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