Athena: integrators/integrate_3d_vl

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_3d_vl_sr.c
00004  *  \brief Integrate MHD equations using 3D version of the directionally
00005  *   unsplit MUSCL-Hancock (VL) integrator. 
00006  *
00007  * PURPOSE: Integrate MHD equations using 3D version of the directionally
00008  *   unsplit MUSCL-Hancock (VL) integrator.  The variables updated are:
00009  *    - U.[d,M1,M2,M3,E,B1c,B2c,B3c,s] -- where U is of type ConsS
00010  *    - B1i, B2i, B3i  -- interface magnetic field
00011  *   Also adds gravitational source terms, self-gravity, and the H-correction
00012  *   of Sanders et al.
00013  *   - For adb hydro, requires (9*Cons1DS + 3*Real + 1*ConsS) = 53 3D arrays
00014  *   - For adb mhd, requires   (9*Cons1DS + 9*Real + 1*ConsS) = 80 3D arrays
00015  *
00016  * REFERENCE: 
00017  * - J.M Stone & T.A. Gardiner, "A simple, unsplit Godunov method
00018  *   for multidimensional MHD", NewA 14, 139 (2009)
00019  *
00020  * - R. Sanders, E. Morano, & M.-C. Druguet, "Multidimensinal dissipation for
00021  *   upwind schemes: stability and applications to gas dynamics", JCP, 145, 511
00022  *   (1998)
00023  *
00024  * CONTAINS PUBLIC FUNCTIONS: 
00025  * - integrate_3d_vl()
00026  * - integrate_destruct_3d()
00027  * - integrate_init_3d() */
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 /* MHD */
00047 
00048 /* The L/R states of primitive variables and fluxes at each cell face */
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 /* The interface magnetic fields and emfs */
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 /* MHD */
00072 
00073 /* 1D scratch vectors used by lr_states and flux functions */
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 /* primitive variables at t^{n} computed in predict step */
00079 static PrimS ***W=NULL;
00080 
00081 /* conserved & primitive variables at t^{n+1/2} computed in predict step */
00082 static ConsS ***Uhalf=NULL;
00083 static PrimS ***Whalf=NULL;
00084 
00085 /* variables needed for H-correction of Sanders et al (1998) */
00086 extern Real etah;
00087 #ifdef H_CORRECTION
00088 static Real ***eta1=NULL, ***eta2=NULL, ***eta3=NULL;
00089 #endif
00090 
00091 
00092 /*==============================================================================
00093  * PRIVATE FUNCTION PROTOTYPES: 
00094  *   integrate_emf1_corner() - upwind CT method of GS (2005) for emf1
00095  *   integrate_emf2_corner() - upwind CT method of GS (2005) for emf2 
00096  *   integrate_emf3_corner() - upwind CT method of GS (2005) for emf3
00097  *   FixCell() - apply first-order correction to one cell
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 /*=========================== PUBLIC FUNCTIONS ===============================*/
00109 /*----------------------------------------------------------------------------*/
00110 /*! \fn void integrate_3d_vl(DomainS *pD) 
00111  *  \brief 3D van Leer unsplit integrator for MHD. 
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 /* Set etah=0 so first calls to flux functions do not use H-correction */
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 /* MHD */
00167       }
00168     }
00169   }
00170 
00171 /*=== STEP 1: Compute first-order fluxes at t^{n} in x1-direction ============*/
00172 /* No source terms are needed since there is no temporal evolution */
00173 
00174 /*--- Step 1a ------------------------------------------------------------------
00175  * Load 1D vector of primitive variables;
00176  * W1d = (d, V1, V2, V3, P, B2c, B3c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 1b ------------------------------------------------------------------
00201  * Compute first-order L/R states */
00202     for (i=il; i<=ie+nghost; i++) {
00203       Wl[i] = W1d[i-1];
00204       Wr[i] = W1d[i  ];
00205 
00206 /* Compute U from W in case Pfloor used in Cons1D_to_Prim1D */
00207       Ul[i] = Prim1D_to_Cons1D(&Wl[i], &Bxc[i-1]);
00208       Ur[i] = Prim1D_to_Cons1D(&Wr[i], &Bxc[i  ]);
00209     }
00210 
00211 /*--- Step 1c ------------------------------------------------------------------
00212  * No source terms needed.
00213  */
00214 
00215 /*--- Step 1d ------------------------------------------------------------------
00216  * Compute flux in x1-direction */
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 /*=== STEP 2: Compute first-order fluxes at t^{n} in x2-direction ============*/
00228 /* No source terms are needed since there is no temporal evolution */
00229 
00230 /*--- Step 2a ------------------------------------------------------------------
00231  * Load 1D vector of primitive variables;
00232  * W1d = (d, V2, V3, V1, P, B3c, B1c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 2b ------------------------------------------------------------------
00257  * Compute first-order L/R states */
00258       for (j=jl; j<=je+nghost; j++) {
00259         Wl[j] = W1d[j-1];
00260         Wr[j] = W1d[j  ];
00261 
00262 /* Compute U from W in case Pfloor used in Cons1D_to_Prim1D */
00263         Ul[j] = Prim1D_to_Cons1D(&Wl[j], &Bxc[j-1]);
00264         Ur[j] = Prim1D_to_Cons1D(&Wr[j], &Bxc[j  ]);
00265       }
00266 
00267 /*--- Step 2c ------------------------------------------------------------------
00268  * No source terms needed
00269  */
00270 
00271 /*--- Step 2d ------------------------------------------------------------------
00272  * Compute flux in x2-direction */
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 /*=== STEP 3: Compute first-order fluxes at t^{n} in x3-direction ============*/
00285 /* No source terms are needed since there is no temporal evolution */
00286 
00287 /*--- Step 3a ------------------------------------------------------------------
00288  * Load 1D vector of primitive variables;
00289  * W1d = (d, V3, V1, V2, P, B1c, B2c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 3b ------------------------------------------------------------------
00314  * Compute first-order L/R states */      
00315       for (k=kl; k<=ke+nghost; k++) { 
00316         Wl[k] = W1d[k-1];
00317         Wr[k] = W1d[k  ]; 
00318 
00319 /* Compute U from W in case Pfloor used in Cons1D_to_Prim1D */
00320         Ul[k] = Prim1D_to_Cons1D(&Wl[k], &Bxc[k-1]);
00321         Ur[k] = Prim1D_to_Cons1D(&Wr[k], &Bxc[k  ]);
00322       }
00323 
00324 /*--- Step 3c ------------------------------------------------------------------
00325  * No source terms needed. 
00326  */
00327 
00328 /*--- Step 3d ------------------------------------------------------------------
00329  * Compute flux in x1-direction */
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 /*=== STEP 4:  Update face-centered B for 0.5*dt =============================*/
00341 
00342 /*--- Step 4a ------------------------------------------------------------------
00343  * Calculate the cell centered value of emf1,2,3 at t^{n} and integrate
00344  * to corner.
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 /*--- Step 4b ------------------------------------------------------------------
00365  * Update the interface magnetic fields using CT for a half time step.
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 /*--- Step 4c ------------------------------------------------------------------
00394  * Compute cell-centered magnetic fields at half-timestep from average of
00395  * face-centered fields.
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 /* MHD */
00408 
00409 /*=== STEP 5: Update cell-centered variables to half-timestep ================*/
00410 
00411 /*--- Step 5a ------------------------------------------------------------------
00412  * Update cell-centered variables to half-timestep using x1-fluxes
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 /* BAROTROPIC */
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 /*--- Step 5b ------------------------------------------------------------------
00438  * Update cell-centered variables to half-timestep using x2-fluxes
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 /* BAROTROPIC */
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 /*--- Step 5c ------------------------------------------------------------------
00464  * Update cell-centered variables to half-timestep using x3-fluxes
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 /* BAROTROPIC */
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 /*--- Step 5d ------------------------------------------------------------------
00491  * With first-order flux correction, save predict fluxes and emf3
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 /*=== STEP 6: Add source terms to predict values at half-timestep ============*/
00564 
00565 /*--- Step 6a ------------------------------------------------------------------
00566  * Add source terms from a static gravitational potential for 0.5*dt to predict
00567  * step.  To improve conservation of total energy, we average the energy
00568  * source term computed at cell faces.
00569  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
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 /*=== STEP 7: Conserved->Primitive variable inversion at t^{n+1/2} ===========*/
00608         
00609 /* Invert conserved variables at t^{n+1/2} to primitive variables. With FOFC, 
00610  * check if cell-centered d < 0, P< 0, or v^2 > 1. With Entropy fix, correct
00611  * by computing new primitive state using the entropy variable, otherwise
00612  * correct by switching back to values at beginning of step, rendering update
00613  * first order in time for that cell.
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 /* USE_ENTROPY_FIX */
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 /* USE_ENTROPY_FIX */
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 /* USE_ENTROPY_FIX */
00687   }
00688 #endif 
00689 
00690 /*=== STEP 8: Compute second-order L/R x1-interface states ===================*/
00691 
00692 /*--- Step 8a ------------------------------------------------------------------
00693  * Load 1D vector of primitve variables;
00694  * W1d = (d, V1, V2, V3, P, B2c, B3c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 8b ------------------------------------------------------------------
00718  * Compute L and R states at x1-interfaces, store in 3D array
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 /*=== STEP 9: Compute second-order L/R x2-interface states ===================*/
00746 
00747 /*--- Step 9a ------------------------------------------------------------------
00748  * Load 1D vector of conserved variables;
00749  * W1d = (d, V2, V3, V1, P, B3c, B1c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 9b ------------------------------------------------------------------
00773  * Compute L and R states at x2-interfaces, store in 3D array
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 /*=== STEP 10: Compute second-order L/R x3-interface states ==================*/
00801 
00802 /*--- Step 9a ------------------------------------------------------------------
00803  * Load 1D vector of conserved variables;
00804  * W1d = (d, V3, V1, V2, P, B1c, B2c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 9b ------------------------------------------------------------------
00828  * Compute L and R states at x3-interfaces, store in 3D array
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 /*=== STEP 11: Compute 3D x1-Flux, x2-Flux, x3-Flux ==========================*/
00855 
00856 /*--- Step 11a -----------------------------------------------------------------
00857  * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
00858  *  al. (1998)) for H-correction, if needed.
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 /* H_CORRECTION */
00907 
00908 /*--- Step 11b -----------------------------------------------------------------
00909  * Compute second-order fluxes in x1-direction
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 /* H_CORRECTION */
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 /* revert to predictor flux if this flux Nan'ed */
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 /*--- Step 11c -----------------------------------------------------------------
00978  * Compute second-order fluxes in x2-direction
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 /* H_CORRECTION */
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 /* revert to predictor flux if this flux NaN'ed */
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 /*--- Step 11d -----------------------------------------------------------------
01048  * Compute second-order fluxes in x3-direction
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 /* H_CORRECTION */
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 /* revert to predictor flux if this flux NaN'ed */
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 /*=== STEP 12: Update face-centered B for a full timestep ====================*/
01117         
01118 /*--- Step 12a -----------------------------------------------------------------
01119  * Calculate the cell centered value of emf1,2,3 at the half-time-step.
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 /*--- Step 12b -----------------------------------------------------------------
01138  * Integrate emf*^{n+1/2} to the grid cell corners
01139  */
01140 
01141 #ifdef MHD
01142   integrate_emf1_corner(pG);
01143   integrate_emf2_corner(pG);
01144   integrate_emf3_corner(pG);
01145 
01146 /*--- Step 12c -----------------------------------------------------------------
01147  * Update the interface magnetic fields using CT for a full time step.
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 /*--- Step 12d -----------------------------------------------------------------
01179  * Set cell centered magnetic fields to average of updated face centered fields.
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 /*=== STEP 13: Add source terms for a full timestep using n+1/2 states =======*/
01194        
01195 /*--- Step 13a -----------------------------------------------------------------
01196  * Add gravitational source terms due to a Static Potential
01197  * To improve conservation of total energy, we average the energy
01198  * source term computed at cell faces.
01199  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
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 /*=== STEP 14: Update cell-centered values for a full timestep ===============*/
01238 
01239 /*--- Step 14a -----------------------------------------------------------------
01240  * Update cell-centered variables in pG using 3D x1-Fluxes
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 /* BAROTROPIC */
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 /*--- Step 14b -----------------------------------------------------------------
01266  * Update cell-centered variables in pG using 3D x2-Fluxes
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 /* BAROTROPIC */
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 /*--- Step 14c -----------------------------------------------------------------
01292  * Update cell-centered variables in pG using 3D x3-Fluxes
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 /* BAROTROPIC */
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 /*=== STEP 15: First-order flux correction ===================================*/
01319 /*--- Step 15a -----------------------------------------------------------------
01320  * If cell-centered d or P have gone negative, or if v^2 > 1 in SR, correct
01321  * by using 1st order predictor fluxes */
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 /*--- Step 15b -----------------------------------------------------------------
01364  * In SR the first-order flux correction can fail to fix an unphysical state.
01365  * We must fix these cells in order to avoid NaN's at the next timestep,
01366  * particuarly if v^2 > 1. We have 2 approaches; firstly, we use the entropy
01367  * equation (which we have not applied a 1st order flux correction to) to
01368  * calculate the pressure and the Lorentz factor of the gas. If this produces
01369  * and unphysical state, then we floor the pressure and iterate on v^2 until
01370  * v^2 < 1. Possibly could improved by averaging density and pressure from
01371  * adjacent cells and then calculating pressure. */
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 /* USE_ENTROPY_FIX */
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 /* MHD */
01442 #endif /* FIRST_ORDER_FLUX_CORRECTION */
01443 
01444 
01445 #ifdef STATIC_MESH_REFINEMENT
01446 /*=== STEP 16: With SMR, store fluxes at fine/coarse boundaries ==============*/
01447 
01448 /*--- Step 16a -----------------------------------------------------------------
01449  * store fluxes at boundaries of child grids.
01450  */
01451 
01452   for (ncg=0; ncg<pG->NCGrid; ncg++) {
01453 
01454 /* x1-boundaries of child Grids (interior to THIS Grid) */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01498       }
01499     }
01500 
01501 /* x2-boundaries of child Grids (interior to THIS Grid) */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01545       }
01546     }
01547 
01548 /* x3-boundaries of child Grids (interior to THIS Grid) */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01592       }
01593     }
01594   } /* end loop over child Grids */
01595 
01596 /*--- Step 16b -----------------------------------------------------------------
01597  * store fluxes at boundaries with parent grids.
01598  */
01599 
01600   for (npg=0; npg<pG->NPGrid; npg++) {
01601 
01602 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01646       }
01647     }
01648 
01649 /* x2-boundaries of parent Grids (at boundaries of THIS Grid)  */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01693       }
01694     }
01695 
01696 /* x3-boundaries of parent Grids (at boundaries of THIS Grid)  */
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 /* BAROTROPIC */
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 /* MHD */
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 /* MHD */
01740       }
01741     }
01742   }
01743 
01744 #endif /* STATIC_MESH_REFINEMENT */
01745 
01746   return;
01747 }
01748 
01749 
01750 /*----------------------------------------------------------------------------*/
01751 /*! \fn void integrate_init_3d(MeshS *pM)
01752  *  \brief Allocate temporary integration arrays */
01753 void integrate_init_3d(MeshS *pM)
01754 {
01755   int nmax,size1=0,size2=0,size3=0,nl,nd;
01756 
01757 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
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 /* MHD */
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 /* H_CORRECTION */
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 /* MHD */
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 /* FIRST_ORDER_FLUX_CORRECTION */
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 /*! \fn void integrate_destruct_3d(void)
01893  *  \brief Free temporary integration arrays */
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 /* MHD */
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 /* H_CORRECTION */
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 /* MHD */
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 /* FIRST_ORDER_FLUX_CORRECTION */
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 /*=========================== PRIVATE FUNCTIONS ==============================*/
01965 
01966 /*----------------------------------------------------------------------------*/
01967 
01968 #ifdef MHD
01969 /*! \fn static void integrate_emf1_corner(const GridS *pG)
01970  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
01971  *
01972  *  Note:
01973  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
01974  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
01975  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
01976  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
01977  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
01978  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX
01979  *- first_order_correction() - Added by N. Lemaster to run supersonic turbulence
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 /* NOTE: The x2-Flux of By is -E1. */
01996 /*       The x3-Flux of Bz is +E1. */
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 /*! \fn static void integrate_emf2_corner(const GridS *pG)
02045  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
02046  *
02047  *  Note:
02048  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
02049  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
02050  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
02051  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
02052  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
02053  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX
02054  *- first_order_correction() - Added by N. Lemaster to run supersonic turbulence
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 /* NOTE: The x1-Flux of Bz is +E2. */
02070 /*       The x3-Flux of By is -E2. */
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 /*! \fn static void integrate_emf3_corner(const GridS *pG)
02120  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
02121  *
02122  *  Note:
02123  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
02124  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
02125  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
02126  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
02127  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
02128  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX
02129  *- first_order_correction() - Added by N. Lemaster to run supersonic turbulence
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 /* NOTE: The x1-Flux of By is -E3. */
02144 /*       The x2-Flux of Bx is +E3. */
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 /* MHD */
02191 
02192 #ifdef FIRST_ORDER_FLUX_CORRECTION
02193 /*----------------------------------------------------------------------------*/
02194 /*! \fn static void FixCell(GridS *pG, Int3Vect ix)
02195  *  \brief Uses first order fluxes to fix negative d,P or superluminal v
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 /* Compute difference of predictor and corrector fluxes at cell faces */
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 /* BAROTROPIC */
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 /* MHD */
02268 
02269 /* Use flux differences to correct bad cell-centered quantities */
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 /* BAROTROPIC */
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 /* BAROTROPIC */
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 /* BAROTROPIC */
02294 
02295 /* Use flux differences to correct cell-centered values at i-1 and i+1 */
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 /* BAROTROPIC */
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 /* BAROTROPIC */
02315   }
02316 
02317 /* Use flux differences to correct cell-centered values at j-1 and j+1 */
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 /* BAROTROPIC */
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 /* BAROTROPIC */
02337   }
02338 
02339 /* Use flux differences to correct cell-centered values at k-1 and k+1 */
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 /* BAROTROPIC */
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 /* BAROTROPIC */
02359   }
02360 
02361 /* Use emf differences to correct face-centered fields on edges of bad cell */
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 /* Use emf differences to correct face-centered fields around bad cell */
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 /* Compute new cell-centered fields */
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 /* MHD */
02425 
02426 /* With SMR, replace higher-order fluxes with predict fluxes in case they are
02427  * used at fine/coarse grid boundaries */
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 /* STATIC_MESH_REFINEMENT */
02433 
02434 }
02435 #endif /* FIRST_ORDER_FLUX_CORRECTION */
02436 
02437 #endif /* VL_INTEGRATOR */
02438 
02439 #endif /* SPECIAL_RELATIVITY */
integrators/integrate_3d_vl_sr.c
