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 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <string.h>
00028 #include <stdarg.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033
00034 #ifdef STATIC_MESH_REFINEMENT
00035
00036 static double **send_bufP= NULL;
00037 static double **send_bufRC=NULL;
00038 static double ***recv_bufP= NULL;
00039 #ifdef MPI_PARALLEL
00040 static double ***recv_bufRC=NULL;
00041 static MPI_Request ***recv_rq=NULL;
00042 static MPI_Request **send_rq=NULL;
00043 #endif
00044 static int maxND, *start_addrP;
00045
00046 static ConsS ***GZ[3];
00047 #ifdef MHD
00048 Real **SMRemf1, **SMRemf2, **SMRemf3;
00049 Real3Vect ***BFld[3];
00050 #endif
00051
00052
00053
00054
00055
00056
00057
00058
00059 void ProCon(const ConsS Uim1,const ConsS Ui, const ConsS Uip1,
00060 const ConsS Ujm1,const ConsS Ujp1,
00061 const ConsS Ukm1,const ConsS Ukp1, ConsS PCon[][2][2]);
00062 #ifdef MHD
00063 void ProFld(Real3Vect BGZ[][3][3], Real3Vect PFld[][3][3],
00064 const Real dx1c, const Real dx2c, const Real dx3c);
00065 #endif
00066 #ifndef FIRST_ORDER
00067 static Real mcd_slope(const Real vl, const Real vc, const Real vr);
00068 #endif
00069
00070
00071
00072
00073
00074
00075
00076
00077 void RestrictCorrect(MeshS *pM)
00078 {
00079 GridS *pG;
00080 int nl,nd,ncg,dim,nDim,npg,rbufN,start_addr,cnt,nCons,nFlx;
00081 int i,ii,ics,ice,ips,ipe;
00082 int j,jj,jcs,jce,jps,jpe;
00083 int k,kk,kcs,kce,kps,kpe;
00084 Real q1,q2,q3,fact;
00085 double *pRcv,*pSnd;
00086 GridOvrlpS *pCO, *pPO;
00087 #if (NSCALARS > 0)
00088 int n;
00089 #endif
00090 #ifdef MHD
00091 int ib,jb,kb,nFld=0;
00092 #endif
00093 #ifdef MPI_PARALLEL
00094 int ierr,mAddress,mIndex,mCount;
00095 #endif
00096
00097
00098 nDim=1;
00099 for (i=1; i<3; i++) if (pM->Nx[i]>1) nDim++;
00100
00101
00102
00103 for (nl=(pM->NLevels)-1; nl>=0; nl--){
00104
00105 #ifdef MPI_PARALLEL
00106
00107
00108
00109
00110 if (nl>0) {
00111 for (nd=0; nd<(pM->DomainsPerLevel[nl-1]); nd++){
00112 if (pM->Domain[nl-1][nd].Grid != NULL) {
00113 pG=pM->Domain[nl-1][nd].Grid;
00114
00115
00116
00117
00118 mAddress = 0;
00119 rbufN = ((nl-1) % 2);
00120 for (ncg=(pG->NmyCGrid); ncg<(pG->NCGrid); ncg++){
00121 mIndex = ncg - pG->NmyCGrid;
00122 ierr = MPI_Irecv(&(recv_bufRC[rbufN][nd][mAddress]),
00123 pG->CGrid[ncg].nWordsRC, MPI_DOUBLE, pG->CGrid[ncg].ID,
00124 pG->CGrid[ncg].DomN, pM->Domain[nl-1][nd].Comm_Children,
00125 &(recv_rq[nl-1][nd][mIndex]));
00126 mAddress += pG->CGrid[ncg].nWordsRC;
00127 }
00128
00129 }
00130 }
00131 }
00132 #endif
00133
00134
00135
00136
00137
00138 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00139
00140 if (pM->Domain[nl][nd].Grid != NULL) {
00141 pG=pM->Domain[nl][nd].Grid;
00142 rbufN = (nl % 2);
00143
00144 for (ncg=0; ncg<(pG->NCGrid); ncg++){
00145
00146
00147
00148
00149
00150
00151 if (ncg < pG->NmyCGrid) {
00152 pCO=(GridOvrlpS*)&(pG->CGrid[ncg]);
00153 pRcv = (double*)&(send_bufRC[pCO->DomN][0]);
00154 } else {
00155
00156 #ifdef MPI_PARALLEL
00157
00158
00159
00160
00161 mCount = pG->NCGrid - pG->NmyCGrid;
00162 ierr = MPI_Waitany(mCount,recv_rq[nl][nd],&mIndex,MPI_STATUS_IGNORE);
00163 if(mIndex == MPI_UNDEFINED){
00164 ath_error("[RestCorr]: Invalid request index nl=%i nd=%i\n",nl,nd);
00165 }
00166
00167
00168 mAddress = 0;
00169 mIndex += pG->NmyCGrid;
00170 for (i=pG->NmyCGrid; i<mIndex; i++) mAddress += pG->CGrid[i].nWordsRC;
00171 pCO=(GridOvrlpS*)&(pG->CGrid[mIndex]);
00172 pRcv = (double*)&(recv_bufRC[rbufN][nd][mAddress]);
00173 #else
00174
00175
00176 ath_error("[RestCorr]: no Child grid on Domain[%d][%d]\n",nl,nd);
00177 #endif
00178 }
00179
00180
00181
00182 ics = pCO->ijks[0];
00183 ice = pCO->ijke[0];
00184 jcs = pCO->ijks[1];
00185 jce = pCO->ijke[1];
00186 kcs = pCO->ijks[2];
00187 kce = pCO->ijke[2];
00188
00189
00190
00191 for (k=kcs; k<=kce; k++) {
00192 for (j=jcs; j<=jce; j++) {
00193 for (i=ics; i<=ice; i++) {
00194 pG->U[k][j][i].d = *(pRcv++);
00195 pG->U[k][j][i].M1 = *(pRcv++);
00196 pG->U[k][j][i].M2 = *(pRcv++);
00197 pG->U[k][j][i].M3 = *(pRcv++);
00198 #ifndef BAROTROPIC
00199 pG->U[k][j][i].E = *(pRcv++);
00200 #endif
00201 #ifdef MHD
00202 pG->U[k][j][i].B1c = *(pRcv++);
00203 pG->U[k][j][i].B2c = *(pRcv++);
00204 pG->U[k][j][i].B3c = *(pRcv++);
00205 #endif
00206 #if (NSCALARS > 0)
00207 for (n=0; n<NSCALARS; n++) pG->U[k][j][i].s[n] = *(pRcv++);
00208 #endif
00209 }
00210 }
00211 }
00212
00213 #ifdef MHD
00214
00215
00216
00217
00218
00219 if (nDim == 1){
00220 for (i=ics; i<=ice; i++) {
00221 pG->B1i[kcs][jcs][i] = pG->U[kcs][jcs][i].B1c;
00222 pG->B2i[kcs][jcs][i] = pG->U[kcs][jcs][i].B2c;
00223 pG->B3i[kcs][jcs][i] = pG->U[kcs][jcs][i].B3c;
00224 }
00225 } else {
00226
00227 if (nDim == 2){
00228
00229 for (j=jcs ; j<=jce; j++) {
00230 for (i=ics+1; i<=ice; i++) {
00231
00232
00233 if (i==ics+1){
00234 if (pCO->myFlx[0] == NULL) {pG->B1i[kcs][j][ics] = *(pRcv++);}
00235 else {pRcv++;}
00236 }
00237
00238 pG->B1i[kcs][j][i] = *(pRcv++);
00239
00240
00241
00242 if (i==ice){
00243 if (pCO->myFlx[1] == NULL) {pG->B1i[kcs][j][ice+1] = *(pRcv++);}
00244 else {pRcv++;}
00245 }
00246 }}
00247
00248
00249
00250
00251 for (i=ics; i<=ice; i++) {
00252 if (pCO->myFlx[2] == NULL) {pG->B2i[kcs][jcs][i] = *(pRcv++);}
00253 else {pRcv++;}
00254 }
00255
00256 for (j=jcs+1; j<=jce; j++) {
00257 for (i=ics ; i<=ice; i++) {
00258 pG->B2i[kcs][j][i] = *(pRcv++);
00259 }}
00260
00261
00262
00263 for (i=ics; i<=ice; i++) {
00264 if (pCO->myFlx[3] == NULL) {pG->B2i[kcs][jce+1][i] = *(pRcv++);}
00265 else {pRcv++;}
00266 }
00267
00268
00269 for (j=jcs; j<=jce; j++) {
00270 for (i=ics; i<=ice; i++) {
00271 pG->U[kcs][j][i].B1c=0.5*(pG->B1i[kcs][j][i] +pG->B1i[kcs][j][i+1]);
00272 pG->U[kcs][j][i].B2c=0.5*(pG->B2i[kcs][j][i] +pG->B2i[kcs][j+1][i]);
00273 pG->B3i[kcs][j][i] = pG->U[kcs][j][i].B3c;
00274 }}
00275
00276 } else {
00277
00278 for (k=kcs ; k<=kce; k++) {
00279 for (j=jcs ; j<=jce; j++) {
00280 for (i=ics+1; i<=ice; i++) {
00281
00282
00283 if (i==ics+1){
00284 if (pCO->myFlx[0] == NULL) {pG->B1i[k][j][ics] = *(pRcv++);}
00285 else {pRcv++;}
00286 }
00287
00288 pG->B1i[k][j][i] = *(pRcv++);
00289
00290
00291
00292 if (i==ice){
00293 if (pCO->myFlx[1] == NULL) {pG->B1i[k][j][ice+1] = *(pRcv++);}
00294 else {pRcv++;}
00295 }
00296
00297 }}}
00298
00299
00300 for (k=kcs ; k<=kce; k++) {
00301
00302
00303 for (i=ics; i<=ice; i++) {
00304 if (pCO->myFlx[2] == NULL) {pG->B2i[k][jcs][i] = *(pRcv++);}
00305 else {pRcv++;}
00306 }
00307
00308 for (j=jcs+1; j<=jce; j++) {
00309 for (i=ics ; i<=ice; i++) {
00310 pG->B2i[k][j][i] = *(pRcv++);
00311 }}
00312
00313
00314
00315 for (i=ics; i<=ice; i++) {
00316 if (pCO->myFlx[3] == NULL) {pG->B2i[k][jce+1][i] = *(pRcv++);}
00317 else {pRcv++;}
00318 }
00319 }
00320
00321
00322
00323
00324 for (j=jcs; j<=jce; j++) {
00325 for (i=ics; i<=ice; i++) {
00326 if (pCO->myFlx[4] == NULL) {pG->B3i[kcs][j][i] = *(pRcv++);}
00327 else {pRcv++;}
00328 }}
00329
00330 for (k=kcs+1; k<=kce; k++) {
00331 for (j=jcs ; j<=jce; j++) {
00332 for (i=ics ; i<=ice; i++) {
00333 pG->B3i[k][j][i] = *(pRcv++);
00334 }}}
00335
00336
00337
00338 for (j=jcs; j<=jce; j++) {
00339 for (i=ics; i<=ice; i++) {
00340 if (pCO->myFlx[5] == NULL) {pG->B3i[kce+1][j][i] = *(pRcv++);}
00341 else {pRcv++;}
00342 }}
00343
00344
00345 for (k=kcs; k<=kce; k++) {
00346 for (j=jcs; j<=jce; j++) {
00347 for (i=ics; i<=ice; i++) {
00348 pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i] + pG->B1i[k][j][i+1]);
00349 pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
00350 pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
00351 }}}
00352
00353 }
00354
00355 }
00356 #endif
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 for (dim=0; dim<2; dim++){
00367 if (pCO->myFlx[dim] != NULL) {
00368 if (dim == 0) {i=ics-1; q1=-(pG->dt/pG->dx1);}
00369 if (dim == 1) {i=ice+1; q1= (pG->dt/pG->dx1);}
00370 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00371 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00372 pG->U[k][j][i].d -= q1*(pCO->myFlx[dim][kk][jj].d - *(pRcv++));
00373 pG->U[k][j][i].M1 -= q1*(pCO->myFlx[dim][kk][jj].M1 - *(pRcv++));
00374 pG->U[k][j][i].M2 -= q1*(pCO->myFlx[dim][kk][jj].M2 - *(pRcv++));
00375 pG->U[k][j][i].M3 -= q1*(pCO->myFlx[dim][kk][jj].M3 - *(pRcv++));
00376 #ifndef BAROTROPIC
00377 pG->U[k][j][i].E -= q1*(pCO->myFlx[dim][kk][jj].E - *(pRcv++));
00378 #endif
00379 #ifdef MHD
00380 pG->U[k][j][i].B1c -= q1*(pCO->myFlx[dim][kk][jj].B1c - *(pRcv++));
00381 pG->U[k][j][i].B2c -= q1*(pCO->myFlx[dim][kk][jj].B2c - *(pRcv++));
00382 pG->U[k][j][i].B3c -= q1*(pCO->myFlx[dim][kk][jj].B3c - *(pRcv++));
00383 #endif
00384 #if (NSCALARS > 0)
00385 for (n=0; n<NSCALARS; n++) pG->U[k][j][i].s[n] -=
00386 q1*(pCO->myFlx[dim][kk][jj].s[n] - *(pRcv++));
00387 #endif
00388 }}
00389 }
00390 }
00391
00392
00393
00394 for (dim=2; dim<4; dim++){
00395 if (pCO->myFlx[dim] != NULL) {
00396 if (dim == 2) {j=jcs-1; q2=-(pG->dt/pG->dx2);}
00397 if (dim == 3) {j=jce+1; q2= (pG->dt/pG->dx2);}
00398 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00399 for (i=ics, ii=0; i<=ice; i++, ii++) {
00400 pG->U[k][j][i].d -= q2*(pCO->myFlx[dim][kk][ii].d - *(pRcv++));
00401 pG->U[k][j][i].M1 -= q2*(pCO->myFlx[dim][kk][ii].M1 - *(pRcv++));
00402 pG->U[k][j][i].M2 -= q2*(pCO->myFlx[dim][kk][ii].M2 - *(pRcv++));
00403 pG->U[k][j][i].M3 -= q2*(pCO->myFlx[dim][kk][ii].M3 - *(pRcv++));
00404 #ifndef BAROTROPIC
00405 pG->U[k][j][i].E -= q2*(pCO->myFlx[dim][kk][ii].E - *(pRcv++));
00406 #endif
00407 #ifdef MHD
00408 pG->U[k][j][i].B1c -= q2*(pCO->myFlx[dim][kk][ii].B1c - *(pRcv++));
00409 pG->U[k][j][i].B2c -= q2*(pCO->myFlx[dim][kk][ii].B2c - *(pRcv++));
00410 pG->U[k][j][i].B3c -= q2*(pCO->myFlx[dim][kk][ii].B3c - *(pRcv++));
00411 #endif
00412 #if (NSCALARS > 0)
00413 for (n=0; n<NSCALARS; n++) pG->U[k][j][i].s[n] -=
00414 q2*(pCO->myFlx[dim][kk][ii].s[n] - *(pRcv++));
00415 #endif
00416 }}
00417 }
00418 }
00419
00420
00421
00422 for (dim=4; dim<6; dim++){
00423 if (pCO->myFlx[dim] != NULL) {
00424 if (dim == 4) {k=kcs-1; q3=-(pG->dt/pG->dx3);}
00425 if (dim == 5) {k=kce+1; q3= (pG->dt/pG->dx3);}
00426 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00427 for (i=ics, ii=0; i<=ice; i++, ii++) {
00428 pG->U[k][j][i].d -= q3*(pCO->myFlx[dim][jj][ii].d - *(pRcv++));
00429 pG->U[k][j][i].M1 -= q3*(pCO->myFlx[dim][jj][ii].M1 - *(pRcv++));
00430 pG->U[k][j][i].M2 -= q3*(pCO->myFlx[dim][jj][ii].M2 - *(pRcv++));
00431 pG->U[k][j][i].M3 -= q3*(pCO->myFlx[dim][jj][ii].M3 - *(pRcv++));
00432 #ifndef BAROTROPIC
00433 pG->U[k][j][i].E -= q3*(pCO->myFlx[dim][jj][ii].E - *(pRcv++));
00434 #endif
00435 #ifdef MHD
00436 pG->U[k][j][i].B1c -= q3*(pCO->myFlx[dim][jj][ii].B1c - *(pRcv++));
00437 pG->U[k][j][i].B2c -= q3*(pCO->myFlx[dim][jj][ii].B2c - *(pRcv++));
00438 pG->U[k][j][i].B3c -= q3*(pCO->myFlx[dim][jj][ii].B3c - *(pRcv++));
00439 #endif
00440 #if (NSCALARS > 0)
00441 for (n=0; n<NSCALARS; n++) pG->U[k][j][i].s[n] -=
00442 q3*(pCO->myFlx[dim][jj][ii].s[n] - *(pRcv++));
00443 #endif
00444 }}
00445 }
00446 }
00447
00448
00449
00450
00451
00452
00453 #ifdef MHD
00454 if (nDim > 1) {
00455
00456
00457
00458 for (dim=0; dim<2; dim++){
00459 if (pCO->myEMF3[dim] != NULL) {
00460 if (dim == 0) {
00461 i=ics-1; ib = ics;
00462 q1 = -(pG->dt/pG->dx1);
00463 q2 = -(pG->dt/pG->dx2);
00464 q3 = -(pG->dt/pG->dx3);
00465 }
00466 if (dim == 1) {
00467 i=ice+1; ib = ice+1;
00468 q1 = (pG->dt/pG->dx1);
00469 q2 = -(pG->dt/pG->dx2);
00470 q3 = -(pG->dt/pG->dx3);
00471 }
00472
00473 for (k=kcs, kk=0; k<=kce ; k++, kk++) {
00474 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00475 SMRemf3[kk][jj] = *(pRcv++);
00476 pG->B2i[k][j][i]+= q1*(pCO->myEMF3[dim][kk][jj] -SMRemf3[kk][jj]);
00477 }
00478
00479
00480
00481
00482 jj = jce-jcs+1;
00483 SMRemf3[kk][jj] = *(pRcv++);
00484 if((pCO->myEMF3[3] != NULL) || (jce==pG->je)){
00485 pG->B2i[k][jce+1][i] +=
00486 q1*(pCO->myEMF3[dim][kk][jj] - SMRemf3[kk][jj]);
00487 }
00488 }
00489
00490 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00491 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00492 pG->B1i[k][j][ib] -=
00493 q2*(pCO->myEMF3[dim][kk][jj+1] - SMRemf3[kk][jj+1]) -
00494 q2*(pCO->myEMF3[dim][kk][jj ] - SMRemf3[kk][jj ]);
00495 }}
00496
00497 if (nDim == 3){
00498
00499 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00500 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00501 SMRemf2[kk][jj] = *(pRcv++);
00502 pG->B3i[k][j][i] -= q1*(pCO->myEMF2[dim][kk][jj]-SMRemf2[kk][jj]);
00503 }}
00504
00505
00506
00507
00508 kk = kce-kcs+1;
00509 for (j=jcs, jj=0; j<=jce ; j++, jj++) {
00510 SMRemf2[kk][jj] = *(pRcv++);
00511 if((pCO->myEMF2[5] != NULL) || (kce==pG->ke)){
00512 pG->B3i[kce+1][j][i] -=
00513 q1*(pCO->myEMF2[dim][kk][jj] - SMRemf2[kk][jj]);
00514 }
00515 }
00516
00517 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00518 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00519 pG->B1i[k][j][ib] +=
00520 q3*(pCO->myEMF2[dim][kk+1][jj] - SMRemf2[kk+1][jj]) -
00521 q3*(pCO->myEMF2[dim][kk ][jj] - SMRemf2[kk ][jj]);
00522 }}
00523
00524 for (k=kcs-1; k<=kce+1; k++) {
00525 for (j=jcs; j<=jce; j++) {
00526 pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
00527 }}
00528
00529 }
00530
00531 for (k=kcs; k<=kce; k++) {
00532 for (j=jcs; j<=jce; j++) {
00533 pG->U[k][j][ib ].B1c=0.5*(pG->B1i[k][j][ib]+pG->B1i[k][j][ib+1]);
00534 pG->U[k][j][ib-1].B1c=0.5*(pG->B1i[k][j][ib-1]+pG->B1i[k][j][ib]);
00535 }
00536 for (j=jcs-1; j<=jce+1; j++) {
00537 pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
00538 }
00539 }
00540 }
00541 }
00542
00543
00544
00545 for (dim=2; dim<4; dim++){
00546 if (pCO->myEMF3[dim] != NULL) {
00547 if (dim == 2) {
00548 j=jcs-1; jb=jcs;
00549 q1 = -(pG->dt/pG->dx1);
00550 q2 = -(pG->dt/pG->dx2);
00551 q3 = -(pG->dt/pG->dx3);
00552 }
00553 if (dim == 3) {
00554 j=jce+1; jb=jce+1;
00555 q1 = -(pG->dt/pG->dx1);
00556 q2 = (pG->dt/pG->dx2);
00557 q3 = -(pG->dt/pG->dx3);
00558 }
00559
00560 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00561 for (i=ics, ii=0; i<=ice; i++, ii++) {
00562 SMRemf3[kk][ii] = *(pRcv++);
00563 pG->B1i[k][j][i]-= q2*(pCO->myEMF3[dim][kk][ii] -SMRemf3[kk][ii]);
00564 }
00565
00566
00567
00568
00569 ii = ice-ics+1;
00570 SMRemf3[kk][ii] = *(pRcv++);
00571 if ((pCO->myEMF3[1] != NULL) || (ice==pG->ie)){
00572 pG->B1i[k][j][ice+1] -=
00573 q2*(pCO->myEMF3[dim][kk][ii] - SMRemf3[kk][ii]);
00574 }
00575 }
00576 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00577 for (i=ics, ii=0; i<=ice; i++, ii++) {
00578 pG->B2i[k][jb][i] +=
00579 q1*(pCO->myEMF3[dim][kk][ii+1] - SMRemf3[kk][ii+1]) -
00580 q1*(pCO->myEMF3[dim][kk][ii ] - SMRemf3[kk][ii ]);
00581 }}
00582
00583 if (nDim == 3){
00584
00585 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00586 for (i=ics, ii=0; i<=ice; i++, ii++) {
00587 SMRemf1[kk][ii] = *(pRcv++);
00588 pG->B3i[k][j][i] += q2*(pCO->myEMF1[dim][kk][ii]-SMRemf1[kk][ii]);
00589 }}
00590
00591
00592
00593
00594 kk = kce-kcs+1;
00595 for (i=ics, ii=0; i<=ice ; i++, ii++) {
00596 SMRemf1[kk][ii] = *(pRcv++);
00597 if((pCO->myEMF1[5] != NULL) || (kce==pG->ke)){
00598 pG->B3i[kce+1][j][i] +=
00599 q2*(pCO->myEMF1[dim][kk][ii] - SMRemf1[kk][ii]);
00600 }
00601 }
00602
00603 for (k=kcs, kk=0; k<=kce; k++, kk++) {
00604 for (i=ics, ii=0; i<=ice; i++, ii++) {
00605 pG->B2i[k][jb][i] -=
00606 q3*(pCO->myEMF1[dim][kk+1][ii] - SMRemf1[kk+1][ii]) -
00607 q3*(pCO->myEMF1[dim][kk ][ii] - SMRemf1[kk ][ii]);
00608 }}
00609
00610 for (k=kcs-1; k<=kce+1; k++) {
00611 for (i=ics; i<=ice; i++) {
00612 pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
00613 }}
00614
00615 }
00616
00617 for (k=kcs; k<=kce; k++) {
00618 for (i=ics-1; i<=ice+1; i++) {
00619 pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i] + pG->B1i[k][j][i+1]);
00620 }
00621 for (i=ics; i<=ice; i++) {
00622 pG->U[k][jb ][i].B2c=0.5*(pG->B2i[k][jb][i]+pG->B2i[k][jb+1][i]);
00623 pG->U[k][jb-1][i].B2c=0.5*(pG->B2i[k][jb-1][i]+pG->B2i[k][jb][i]);
00624 }
00625 }
00626 }
00627 }
00628
00629
00630
00631 for (dim=4; dim<6; dim++){
00632 if (pCO->myEMF1[dim] != NULL) {
00633 if (dim == 4) {
00634 k=kcs-1; kb=kcs;
00635 q1 = -(pG->dt/pG->dx1);
00636 q2 = -(pG->dt/pG->dx2);
00637 q3 = -(pG->dt/pG->dx3);
00638 }
00639 if (dim == 5) {
00640 k=kce+1; kb=kce+1;
00641 q1 = -(pG->dt/pG->dx1);
00642 q2 = -(pG->dt/pG->dx2);
00643 q3 = (pG->dt/pG->dx3);
00644 }
00645
00646 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00647 for (i=ics, ii=0; i<=ice; i++, ii++) {
00648 SMRemf1[jj][ii] = *(pRcv++);
00649 pG->B2i[k][j][i] -= q3*(pCO->myEMF1[dim][jj][ii] - SMRemf1[jj][ii]);
00650 }}
00651
00652
00653
00654
00655 jj = jce-jcs+1;
00656 for (i=ics, ii=0; i<=ice ; i++, ii++) {
00657 SMRemf1[jj][ii] = *(pRcv++);
00658 if((pCO->myEMF1[3] != NULL) || (jce==pG->je)){
00659 pG->B2i[k][jce+1][i] -=
00660 q3*(pCO->myEMF1[dim][jj][ii] - SMRemf1[jj][ii]);
00661 }
00662 }
00663
00664 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00665 for (i=ics, ii=0; i<=ice; i++, ii++) {
00666 SMRemf2[jj][ii] = *(pRcv++);
00667 pG->B1i[k][j][i]+= q3*(pCO->myEMF2[dim][jj][ii] -SMRemf2[jj][ii]);
00668 }
00669
00670
00671
00672
00673 ii = ice-ics+1;
00674 SMRemf2[jj][ii] = *(pRcv++);
00675 if ((pCO->myEMF2[1] != NULL) || (ice==pG->ie)){
00676 pG->B1i[k][j][ice+1] +=
00677 q3*(pCO->myEMF2[dim][jj][ii] -SMRemf2[jj][ii]);
00678 }
00679 }
00680
00681 for (j=jcs, jj=0; j<=jce; j++, jj++) {
00682 for (i=ics, ii=0; i<=ice; i++, ii++) {
00683 pG->B3i[kb][j][i] +=
00684 q2*(pCO->myEMF1[dim][jj+1][ii ] - SMRemf1[jj+1][ii ]) -
00685 q2*(pCO->myEMF1[dim][jj ][ii ] - SMRemf1[jj ][ii ]) -
00686 q1*(pCO->myEMF2[dim][jj ][ii+1] - SMRemf2[jj ][ii+1]) +
00687 q1*(pCO->myEMF2[dim][jj ][ii ] - SMRemf2[jj ][ii ]);
00688 }}
00689
00690 for (j=jcs; j<=jce; j++) {
00691 for (i=ics-1; i<=ice+1; i++) {
00692 pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i] + pG->B1i[k][j][i+1]);
00693 }}
00694
00695 for (j=jcs-1; j<=jce+1; j++) {
00696 for (i=ics; i<=ice; i++) {
00697 pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
00698 }}
00699
00700 for (j=jcs; j<=jce; j++) {
00701 for (i=ics; i<=ice; i++) {
00702 pG->U[kb ][j][i].B3c=0.5*(pG->B3i[kb][j][i] + pG->B3i[kb+1][j][i]);
00703 pG->U[kb-1][j][i].B3c=0.5*(pG->B3i[kb-1][j][i] + pG->B3i[kb][j][i]);
00704 }}
00705 }
00706 }
00707
00708 }
00709 #endif
00710
00711 }
00712 }}
00713
00714
00715
00716
00717
00718
00719
00720 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00721
00722 if (pM->Domain[nl][nd].Grid != NULL) {
00723 pG=pM->Domain[nl][nd].Grid;
00724 start_addr=0;
00725
00726 for (npg=0; npg<(pG->NPGrid); npg++){
00727 pPO=(GridOvrlpS*)&(pG->PGrid[npg]);
00728 cnt = 0;
00729
00730
00731
00732 ips = pPO->ijks[0];
00733 ipe = pPO->ijke[0];
00734 jps = pPO->ijks[1];
00735 jpe = pPO->ijke[1];
00736 kps = pPO->ijks[2];
00737 kpe = pPO->ijke[2];
00738
00739
00740
00741
00742 pSnd = (double*)&(send_bufRC[nd][start_addr]);
00743 for (k=kps; k<=kpe; k+=2) {
00744 for (j=jps; j<=jpe; j+=2) {
00745 for (i=ips; i<=ipe; i+=2) {
00746 *(pSnd++) = pG->U[k][j][i].d + pG->U[k][j][i+1].d;
00747 *(pSnd++) = pG->U[k][j][i].M1 + pG->U[k][j][i+1].M1;
00748 *(pSnd++) = pG->U[k][j][i].M2 + pG->U[k][j][i+1].M2;
00749 *(pSnd++) = pG->U[k][j][i].M3 + pG->U[k][j][i+1].M3;
00750 #ifndef BAROTROPIC
00751 *(pSnd++) = pG->U[k][j][i].E + pG->U[k][j][i+1].E;
00752 #endif
00753 #ifdef MHD
00754 *(pSnd++) = pG->U[k][j][i].B1c + pG->U[k][j][i+1].B1c;
00755 *(pSnd++) = pG->U[k][j][i].B2c + pG->U[k][j][i+1].B2c;
00756 *(pSnd++) = pG->U[k][j][i].B3c + pG->U[k][j][i+1].B3c;
00757 #endif
00758 #if (NSCALARS > 0)
00759 for (n=0; n<NSCALARS; n++)
00760 *(pSnd++) = pG->U[k][j][i].s[n] + pG->U[k][j][i+1].s[n];
00761 #endif
00762 }}}
00763 fact = 0.5;
00764 nCons = (ipe-ips+1)*(NVAR)/2;
00765
00766
00767
00768 if (pG->Nx[1] > 1) {
00769 pSnd = (double*)&(send_bufRC[nd][start_addr]);
00770 for (k=kps; k<=kpe; k+=2) {
00771 for (j=jps; j<=jpe; j+=2) {
00772 for (i=ips; i<=ipe; i+=2) {
00773 *(pSnd++) += pG->U[k][j+1][i].d + pG->U[k][j+1][i+1].d;
00774 *(pSnd++) += pG->U[k][j+1][i].M1 + pG->U[k][j+1][i+1].M1;
00775 *(pSnd++) += pG->U[k][j+1][i].M2 + pG->U[k][j+1][i+1].M2;
00776 *(pSnd++) += pG->U[k][j+1][i].M3 + pG->U[k][j+1][i+1].M3;
00777 #ifndef BAROTROPIC
00778 *(pSnd++) += pG->U[k][j+1][i].E + pG->U[k][j+1][i+1].E;
00779 #endif
00780 #ifdef MHD
00781 *(pSnd++) += pG->U[k][j+1][i].B1c + pG->U[k][j+1][i+1].B1c;
00782 *(pSnd++) += pG->U[k][j+1][i].B2c + pG->U[k][j+1][i+1].B2c;
00783 *(pSnd++) += pG->U[k][j+1][i].B3c + pG->U[k][j+1][i+1].B3c;
00784 #endif
00785 #if (NSCALARS > 0)
00786 for (n=0; n<NSCALARS; n++)
00787 *(pSnd++) += pG->U[k][j+1][i].s[n] + pG->U[k][j+1][i+1].s[n];
00788 #endif
00789 }}}
00790 fact = 0.25;
00791 nCons = (jpe-jps+1)*(ipe-ips+1)*(NVAR)/4;
00792 }
00793
00794
00795
00796 if (pG->Nx[2] > 1) {
00797 pSnd = (double*)&(send_bufRC[nd][start_addr]);
00798 for (k=kps; k<=kpe; k+=2) {
00799 for (j=jps; j<=jpe; j+=2) {
00800 for (i=ips; i<=ipe; i+=2) {
00801 *(pSnd++) += pG->U[k+1][j ][i].d + pG->U[k+1][j ][i+1].d +
00802 pG->U[k+1][j+1][i].d + pG->U[k+1][j+1][i+1].d;
00803 *(pSnd++) += pG->U[k+1][j ][i].M1 + pG->U[k+1][j ][i+1].M1 +
00804 pG->U[k+1][j+1][i].M1 + pG->U[k+1][j+1][i+1].M1;
00805 *(pSnd++) += pG->U[k+1][j ][i].M2 + pG->U[k+1][j ][i+1].M2 +
00806 pG->U[k+1][j+1][i].M2 + pG->U[k+1][j+1][i+1].M2;
00807 *(pSnd++) += pG->U[k+1][j ][i].M3 + pG->U[k+1][j ][i+1].M3 +
00808 pG->U[k+1][j+1][i].M3 + pG->U[k+1][j+1][i+1].M3;
00809 #ifndef BAROTROPIC
00810 *(pSnd++) += pG->U[k+1][j ][i].E + pG->U[k+1][j ][i+1].E +
00811 pG->U[k+1][j+1][i].E + pG->U[k+1][j+1][i+1].E;
00812 #endif
00813 #ifdef MHD
00814 *(pSnd++) += pG->U[k+1][j ][i].B1c + pG->U[k+1][j ][i+1].B1c +
00815 pG->U[k+1][j+1][i].B1c + pG->U[k+1][j+1][i+1].B1c;
00816 *(pSnd++) += pG->U[k+1][j ][i].B2c + pG->U[k+1][j ][i+1].B2c +
00817 pG->U[k+1][j+1][i].B2c + pG->U[k+1][j+1][i+1].B2c;
00818 *(pSnd++) += pG->U[k+1][j ][i].B3c + pG->U[k+1][j ][i+1].B3c +
00819 pG->U[k+1][j+1][i].B3c + pG->U[k+1][j+1][i+1].B3c;
00820 #endif
00821 #if (NSCALARS > 0)
00822 for (n=0; n<NSCALARS; n++)
00823 *(pSnd++) += pG->U[k+1][j ][i].s[n] + pG->U[k+1][j ][i+1].s[n] +
00824 pG->U[k+1][j+1][i].s[n] + pG->U[k+1][j+1][i+1].s[n];
00825 #endif
00826 }}}
00827 fact = 0.125;
00828 nCons = (kpe-kps+1)*(jpe-jps+1)*(ipe-ips+1)*(NVAR)/8;
00829 }
00830
00831
00832 pSnd = (double*)&(send_bufRC[nd][start_addr]);
00833 for (i=start_addr; i<(start_addr+nCons); i++) *(pSnd++) *= fact;
00834 cnt = nCons;
00835
00836 #ifdef MHD
00837
00838
00839
00840
00841
00842 if (nDim == 3) {
00843
00844 for (k=kps; k<=kpe ; k+=2) {
00845 for (j=jps; j<=jpe ; j+=2) {
00846 for (i=ips; i<=ipe+1; i+=2) {
00847 *(pSnd++) = 0.25*(pG->B1i[k ][j][i] + pG->B1i[k ][j+1][i]
00848 + pG->B1i[k+1][j][i] + pG->B1i[k+1][j+1][i]);
00849 }}}
00850 for (k=kps; k<=kpe ; k+=2) {
00851 for (j=jps; j<=jpe+1; j+=2) {
00852 for (i=ips; i<=ipe ; i+=2) {
00853 *(pSnd++) = 0.25*(pG->B2i[k ][j][i] + pG->B2i[k ][j][i+1]
00854 + pG->B2i[k+1][j][i] + pG->B2i[k+1][j][i+1]);
00855 }}}
00856 for (k=kps; k<=kpe+1; k+=2) {
00857 for (j=jps; j<=jpe ; j+=2) {
00858 for (i=ips; i<=ipe ; i+=2) {
00859 *(pSnd++) = 0.25*(pG->B3i[k][j ][i] + pG->B3i[k][j ][i+1]
00860 + pG->B3i[k][j+1][i] + pG->B3i[k][j+1][i+1]);
00861 }}}
00862 nFld = ((kpe-kps+1)/2 )*((jpe-jps+1)/2 )*((ipe-ips+1)/2 + 1) +
00863 ((kpe-kps+1)/2 )*((jpe-jps+1)/2 + 1)*((ipe-ips+1)/2 ) +
00864 ((kpe-kps+1)/2 + 1)*((jpe-jps+1)/2 )*((ipe-ips+1)/2 );
00865
00866 } else {
00867
00868 if (nDim == 2) {
00869 for (j=jps; j<=jpe ; j+=2) {
00870 for (i=ips; i<=ipe+1; i+=2) {
00871 *(pSnd++) = 0.5*(pG->B1i[kps][j][i] + pG->B1i[kps][j+1][i]);
00872 }}
00873 for (j=jps; j<=jpe+1; j+=2) {
00874 for (i=ips; i<=ipe ; i+=2) {
00875 *(pSnd++) = 0.5*(pG->B2i[kps][j][i] + pG->B2i[kps][j][i+1]);
00876 }}
00877 nFld = ((jpe-jps+1)/2 )*((ipe-ips+1)/2 + 1) +
00878 ((jpe-jps+1)/2 + 1)*((ipe-ips+1)/2 );
00879 }
00880
00881 }
00882 cnt += nFld;
00883 #endif
00884
00885
00886
00887
00888 for (dim=0; dim<2; dim++){
00889 if (pPO->myFlx[dim] != NULL) {
00890
00891 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
00892
00893 if (nDim == 1) {
00894
00895 *(pSnd++) = pPO->myFlx[dim][kps][jps].d ;
00896 *(pSnd++) = pPO->myFlx[dim][kps][jps].M1;
00897 *(pSnd++) = pPO->myFlx[dim][kps][jps].M2;
00898 *(pSnd++) = pPO->myFlx[dim][kps][jps].M3;
00899 #ifndef BAROTROPIC
00900 *(pSnd++) = pPO->myFlx[dim][kps][jps].E ;
00901 #endif
00902 #ifdef MHD
00903 *(pSnd++) = pPO->myFlx[dim][kps][jps].B1c;
00904 *(pSnd++) = pPO->myFlx[dim][kps][jps].B2c;
00905 *(pSnd++) = pPO->myFlx[dim][kps][jps].B3c;
00906 #endif
00907 #if (NSCALARS > 0)
00908 for (n=0; n<NSCALARS; n++) *(pSnd++) = pPO->myFlx[dim][kps][jps].s[n];
00909 #endif
00910 nFlx = NVAR;
00911 } else {
00912
00913
00914
00915 for (k=0; k<=(kpe-kps); k+=2) {
00916 for (j=0; j<=(jpe-jps); j+=2) {
00917 *(pSnd++) = pPO->myFlx[dim][k][j].d + pPO->myFlx[dim][k][j+1].d;
00918 *(pSnd++) = pPO->myFlx[dim][k][j].M1 + pPO->myFlx[dim][k][j+1].M1;
00919 *(pSnd++) = pPO->myFlx[dim][k][j].M2 + pPO->myFlx[dim][k][j+1].M2;
00920 *(pSnd++) = pPO->myFlx[dim][k][j].M3 + pPO->myFlx[dim][k][j+1].M3;
00921 #ifndef BAROTROPIC
00922 *(pSnd++) = pPO->myFlx[dim][k][j].E + pPO->myFlx[dim][k][j+1].E;
00923 #endif
00924 #ifdef MHD
00925 *(pSnd++) = pPO->myFlx[dim][k][j].B1c + pPO->myFlx[dim][k][j+1].B1c;
00926 *(pSnd++) = pPO->myFlx[dim][k][j].B2c + pPO->myFlx[dim][k][j+1].B2c;
00927 *(pSnd++) = pPO->myFlx[dim][k][j].B3c + pPO->myFlx[dim][k][j+1].B3c;
00928 #endif
00929 #if (NSCALARS > 0)
00930 for (n=0; n<NSCALARS; n++) *(pSnd++) =
00931 pPO->myFlx[dim][k][j].s[n] + pPO->myFlx[dim][k][j+1].s[n];
00932 #endif
00933 }}
00934 fact = 0.5;
00935 nFlx = ((jpe-jps+1)/2)*(NVAR);
00936
00937
00938
00939 if (nDim == 3) {
00940 pSnd = (double*)&(send_bufRC[nd][start_addr+cnt]);
00941 for (k=0; k<=(kpe-kps); k+=2) {
00942 for (j=0; j<=(jpe-jps); j+=2) {
00943 *(pSnd++)+=pPO->myFlx[dim][k+1][j].d +pPO->myFlx[dim][k+1][j+1].d;
00944 *(pSnd++)+=pPO->myFlx[dim][k+1][j].M1 +pPO->myFlx[dim][k+1][j+1].M1;
00945 *(pSnd++)+=pPO->myFlx[dim][k+1][j].M2 +pPO->myFlx[dim][k+1][j+1].M2;
00946 *(pSnd++)+=pPO->myFlx[dim][k+1][j].M3 +pPO->myFlx[dim][k+1][j+1].M3;
00947 #ifndef BAROTROPIC
00948 *(pSnd++)+=pPO->myFlx[dim][k+1][j].E +pPO->myFlx[dim][k+1][j+1].E;
00949 #endif
00950 #ifdef MHD
00951 *(pSnd++)+=pPO->myFlx[dim][k+1][j].B1c +pPO->myFlx[dim][k+1][j+1].B1c;
00952 *(pSnd++)+=pPO->myFlx[dim][k+1][j].B2c +pPO->myFlx[dim][k+1][j+1].B2c;
00953 *(pSnd++)+=pPO->myFlx[dim][k+1][j].B3c +pPO->myFlx[dim][k+1][j+1].B3c;
00954 #endif
00955 #if (NSCALARS > 0)
00956 for (n=0; n<NSCALARS; n++) *(pSnd++) +=
00957 pPO->myFlx[dim][k+1][j].s[n] + pPO->myFlx[dim][k+1][j+1].s[n];
00958 #endif
00959 }}
00960 fact = 0.25;
00961 nFlx = ((kpe-kps+1)*(jpe-jps+1)/4)*(NVAR);
00962 }
00963
00964
00965 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
00966 for (i=(start_addr+cnt); i<(start_addr+cnt+nFlx); i++) *(pSnd++) *=fact;
00967 }
00968 cnt += nFlx;
00969
00970 }}
00971
00972
00973
00974 for (dim=2; dim<4; dim++){
00975 if (pPO->myFlx[dim] != NULL) {
00976
00977 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
00978
00979
00980
00981 for (k=0; k<=(kpe-kps); k+=2) {
00982 for (i=0; i<=(ipe-ips); i+=2) {
00983 *(pSnd++) = pPO->myFlx[dim][k][i].d + pPO->myFlx[dim][k][i+1].d;
00984 *(pSnd++) = pPO->myFlx[dim][k][i].M1 + pPO->myFlx[dim][k][i+1].M1;
00985 *(pSnd++) = pPO->myFlx[dim][k][i].M2 + pPO->myFlx[dim][k][i+1].M2;
00986 *(pSnd++) = pPO->myFlx[dim][k][i].M3 + pPO->myFlx[dim][k][i+1].M3;
00987 #ifndef BAROTROPIC
00988 *(pSnd++) = pPO->myFlx[dim][k][i].E + pPO->myFlx[dim][k][i+1].E;
00989 #endif
00990 #ifdef MHD
00991 *(pSnd++) = pPO->myFlx[dim][k][i].B1c + pPO->myFlx[dim][k][i+1].B1c;
00992 *(pSnd++) = pPO->myFlx[dim][k][i].B2c + pPO->myFlx[dim][k][i+1].B2c;
00993 *(pSnd++) = pPO->myFlx[dim][k][i].B3c + pPO->myFlx[dim][k][i+1].B3c;
00994 #endif
00995 #if (NSCALARS > 0)
00996 for (n=0; n<NSCALARS; n++) *(pSnd++) =
00997 pPO->myFlx[dim][k][i].s[n] + pPO->myFlx[dim][k][i+1].s[n];
00998 #endif
00999 }}
01000 fact = 0.5;
01001 nFlx = ((ipe-ips+1)/2)*(NVAR);
01002
01003
01004
01005 if (nDim == 3) {
01006 pSnd = (double*)&(send_bufRC[nd][start_addr+cnt]);
01007 for (k=0; k<=(kpe-kps); k+=2) {
01008 for (i=0; i<=(ipe-ips); i+=2) {
01009 *(pSnd++) +=pPO->myFlx[dim][k+1][i].d + pPO->myFlx[dim][k+1][i+1].d;
01010 *(pSnd++) +=pPO->myFlx[dim][k+1][i].M1 + pPO->myFlx[dim][k+1][i+1].M1;
01011 *(pSnd++) +=pPO->myFlx[dim][k+1][i].M2 + pPO->myFlx[dim][k+1][i+1].M2;
01012 *(pSnd++) +=pPO->myFlx[dim][k+1][i].M3 + pPO->myFlx[dim][k+1][i+1].M3;
01013 #ifndef BAROTROPIC
01014 *(pSnd++) +=pPO->myFlx[dim][k+1][i].E + pPO->myFlx[dim][k+1][i+1].E;
01015 #endif
01016 #ifdef MHD
01017 *(pSnd++)+= pPO->myFlx[dim][k+1][i].B1c+pPO->myFlx[dim][k+1][i+1].B1c;
01018 *(pSnd++)+= pPO->myFlx[dim][k+1][i].B2c+pPO->myFlx[dim][k+1][i+1].B2c;
01019 *(pSnd++)+= pPO->myFlx[dim][k+1][i].B3c+pPO->myFlx[dim][k+1][i+1].B3c;
01020 #endif
01021 #if (NSCALARS > 0)
01022 for (n=0; n<NSCALARS; n++) *(pSnd++) +=
01023 pPO->myFlx[dim][k+1][i].s[n] + pPO->myFlx[dim][k+1][i+1].s[n];
01024 #endif
01025 }}
01026 fact = 0.25;
01027 nFlx = ((kpe-kps+1)*(ipe-ips+1)/4)*(NVAR);
01028 }
01029
01030
01031 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
01032 for (i=(start_addr+cnt); i<(start_addr+cnt+nFlx); i++) *(pSnd++) *= fact;
01033 cnt += nFlx;
01034
01035 }}
01036
01037
01038
01039 for (dim=4; dim<6; dim++){
01040 if (pPO->myFlx[dim] != NULL) {
01041
01042 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
01043
01044
01045
01046 for (j=0; j<=(jpe-jps); j+=2) {
01047 for (i=0; i<=(ipe-ips); i+=2) {
01048 *(pSnd++) = pPO->myFlx[dim][j][i].d + pPO->myFlx[dim][j][i+1].d;
01049 *(pSnd++) = pPO->myFlx[dim][j][i].M1 + pPO->myFlx[dim][j][i+1].M1;
01050 *(pSnd++) = pPO->myFlx[dim][j][i].M2 + pPO->myFlx[dim][j][i+1].M2;
01051 *(pSnd++) = pPO->myFlx[dim][j][i].M3 + pPO->myFlx[dim][j][i+1].M3;
01052 #ifndef BAROTROPIC
01053 *(pSnd++) = pPO->myFlx[dim][j][i].E + pPO->myFlx[dim][j][i+1].E;
01054 #endif
01055 #ifdef MHD
01056 *(pSnd++) = pPO->myFlx[dim][j][i].B1c + pPO->myFlx[dim][j][i+1].B1c;
01057 *(pSnd++) = pPO->myFlx[dim][j][i].B2c + pPO->myFlx[dim][j][i+1].B2c;
01058 *(pSnd++) = pPO->myFlx[dim][j][i].B3c + pPO->myFlx[dim][j][i+1].B3c;
01059 #endif
01060 #if (NSCALARS > 0)
01061 for (n=0; n<NSCALARS; n++) *(pSnd++) =
01062 pPO->myFlx[dim][j][i].s[n] + pPO->myFlx[dim][j][i+1].s[n];
01063 #endif
01064 }}
01065
01066
01067
01068 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
01069 for (j=0; j<=(jpe-jps); j+=2) {
01070 for (i=0; i<=(ipe-ips); i+=2) {
01071 *(pSnd++) +=pPO->myFlx[dim][j+1][i].d + pPO->myFlx[dim][j+1][i+1].d;
01072 *(pSnd++) +=pPO->myFlx[dim][j+1][i].M1 + pPO->myFlx[dim][j+1][i+1].M1;
01073 *(pSnd++) +=pPO->myFlx[dim][j+1][i].M2 + pPO->myFlx[dim][j+1][i+1].M2;
01074 *(pSnd++) +=pPO->myFlx[dim][j+1][i].M3 + pPO->myFlx[dim][j+1][i+1].M3;
01075 #ifndef BAROTROPIC
01076 *(pSnd++) +=pPO->myFlx[dim][j+1][i].E + pPO->myFlx[dim][j+1][i+1].E;
01077 #endif
01078 #ifdef MHD
01079 *(pSnd++) +=pPO->myFlx[dim][j+1][i].B1c + pPO->myFlx[dim][j+1][i+1].B1c;
01080 *(pSnd++) +=pPO->myFlx[dim][j+1][i].B2c + pPO->myFlx[dim][j+1][i+1].B2c;
01081 *(pSnd++) +=pPO->myFlx[dim][j+1][i].B3c + pPO->myFlx[dim][j+1][i+1].B3c;
01082 #endif
01083 #if (NSCALARS > 0)
01084 for (n=0; n<NSCALARS; n++) *(pSnd++) +=
01085 pPO->myFlx[dim][j+1][i].s[n] + pPO->myFlx[dim][j+1][i+1].s[n];
01086 #endif
01087 }}
01088 fact = 0.25;
01089 nFlx = ((jpe-jps+1)*(ipe-ips+1)/4)*(NVAR);
01090
01091
01092 pSnd = (double*)&(send_bufRC[nd][(start_addr+cnt)]);
01093 for (i=(start_addr+cnt); i<(start_addr+cnt+nFlx); i++) *(pSnd++) *= fact;
01094 cnt += nFlx;
01095
01096 }}
01097
01098
01099
01100
01101
01102
01103
01104 #ifdef MHD
01105 for (dim=0; dim<2; dim++){
01106 if (pPO->myEMF3[dim] != NULL) {
01107
01108
01109
01110 if (pG->Nx[2] == 1) {
01111 for (k=0; k<=(kpe-kps) ; k+=2) {
01112 for (j=0; j<=(jpe-jps)+1; j+=2) {
01113 *(pSnd++) = pPO->myEMF3[dim][k][j];
01114 }}
01115
01116 } else {
01117
01118
01119
01120 for (k=0; k<=(kpe-kps) ; k+=2) {
01121 for (j=0; j<=(jpe-jps)+1; j+=2) {
01122 *(pSnd++) = 0.5*(pPO->myEMF3[dim][k][j]+pPO->myEMF3[dim][k+1][j]);
01123 }}
01124
01125 for (k=0; k<=(kpe-kps)+1; k+=2) {
01126 for (j=0; j<=(jpe-jps) ; j+=2) {
01127 *(pSnd++) = 0.5*(pPO->myEMF2[dim][k][j]+pPO->myEMF2[dim][k][j+1]);
01128 }}
01129 }
01130 }
01131 }
01132
01133
01134
01135 for (dim=2; dim<4; dim++){
01136 if (pPO->myEMF3[dim] != NULL) {
01137
01138
01139
01140 if (pG->Nx[2] == 1) {
01141 for (k=0; k<=(kpe-kps) ; k+=2) {
01142 for (i=0; i<=(ipe-ips)+1; i+=2) {
01143 *(pSnd++) = pPO->myEMF3[dim][k][i];
01144 }}
01145
01146 } else {
01147
01148
01149
01150 for (k=0; k<=(kpe-kps) ; k+=2) {
01151 for (i=0; i<=(ipe-ips)+1; i+=2) {
01152 *(pSnd++) = 0.5*(pPO->myEMF3[dim][k][i]+pPO->myEMF3[dim][k+1][i]);
01153 }}
01154
01155 for (k=0; k<=(kpe-kps)+1; k+=2) {
01156 for (i=0; i<=(ipe-ips) ; i+=2) {
01157 *(pSnd++) = 0.5*(pPO->myEMF1[dim][k][i]+pPO->myEMF1[dim][k][i+1]);
01158 }}
01159 }
01160 }
01161 }
01162
01163
01164
01165
01166 for (dim=4; dim<6; dim++){
01167 if (pPO->myEMF1[dim] != NULL) {
01168
01169
01170
01171 for (j=0; j<=(jpe-jps)+1; j+=2) {
01172 for (i=0; i<=(ipe-ips) ; i+=2) {
01173 *(pSnd++) = 0.5*(pPO->myEMF1[dim][j][i] + pPO->myEMF1[dim][j][i+1]);
01174 }}
01175
01176 for (j=0; j<=(jpe-jps) ; j+=2) {
01177 for (i=0; i<=(ipe-ips)+1; i+=2) {
01178 *(pSnd++) = 0.5*(pPO->myEMF2[dim][j][i] + pPO->myEMF2[dim][j+1][i]);
01179 }}
01180 }
01181 }
01182
01183 #endif
01184
01185 #ifdef MPI_PARALLEL
01186
01187
01188
01189 if (npg >= pG->NmyPGrid){
01190 mIndex = npg - pG->NmyPGrid;
01191 ierr = MPI_Isend(&(send_bufRC[nd][start_addr]), pG->PGrid[npg].nWordsRC,
01192 MPI_DOUBLE, pG->PGrid[npg].ID, nd, pM->Domain[nl][nd].Comm_Parent,
01193 &(send_rq[nd][mIndex]));
01194 }
01195 #endif
01196
01197 start_addr += pG->PGrid[npg].nWordsRC;
01198
01199 }
01200 }}
01201
01202 #ifdef MPI_PARALLEL
01203
01204
01205
01206
01207 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01208 if (pM->Domain[nl][nd].Grid != NULL) {
01209 pG=pM->Domain[nl][nd].Grid;
01210
01211 if (pG->NPGrid > pG->NmyPGrid) {
01212 mCount = pG->NPGrid - pG->NmyPGrid;
01213 ierr = MPI_Waitall(mCount, send_rq[nd], MPI_STATUS_IGNORE);
01214 }
01215 }
01216 }
01217 #endif
01218
01219 }
01220 }
01221
01222
01223
01224
01225
01226
01227 void Prolongate(MeshS *pM)
01228 {
01229 GridS *pG;
01230 int nDim,nl,nd,ncg,dim,npg,rbufN,id,l,m,n,mend,nend;
01231 int i,ii,ics,ice,ips,ipe,igzs,igze;
01232 int j,jj,jcs,jce,jps,jpe,jgzs,jgze;
01233 int k,kk,kcs,kce,kps,kpe,kgzs,kgze;
01234 int ngz1,ngz2,ngz3;
01235 double *pRcv,*pSnd;
01236 GridOvrlpS *pCO, *pPO;
01237 ConsS ProlongedC[2][2][2];
01238 #if (NSCALARS > 0)
01239 int n;
01240 #endif
01241 #ifdef MHD
01242 Real3Vect BGZ[3][3][3], ProlongedF[3][3][3];
01243 #endif
01244 #ifdef MPI_PARALLEL
01245 int ierr,mAddress,mIndex,mCount;
01246 #endif
01247
01248
01249 nDim=1;
01250 for (dim=1; dim<3; dim++) if (pM->Nx[dim]>1) nDim++;
01251
01252
01253
01254 for (nl=0; nl<(pM->NLevels); nl++){
01255
01256 #ifdef MPI_PARALLEL
01257
01258
01259
01260
01261 if (nl<(pM->NLevels)-1) {
01262 for (nd=0; nd<(pM->DomainsPerLevel[nl+1]); nd++){
01263 if (pM->Domain[nl+1][nd].Grid != NULL) {
01264 pG=pM->Domain[nl+1][nd].Grid;
01265
01266
01267
01268
01269
01270
01271 mAddress = 0;
01272 rbufN = ((nl+1) % 2);
01273 if (pG->NmyPGrid > 0) mAddress = pG->PGrid[0].nWordsP;
01274
01275 for (npg=(pG->NmyPGrid); npg<(pG->NPGrid); npg++){
01276 mIndex = npg - pG->NmyPGrid;
01277 ierr = MPI_Irecv(&(recv_bufP[rbufN][nd][mAddress]),
01278 pG->PGrid[npg].nWordsP, MPI_DOUBLE, pG->PGrid[npg].ID,
01279 pG->PGrid[npg].DomN, pM->Domain[nl+1][nd].Comm_Parent,
01280 &(recv_rq[nl+1][nd][mIndex]));
01281 mAddress += pG->PGrid[npg].nWordsP;
01282 }
01283
01284 }
01285 }
01286 }
01287 #endif
01288
01289
01290
01291
01292 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01293
01294 if (pM->Domain[nl][nd].Grid != NULL) {
01295 pG=pM->Domain[nl][nd].Grid;
01296 for(i=0; i<maxND; i++) start_addrP[i] = 0;
01297
01298 for (ncg=0; ncg<(pG->NCGrid); ncg++){
01299 pCO=(GridOvrlpS*)&(pG->CGrid[ncg]);
01300
01301
01302
01303 pSnd = (double*)&(send_bufP[pCO->DomN][start_addrP[pCO->DomN]]);
01304
01305 for (dim=0; dim<(2*nDim); dim++){
01306 if (pCO->myFlx[dim] != NULL) {
01307
01308
01309
01310 ics = pCO->ijks[0] - (nghost/2) - 1;
01311 ice = pCO->ijke[0] + (nghost/2) + 1;
01312 if (pG->Nx[1] > 1) {
01313 jcs = pCO->ijks[1] - (nghost/2) - 1;
01314 jce = pCO->ijke[1] + (nghost/2) + 1;
01315 } else {
01316 jcs = pCO->ijks[1];
01317 jce = pCO->ijke[1];
01318 }
01319 if (pG->Nx[2] > 1) {
01320 kcs = pCO->ijks[2] - (nghost/2) - 1;
01321 kce = pCO->ijke[2] + (nghost/2) + 1;
01322 } else {
01323 kcs = pCO->ijks[2];
01324 kce = pCO->ijke[2];
01325 }
01326 if (dim == 0) ice = pCO->ijks[0];
01327 if (dim == 1) ics = pCO->ijke[0];
01328 if (dim == 2) jce = pCO->ijks[1];
01329 if (dim == 3) jcs = pCO->ijke[1];
01330 if (dim == 4) kce = pCO->ijks[2];
01331 if (dim == 5) kcs = pCO->ijke[2];
01332
01333
01334
01335
01336 for (k=kcs; k<=kce; k++) {
01337 for (j=jcs; j<=jce; j++) {
01338 for (i=ics; i<=ice; i++) {
01339 *(pSnd++) = pG->U[k][j][i].d;
01340 *(pSnd++) = pG->U[k][j][i].M1;
01341 *(pSnd++) = pG->U[k][j][i].M2;
01342 *(pSnd++) = pG->U[k][j][i].M3;
01343 #ifndef BAROTROPIC
01344 *(pSnd++) = pG->U[k][j][i].E;
01345 #endif
01346 #ifdef MHD
01347 *(pSnd++) = pG->U[k][j][i].B1c;
01348 *(pSnd++) = pG->U[k][j][i].B2c;
01349 *(pSnd++) = pG->U[k][j][i].B3c;
01350 *(pSnd++) = pG->B1i[k][j][i];
01351 *(pSnd++) = pG->B2i[k][j][i];
01352 *(pSnd++) = pG->B3i[k][j][i];
01353 #endif
01354 #if (NSCALARS > 0)
01355 for (n=0; n<NSCALARS; n++) {
01356 *(pSnd++) = pG->U[k][j][i].s[n];
01357 }
01358 #endif
01359 }}}
01360 }
01361 }
01362
01363
01364
01365 #ifdef MPI_PARALLEL
01366 if (ncg >= pG->NmyCGrid) {
01367 mIndex = ncg - pG->NmyCGrid;
01368 ierr = MPI_Isend(&(send_bufP[pCO->DomN][start_addrP[pCO->DomN]]),
01369 pG->CGrid[ncg].nWordsP, MPI_DOUBLE, pG->CGrid[ncg].ID, nd,
01370 pM->Domain[nl][nd].Comm_Children, &(send_rq[nd][mIndex]));
01371 }
01372 #endif
01373
01374 start_addrP[pCO->DomN] += pG->CGrid[ncg].nWordsP;
01375
01376 }
01377 }}
01378
01379
01380
01381
01382
01383
01384 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01385
01386 if (pM->Domain[nl][nd].Grid != NULL) {
01387 pG=pM->Domain[nl][nd].Grid;
01388 rbufN = (nl % 2);
01389
01390 for (npg=0; npg<(pG->NPGrid); npg++){
01391
01392
01393
01394 if (npg < pG->NmyPGrid) {
01395 pPO = (GridOvrlpS*)&(pG->PGrid[npg]);
01396 pRcv = (double*)&(recv_bufP[rbufN][nd][0]);
01397 } else {
01398
01399 #ifdef MPI_PARALLEL
01400
01401
01402
01403 mCount = pG->NPGrid - pG->NmyPGrid;
01404 ierr = MPI_Waitany(mCount,recv_rq[nl][nd],&mIndex,MPI_STATUS_IGNORE);
01405 if(mIndex == MPI_UNDEFINED){
01406 ath_error("[Prolong]: Invalid request index nl=%i nd=%i\n",nl,nd);
01407 }
01408
01409
01410
01411
01412 mAddress = 0;
01413 mIndex += pG->NmyPGrid;
01414 for (i=0; i<mIndex; i++) mAddress += pG->PGrid[i].nWordsP;
01415 pPO = (GridOvrlpS*)&(pG->PGrid[mIndex]);
01416 pRcv = (double*)&(recv_bufP[rbufN][nd][mAddress]);
01417
01418 #else
01419
01420
01421 ath_error("[Prolong]: no Parent Grid on Domain[%d][%d]\n",nl,nd);
01422 #endif
01423 }
01424
01425
01426
01427
01428 for (dim=0; dim<(2*nDim); dim++){
01429 if (pPO->myFlx[dim] != NULL) {
01430
01431
01432
01433
01434
01435 if (dim == 0 || dim == 1) {
01436 ngz1 = (nghost/2) + 2;
01437 id = 0;
01438 } else {
01439 ngz1 = (pPO->ijke[0] - pPO->ijks[0] + 1)/2 + nghost + 2;
01440 }
01441
01442 if (dim == 2 || dim == 3) {
01443 ngz2 = (nghost/2) + 2;
01444 id = 1;
01445 } else {
01446 ngz2 = (pPO->ijke[1] - pPO->ijks[1] + 1)/2 + nghost + 2;
01447 }
01448
01449 if (dim == 4 || dim == 5) {
01450 ngz3 = (nghost/2) + 2;
01451 id = 2;
01452 } else {
01453 ngz3 = (pPO->ijke[2] - pPO->ijks[2] + 1)/2 + nghost + 2;
01454 }
01455
01456 igzs = 0;
01457 igze = ngz1-1;
01458 if (pG->Nx[1] > 1) {
01459 jgzs = 0;
01460 jgze = ngz2-1;
01461 mend = 1;
01462 } else {
01463 ngz2 = 1;
01464 jgzs = 1;
01465 jgze = 1;
01466 mend = 0;
01467 }
01468 if (pG->Nx[2] > 1) {
01469 kgzs = 0;
01470 kgze = ngz3-1;
01471 nend = 1;
01472 } else {
01473 ngz3 = 1;
01474 kgzs = 1;
01475 kgze = 1;
01476 nend = 0;
01477 }
01478
01479
01480
01481 for (k=kgzs; k<=kgze; k++) {
01482 for (j=jgzs; j<=jgze; j++) {
01483 for (i=igzs; i<=igze; i++) {
01484 GZ[id][k][j][i].d = *(pRcv++);
01485 GZ[id][k][j][i].M1 = *(pRcv++);
01486 GZ[id][k][j][i].M2 = *(pRcv++);
01487 GZ[id][k][j][i].M3 = *(pRcv++);
01488 #ifndef BAROTROPIC
01489 GZ[id][k][j][i].E = *(pRcv++);
01490 #endif
01491 #ifdef MHD
01492 GZ[id][k][j][i].B1c = *(pRcv++);
01493 GZ[id][k][j][i].B2c = *(pRcv++);
01494 GZ[id][k][j][i].B3c = *(pRcv++);
01495 BFld[id][k][j][i].x = *(pRcv++);
01496 BFld[id][k][j][i].y = *(pRcv++);
01497 BFld[id][k][j][i].z = *(pRcv++);
01498 #endif
01499 #if (NSCALARS > 0)
01500 for (n=0; n<NSCALARS; n++) {
01501 GZ[id][k][j][i].s[n] = *(pRcv++);
01502 }
01503 #endif
01504 }}}
01505
01506
01507
01508 if (nDim == 1) {
01509 for (i=igzs; i<=igze; i++) {
01510 GZ[id][1][0][i] = GZ[id][1][1][i];
01511 GZ[id][1][2][i] = GZ[id][1][1][i];
01512 GZ[id][0][1][i] = GZ[id][1][1][i];
01513 GZ[id][2][1][i] = GZ[id][1][1][i];
01514 }
01515 }
01516
01517 if (nDim == 2) {
01518 for (j=jgzs; j<=jgze; j++) {
01519 for (i=igzs; i<=igze; i++) {
01520 GZ[id][0][j][i] = GZ[id][1][j][i];
01521 GZ[id][2][j][i] = GZ[id][1][j][i];
01522 }}
01523 #ifdef MHD
01524 for (j=jgzs; j<=jgze; j++) {
01525 for (i=igzs; i<=igze; i++) {
01526 BFld[id][0][j][i] = BFld[id][1][j][i];
01527 BFld[id][2][j][i] = BFld[id][1][j][i];
01528 }}
01529 #endif
01530 }
01531
01532
01533
01534
01535 ips = pPO->ijks[0] - nghost;
01536 ipe = pPO->ijke[0] + nghost;
01537 if (pG->Nx[1] > 1) {
01538 jps = pPO->ijks[1] - nghost;
01539 jpe = pPO->ijke[1] + nghost;
01540 } else {
01541 jps = pPO->ijks[1];
01542 jpe = pPO->ijke[1];
01543 }
01544 if (pG->Nx[2] > 1) {
01545 kps = pPO->ijks[2] - nghost;
01546 kpe = pPO->ijke[2] + nghost;
01547 } else {
01548 kps = pPO->ijks[2];
01549 kpe = pPO->ijke[2];
01550 }
01551 if (dim == 0) {ipe = pPO->ijks[0] - 1;}
01552 if (dim == 1) {ips = pPO->ijke[0] + 1;}
01553 if (dim == 2) {jpe = pPO->ijks[1] - 1;}
01554 if (dim == 3) {jps = pPO->ijke[1] + 1;}
01555 if (dim == 4) {kpe = pPO->ijks[2] - 1;}
01556 if (dim == 5) {kps = pPO->ijke[2] + 1;}
01557
01558
01559
01560 for (k=kps, kk=1; k<=kpe; k+=2, kk++) {
01561 for (j=jps, jj=1; j<=jpe; j+=2, jj++) {
01562 for (i=ips, ii=1; i<=ipe; i+=2, ii++) {
01563
01564 ProCon(GZ[id][kk][jj][ii-1],GZ[id][kk][jj][ii],GZ[id][kk][jj][ii+1],
01565 GZ[id][kk][jj-1][ii], GZ[id][kk][jj+1][ii],
01566 GZ[id][kk-1][jj][ii], GZ[id][kk+1][jj][ii],
01567 ProlongedC);
01568
01569
01570
01571 for (n=0; n<=nend; n++) {
01572 for (m=0; m<=mend; m++) {
01573 for (l=0; l<=1; l++) {
01574 pG->U[k+n][j+m][i+l].d = ProlongedC[n][m][l].d;
01575 pG->U[k+n][j+m][i+l].M1 = ProlongedC[n][m][l].M1;
01576 pG->U[k+n][j+m][i+l].M2 = ProlongedC[n][m][l].M2;
01577 pG->U[k+n][j+m][i+l].M3 = ProlongedC[n][m][l].M3;
01578 #ifndef BAROTROPIC
01579 pG->U[k+n][j+m][i+l].E = ProlongedC[n][m][l].E;
01580 #endif
01581 #ifdef MHD
01582 pG->U[k+n][j+m][i+l].B1c = ProlongedC[n][m][l].B1c;
01583 pG->U[k+n][j+m][i+l].B2c = ProlongedC[n][m][l].B2c;
01584 pG->U[k+n][j+m][i+l].B3c = ProlongedC[n][m][l].B3c;
01585 #endif
01586 #if (NSCALARS > 0)
01587 for (ns=0; ns<NSCALARS; ns++)
01588 pG->U[k+n][j+m][i+l].s[ns] = ProlongedC[n][m][l].s[ns];
01589 #endif
01590 }}}
01591
01592 #ifdef MHD
01593
01594
01595
01596 if (nDim == 1) {
01597 for (l=0; l<=1; l++) {
01598 pG->B1i[k][j][i+l] = pG->U[k][j][i+l].B1c;
01599 pG->B2i[k][j][i+l] = pG->U[k][j][i+l].B2c;
01600 pG->B3i[k][j][i+l] = pG->U[k][j][i+l].B3c;
01601 }
01602 } else {
01603 for (n=0; n<3; n++) {
01604 for (m=0; m<3; m++) {
01605 for (l=0; l<3; l++) {
01606 ProlongedF[n][m][l].x = 0.0;
01607 ProlongedF[n][m][l].y = 0.0;
01608 ProlongedF[n][m][l].z = 0.0;
01609 }}}
01610 }
01611
01612
01613
01614 if (nDim == 2 || nDim ==3) {
01615 for (n=0; n<3; n++) {
01616 for (m=0; m<3; m++) {
01617 for (l=0; l<3; l++) {
01618 BGZ[n][m][l].x = BFld[id][kk+(n-1)][jj+(m-1)][ii+(l-1)].x;
01619 BGZ[n][m][l].y = BFld[id][kk+(n-1)][jj+(m-1)][ii+(l-1)].y;
01620 BGZ[n][m][l].z = BFld[id][kk+(n-1)][jj+(m-1)][ii+(l-1)].z;
01621 }}}
01622
01623
01624
01625
01626
01627
01628
01629 if ((dim == 0) &&
01630 (i == (ipe-1)) &&
01631 ((j >= (jps+nghost)) || (pPO->myFlx[2]==NULL)) &&
01632 ((j < (jpe-nghost)) || (pPO->myFlx[3]==NULL)) ){
01633 ProlongedF[0][0][2].x = pG->B1i[k][j ][i+2];
01634 ProlongedF[0][1][2].x = pG->B1i[k][j+1][i+2];
01635 ProlongedF[1][0][2].x = pG->B1i[k][j ][i+2];
01636 ProlongedF[1][1][2].x = pG->B1i[k][j+1][i+2];
01637 if ((nDim == 3) &&
01638 ((k >= (kps+nghost)) || (pPO->myFlx[4]==NULL)) &&
01639 ((k < (kpe-nghost)) || (pPO->myFlx[5]==NULL)) ){
01640 ProlongedF[1][0][2].x = pG->B1i[k+1][j ][i+2];
01641 ProlongedF[1][1][2].x = pG->B1i[k+1][j+1][i+2];
01642 }
01643 }
01644
01645
01646 if ((dim == 1) &&
01647 (i == ips) &&
01648 ((j >= (jps+nghost)) || (pPO->myFlx[2]==NULL)) &&
01649 ((j < (jpe-nghost)) || (pPO->myFlx[3]==NULL)) ){
01650 ProlongedF[0][0][0].x = pG->B1i[k][j ][i];
01651 ProlongedF[0][1][0].x = pG->B1i[k][j+1][i];
01652 ProlongedF[1][0][0].x = pG->B1i[k][j ][i];
01653 ProlongedF[1][1][0].x = pG->B1i[k][j+1][i];
01654 if ((nDim == 3) &&
01655 ((k >= (kps+nghost)) || (pPO->myFlx[4]==NULL)) &&
01656 ((k < (kpe-nghost)) || (pPO->myFlx[5]==NULL)) ){
01657 ProlongedF[1][0][0].x = pG->B1i[k+1][j ][i];
01658 ProlongedF[1][1][0].x = pG->B1i[k+1][j+1][i];
01659 }
01660 }
01661
01662
01663 if ((dim == 2) &&
01664 (j == (jpe-1)) &&
01665 ((i >= (ips+nghost)) || (pPO->myFlx[0]==NULL)) &&
01666 ((i < (ipe-nghost)) || (pPO->myFlx[1]==NULL)) ){
01667 ProlongedF[0][2][0].y = pG->B2i[k][j+2][i ];
01668 ProlongedF[0][2][1].y = pG->B2i[k][j+2][i+1];
01669 ProlongedF[1][2][0].y = pG->B2i[k][j+2][i ];
01670 ProlongedF[1][2][1].y = pG->B2i[k][j+2][i+1];
01671 if ((nDim == 3) &&
01672 ((k >= (kps+nghost)) || (pPO->myFlx[4]==NULL)) &&
01673 ((k < (kpe-nghost)) || (pPO->myFlx[5]==NULL)) ){
01674 ProlongedF[1][2][0].y = pG->B2i[k+1][j+2][i ];
01675 ProlongedF[1][2][1].y = pG->B2i[k+1][j+2][i+1];
01676 }
01677 }
01678
01679
01680 if ((dim == 3) &&
01681 (j == jps) &&
01682 ((i >= (ips+nghost)) || (pPO->myFlx[0]==NULL)) &&
01683 ((i < (ipe-nghost)) || (pPO->myFlx[1]==NULL)) ){
01684 ProlongedF[0][0][0].y = pG->B2i[k][j][i ];
01685 ProlongedF[0][0][1].y = pG->B2i[k][j][i+1];
01686 ProlongedF[1][0][0].y = pG->B2i[k][j][i ];
01687 ProlongedF[1][0][1].y = pG->B2i[k][j][i+1];
01688 if ((nDim == 3) &&
01689 ((k >= (kps+nghost)) || (pPO->myFlx[4]==NULL)) &&
01690 ((k < (kpe-nghost)) || (pPO->myFlx[5]==NULL)) ){
01691 ProlongedF[1][0][0].y = pG->B2i[k+1][j][i ];
01692 ProlongedF[1][0][1].y = pG->B2i[k+1][j][i+1];
01693 }
01694 }
01695
01696
01697 if ((dim == 4) &&
01698 (k == (kpe-1)) &&
01699 ((i >= (ips+nghost)) || (pPO->myFlx[0]==NULL)) &&
01700 ((i < (ipe-nghost)) || (pPO->myFlx[1]==NULL)) &&
01701 ((j >= (jps+nghost)) || (pPO->myFlx[2]==NULL)) &&
01702 ((j < (jpe-nghost)) || (pPO->myFlx[3]==NULL)) ){
01703 ProlongedF[2][0][0].z = pG->B3i[k+2][j ][i ];
01704 ProlongedF[2][0][1].z = pG->B3i[k+2][j ][i+1];
01705 ProlongedF[2][1][0].z = pG->B3i[k+2][j+1][i ];
01706 ProlongedF[2][1][1].z = pG->B3i[k+2][j+1][i+1];
01707 }
01708
01709
01710 if ((dim == 5) &&
01711 (k == kps) &&
01712 ((i >= (ips+nghost)) || (pPO->myFlx[0]==NULL)) &&
01713 ((i < (ipe-nghost)) || (pPO->myFlx[1]==NULL)) &&
01714 ((j >= (jps+nghost)) || (pPO->myFlx[2]==NULL)) &&
01715 ((j < (jpe-nghost)) || (pPO->myFlx[3]==NULL)) ){
01716 ProlongedF[0][0][0].z = pG->B3i[k][j ][i ];
01717 ProlongedF[0][0][1].z = pG->B3i[k][j ][i+1];
01718 ProlongedF[0][1][0].z = pG->B3i[k][j+1][i ];
01719 ProlongedF[0][1][1].z = pG->B3i[k][j+1][i+1];
01720 }
01721
01722 ProFld(BGZ, ProlongedF, pG->dx1, pG->dx2, pG->dx3);
01723
01724 for (n=0; n<=nend; n++) {
01725 for (m=0; m<=mend; m++) {
01726 for (l=0; l<=1; l++) {
01727 if (dim != 1 || (i+l) != ips)
01728 pG->B1i[k+n][j+m][i+l] = ProlongedF[n][m][l].x;
01729 if (dim != 3 || (j+m) != jps)
01730 pG->B2i[k+n][j+m][i+l] = ProlongedF[n][m][l].y;
01731 if (dim != 5 || (k+n) != kps)
01732 pG->B3i[k+n][j+m][i+l] = ProlongedF[n][m][l].z;
01733
01734 pG->U[k+n][j+m][i+l].B1c =
01735 0.5*(ProlongedF[n][m][l].x + ProlongedF[n][m][l+1].x);
01736 pG->U[k+n][j+m][i+l].B2c =
01737 0.5*(ProlongedF[n][m][l].y + ProlongedF[n][m+1][l].y);
01738 pG->U[k+n][j+m][i+l].B3c =
01739 0.5*(ProlongedF[n][m][l].z + ProlongedF[n+1][m][l].z);
01740 }}}
01741 }
01742
01743 #endif
01744 }}}
01745
01746 }
01747 }
01748
01749 }
01750 }}
01751
01752
01753
01754
01755
01756
01757 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01758 if (pM->Domain[nl][nd].Grid != NULL) {
01759 pG=pM->Domain[nl][nd].Grid;
01760 rbufN = ((nl+1) % 2);
01761
01762
01763
01764
01765 for (ncg=0; ncg<(pG->NmyCGrid); ncg++){
01766 pCO=(GridOvrlpS*)&(pG->CGrid[ncg]);
01767
01768 for (i=0; i<pCO->nWordsP; i++) {
01769 recv_bufP[rbufN][pCO->DomN][i]=send_bufP[pCO->DomN][i];
01770 }
01771 }
01772 }
01773 }
01774
01775 #ifdef MPI_PARALLEL
01776
01777
01778 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01779 if (pM->Domain[nl][nd].Grid != NULL) {
01780 pG=pM->Domain[nl][nd].Grid;
01781
01782 if (pG->NCGrid > pG->NmyCGrid) {
01783 mCount = pG->NCGrid - pG->NmyCGrid;
01784 ierr = MPI_Waitall(mCount, send_rq[nd], MPI_STATUS_IGNORE);
01785 }
01786 }
01787 }
01788 #endif
01789
01790 }
01791 }
01792
01793
01794
01795
01796
01797
01798
01799 void SMR_init(MeshS *pM)
01800 {
01801 int nl,nd,sendRC,recvRC,sendP,recvP,npg,ncg;
01802 int max_sendRC=1,max_recvRC=1,max_sendP=1,max_recvP=1;
01803 int max1=0,max2=0,max3=0,maxCG=1;
01804 #ifdef MHD
01805 int ngh1;
01806 #endif
01807 GridS *pG;
01808
01809 maxND=1;
01810 for (nl=0; nl<(pM->NLevels); nl++) maxND=MAX(maxND,pM->DomainsPerLevel[nl]);
01811 if((start_addrP = (int*)calloc_1d_array(maxND,sizeof(int))) == NULL)
01812 ath_error("[SMR_init]:Failed to allocate start_addrP\n");
01813
01814
01815
01816
01817 for (nl=0; nl<(pM->NLevels); nl++){
01818 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01819 sendRC=0;
01820 recvRC=0;
01821 sendP =0;
01822 recvP =0;
01823
01824 if (pM->Domain[nl][nd].Grid != NULL) {
01825 pG=pM->Domain[nl][nd].Grid;
01826
01827 for (npg=0; npg<pG->NPGrid; npg++){
01828 sendRC += pG->PGrid[npg].nWordsRC;
01829 recvP += pG->PGrid[npg].nWordsP;
01830 }
01831 for (ncg=0; ncg<pG->NCGrid; ncg++){
01832 recvRC += pG->CGrid[ncg].nWordsRC;
01833 sendP += pG->CGrid[ncg].nWordsP;
01834 }
01835
01836 max_sendRC = MAX(max_sendRC,sendRC);
01837 max_recvRC = MAX(max_recvRC,recvRC);
01838 max_sendP = MAX(max_sendP ,sendP );
01839 max_recvP = MAX(max_recvP ,recvP );
01840 max1 = MAX(max1,(pG->Nx[0]+1));
01841 max2 = MAX(max2,(pG->Nx[1]+1));
01842 max3 = MAX(max3,(pG->Nx[2]+1));
01843 maxCG = MAX(maxCG,pG->NCGrid);
01844 }
01845 }
01846 }
01847
01848
01849
01850 if((send_bufRC =
01851 (double**)calloc_2d_array(maxND,max_sendRC,sizeof(double))) == NULL)
01852 ath_error("[SMR_init]:Failed to allocate send_bufRC\n");
01853
01854 #ifdef MPI_PARALLEL
01855 if((recv_bufRC =
01856 (double***)calloc_3d_array(2,maxND,max_recvRC,sizeof(double))) == NULL)
01857 ath_error("[SMR_init]: Failed to allocate recv_bufRC\n");
01858 if((recv_rq = (MPI_Request***)
01859 calloc_3d_array(pM->NLevels,maxND,maxCG,sizeof(MPI_Request))) == NULL)
01860 ath_error("[SMR_init]: Failed to allocate recv MPI_Request array\n");
01861 if((send_rq = (MPI_Request**)
01862 calloc_2d_array(maxND,maxCG,sizeof(MPI_Request))) == NULL)
01863 ath_error("[SMR_init]: Failed to allocate send MPI_Request array\n");
01864 #endif
01865
01866 #ifdef MHD
01867 if((SMRemf1 =
01868 (Real**)calloc_2d_array(MAX(max2,max3),max1,sizeof(Real))) == NULL)
01869 ath_error("[smr_init]Failed to calloc_2d_array for SMRemf1\n");;
01870 if((SMRemf2 =
01871 (Real**)calloc_2d_array(MAX(max2,max3),MAX(max1,max2),sizeof(Real))) ==NULL)
01872 ath_error("[smr_init]Failed to calloc_2d_array for SMRemf2\n");;
01873 if((SMRemf3 =
01874 (Real**)calloc_2d_array(max3,MAX(max1,max2),sizeof(Real))) == NULL)
01875 ath_error("[smr_init]Failed to calloc_2d_array for SMRemf3\n");;
01876 #endif
01877
01878
01879
01880 if((send_bufP =
01881 (double**)calloc_2d_array(maxND,max_sendP,sizeof(double))) == NULL)
01882 ath_error("[SMR_init]:Failed to allocate send_bufP\n");
01883
01884 if((recv_bufP =
01885 (double***)calloc_3d_array(2,maxND,max_recvP,sizeof(double))) == NULL)
01886 ath_error("[SMR_init]: Failed to allocate recv_bufP\n");
01887
01888 max1 += 2*nghost;
01889 max2 += 2*nghost;
01890 max3 += 2*nghost;
01891
01892 if((GZ[0]=(ConsS***)calloc_3d_array(max3,max2,nghost,sizeof(ConsS)))
01893 ==NULL) ath_error("[SMR_init]:Failed to allocate GZ[0]C\n");
01894 if((GZ[1]=(ConsS***)calloc_3d_array(max3,nghost,max1,sizeof(ConsS)))
01895 ==NULL) ath_error("[SMR_init]:Failed to allocate GZ[1]C\n");
01896 if((GZ[2]=(ConsS***)calloc_3d_array(nghost,max2,max1,sizeof(ConsS)))
01897 ==NULL) ath_error("[SMR_init]:Failed to allocate GZ[2]C\n");
01898 #ifdef MHD
01899 ngh1 = nghost + 1;
01900 if((BFld[0]=(Real3Vect***)calloc_3d_array(max3,max2,ngh1,sizeof(Real3Vect)))
01901 ==NULL) ath_error("[SMR_init]:Failed to allocate BFld[0]C\n");
01902 if((BFld[1]=(Real3Vect***)calloc_3d_array(max3,ngh1,max1,sizeof(Real3Vect)))
01903 ==NULL) ath_error("[SMR_init]:Failed to allocate BFld[1]C\n");
01904 if((BFld[2]=(Real3Vect***)calloc_3d_array(ngh1,max2,max1,sizeof(Real3Vect)))
01905 ==NULL) ath_error("[SMR_init]:Failed to allocate BFld[2]C\n");
01906 #endif
01907
01908 return;
01909 }
01910
01911
01912
01913
01914
01915
01916
01917
01918 void ProCon(const ConsS Uim1,const ConsS Ui, const ConsS Uip1,
01919 const ConsS Ujm1,const ConsS Ujp1,
01920 const ConsS Ukm1,const ConsS Ukp1, ConsS PCon[][2][2])
01921 {
01922 int i,j,k;
01923 Real dq1,dq2,dq3,Pim1,Pi,Pip1,Pjm1,Pjp1,Pkm1,Pkp1;
01924 #ifdef SPECIAL_RELATIVITY
01925 PrimS W;
01926 Real Vsq;
01927 int fail,dfail,Pfail,Vfail;
01928 #endif
01929 #if (NSCALARS > 0)
01930 int n;
01931 #endif
01932
01933
01934 #ifdef FIRST_ORDER
01935
01936 for (k=0; k<2; k++){
01937 for (j=0; j<2; j++){
01938 for (i=0; i<2; i++){
01939 PCon[k][j][i].d = Ui.d;
01940 PCon[k][j][i].M1 = Ui.M1;
01941 PCon[k][j][i].M2 = Ui.M2;
01942 PCon[k][j][i].M3 = Ui.M3;
01943 #ifndef BAROTROPIC
01944 PCon[k][j][i].E = Ui.E;
01945 #endif
01946 #ifdef MHD
01947 PCon[k][j][i].B1c = Ui.B1c;
01948 PCon[k][j][i].B2c = Ui.B2c;
01949 PCon[k][j][i].B3c = Ui.B3c;
01950 #endif
01951 #if (NSCALARS > 0)
01952 for (n=0; n<NSCALARS; n++) PCon[k][j][i].s[n] = Ui.s[n];
01953 #endif
01954 }}}
01955
01956
01957 #else
01958
01959
01960 dq1 = mcd_slope(Uim1.d, Ui.d, Uip1.d);
01961 dq2 = mcd_slope(Ujm1.d, Ui.d, Ujp1.d);
01962 dq3 = mcd_slope(Ukm1.d, Ui.d, Ukp1.d);
01963 for (k=0; k<2; k++){
01964 for (j=0; j<2; j++){
01965 for (i=0; i<2; i++){
01966 PCon[k][j][i].d = Ui.d
01967 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
01968 }}}
01969
01970
01971 dq1 = mcd_slope(Uim1.M1, Ui.M1, Uip1.M1);
01972 dq2 = mcd_slope(Ujm1.M1, Ui.M1, Ujp1.M1);
01973 dq3 = mcd_slope(Ukm1.M1, Ui.M1, Ukp1.M1);
01974 for (k=0; k<2; k++){
01975 for (j=0; j<2; j++){
01976 for (i=0; i<2; i++){
01977 PCon[k][j][i].M1 = Ui.M1
01978 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
01979 }}}
01980
01981
01982 dq1 = mcd_slope(Uim1.M2, Ui.M2, Uip1.M2);
01983 dq2 = mcd_slope(Ujm1.M2, Ui.M2, Ujp1.M2);
01984 dq3 = mcd_slope(Ukm1.M2, Ui.M2, Ukp1.M2);
01985 for (k=0; k<2; k++){
01986 for (j=0; j<2; j++){
01987 for (i=0; i<2; i++){
01988 PCon[k][j][i].M2 = Ui.M2
01989 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
01990 }}}
01991
01992
01993 dq1 = mcd_slope(Uim1.M3, Ui.M3, Uip1.M3);
01994 dq2 = mcd_slope(Ujm1.M3, Ui.M3, Ujp1.M3);
01995 dq3 = mcd_slope(Ukm1.M3, Ui.M3, Ukp1.M3);
01996 for (k=0; k<2; k++){
01997 for (j=0; j<2; j++){
01998 for (i=0; i<2; i++){
01999 PCon[k][j][i].M3 = Ui.M3
02000 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02001 }}}
02002
02003 #ifdef MHD
02004
02005 dq1 = mcd_slope(Uim1.B1c, Ui.B1c, Uip1.B1c);
02006 dq2 = mcd_slope(Ujm1.B1c, Ui.B1c, Ujp1.B1c);
02007 dq3 = mcd_slope(Ukm1.B1c, Ui.B1c, Ukp1.B1c);
02008 for (k=0; k<2; k++){
02009 for (j=0; j<2; j++){
02010 for (i=0; i<2; i++){
02011 PCon[k][j][i].B1c = Ui.B1c
02012 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02013 }}}
02014
02015
02016 dq1 = mcd_slope(Uim1.B2c, Ui.B2c, Uip1.B2c);
02017 dq2 = mcd_slope(Ujm1.B2c, Ui.B2c, Ujp1.B2c);
02018 dq3 = mcd_slope(Ukm1.B2c, Ui.B2c, Ukp1.B2c);
02019 for (k=0; k<2; k++){
02020 for (j=0; j<2; j++){
02021 for (i=0; i<2; i++){
02022 PCon[k][j][i].B2c = Ui.B2c
02023 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02024 }}}
02025
02026
02027 dq1 = mcd_slope(Uim1.B3c, Ui.B3c, Uip1.B3c);
02028 dq2 = mcd_slope(Ujm1.B3c, Ui.B3c, Ujp1.B3c);
02029 dq3 = mcd_slope(Ukm1.B3c, Ui.B3c, Ukp1.B3c);
02030 for (k=0; k<2; k++){
02031 for (j=0; j<2; j++){
02032 for (i=0; i<2; i++){
02033 PCon[k][j][i].B3c = Ui.B3c
02034 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02035 }}}
02036 #endif
02037
02038 #ifndef BAROTROPIC
02039 #ifdef SPECIAL_RELATIVITY
02040
02041
02042
02043
02044
02045 dq1 = mcd_slope(Uim1.E, Ui.E, Uip1.E);
02046 dq2 = mcd_slope(Ujm1.E, Ui.E, Ujp1.E);
02047 dq3 = mcd_slope(Ukm1.E, Ui.E, Ukp1.E);
02048 for (k=0; k<2; k++){
02049 for (j=0; j<2; j++){
02050 for (i=0; i<2; i++){
02051 PCon[k][j][i].E = Ui.E
02052 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02053 }}}
02054
02055 #else
02056
02057
02058
02059 Pi = Ui.E - 0.5*(SQR(Ui.M1 ) + SQR(Ui.M2 ) + SQR(Ui.M3 ))/Ui.d;
02060 Pim1 = Uim1.E - 0.5*(SQR(Uim1.M1) + SQR(Uim1.M2) + SQR(Uim1.M3))/Uim1.d;
02061 Pip1 = Uip1.E - 0.5*(SQR(Uip1.M1) + SQR(Uip1.M2) + SQR(Uip1.M3))/Uip1.d;
02062 #ifdef MHD
02063 Pi -= 0.5*(SQR(Ui.B1c ) + SQR(Ui.B2c ) + SQR(Ui.B3c ));
02064 Pim1 -= 0.5*(SQR(Uim1.B1c) + SQR(Uim1.B2c) + SQR(Uim1.B3c));
02065 Pip1 -= 0.5*(SQR(Uip1.B1c) + SQR(Uip1.B2c) + SQR(Uip1.B3c));
02066 #endif
02067 dq1 = mcd_slope(Pim1, Pi, Pip1);
02068
02069 Pjm1 = Ujm1.E - 0.5*(SQR(Ujm1.M1) + SQR(Ujm1.M2) + SQR(Ujm1.M3))/Ujm1.d;
02070 Pjp1 = Ujp1.E - 0.5*(SQR(Ujp1.M1) + SQR(Ujp1.M2) + SQR(Ujp1.M3))/Ujp1.d;
02071 #ifdef MHD
02072 Pjm1 -= 0.5*(SQR(Ujm1.B1c) + SQR(Ujm1.B2c) + SQR(Ujm1.B3c));
02073 Pjp1 -= 0.5*(SQR(Ujp1.B1c) + SQR(Ujp1.B2c) + SQR(Ujp1.B3c));
02074 #endif
02075 dq2 = mcd_slope(Pjm1, Pi, Pjp1);
02076
02077 Pkm1 = Ukm1.E - 0.5*(SQR(Ukm1.M1) + SQR(Ukm1.M2) + SQR(Ukm1.M3))/Ukm1.d;
02078 Pkp1 = Ukp1.E - 0.5*(SQR(Ukp1.M1) + SQR(Ukp1.M2) + SQR(Ukp1.M3))/Ukp1.d;
02079 #ifdef MHD
02080 Pkm1 -= 0.5*(SQR(Ukm1.B1c) + SQR(Ukm1.B2c) + SQR(Ukm1.B3c));
02081 Pkp1 -= 0.5*(SQR(Ukp1.B1c) + SQR(Ukp1.B2c) + SQR(Ukp1.B3c));
02082 #endif
02083 dq3 = mcd_slope(Pkm1, Pi, Pkp1);
02084
02085 for (k=0; k<2; k++){
02086 for (j=0; j<2; j++){
02087 for (i=0; i<2; i++){
02088 PCon[k][j][i].E = Pi
02089 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;
02090 PCon[k][j][i].E += 0.5*(SQR(PCon[k][j][i].M1) + SQR(PCon[k][j][i].M2) +
02091 SQR(PCon[k][j][i].M3))/PCon[k][j][i].d;
02092 #ifdef MHD
02093 PCon[k][j][i].E += 0.5*(SQR(PCon[k][j][i].B1c) + SQR(PCon[k][j][i].B2c) +
02094 SQR(PCon[k][j][i].B3c));
02095 #endif
02096 }}}
02097
02098 #endif
02099 #endif
02100
02101 #if (NSCALARS > 0)
02102
02103 for (n=0; n<NSCALARS; n++) {
02104 dq1 = mcd_slope(Uim1.s[n], Ui.s[n], Uip1.s[n]);
02105 dq2 = mcd_slope(Ujm1.s[n], Ui.s[n], Ujp1.s[n]);
02106 dq3 = mcd_slope(Ukm1.s[n], Ui.s[n], Ukp1.s[n]);
02107 for (k=0; k<2; k++){
02108 for (j=0; j<2; j++){
02109 for (i=0; i<2; i++){
02110 PCon[k][j][i].s[n] = Ui.s[n]
02111 + (0.5*i - 0.25)*dq1 + (0.5*j - 0.25)*dq2 + (0.5*k - 0.25)*dq3;;
02112 }}}
02113 }
02114 #endif
02115
02116 #ifdef SPECIAL_RELATIVITY
02117
02118
02119 dfail = 0;
02120 Pfail = 0;
02121 Vfail = 0;
02122 fail = 0;
02123 for (k=0; k<2; k++){
02124 for (j=0; j<2; j++){
02125 for (i=0; i<2; i++){
02126 W = check_Prim(&(PCon[k][j][i]));
02127 Vsq = SQR(W.V1) + SQR(W.V2) + SQR(W.V3);
02128 if (W.d < 0.0){
02129 dfail++;
02130 fail = 1;
02131 }
02132 if (W.P < 0.0){
02133 Pfail++;
02134 fail = 1;
02135 }
02136 if (Vsq > 1.0){
02137 Vfail++;
02138 fail = 1;
02139 }
02140 }}}
02141
02142
02143
02144 if (fail) {
02145
02146 for (k=0; k<2; k++){
02147 for (j=0; j<2; j++){
02148 for (i=0; i<2; i++){
02149 PCon[k][j][i].d = Ui.d;
02150 PCon[k][j][i].M1 = Ui.M1;
02151 PCon[k][j][i].M2 = Ui.M2;
02152 PCon[k][j][i].M3 = Ui.M3;
02153 #ifndef BAROTROPIC
02154 PCon[k][j][i].E = Ui.E;
02155 #endif
02156 #ifdef MHD
02157 PCon[k][j][i].B1c = Ui.B1c;
02158 PCon[k][j][i].B2c = Ui.B2c;
02159 PCon[k][j][i].B3c = Ui.B3c;
02160 #endif
02161 #if (NSCALARS > 0)
02162 for (n=0; n<NSCALARS; n++) PCon[k][j][i].s[n] = Ui.s[n];
02163 #endif
02164 }}}
02165 }
02166
02167 dfail = 0;
02168 Pfail = 0;
02169 Vfail = 0;
02170 fail = 0;
02171
02172 #endif SPECIAL_RELATIVITY
02173
02174 #endif
02175 }
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185 #ifdef MHD
02186 void ProFld(Real3Vect BGZ[][3][3], Real3Vect PFld[][3][3],
02187 const Real dx1c, const Real dx2c, const Real dx3c)
02188 {
02189 int i,j,k;
02190 Real dBdx,dBdy,dBdz,Uxx,Vyy,Wzz,Uxyz,Vxyz,Wxyz;
02191
02192
02193
02194 if (PFld[0][0][0].x == 0.0) {
02195 dBdy = mcd_slope(BGZ[1][0][1].x, BGZ[1][1][1].x, BGZ[1][2][1].x);
02196 dBdz = mcd_slope(BGZ[0][1][1].x, BGZ[1][1][1].x, BGZ[2][1][1].x);
02197
02198 PFld[0][0][0].x = BGZ[1][1][1].x - 0.25*dBdy - 0.25*dBdz;
02199 PFld[0][1][0].x = BGZ[1][1][1].x + 0.25*dBdy - 0.25*dBdz;
02200 PFld[1][0][0].x = BGZ[1][1][1].x - 0.25*dBdy + 0.25*dBdz;
02201 PFld[1][1][0].x = BGZ[1][1][1].x + 0.25*dBdy + 0.25*dBdz;
02202 }
02203
02204
02205
02206 if (PFld[0][0][2].x == 0.0) {
02207 dBdy = mcd_slope(BGZ[1][0][2].x, BGZ[1][1][2].x, BGZ[1][2][2].x);
02208 dBdz = mcd_slope(BGZ[0][1][2].x, BGZ[1][1][2].x, BGZ[2][1][2].x);
02209
02210 PFld[0][0][2].x = BGZ[1][1][2].x - 0.25*dBdy - 0.25*dBdz;
02211 PFld[0][1][2].x = BGZ[1][1][2].x + 0.25*dBdy - 0.25*dBdz;
02212 PFld[1][0][2].x = BGZ[1][1][2].x - 0.25*dBdy + 0.25*dBdz;
02213 PFld[1][1][2].x = BGZ[1][1][2].x + 0.25*dBdy + 0.25*dBdz;
02214 }
02215
02216
02217
02218 if (PFld[0][0][0].y == 0.0) {
02219 dBdx = mcd_slope(BGZ[1][1][0].y, BGZ[1][1][1].y, BGZ[1][1][2].y);
02220 dBdz = mcd_slope(BGZ[0][1][1].y, BGZ[1][1][1].y, BGZ[2][1][1].y);
02221
02222 PFld[0][0][0].y = BGZ[1][1][1].y - 0.25*dBdx - 0.25*dBdz;
02223 PFld[0][0][1].y = BGZ[1][1][1].y + 0.25*dBdx - 0.25*dBdz;
02224 PFld[1][0][0].y = BGZ[1][1][1].y - 0.25*dBdx + 0.25*dBdz;
02225 PFld[1][0][1].y = BGZ[1][1][1].y + 0.25*dBdx + 0.25*dBdz;
02226 }
02227
02228
02229
02230 if (PFld[0][2][0].y == 0.0) {
02231 dBdx = mcd_slope(BGZ[1][2][0].y, BGZ[1][2][1].y, BGZ[1][2][2].y);
02232 dBdz = mcd_slope(BGZ[0][2][1].y, BGZ[1][2][1].y, BGZ[2][2][1].y);
02233
02234 PFld[0][2][0].y = BGZ[1][2][1].y - 0.25*dBdx - 0.25*dBdz;
02235 PFld[0][2][1].y = BGZ[1][2][1].y + 0.25*dBdx - 0.25*dBdz;
02236 PFld[1][2][0].y = BGZ[1][2][1].y - 0.25*dBdx + 0.25*dBdz;
02237 PFld[1][2][1].y = BGZ[1][2][1].y + 0.25*dBdx + 0.25*dBdz;
02238 }
02239
02240
02241
02242 if (PFld[0][0][0].z == 0.0) {
02243 dBdx = mcd_slope(BGZ[1][1][0].z, BGZ[1][1][1].z, BGZ[1][1][2].z);
02244 dBdy = mcd_slope(BGZ[1][0][1].z, BGZ[1][1][1].z, BGZ[1][2][1].z);
02245
02246 PFld[0][0][0].z = BGZ[1][1][1].z - 0.25*dBdx - 0.25*dBdy;
02247 PFld[0][0][1].z = BGZ[1][1][1].z + 0.25*dBdx - 0.25*dBdy;
02248 PFld[0][1][0].z = BGZ[1][1][1].z - 0.25*dBdx + 0.25*dBdy;
02249 PFld[0][1][1].z = BGZ[1][1][1].z + 0.25*dBdx + 0.25*dBdy;
02250 }
02251
02252
02253
02254 if (PFld[2][0][0].z == 0.0) {
02255 dBdx = mcd_slope(BGZ[2][1][0].z, BGZ[2][1][1].z, BGZ[2][1][2].z);
02256 dBdy = mcd_slope(BGZ[2][0][1].z, BGZ[2][1][1].z, BGZ[2][2][1].z);
02257
02258 PFld[2][0][0].z = BGZ[2][1][1].z - 0.25*dBdx - 0.25*dBdy;
02259 PFld[2][0][1].z = BGZ[2][1][1].z + 0.25*dBdx - 0.25*dBdy;
02260 PFld[2][1][0].z = BGZ[2][1][1].z - 0.25*dBdx + 0.25*dBdy;
02261 PFld[2][1][1].z = BGZ[2][1][1].z + 0.25*dBdx + 0.25*dBdy;
02262 }
02263
02264
02265
02266
02267
02268 Uxx = Vyy = Wzz = 0.0;
02269 Uxyz = Vxyz = Wxyz = 0.0;
02270 for(j=0; j<2; j++){
02271 for(i=0; i<2; i++){
02272 Uxx += (2*i-1)*((2*j-1)*dx3c*(PFld[0][2*j][i].y + PFld[1][2*j][i].y) +
02273 dx2c*(PFld[2][j ][i].z - PFld[0][j ][i].z) );
02274
02275 Vyy += (2*j-1)*( dx1c*(PFld[2][j][i ].z - PFld[0][j][i ].z) +
02276 (2*i-1)*dx3c*(PFld[0][j][2*i].x + PFld[1][j][2*i].x) );
02277
02278 Wzz += ((2*i-1)*dx2c*(PFld[1][j][2*i].x - PFld[0][j][2*i].x) +
02279 (2*j-1)*dx1c*(PFld[1][2*j][i].y - PFld[0][2*j][i].y) );
02280
02281 Uxyz += (2*i-1)*(2*j-1)*(PFld[1][j][2*i].x - PFld[0][j][2*i].x);
02282 Vxyz += (2*i-1)*(2*j-1)*(PFld[1][2*j][i].y - PFld[0][2*j][i].y);
02283 Wxyz += (2*i-1)*(2*j-1)*(PFld[2][j][i].z - PFld[0][j][i].z);
02284 }}
02285
02286
02287
02288 Uxx *= 0.125*dx1c;
02289 Vyy *= 0.125*dx2c;
02290 Wzz *= 0.125*dx3c;
02291 Uxyz *= 0.125*dx2c*dx3c/(dx2c*dx2c + dx3c*dx3c);
02292 Vxyz *= 0.125*dx1c*dx3c/(dx1c*dx1c + dx3c*dx3c);
02293 Wxyz *= 0.125*dx1c*dx2c/(dx1c*dx1c + dx2c*dx2c);
02294
02295
02296
02297 for(k=0; k<2; k++){
02298 for(j=0; j<2; j++){
02299 PFld[k][j][1].x = 0.5*(PFld[k][j][0].x + PFld[k][j][2].x) + Uxx/(dx2c*dx3c)
02300 + (2*k-1)*(dx3c/dx2c)*Vxyz + (2*j-1)*(dx2c/dx3c)*Wxyz;
02301 }}
02302
02303
02304
02305 for(k=0; k<2; k++){
02306 for(i=0; i<2; i++){
02307 PFld[k][1][i].y = 0.5*(PFld[k][0][i].y + PFld[k][2][i].y) + Vyy/(dx3c*dx1c)
02308 + (2*i-1)*(dx1c/dx3c)*Wxyz + (2*k-1)*(dx3c/dx1c)*Uxyz;
02309 }}
02310
02311
02312
02313 for(j=0; j<2; j++){
02314 for(i=0; i<2; i++){
02315 PFld[1][j][i].z = 0.5*(PFld[0][j][i].z + PFld[2][j][i].z) + Wzz/(dx1c*dx2c)
02316 + (2*j-1)*(dx2c/dx1c)*Uxyz + (2*i-1)*(dx1c/dx2c)*Vxyz;
02317 }}
02318
02319 }
02320 #endif
02321
02322
02323
02324
02325
02326
02327 #ifndef FIRST_ORDER
02328 static Real mcd_slope(const Real vl, const Real vc, const Real vr){
02329
02330 Real dvl = (vc - vl), dvr = (vr - vc);
02331 Real dv, dvm;
02332
02333 if(dvl > 0.0 && dvr > 0.0){
02334 dv = 2.0*(dvl < dvr ? dvl : dvr);
02335 dvm = 0.5*(dvl + dvr);
02336 return (dvm < dv ? dvm : dv);
02337 }
02338 else if(dvl < 0.0 && dvr < 0.0){
02339 dv = 2.0*(dvl > dvr ? dvl : dvr);
02340 dvm = 0.5*(dvl + dvr);
02341 return (dvm > dv ? dvm : dv);
02342 }
02343
02344 return 0.0;
02345 }
02346 #endif
02347
02348 #endif