Athena: integrators/integrate_2d

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_2d_ctu.c
00004  *  \brief Integrate MHD equations in 2D using the directionally unsplit CTU
00005  *   method of Colella (1990). 
00006  *
00007  * PURPOSE: Integrate MHD equations in 2D using the directionally unsplit CTU
00008  *   method of Colella (1990).  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, optically-thin cooling,
00012  *   and H-correction of Sanders et al.
00013  *
00014  * REFERENCES:
00015  * - P. Colella, "Multidimensional upwind methods for hyperbolic conservation
00016  *   laws", JCP, 87, 171 (1990)
00017  *
00018  * - T. Gardiner & J.M. Stone, "An unsplit Godunov method for ideal MHD via
00019  *   constrained transport", JCP, 205, 509 (2005)
00020  *
00021  * - R. Sanders, E. Morano, & M.-C. Druguet, "Multidimensinal dissipation for
00022  *   upwind schemes: stability and applications to gas dynamics", JCP, 145, 511
00023  *   (1998)
00024  *
00025  * CONTAINS PUBLIC FUNCTIONS: 
00026  * - integrate_2d_ctu()
00027  * - integrate_init_2d()
00028  * - integrate_destruct_2d() */
00029 /*============================================================================*/
00030 
00031 #include <math.h>
00032 #include <stdio.h>
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040 #ifdef PARTICLES
00041 #include "../particles/particle.h"
00042 #endif
00043 
00044 #if defined(CTU_INTEGRATOR)
00045 #ifdef SPECIAL_RELATIVITY
00046 #error : The CTU integrator cannot be used for special relativity.
00047 #endif /* SPECIAL_RELATIVITY */
00048 
00049 /* The L/R states of conserved variables and fluxes at each cell face */
00050 static Cons1DS **Ul_x1Face=NULL, **Ur_x1Face=NULL;
00051 static Cons1DS **Ul_x2Face=NULL, **Ur_x2Face=NULL;
00052 static Cons1DS **x1Flux=NULL, **x2Flux=NULL;
00053 
00054 /* The interface magnetic fields and emfs */
00055 #ifdef MHD
00056 static Real **B1_x1Face=NULL, **B2_x2Face=NULL;
00057 static Real **emf3=NULL, **emf3_cc=NULL;
00058 #endif /* MHD */
00059 
00060 /* 1D scratch vectors used by lr_states and flux functions */
00061 static Real *Bxc=NULL, *Bxi=NULL;
00062 static Prim1DS *W=NULL, *Wl=NULL, *Wr=NULL;
00063 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00064 
00065 /* density and Pressure at t^{n+1/2} needed by MHD, cooling, and gravity */
00066 static Real **dhalf = NULL,**phalf = NULL;
00067 
00068 /* variables needed for H-correction of Sanders et al (1998) */
00069 extern Real etah;
00070 #ifdef H_CORRECTION
00071 static Real **eta1=NULL, **eta2=NULL;
00072 #endif
00073 
00074 /* VARIABLES NEEDED FOR CYLINDRICAL COORDINATES */
00075 #ifdef CYLINDRICAL
00076 static Real **geom_src=NULL;
00077 #endif
00078 
00079 /*==============================================================================
00080  * PRIVATE FUNCTION PROTOTYPES: 
00081  *   integrate_emf3_corner() - the upwind CT method in Gardiner & Stone (2005) 
00082  *============================================================================*/
00083 
00084 #ifdef MHD
00085 static void integrate_emf3_corner(GridS *pG);
00086 #endif
00087 
00088 /*=========================== PUBLIC FUNCTIONS ===============================*/
00089 /*----------------------------------------------------------------------------*/
00090 /*! \fn void integrate_2d_ctu(DomainS *pD)
00091  *  \brief CTU integrator in 2D.
00092  *
00093  *   The numbering of steps follows the numbering in the 3D version.
00094  *   NOT ALL STEPS ARE NEEDED IN 2D.
00095  */
00096 
00097 void integrate_2d_ctu(DomainS *pD)
00098 {
00099   GridS *pG=(pD->Grid);
00100   Real dtodx1 = pG->dt/pG->dx1, dtodx2 = pG->dt/pG->dx2;
00101   Real hdtodx1 = 0.5*dtodx1, hdtodx2 = 0.5*dtodx2;
00102   Real hdt = 0.5*pG->dt, dx2 = pG->dx2;
00103   int i,il,iu,is=pG->is, ie=pG->ie;
00104   int j,jl,ju,js=pG->js, je=pG->je;
00105   int ks=pG->ks;
00106   Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,M1h,M2h,M3h,Bx=0.0;
00107 #ifndef BAROTROPIC
00108   Real coolfl,coolfr,coolf,Eh=0.0;
00109 #endif
00110 #ifdef MHD
00111   Real MHD_src,dbx,dby,B1,B2,B3,V3;
00112   Real B1ch, B2ch, B3ch;
00113 #endif
00114 #ifdef H_CORRECTION
00115   Real cfr,cfl,lambdar,lambdal;
00116 #endif
00117 #if (NSCALARS > 0)
00118   int n;
00119 #endif
00120 #ifdef SELF_GRAVITY
00121   Real gxl,gxr,gyl,gyr,flux_m1l,flux_m1r,flux_m2l,flux_m2r;
00122 #endif
00123 #ifdef FEEDBACK
00124   Real dt1 = 1.0/pG->dt;
00125 #endif
00126 #ifdef SHEARING_BOX
00127 /* in XY shearing sheet 2=phi; in XZ shearing sheet 2=Z and 3=phi */
00128   Real Vphi, Mphi;
00129   Real M1n, dM2n, dM3n;   /* M1, dM2/3=(My+d*q*Omega_0*x) at time n */
00130   Real M1e, dM2e, dM3e;   /* M1, dM2/3 evolved by dt/2  */
00131   Real flx1_dM2, frx1_dM2, flx2_dM2, frx2_dM2;
00132   Real flx1_dM3, frx1_dM3, flx2_dM3, frx2_dM3;
00133   Real fact, qom, om_dt = Omega_0*pG->dt;
00134 #endif /* SHEARING_BOX */
00135 #ifdef STATIC_MESH_REFINEMENT
00136   int ncg,npg,dim;
00137   int ii,ics,ice,jj,jcs,jce,ips,ipe,jps,jpe;
00138 #endif
00139 
00140   /* VARIABLES NEEDED FOR CYLINDRICAL COORDINATES */
00141 #ifdef CYLINDRICAL
00142 #ifndef ISOTHERMAL
00143   Real Pavgh;
00144 #endif
00145   Real g,gl,gr,rinv;
00146   Real geom_src_d,geom_src_Vx,geom_src_Vy,geom_src_P,geom_src_By,geom_src_Bz;
00147   const Real *r=pG->r, *ri=pG->ri;
00148 #endif /* CYLINDRICAL */
00149   Real lsf=1.0, rsf=1.0;
00150 
00151 /* With particles, one more ghost cell must be updated in predict step */
00152 #ifdef PARTICLES
00153   Real d1;
00154   il = is - 3;
00155   iu = ie + 3;
00156   jl = js - 3;
00157   ju = je + 3;
00158 #else
00159   il = is - 2;
00160   iu = ie + 2;
00161   jl = js - 2;
00162   ju = je + 2;
00163 #endif
00164 
00165 /* Set etah=0 so first calls to flux functions do not use H-correction */
00166   etah = 0.0;
00167 
00168 /* Compute predictor feedback from particle drag */
00169 #ifdef FEEDBACK
00170   feedback_predictor(pG);
00171 #endif
00172 
00173 /*=== STEP 1: Compute L/R x1-interface states and 1D x1-Fluxes ===============*/
00174 
00175 /*--- Step 1a ------------------------------------------------------------------
00176  * Load 1D vector of conserved variables;
00177  * U1d = (d, M1, M2, M3, E, B2c, B3c, s[n])
00178  */
00179 
00180   for (j=jl; j<=ju; j++) {
00181     for (i=is-nghost; i<=ie+nghost; i++) {
00182       U1d[i].d  = pG->U[ks][j][i].d;
00183       U1d[i].Mx = pG->U[ks][j][i].M1;
00184       U1d[i].My = pG->U[ks][j][i].M2;
00185       U1d[i].Mz = pG->U[ks][j][i].M3;
00186 #ifndef BAROTROPIC
00187       U1d[i].E  = pG->U[ks][j][i].E;
00188 #endif /* BAROTROPIC */
00189 #ifdef MHD
00190       U1d[i].By = pG->U[ks][j][i].B2c;
00191       U1d[i].Bz = pG->U[ks][j][i].B3c;
00192       Bxc[i] = pG->U[ks][j][i].B1c;
00193       Bxi[i] = pG->B1i[ks][j][i];
00194       B1_x1Face[j][i] = pG->B1i[ks][j][i];
00195 #endif /* MHD */
00196 #if (NSCALARS > 0)
00197       for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[ks][j][i].s[n];
00198 #endif
00199     }
00200 
00201 /*--- Step 1b ------------------------------------------------------------------
00202  * Compute L and R states at X1-interfaces, add "MHD source terms" for 0.5*dt
00203  */
00204 
00205     for (i=is-nghost; i<=ie+nghost; i++) {
00206       W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00207 
00208       /* CALCULATE THE CELL-CENTERED GEOMETRIC SOURCE VECTOR NOW USING U^{n}
00209       * THIS WILL BE USED AT THE END OF THE INTEGRATION STEP AS A SOURCE TERM
00210       * FOR THE CELL-CENTERED CONSERVED VARIABLES (STEPS 6D,8B) */
00211 #ifdef CYLINDRICAL 
00212       geom_src[j][i]  = W[i].d*SQR(W[i].Vy);
00213 #ifdef MHD
00214       geom_src[j][i] += 0.5*(SQR(Bxc[i]) - SQR(W[i].By) + SQR(W[i].Bz));
00215 #endif
00216 #ifdef ISOTHERMAL
00217       geom_src[j][i] += Iso_csound2*W[i].d;
00218 #else
00219       geom_src[j][i] += W[i].P;
00220 #endif
00221       geom_src[j][i] /= x1vc(pG,i);
00222 #endif /* CYLINDRICAL */
00223     }
00224 
00225     lr_states(pG,W,Bxc,pG->dt,pG->dx1,il+1,iu-1,Wl,Wr,1);
00226 
00227 #ifdef MHD
00228     for (i=il+1; i<=iu; i++) {
00229 #ifdef CYLINDRICAL
00230       rsf = ri[i]/r[i-1];  lsf = ri[i-1]/r[i-1];
00231 #endif
00232       MHD_src = (pG->U[ks][j][i-1].M2/pG->U[ks][j][i-1].d)*
00233                (rsf*pG->B1i[ks][j][i] - lsf*pG->B1i[ks][j][i-1])/pG->dx1;
00234       Wl[i].By += hdt*MHD_src;
00235 
00236 #ifdef CYLINDRICAL
00237       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00238 #endif
00239       MHD_src = (pG->U[ks][j][i].M2/pG->U[ks][j][i].d)*
00240                (rsf*pG->B1i[ks][j][i+1] - lsf*pG->B1i[ks][j][i])/pG->dx1;
00241       Wr[i].By += hdt*MHD_src;
00242     }
00243 #endif /* MHD */
00244 
00245 /*--- Step 1c ------------------------------------------------------------------
00246  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00247  */
00248 
00249     if (StaticGravPot != NULL){
00250       for (i=il+1; i<=iu; i++) {
00251         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00252 #ifdef CYLINDRICAL
00253         gl = (*x1GravAcc)(x1vc(pG,i-1),x2,x3);
00254         gr = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00255 
00256         /* APPLY GRAV. SOURCE TERMS TO VELOCITY USING ACCELERATION FOR (dt/2) */
00257         Wl[i].Vx -= hdt*gl;
00258         Wr[i].Vx -= hdt*gr;
00259 #else
00260         phicr = (*StaticGravPot)( x1             ,x2,x3);
00261         phicl = (*StaticGravPot)((x1-    pG->dx1),x2,x3);
00262         phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00263 
00264         Wl[i].Vx -= dtodx1*(phifc - phicl);
00265         Wr[i].Vx -= dtodx1*(phicr - phifc);
00266 #endif
00267       }
00268     }
00269 
00270 /*--- Step 1c (cont) -----------------------------------------------------------
00271  * Add source terms for self-gravity for 0.5*dt to L/R states
00272  */
00273 
00274 #ifdef SELF_GRAVITY
00275     for (i=il+1; i<=iu; i++) {
00276       Wl[i].Vx -= hdtodx1*(pG->Phi[ks][j][i] - pG->Phi[ks][j][i-1]);
00277       Wr[i].Vx -= hdtodx1*(pG->Phi[ks][j][i] - pG->Phi[ks][j][i-1]);
00278     }
00279 #endif
00280 
00281 /*--- Step 1c (cont) -----------------------------------------------------------
00282  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00283  */
00284 
00285 #ifndef BAROTROPIC
00286     if (CoolingFunc != NULL){
00287       for (i=il+1; i<=iu; i++) {
00288         coolfl = (*CoolingFunc)(Wl[i].d,Wl[i].P,(0.5*pG->dt));
00289         coolfr = (*CoolingFunc)(Wr[i].d,Wr[i].P,(0.5*pG->dt));
00290 
00291         Wl[i].P -= 0.5*pG->dt*Gamma_1*coolfl;
00292         Wr[i].P -= 0.5*pG->dt*Gamma_1*coolfr;
00293       }
00294     }
00295 #endif /* BAROTROPIC */
00296 
00297 /*--- Step 1c (cont) -----------------------------------------------------------
00298  * Add source terms for shearing box (Coriolis forces) for 0.5*dt to L/R states
00299  * starting with tidal gravity terms added through the ShearingBoxPot
00300  *    Vx source term = (dt/2)*( 2 Omega_0 Vy)
00301  *    Vy source term = (dt/2)*(-2 Omega_0 Vx)
00302  *    Vy source term = (dt/2)*((q-2) Omega_0 Vx) (with FARGO)
00303  * (x1,x2,x3) in code = (X,Z,Y) in 2D shearing sheet
00304  */
00305 
00306 #ifdef SHEARING_BOX
00307     if (ShearingBoxPot != NULL){
00308       for (i=il+1; i<=iu; i++) {
00309         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00310         phicr = (*ShearingBoxPot)( x1             ,x2,x3);
00311         phicl = (*ShearingBoxPot)((x1-    pG->dx1),x2,x3);
00312         phifc = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00313 
00314         Wl[i].Vx -= dtodx1*(phifc - phicl);
00315         Wr[i].Vx -= dtodx1*(phicr - phifc);
00316       }
00317     }
00318 
00319     if (ShBoxCoord == xz){
00320       for (i=il+1; i<=iu; i++) {
00321         Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vz;
00322         Wr[i].Vx += pG->dt*Omega_0*W[i].Vz;
00323 #ifdef FARGO
00324         Wl[i].Vz += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00325         Wr[i].Vz += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00326 #else
00327         Wl[i].Vz -= pG->dt*Omega_0*W[i-1].Vx;
00328         Wr[i].Vz -= pG->dt*Omega_0*W[i].Vx;
00329 #endif
00330       }
00331     }
00332 
00333     if (ShBoxCoord == xy) {
00334       for (i=il+1; i<=iu; i++) {
00335         Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vy;
00336         Wr[i].Vx += pG->dt*Omega_0*W[i].Vy;
00337 #ifdef FARGO
00338         Wl[i].Vy += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00339         Wr[i].Vy += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00340 #else
00341         Wl[i].Vy -= pG->dt*Omega_0*W[i-1].Vx;
00342         Wr[i].Vy -= pG->dt*Omega_0*W[i].Vx;
00343 #endif
00344       }
00345     }
00346 #endif /* SHEARING_BOX */
00347 
00348 /*--- Step 1c (cont) -----------------------------------------------------------
00349  * Add source terms for particle feedback for 0.5*dt to L/R states
00350  */
00351 
00352 #ifdef FEEDBACK
00353     for (i=il+1; i<=iu; i++) {
00354       d1 = 1.0/W[i-1].d;
00355       Wl[i].Vx -= pG->Coup[ks][j][i-1].fb1*d1;
00356       Wl[i].Vy -= pG->Coup[ks][j][i-1].fb2*d1;
00357       Wl[i].Vz -= pG->Coup[ks][j][i-1].fb3*d1;
00358 
00359       d1 = 1.0/W[i].d;
00360       Wr[i].Vx -= pG->Coup[ks][j][i].fb1*d1;
00361       Wr[i].Vy -= pG->Coup[ks][j][i].fb2*d1;
00362       Wr[i].Vz -= pG->Coup[ks][j][i].fb3*d1;
00363 
00364 #ifndef BAROTROPIC
00365       Wl[i].P += pG->Coup[ks][j][i-1].Eloss*Gamma_1;
00366       Wr[i].P += pG->Coup[ks][j][i].Eloss*Gamma_1;
00367 #endif
00368     }
00369 #endif /* FEEDBACK */
00370 
00371 /*--- Step 1c (cont) -----------------------------------------------------------
00372  * ADD THE GEOMETRIC SOURCE-TERMS NOW USING CELL-CENTERED PRIMITIVE
00373  * VARIABLES AT TIME t^n
00374  */
00375 #ifdef CYLINDRICAL
00376       for (i=il+1; i<=iu; i++) {
00377 
00378         // LEFT STATE GEOMETRIC SOURCE TERM (USES W[i-1])
00379         rinv = 1.0/x1vc(pG,i-1);
00380         geom_src_d  = -W[i-1].d*W[i-1].Vx*rinv;
00381         geom_src_Vx =  SQR(W[i-1].Vy);
00382         geom_src_Vy = -W[i-1].Vx*W[i-1].Vy;
00383 #ifdef MHD
00384         geom_src_Vx -= SQR(W[i-1].By)/W[i-1].d;
00385         geom_src_Vy += Bxc[i-1]*W[i-1].By/W[i-1].d;
00386         geom_src_By =  -W[i-1].Vy*Bxc[i-1]*rinv;
00387         geom_src_Bz =  -W[i-1].Vx*W[i-1].Bz*rinv;
00388 #endif /* MHD */
00389         geom_src_Vx *= rinv;
00390         geom_src_Vy *= rinv;
00391 #ifndef ISOTHERMAL
00392         geom_src_P  = -Gamma*W[i-1].P*W[i-1].Vx*rinv;
00393 #endif /* ISOTHERMAL */
00394 
00395         // ADD SOURCE TERM TO LEFT STATE
00396         Wl[i].d  += hdt*geom_src_d;
00397         Wl[i].Vx += hdt*geom_src_Vx;
00398         Wl[i].Vy += hdt*geom_src_Vy;
00399 #ifdef MHD
00400         Wl[i].By += hdt*geom_src_By;
00401         Wl[i].Bz += hdt*geom_src_Bz;
00402 #endif /* MHD */
00403 #ifndef ISOTHERMAL
00404         Wl[i].P  += hdt*geom_src_P;
00405 #endif /* ISOTHERMAL */
00406 
00407         // RIGHT STATE GEOMETRIC SOURCE TERM (USES W[i])
00408         rinv = 1.0/x1vc(pG,i);
00409         geom_src_d  = -W[i].d*W[i].Vx*rinv;
00410         geom_src_Vx =  SQR(W[i].Vy);
00411         geom_src_Vy = -W[i].Vx*W[i].Vy;
00412 #ifdef MHD
00413         geom_src_Vx -= SQR(W[i].By)/W[i].d;
00414         geom_src_Vy += Bxc[i]*W[i].By/W[i].d;
00415         geom_src_By =  -W[i].Vy*Bxc[i]*rinv;
00416         geom_src_Bz =  -W[i].Vx*W[i].Bz*rinv;
00417 #endif /* MHD */
00418         geom_src_Vx *= rinv;
00419         geom_src_Vy *= rinv;
00420 #ifndef ISOTHERMAL
00421         geom_src_P  = -Gamma*W[i].P*W[i].Vx*rinv;
00422 #endif /* ISOTHERMAL */
00423 
00424         // ADD SOURCE TERM TO RIGHT STATE
00425         Wr[i].d  += hdt*geom_src_d;
00426         Wr[i].Vx += hdt*geom_src_Vx;
00427         Wr[i].Vy += hdt*geom_src_Vy;
00428 #ifdef MHD
00429         Wr[i].By += hdt*geom_src_By;
00430         Wr[i].Bz += hdt*geom_src_Bz;
00431 #endif /* MHD */
00432 #ifndef ISOTHERMAL
00433         Wr[i].P  += hdt*geom_src_P;
00434 #endif /* ISOTHERMAL */
00435       }
00436 #endif /* CYLINDRICAL */
00437 
00438 /*--- Step 1d ------------------------------------------------------------------
00439  * Compute 1D fluxes in x1-direction, storing into 2D array
00440  */
00441 
00442     for (i=il+1; i<=iu; i++) {
00443       Ul_x1Face[j][i] = Prim1D_to_Cons1D(&Wl[i],&Bxi[i]);
00444       Ur_x1Face[j][i] = Prim1D_to_Cons1D(&Wr[i],&Bxi[i]);
00445 
00446 #ifdef MHD
00447       Bx = B1_x1Face[j][i];
00448 #endif
00449       fluxes(Ul_x1Face[j][i],Ur_x1Face[j][i],Wl[i],Wr[i],Bx,&x1Flux[j][i]);
00450     }
00451   }
00452 
00453 /*=== STEP 2: Compute L/R x2-interface states and 1D x2-Fluxes ===============*/
00454 
00455 /*--- Step 2a ------------------------------------------------------------------
00456  * Load 1D vector of conserved variables;
00457  * U1d = (d, M2, M3, M1, E, B3c, B1c, s[n])
00458  */
00459 
00460   for (i=il; i<=iu; i++) {
00461 #ifdef CYLINDRICAL
00462     dx2 = r[i]*pG->dx2;
00463     dtodx2 = pG->dt/dx2;
00464     hdtodx2 = 0.5*dtodx2;
00465 #endif
00466     for (j=js-nghost; j<=je+nghost; j++) {
00467       U1d[j].d  = pG->U[ks][j][i].d;
00468       U1d[j].Mx = pG->U[ks][j][i].M2;
00469       U1d[j].My = pG->U[ks][j][i].M3;
00470       U1d[j].Mz = pG->U[ks][j][i].M1;
00471 #ifndef BAROTROPIC
00472       U1d[j].E  = pG->U[ks][j][i].E;
00473 #endif /* BAROTROPIC */
00474 #ifdef MHD
00475       U1d[j].By = pG->U[ks][j][i].B3c;
00476       U1d[j].Bz = pG->U[ks][j][i].B1c;
00477       Bxc[j] = pG->U[ks][j][i].B2c;
00478       Bxi[j] = pG->B2i[ks][j][i];
00479       B2_x2Face[j][i] = pG->B2i[ks][j][i];
00480 #endif /* MHD */
00481 #if (NSCALARS > 0)
00482       for (n=0; n<NSCALARS; n++) U1d[j].s[n] = pG->U[ks][j][i].s[n];
00483 #endif
00484     }
00485 
00486 /*--- Step 2b ------------------------------------------------------------------
00487  * Compute L and R states at X2-interfaces, add "MHD source terms" for 0.5*dt
00488  */
00489 
00490     for (j=js-nghost; j<=je+nghost; j++) {
00491       W[j] = Cons1D_to_Prim1D(&U1d[j],&Bxc[j]);
00492     }
00493 
00494     lr_states(pG,W,Bxc,pG->dt,dx2,jl+1,ju-1,Wl,Wr,2);
00495 
00496 #ifdef MHD
00497     for (j=jl+1; j<=ju; j++) {
00498       MHD_src = (pG->U[ks][j-1][i].M1/pG->U[ks][j-1][i].d)*
00499         (pG->B2i[ks][j][i] - pG->B2i[ks][j-1][i])/dx2;
00500       Wl[j].Bz += hdt*MHD_src;
00501 
00502       MHD_src = (pG->U[ks][j][i].M1/pG->U[ks][j][i].d)*
00503         (pG->B2i[ks][j+1][i] - pG->B2i[ks][j][i])/dx2;
00504       Wr[j].Bz += hdt*MHD_src;
00505     }
00506 #endif
00507 
00508 /*--- Step 2c ------------------------------------------------------------------
00509  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00510  */
00511   
00512     if (StaticGravPot != NULL){
00513       for (j=jl+1; j<=ju; j++) {
00514         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00515         phicr = (*StaticGravPot)(x1, x2             ,x3);
00516         phicl = (*StaticGravPot)(x1,(x2-    pG->dx2),x3);
00517         phifc = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00518 
00519         Wl[j].Vx -= dtodx2*(phifc - phicl);
00520         Wr[j].Vx -= dtodx2*(phicr - phifc);
00521       }
00522     }
00523 
00524 /*--- Step 2c (cont) -----------------------------------------------------------
00525  * Add source terms for self-gravity for 0.5*dt to L/R states
00526  */
00527 
00528 #ifdef SELF_GRAVITY
00529     for (j=jl+1; j<=ju; j++) {
00530       Wl[j].Vx -= hdtodx2*(pG->Phi[ks][j][i] - pG->Phi[ks][j-1][i]);
00531       Wr[j].Vx -= hdtodx2*(pG->Phi[ks][j][i] - pG->Phi[ks][j-1][i]);
00532     }
00533 #endif
00534 
00535 /*--- Step 2c (cont) -----------------------------------------------------------
00536  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00537  */
00538 
00539 #ifndef BAROTROPIC
00540     if (CoolingFunc != NULL){
00541       for (j=jl+1; j<=ju; j++) {
00542         coolfl = (*CoolingFunc)(Wl[j].d,Wl[j].P,(0.5*pG->dt));
00543         coolfr = (*CoolingFunc)(Wr[j].d,Wr[j].P,(0.5*pG->dt));
00544 
00545         Wl[j].P -= 0.5*pG->dt*Gamma_1*coolfl;
00546         Wr[j].P -= 0.5*pG->dt*Gamma_1*coolfr;
00547       }
00548     }
00549 #endif /* BAROTROPIC */
00550 
00551 /*--- Step 2c (cont) -----------------------------------------------------------
00552  * Add source terms for particle feedback for 0.5*dt to L/R states
00553  */
00554 
00555 #ifdef FEEDBACK
00556    for (j=jl+1; j<=ju; j++) {
00557       d1 = 1.0/W[j-1].d;
00558       Wl[j].Vx -= pG->Coup[ks][j-1][i].fb2*d1;
00559       Wl[j].Vy -= pG->Coup[ks][j-1][i].fb3*d1;
00560       Wl[j].Vz -= pG->Coup[ks][j-1][i].fb1*d1;
00561 
00562       d1 = 1.0/W[j].d;
00563       Wr[j].Vx -= pG->Coup[ks][j][i].fb2*d1;
00564       Wr[j].Vy -= pG->Coup[ks][j][i].fb3*d1;
00565       Wr[j].Vz -= pG->Coup[ks][j][i].fb1*d1;
00566 
00567 #ifndef BAROTROPIC
00568       Wl[i].P += pG->Coup[ks][j-1][i].Eloss*Gamma_1;
00569       Wr[i].P += pG->Coup[ks][j][i].Eloss*Gamma_1;
00570 #endif
00571     }
00572 #endif /* FEEDBACK */
00573 
00574 /*--- Step 2d ------------------------------------------------------------------
00575  * Compute 1D fluxes in x2-direction, storing into 2D array
00576  */
00577 
00578     for (j=jl+1; j<=ju; j++) {
00579       Ul_x2Face[j][i] = Prim1D_to_Cons1D(&Wl[j],&Bxi[j]);
00580       Ur_x2Face[j][i] = Prim1D_to_Cons1D(&Wr[j],&Bxi[j]);
00581 #ifdef MHD
00582       Bx = B2_x2Face[j][i];
00583 #endif
00584       fluxes(Ul_x2Face[j][i],Ur_x2Face[j][i],Wl[j],Wr[j],Bx,&x2Flux[j][i]);
00585     }
00586   }
00587 
00588 /*=== STEP 3: Not needed in 2D ===*/
00589 
00590 /*=== STEP 4:  Update face-centered B for 0.5*dt =============================*/
00591 
00592 /*--- Step 4a ------------------------------------------------------------------
00593  * Calculate the cell centered value of emf_3 at t^{n} and integrate
00594  * to corner.
00595  */
00596 
00597 #ifdef MHD
00598   for (j=jl; j<=ju; j++) {
00599     for (i=il; i<=iu; i++) {
00600       emf3_cc[j][i] =
00601         (pG->U[ks][j][i].B1c*pG->U[ks][j][i].M2 -
00602          pG->U[ks][j][i].B2c*pG->U[ks][j][i].M1 )/pG->U[ks][j][i].d;
00603     }
00604   }
00605   integrate_emf3_corner(pG);
00606 
00607 /*--- Step 4b ------------------------------------------------------------------
00608  * Update the interface magnetic fields using CT for a half time step.
00609  */
00610 
00611   for (j=jl+1; j<=ju-1; j++) {
00612     for (i=il+1; i<=iu-1; i++) {
00613 #ifdef CYLINDRICAL
00614       hdtodx2 = hdt/(ri[i]*pG->dx2);
00615 #endif
00616       B1_x1Face[j][i] -= hdtodx2*(emf3[j+1][i  ] - emf3[j][i]);
00617       B2_x2Face[j][i] += hdtodx1*(emf3[j  ][i+1] - emf3[j][i]);
00618     }
00619 #ifdef CYLINDRICAL
00620     hdtodx2 = hdt/(ri[iu]*pG->dx2);
00621 #endif
00622     B1_x1Face[j][iu] -= hdtodx2*(emf3[j+1][iu] - emf3[j][iu]);
00623   }
00624   for (i=il+1; i<=iu-1; i++) {
00625     B2_x2Face[ju][i] += hdtodx1*(emf3[ju][i+1] - emf3[ju][i]);
00626   }
00627 #endif /* MHD */
00628 
00629 /*=== STEP 5: Correct x1-interface states with transverse flux gradients =====*/
00630 
00631 /*--- Step 5a ------------------------------------------------------------------
00632  * Correct x1-interface states using x2-fluxes computed in Step 2d.
00633  * Since the fluxes come from an x2-sweep, (x,y,z) on RHS -> (z,x,y) on LHS
00634  */
00635 
00636   for (j=jl+1; j<=ju-1; j++) {
00637     for (i=il+1; i<=iu; i++) {
00638 #ifdef CYLINDRICAL
00639       hdtodx2 = hdt/(r[i-1]*pG->dx2);
00640 #endif
00641       Ul_x1Face[j][i].d  -= hdtodx2*(x2Flux[j+1][i-1].d  - x2Flux[j][i-1].d );
00642       Ul_x1Face[j][i].Mx -= hdtodx2*(x2Flux[j+1][i-1].Mz - x2Flux[j][i-1].Mz);
00643       Ul_x1Face[j][i].My -= hdtodx2*(x2Flux[j+1][i-1].Mx - x2Flux[j][i-1].Mx);
00644       Ul_x1Face[j][i].Mz -= hdtodx2*(x2Flux[j+1][i-1].My - x2Flux[j][i-1].My);
00645 #ifndef BAROTROPIC
00646       Ul_x1Face[j][i].E  -= hdtodx2*(x2Flux[j+1][i-1].E  - x2Flux[j][i-1].E );
00647 #endif /* BAROTROPIC */
00648 #ifdef MHD
00649       Ul_x1Face[j][i].Bz -= hdtodx2*(x2Flux[j+1][i-1].By - x2Flux[j][i-1].By);
00650 #endif
00651 #if (NSCALARS > 0)
00652       for (n=0; n<NSCALARS; n++) {
00653       Ul_x1Face[j][i].s[n]-=hdtodx2*(x2Flux[j+1][i-1].s[n]-x2Flux[j][i-1].s[n]);
00654       }
00655 #endif
00656 
00657 #ifdef CYLINDRICAL
00658       hdtodx2 = hdt/(r[i]*pG->dx2);
00659 #endif
00660       Ur_x1Face[j][i].d  -= hdtodx2*(x2Flux[j+1][i  ].d  - x2Flux[j][i  ].d );
00661       Ur_x1Face[j][i].Mx -= hdtodx2*(x2Flux[j+1][i  ].Mz - x2Flux[j][i  ].Mz);
00662       Ur_x1Face[j][i].My -= hdtodx2*(x2Flux[j+1][i  ].Mx - x2Flux[j][i  ].Mx);
00663       Ur_x1Face[j][i].Mz -= hdtodx2*(x2Flux[j+1][i  ].My - x2Flux[j][i  ].My);
00664 #ifndef BAROTROPIC
00665       Ur_x1Face[j][i].E  -= hdtodx2*(x2Flux[j+1][i  ].E  - x2Flux[j][i  ].E );
00666 #endif /* BAROTROPIC */
00667 #ifdef MHD
00668       Ur_x1Face[j][i].Bz -= hdtodx2*(x2Flux[j+1][i  ].By - x2Flux[j][i  ].By);
00669 #endif
00670 #if (NSCALARS > 0)
00671       for (n=0; n<NSCALARS; n++) {
00672       Ur_x1Face[j][i].s[n]-=hdtodx2*(x2Flux[j+1][i  ].s[n]-x2Flux[j][i  ].s[n]);
00673       }
00674 #endif
00675     }
00676   }
00677 
00678 /*--- Step 5b: Not needed in 2D ---*/
00679 /*--- Step 5c ------------------------------------------------------------------
00680  * Add the "MHD source terms" from the x2-flux gradient to the conservative
00681  * variables on the x1Face..
00682  */
00683 
00684 #ifdef MHD
00685   for (j=jl+1; j<=ju-1; j++) {
00686     for (i=il+1; i<=iu; i++) {
00687 #ifdef CYLINDRICAL
00688       dbx = (ri[i]*pG->B1i[ks][j][i] - ri[i-1]*pG->B1i[ks][j][i-1])/r[i-1];
00689 #else
00690       dbx = pG->B1i[ks][j][i] - pG->B1i[ks][j][i-1];
00691 #endif
00692       B1 = pG->U[ks][j][i-1].B1c;
00693       B2 = pG->U[ks][j][i-1].B2c;
00694       B3 = pG->U[ks][j][i-1].B3c;
00695       V3 = pG->U[ks][j][i-1].M3/pG->U[ks][j][i-1].d;
00696 
00697       Ul_x1Face[j][i].Mx += hdtodx1*B1*dbx;
00698       Ul_x1Face[j][i].My += hdtodx1*B2*dbx;
00699       Ul_x1Face[j][i].Mz += hdtodx1*B3*dbx;
00700       Ul_x1Face[j][i].Bz += hdtodx1*V3*dbx;
00701 #ifndef BAROTROPIC
00702       Ul_x1Face[j][i].E  += hdtodx1*B3*V3*dbx;
00703 #endif /* BAROTROPIC */
00704 
00705 #ifdef CYLINDRICAL
00706       dbx = (ri[i+1]*pG->B1i[ks][j][i+1] - ri[i]*pG->B1i[ks][j][i])/r[i];
00707 #else
00708       dbx = pG->B1i[ks][j][i+1] - pG->B1i[ks][j][i];
00709 #endif
00710       B1 = pG->U[ks][j][i].B1c;
00711       B2 = pG->U[ks][j][i].B2c;
00712       B3 = pG->U[ks][j][i].B3c;
00713       V3 = pG->U[ks][j][i].M3/pG->U[ks][j][i].d;
00714 
00715       Ur_x1Face[j][i].Mx += hdtodx1*B1*dbx;
00716       Ur_x1Face[j][i].My += hdtodx1*B2*dbx;
00717       Ur_x1Face[j][i].Mz += hdtodx1*B3*dbx;
00718       Ur_x1Face[j][i].Bz += hdtodx1*V3*dbx;
00719 #ifndef BAROTROPIC
00720       Ur_x1Face[j][i].E  += hdtodx1*B3*V3*dbx;
00721 #endif /* BAROTROPIC */
00722     }
00723   }
00724 #endif /* MHD */
00725 
00726 /*--- Step 5d ------------------------------------------------------------------
00727  * Add source terms for a static gravitational potential arising from x2-Flux
00728  * gradients.  To improve conservation of total energy, average
00729  * the energy source term computed at cell faces.
00730  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00731  */
00732 
00733   if (StaticGravPot != NULL){
00734     for (j=jl+1; j<=ju-1; j++) {
00735       for (i=il+1; i<=iu; i++) {
00736         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00737         phic = (*StaticGravPot)(x1, x2             ,x3);
00738         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00739         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00740 
00741 #ifdef CYLINDRICAL
00742         hdtodx2 = hdt/(r[i]*pG->dx2);
00743 #endif
00744         Ur_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00745 #ifndef BAROTROPIC
00746         Ur_x1Face[j][i].E -= hdtodx2*(x2Flux[j  ][i  ].d*(phic - phil) +
00747                                       x2Flux[j+1][i  ].d*(phir - phic));
00748 #endif
00749 
00750         phic = (*StaticGravPot)((x1-pG->dx1), x2             ,x3);
00751         phir = (*StaticGravPot)((x1-pG->dx1),(x2+0.5*pG->dx2),x3);
00752         phil = (*StaticGravPot)((x1-pG->dx1),(x2-0.5*pG->dx2),x3);
00753 
00754 #ifdef CYLINDRICAL
00755         hdtodx2 = hdt/(r[i-1]*pG->dx2);
00756 #endif
00757         Ul_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i-1].d;
00758 #ifndef BAROTROPIC
00759         Ul_x1Face[j][i].E -= hdtodx2*(x2Flux[j  ][i-1].d*(phic - phil) +
00760                                       x2Flux[j+1][i-1].d*(phir - phic));
00761 #endif
00762       }
00763     }
00764   }
00765 
00766 /*--- Step 5d (cont) -----------------------------------------------------------
00767  * Add source terms for self gravity arising from x2-Flux gradient
00768  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00769  */
00770 
00771 #ifdef SELF_GRAVITY
00772   for (j=jl+1; j<=ju-1; j++) {
00773     for (i=il+1; i<=iu; i++) {
00774       phic = pG->Phi[ks][j][i];
00775       phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j+1][i]);
00776       phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j-1][i]);
00777 
00778       Ur_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00779 #ifndef BAROTROPIC
00780       Ur_x1Face[j][i].E -= hdtodx2*(x2Flux[j  ][i  ].d*(phic - phil) +
00781                                     x2Flux[j+1][i  ].d*(phir - phic));
00782 #endif
00783 
00784       phic = pG->Phi[ks][j][i-1];
00785       phir = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j+1][i-1]);
00786       phil = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j-1][i-1]);
00787 
00788       Ul_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i-1].d;
00789 #ifndef BAROTROPIC
00790       Ul_x1Face[j][i].E -= hdtodx2*(x2Flux[j  ][i-1].d*(phic - phil) +
00791                                     x2Flux[j+1][i-1].d*(phir - phic));
00792 #endif
00793     }
00794   }
00795 #endif /* SELF_GRAVITY */
00796 
00797 /*=== STEP 6: Correct x2-interface states with transverse flux gradients =====*/
00798 
00799 /*--- Step 6a ------------------------------------------------------------------
00800  * Correct x2-interface states using x1-fluxes computed in Step 1d.
00801  * Since the fluxes come from an x1-sweep, (x,y,z) on RHS -> (y,z,x) on LHS
00802  */
00803 
00804   for (j=jl+1; j<=ju; j++) {
00805     for (i=il+1; i<=iu-1; i++) {
00806 #ifdef CYLINDRICAL
00807       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00808 #endif
00809       Ul_x2Face[j][i].d  -= hdtodx1*(rsf*x1Flux[j-1][i+1].d  - lsf*x1Flux[j-1][i].d );
00810       Ul_x2Face[j][i].Mx -= hdtodx1*(SQR(rsf)*x1Flux[j-1][i+1].My - SQR(lsf)*x1Flux[j-1][i].My);
00811       Ul_x2Face[j][i].My -= hdtodx1*(rsf*x1Flux[j-1][i+1].Mz - lsf*x1Flux[j-1][i].Mz);
00812       Ul_x2Face[j][i].Mz -= hdtodx1*(rsf*x1Flux[j-1][i+1].Mx - lsf*x1Flux[j-1][i].Mx);
00813 #ifndef BAROTROPIC
00814       Ul_x2Face[j][i].E  -= hdtodx1*(rsf*x1Flux[j-1][i+1].E  - lsf*x1Flux[j-1][i].E );
00815 #endif /* BAROTROPIC */
00816 #ifdef MHD
00817       Ul_x2Face[j][i].By -= hdtodx1*(rsf*x1Flux[j-1][i+1].Bz - lsf*x1Flux[j-1][i].Bz);
00818 #endif
00819 
00820       Ur_x2Face[j][i].d  -= hdtodx1*(rsf*x1Flux[j  ][i+1].d  - lsf*x1Flux[j  ][i].d );
00821       Ur_x2Face[j][i].Mx -= hdtodx1*(SQR(rsf)*x1Flux[j  ][i+1].My - SQR(lsf)*x1Flux[j  ][i].My);
00822       Ur_x2Face[j][i].My -= hdtodx1*(rsf*x1Flux[j  ][i+1].Mz - lsf*x1Flux[j  ][i].Mz);
00823       Ur_x2Face[j][i].Mz -= hdtodx1*(rsf*x1Flux[j  ][i+1].Mx - lsf*x1Flux[j  ][i].Mx);
00824 #ifndef BAROTROPIC
00825       Ur_x2Face[j][i].E  -= hdtodx1*(rsf*x1Flux[j  ][i+1].E  - lsf*x1Flux[j  ][i].E );
00826 #endif /* BAROTROPIC */
00827 #ifdef MHD
00828       Ur_x2Face[j][i].By -= hdtodx1*(rsf*x1Flux[j  ][i+1].Bz - lsf*x1Flux[j  ][i].Bz);
00829 #endif
00830 #if (NSCALARS > 0)
00831       for (n=0; n<NSCALARS; n++) {
00832       Ul_x2Face[j][i].s[n]-=hdtodx1*(rsf*x1Flux[j-1][i+1].s[n]-lsf*x1Flux[j-1][i].s[n]);
00833       Ur_x2Face[j][i].s[n]-=hdtodx1*(rsf*x1Flux[j  ][i+1].s[n]-lsf*x1Flux[j  ][i].s[n]);
00834       }
00835 #endif
00836     }
00837   }
00838 
00839 /*--- Step 6b: Not needed in 2D ---*/
00840 /*--- Step 6c ------------------------------------------------------------------
00841  * Add the "MHD source terms" from the x1-flux-gradients to the
00842  * conservative variables on the x2Face.
00843  */
00844 
00845 #ifdef MHD
00846   for (j=jl+1; j<=ju; j++) {
00847     for (i=il+1; i<=iu-1; i++) {
00848 #ifdef CYLINDRICAL
00849       hdtodx2 = hdt/(r[i]*pG->dx2);
00850 #endif
00851       dby = pG->B2i[ks][j][i] - pG->B2i[ks][j-1][i];
00852       B1 = pG->U[ks][j-1][i].B1c;
00853       B2 = pG->U[ks][j-1][i].B2c;
00854       B3 = pG->U[ks][j-1][i].B3c;
00855       V3 = pG->U[ks][j-1][i].M3/pG->U[ks][j-1][i].d;
00856 
00857       Ul_x2Face[j][i].Mz += hdtodx2*B1*dby;
00858       Ul_x2Face[j][i].Mx += hdtodx2*B2*dby;
00859       Ul_x2Face[j][i].My += hdtodx2*B3*dby;
00860       Ul_x2Face[j][i].By += hdtodx2*V3*dby;
00861 #ifndef BAROTROPIC
00862       Ul_x2Face[j][i].E  += hdtodx2*B3*V3*dby;
00863 #endif /* BAROTROPIC */
00864 
00865       dby = pG->B2i[ks][j+1][i] - pG->B2i[ks][j][i];
00866       B1 = pG->U[ks][j][i].B1c;
00867       B2 = pG->U[ks][j][i].B2c;
00868       B3 = pG->U[ks][j][i].B3c;
00869       V3 = pG->U[ks][j][i].M3/pG->U[ks][j][i].d;
00870 
00871       Ur_x2Face[j][i].Mz += hdtodx2*B1*dby;
00872       Ur_x2Face[j][i].Mx += hdtodx2*B2*dby;
00873       Ur_x2Face[j][i].My += hdtodx2*B3*dby;
00874       Ur_x2Face[j][i].By += hdtodx2*V3*dby;
00875 #ifndef BAROTROPIC
00876       Ur_x2Face[j][i].E  += hdtodx2*B3*V3*dby;
00877 #endif /* BAROTROPIC */
00878     }
00879   }
00880 #endif /* MHD */
00881 
00882 /*--- Step 6d ------------------------------------------------------------------
00883  * Add source terms for a static gravitational potential arising from x1-Flux
00884  * gradients. To improve conservation of total energy, average the energy
00885  * source term computed at cell faces.
00886  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00887  */
00888 
00889   if (StaticGravPot != NULL){
00890     for (j=jl+1; j<=ju; j++) {
00891       for (i=il+1; i<=iu-1; i++) {
00892         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00893         phic = (*StaticGravPot)((x1            ),x2,x3);
00894         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00895         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00896 
00897 #ifdef CYLINDRICAL
00898         g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00899         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00900         Ur_x2Face[j][i].Mz -= hdt*pG->U[ks][j][i].d*g;
00901 #else
00902         Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00903 #endif
00904 #ifndef BAROTROPIC
00905         Ur_x2Face[j][i].E -= hdtodx1*(lsf*x1Flux[j  ][i  ].d*(phic - phil) +
00906                                       rsf*x1Flux[j  ][i+1].d*(phir - phic));
00907 #endif
00908 
00909         phic = (*StaticGravPot)((x1            ),(x2-pG->dx2),x3);
00910         phir = (*StaticGravPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
00911         phil = (*StaticGravPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
00912 
00913 #ifdef CYLINDRICAL
00914         g = (*x1GravAcc)(x1vc(pG,i),(x2-pG->dx2),x3);
00915         Ul_x2Face[j][i].Mz -= hdt*pG->U[ks][j-1][i].d*g;
00916 #else
00917         Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00918 #endif
00919 #ifndef BAROTROPIC
00920         Ul_x2Face[j][i].E -= hdtodx1*(lsf*x1Flux[j-1][i  ].d*(phic - phil) +
00921                                       rsf*x1Flux[j-1][i+1].d*(phir - phic));
00922 #endif
00923       }
00924     }
00925   }
00926 
00927 /*--- Step 6d (cont) -----------------------------------------------------------
00928  * Add source terms for self gravity arising from x1-Flux gradients
00929  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00930  */
00931 
00932 #ifdef SELF_GRAVITY
00933   for (j=jl+1; j<=ju; j++) {
00934     for (i=il+1; i<=iu-1; i++) {
00935       phic = pG->Phi[ks][j][i];
00936       phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i+1]);
00937       phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i-1]);
00938 
00939       Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00940 #ifndef BAROTROPIC
00941       Ur_x2Face[j][i].E -= hdtodx1*(x1Flux[j  ][i  ].d*(phic - phil) +
00942                                     x1Flux[j  ][i+1].d*(phir - phic));
00943 #endif
00944 
00945       phic = pG->Phi[ks][j-1][i];
00946       phir = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j-1][i+1]);
00947       phil = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j-1][i-1]);
00948 
00949       Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00950 #ifndef BAROTROPIC
00951       Ul_x2Face[j][i].E -= hdtodx1*(x1Flux[j-1][i  ].d*(phic - phil) +
00952                                     x1Flux[j-1][i+1].d*(phir - phic));
00953 #endif
00954     }
00955   }
00956 
00957 #endif /* SELF_GRAVITY */
00958 
00959 /*--- Step 6d (cont) -----------------------------------------------------------
00960  * Add source terms for shearing box (Coriolis forces) for 0.5*dt arising from
00961  * x1-Flux gradient.  The tidal gravity terms are added via ShearingBoxPot
00962  *    Vx source term is (dt/2)( 2 Omega_0 Vy)
00963  *    Vy source term is (dt/2)(-2 Omega_0 Vx)
00964  *    Vy source term is (dt/2)((q-2) Omega_0 Vx) (with FARGO)
00965  * (x1,x2,x3) in code = (X,Z,Y) in shearing sheet
00966  */
00967 
00968 #ifdef SHEARING_BOX
00969   if (ShearingBoxPot != NULL){
00970     for (j=jl+1; j<=ju; j++) {
00971       for (i=il+1; i<=iu-1; i++) {
00972         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00973         phic = (*ShearingBoxPot)((x1            ),x2,x3);
00974         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
00975         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00976 
00977         Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00978 #ifndef BAROTROPIC
00979         Ur_x2Face[j][i].E -= hdtodx1*(x1Flux[j  ][i  ].d*(phic - phil) +
00980                                       x1Flux[j  ][i+1].d*(phir - phic));
00981 #endif
00982 
00983         phic = (*ShearingBoxPot)((x1            ),(x2-pG->dx2),x3);
00984         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
00985         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
00986 
00987         Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00988 #ifndef BAROTROPIC
00989         Ul_x2Face[j][i].E -= hdtodx1*(x1Flux[j-1][i  ].d*(phic - phil) +
00990                                       x1Flux[j-1][i+1].d*(phir - phic));
00991 #endif
00992       }
00993     }
00994   }
00995 
00996   if (ShBoxCoord == xz){
00997     for (j=jl+1; j<=ju; j++) {
00998       for (i=il+1; i<=iu-1; i++) {
00999         Ur_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j][i].M3;
01000         Ul_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j-1][i].M3;
01001 #ifdef FARGO
01002         Ur_x2Face[j][i].My += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01003         Ul_x2Face[j][i].My += hdt*(qshear-2.)*Omega_0*pG->U[ks][j-1][i].M1;
01004 #else
01005         Ur_x2Face[j][i].My -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01006         Ul_x2Face[j][i].My -= pG->dt*Omega_0*pG->U[ks][j-1][i].M1;
01007 #endif
01008       }
01009     }
01010   }
01011 
01012   if (ShBoxCoord == xy){
01013     for (j=jl+1; j<=ju; j++) {
01014       for (i=il+1; i<=iu-1; i++) {
01015         Ur_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j][i].M2;
01016         Ul_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j-1][i].M2;
01017 #ifdef FARGO
01018         Ur_x2Face[j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01019         Ul_x2Face[j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[ks][j-1][i].M1;
01020 #else
01021         Ur_x2Face[j][i].Mx -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01022         Ul_x2Face[j][i].Mx -= pG->dt*Omega_0*pG->U[ks][j-1][i].M1;
01023 #endif
01024       }
01025     }
01026   }
01027 #endif /* SHEARING_BOX */
01028 
01029 /*--- Step 6d (cont) -----------------------------------------------------------
01030  * ADD THE GEOMETRIC SOURCE TERMS IN THE x1-DIRECTION FOR dt/2
01031  */
01032 #ifdef CYLINDRICAL
01033   for (j=js-1; j<=ju; j++) {
01034     for (i=is-1; i<=ie+1; i++) {
01035       Ur_x2Face[j][i].Mz += hdt*geom_src[j  ][i];
01036       Ul_x2Face[j][i].Mz += hdt*geom_src[j-1][i];
01037     }
01038   }
01039 #endif
01040 
01041 /*=== STEP 7: Not needed in 2D ===*/
01042 
01043 /*=== STEP 8: Compute cell-centered values at n+1/2 ==========================*/
01044 
01045 /*--- Step 8a ------------------------------------------------------------------
01046  * Calculate d^{n+1/2} (needed with static potential, cooling, or MHD)
01047  */
01048 
01049 #ifndef MHD
01050 #ifndef PARTICLES
01051   if ((StaticGravPot != NULL) || (CoolingFunc != NULL)) 
01052 #endif
01053 #endif
01054   {
01055     for (j=jl+1; j<=ju-1; j++) {
01056       for (i=il+1; i<=iu-1; i++) {
01057 #ifdef CYLINDRICAL
01058         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01059         hdtodx2 = hdt/(r[i]*pG->dx2);
01060 #endif
01061         dhalf[j][i] = pG->U[ks][j][i].d
01062           - hdtodx1*(rsf*x1Flux[j  ][i+1].d - lsf*x1Flux[j][i].d)
01063           - hdtodx2*(    x2Flux[j+1][i  ].d -     x2Flux[j][i].d);
01064 #ifdef PARTICLES
01065         pG->Coup[ks][j][i].grid_d = dhalf[j][i];
01066 #endif
01067       }
01068     }
01069   }
01070 
01071 /*--- Step 8b ------------------------------------------------------------------
01072  * Calculate P^{n+1/2} (needed with cooling), and cell centered emf3_cc^{n+1/2}
01073  */
01074 
01075 #ifndef MHD
01076 #ifndef PARTICLES
01077   if (CoolingFunc != NULL) 
01078 #endif /* PARTICLES */
01079 #endif /* MHD */
01080   {
01081   for (j=jl+1; j<=ju-1; j++) {
01082     for (i=il+1; i<=iu-1; i++) {
01083 #ifdef CYLINDRICAL
01084       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01085       hdtodx2 = hdt/(r[i]*pG->dx2);
01086 #endif
01087       M1h = pG->U[ks][j][i].M1
01088         - hdtodx1*(rsf*x1Flux[j][i+1].Mx - lsf*x1Flux[j][i].Mx)
01089         - hdtodx2*(    x2Flux[j+1][i].Mz -     x2Flux[j][i].Mz);
01090 
01091       M2h = pG->U[ks][j][i].M2
01092         - hdtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My)
01093         - hdtodx2*(         x2Flux[j+1][i].Mx -          x2Flux[j][i].Mx);
01094 
01095       M3h = pG->U[ks][j][i].M3
01096         - hdtodx1*(rsf*x1Flux[j][i+1].Mz - lsf*x1Flux[j][i].Mz)
01097         - hdtodx2*(    x2Flux[j+1][i].My -     x2Flux[j][i].My);
01098 
01099 #ifndef BAROTROPIC
01100       Eh = pG->U[ks][j][i].E
01101         - hdtodx1*(rsf*x1Flux[j][i+1].E - lsf*x1Flux[j][i].E)
01102         - hdtodx2*(    x2Flux[j+1][i].E -     x2Flux[j][i].E);
01103 #endif
01104 
01105 /* Add source terms for fixed gravitational potential */
01106       if (StaticGravPot != NULL){
01107         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01108 #ifdef CYLINDRICAL
01109         g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01110         M1h -= hdt*pG->U[ks][j][i].d*g;
01111 #else
01112         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01113         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01114         M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01115 #endif
01116         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01117         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01118         M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01119       }
01120 
01121 /* Add source terms due to self-gravity  */
01122 #ifdef SELF_GRAVITY
01123       phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i+1]);
01124       phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i-1]);
01125       M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01126 
01127       phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j+1][i]);
01128       phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j-1][i]);
01129       M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01130 #endif /* SELF_GRAVITY */
01131 
01132 /* Add the tidal gravity and Coriolis terms for shearing box. */
01133 #ifdef SHEARING_BOX
01134       if (ShearingBoxPot != NULL){
01135         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01136         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01137         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01138         M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01139 
01140         phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
01141         phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
01142         M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01143       }
01144 
01145       if (ShBoxCoord == xy) M1h += pG->dt*Omega_0*pG->U[ks][j][i].M2;
01146       if (ShBoxCoord == xz) M1h += pG->dt*Omega_0*pG->U[ks][j][i].M3;
01147 #ifdef FARGO
01148       if (ShBoxCoord == xy) M2h += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01149       if (ShBoxCoord == xz) M3h += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01150 #else
01151       if (ShBoxCoord == xy) M2h -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01152       if (ShBoxCoord == xz) M3h -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01153 #endif
01154 #endif /* SHEARING_BOX */
01155 
01156 /* Add the particle feedback terms */
01157 #ifdef FEEDBACK
01158       M1h -= pG->Coup[ks][j][i].fb1;
01159       M2h -= pG->Coup[ks][j][i].fb2;
01160       M3h -= pG->Coup[ks][j][i].fb3;
01161 #endif /* FEEDBACK */
01162 
01163 
01164 /* Add the geometric source term */
01165 #ifdef CYLINDRICAL
01166       M1h += hdt*geom_src[j][i];
01167 #endif
01168 
01169 #ifndef BAROTROPIC
01170       phalf[j][i] = Eh - 0.5*(M1h*M1h + M2h*M2h + M3h*M3h)/dhalf[j][i];
01171 #endif
01172 
01173 #ifdef MHD
01174       B1ch = 0.5*(lsf*B1_x1Face[j][i] + rsf*B1_x1Face[j][i+1]);
01175       B2ch = 0.5*(    B2_x2Face[j][i] +     B2_x2Face[j+1][i]);
01176       B3ch = pG->U[ks][j][i].B3c 
01177         - hdtodx1*(rsf*x1Flux[j][i+1].Bz - lsf*x1Flux[j][i].Bz)
01178         - hdtodx2*(    x2Flux[j+1][i].By -     x2Flux[j][i].By);
01179       emf3_cc[j][i] = (B1ch*M2h - B2ch*M1h)/dhalf[j][i];
01180 
01181 #ifndef BAROTROPIC
01182       phalf[j][i] -= 0.5*(B1ch*B1ch + B2ch*B2ch + B3ch*B3ch);
01183 #endif
01184 #endif /* MHD */
01185 
01186 #ifndef BAROTROPIC
01187       phalf[j][i] *= Gamma_1;
01188 #endif
01189 
01190 #ifdef PARTICLES
01191       d1 = 1.0/dhalf[j][i];
01192       pG->Coup[ks][j][i].grid_v1 = M1h*d1;
01193       pG->Coup[ks][j][i].grid_v2 = M2h*d1;
01194       pG->Coup[ks][j][i].grid_v3 = M3h*d1;
01195 #ifndef BAROTROPIC
01196       pG->Coup[ks][j][i].grid_cs = sqrt(Gamma*phalf[j][i]*d1);
01197 #endif  /* BAROTROPIC */
01198 #endif /* PARTICLES */
01199     }
01200   }
01201   }
01202 
01203 /*=== STEP 8.5: Integrate the particles, compute the feedback ================*/
01204 
01205 #ifdef PARTICLES
01206   Integrate_Particles(pG,pD);
01207 #ifdef FEEDBACK
01208   exchange_feedback(pG, pD);
01209 #endif
01210 #endif
01211 
01212 /*=== STEP 9: Compute 2D x1-Flux, x2-Flux ====================================*/
01213 
01214 /*--- Step 9a ------------------------------------------------------------------
01215  * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
01216  *  al. (1998)) for H-correction
01217  */
01218 
01219 #ifdef H_CORRECTION
01220   for (j=js-1; j<=je+1; j++) {
01221     for (i=is-1; i<=ie+2; i++) {
01222 #ifdef MHD
01223       Bx = B1_x1Face[j][i];
01224 #endif
01225       cfr = cfast(&(Ur_x1Face[j][i]),&Bx);
01226       cfl = cfast(&(Ul_x1Face[j][i]),&Bx);
01227       lambdar = Ur_x1Face[j][i].Mx/Ur_x1Face[j][i].d + cfr;
01228       lambdal = Ul_x1Face[j][i].Mx/Ul_x1Face[j][i].d - cfl;
01229       eta1[j][i] = 0.5*fabs(lambdar - lambdal);
01230     }
01231   }
01232 
01233   for (j=js-1; j<=je+2; j++) {
01234     for (i=is-1; i<=ie+1; i++) {
01235 #ifdef MHD
01236       Bx = B2_x2Face[j][i];
01237 #endif
01238       cfr = cfast(&(Ur_x2Face[j][i]),&Bx);
01239       cfl = cfast(&(Ul_x2Face[j][i]),&Bx);
01240       lambdar = Ur_x2Face[j][i].Mx/Ur_x2Face[j][i].d + cfr;
01241       lambdal = Ul_x2Face[j][i].Mx/Ul_x2Face[j][i].d - cfl;
01242       eta2[j][i] = 0.5*fabs(lambdar - lambdal);
01243     }
01244   }
01245 #endif /* H_CORRECTION */
01246 
01247 /*--- Step 9b ------------------------------------------------------------------
01248  * Compute 2D x1-fluxes from corrected L/R states.
01249  */
01250 
01251   for (j=js-1; j<=je+1; j++) {
01252     for (i=is; i<=ie+1; i++) {
01253 #ifdef H_CORRECTION
01254       etah = MAX(eta2[j][i-1],eta2[j][i]);
01255       etah = MAX(etah,eta2[j+1][i-1]);
01256       etah = MAX(etah,eta2[j+1][i  ]);
01257       etah = MAX(etah,eta1[j  ][i  ]);
01258 #endif /* H_CORRECTION */
01259 #ifdef MHD
01260       Bx = B1_x1Face[j][i];
01261 #endif
01262       Wl[i] = Cons1D_to_Prim1D(&Ul_x1Face[j][i],&Bx);
01263       Wr[i] = Cons1D_to_Prim1D(&Ur_x1Face[j][i],&Bx);
01264 
01265       fluxes(Ul_x1Face[j][i],Ur_x1Face[j][i],Wl[i],Wr[i],Bx,&x1Flux[j][i]);
01266     }
01267   }
01268 
01269 /*--- Step 9c ------------------------------------------------------------------
01270  * Compute 2D x2-fluxes from corrected L/R states.
01271  */
01272 
01273   for (j=js; j<=je+1; j++) {
01274     for (i=is-1; i<=ie+1; i++) {
01275 #ifdef H_CORRECTION
01276       etah = MAX(eta1[j-1][i],eta1[j][i]);
01277       etah = MAX(etah,eta1[j-1][i+1]);
01278       etah = MAX(etah,eta1[j  ][i+1]);
01279       etah = MAX(etah,eta2[j  ][i  ]);
01280 #endif /* H_CORRECTION */
01281 #ifdef MHD
01282       Bx = B2_x2Face[j][i];
01283 #endif
01284       Wl[i] = Cons1D_to_Prim1D(&Ul_x2Face[j][i],&Bx);
01285       Wr[i] = Cons1D_to_Prim1D(&Ur_x2Face[j][i],&Bx);
01286 
01287       fluxes(Ul_x2Face[j][i],Ur_x2Face[j][i],Wl[i],Wr[i],Bx,&x2Flux[j][i]);
01288     }
01289   }
01290 
01291 /*=== STEP 10: Update face-centered B for a full timestep ====================*/
01292 
01293 /*--- Step 10a -----------------------------------------------------------------
01294  * Integrate emf*^{n+1/2} to the grid cell corners
01295  */
01296 
01297 #ifdef MHD
01298   integrate_emf3_corner(pG);
01299 
01300 /*--- Step 10b -----------------------------------------------------------------
01301  * Update the interface magnetic fields using CT for a full time step.
01302  */
01303 
01304   for (j=js; j<=je; j++) {
01305     for (i=is; i<=ie; i++) {
01306 #ifdef CYLINDRICAL
01307       dtodx2 = pG->dt/(ri[i]*pG->dx2);
01308 #endif
01309       pG->B1i[ks][j][i] -= dtodx2*(emf3[j+1][i  ] - emf3[j][i]);
01310       pG->B2i[ks][j][i] += dtodx1*(emf3[j  ][i+1] - emf3[j][i]);
01311     }
01312 #ifdef CYLINDRICAL
01313     dtodx2 = pG->dt/(ri[ie+1]*pG->dx2);
01314 #endif
01315     pG->B1i[ks][j][ie+1] -= dtodx2*(emf3[j+1][ie+1] - emf3[j][ie+1]);
01316   }
01317   for (i=is; i<=ie; i++) {
01318     pG->B2i[ks][je+1][i] += dtodx1*(emf3[je+1][i+1] - emf3[je+1][i]);
01319   }
01320 #endif /* MHD */
01321 
01322 /*=== STEP 11: Add source terms for a full timestep using n+1/2 states =======*/
01323 
01324 /*--- Step 11a -----------------------------------------------------------------
01325  * ADD GEOMETRIC SOURCE TERMS.
01326  */
01327 
01328 #ifdef CYLINDRICAL
01329   for (j=js; j<=je; j++) {
01330     for (i=is; i<=ie; i++) {
01331       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01332       hdtodx2 = hdt/(r[i]*pG->dx2);
01333 
01334       /* CALCULATE d AT TIME n+1/2 */
01335       dhalf[j][i] = pG->U[ks][j][i].d
01336         - hdtodx1*(rsf*x1Flux[j  ][i+1].d - lsf*x1Flux[j][i].d)
01337         - hdtodx2*(    x2Flux[j+1][i  ].d -     x2Flux[j][i].d);
01338 
01339       /* CALCULATE M2 AT TIME n+1/2 */
01340       M2h = pG->U[ks][j][i].M2
01341         - hdtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My)
01342         - hdtodx2*(         x2Flux[j+1][i].Mx -          x2Flux[j][i].Mx);
01343 
01344       /* ADD SOURCE TERM FOR FIXED GRAVITATIONAL POTENTIAL FOR 0.5*dt */
01345       if (StaticGravPot != NULL){
01346         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01347         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01348         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01349         M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01350       }
01351 
01352       /* COMPUTE GEOMETRIC SOURCE TERM AT TIME n+1/2 */
01353       geom_src[j][i] = SQR(M2h)/dhalf[j][i];
01354 #ifdef MHD
01355       B2ch = 0.5*(B2_x2Face[j][i] + B2_x2Face[j+1][i]);
01356       geom_src[j][i] -= SQR(B2ch);
01357 #endif
01358 #ifdef ISOTHERMAL
01359       geom_src[j][i] += Iso_csound2*dhalf[j][i];
01360 #ifdef MHD
01361       B1ch = 0.5*(rsf*B1_x1Face[j][i+1] + lsf*B1_x1Face[j][i]);
01362       B3ch = pG->U[ks][j][i].B3c
01363         - hdtodx1*(rsf*x1Flux[j  ][i+1].Bz - lsf*x1Flux[j][i].Bz)
01364         - hdtodx2*(    x2Flux[j+1][i  ].By -     x2Flux[j][i].By);
01365       geom_src[j][i] += 0.5*(SQR(B1ch)+SQR(B2ch)+SQR(B3ch));
01366 #endif
01367 #else /* ISOTHERMAL */
01368       Pavgh = 0.5*(lsf*x1Flux[j][i].Pflux + rsf*x1Flux[j][i+1].Pflux);
01369       geom_src[j][i] += Pavgh;
01370 #endif
01371       geom_src[j][i] /= x1vc(pG,i);
01372 
01373       /* ADD TIME-CENTERED GEOMETRIC SOURCE TERM FOR FULL dt */
01374       pG->U[ks][j][i].M1 += pG->dt*geom_src[j][i];
01375     }
01376   }
01377 #endif /* CYLINDRICAL */
01378 
01379 /*--- Step 11a -----------------------------------------------------------------
01380  * Add gravitational (or shearing box) source terms.
01381  * Remember with the 2D shearing box (1,2,3) = (x,z,y)
01382  *   A Crank-Nicholson update is used for shearing box terms.
01383  *   The energy source terms computed at cell faces are averaged to improve
01384  * conservation of total energy.
01385  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
01386  */
01387 
01388 #ifdef SHEARING_BOX
01389   fact = om_dt/(2. + (2.-qshear)*om_dt*om_dt);
01390   qom = qshear*Omega_0;
01391   for(j=js; j<=je; j++) {
01392     for(i=is; i<=ie; i++) {
01393       cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01394 
01395 /* Store the current state */
01396       M1n  = pG->U[ks][j][i].M1;
01397 #ifdef FARGO
01398       if (ShBoxCoord==xy) dM2n = pG->U[ks][j][i].M2;
01399       if (ShBoxCoord==xz) dM3n = pG->U[ks][j][i].M3;
01400 #else
01401       if (ShBoxCoord==xy) dM2n = pG->U[ks][j][i].M2 + qom*x1*pG->U[ks][j][i].d;
01402       if (ShBoxCoord==xz) dM3n = pG->U[ks][j][i].M3 + qom*x1*pG->U[ks][j][i].d;
01403 #endif
01404 
01405 /* Calculate the flux for the y-momentum fluctuation (M3 in 2D) */
01406       if (ShBoxCoord==xy){
01407         frx1_dM2 = x1Flux[j][i+1].My;
01408         flx1_dM2 = x1Flux[j][i  ].My;
01409         frx2_dM2 = x2Flux[j+1][i].Mx;
01410         flx2_dM2 = x2Flux[j  ][i].Mx;
01411       }
01412       if (ShBoxCoord==xz){
01413         frx1_dM3 = x1Flux[j][i+1].Mz;
01414         flx1_dM3 = x1Flux[j][i  ].Mz;
01415         frx2_dM3 = x2Flux[j+1][i].My;
01416         flx2_dM3 = x2Flux[j  ][i].My;
01417       }
01418 #ifndef FARGO
01419       if (ShBoxCoord==xy){
01420         frx1_dM2 += qom*(x1+0.5*pG->dx1)*x1Flux[j][i+1].d;
01421         flx1_dM2 += qom*(x1-0.5*pG->dx1)*x1Flux[j][i  ].d;
01422         frx2_dM2 += qom*(x1            )*x2Flux[j+1][i].d;
01423         flx2_dM2 += qom*(x1            )*x2Flux[j  ][i].d;
01424       }
01425       if (ShBoxCoord==xz){
01426         frx1_dM3 += qom*(x1+0.5*pG->dx1)*x1Flux[j][i+1].d;
01427         flx1_dM3 += qom*(x1-0.5*pG->dx1)*x1Flux[j][i  ].d;
01428         frx2_dM3 += qom*(x1            )*x2Flux[j+1][i].d;
01429         flx2_dM3 += qom*(x1            )*x2Flux[j  ][i].d;
01430       }
01431 #endif
01432 
01433 /* evolve M1n and dM3n by dt/2 using flux gradients */
01434       M1e = M1n - hdtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx)
01435                 - hdtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
01436 
01437       if (ShBoxCoord==xy){
01438         dM2e = dM2n - hdtodx1*(frx1_dM2 - flx1_dM2)
01439                     - hdtodx2*(frx2_dM2 - flx2_dM2);
01440       }
01441       if (ShBoxCoord==xz){
01442         dM3e = dM3n - hdtodx1*(frx1_dM3 - flx1_dM3)
01443                     - hdtodx2*(frx2_dM3 - flx2_dM3);
01444       }
01445 
01446 #ifdef FEEDBACK
01447       M1e -= 0.5*pG->Coup[ks][j][i].fb1;
01448       if (ShBoxCoord==xy) dM2e -= 0.5*pG->Coup[ks][j][i].fb2;
01449       if (ShBoxCoord==xz) dM3e -= 0.5*pG->Coup[ks][j][i].fb3;
01450 #endif
01451 
01452 /* Update the 1- and 2-momentum (or 1- and 3-momentum in XZ 2D shearing box)
01453  * for the Coriolis and tidal potential source terms using a Crank-Nicholson
01454  * discretization for the momentum fluctuation equation. */
01455 
01456       if (ShBoxCoord==xy){
01457         pG->U[ks][j][i].M1 += (4.0*dM2e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
01458         pG->U[ks][j][i].M2 += 2.0*(qshear-2.)*(M1e + om_dt*dM2e)*fact;
01459 #ifndef FARGO
01460         pG->U[ks][j][i].M2 -=0.5*qshear*om_dt*(x1Flux[j][i].d+x1Flux[j][i+1].d);
01461 #endif
01462       }
01463       if (ShBoxCoord==xz){
01464         pG->U[ks][j][i].M1 += (4.0*dM3e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
01465         pG->U[ks][j][i].M3 += 2.0*(qshear-2.)*(M1e + om_dt*dM3e)*fact;
01466 #ifndef FARGO
01467         pG->U[ks][j][i].M3 -=0.5*qshear*om_dt*(x1Flux[j][i].d+x1Flux[j][i+1].d);
01468 #endif
01469       }
01470 
01471 /* Update the energy for a fixed potential.
01472  * This update is identical to non-SHEARING_BOX below  */
01473 
01474       phic = (*ShearingBoxPot)((x1            ),x2,x3);
01475       phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01476       phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01477 #ifndef BAROTROPIC
01478       pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i  ].d*(phic - phil) +
01479                                    x1Flux[j][i+1].d*(phir - phic));
01480 #endif
01481 
01482       phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
01483       phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
01484 #ifndef BAROTROPIC
01485       pG->U[ks][j][i].E -= dtodx2*(x2Flux[j  ][i].d*(phic - phil) +
01486                                    x2Flux[j+1][i].d*(phir - phic));
01487 #endif
01488     }
01489   }
01490 
01491 #endif /* SHEARING_BOX */
01492 
01493   if (StaticGravPot != NULL){
01494     for (j=js; j<=je; j++) {
01495       for (i=is; i<=ie; i++) {
01496         cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01497         phic = (*StaticGravPot)((x1            ),x2,x3);
01498         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01499         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01500 
01501 #ifdef CYLINDRICAL
01502         g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01503         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01504         dtodx2 = pG->dt/(r[i]*pG->dx2);
01505         pG->U[ks][j][i].M1 -= pG->dt*dhalf[j][i]*g;
01506 #else
01507         pG->U[ks][j][i].M1 -= dtodx1*dhalf[j][i]*(phir-phil);
01508 #endif
01509 
01510 #ifndef BAROTROPIC
01511         pG->U[ks][j][i].E -= dtodx1*(lsf*x1Flux[j][i  ].d*(phic - phil) +
01512                                      rsf*x1Flux[j][i+1].d*(phir - phic));
01513 #endif
01514         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01515         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01516 
01517         pG->U[ks][j][i].M2 -= dtodx2*dhalf[j][i]*(phir-phil);
01518 
01519 #ifndef BAROTROPIC
01520         pG->U[ks][j][i].E -= dtodx2*(x2Flux[j  ][i].d*(phic - phil) +
01521                                      x2Flux[j+1][i].d*(phir - phic));
01522 #endif
01523       }
01524     }
01525   }
01526 
01527 
01528 /*--- Step 11b -----------------------------------------------------------------
01529  * Add source terms for self-gravity.
01530  * A flux correction using Phi^{n+1} in the main loop is required to make
01531  * the source terms 2nd order: see selfg_flux_correction().
01532  */
01533 
01534 #ifdef SELF_GRAVITY
01535 /* Add fluxes and source terms due to (d/dx1) terms  */
01536 
01537   for (j=js; j<=je; j++){
01538     for (i=is; i<=ie; i++){
01539       phic = pG->Phi[ks][j][i];
01540       phil = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j][i  ]);
01541       phir = 0.5*(pG->Phi[ks][j][i  ] + pG->Phi[ks][j][i+1]);
01542 
01543 /* gx and gy centered at L and R x1-faces */
01544       gxl = (pG->Phi[ks][j][i-1] - pG->Phi[ks][j][i  ])/(pG->dx1);
01545       gxr = (pG->Phi[ks][j][i  ] - pG->Phi[ks][j][i+1])/(pG->dx1);
01546 
01547       gyl = 0.25*((pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j+1][i-1]) +
01548                   (pG->Phi[ks][j-1][i  ] - pG->Phi[ks][j+1][i  ]) )/(pG->dx2);
01549       gyr = 0.25*((pG->Phi[ks][j-1][i  ] - pG->Phi[ks][j+1][i  ]) +
01550                   (pG->Phi[ks][j-1][i+1] - pG->Phi[ks][j+1][i+1]) )/(pG->dx2);
01551 
01552 /* momentum fluxes in x1.  2nd term is needed only if Jean's swindle used */
01553       flux_m1l = 0.5*(gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
01554       flux_m1r = 0.5*(gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
01555 
01556       flux_m2l = gxl*gyl/four_pi_G;
01557       flux_m2r = gxr*gyr/four_pi_G;
01558 
01559 /* Update momenta and energy with d/dx1 terms  */
01560       pG->U[ks][j][i].M1 -= dtodx1*(flux_m1r - flux_m1l);
01561       pG->U[ks][j][i].M2 -= dtodx1*(flux_m2r - flux_m2l);
01562 #ifndef BAROTROPIC
01563       pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i  ].d*(phic - phil) +
01564                                    x1Flux[j][i+1].d*(phir - phic));
01565 #endif
01566     }
01567   }
01568 
01569 /* Add fluxes and source terms due to (d/dx2) terms  */
01570 
01571   for (j=js; j<=je; j++){
01572     for (i=is; i<=ie; i++){
01573       phic = pG->Phi[ks][j][i];
01574       phil = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j  ][i]);
01575       phir = 0.5*(pG->Phi[ks][j  ][i] + pG->Phi[ks][j+1][i]);
01576 
01577 /* gx and gy centered at L and R x2-faces */
01578       gxl = 0.25*((pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j-1][i+1]) +
01579                   (pG->Phi[ks][j  ][i-1] - pG->Phi[ks][j  ][i+1]) )/(pG->dx1);
01580       gxr = 0.25*((pG->Phi[ks][j  ][i-1] - pG->Phi[ks][j  ][i+1]) +
01581                   (pG->Phi[ks][j+1][i-1] - pG->Phi[ks][j+1][i+1]) )/(pG->dx1);
01582 
01583       gyl = (pG->Phi[ks][j-1][i] - pG->Phi[ks][j  ][i])/(pG->dx2);
01584       gyr = (pG->Phi[ks][j  ][i] - pG->Phi[ks][j+1][i])/(pG->dx2);
01585 
01586 /* momentum fluxes in x2.  2nd term is needed only if Jean's swindle used */
01587       flux_m1l = gyl*gxl/four_pi_G;
01588       flux_m1r = gyr*gxr/four_pi_G;
01589 
01590       flux_m2l = 0.5*(gyl*gyl-gxl*gxl)/four_pi_G + grav_mean_rho*phil;
01591       flux_m2r = 0.5*(gyr*gyr-gxr*gxr)/four_pi_G + grav_mean_rho*phir;
01592 
01593 /* Update momenta and energy with d/dx2 terms  */
01594       pG->U[ks][j][i].M1 -= dtodx2*(flux_m1r - flux_m1l);
01595       pG->U[ks][j][i].M2 -= dtodx2*(flux_m2r - flux_m2l);
01596 #ifndef BAROTROPIC
01597       pG->U[ks][j][i].E -= dtodx2*(x2Flux[j  ][i].d*(phic - phil) +
01598                                    x2Flux[j+1][i].d*(phir - phic));
01599 #endif
01600     }
01601   }
01602 
01603 /* Save mass fluxes in Grid structure for source term correction in main loop */
01604   for (j=js; j<=je+1; j++) {
01605     for (i=is; i<=ie+1; i++) {
01606       pG->x1MassFlux[ks][j][i] = x1Flux[j][i].d;
01607       pG->x2MassFlux[ks][j][i] = x2Flux[j][i].d;
01608     }
01609   }
01610 #endif /* SELF_GRAVITY */
01611 
01612 /*--- Step 11c -----------------------------------------------------------------
01613  * Add source terms for optically thin cooling
01614  */
01615 
01616 #ifndef BAROTROPIC
01617   if (CoolingFunc != NULL){
01618     for (j=js; j<=je; j++){
01619       for (i=is; i<=ie; i++){
01620         coolf = (*CoolingFunc)(dhalf[j][i],phalf[j][i],pG->dt);
01621         pG->U[ks][j][i].E -= pG->dt*coolf;
01622       }
01623     }
01624   }
01625 #endif /* BAROTROPIC */
01626 
01627 /*--- Step 11d -----------------------------------------------------------------
01628  * Add source terms for particle feedback
01629  */
01630 
01631 #ifdef FEEDBACK
01632   for (j=js; j<=je; j++)
01633     for (i=is; i<=ie; i++) {
01634       pG->U[ks][j][i].M1 -= pG->Coup[ks][j][i].fb1;
01635       pG->U[ks][j][i].M2 -= pG->Coup[ks][j][i].fb2;
01636       pG->U[ks][j][i].M3 -= pG->Coup[ks][j][i].fb3;
01637 #ifndef BAROTROPIC
01638       pG->U[ks][j][i].E += pG->Coup[ks][j][i].Eloss;
01639       pG->Coup[ks][j][i].Eloss *= dt1; /* for history output purpose */
01640       pG->Eloss[ks][j][i] *= dt1; /* for history output purpose */
01641 #endif
01642     }
01643 #endif
01644 
01645 /*=== STEP 12: Update cell-centered values for a full timestep ===============*/
01646 
01647 /*--- Step 12a -----------------------------------------------------------------
01648  * Update cell-centered variables in pG using 2D x1-fluxes
01649  */
01650 
01651   for (j=js; j<=je; j++) {
01652     for (i=is; i<=ie; i++) {
01653 #ifdef CYLINDRICAL
01654       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01655 #endif
01656       pG->U[ks][j][i].d  -= dtodx1*(rsf*x1Flux[j][i+1].d  - lsf*x1Flux[j][i].d );
01657       pG->U[ks][j][i].M1 -= dtodx1*(rsf*x1Flux[j][i+1].Mx - lsf*x1Flux[j][i].Mx);
01658       pG->U[ks][j][i].M2 -= dtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My);
01659       pG->U[ks][j][i].M3 -= dtodx1*(rsf*x1Flux[j][i+1].Mz - lsf*x1Flux[j][i].Mz);
01660 #ifndef BAROTROPIC
01661       pG->U[ks][j][i].E  -= dtodx1*(rsf*x1Flux[j][i+1].E  - lsf*x1Flux[j][i].E );
01662 #endif /* BAROTROPIC */
01663 #ifdef MHD
01664       pG->U[ks][j][i].B2c -= dtodx1*(x1Flux[j][i+1].By - x1Flux[j][i].By);
01665       pG->U[ks][j][i].B3c -= dtodx1*(rsf*x1Flux[j][i+1].Bz - lsf*x1Flux[j][i].Bz);
01666 #endif /* MHD */
01667 #if (NSCALARS > 0)
01668       for (n=0; n<NSCALARS; n++)
01669         pG->U[ks][j][i].s[n] -= dtodx1*(rsf*x1Flux[j][i+1].s[n] 
01670                                       - lsf*x1Flux[j][i  ].s[n]);
01671 #endif
01672     }
01673   }
01674 
01675 /*--- Step 12b -----------------------------------------------------------------
01676  * Update cell-centered variables in pG using 2D x2-fluxes
01677  */
01678 
01679   for (j=js; j<=je; j++) {
01680     for (i=is; i<=ie; i++) {
01681 #ifdef CYLINDRICAL
01682       dtodx2 = pG->dt/(r[i]*pG->dx2);
01683 #endif
01684       pG->U[ks][j][i].d  -= dtodx2*(x2Flux[j+1][i].d  - x2Flux[j][i].d );
01685       pG->U[ks][j][i].M1 -= dtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
01686       pG->U[ks][j][i].M2 -= dtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
01687       pG->U[ks][j][i].M3 -= dtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
01688 #ifndef BAROTROPIC
01689       pG->U[ks][j][i].E  -= dtodx2*(x2Flux[j+1][i].E  - x2Flux[j][i].E );
01690 #endif /* BAROTROPIC */
01691 #ifdef MHD
01692       pG->U[ks][j][i].B3c -= dtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
01693       pG->U[ks][j][i].B1c -= dtodx2*(x2Flux[j+1][i].Bz - x2Flux[j][i].Bz);
01694 #endif /* MHD */
01695 #if (NSCALARS > 0)
01696       for (n=0; n<NSCALARS; n++)
01697         pG->U[ks][j][i].s[n] -= dtodx2*(x2Flux[j+1][i].s[n] 
01698                                          - x2Flux[j  ][i].s[n]);
01699 #endif
01700     }
01701   }
01702 
01703 /*--- Step 12c: Not needed in 2D ---*/
01704 /*--- Step 12d -----------------------------------------------------------------
01705  * Set cell centered magnetic fields to average of updated face centered fields.
01706  */
01707 
01708 #ifdef MHD
01709   for (j=js; j<=je; j++) {
01710     for (i=is; i<=ie; i++) {
01711 #ifdef CYLINDRICAL
01712       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01713 #endif
01714       pG->U[ks][j][i].B1c =0.5*(lsf*pG->B1i[ks][j][i] + rsf*pG->B1i[ks][j][i+1]);
01715       pG->U[ks][j][i].B2c =0.5*(    pG->B2i[ks][j][i] +     pG->B2i[ks][j+1][i]);
01716 /* Set the 3-interface magnetic field equal to the cell center field. */
01717       pG->B3i[ks][j][i] = pG->U[ks][j][i].B3c;
01718     }
01719   }
01720 #endif /* MHD */
01721 
01722 #ifdef STATIC_MESH_REFINEMENT
01723 /*--- Step 12e -----------------------------------------------------------------
01724  * With SMR, store fluxes at boundaries of child and parent grids.
01725  */
01726 
01727   for (ncg=0; ncg<pG->NCGrid; ncg++) {
01728 
01729 /* x1-boundaries of child Grids (interior to THIS Grid) */
01730 
01731     for (dim=0; dim<2; dim++){
01732       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01733 
01734         if (dim==0) i = pG->CGrid[ncg].ijks[0];
01735         if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
01736         jcs = pG->CGrid[ncg].ijks[1];
01737         jce = pG->CGrid[ncg].ijke[1];
01738 
01739         for (j=jcs, jj=0; j<=jce; j++, jj++){
01740           pG->CGrid[ncg].myFlx[dim][ks][jj].d  = x1Flux[j][i].d; 
01741           pG->CGrid[ncg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx; 
01742           pG->CGrid[ncg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01743           pG->CGrid[ncg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz; 
01744 #ifndef BAROTROPIC
01745           pG->CGrid[ncg].myFlx[dim][ks][jj].E  = x1Flux[j][i].E; 
01746 #endif /* BAROTROPIC */
01747 #ifdef MHD
01748           pG->CGrid[ncg].myFlx[dim][ks][jj].B1c = 0.0;
01749           pG->CGrid[ncg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By; 
01750           pG->CGrid[ncg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz; 
01751 #endif /* MHD */
01752 #if (NSCALARS > 0)
01753           for (n=0; n<NSCALARS; n++)
01754             pG->CGrid[ncg].myFlx[dim][ks][jj].s[n]  = x1Flux[j][i].s[n]; 
01755 #endif
01756         }
01757 #ifdef MHD
01758         for (j=jcs, jj=0; j<=jce+1; j++, jj++){
01759           pG->CGrid[ncg].myEMF3[dim][ks][jj] = emf3[j][i];
01760         }
01761 #endif /* MHD */
01762       }
01763     }
01764 
01765 /* x2-boundaries of child Grids (interior to THIS Grid) */
01766 
01767     for (dim=2; dim<4; dim++){
01768       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01769 
01770         ics = pG->CGrid[ncg].ijks[0];
01771         ice = pG->CGrid[ncg].ijke[0];
01772         if (dim==2) j = pG->CGrid[ncg].ijks[1];
01773         if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
01774 
01775         for (i=ics, ii=0; i<=ice; i++, ii++){
01776           pG->CGrid[ncg].myFlx[dim][ks][ii].d  = x2Flux[j][i].d; 
01777           pG->CGrid[ncg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz; 
01778           pG->CGrid[ncg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01779           pG->CGrid[ncg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My; 
01780 #ifndef BAROTROPIC
01781           pG->CGrid[ncg].myFlx[dim][ks][ii].E  = x2Flux[j][i].E; 
01782 #endif /* BAROTROPIC */
01783 #ifdef MHD
01784           pG->CGrid[ncg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz; 
01785           pG->CGrid[ncg].myFlx[dim][ks][ii].B2c = 0.0;
01786           pG->CGrid[ncg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By; 
01787 #endif /* MHD */
01788 #if (NSCALARS > 0)
01789           for (n=0; n<NSCALARS; n++)
01790             pG->CGrid[ncg].myFlx[dim][ks][ii].s[n]  = x2Flux[j][i].s[n]; 
01791 #endif
01792         }
01793 #ifdef MHD
01794         for (i=ics, ii=0; i<=ice+1; i++, ii++){
01795           pG->CGrid[ncg].myEMF3[dim][ks][ii] = emf3[j][i];
01796         }
01797 #endif /* MHD */
01798       }
01799     }
01800   }
01801 
01802   for (npg=0; npg<pG->NPGrid; npg++) {
01803 
01804 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
01805 
01806     for (dim=0; dim<2; dim++){
01807       if (pG->PGrid[npg].myFlx[dim] != NULL) {
01808 
01809         if (dim==0) i = pG->PGrid[npg].ijks[0];
01810         if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
01811         jps = pG->PGrid[npg].ijks[1];
01812         jpe = pG->PGrid[npg].ijke[1];
01813 
01814         for (j=jps, jj=0; j<=jpe; j++, jj++){
01815           pG->PGrid[npg].myFlx[dim][ks][jj].d  = x1Flux[j][i].d; 
01816           pG->PGrid[npg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx; 
01817           pG->PGrid[npg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01818           pG->PGrid[npg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz; 
01819 #ifndef BAROTROPIC
01820           pG->PGrid[npg].myFlx[dim][ks][jj].E  = x1Flux[j][i].E; 
01821 #endif /* BAROTROPIC */
01822 #ifdef MHD
01823           pG->PGrid[npg].myFlx[dim][ks][jj].B1c = 0.0;
01824           pG->PGrid[npg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By; 
01825           pG->PGrid[npg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz; 
01826 #endif /* MHD */
01827 #if (NSCALARS > 0)
01828           for (n=0; n<NSCALARS; n++)
01829             pG->PGrid[npg].myFlx[dim][ks][jj].s[n]  = x1Flux[j][i].s[n]; 
01830 #endif
01831         }
01832 #ifdef MHD
01833         for (j=jps, jj=0; j<=jpe+1; j++, jj++){
01834           pG->PGrid[npg].myEMF3[dim][ks][jj] = emf3[j][i];
01835         }
01836 #endif /* MHD */
01837       }
01838     }
01839 
01840 /* x2-boundaries of parent Grids (at boundaries of THIS Grid)  */
01841 
01842     for (dim=2; dim<4; dim++){
01843       if (pG->PGrid[npg].myFlx[dim] != NULL) {
01844 
01845         ips = pG->PGrid[npg].ijks[0];
01846         ipe = pG->PGrid[npg].ijke[0];
01847         if (dim==2) j = pG->PGrid[npg].ijks[1];
01848         if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
01849 
01850         for (i=ips, ii=0; i<=ipe; i++, ii++){
01851           pG->PGrid[npg].myFlx[dim][ks][ii].d  = x2Flux[j][i].d; 
01852           pG->PGrid[npg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz; 
01853           pG->PGrid[npg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01854           pG->PGrid[npg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My; 
01855 #ifndef BAROTROPIC
01856           pG->PGrid[npg].myFlx[dim][ks][ii].E  = x2Flux[j][i].E; 
01857 #endif /* BAROTROPIC */
01858 #ifdef MHD
01859           pG->PGrid[npg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz; 
01860           pG->PGrid[npg].myFlx[dim][ks][ii].B2c = 0.0;
01861           pG->PGrid[npg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By; 
01862 #endif /* MHD */
01863 #if (NSCALARS > 0)
01864           for (n=0; n<NSCALARS; n++)
01865             pG->PGrid[npg].myFlx[dim][ks][ii].s[n]  = x2Flux[j][i].s[n]; 
01866 #endif
01867         }
01868 #ifdef MHD
01869         for (i=ips, ii=0; i<=ipe+1; i++, ii++){
01870           pG->PGrid[npg].myEMF3[dim][ks][ii] = emf3[j][i];
01871         }
01872 #endif /* MHD */
01873       }
01874     }
01875   }
01876 
01877 #endif /* STATIC_MESH_REFINEMENT */
01878 
01879   return;
01880 }
01881 
01882 /*----------------------------------------------------------------------------*/
01883 /*! \fn void integrate_init_2d(MeshS *pM)
01884  *  \brief Allocate temporary integration arrays */
01885 void integrate_init_2d(MeshS *pM)
01886 {
01887   int nmax,size1=0,size2=0,nl,nd;
01888 
01889 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2 */
01890   for (nl=0; nl<(pM->NLevels); nl++){
01891     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01892       if (pM->Domain[nl][nd].Grid != NULL) {
01893         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01894           size1 = pM->Domain[nl][nd].Grid->Nx[0];
01895         }
01896         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01897           size2 = pM->Domain[nl][nd].Grid->Nx[1];
01898         }
01899       }
01900     }
01901   }
01902 
01903   size1 = size1 + 2*nghost;
01904   size2 = size2 + 2*nghost;
01905   nmax = MAX(size1,size2);
01906 
01907 #ifdef MHD
01908   if ((emf3 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01909     goto on_error;
01910   if ((emf3_cc = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01911     goto on_error;
01912 #endif /* MHD */
01913 #ifdef H_CORRECTION
01914   if ((eta1 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01915     goto on_error;
01916   if ((eta2 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01917     goto on_error;
01918 #endif /* H_CORRECTION */
01919 
01920   if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01921   if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01922 
01923 #ifdef MHD
01924   if ((B1_x1Face = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01925     goto on_error;
01926   if ((B2_x2Face = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01927     goto on_error;
01928 #endif /* MHD */
01929 
01930   if ((U1d= (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01931   if ((Ul = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01932   if ((Ur = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01933   if ((W  = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01934   if ((Wl = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01935   if ((Wr = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01936 
01937   if ((Ul_x1Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01938     goto on_error;
01939   if ((Ur_x1Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01940     goto on_error;
01941   if ((Ul_x2Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01942     goto on_error;
01943   if ((Ur_x2Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01944     goto on_error;
01945 
01946   if ((x1Flux   =(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01947     goto on_error;
01948   if ((x2Flux   =(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01949     goto on_error;
01950 
01951 #ifndef CYLINDRICAL
01952 #ifndef MHD
01953 #ifndef PARTICLES
01954   if((StaticGravPot != NULL) || (CoolingFunc != NULL))
01955 #endif
01956 #endif
01957 #endif
01958   {
01959   if ((dhalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01960     goto on_error;
01961   if ((phalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01962     goto on_error;
01963   }
01964 
01965   /* DATA STRUCTURES FOR CYLINDRICAL COORDINATES */
01966 #ifdef CYLINDRICAL
01967   if ((geom_src = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL) 
01968     goto on_error;
01969 #endif
01970 
01971   return;
01972 
01973   on_error:
01974     integrate_destruct();
01975     ath_error("[integrate_init]: malloc returned a NULL pointer\n");
01976 }
01977 
01978 /*----------------------------------------------------------------------------*/
01979 /*! \fn void integrate_destruct_2d(void)
01980  *  \brief Free temporary integration arrays */
01981 void integrate_destruct_2d(void)
01982 {
01983 #ifdef MHD
01984   if (emf3    != NULL) free_2d_array(emf3);
01985   if (emf3_cc != NULL) free_2d_array(emf3_cc);
01986 #endif /* MHD */
01987 #ifdef H_CORRECTION
01988   if (eta1 != NULL) free_2d_array(eta1);
01989   if (eta2 != NULL) free_2d_array(eta2);
01990 #endif /* H_CORRECTION */
01991   if (Bxc != NULL) free(Bxc);
01992   if (Bxi != NULL) free(Bxi);
01993 #ifdef MHD
01994   if (B1_x1Face != NULL) free_2d_array(B1_x1Face);
01995   if (B2_x2Face != NULL) free_2d_array(B2_x2Face);
01996 #endif /* MHD */
01997 
01998   if (U1d      != NULL) free(U1d);
01999   if (Ul       != NULL) free(Ul);
02000   if (Ur       != NULL) free(Ur);
02001   if (W        != NULL) free(W);
02002   if (Wl       != NULL) free(Wl);
02003   if (Wr       != NULL) free(Wr);
02004 
02005   if (Ul_x1Face != NULL) free_2d_array(Ul_x1Face);
02006   if (Ur_x1Face != NULL) free_2d_array(Ur_x1Face);
02007   if (Ul_x2Face != NULL) free_2d_array(Ul_x2Face);
02008   if (Ur_x2Face != NULL) free_2d_array(Ur_x2Face);
02009   if (x1Flux    != NULL) free_2d_array(x1Flux);
02010   if (x2Flux    != NULL) free_2d_array(x2Flux);
02011   if (dhalf     != NULL) free_2d_array(dhalf);
02012   if (phalf     != NULL) free_2d_array(phalf);
02013 
02014   /* DATA STRUCTURES FOR CYLINDRICAL COORDINATES */
02015 #ifdef CYLINDRICAL
02016   if (geom_src  != NULL) free_2d_array(geom_src);
02017 #endif
02018 
02019   return;
02020 }
02021 
02022 /*=========================== PRIVATE FUNCTIONS ==============================*/
02023 
02024 /*----------------------------------------------------------------------------*/
02025 /*! \fn static void integrate_emf3_corner(GridS *pG)
02026  *  \brief Corner EMF terms in the upwind CT method */
02027 #ifdef MHD
02028 static void integrate_emf3_corner(GridS *pG)
02029 {
02030   int i,is,ie,j,js,je;
02031   Real emf_l1, emf_r1, emf_l2, emf_r2;
02032   Real rsf=1.0,lsf=1.0;
02033 
02034   is = pG->is;   ie = pG->ie;
02035   js = pG->js;   je = pG->je;
02036 
02037 /* NOTE: The x1-Flux of B2 is -E3.  The x2-Flux of B1 is +E3. */
02038   for (j=js-1; j<=je+2; j++) {
02039     for (i=is-1; i<=ie+2; i++) {
02040 #ifdef CYLINDRICAL
02041       rsf = pG->ri[i]/pG->r[i];  lsf = pG->ri[i]/pG->r[i-1];
02042 #endif
02043       if (x1Flux[j-1][i].d > 0.0) {
02044         emf_l2 = -x1Flux[j-1][i].By
02045           + (x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1])*lsf;
02046       }
02047       else if (x1Flux[j-1][i].d < 0.0) {
02048         emf_l2 = -x1Flux[j-1][i].By
02049           + (x2Flux[j][i].Bz - emf3_cc[j-1][i])*rsf;
02050 
02051       } else {
02052         emf_l2 = -x1Flux[j-1][i].By
02053           + 0.5*((x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1])*lsf + 
02054                  (x2Flux[j][i  ].Bz - emf3_cc[j-1][i  ])*rsf );
02055       }
02056 
02057       if (x1Flux[j][i].d > 0.0) {
02058         emf_r2 = -x1Flux[j][i].By 
02059           + (x2Flux[j][i-1].Bz - emf3_cc[j][i-1])*lsf;
02060       }
02061       else if (x1Flux[j][i].d < 0.0) {
02062         emf_r2 = -x1Flux[j][i].By 
02063           + (x2Flux[j][i].Bz - emf3_cc[j][i])*rsf;
02064 
02065       } else {
02066         emf_r2 = -x1Flux[j][i].By 
02067           + 0.5*((x2Flux[j][i-1].Bz - emf3_cc[j][i-1])*lsf + 
02068                  (x2Flux[j][i  ].Bz - emf3_cc[j][i  ])*rsf );
02069       }
02070 
02071       if (x2Flux[j][i-1].d > 0.0) {
02072         emf_l1 = x2Flux[j][i-1].Bz
02073           + (-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]);
02074       }
02075       else if (x2Flux[j][i-1].d < 0.0) {
02076         emf_l1 = x2Flux[j][i-1].Bz 
02077           + (-x1Flux[j][i].By - emf3_cc[j][i-1]);
02078       } else {
02079         emf_l1 = x2Flux[j][i-1].Bz
02080           + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]
02081                  -x1Flux[j  ][i].By - emf3_cc[j  ][i-1] );
02082       }
02083 
02084       if (x2Flux[j][i].d > 0.0) {
02085         emf_r1 = x2Flux[j][i].Bz
02086           + (-x1Flux[j-1][i].By - emf3_cc[j-1][i]);
02087       }
02088       else if (x2Flux[j][i].d < 0.0) {
02089         emf_r1 = x2Flux[j][i].Bz
02090           + (-x1Flux[j][i].By - emf3_cc[j][i]);
02091       } else {
02092         emf_r1 = x2Flux[j][i].Bz
02093           + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i]
02094                  -x1Flux[j  ][i].By - emf3_cc[j  ][i] );
02095       }
02096 
02097       emf3[j][i] = 0.25*(emf_l1 + emf_r1 + emf_l2 + emf_r2);
02098     }
02099   }
02100 
02101   return;
02102 }
02103 #endif /* MHD */
02104 
02105 #endif /* CTU_INTEGRATOR */
integrators/integrate_2d_ctu.c
