00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <math.h>
00034 #include "../defs.h"
00035 #include "../athena.h"
00036 #include "../prototypes.h"
00037 #include "prototypes.h"
00038 #include "particle.h"
00039 #include "../globals.h"
00040
00041
00042 #ifdef FEEDBACK
00043
00044
00045 #define NVAR_F 4
00046
00047 #define NGF 2
00048
00049 #ifdef MPI_PARALLEL
00050
00051
00052
00053 static double *send_buf = NULL, *recv_buf = NULL;
00054
00055 #endif
00056
00057
00058 static int my_iproc, my_jproc, my_kproc;
00059
00060 static int il,iu, jl,ju, kl,ku;
00061
00062 static Real x1min,x1max,x2min,x2max,x3min,x3max;
00063
00064 static Real Lx1, Lx2, Lx3;
00065
00066
00067 static VBCFun_t apply_ix1 = NULL, apply_ox1 = NULL;
00068 static VBCFun_t apply_ix2 = NULL, apply_ox2 = NULL;
00069 static VBCFun_t apply_ix3 = NULL, apply_ox3 = NULL;
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 static void reflect_ix1_feedback(Grid *pG);
00084 static void reflect_ox1_feedback(Grid *pG);
00085 static void reflect_ix2_feedback(Grid *pG);
00086 static void reflect_ox2_feedback(Grid *pG);
00087 static void reflect_ix3_feedback(Grid *pG);
00088 static void reflect_ox3_feedback(Grid *pG);
00089
00090 static void outflow_feedback(Grid *pG);
00091
00092 static void periodic_ix1_feedback(Grid *pG);
00093 static void periodic_ox1_feedback(Grid *pG);
00094 static void periodic_ix2_feedback(Grid *pG);
00095 static void periodic_ox2_feedback(Grid *pG);
00096 static void periodic_ix3_feedback(Grid *pG);
00097 static void periodic_ox3_feedback(Grid *pG);
00098
00099 #ifdef MPI_PARALLEL
00100 static void send_ix1_feedback(Grid *pG);
00101 static void send_ox1_feedback(Grid *pG);
00102 static void send_ix2_feedback(Grid *pG);
00103 static void send_ox2_feedback(Grid *pG);
00104 static void send_ix3_feedback(Grid *pG);
00105 static void send_ox3_feedback(Grid *pG);
00106
00107 static void recv_ix1_feedback(Grid *pG, MPI_Request *prq);
00108 static void recv_ox1_feedback(Grid *pG, MPI_Request *prq);
00109 static void recv_ix2_feedback(Grid *pG, MPI_Request *prq);
00110 static void recv_ox2_feedback(Grid *pG, MPI_Request *prq);
00111 static void recv_ix3_feedback(Grid *pG, MPI_Request *prq);
00112 static void recv_ox3_feedback(Grid *pG, MPI_Request *prq);
00113 #endif
00114
00115 #ifdef SHEARING_BOX
00116 static void shearingbox_ix1_feedback(Grid *pG, Domain *pD);
00117 static void shearingbox_ox1_feedback(Grid *pG, Domain *pD);
00118 #endif
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 void exchange_feedback(Grid *pG, Domain *pD)
00138 {
00139 #ifdef MPI_PARALLEL
00140 int cnt1, cnt2, cnt3, cnt, err;
00141 MPI_Request rq;
00142 #endif
00143
00144
00145
00146
00147 if (pG->Nx3 > 1){
00148
00149 #ifdef MPI_PARALLEL
00150 cnt1 = pG->Nx1 > 1 ? pG->Nx1 + 2*NGF : 1;
00151 cnt2 = pG->Nx2 > 1 ? pG->Nx2 + 2*NGF : 1;
00152 cnt = NGF*cnt1*cnt2*NVAR_F;
00153
00154
00155 if (pG->rx3_id >= 0 && pG->lx3_id >= 0) {
00156
00157 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx3_id,
00158 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00159 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00160
00161 send_ox3_feedback(pG);
00162 recv_ix3_feedback(pG, &rq);
00163
00164
00165 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx3_id,
00166 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00167 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00168
00169 send_ix3_feedback(pG);
00170 recv_ox3_feedback(pG, &rq);
00171 }
00172
00173
00174 if (pG->rx3_id >= 0 && pG->lx3_id < 0) {
00175
00176 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx3_id,
00177 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00178 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00179
00180 send_ox3_feedback(pG);
00181 (*apply_ix3) (pG);
00182 recv_ox3_feedback(pG, &rq);
00183 }
00184
00185
00186 if (pG->rx3_id < 0 && pG->lx3_id >= 0) {
00187
00188 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx3_id,
00189 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00190 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00191
00192 send_ix3_feedback(pG);
00193 (*apply_ox3) (pG);
00194 recv_ix3_feedback(pG, &rq);
00195 }
00196 #endif
00197
00198
00199 if (pG->rx3_id < 0 && pG->lx3_id < 0) {
00200 (*apply_ix3)(pG);
00201 (*apply_ox3)(pG);
00202 }
00203
00204 }
00205
00206
00207
00208
00209 if (pG->Nx2 > 1){
00210
00211 #ifdef MPI_PARALLEL
00212 cnt1 = pG->Nx1 > 1 ? pG->Nx1 + 2*NGF : 1;
00213 cnt3 = pG->Nx3 > 1 ? pG->Nx3 : 1;
00214 cnt = NGF*cnt1*cnt3*NVAR_F;
00215
00216
00217 if (pG->rx2_id >= 0 && pG->lx2_id >= 0) {
00218
00219 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx2_id,
00220 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00221 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00222
00223 send_ox2_feedback(pG);
00224 recv_ix2_feedback(pG, &rq);
00225
00226
00227 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx2_id,
00228 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00229 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00230
00231 send_ix2_feedback(pG);
00232 recv_ox2_feedback(pG, &rq);
00233 }
00234
00235
00236 if (pG->rx2_id >= 0 && pG->lx2_id < 0) {
00237
00238 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx2_id,
00239 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00240 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00241
00242 send_ox2_feedback(pG);
00243 (*apply_ix2) (pG);
00244 recv_ox2_feedback(pG, &rq);
00245 }
00246
00247
00248 if (pG->rx2_id < 0 && pG->lx2_id >= 0) {
00249
00250 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx2_id,
00251 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00252 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00253
00254 send_ix2_feedback(pG);
00255 (*apply_ox2) (pG);
00256 recv_ix2_feedback(pG, &rq);
00257 }
00258 #endif
00259
00260
00261 if (pG->rx2_id < 0 && pG->lx2_id < 0) {
00262 (*apply_ix2)(pG);
00263 (*apply_ox2)(pG);
00264 }
00265
00266 }
00267
00268
00269
00270 #ifdef SHEARING_BOX
00271 #ifndef FARGO
00272 if (pG->Nx3 > 1)
00273 {
00274 if (my_iproc == 0)
00275 shearingbox_ix1_feedback(pG, pD);
00276
00277 if (my_iproc == pD->NGrid_x1-1)
00278 shearingbox_ox1_feedback(pG, pD);
00279 }
00280 #endif
00281 #endif
00282
00283
00284
00285
00286 if (pG->Nx1 > 1){
00287
00288 #ifdef MPI_PARALLEL
00289 cnt2 = pG->Nx2 > 1 ? pG->Nx2 : 1;
00290 cnt3 = pG->Nx3 > 1 ? pG->Nx3 : 1;
00291 cnt = NGF*cnt2*cnt3*NVAR_F;
00292
00293
00294 if (pG->rx1_id >= 0 && pG->lx1_id >= 0) {
00295
00296 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx1_id,
00297 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00298 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00299
00300 send_ox1_feedback(pG);
00301 recv_ix1_feedback(pG, &rq);
00302
00303
00304 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx1_id,
00305 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00306 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00307
00308 send_ix1_feedback(pG);
00309 recv_ox1_feedback(pG, &rq);
00310 }
00311
00312
00313 if (pG->rx1_id >= 0 && pG->lx1_id < 0) {
00314
00315 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->rx1_id,
00316 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00317 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00318
00319 send_ox1_feedback(pG);
00320 (*apply_ix1) (pG);
00321 recv_ox1_feedback(pG, &rq);
00322 }
00323
00324
00325 if (pG->rx1_id < 0 && pG->lx1_id >= 0) {
00326
00327 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pG->lx1_id,
00328 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00329 if(err) ath_error("[exchange_feedback]: MPI_Irecv error = %d\n",err);
00330
00331 send_ix1_feedback(pG);
00332 (*apply_ox1) (pG);
00333 recv_ix1_feedback(pG, &rq);
00334 }
00335 #endif
00336
00337
00338 if (pG->rx1_id < 0 && pG->lx1_id < 0) {
00339 (*apply_ix1)(pG);
00340 (*apply_ox1)(pG);
00341 }
00342
00343 }
00344
00345 return;
00346
00347 }
00348
00349
00350
00351
00352
00353
00354 void exchange_feedback_init(Grid *pG, Domain *pD)
00355 {
00356 int ibc_x1, obc_x1;
00357 int ibc_x2, obc_x2;
00358 int ibc_x3, obc_x3;
00359 #ifdef MPI_PARALLEL
00360 int i,j,k;
00361 int x1cnt, x2cnt, x3cnt;
00362 int nx1t, nx2t, nx3t, size;
00363 #endif
00364
00365
00366 if (pG->Nx1 > 1) {
00367 il = pG->is - NGF; iu = pG->ie + NGF;
00368 } else {
00369 il = pG->is; iu = pG->ie;
00370 }
00371
00372 if (pG->Nx2 > 1) {
00373 jl = pG->js - NGF; ju = pG->je + NGF;
00374 } else {
00375 jl = pG->js; ju = pG->je;
00376 }
00377
00378 if (pG->Nx3 > 1) {
00379 kl = pG->ks - NGF; ku = pG->ke + NGF;
00380 } else {
00381 kl = pG->ks; ku = pG->ke;
00382 }
00383
00384
00385 x1min = par_getd("grid","x1min");
00386 x1max = par_getd("grid","x1max");
00387 Lx1 = x1max - x1min;
00388
00389 x2min = par_getd("grid","x2min");
00390 x2max = par_getd("grid","x2max");
00391 Lx2 = x2max - x2min;
00392
00393 x3min = par_getd("grid","x3min");
00394 x3max = par_getd("grid","x3max");
00395 Lx3 = x3max - x3min;
00396
00397 get_myGridIndex(pD, pG->my_id, &my_iproc, &my_jproc, &my_kproc);
00398
00399
00400
00401 if(pG->Nx1 > 1) {
00402 if(apply_ix1 == NULL){
00403
00404 ibc_x1 = par_geti("grid","ibc_x1");
00405 switch(ibc_x1){
00406
00407 case 1:
00408 apply_ix1 = reflect_ix1_feedback;
00409 break;
00410 case 5:
00411 apply_ix1 = reflect_ix1_feedback;
00412 break;
00413
00414 case 2:
00415 apply_ix1 = outflow_feedback;
00416 break;
00417
00418 case 4:
00419 apply_ix1 = periodic_ix1_feedback;
00420 break;
00421
00422 default:
00423 ath_perr(-1,"[exchange_feedback_init]: ibc_x1 = %d unknown\n",ibc_x1);
00424 exit(EXIT_FAILURE);
00425 }
00426 }
00427
00428 if(apply_ox1 == NULL){
00429
00430 obc_x1 = par_geti("grid","obc_x1");
00431 switch(obc_x1){
00432
00433 case 1:
00434 apply_ox1 = reflect_ox1_feedback;
00435 break;
00436 case 5:
00437 apply_ox1 = reflect_ox1_feedback;
00438 break;
00439
00440 case 2:
00441 apply_ox1 = outflow_feedback;
00442 break;
00443
00444 case 4:
00445 apply_ox1 = periodic_ox1_feedback;
00446 break;
00447
00448 default:
00449 ath_perr(-1,"[exchange_feedback_init]: obc_x1 = %d unknown\n",obc_x1);
00450 exit(EXIT_FAILURE);
00451 }
00452 }
00453 }
00454
00455
00456
00457 if(pG->Nx2 > 1) {
00458 if(apply_ix2 == NULL){
00459
00460 ibc_x2 = par_geti("grid","ibc_x2");
00461 switch(ibc_x2){
00462
00463 case 1:
00464 apply_ix2 = reflect_ix2_feedback;
00465 break;
00466 case 5:
00467 apply_ix2 = reflect_ix2_feedback;
00468 break;
00469
00470 case 2:
00471 apply_ix2 = outflow_feedback;
00472 break;
00473
00474 case 4:
00475 apply_ix2 = periodic_ix2_feedback;
00476 break;
00477
00478 default:
00479 ath_perr(-1,"[exchange_feedback_init]: ibc_x2 = %d unknown\n",ibc_x2);
00480 exit(EXIT_FAILURE);
00481 }
00482 }
00483
00484 if(apply_ox2 == NULL){
00485
00486 obc_x2 = par_geti("grid","obc_x2");
00487 switch(obc_x2){
00488
00489 case 1:
00490 apply_ox2 = reflect_ox2_feedback;
00491 break;
00492 case 5:
00493 apply_ox2 = reflect_ox2_feedback;
00494 break;
00495
00496 case 2:
00497 apply_ox2 = outflow_feedback;
00498 break;
00499
00500 case 4:
00501 apply_ox2 = periodic_ox2_feedback;
00502 break;
00503
00504 default:
00505 ath_perr(-1,"[exchange_feedback_init]: obc_x2 = %d unknown\n",obc_x2);
00506 exit(EXIT_FAILURE);
00507 }
00508 }
00509 }
00510
00511
00512
00513 if(pG->Nx3 > 1) {
00514 if(apply_ix3 == NULL){
00515
00516 ibc_x3 = par_geti("grid","ibc_x3");
00517 switch(ibc_x3){
00518
00519 case 1:
00520 apply_ix3 = reflect_ix3_feedback;
00521 break;
00522 case 5:
00523 apply_ix3 = reflect_ix3_feedback;
00524 break;
00525
00526 case 2:
00527 apply_ix3 = outflow_feedback;
00528 break;
00529
00530 case 4:
00531 apply_ix3 = periodic_ix3_feedback;
00532 break;
00533
00534 default:
00535 ath_perr(-1,"[exchange_feedback_init]: ibc_x3 = %d unknown\n",ibc_x3);
00536 exit(EXIT_FAILURE);
00537 }
00538 }
00539
00540 if(apply_ox3 == NULL){
00541
00542 obc_x3 = par_geti("grid","obc_x3");
00543 switch(obc_x3){
00544
00545 case 1:
00546 apply_ox3 = reflect_ox3_feedback;
00547 break;
00548 case 5:
00549 apply_ox3 = reflect_ox3_feedback;
00550 break;
00551
00552 case 2:
00553 apply_ox3 = outflow_feedback;
00554 break;
00555
00556 case 4:
00557 apply_ox3 = periodic_ox3_feedback;
00558 break;
00559
00560 default:
00561 ath_perr(-1,"[exchange_feedback_init]: obc_x3 = %d unknown\n",obc_x3);
00562 exit(EXIT_FAILURE);
00563 }
00564 }
00565 }
00566
00567 #ifdef MPI_PARALLEL
00568 x1cnt = x2cnt = x3cnt = 0;
00569
00570 for (k=0; k<(pD->NGrid_x3); k++){
00571 for (j=0; j<(pD->NGrid_x2); j++){
00572 for (i=0; i<(pD->NGrid_x1); i++){
00573 if(pD->NGrid_x3 > 1){
00574 nx1t = pD->GridArray[k][j][i].ige - pD->GridArray[k][j][i].igs + 1;
00575 if(nx1t > 1) nx1t += 2*NGF;
00576
00577 nx2t = pD->GridArray[k][j][i].jge - pD->GridArray[k][j][i].jgs + 1;
00578 if(nx2t > 1) nx2t += 2*NGF;
00579
00580 x3cnt = nx1t*nx2t > x3cnt ? nx1t*nx2t : x3cnt;
00581 }
00582
00583 if(pD->NGrid_x2 > 1){
00584 nx1t = pD->GridArray[k][j][i].ige - pD->GridArray[k][j][i].igs + 1;
00585 if(nx1t > 1) nx1t += 2*NGF;
00586
00587 nx3t = pD->GridArray[k][j][i].kge - pD->GridArray[k][j][i].kgs + 1;
00588 if(nx3t > 1) nx3t += 1;
00589
00590 x2cnt = nx1t*nx3t > x2cnt ? nx1t*nx3t : x2cnt;
00591 }
00592
00593 if(pD->NGrid_x1 > 1){
00594 nx2t = pD->GridArray[k][j][i].jge - pD->GridArray[k][j][i].jgs + 1;
00595 if(nx2t > 1) nx2t += 1;
00596
00597 nx3t = pD->GridArray[k][j][i].kge - pD->GridArray[k][j][i].kgs + 1;
00598 if(nx3t > 1) nx3t += 1;
00599
00600 x1cnt = nx2t*nx3t > x1cnt ? nx2t*nx3t : x1cnt;
00601 }
00602
00603 }
00604 }
00605 }
00606
00607 size = x1cnt > x2cnt ? x1cnt : x2cnt;
00608 size = x3cnt > size ? x3cnt : size;
00609
00610 size *= NGF;
00611
00612 if (size > 0) {
00613 if((send_buf = (double*)malloc(size*NVAR_F*sizeof(double))) == NULL)
00614 ath_error("[exchange_feedback_init]: Failed to allocate send buffer\n");
00615
00616 if((recv_buf = (double*)malloc(size*NVAR_F*sizeof(double))) == NULL)
00617 ath_error("[exchange_feedback_init]: Failed to allocate recv buffer\n");
00618 }
00619 #endif
00620
00621 return;
00622 }
00623
00624
00625
00626
00627
00628
00629
00630 void exchange_feedback_fun(enum Direction dir, VBCFun_t prob_bc)
00631 {
00632 switch(dir){
00633 case left_x1:
00634 apply_ix1 = prob_bc;
00635 break;
00636 case right_x1:
00637 apply_ox1 = prob_bc;
00638 break;
00639 case left_x2:
00640 apply_ix2 = prob_bc;
00641 break;
00642 case right_x2:
00643 apply_ox2 = prob_bc;
00644 break;
00645 case left_x3:
00646 apply_ix3 = prob_bc;
00647 break;
00648 case right_x3:
00649 apply_ox3 = prob_bc;
00650 break;
00651 default:
00652 ath_perr(-1,"[set_bvals_particle_fun]: Unknown direction = %d\n",dir);
00653 exit(EXIT_FAILURE);
00654 }
00655 return;
00656 }
00657
00658
00659
00660 void exchange_feedback_destruct(Grid *pG, Domain *pD)
00661 {
00662 apply_ix1 = NULL;
00663 apply_ox1 = NULL;
00664 apply_ix2 = NULL;
00665 apply_ox2 = NULL;
00666 apply_ix3 = NULL;
00667 apply_ox3 = NULL;
00668 #ifdef MPI_PARALLEL
00669 free(send_buf);
00670 free(recv_buf);
00671 #endif
00672 return;
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690 static void reflect_ix3_feedback(Grid *pG)
00691 {
00692 int kr;
00693 int i,j,k;
00694
00695 for (k=nghost-NGF; k<nghost; k++) {
00696 kr = 2*nghost-k-1;
00697 for (j=jl; j<=ju; j++) {
00698 for (i=il; i<=iu; i++) {
00699 pG->Coup[kr][j][i].fb1 += pG->Coup[k][j][i].fb1;
00700 pG->Coup[kr][j][i].fb2 += pG->Coup[k][j][i].fb2;
00701 pG->Coup[kr][j][i].fb3 -= pG->Coup[k][j][i].fb3;
00702
00703 pG->Coup[kr][j][i].Eloss += pG->Coup[k][j][i].Eloss;
00704 }
00705 }
00706 }
00707
00708 return;
00709 }
00710
00711
00712
00713
00714
00715
00716 static void reflect_ox3_feedback(Grid *pG)
00717 {
00718 int kr;
00719 int i,j,k;
00720
00721 for (k=pG->ke+1; k<=pG->ke+NGF; k++) {
00722 kr = 2*pG->ke-k+1;
00723 for (j=jl; j<=ju; j++) {
00724 for (i=il; i<=iu; i++) {
00725 pG->Coup[kr][j][i].fb1 += pG->Coup[k][j][i].fb1;
00726 pG->Coup[kr][j][i].fb2 += pG->Coup[k][j][i].fb2;
00727 pG->Coup[kr][j][i].fb3 -= pG->Coup[k][j][i].fb3;
00728
00729 pG->Coup[kr][j][i].Eloss += pG->Coup[k][j][i].Eloss;
00730 }
00731 }
00732 }
00733
00734 return;
00735 }
00736
00737
00738
00739
00740
00741 static void reflect_ix2_feedback(Grid *pG)
00742 {
00743 int jr;
00744 int i,j,k;
00745
00746 for (k=pG->ks; k<=pG->ke; k++) {
00747 for (j=nghost-NGF; j<nghost; j++) {
00748 jr = 2*nghost-j-1;
00749 for (i=il; i<=iu; i++) {
00750 pG->Coup[k][jr][i].fb1 += pG->Coup[k][j][i].fb1;
00751 pG->Coup[k][jr][i].fb2 -= pG->Coup[k][j][i].fb2;
00752 pG->Coup[k][jr][i].fb3 += pG->Coup[k][j][i].fb3;
00753
00754 pG->Coup[k][jr][i].Eloss += pG->Coup[k][j][i].Eloss;
00755 }
00756 }
00757 }
00758
00759 return;
00760 }
00761
00762
00763
00764
00765
00766
00767 static void reflect_ox2_feedback(Grid *pG)
00768 {
00769 int jr;
00770 int i,j,k;
00771
00772 for (k=pG->ks; k<=pG->ke; k++) {
00773 for (j=pG->je+1; j<=pG->je+NGF; j++) {
00774 jr = 2*pG->je-j+1;
00775 for (i=il; i<=iu; i++) {
00776 pG->Coup[k][jr][i].fb1 += pG->Coup[k][j][i].fb1;
00777 pG->Coup[k][jr][i].fb2 -= pG->Coup[k][j][i].fb2;
00778 pG->Coup[k][jr][i].fb3 += pG->Coup[k][j][i].fb3;
00779
00780 pG->Coup[k][jr][i].Eloss += pG->Coup[k][j][i].Eloss;
00781 }
00782 }
00783 }
00784
00785 return;
00786 }
00787
00788
00789
00790
00791
00792
00793 static void reflect_ix1_feedback(Grid *pG)
00794 {
00795 int ir;
00796 int i,j,k;
00797
00798 for (k=pG->ks; k<=pG->ke; k++) {
00799 for (j=pG->js; j<=pG->je; j++) {
00800 for (i=nghost-NGF; i<nghost; i++) {
00801 ir = 2*nghost-i-1;
00802 pG->Coup[k][j][ir].fb1 -= pG->Coup[k][j][i].fb1;
00803 pG->Coup[k][j][ir].fb2 += pG->Coup[k][j][i].fb2;
00804 pG->Coup[k][j][ir].fb3 += pG->Coup[k][j][i].fb3;
00805
00806 pG->Coup[k][j][ir].Eloss += pG->Coup[k][j][i].Eloss;
00807 }
00808 }
00809 }
00810
00811 return;
00812 }
00813
00814
00815
00816
00817
00818
00819 static void reflect_ox1_feedback(Grid *pG)
00820 {
00821 int ir;
00822 int i,j,k;
00823
00824 for (k=pG->ks; k<=pG->ke; k++) {
00825 for (j=pG->js; j<=pG->je; j++) {
00826 for (i=pG->ie+1; i<=pG->ie+NGF; i++) {
00827 ir = 2*pG->ie-i+1;
00828 pG->Coup[k][j][ir].fb1 -= pG->Coup[k][j][i].fb1;
00829 pG->Coup[k][j][ir].fb2 += pG->Coup[k][j][i].fb2;
00830 pG->Coup[k][j][ir].fb3 += pG->Coup[k][j][i].fb3;
00831
00832 pG->Coup[k][j][ir].Eloss += pG->Coup[k][j][i].Eloss;
00833 }
00834 }
00835 }
00836
00837 return;
00838 }
00839
00840
00841
00842
00843
00844
00845 static void outflow_feedback(Grid *pG)
00846 {
00847 return;
00848 }
00849
00850
00851
00852
00853
00854 static void periodic_ix3_feedback(Grid *pG)
00855 {
00856 int dk = pG->ke - pG->ks + 1;
00857 int i,j,k;
00858
00859 for (k=nghost; k<nghost+NGF; k++) {
00860 for (j=jl; j<=ju; j++) {
00861 for (i=il; i<=iu; i++) {
00862 pG->Coup[k][j][i].fb1 += pG->Coup[k+dk][j][i].fb1;
00863 pG->Coup[k][j][i].fb2 += pG->Coup[k+dk][j][i].fb2;
00864 pG->Coup[k][j][i].fb3 += pG->Coup[k+dk][j][i].fb3;
00865
00866 pG->Coup[k][j][i].Eloss += pG->Coup[k+dk][j][i].Eloss;
00867 }
00868 }
00869 }
00870
00871 return;
00872 }
00873
00874
00875
00876
00877
00878
00879 static void periodic_ox3_feedback(Grid *pG)
00880 {
00881 int dk = pG->ke - pG->ks + 1;
00882 int i,j,k;
00883
00884 for (k=nghost-NGF; k<nghost; k++) {
00885 for (j=jl; j<=ju; j++) {
00886 for (i=il; i<=iu; i++) {
00887 pG->Coup[k+dk][j][i].fb1 += pG->Coup[k][j][i].fb1;
00888 pG->Coup[k+dk][j][i].fb2 += pG->Coup[k][j][i].fb2;
00889 pG->Coup[k+dk][j][i].fb3 += pG->Coup[k][j][i].fb3;
00890
00891 pG->Coup[k+dk][j][i].Eloss += pG->Coup[k][j][i].Eloss;
00892 }
00893 }
00894 }
00895
00896 return;
00897 }
00898
00899
00900
00901
00902
00903
00904 static void periodic_ix2_feedback(Grid *pG)
00905 {
00906 int dj = pG->je - pG->js + 1;
00907 int i,j,k;
00908
00909 for (k=pG->ks; k<=pG->ke; k++) {
00910 for (j=nghost; j<nghost+NGF; j++) {
00911 for (i=il; i<=iu; i++) {
00912 pG->Coup[k][j][i].fb1 += pG->Coup[k][j+dj][i].fb1;
00913 pG->Coup[k][j][i].fb2 += pG->Coup[k][j+dj][i].fb2;
00914 pG->Coup[k][j][i].fb3 += pG->Coup[k][j+dj][i].fb3;
00915
00916 pG->Coup[k][j][i].Eloss += pG->Coup[k][j+dj][i].Eloss;
00917 }
00918 }
00919 }
00920
00921 return;
00922 }
00923
00924
00925
00926
00927
00928
00929 static void periodic_ox2_feedback(Grid *pG)
00930 {
00931 int dj = pG->je - pG->js + 1;
00932 int i,j,k;
00933
00934 for (k=pG->ks; k<=pG->ke; k++) {
00935 for (j=nghost-NGF; j<nghost; j++) {
00936 for (i=il; i<=iu; i++) {
00937 pG->Coup[k][j+dj][i].fb1 += pG->Coup[k][j][i].fb1;
00938 pG->Coup[k][j+dj][i].fb2 += pG->Coup[k][j][i].fb2;
00939 pG->Coup[k][j+dj][i].fb3 += pG->Coup[k][j][i].fb3;
00940
00941 pG->Coup[k][j+dj][i].Eloss += pG->Coup[k][j][i].Eloss;
00942 }
00943 }
00944 }
00945
00946 return;
00947 }
00948
00949
00950
00951
00952
00953
00954 static void periodic_ix1_feedback(Grid *pG)
00955 {
00956 int di = pG->ie - pG->is + 1;
00957 int i,j,k;
00958
00959 for (k=pG->ks; k<=pG->ke; k++) {
00960 for (j=pG->js; j<=pG->je; j++) {
00961 for (i=nghost; i<nghost+NGF; i++) {
00962 pG->Coup[k][j][i].fb1 += pG->Coup[k][j][i+di].fb1;
00963 pG->Coup[k][j][i].fb2 += pG->Coup[k][j][i+di].fb2;
00964 pG->Coup[k][j][i].fb3 += pG->Coup[k][j][i+di].fb3;
00965
00966 pG->Coup[k][j][i].Eloss += pG->Coup[k][j][i+di].Eloss;
00967 }
00968 }
00969 }
00970
00971 return;
00972 }
00973
00974
00975
00976
00977
00978
00979 static void periodic_ox1_feedback(Grid *pG)
00980 {
00981 int di = pG->ie - pG->is + 1;
00982 int i,j,k;
00983
00984 for (k=pG->ks; k<=pG->ke; k++) {
00985 for (j=pG->js; j<=pG->je; j++) {
00986 for (i=nghost-NGF; i<nghost; i++) {
00987 pG->Coup[k][j][i+di].fb1 += pG->Coup[k][j][i].fb1;
00988 pG->Coup[k][j][i+di].fb2 += pG->Coup[k][j][i].fb2;
00989 pG->Coup[k][j][i+di].fb3 += pG->Coup[k][j][i].fb3;
00990
00991 pG->Coup[k][j][i+di].Eloss += pG->Coup[k][j][i].Eloss;
00992 }
00993 }
00994 }
00995
00996 return;
00997 }
00998
00999 #ifdef MPI_PARALLEL
01000
01001
01002
01003
01004
01005
01006 static void send_ix3_feedback(Grid *pG)
01007 {
01008 int i,j,k,cnt,err;
01009 GPCouple *pq;
01010 double *pd = send_buf;
01011
01012
01013 cnt = (iu-il+1)*(ju-jl+1)*NGF*NVAR_F;
01014
01015 for (k=nghost-NGF; k<nghost; k++) {
01016 for (j=jl; j<=ju; j++) {
01017 for (i=il; i<=iu; i++) {
01018
01019 pq = &(pG->Coup[k][j][i]);
01020
01021 *(pd++) = pq->fb1;
01022 *(pd++) = pq->fb2;
01023 *(pd++) = pq->fb3;
01024
01025 *(pd++) = pq->Eloss;
01026 }
01027 }
01028 }
01029
01030
01031 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->lx3_id,
01032 boundary_cells_tag, MPI_COMM_WORLD);
01033 if(err) ath_error("[send_ix3_feedback]: MPI_Send error = %d\n",err);
01034
01035 return;
01036 }
01037
01038
01039
01040
01041
01042
01043 static void send_ox3_feedback(Grid *pG)
01044 {
01045 int i,j,k,cnt,err;
01046 GPCouple *pq;
01047 double *pd = send_buf;
01048
01049
01050 cnt = (iu-il+1)*(ju-jl+1)*NGF*NVAR_F;
01051
01052 for (k=pG->ke+1; k<=pG->ke+NGF; k++) {
01053 for (j=jl; j<=ju; j++) {
01054 for (i=il; i<=iu; i++) {
01055
01056 pq = &(pG->Coup[k][j][i]);
01057
01058 *(pd++) = pq->fb1;
01059 *(pd++) = pq->fb2;
01060 *(pd++) = pq->fb3;
01061
01062 *(pd++) = pq->Eloss;
01063 }
01064 }
01065 }
01066
01067
01068 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->rx3_id,
01069 boundary_cells_tag, MPI_COMM_WORLD);
01070 if(err) ath_error("[send_ox3_feedback]: MPI_Send error = %d\n",err);
01071
01072 return;
01073 }
01074
01075
01076
01077
01078
01079
01080 static void send_ix2_feedback(Grid *pG)
01081 {
01082 int i,j,k,cnt,err;
01083 GPCouple *pq;
01084 double *pd = send_buf;
01085
01086
01087 cnt = (iu-il+1)*(pG->ke-pG->ks+1)*NGF*NVAR_F;
01088
01089 for (k=pG->ks; k<=pG->ke; k++) {
01090 for (j=nghost-NGF; j<nghost; j++) {
01091 for (i=il; i<=iu; i++) {
01092
01093 pq = &(pG->Coup[k][j][i]);
01094
01095 *(pd++) = pq->fb1;
01096 *(pd++) = pq->fb2;
01097 *(pd++) = pq->fb3;
01098
01099 *(pd++) = pq->Eloss;
01100 }
01101 }
01102 }
01103
01104
01105 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->lx2_id,
01106 boundary_cells_tag, MPI_COMM_WORLD);
01107 if(err) ath_error("[send_ix2_feedback]: MPI_Send error = %d\n",err);
01108
01109 return;
01110 }
01111
01112
01113
01114
01115
01116
01117 static void send_ox2_feedback(Grid *pG)
01118 {
01119 int i,j,k,cnt,err;
01120 GPCouple *pq;
01121 double *pd = send_buf;
01122
01123
01124 cnt = (iu-il+1)*(pG->ke-pG->ks+1)*NGF*NVAR_F;
01125
01126 for (k=pG->ks; k<=pG->ke; k++) {
01127 for (j=pG->je+1; j<=pG->je+NGF; j++) {
01128 for (i=il; i<=iu; i++) {
01129
01130 pq = &(pG->Coup[k][j][i]);
01131
01132 *(pd++) = pq->fb1;
01133 *(pd++) = pq->fb2;
01134 *(pd++) = pq->fb3;
01135
01136 *(pd++) = pq->Eloss;
01137 }
01138 }
01139 }
01140
01141
01142 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->rx2_id,
01143 boundary_cells_tag, MPI_COMM_WORLD);
01144 if(err) ath_error("[send_ox2_feedback]: MPI_Send error = %d\n",err);
01145
01146 return;
01147 }
01148
01149
01150
01151
01152
01153
01154 static void send_ix1_feedback(Grid *pG)
01155 {
01156 int i,j,k,cnt,err;
01157 GPCouple *pq;
01158 double *pd = send_buf;
01159
01160
01161 cnt = (pG->je-pG->js+1)*(pG->ke-pG->ks+1)*NGF*NVAR_F;
01162
01163 for (k=pG->ks; k<=pG->ke; k++) {
01164 for (j=pG->js; j<=pG->je; j++) {
01165 for (i=nghost-NGF; i<nghost; i++) {
01166
01167 pq = &(pG->Coup[k][j][i]);
01168
01169 *(pd++) = pq->fb1;
01170 *(pd++) = pq->fb2;
01171 *(pd++) = pq->fb3;
01172
01173 *(pd++) = pq->Eloss;
01174 }
01175 }
01176 }
01177
01178
01179 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->lx1_id,
01180 boundary_cells_tag, MPI_COMM_WORLD);
01181 if(err) ath_error("[send_ix1_feedback]: MPI_Send error = %d\n",err);
01182
01183 return;
01184 }
01185
01186
01187
01188
01189
01190
01191 static void send_ox1_feedback(Grid *pG)
01192 {
01193 int i,j,k,cnt,err;
01194 GPCouple *pq;
01195 double *pd = send_buf;
01196
01197
01198 cnt = (pG->je-pG->js+1)*(pG->ke-pG->ks+1)*NGF*NVAR_F;
01199
01200 for (k=pG->ks; k<=pG->ke; k++) {
01201 for (j=pG->js; j<=pG->je; j++) {
01202 for (i=pG->ie+1; i<=pG->ie+NGF; i++) {
01203
01204 pq = &(pG->Coup[k][j][i]);
01205
01206 *(pd++) = pq->fb1;
01207 *(pd++) = pq->fb2;
01208 *(pd++) = pq->fb3;
01209
01210 *(pd++) = pq->Eloss;
01211 }
01212 }
01213 }
01214
01215
01216 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pG->rx1_id,
01217 boundary_cells_tag, MPI_COMM_WORLD);
01218 if(err) ath_error("[send_ox1_feedback]: MPI_Send error = %d\n",err);
01219
01220 return;
01221 }
01222
01223
01224
01225
01226
01227
01228 static void recv_ix3_feedback(Grid *pG, MPI_Request *prq)
01229 {
01230 int i,j,k,err;
01231 GPCouple *pq;
01232 MPI_Status stat;
01233 double *pd = recv_buf;
01234
01235
01236 err = MPI_Wait(prq, &stat);
01237 if(err) ath_error("[receive_ix3_feedback]: MPI_Wait error = %d\n",err);
01238
01239
01240 for (k=pG->ks; k<pG->ks+NGF; k++){
01241 for (j=jl; j<=ju; j++){
01242 for (i=il; i<=iu; i++){
01243
01244 pq = &(pG->Coup[k][j][i]);
01245
01246 pq->fb1 += *(pd++);
01247 pq->fb2 += *(pd++);
01248 pq->fb3 += *(pd++);
01249
01250 pq->Eloss += *(pd++);
01251 }
01252 }
01253 }
01254
01255 return;
01256 }
01257
01258
01259
01260
01261
01262
01263 static void recv_ox3_feedback(Grid *pG, MPI_Request *prq)
01264 {
01265 int i,j,k,err;
01266 GPCouple *pq;
01267 MPI_Status stat;
01268 double *pd = recv_buf;
01269
01270
01271 err = MPI_Wait(prq, &stat);
01272 if(err) ath_error("[receive_ox3_feedback]: MPI_Wait error = %d\n",err);
01273
01274
01275 for (k=pG->ke-NGF+1; k<=pG->ke; k++){
01276 for (j=jl; j<=ju; j++){
01277 for (i=il; i<=iu; i++){
01278
01279 pq = &(pG->Coup[k][j][i]);
01280
01281 pq->fb1 += *(pd++);
01282 pq->fb2 += *(pd++);
01283 pq->fb3 += *(pd++);
01284
01285 pq->Eloss += *(pd++);
01286 }
01287 }
01288 }
01289
01290 return;
01291 }
01292
01293
01294
01295
01296
01297
01298 static void recv_ix2_feedback(Grid *pG, MPI_Request *prq)
01299 {
01300 int i,j,k,err;
01301 GPCouple *pq;
01302 MPI_Status stat;
01303 double *pd = recv_buf;
01304
01305
01306 err = MPI_Wait(prq, &stat);
01307 if(err) ath_error("[receive_ix2_feedback]: MPI_Wait error = %d\n",err);
01308
01309
01310 for (k=pG->ks; k<=pG->ke; k++){
01311 for (j=pG->js; j<pG->js+NGF; j++){
01312 for (i=il; i<=iu; i++){
01313
01314 pq = &(pG->Coup[k][j][i]);
01315
01316 pq->fb1 += *(pd++);
01317 pq->fb2 += *(pd++);
01318 pq->fb3 += *(pd++);
01319
01320 pq->Eloss += *(pd++);
01321 }
01322 }
01323 }
01324
01325 return;
01326 }
01327
01328
01329
01330
01331
01332
01333 static void recv_ox2_feedback(Grid *pG, MPI_Request *prq)
01334 {
01335 int i,j,k,err;
01336 GPCouple *pq;
01337 MPI_Status stat;
01338 double *pd = recv_buf;
01339
01340
01341 err = MPI_Wait(prq, &stat);
01342 if(err) ath_error("[receive_ox2_feedback]: MPI_Wait error = %d\n",err);
01343
01344
01345 for (k=pG->ks; k<=pG->ke; k++){
01346 for (j=pG->je-NGF+1; j<=pG->je; j++){
01347 for (i=il; i<=iu; i++){
01348
01349 pq = &(pG->Coup[k][j][i]);
01350
01351 pq->fb1 += *(pd++);
01352 pq->fb2 += *(pd++);
01353 pq->fb3 += *(pd++);
01354
01355 pq->Eloss += *(pd++);
01356 }
01357 }
01358 }
01359
01360 return;
01361 }
01362
01363
01364
01365
01366
01367
01368 static void recv_ix1_feedback(Grid *pG, MPI_Request *prq)
01369 {
01370 int i,j,k,err;
01371 GPCouple *pq;
01372 MPI_Status stat;
01373 double *pd = recv_buf;
01374
01375
01376 err = MPI_Wait(prq, &stat);
01377 if(err) ath_error("[receive_ix1_feedback]: MPI_Wait error = %d\n",err);
01378
01379
01380 for (k=pG->ks; k<=pG->ke; k++){
01381 for (j=pG->js; j<=pG->je; j++){
01382 for (i=pG->is; i<pG->is+NGF; i++){
01383
01384 pq = &(pG->Coup[k][j][i]);
01385
01386 pq->fb1 += *(pd++);
01387 pq->fb2 += *(pd++);
01388 pq->fb3 += *(pd++);
01389
01390 pq->Eloss += *(pd++);
01391 }
01392 }
01393 }
01394
01395 return;
01396 }
01397
01398
01399
01400
01401
01402
01403 static void recv_ox1_feedback(Grid *pG, MPI_Request *prq)
01404 {
01405 int i,j,k,err;
01406 GPCouple *pq;
01407 MPI_Status stat;
01408 double *pd = recv_buf;
01409
01410
01411 err = MPI_Wait(prq, &stat);
01412 if(err) ath_error("[receive_ox1_feedback]: MPI_Wait error = %d\n",err);
01413
01414
01415 for (k=pG->ks; k<=pG->ke; k++){
01416 for (j=pG->js; j<=pG->je; j++){
01417 for (i=pG->ie-NGF+1; i<=pG->ie; i++){
01418
01419 pq = &(pG->Coup[k][j][i]);
01420
01421 pq->fb1 += *(pd++);
01422 pq->fb2 += *(pd++);
01423 pq->fb3 += *(pd++);
01424
01425 pq->Eloss += *(pd++);
01426 }
01427 }
01428 }
01429
01430 return;
01431 }
01432
01433 #endif
01434
01435 #ifdef SHEARING_BOX
01436
01437
01438
01439
01440
01441
01442
01443
01444 static void shearingbox_ix1_feedback(Grid *pG, Domain *pD)
01445 {
01446 return;
01447 }
01448
01449
01450
01451
01452
01453
01454
01455
01456 static void shearingbox_ox1_feedback(Grid *pG, Domain *pD)
01457 {
01458 return;
01459 }
01460
01461 #endif
01462
01463 #endif