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