Athena: integrators/integrate_2d_vl

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_2d_vl_sr.c
00004  *  \brief Integrate MHD equations using 2D version of the directionally
00005  *   unsplit MUSCL-Hancock (VL) integrator.
00006  *
00007  * PURPOSE: Integrate MHD equations using 2D 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  -- interface magnetic field
00011  *   Also adds gravitational source terms, self-gravity, and the H-correction
00012  *   of Sanders et al.
00013  *
00014  * REFERENCE: 
00015  * - J.M Stone & T.A. Gardiner, "A simple, unsplit Godunov method
00016  *   for multidimensional MHD", NewA 14, 139 (2009)
00017  *
00018  * - R. Sanders, E. Morano, & M.-C. Druguet, "Multidimensional dissipation for
00019  *   upwind schemes: stability and applications to gas dynamics", JCP, 145, 511
00020  *   (1998)
00021  *
00022  * CONTAINS PUBLIC FUNCTIONS: 
00023  * - integrate_2d_vl()
00024  * - integrate_destruct_2d()
00025  * - integrate_init_2d() */
00026 /*============================================================================*/
00027 
00028 #include <math.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include "../defs.h"
00033 #include "../athena.h"
00034 #include "../globals.h"
00035 #include "prototypes.h"
00036 #include "../prototypes.h"
00037 
00038 #ifdef SPECIAL_RELATIVITY
00039 
00040 #if defined(VL_INTEGRATOR) && defined(CARTESIAN)
00041 
00042 #ifdef MHD
00043 #define USE_ENTROPY_FIX
00044 #endif /* MHD */
00045 
00046 /* The L/R states of primitive variables and fluxes at each cell face */
00047 static Prim1DS **Wl_x1Face=NULL, **Wr_x1Face=NULL;
00048 static Prim1DS **Wl_x2Face=NULL, **Wr_x2Face=NULL;
00049 static Cons1DS **x1Flux=NULL, **x2Flux=NULL;
00050 #ifdef FIRST_ORDER_FLUX_CORRECTION
00051 static Cons1DS **x1FluxP=NULL, **x2FluxP=NULL;
00052 static Real **emf3P=NULL;
00053 #endif
00054 #ifdef USE_ENTROPY_FIX
00055 static Real **S=NULL,**Shalf=NULL;
00056 static Real **x1FluxS=NULL,**x2FluxS=NULL;
00057 #ifdef FIRST_ORDER_FLUX_CORRECTION
00058 static Real **x1FluxSP=NULL,**x2FluxSP=NULL;
00059 #endif
00060 #endif
00061 
00062 /* The interface magnetic fields and emfs */
00063 #ifdef MHD
00064 static Real **B1_x1Face=NULL, **B2_x2Face=NULL;
00065 static Real **emf3=NULL, **emf3_cc=NULL;
00066 #endif /* MHD */
00067 
00068 /* 1D scratch vectors used by lr_states and flux functions */
00069 static Real *Bxc=NULL, *Bxi=NULL;
00070 static Prim1DS *W1d=NULL, *Wl=NULL, *Wr=NULL;
00071 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00072 
00073 /* primitive variables at t^{n} */
00074 static PrimS **W=NULL;
00075 
00076 /* conserved and primitive variables at t^{n+1/2} computed in predict step */
00077 static ConsS **Uhalf=NULL;
00078 static PrimS **Whalf=NULL;
00079 
00080 /* variables needed for H-correction of Sanders et al (1998) */
00081 extern Real etah;
00082 #ifdef H_CORRECTION
00083 static Real **eta1=NULL, **eta2=NULL;
00084 #endif
00085 
00086 /*==============================================================================
00087  * PRIVATE FUNCTION PROTOTYPES: 
00088  *   integrate_emf3_corner() - the upwind CT method of Gardiner & Stone (2005) 
00089  *   FixCell() - apply first-order correction to one cell
00090  *============================================================================*/
00091 #ifdef MHD
00092 static void integrate_emf3_corner(GridS *pG);
00093 #endif
00094 #ifdef FIRST_ORDER_FLUX_CORRECTION
00095 static void FixCell(GridS *pG, Int3Vect);
00096 #endif
00097 
00098 /*=========================== PUBLIC FUNCTIONS ===============================*/
00099 /*----------------------------------------------------------------------------*/
00100 /*! \fn void integrate_2d_vl(DomainS *pD)
00101  *  \brief van Leer unsplit integrator in 2D. 
00102  *
00103  *   The numbering of steps follows the numbering in the 3D version.
00104  *   NOT ALL STEPS ARE NEEDED IN 2D.
00105  */
00106 void integrate_2d_vl(DomainS *pD)
00107 {
00108   GridS *pG=(pD->Grid);
00109   Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2;
00110   Real hdtodx1 = 0.5*dtodx1, hdtodx2 = 0.5*dtodx2;
00111   int i, is = pG->is, ie = pG->ie;
00112   int j, js = pG->js, je = pG->je;
00113   int ks = pG->ks;
00114   int cart_x1 = 1, cart_x2 = 2, cart_x3 = 3;
00115   Real pmhalf,pmnew;
00116   Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,Bx;
00117 #if (NSCALARS > 0)
00118   int n;
00119 #endif
00120 #ifdef SELF_GRAVITY
00121   Real gxl,gxr,gyl,gyr,flx_m1l,flx_m1r,flx_m2l,flx_m2r;
00122 #endif
00123 #ifdef H_CORRECTION
00124   Real cfr,cfl,lambdar,lambdal;
00125 #endif
00126 #ifdef STATIC_MESH_REFINEMENT
00127   int ncg,npg,dim;
00128   int ii,ics,ice,jj,jcs,jce,ips,ipe,jps,jpe;
00129 #endif
00130 #ifdef FIRST_ORDER_FLUX_CORRECTION
00131   ConsS Ucheck;
00132   PrimS Wcheck;
00133   int flag_cell=0,negd=0,negP=0,superl=0,NaNFlux=0;
00134   int entropy,final,fail;
00135   Real Vsq;
00136   Int3Vect BadCell;
00137 #endif
00138   int il=is-(nghost-1), iu=ie+(nghost-1);
00139   int jl=js-(nghost-1), ju=je+(nghost-1);
00140 
00141 /* Set etah=0 so first calls to flux functions do not use H-correction */
00142   etah = 0.0;
00143 
00144   for (j=js-nghost; j<=je+nghost; j++) {
00145     for (i=is-nghost; i<=ie+nghost; i++) {
00146       Uhalf[j][i] = pG->U[ks][j][i];
00147       W[j][i] = Cons_to_Prim(&(pG->U[ks][j][i]));
00148 #ifdef MHD
00149       B1_x1Face[j][i] = pG->B1i[ks][j][i]; 
00150       B2_x2Face[j][i] = pG->B2i[ks][j][i]; 
00151 #endif /* MHD */
00152 #ifdef USE_ENTROPY_FIX
00153       S[j][i] = W[j][i].P * pow(W[j][i].d,1.0-Gamma);
00154       S[j][i]*= pG->U[ks][j][i].d / W[j][i].d;
00155       Shalf[j][i] = S[j][i];
00156 #endif
00157     }
00158   }
00159 
00160 /*=== STEP 1: Compute first-order fluxes at t^{n} in x1-direction ============*/
00161 /* No source terms are needed since there is no temporal evolution */
00162 
00163 /*--- Step 1a ------------------------------------------------------------------
00164  * Load 1D vector of conserved variables;
00165  * W1d = (d, V1, V2, V3, P, B2c, B3c, s[n])
00166  */
00167 
00168   for (j=js-nghost; j<=je+nghost; j++) {
00169     for (i=is-nghost; i<=ie+nghost; i++) {
00170       W1d[i].d  = W[j][i].d;
00171       W1d[i].Vx = W[j][i].V1;
00172       W1d[i].Vy = W[j][i].V2;
00173       W1d[i].Vz = W[j][i].V3;
00174       W1d[i].P  = W[j][i].P;
00175 #ifdef MHD
00176       W1d[i].By = W[j][i].B2c;
00177       W1d[i].Bz = W[j][i].B3c;
00178       Bxc[i] = W[j][i].B1c;
00179       Bxi[i] = pG->B1i[ks][j][i];
00180 #endif /* MHD */
00181 #if (NSCALARS > 0)
00182       for (n=0; n<NSCALARS; n++) W1d[i].s[n] = W[j][i].s[n];
00183 #endif
00184     }
00185 
00186 /*--- Step 1b ------------------------------------------------------------------
00187  * Compute first-order L/R states */
00188     for (i=il; i<=ie+nghost; i++) {
00189       Wl[i] = W1d[i-1];
00190       Wr[i] = W1d[i  ];
00191 
00192 /* Compute U from W in case Pfloor used in Cons1D_to_Prim1D */
00193       Ul[i] = Prim1D_to_Cons1D(&Wl[i], &Bxc[i-1]);
00194       Ur[i] = Prim1D_to_Cons1D(&Wr[i], &Bxc[i  ]);
00195     }
00196 
00197 /*--- Step 1c ------------------------------------------------------------------
00198  * No source terms needed */
00199 
00200 /*--- Step 1d ------------------------------------------------------------------
00201  * Compute flux in x1-direction */
00202 
00203     for (i=il; i<=ie+nghost; i++) {
00204       fluxes(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1Flux[j][i]);
00205 #ifdef USE_ENTROPY_FIX
00206       entropy_flux(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1FluxS[j][i]);
00207 #endif
00208     }
00209   }
00210 
00211 /*=== STEP 2: Compute first-order fluxes at t^{n} in x2-direction ============*/
00212 /* No source terms are needed since there is no temporal evolution */
00213 
00214 /*--- Step 2a ------------------------------------------------------------------
00215  * Load 1D vector of primitive variables;
00216  * W1d = (d, V2, V3, V1, P, B3c, B1c, s[n])
00217  */
00218 
00219   for (i=is-nghost; i<=ie+nghost; i++) {
00220     for (j=js-nghost; j<=je+nghost; j++) {
00221       W1d[j].d  = W[j][i].d;
00222       W1d[j].Vx = W[j][i].V2;
00223       W1d[j].Vy = W[j][i].V3;
00224       W1d[j].Vz = W[j][i].V1;
00225       W1d[j].P  = W[j][i].P;
00226 #ifdef MHD
00227       W1d[j].By = W[j][i].B3c;
00228       W1d[j].Bz = W[j][i].B1c;
00229       Bxc[j] = W[j][i].B2c;
00230       Bxi[j] = pG->B2i[ks][j][i];
00231 #endif /* MHD */
00232 #if (NSCALARS > 0)
00233       for (n=0; n<NSCALARS; n++) W1d[j].s[n] = W[j][i].s[n];
00234 #endif
00235     }
00236 
00237 /*--- Step 2b ------------------------------------------------------------------
00238  * Compute first-order L/R states */
00239 
00240     for (j=jl; j<=je+nghost; j++) {
00241       Wl[j] = W1d[j-1];
00242       Wr[j] = W1d[j  ];
00243 
00244 /* Compute U from W in case Pfloor used in Cons1D_to_Prim1D */
00245       Ul[j] = Prim1D_to_Cons1D(&Wl[j], &Bxc[j-1]);
00246       Ur[j] = Prim1D_to_Cons1D(&Wr[j], &Bxc[j  ]);
00247     }
00248 
00249 /*--- Step 2c ------------------------------------------------------------------
00250  * No source terms needed */
00251 
00252 /*--- Step 2d ------------------------------------------------------------------
00253  * Compute flux in x2-direction */
00254 
00255     for (j=jl; j<=je+nghost; j++) {
00256       fluxes(Ul[j],Ur[j],Wl[j],Wr[j],Bxi[j],&x2Flux[j][i]);
00257 #ifdef USE_ENTROPY_FIX
00258       entropy_flux(Ul[j],Ur[j],Wl[j],Wr[j],Bxi[j],&x2FluxS[j][i]);
00259 #endif
00260     }
00261   }
00262 
00263 /*=== STEP 3: Not needed in 2D ===*/
00264 
00265 /*=== STEP 4:  Update face-centered B for 0.5*dt =============================*/
00266 
00267 /*--- Step 4a ------------------------------------------------------------------
00268  * Calculate the cell centered value of emf1,2,3 at t^{n} and integrate
00269  * to corner.
00270  */
00271 
00272 #ifdef MHD
00273   for (j=js-nghost; j<=je+nghost; j++) {
00274     for (i=is-nghost; i<=ie+nghost; i++) {
00275       emf3_cc[j][i] = (W[j][i].B1c*W[j][i].V2 - W[j][i].B2c*W[j][i].V1);
00276     }
00277   }
00278   integrate_emf3_corner(pG);
00279 
00280 /*--- Step 4b ------------------------------------------------------------------
00281  * Update the interface magnetic fields using CT for a half time step.
00282  */
00283 
00284   for (j=jl; j<=ju; j++) {
00285     for (i=il; i<=iu; i++) {
00286       B1_x1Face[j][i] -= hdtodx2*(emf3[j+1][i  ] - emf3[j][i]);
00287       B2_x2Face[j][i] += hdtodx1*(emf3[j  ][i+1] - emf3[j][i]);
00288     }
00289     B1_x1Face[j][iu+1] -= hdtodx2*(emf3[j+1][iu+1]-emf3[j][iu+1]);
00290   }
00291   for (i=il; i<=iu; i++) {
00292     B2_x2Face[ju+1][i] += hdtodx1*(emf3[ju+1][i+1]-emf3[ju+1][i]);
00293   }
00294 
00295 /*--- Step 4c ------------------------------------------------------------------
00296  * Compute cell-centered magnetic fields at half-timestep from average of
00297  * face-centered fields.
00298  */
00299 
00300   for (j=jl; j<=ju; j++) {
00301     for (i=il; i<=iu; i++) {
00302       Uhalf[j][i].B1c = 0.5*(B1_x1Face[j][i] + B1_x1Face[j][i+1]);
00303       Uhalf[j][i].B2c = 0.5*(B2_x2Face[j][i] + B2_x2Face[j+1][i]);
00304     }
00305   }
00306 #endif /* MHD */
00307 
00308 /*=== STEP 5: Update cell-centered variables to half-timestep ================*/
00309 
00310 /*--- Step 5a ------------------------------------------------------------------
00311  * Update cell-centered variables (including B3c) to half-timestep with x1Flux
00312  */
00313 
00314   for (j=jl; j<=ju; j++) {
00315     for (i=il; i<=iu; i++) {
00316       Uhalf[j][i].d   -= hdtodx1*(x1Flux[j][i+1].d  - x1Flux[j][i].d );
00317       Uhalf[j][i].M1  -= hdtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx);
00318       Uhalf[j][i].M2  -= hdtodx1*(x1Flux[j][i+1].My - x1Flux[j][i].My);
00319       Uhalf[j][i].M3  -= hdtodx1*(x1Flux[j][i+1].Mz - x1Flux[j][i].Mz);
00320 #ifndef BAROTROPIC
00321       Uhalf[j][i].E   -= hdtodx1*(x1Flux[j][i+1].E  - x1Flux[j][i].E );
00322 #endif /* BAROTROPIC */
00323 #ifdef MHD
00324       Uhalf[j][i].B3c -= hdtodx1*(x1Flux[j][i+1].Bz - x1Flux[j][i].Bz);
00325 #endif /* MHD */
00326 #if (NSCALARS > 0)
00327       for (n=0; n<NSCALARS; n++)
00328         Uhalf[j][i].s[n] -= hdtodx1*(x1Flux[j][i+1].s[n] - x1Flux[j][i].s[n]);
00329 #endif
00330 #ifdef USE_ENTROPY_FIX
00331       Shalf[j][i]   -= hdtodx1*(x1FluxS[j][i+1]  - x1FluxS[j][i]);
00332 #endif
00333     }
00334   }
00335 
00336 /*--- Step 5b ------------------------------------------------------------------
00337  * Update cell-centered variables (including B3c) to half-timestep with x2Flux
00338  */
00339 
00340   for (j=jl; j<=ju; j++) {
00341     for (i=il; i<=iu; i++) {
00342       Uhalf[j][i].d   -= hdtodx2*(x2Flux[j+1][i].d  - x2Flux[j][i].d );
00343       Uhalf[j][i].M1  -= hdtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
00344       Uhalf[j][i].M2  -= hdtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
00345       Uhalf[j][i].M3  -= hdtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
00346 #ifndef BAROTROPIC
00347       Uhalf[j][i].E   -= hdtodx2*(x2Flux[j+1][i].E  - x2Flux[j][i].E );
00348 #endif /* BAROTROPIC */
00349 #ifdef MHD
00350       Uhalf[j][i].B3c -= hdtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
00351 #endif /* MHD */
00352 #if (NSCALARS > 0)
00353       for (n=0; n<NSCALARS; n++)
00354         Uhalf[j][i].s[n] -= hdtodx2*(x2Flux[j+1][i].s[n] - x2Flux[j][i].s[n]);
00355 #endif
00356 #ifdef USE_ENTROPY_FIX
00357       Shalf[j][i]   -= hdtodx2*(x2FluxS[j+1][i]  - x2FluxS[j][i] );
00358 #endif
00359     }
00360   }
00361 
00362 #ifdef FIRST_ORDER_FLUX_CORRECTION
00363 /*--- Step 5d ------------------------------------------------------------------
00364  * With first-order flux correction, check predict fluxes & emf3 for NaN's
00365  * and save predict fluxes & emf3 for first-order flux correction
00366  */
00367 
00368   NaNFlux = 0;
00369   for (j=js; j<=je+1; j++) {
00370     for (i=is; i<=ie+1; i++) {
00371       if ((x1Flux[j][i].d  != x1Flux[j][i].d)  ||
00372 #ifndef BAROTROPIC
00373           (x1Flux[j][i].E  != x1Flux[j][i].E)  ||
00374 #endif
00375 #ifdef MHD
00376           (x1Flux[j][i].By != x1Flux[j][i].By) ||
00377           (x1Flux[j][i].Bz != x1Flux[j][i].Bz) ||
00378 #endif
00379           (x1Flux[j][i].Mx != x1Flux[j][i].Mx) ||
00380           (x1Flux[j][i].My != x1Flux[j][i].My) ||
00381           (x1Flux[j][i].Mz != x1Flux[j][i].Mz)) {
00382         NaNFlux++;
00383       }          
00384       if ((x2Flux[j][i].d  != x2Flux[j][i].d)  ||
00385 #ifndef BAROTROPIC
00386           (x2Flux[j][i].E  != x2Flux[j][i].E)  ||
00387 #endif
00388 #ifdef MHD
00389           (x2Flux[j][i].By != x2Flux[j][i].By) ||
00390           (x2Flux[j][i].Bz != x2Flux[j][i].Bz) ||
00391 #endif
00392           (x2Flux[j][i].Mx != x2Flux[j][i].Mx) ||
00393           (x2Flux[j][i].My != x2Flux[j][i].My) ||
00394           (x2Flux[j][i].Mz != x2Flux[j][i].Mz)) {
00395         NaNFlux++;
00396       }
00397       x1FluxP[j][i] = x1Flux[j][i];
00398       x2FluxP[j][i] = x2Flux[j][i];
00399 #ifdef MHD
00400       emf3P[j][i] = emf3[j][i];
00401 #endif
00402 #ifdef USE_ENTROPY_FIX
00403       x1FluxSP[j][i] = x1FluxS[j][i];
00404       x2FluxSP[j][i] = x2FluxS[j][i];
00405 #endif
00406     }
00407   }
00408   if (NaNFlux != 0) {
00409     ath_error("[Step5d] %i first-order fluxes are NaN!\n",NaNFlux);
00410   }
00411 #endif
00412 
00413 /*=== STEP 6: Add source terms to predict values at half-timestep ============*/
00414 
00415 /*--- Step 6a ------------------------------------------------------------------
00416  * Add source terms from a static gravitational potential for 0.5*dt to predict
00417  * step.  To improve conservation of total energy, we average the energy
00418  * source term computed at cell faces.
00419  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00420  */
00421 
00422   if (StaticGravPot != NULL){
00423     for (j=jl; j<=ju; j++) {
00424       for (i=il; i<=iu; i++) {
00425         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00426         phic = (*StaticGravPot)( x1,             x2,x3);
00427         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00428         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00429 
00430         Uhalf[j][i].M1 -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00431 #ifndef BAROTROPIC
00432         Uhalf[j][i].E -= hdtodx1*(x1Flux[j][i  ].d*(phic - phil)
00433                            + x1Flux[j][i+1].d*(phir - phic));
00434 #endif
00435         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00436         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00437 
00438         Uhalf[j][i].M2 -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00439 #ifndef BAROTROPIC
00440         Uhalf[j][i].E -= hdtodx2*(x2Flux[j  ][i].d*(phic - phil)
00441                                 + x2Flux[j+1][i].d*(phir - phic));
00442 #endif
00443       }
00444     }
00445   }
00446 
00447 
00448 /*=== STEP 7: Conserved->Primitive variable inversion at t^{n+1/2} ===========*/
00449         
00450 /* Invert conserved variables at t^{n+1/2} to primitive variables. With FOFC, 
00451  * check if cell-centered d < 0, P< 0, or v^2 > 1. With Entropy fix, correct
00452  * by computing new primitive state using the entropy variable, otherwise
00453  * correct by switching back to values at beginning of step, rendering update
00454  * first order in time for that cell.
00455  */
00456 
00457 #ifdef FIRST_ORDER_FLUX_CORRECTION        
00458   negd = 0;
00459   negP = 0;
00460   superl = 0;
00461   flag_cell = 0;
00462 #endif
00463   for (j=js-nghost; j<=je+nghost; j++) {
00464     for (i=is-nghost; i<=ie+nghost; i++) {
00465       Whalf[j][i] = check_Prim(&(Uhalf[j][i]));
00466 #ifdef FIRST_ORDER_FLUX_CORRECTION   
00467       if (Whalf[j][i].d < 0.0) {
00468         flag_cell = 1;
00469         negd++;
00470       }
00471       if (Whalf[j][i].P < 0.0) {
00472         flag_cell = 1;
00473         negP++;
00474       }
00475       Vsq = SQR(Whalf[j][i].V1) + SQR(Whalf[j][i].V2) + SQR(Whalf[j][i].V3);
00476       if (Vsq > 1.0) {
00477         flag_cell = 1;
00478         superl++;
00479       }
00480       if (flag_cell != 0) {
00481 #ifdef USE_ENTROPY_FIX
00482         Wcheck = entropy_fix (&(Uhalf[j][i]),&(Shalf[j][i]));
00483         Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00484         if (Wcheck.d > 0.0 && Wcheck.P > 0.0 && Vsq < 1.0){
00485           entropy++;
00486           Whalf[j][i].d = Wcheck.d;
00487           Whalf[j][i].P = Wcheck.P;
00488           Whalf[j][i].V1 = Wcheck.V1;
00489           Whalf[j][i].V2 = Wcheck.V2;
00490           Whalf[j][i].V3 = Wcheck.V3;
00491           flag_cell=0;
00492         } else {
00493 #endif /* USE_ENTROPY_FIX */
00494           Whalf[j][i].d = W[j][i].d;
00495           Whalf[j][i].V1 = W[j][i].V1;
00496           Whalf[j][i].V2 = W[j][i].V2;
00497           Whalf[j][i].V3 = W[j][i].V3;
00498           Whalf[j][i].P = W[j][i].P;
00499           flag_cell=0;
00500 #ifdef USE_ENTROPY_FIX
00501         }
00502 #endif /* USE_ENTROPY_FIX */
00503       }
00504 #endif
00505     }
00506   }
00507 
00508 #ifdef FIRST_ORDER_FLUX_CORRECTION  
00509   if (negd > 0 || negP > 0 || superl > 0){
00510     printf("[Step7]: %i cells had d<0; %i cells had P<0\n",negd,negP);
00511     printf("[Step7]: %i cells had v>1 at t_half\n",superl);
00512 #ifdef USE_ENTROPY_FIX
00513     printf("[Step7]: %i cells fixed using entropy at t_half\n",entropy);
00514 #endif /* USE_ENTROPY_FIX */
00515   }
00516 #endif 
00517 
00518 /*=== STEP 8: Compute second-order L/R x1-interface states ===================*/
00519 
00520 /*--- Step 8a ------------------------------------------------------------------
00521  * Load 1D vector of primitive variables;
00522  * W1d = (d, V1, V2, V3, P, B2c, B3c, s[n])
00523  */
00524 
00525   for (j=js-1; j<=je+1; j++) {
00526     for (i=il; i<=iu; i++) {
00527       W1d[i].d  = Whalf[j][i].d;
00528       W1d[i].Vx = Whalf[j][i].V1;
00529       W1d[i].Vy = Whalf[j][i].V2;
00530       W1d[i].Vz = Whalf[j][i].V3;
00531 #ifndef BAROTROPIC
00532       W1d[i].P  = Whalf[j][i].P;
00533 #endif /* BAROTROPIC */
00534 #ifdef MHD
00535       W1d[i].By = Whalf[j][i].B2c;
00536       W1d[i].Bz = Whalf[j][i].B3c;
00537       Bxc[i] = Whalf[j][i].B1c;
00538 #endif /* MHD */
00539 #if (NSCALARS > 0)
00540       for (n=0; n<NSCALARS; n++) W1d[i].s[n] = Whalf[j][i].s[n];
00541 #endif /* NSCALARS */
00542     }
00543 
00544 /*--- Step 8b ------------------------------------------------------------------
00545  * Compute L/R states on x1-interfaces, store into arrays
00546  */
00547 
00548     lr_states(pG,W1d,Bxc,pG->dt,pG->dx1,is,ie,Wl,Wr,cart_x1);
00549 
00550 #ifdef FIRST_ORDER_FLUX_CORRECTION
00551     for (i=il; i<=iu; i++) {
00552       Vsq = SQR(Wl[i].Vx) + SQR(Wl[i].Vy) + SQR(Wl[i].Vz);
00553       if (Vsq > 1.0){
00554         Wl[i] = W1d[i];
00555         Wr[i] = W1d[i];
00556       }
00557       Vsq = SQR(Wr[i].Vx) + SQR(Wr[i].Vy) + SQR(Wr[i].Vz);
00558       if (Vsq > 1.0){
00559         Wl[i] = W1d[i];
00560         Wr[i] = W1d[i];
00561       }
00562     }
00563 #endif
00564 
00565     for (i=is; i<=ie+1; i++) {
00566       Wl_x1Face[j][i] = Wl[i];
00567       Wr_x1Face[j][i] = Wr[i];
00568     }
00569   }
00570 
00571 
00572 
00573 /*=== STEP 9: Compute second-order L/R x2-interface states ===================*/
00574 
00575 /*--- Step 9a ------------------------------------------------------------------
00576  * Load 1D vector of primitive variables;
00577  * W1d = (d, V2, V3, V1, P, B3c, B1c, s[n])
00578  */
00579 
00580   for (i=is-1; i<=ie+1; i++) {
00581     for (j=jl; j<=ju; j++) {
00582       W1d[j].d  = Whalf[j][i].d;
00583       W1d[j].Vx = Whalf[j][i].V2;
00584       W1d[j].Vy = Whalf[j][i].V3;
00585       W1d[j].Vz = Whalf[j][i].V1;
00586 #ifndef BAROTROPIC
00587       W1d[j].P  = Whalf[j][i].P;
00588 #endif /* BAROTROPIC */
00589 #ifdef MHD
00590       W1d[j].By = Whalf[j][i].B3c;
00591       W1d[j].Bz = Whalf[j][i].B1c;
00592       Bxc[j] = Whalf[j][i].B2c;
00593 #endif /* MHD */
00594 #if (NSCALARS > 0)
00595       for (n=0; n<NSCALARS; n++) W1d[j].s[n] = Whalf[j][i].s[n];
00596 #endif /* NSCALARS */
00597     }
00598 
00599 /*--- Step 9b ------------------------------------------------------------------
00600  * Compute L/R states on x1-interfaces, store into arrays
00601  */
00602 
00603     lr_states(pG,W1d,Bxc,pG->dt,pG->dx2,js,je,Wl,Wr,cart_x2);
00604 
00605 #ifdef FIRST_ORDER_FLUX_CORRECTION
00606     for (j=jl; j<=ju; j++) {
00607       Vsq = SQR(Wl[j].Vx) + SQR(Wl[j].Vy) + SQR(Wl[j].Vz);
00608       if (Vsq > 1.0){
00609         Wl[j] = W1d[j];
00610         Wr[j] = W1d[j];
00611       }
00612       Vsq = SQR(Wr[j].Vx) + SQR(Wr[j].Vy) + SQR(Wr[j].Vz);
00613       if (Vsq > 1.0){
00614         Wl[j] = W1d[j];
00615         Wr[j] = W1d[j];
00616       }
00617     }
00618 #endif
00619 
00620     for (j=js; j<=je+1; j++) {
00621       Wl_x2Face[j][i] = Wl[j];
00622       Wr_x2Face[j][i] = Wr[j];
00623     }
00624   }
00625 
00626 /*=== STEP 10: Not needed in 2D ===*/
00627 
00628 /*=== STEP 11: Compute 2D x1-Flux, x2-Flux, ==================================*/
00629 
00630 /*--- Step 11a -----------------------------------------------------------------
00631  * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
00632  *  al. (1998)) for H-correction, if needed.
00633  */
00634 
00635 #ifdef H_CORRECTION
00636   for (j=js-1; j<=je+1; j++) {
00637     for (i=is-1; i<=iu; i++) {
00638 #ifdef MHD
00639       Bx = B1_x1Face[j][i];
00640 #endif
00641       cfr = cfast(&(Ur_x1Face[j][i]),&Bx);
00642       cfl = cfast(&(Ul_x1Face[j][i]),&Bx);
00643       lambdar = Ur_x1Face[j][i].Mx/Ur_x1Face[j][i].d + cfr;
00644       lambdal = Ul_x1Face[j][i].Mx/Ul_x1Face[j][i].d - cfl;
00645       eta1[j][i] = 0.5*fabs(lambdar - lambdal);
00646     }
00647   }
00648 
00649   for (j=js-1; j<=ju; j++) {
00650     for (i=is-1; i<=ie+1; i++) {
00651 #ifdef MHD
00652       Bx = B2_x2Face[j][i];
00653 #endif
00654       cfr = cfast(&(Ur_x2Face[j][i]),&Bx);
00655       cfl = cfast(&(Ul_x2Face[j][i]),&Bx);
00656       lambdar = Ur_x2Face[j][i].Mx/Ur_x2Face[j][i].d + cfr;
00657       lambdal = Ul_x2Face[j][i].Mx/Ul_x2Face[j][i].d - cfl;
00658       eta2[j][i] = 0.5*fabs(lambdar - lambdal);
00659     }
00660   }
00661 #endif /* H_CORRECTION */
00662 
00663 /*--- Step 11b -----------------------------------------------------------------
00664  * Compute second-order fluxes in x1-direction
00665  */
00666 
00667   for (j=js-1; j<=je+1; j++) {
00668     for (i=is; i<=ie+1; i++) {
00669 #ifdef H_CORRECTION
00670       etah = MAX(eta2[j][i-1],eta2[j][i]);
00671       etah = MAX(etah,eta2[j+1][i-1]);
00672       etah = MAX(etah,eta2[j+1][i  ]);
00673       etah = MAX(etah,eta1[j  ][i  ]);
00674 #endif /* H_CORRECTION */
00675 #ifdef MHD
00676       Bx = B1_x1Face[j][i];
00677 #endif
00678       Ul[i] = Prim1D_to_Cons1D(&Wl_x1Face[j][i],&Bx);
00679       Ur[i] = Prim1D_to_Cons1D(&Wr_x1Face[j][i],&Bx);
00680 
00681       fluxes(Ul[i],Ur[i],Wl_x1Face[j][i],Wr_x1Face[j][i],Bx,&x1Flux[j][i]);
00682 
00683 #ifdef USE_ENTROPY_FIX
00684       entropy_flux(Ul[i],          Ur[i],
00685                    Wl_x1Face[j][i],Wr_x1Face[j][i],
00686                    Bx,             &x1FluxS[j][i]);
00687 #endif
00688 
00689 #ifdef FIRST_ORDER_FLUX_CORRECTION
00690 /* revert to predictor flux if this flux Nan'ed */
00691       if ((x1Flux[j][i].d  != x1Flux[j][i].d)  ||
00692 #ifndef BAROTROPIC
00693           (x1Flux[j][i].E  != x1Flux[j][i].E)  ||
00694 #endif
00695 #ifdef MHD
00696           (x1Flux[j][i].By != x1Flux[j][i].By) ||
00697           (x1Flux[j][i].Bz != x1Flux[j][i].Bz) ||
00698 #endif
00699           (x1Flux[j][i].Mx != x1Flux[j][i].Mx) ||
00700           (x1Flux[j][i].My != x1Flux[j][i].My) ||
00701           (x1Flux[j][i].Mz != x1Flux[j][i].Mz)) {
00702         x1Flux[j][i] = x1FluxP[j][i];
00703         NaNFlux++;
00704       }
00705 #endif
00706 
00707     }
00708   }
00709 
00710 /*--- Step 11c -----------------------------------------------------------------
00711  * Compute second-order fluxes in x2-direction
00712  */
00713 
00714   for (j=js; j<=je+1; j++) {
00715     for (i=is-1; i<=ie+1; i++) {
00716 #ifdef H_CORRECTION
00717       etah = MAX(eta1[j-1][i],eta1[j][i]);
00718       etah = MAX(etah,eta1[j-1][i+1]);
00719       etah = MAX(etah,eta1[j  ][i+1]);
00720       etah = MAX(etah,eta2[j  ][i]);
00721 #endif /* H_CORRECTION */
00722 #ifdef MHD
00723       Bx = B2_x2Face[j][i];
00724 #endif
00725       Ul[i] = Prim1D_to_Cons1D(&Wl_x2Face[j][i],&Bx);
00726       Ur[i] = Prim1D_to_Cons1D(&Wr_x2Face[j][i],&Bx);
00727 
00728       fluxes(Ul[i],Ur[i],Wl_x2Face[j][i],Wr_x2Face[j][i],Bx,&x2Flux[j][i]);
00729 
00730 #ifdef USE_ENTROPY_FIX
00731       entropy_flux(Ul[i],          Ur[i],
00732                    Wl_x2Face[j][i],Wr_x2Face[j][i],
00733                    Bx,             &x2FluxS[j][i]);
00734 #endif
00735 
00736 #ifdef FIRST_ORDER_FLUX_CORRECTION
00737 /* revert to predictor flux if this flux NaN'ed */
00738       if ((x2Flux[j][i].d  != x2Flux[j][i].d)  ||
00739 #ifndef BAROTROPIC
00740           (x2Flux[j][i].E  != x2Flux[j][i].E)  ||
00741 #endif
00742 #ifdef MHD
00743           (x2Flux[j][i].By != x2Flux[j][i].By) ||
00744           (x2Flux[j][i].Bz != x2Flux[j][i].Bz) ||
00745 #endif
00746           (x2Flux[j][i].Mx != x2Flux[j][i].Mx) ||
00747           (x2Flux[j][i].My != x2Flux[j][i].My) ||
00748           (x2Flux[j][i].Mz != x2Flux[j][i].Mz)) {
00749         x2Flux[j][i] = x2FluxP[j][i];
00750         NaNFlux++;
00751       }
00752 #endif
00753 
00754     }
00755   }
00756 
00757 #ifdef FIRST_ORDER_FLUX_CORRECTION
00758   if (NaNFlux != 0) {
00759     printf("[Step11] %i second-order fluxes replaced\n",NaNFlux);
00760     NaNFlux=0;
00761   }
00762 #endif
00763   
00764 
00765 /*=== STEP 12: Update face-centered B for a full timestep ====================*/
00766         
00767 /*--- Step 12a -----------------------------------------------------------------
00768  * Calculate the cell centered value of emf1,2,3 at the half-time-step.
00769  */
00770 
00771 #ifdef MHD
00772   for (j=js-1; j<=je+1; j++) {
00773     for (i=is-1; i<=ie+1; i++) {
00774       emf3_cc[j][i] = (Whalf[j][i].B1c*Whalf[j][i].V2 - 
00775                        Whalf[j][i].B2c*Whalf[j][i].V1);
00776     }
00777   }
00778 
00779 /*--- Step 12b -----------------------------------------------------------------
00780  * Integrate emf*^{n+1/2} to the grid cell corners
00781  */
00782 
00783   integrate_emf3_corner(pG);
00784 
00785 /*--- Step 12c -----------------------------------------------------------------
00786  * Update the interface magnetic fields using CT for a full time step.
00787  */
00788 
00789   for (j=js; j<=je; j++) {
00790     for (i=is; i<=ie; i++) {
00791       pG->B1i[ks][j][i] -= dtodx2*(emf3[j+1][i  ] - emf3[j][i]);
00792       pG->B2i[ks][j][i] += dtodx1*(emf3[j  ][i+1] - emf3[j][i]);
00793     }
00794     pG->B1i[ks][j][ie+1] -= dtodx2*(emf3[j+1][ie+1] - emf3[j][ie+1]);
00795   }
00796   for (i=is; i<=ie; i++) {
00797     pG->B2i[ks][je+1][i] += dtodx1*(emf3[je+1][i+1] - emf3[je+1][i]);
00798   }
00799 
00800 /*--- Step 12d -----------------------------------------------------------------
00801  * Set cell centered magnetic fields to average of updated face centered fields.
00802  */
00803 
00804   for (j=js; j<=je; j++) {
00805     for (i=is; i<=ie; i++) {
00806       pG->U[ks][j][i].B1c = 0.5*(pG->B1i[ks][j][i]+pG->B1i[ks][j][i+1]);
00807       pG->U[ks][j][i].B2c = 0.5*(pG->B2i[ks][j][i]+pG->B2i[ks][j+1][i]);
00808     }
00809   }
00810 #endif /* MHD */
00811 
00812 /*=== STEP 13: Add source terms for a full timestep using n+1/2 states =======*/
00813        
00814 /*--- Step 13a -----------------------------------------------------------------
00815  * Add gravitational source terms due to a Static Potential
00816  * To improve conservation of total energy, we average the energy
00817  * source term computed at cell faces.
00818  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
00819  */
00820 
00821   if (StaticGravPot != NULL){
00822     for (j=js; j<=je; j++) {
00823       for (i=is; i<=ie; i++) {
00824         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00825         phic = (*StaticGravPot)( x1,             x2,x3);
00826         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00827         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00828 
00829         pG->U[ks][j][i].M1 -= dtodx1*(phir-phil)*Uhalf[j][i].d;
00830 #ifndef BAROTROPIC
00831         pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i  ].d*(phic - phil)
00832                                    + x1Flux[j][i+1].d*(phir - phic));
00833 #endif
00834         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00835         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00836 
00837         pG->U[ks][j][i].M2 -= dtodx2*(phir-phil)*Uhalf[j][i].d;
00838 #ifndef BAROTROPIC
00839         pG->U[ks][j][i].E -= dtodx2*(x2Flux[j  ][i].d*(phic - phil)
00840                                    + x2Flux[j+1][i].d*(phir - phic));
00841 #endif
00842       }
00843     }
00844   }
00845 
00846 /*=== STEP 14: Update cell-centered values for a full timestep ===============*/
00847 
00848 /*--- Step 14a -----------------------------------------------------------------
00849  * Update cell-centered variables in pG (including B3c) using 2D x1-Fluxes
00850  */
00851 
00852   for (j=js; j<=je; j++) {
00853     for (i=is; i<=ie; i++) {
00854       pG->U[ks][j][i].d   -= dtodx1*(x1Flux[j][i+1].d  - x1Flux[j][i].d );
00855       pG->U[ks][j][i].M1  -= dtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx);
00856       pG->U[ks][j][i].M2  -= dtodx1*(x1Flux[j][i+1].My - x1Flux[j][i].My);
00857       pG->U[ks][j][i].M3  -= dtodx1*(x1Flux[j][i+1].Mz - x1Flux[j][i].Mz);
00858 #ifndef BAROTROPIC
00859       pG->U[ks][j][i].E   -= dtodx1*(x1Flux[j][i+1].E  - x1Flux[j][i].E );
00860 #endif /* BAROTROPIC */
00861 #ifdef MHD
00862       pG->U[ks][j][i].B3c -= dtodx1*(x1Flux[j][i+1].Bz - x1Flux[j][i].Bz);
00863 #endif /* MHD */
00864 #if (NSCALARS > 0)
00865       for (n=0; n<NSCALARS; n++)
00866         pG->U[ks][j][i].s[n] -= dtodx1*(x1Flux[j][i+1].s[n]
00867                                       - x1Flux[j][i  ].s[n]);
00868 #endif
00869 #ifdef USE_ENTROPY_FIX
00870       S[j][i] -= dtodx1*(x1FluxS[j][i+1]  - x1FluxS[j][i]);
00871 #endif
00872     }
00873   }
00874 
00875 /*--- Step 14b -----------------------------------------------------------------
00876  * Update cell-centered variables in pG (including B3c) using 2D x2-Fluxes
00877  */
00878 
00879   for (j=js; j<=je; j++) {
00880     for (i=is; i<=ie; i++) {
00881       pG->U[ks][j][i].d   -= dtodx2*(x2Flux[j+1][i].d  - x2Flux[j][i].d );
00882       pG->U[ks][j][i].M1  -= dtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
00883       pG->U[ks][j][i].M2  -= dtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
00884       pG->U[ks][j][i].M3  -= dtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
00885 #ifndef BAROTROPIC
00886       pG->U[ks][j][i].E   -= dtodx2*(x2Flux[j+1][i].E  - x2Flux[j][i].E );
00887 #endif /* BAROTROPIC */
00888 #ifdef MHD
00889       pG->U[ks][j][i].B3c -= dtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
00890 #endif /* MHD */
00891 #if (NSCALARS > 0)
00892       for (n=0; n<NSCALARS; n++)
00893         pG->U[ks][j][i].s[n] -= dtodx2*(x2Flux[j+1][i].s[n]
00894                                       - x2Flux[j  ][i].s[n]);
00895 #endif
00896 #ifdef USE_ENTROPY_FIX
00897       S[j][i] -= dtodx2*(x2FluxS[j+1][i]  - x2FluxS[j][i]);
00898 #endif
00899     }
00900   }  
00901 
00902 #ifdef FIRST_ORDER_FLUX_CORRECTION
00903 /*=== STEP 15: First-order flux correction ===================================*/
00904 
00905 /*--- Step 15a -----------------------------------------------------------------
00906  * If cell-centered d or P have gone negative, or if v^2 > 1 in SR, correct
00907  * by using 1st order predictor fluxes
00908  */
00909 
00910   for (j=js; j<=je; j++) {
00911     for (i=is; i<=ie; i++) {
00912       Wcheck = check_Prim(&(pG->U[ks][j][i]));
00913       if (Wcheck.d < 0.0) {
00914         flag_cell = 1;
00915         BadCell.i = i;
00916         BadCell.j = j;
00917         BadCell.k = ks;
00918         negd++;
00919       }
00920       if (Wcheck.P < 0.0) {
00921         flag_cell = 1;
00922         BadCell.i = i;
00923         BadCell.j = j;
00924         BadCell.k = ks;
00925         negP++;
00926       }
00927       Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00928       if (Vsq > 1.0) {
00929         flag_cell = 1;
00930         BadCell.i = i;
00931         BadCell.j = j;
00932         BadCell.k = ks;
00933         superl++;
00934       }
00935       if (flag_cell != 0) {
00936         FixCell(pG, BadCell);
00937         flag_cell=0;
00938       }
00939     }
00940   }
00941 
00942   if (negd > 0 || negP > 0 || superl > 0){
00943     printf("[Step15a]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00944     printf("[Step15a]: %i cells had v>1 at 1st order correction\n",superl);
00945   }
00946 
00947 /*--- Step 15b -----------------------------------------------------------------
00948  * In SR the first-order flux correction can fail to fix an unphysical state.
00949  * We must fix these cells in order to avoid NaN's at the next timestep,
00950  * particuarly if v^2 > 1. We have 2 approaches; firstly, we use the entropy
00951  * equation (which we have not applied a 1st order flux correction to) to
00952  * calculate the pressure and the Lorentz factor of the gas. If this produces
00953  * and unphysical state, then we floor the pressure and iterate on v^2 until
00954  * v^2 < 1. Possibly could improved by averaging density and pressure from
00955  * adjacent cells and then calculating pressure.
00956  */
00957 
00958 #ifdef MHD
00959   negd = 0;
00960   negP = 0;
00961   superl = 0;
00962   entropy = 0;
00963   final = 0;
00964   fail = 0;
00965   for (j=js; j<=je; j++) {
00966     for (i=is; i<=ie; i++) {
00967       flag_cell=0;
00968       Wcheck = check_Prim(&(pG->U[ks][j][i]));
00969       if (Wcheck.d < 0.0) {
00970         flag_cell = 1;
00971         negd++;
00972       }
00973       if (Wcheck.P < 0.0) {
00974         flag_cell = 1;
00975         negP++;
00976       }
00977       Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00978       if (Vsq > 1.0) {
00979         flag_cell = 1;
00980         superl++;
00981       }
00982 #ifdef USE_ENTROPY_FIX
00983       if (flag_cell != 0) {
00984         entropy++;
00985         Wcheck = entropy_fix (&(pG->U[ks][j][i]),&(S[j][i]));
00986         Ucheck = Prim_to_Cons(&Wcheck);
00987         Wcheck = check_Prim(&Ucheck);
00988         Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00989         if (Wcheck.d > 0.0 && Wcheck.P > 0.0 && Vsq < 1.0){
00990           pG->U[ks][j][i].d = Ucheck.d;
00991           pG->U[ks][j][i].M1 = Ucheck.M1;
00992           pG->U[ks][j][i].M2 = Ucheck.M2;
00993           pG->U[ks][j][i].M3 = Ucheck.M3;
00994           pG->U[ks][j][i].E = Ucheck.E;
00995           flag_cell = 0;
00996         }
00997       }
00998 #endif
00999       if (flag_cell != 0) {
01000         final++;
01001         Wcheck = fix_vsq (&(pG->U[ks][j][i]));
01002         Ucheck = Prim_to_Cons(&Wcheck);
01003         pG->U[ks][j][i].d = Ucheck.d;
01004         pG->U[ks][j][i].M1 = Ucheck.M1;
01005         pG->U[ks][j][i].M2 = Ucheck.M2;
01006         pG->U[ks][j][i].M3 = Ucheck.M3;
01007         pG->U[ks][j][i].E = Ucheck.E;
01008         Wcheck = check_Prim(&(pG->U[ks][j][i]));
01009         Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
01010         if (Wcheck.d < 0.0 || Wcheck.P < 0.0 || Vsq > 1.0){
01011           fail++;
01012         }  
01013       }
01014     }
01015   }
01016 
01017   if (negd > 0 || negP > 0 || superl > 0) {
01018     printf("[Step15b]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
01019     printf("[Step15b]: %i cells had v>1 after 1st order correction\n",superl);
01020     printf("[Step15b]: %i cells required an entropy fix\n",entropy);
01021     printf("[Step15b]: %i cells required a  final fix\n",final);
01022     printf("[Step15b]: %i cells had an unphysical state\n",fail);
01023   }
01024 #endif /* MHD */
01025 #endif /* FIRST_ORDER_FLUX_CORRECTION */
01026 
01027 #ifdef STATIC_MESH_REFINEMENT
01028 /*=== STEP 16: With SMR, store fluxes at fine/coarse boundaries ==============*/
01029 
01030 /*--- Step 16a -----------------------------------------------------------------
01031  * store fluxes at boundaries of child grids.
01032  */
01033 
01034   for (ncg=0; ncg<pG->NCGrid; ncg++) {
01035 
01036 /* x1-boundaries of child Grids (interior to THIS Grid) */
01037 
01038     for (dim=0; dim<2; dim++){
01039       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01040 
01041         if (dim==0) i = pG->CGrid[ncg].ijks[0];
01042         if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
01043         jcs = pG->CGrid[ncg].ijks[1];
01044         jce = pG->CGrid[ncg].ijke[1];
01045 
01046         for (j=jcs, jj=0; j<=jce; j++, jj++){
01047           pG->CGrid[ncg].myFlx[dim][ks][jj].d  = x1Flux[j][i].d;
01048           pG->CGrid[ncg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01049           pG->CGrid[ncg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01050           pG->CGrid[ncg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01051 #ifndef BAROTROPIC
01052           pG->CGrid[ncg].myFlx[dim][ks][jj].E  = x1Flux[j][i].E;
01053 #endif /* BAROTROPIC */
01054 #ifdef MHD
01055           pG->CGrid[ncg].myFlx[dim][ks][jj].B1c = 0.0;
01056           pG->CGrid[ncg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01057           pG->CGrid[ncg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01058 #endif /* MHD */
01059 #if (NSCALARS > 0)
01060           for (n=0; n<NSCALARS; n++)
01061             pG->CGrid[ncg].myFlx[dim][ks][jj].s[n]  = x1Flux[j][i].s[n];
01062 #endif
01063         }
01064 #ifdef MHD
01065         for (j=jcs, jj=0; j<=jce+1; j++, jj++){
01066           pG->CGrid[ncg].myEMF3[dim][ks][jj] = emf3[j][i];
01067         }
01068 #endif /* MHD */
01069       }
01070     }
01071 
01072 /* x2-boundaries of child Grids (interior to THIS Grid) */
01073 
01074     for (dim=2; dim<4; dim++){
01075       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01076 
01077         ics = pG->CGrid[ncg].ijks[0];
01078         ice = pG->CGrid[ncg].ijke[0];
01079         if (dim==2) j = pG->CGrid[ncg].ijks[1];
01080         if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
01081 
01082         for (i=ics, ii=0; i<=ice; i++, ii++){
01083           pG->CGrid[ncg].myFlx[dim][ks][ii].d  = x2Flux[j][i].d;
01084           pG->CGrid[ncg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01085           pG->CGrid[ncg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01086           pG->CGrid[ncg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01087 #ifndef BAROTROPIC
01088           pG->CGrid[ncg].myFlx[dim][ks][ii].E  = x2Flux[j][i].E;
01089 #endif /* BAROTROPIC */
01090 #ifdef MHD
01091           pG->CGrid[ncg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01092           pG->CGrid[ncg].myFlx[dim][ks][ii].B2c = 0.0;
01093           pG->CGrid[ncg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01094 #endif /* MHD */
01095 #if (NSCALARS > 0)
01096           for (n=0; n<NSCALARS; n++)
01097             pG->CGrid[ncg].myFlx[dim][ks][ii].s[n]  = x2Flux[j][i].s[n];
01098 #endif
01099         }
01100 #ifdef MHD
01101         for (i=ics, ii=0; i<=ice+1; i++, ii++){
01102           pG->CGrid[ncg].myEMF3[dim][ks][ii] = emf3[j][i];
01103         }
01104 #endif /* MHD */
01105       }
01106     }
01107   }
01108 
01109 /*--- Step 16b -----------------------------------------------------------------
01110  * store fluxes at boundaries with parent grids.
01111  */
01112 
01113   for (npg=0; npg<pG->NPGrid; npg++) {
01114 
01115 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
01116 
01117     for (dim=0; dim<2; dim++){
01118       if (pG->PGrid[npg].myFlx[dim] != NULL) {
01119 
01120         if (dim==0) i = pG->PGrid[npg].ijks[0];
01121         if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
01122         jps = pG->PGrid[npg].ijks[1];
01123         jpe = pG->PGrid[npg].ijke[1];
01124 
01125         for (j=jps, jj=0; j<=jpe; j++, jj++){
01126           pG->PGrid[npg].myFlx[dim][ks][jj].d  = x1Flux[j][i].d;
01127           pG->PGrid[npg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01128           pG->PGrid[npg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01129           pG->PGrid[npg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01130 #ifndef BAROTROPIC
01131           pG->PGrid[npg].myFlx[dim][ks][jj].E  = x1Flux[j][i].E;
01132 #endif /* BAROTROPIC */
01133 #ifdef MHD
01134           pG->PGrid[npg].myFlx[dim][ks][jj].B1c = 0.0;
01135           pG->PGrid[npg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01136           pG->PGrid[npg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01137 #endif /* MHD */
01138 #if (NSCALARS > 0)
01139           for (n=0; n<NSCALARS; n++)
01140             pG->PGrid[npg].myFlx[dim][ks][jj].s[n]  = x1Flux[j][i].s[n];
01141 #endif
01142         }
01143 #ifdef MHD
01144         for (j=jps, jj=0; j<=jpe+1; j++, jj++){
01145           pG->PGrid[npg].myEMF3[dim][ks][jj] = emf3[j][i];
01146         }
01147 #endif /* MHD */
01148       }
01149     }
01150 
01151 /* x2-boundaries of parent Grids (at boundaries of THIS Grid)  */
01152 
01153     for (dim=2; dim<4; dim++){
01154       if (pG->PGrid[npg].myFlx[dim] != NULL) {
01155 
01156         ips = pG->PGrid[npg].ijks[0];
01157         ipe = pG->PGrid[npg].ijke[0];
01158         if (dim==2) j = pG->PGrid[npg].ijks[1];
01159         if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
01160 
01161         for (i=ips, ii=0; i<=ipe; i++, ii++){
01162           pG->PGrid[npg].myFlx[dim][ks][ii].d  = x2Flux[j][i].d;
01163           pG->PGrid[npg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01164           pG->PGrid[npg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01165           pG->PGrid[npg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01166 #ifndef BAROTROPIC
01167           pG->PGrid[npg].myFlx[dim][ks][ii].E  = x2Flux[j][i].E;
01168 #endif /* BAROTROPIC */
01169 #ifdef MHD
01170           pG->PGrid[npg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01171           pG->PGrid[npg].myFlx[dim][ks][ii].B2c = 0.0;
01172           pG->PGrid[npg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01173 #endif /* MHD */
01174 #if (NSCALARS > 0)
01175           for (n=0; n<NSCALARS; n++)
01176             pG->PGrid[npg].myFlx[dim][ks][ii].s[n]  = x2Flux[j][i].s[n];
01177 #endif
01178         }
01179 #ifdef MHD
01180         for (i=ips, ii=0; i<=ipe+1; i++, ii++){
01181           pG->PGrid[npg].myEMF3[dim][ks][ii] = emf3[j][i];
01182         }
01183 #endif /* MHD */
01184       }
01185     }
01186   }
01187 
01188 #endif /* STATIC_MESH_REFINEMENT */
01189 
01190   return;
01191 }
01192 
01193 /*----------------------------------------------------------------------------*/
01194 /*! \fn void integrate_init_2d(MeshS *pM)
01195  *  \brief Allocate temporary integration arrays */
01196 void integrate_init_2d(MeshS *pM)
01197 {
01198   int nmax,size1=0,size2=0,nl,nd;
01199 
01200 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2 */
01201   for (nl=0; nl<(pM->NLevels); nl++){
01202     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01203       if (pM->Domain[nl][nd].Grid != NULL) {
01204         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01205           size1 = pM->Domain[nl][nd].Grid->Nx[0];
01206         }
01207         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01208           size2 = pM->Domain[nl][nd].Grid->Nx[1];
01209         }
01210       }
01211     }
01212   }
01213 
01214   size1 = size1 + 2*nghost;
01215   size2 = size2 + 2*nghost;
01216   nmax = MAX(size1,size2);
01217 
01218 #ifdef MHD
01219   if ((emf3 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01220     goto on_error;
01221   if ((emf3_cc = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01222     goto on_error;
01223 #endif /* MHD */
01224 #ifdef H_CORRECTION
01225   if ((eta1 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01226     goto on_error;
01227   if ((eta2 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01228     goto on_error;
01229 #endif /* H_CORRECTION */
01230   if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01231   if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01232 #ifdef MHD
01233   if ((B1_x1Face = (Real**)calloc_2d_array(size2,size1, sizeof(Real))) == NULL)
01234     goto on_error;
01235   if ((B2_x2Face = (Real**)calloc_2d_array(size2,size1, sizeof(Real))) == NULL)
01236     goto on_error;
01237 #endif /* MHD */
01238 
01239   if ((U1d= (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01240   if ((Ul = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01241   if ((Ur = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01242   if ((W1d= (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01243   if ((Wl = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01244   if ((Wr = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01245 
01246   if ((Wl_x1Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01247     == NULL) goto on_error;
01248   if ((Wr_x1Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01249     == NULL) goto on_error;
01250   if ((Wl_x2Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01251     == NULL) goto on_error;
01252   if ((Wr_x2Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01253     == NULL) goto on_error;
01254   if ((x1Flux    = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS))) 
01255     == NULL) goto on_error;
01256   if ((x2Flux    = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS))) 
01257     == NULL) goto on_error;
01258 #ifdef FIRST_ORDER_FLUX_CORRECTION
01259   if ((x1FluxP = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS))) 
01260     == NULL) goto on_error;
01261   if ((x2FluxP = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS))) 
01262     == NULL) goto on_error;
01263 #ifdef MHD
01264   if ((emf3P = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01265     goto on_error;
01266 #endif
01267 #endif /* FIRST_ORDER_FLUX_CORRECTION */
01268 #ifdef USE_ENTROPY_FIX
01269   if ((S = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01270     goto on_error;
01271   if ((Shalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01272     goto on_error;
01273   if ((x1FluxS = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01274     goto on_error;
01275   if ((x2FluxS = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01276     goto on_error;
01277 #ifdef FIRST_ORDER_FLUX_CORRECTION
01278   if ((x1FluxSP = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01279     goto on_error;
01280   if ((x2FluxSP = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01281     goto on_error;
01282 #endif
01283 #endif
01284 
01285   if ((Uhalf = (ConsS**)calloc_2d_array(size2,size1,sizeof(ConsS)))==NULL)
01286     goto on_error;
01287   if ((Whalf = (PrimS**)calloc_2d_array(size2,size1,sizeof(PrimS)))==NULL)
01288     goto on_error;
01289   if ((W     = (PrimS**)calloc_2d_array(size2,size1,sizeof(PrimS)))==NULL)
01290     goto on_error;
01291 
01292   return;
01293 
01294   on_error:
01295     integrate_destruct();
01296     ath_error("[integrate_init]: malloc returned a NULL pointer\n");
01297 }
01298 
01299 /*----------------------------------------------------------------------------*/
01300 /*! \fn void integrate_destruct_2d(void)
01301  *  \brief Free temporary integration arrays */
01302 void integrate_destruct_2d(void)
01303 {
01304 #ifdef MHD
01305   if (emf3    != NULL) free_2d_array(emf3);
01306   if (emf3_cc != NULL) free_2d_array(emf3_cc);
01307 #endif /* MHD */
01308 #ifdef H_CORRECTION
01309   if (eta1 != NULL) free_2d_array(eta1);
01310   if (eta2 != NULL) free_2d_array(eta2);
01311 #endif /* H_CORRECTION */
01312   if (Bxc != NULL) free(Bxc);
01313   if (Bxi != NULL) free(Bxi);
01314 #ifdef MHD
01315   if (B1_x1Face != NULL) free_2d_array(B1_x1Face);
01316   if (B2_x2Face != NULL) free_2d_array(B2_x2Face);
01317 #endif /* MHD */
01318 
01319   if (U1d != NULL) free(U1d);
01320   if (Ul  != NULL) free(Ul);
01321   if (Ur  != NULL) free(Ur);
01322   if (W1d != NULL) free(W1d);
01323   if (Wl  != NULL) free(Wl);
01324   if (Wr  != NULL) free(Wr);
01325 
01326   if (Wl_x1Face != NULL) free_2d_array(Wl_x1Face);
01327   if (Wr_x1Face != NULL) free_2d_array(Wr_x1Face);
01328   if (Wl_x2Face != NULL) free_2d_array(Wl_x2Face);
01329   if (Wr_x2Face != NULL) free_2d_array(Wr_x2Face);
01330   if (x1Flux    != NULL) free_2d_array(x1Flux);
01331   if (x2Flux    != NULL) free_2d_array(x2Flux);
01332 #ifdef FIRST_ORDER_FLUX_CORRECTION
01333   if (x1FluxP   != NULL) free_2d_array(x1FluxP);
01334   if (x2FluxP   != NULL) free_2d_array(x2FluxP);
01335 #ifdef MHD
01336   if (emf3P     != NULL) free_2d_array(emf3P);
01337 #endif
01338 #endif
01339 #ifdef USE_ENTROPY_FIX
01340   if (S         != NULL) free_2d_array(S);
01341   if (Shalf     != NULL) free_2d_array(Shalf);
01342   if (x1FluxS   != NULL) free_2d_array(x1FluxS);
01343   if (x2FluxS   != NULL) free_2d_array(x2FluxS);
01344 #ifdef FIRST_ORDER_FLUX_CORRECTION
01345   if (x1FluxSP  != NULL) free_2d_array(x1FluxSP);
01346   if (x2FluxSP  != NULL) free_2d_array(x2FluxSP);
01347 #endif
01348 #endif
01349 
01350   if (Uhalf    != NULL) free_2d_array(Uhalf);
01351   if (Whalf    != NULL) free_2d_array(Whalf);
01352   if (W        != NULL) free_2d_array(W);
01353 
01354   return;
01355 }
01356 
01357 /*=========================== PRIVATE FUNCTIONS ==============================*/
01358 
01359 /*----------------------------------------------------------------------------*/
01360 /*! \fn static void integrate_emf3_corner(GridS *pG)
01361  *  \brief Integrates face centered B-fluxes to compute corner EMFs. 
01362  *
01363  *   Note:
01364  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
01365  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
01366  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
01367  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
01368  */ 
01369 
01370 #ifdef MHD
01371 static void integrate_emf3_corner(GridS *pG)
01372 {
01373   int i,il,iu,j,jl,ju;
01374   Real emf_l1, emf_r1, emf_l2, emf_r2;
01375 
01376   il = pG->is-(nghost-1);   iu = pG->ie+(nghost-1);
01377   jl = pG->js-(nghost-1);   ju = pG->je+(nghost-1);
01378 
01379 /* NOTE: The x1-Flux of B2 is -E3.  The x2-Flux of B1 is +E3. */
01380   for (j=jl; j<=ju+1; j++) {
01381     for (i=il; i<=iu+1; i++) {
01382       if (x1Flux[j-1][i].d > 0.0) {
01383         emf_l2 = -x1Flux[j-1][i].By
01384           + (x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1]);
01385       }
01386       else if (x1Flux[j-1][i].d < 0.0) {
01387         emf_l2 = -x1Flux[j-1][i].By
01388           + (x2Flux[j][i].Bz - emf3_cc[j-1][i]);
01389 
01390       } else {
01391         emf_l2 = -x1Flux[j-1][i].By
01392           + 0.5*(x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1] +
01393                  x2Flux[j][i  ].Bz - emf3_cc[j-1][i  ] );
01394       }
01395 
01396       if (x1Flux[j][i].d > 0.0) {
01397         emf_r2 = -x1Flux[j][i].By
01398           + (x2Flux[j][i-1].Bz - emf3_cc[j][i-1]);
01399       }
01400       else if (x1Flux[j][i].d < 0.0) {
01401         emf_r2 = -x1Flux[j][i].By
01402           + (x2Flux[j][i].Bz - emf3_cc[j][i]);
01403 
01404       } else {
01405         emf_r2 = -x1Flux[j][i].By
01406           + 0.5*(x2Flux[j][i-1].Bz - emf3_cc[j][i-1] +
01407                  x2Flux[j][i  ].Bz - emf3_cc[j][i  ] );
01408       }
01409 
01410       if (x2Flux[j][i-1].d > 0.0) {
01411         emf_l1 = x2Flux[j][i-1].Bz
01412           + (-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]);
01413       }
01414       else if (x2Flux[j][i-1].d < 0.0) {
01415         emf_l1 = x2Flux[j][i-1].Bz
01416           + (-x1Flux[j][i].By - emf3_cc[j][i-1]);
01417       } else {
01418         emf_l1 = x2Flux[j][i-1].Bz
01419           + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]
01420                  -x1Flux[j  ][i].By - emf3_cc[j  ][i-1] );
01421       }
01422 
01423       if (x2Flux[j][i].d > 0.0) {
01424         emf_r1 = x2Flux[j][i].Bz
01425           + (-x1Flux[j-1][i].By - emf3_cc[j-1][i]);
01426       }
01427       else if (x2Flux[j][i].d < 0.0) {
01428         emf_r1 = x2Flux[j][i].Bz
01429           + (-x1Flux[j][i].By - emf3_cc[j][i]);
01430       } else {
01431         emf_r1 = x2Flux[j][i].Bz
01432           + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i]
01433                  -x1Flux[j  ][i].By - emf3_cc[j  ][i] );
01434       }
01435 
01436       emf3[j][i] = 0.25*(emf_l1 + emf_r1 + emf_l2 + emf_r2);
01437     }
01438   }
01439 
01440   return;
01441 }
01442 #endif /* MHD */
01443 
01444 #ifdef FIRST_ORDER_FLUX_CORRECTION
01445 /*----------------------------------------------------------------------------*/
01446 /*! \fn static void FixCell(GridS *pG, Int3Vect ix)
01447  *  \brief Uses first order fluxes to fix negative d,P or superluminal v
01448  */ 
01449 static void FixCell(GridS *pG, Int3Vect ix)
01450 {
01451   int ks=pG->ks;
01452   Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2;
01453   Cons1DS x1FD_i, x1FD_ip1, x2FD_j, x2FD_jp1;
01454 #ifdef MHD
01455   int i,j;
01456   Real emf3D_ji, emf3D_jip1, emf3D_jp1i, emf3D_jp1ip1;
01457 #endif
01458 #ifdef USE_ENTROPY_FIX
01459   Real x1SFD_i, x1SFD_ip1, x2SFD_j, x2SFD_jp1;
01460 #endif
01461 
01462 
01463 /* Compute difference of predictor and corrector fluxes at cell faces */
01464 
01465   x1FD_i.d = x1Flux[ix.j][ix.i].d - x1FluxP[ix.j][ix.i].d;
01466   x2FD_j.d = x2Flux[ix.j][ix.i].d - x2FluxP[ix.j][ix.i].d;
01467 
01468   x1FD_ip1.d = x1Flux[ix.j][ix.i+1].d - x1FluxP[ix.j][ix.i+1].d;
01469   x2FD_jp1.d = x2Flux[ix.j+1][ix.i].d - x2FluxP[ix.j+1][ix.i].d;
01470 
01471   x1FD_i.Mx = x1Flux[ix.j][ix.i].Mx - x1FluxP[ix.j][ix.i].Mx;
01472   x2FD_j.Mx = x2Flux[ix.j][ix.i].Mx - x2FluxP[ix.j][ix.i].Mx;
01473 
01474   x1FD_ip1.Mx = x1Flux[ix.j][ix.i+1].Mx - x1FluxP[ix.j][ix.i+1].Mx;
01475   x2FD_jp1.Mx = x2Flux[ix.j+1][ix.i].Mx - x2FluxP[ix.j+1][ix.i].Mx;
01476 
01477   x1FD_i.My = x1Flux[ix.j][ix.i].My - x1FluxP[ix.j][ix.i].My;
01478   x2FD_j.My = x2Flux[ix.j][ix.i].My - x2FluxP[ix.j][ix.i].My;
01479 
01480   x1FD_ip1.My = x1Flux[ix.j][ix.i+1].My - x1FluxP[ix.j][ix.i+1].My;
01481   x2FD_jp1.My = x2Flux[ix.j+1][ix.i].My - x2FluxP[ix.j+1][ix.i].My;
01482 
01483   x1FD_i.Mz = x1Flux[ix.j][ix.i].Mz - x1FluxP[ix.j][ix.i].Mz;
01484   x2FD_j.Mz = x2Flux[ix.j][ix.i].Mz - x2FluxP[ix.j][ix.i].Mz;
01485 
01486   x1FD_ip1.Mz = x1Flux[ix.j][ix.i+1].Mz - x1FluxP[ix.j][ix.i+1].Mz;
01487   x2FD_jp1.Mz = x2Flux[ix.j+1][ix.i].Mz - x2FluxP[ix.j+1][ix.i].Mz;
01488 
01489 #ifndef BAROTROPIC
01490   x1FD_i.E = x1Flux[ix.j][ix.i].E - x1FluxP[ix.j][ix.i].E;
01491   x2FD_j.E = x2Flux[ix.j][ix.i].E - x2FluxP[ix.j][ix.i].E;
01492 
01493   x1FD_ip1.E = x1Flux[ix.j][ix.i+1].E - x1FluxP[ix.j][ix.i+1].E;
01494   x2FD_jp1.E = x2Flux[ix.j+1][ix.i].E - x2FluxP[ix.j+1][ix.i].E;
01495 #endif /* BAROTROPIC */
01496 
01497 #ifdef MHD
01498   x1FD_i.Bz = x1Flux[ix.j][ix.i].Bz - x1FluxP[ix.j][ix.i].Bz;
01499   x2FD_j.By = x2Flux[ix.j][ix.i].By - x2FluxP[ix.j][ix.i].By;
01500 
01501   x1FD_ip1.Bz = x1Flux[ix.j][ix.i+1].Bz - x1FluxP[ix.j][ix.i+1].Bz;
01502   x2FD_jp1.By = x2Flux[ix.j+1][ix.i].By - x2FluxP[ix.j+1][ix.i].By;
01503 
01504   emf3D_ji     = emf3[ix.j  ][ix.i  ] - emf3P[ix.j  ][ix.i  ];
01505   emf3D_jip1   = emf3[ix.j  ][ix.i+1] - emf3P[ix.j  ][ix.i+1];
01506   emf3D_jp1i   = emf3[ix.j+1][ix.i  ] - emf3P[ix.j+1][ix.i  ];
01507   emf3D_jp1ip1 = emf3[ix.j+1][ix.i+1] - emf3P[ix.j+1][ix.i+1];
01508 #endif /* MHD */
01509 
01510   /*#ifdef USE_ENTROPY_FIX
01511   x1SFD_i = x1FluxS[ix.j][ix.i] - x1FluxSP[ix.j][ix.i];
01512   x2SFD_j = x2FluxS[ix.j][ix.i] - x2FluxSP[ix.j][ix.i];
01513 
01514   x1SFD_ip1 = x1FluxS[ix.j][ix.i+1] - x1FluxSP[ix.j][ix.i+1];
01515   x2SFD_jp1 = x2FluxS[ix.j+1][ix.i] - x2FluxSP[ix.j+1][ix.i];
01516   #endif*/
01517 
01518 /* Use flux differences to correct bad cell-centered quantities */
01519 
01520   pG->U[ks][ix.j][ix.i].d  += dtodx1*(x1FD_ip1.d  - x1FD_i.d );
01521   pG->U[ks][ix.j][ix.i].M1 += dtodx1*(x1FD_ip1.Mx - x1FD_i.Mx);
01522   pG->U[ks][ix.j][ix.i].M2 += dtodx1*(x1FD_ip1.My - x1FD_i.My);
01523   pG->U[ks][ix.j][ix.i].M3 += dtodx1*(x1FD_ip1.Mz - x1FD_i.Mz);
01524 #ifndef BAROTROPIC
01525   pG->U[ks][ix.j][ix.i].E  += dtodx1*(x1FD_ip1.E  - x1FD_i.E );
01526 #endif /* BAROTROPIC */
01527 #ifdef MHD
01528   pG->U[ks][ix.j][ix.i].B3c += dtodx1*(x1FD_ip1.Bz - x1FD_i.Bz);
01529 #endif /* MHD */
01530   /*#ifdef USE_ENTROPY_FIX
01531   S[ix.j][ix.i]  += dtodx1*(x1SFD_ip1  - x1SFD_i );
01532   #endif*/
01533 
01534   pG->U[ks][ix.j][ix.i].d  += dtodx2*(x2FD_jp1.d  - x2FD_j.d );
01535   pG->U[ks][ix.j][ix.i].M1 += dtodx2*(x2FD_jp1.Mz - x2FD_j.Mz);
01536   pG->U[ks][ix.j][ix.i].M2 += dtodx2*(x2FD_jp1.Mx - x2FD_j.Mx);
01537   pG->U[ks][ix.j][ix.i].M3 += dtodx2*(x2FD_jp1.My - x2FD_j.My);
01538 #ifndef BAROTROPIC
01539   pG->U[ks][ix.j][ix.i].E  += dtodx2*(x2FD_jp1.E  - x2FD_j.E );
01540 #endif /* BAROTROPIC */
01541 #ifdef MHD
01542   pG->U[ks][ix.j][ix.i].B3c += dtodx2*(x2FD_jp1.By - x2FD_j.By);
01543 #endif /* MHD */
01544   /*#ifdef USE_ENTROPY_FIX
01545   S[ix.j][ix.i]  += dtodx2*(x2SFD_jp1  - x2SFD_j );
01546   #endif*/
01547 
01548 /* Use flux differences to correct cell-centered values at i-1 and i+1 */
01549 
01550   if (ix.i > pG->is) {
01551     pG->U[ks][ix.j][ix.i-1].d  += dtodx1*(x1FD_i.d );
01552     pG->U[ks][ix.j][ix.i-1].M1 += dtodx1*(x1FD_i.Mx);
01553     pG->U[ks][ix.j][ix.i-1].M2 += dtodx1*(x1FD_i.My);
01554     pG->U[ks][ix.j][ix.i-1].M3 += dtodx1*(x1FD_i.Mz);
01555 #ifndef BAROTROPIC
01556     pG->U[ks][ix.j][ix.i-1].E  += dtodx1*(x1FD_i.E );
01557 #endif /* BAROTROPIC */
01558 #ifdef MHD
01559     pG->U[ks][ix.j][ix.i-1].B3c += dtodx1*(x1FD_i.Bz);
01560 #endif /* MHD */
01561     /*#ifdef USE_ENTROPY_FIX
01562     S[ix.j][ix.i-1]  += dtodx1*(x1SFD_i );
01563     #endif*/
01564   }
01565 
01566   if (ix.i < pG->ie) {
01567     pG->U[ks][ix.j][ix.i+1].d  -= dtodx1*(x1FD_ip1.d );
01568     pG->U[ks][ix.j][ix.i+1].M1 -= dtodx1*(x1FD_ip1.Mx);
01569     pG->U[ks][ix.j][ix.i+1].M2 -= dtodx1*(x1FD_ip1.My);
01570     pG->U[ks][ix.j][ix.i+1].M3 -= dtodx1*(x1FD_ip1.Mz);
01571 #ifndef BAROTROPIC
01572     pG->U[ks][ix.j][ix.i+1].E  -= dtodx1*(x1FD_ip1.E );
01573 #endif /* BAROTROPIC */
01574 #ifdef MHD
01575     pG->U[ks][ix.j][ix.i+1].B3c -= dtodx1*(x1FD_ip1.Bz);
01576 #endif /* MHD */
01577     /*#ifdef USE_ENTROPY_FIX
01578     S[ix.j][ix.i+1]  -= dtodx1*(x1SFD_ip1 );
01579     #endif*/
01580   }
01581 
01582 /* Use flux differences to correct cell-centered values at j-1 and j+1 */
01583 
01584   if (ix.j > pG->js) {
01585     pG->U[ks][ix.j-1][ix.i].d  += dtodx2*(x2FD_j.d );
01586     pG->U[ks][ix.j-1][ix.i].M1 += dtodx2*(x2FD_j.Mz);
01587     pG->U[ks][ix.j-1][ix.i].M2 += dtodx2*(x2FD_j.Mx);
01588     pG->U[ks][ix.j-1][ix.i].M3 += dtodx2*(x2FD_j.My);
01589 #ifndef BAROTROPIC
01590     pG->U[ks][ix.j-1][ix.i].E  += dtodx2*(x2FD_j.E );
01591 #endif /* BAROTROPIC */
01592 #ifdef MHD
01593     pG->U[ks][ix.j-1][ix.i].B3c += dtodx2*(x2FD_j.By);
01594 #endif /* MHD */
01595     /*#ifdef USE_ENTROPY_FIX
01596    S[ix.j-1][ix.i]  += dtodx2*(x2SFD_j );
01597    #endif*/
01598   }
01599 
01600   if (ix.j < pG->je) {
01601     pG->U[ks][ix.j+1][ix.i].d  -= dtodx2*(x2FD_jp1.d );
01602     pG->U[ks][ix.j+1][ix.i].M1 -= dtodx2*(x2FD_jp1.Mz);
01603     pG->U[ks][ix.j+1][ix.i].M2 -= dtodx2*(x2FD_jp1.Mx);
01604     pG->U[ks][ix.j+1][ix.i].M3 -= dtodx2*(x2FD_jp1.My);
01605 #ifndef BAROTROPIC
01606     pG->U[ks][ix.j+1][ix.i].E  -= dtodx2*(x2FD_jp1.E );
01607 #endif /* BAROTROPIC */
01608 #ifdef MHD
01609     pG->U[ks][ix.j+1][ix.i].B3c -= dtodx2*(x2FD_jp1.By);
01610 #endif /* MHD */
01611     /*#ifdef USE_ENTROPY_FIX
01612     S[ix.j+1][ix.i]  -= dtodx2*(x2SFD_jp1 );
01613     #endif*/
01614   }
01615 
01616 /* Use emf3 difference to correct face-centered fields on edges of bad cell */
01617 #ifdef MHD
01618   pG->B1i[ks][ix.j][ix.i  ] += dtodx2*(emf3D_jp1i   - emf3D_ji);
01619   pG->B1i[ks][ix.j][ix.i+1] += dtodx2*(emf3D_jp1ip1 - emf3D_jip1);
01620   pG->B2i[ks][ix.j  ][ix.i] -= dtodx1*(emf3D_jip1   - emf3D_ji);
01621   pG->B2i[ks][ix.j+1][ix.i] -= dtodx1*(emf3D_jp1ip1 - emf3D_jp1i);
01622 
01623 /* Use emf3 difference to correct face-centered fields around bad cell */
01624   if (ix.j > pG->js) {
01625     pG->B1i[ks][ix.j-1][ix.i  ] -= dtodx2*(emf3D_ji);
01626     pG->B1i[ks][ix.j-1][ix.i+1] -= dtodx2*(emf3D_jip1);
01627   }
01628   if (ix.j < pG->je) {
01629     pG->B1i[ks][ix.j+1][ix.i  ] += dtodx2*(emf3D_jp1i);
01630     pG->B1i[ks][ix.j+1][ix.i+1] += dtodx2*(emf3D_jp1ip1);
01631   }
01632 
01633   if (ix.i > pG->is) {
01634     pG->B2i[ks][ix.j  ][ix.i-1] += dtodx1*(emf3D_ji);
01635     pG->B2i[ks][ix.j+1][ix.i-1] += dtodx1*(emf3D_jp1i);
01636   }
01637   if (ix.i < pG->ie) {
01638     pG->B2i[ks][ix.j  ][ix.i+1] -= dtodx1*(emf3D_jip1);
01639     pG->B2i[ks][ix.j+1][ix.i+1] -= dtodx1*(emf3D_jp1ip1);
01640   }
01641 
01642 /* compute new cell-centered fields */
01643   for (j=(ix.j-1); j<=(ix.j+1); j++) {
01644   for (i=(ix.i-1); i<=(ix.i+1); i++) {
01645     pG->U[ks][j][i].B1c = 0.5*(pG->B1i[ks][j][i] + pG->B1i[ks][j][i+1]);
01646     pG->U[ks][j][i].B2c = 0.5*(pG->B2i[ks][j][i] + pG->B2i[ks][j+1][i]);
01647   }}
01648 #endif /* MHD */
01649 
01650 #ifdef STATIC_MESH_REFINEMENT
01651 /* With SMR, replace higher-order fluxes with predict fluxes in case they are
01652  * used at fine/coarse grid boundaries */ 
01653 
01654   x1Flux[ix.j][ix.i] = x1FluxP[ix.j][ix.i];
01655   x2Flux[ix.j][ix.i] = x2FluxP[ix.j][ix.i];
01656 #endif /* STATIC_MESH_REFINEMENT */
01657 
01658 }
01659 #endif /* FIRST_ORDER_FLUX_CORRECTION */
01660 
01661 #endif /* VL_INTEGRATOR */
01662 
01663 #endif /* SPECIAL_RELATIVITY */
integrators/integrate_2d_vl_sr.c
