Athena: integrators/integrate_2d

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