Athena: integrators/integrate_1d_vl

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_1d_vl_sr.c
00004  *  \brief Integrate SRMHD equations using 1D version of MUSCL-Hancock (VL)
00005  *   integrator.
00006  *
00007  * PURPOSE: Integrate SRMHD equations using 1D version of MUSCL-Hancock (VL)
00008  *   integrator.  Updates U.[d,M1,M2,M3,E,B2c,B3c,s] in Grid structure.
00009  *   Adds gravitational source terms, self-gravity.
00010  *        
00011  * CONTAINS PUBLIC FUNCTIONS: 
00012  * - integrate_1d_vl()
00013  * - integrate_init_1d()
00014  * - integrate_destruct_1d() */
00015 /*============================================================================*/
00016 
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include "../defs.h"
00022 #include "../athena.h"
00023 #include "../globals.h"
00024 #include "prototypes.h"
00025 #include "../prototypes.h"
00026 
00027 #ifdef SPECIAL_RELATIVITY
00028 
00029 #if defined(VL_INTEGRATOR) && defined(CARTESIAN)
00030 
00031 /* The L/R states of primitive variables and fluxes at each cell face */
00032 static Prim1DS *Wl_x1Face=NULL, *Wr_x1Face=NULL;
00033 static Cons1DS *x1Flux=NULL;
00034 
00035 /* 1D scratch vectors used by lr_states and flux functions */
00036 static Real *Bxc=NULL, *Bxi=NULL;
00037 static Prim1DS *W1d=NULL, *Wl=NULL, *Wr=NULL;
00038 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00039 #ifdef FIRST_ORDER_FLUX_CORRECTION
00040 static Cons1DS *x1FluxP=NULL;
00041 #endif
00042 
00043 /* conserved & primitive variables at t^{n} */
00044 static PrimS *W=NULL;
00045 
00046 /* conserved variables at t^{n+1/2} computed in predict step */
00047 static ConsS *Uhalf=NULL;
00048 static PrimS *Whalf=NULL;
00049 
00050 #ifdef FIRST_ORDER_FLUX_CORRECTION
00051 static void FixCell(GridS *pG, Int3Vect);
00052 #endif
00053 
00054 
00055 /*=========================== PUBLIC FUNCTIONS ===============================*/
00056 /*----------------------------------------------------------------------------*/
00057 /*! \fn void integrate_1d_vl(DomainS *pD) 
00058  *  \brief 1D version of van Leer unsplit integrator for MHD. 
00059  *
00060  *   The numbering of steps follows the numbering in the 3D version.
00061  *   NOT ALL STEPS ARE NEEDED IN 1D.
00062  */
00063 void integrate_1d_vl(DomainS *pD)
00064 {
00065   GridS *pG=(pD->Grid);
00066   ConsS U;
00067   Real dtodx1=pG->dt/pG->dx1, hdtodx1=0.5*pG->dt/pG->dx1;
00068   int i, is = pG->is, ie = pG->ie;
00069   int js = pG->js;
00070   int ks = pG->ks;
00071   int cart_x1 = 1, cart_x2 = 2, cart_x3 = 3;
00072   Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic;
00073 #if (NSCALARS > 0)
00074   int n;
00075 #endif
00076 #ifdef SELF_GRAVITY
00077   Real gxl,gxr,flx_m1l,flx_m1r;
00078 #endif
00079 #ifdef STATIC_MESH_REFINEMENT
00080   int ncg,npg,dim;
00081 #endif
00082 #ifdef FIRST_ORDER_FLUX_CORRECTION
00083   int flag_cell=0,negd=0,negP=0,superl=0,NaNFlux=0;
00084   int fail=0,final=0;
00085   Real Vsq,Bx;
00086   PrimS Wcheck;
00087   Int3Vect BadCell;
00088 #endif
00089 
00090 
00091   int il=is-(nghost-1), iu=ie+(nghost-1);
00092 
00093   for (i=is-nghost; i<=ie+nghost; i++) {
00094     Uhalf[i] = pG->U[ks][js][i];
00095     W[i] = Cons_to_Prim(&(pG->U[ks][js][i]));
00096   }
00097 
00098 /*=== STEP 1: Compute first-order fluxes at t^{n} in x1-direction ============*/
00099 /* No source terms are needed since there is no temporal evolution */
00100 
00101 /*--- Step 1a ------------------------------------------------------------------
00102  * Load 1D vector of primitive variables;  
00103  * W1d = (d, V1, V2, V3, P, B2c, B3c, s[n])
00104  */
00105 
00106   for (i=is-nghost; i<=ie+nghost; i++) {
00107     W1d[i].d  = W[i].d;
00108     W1d[i].Vx = W[i].V1;
00109     W1d[i].Vy = W[i].V2;
00110     W1d[i].Vz = W[i].V3;
00111     W1d[i].P  = W[i].P;
00112 #ifdef MHD
00113     W1d[i].By = W[i].B2c;
00114     W1d[i].Bz = W[i].B3c;
00115     Bxc[i] = W[i].B1c;
00116     Bxi[i] = pG->B1i[ks][js][i];
00117 #endif /* MHD */
00118 #if (NSCALARS > 0)
00119     for (n=0; n<NSCALARS; n++) W1d[i].s[n] = W[i].s[n];
00120 #endif
00121   }
00122 
00123 /*--- Step 1b ------------------------------------------------------------------
00124  * Compute first-order L/R states */
00125 
00126 /* Ensure that W & U are consistent */
00127   for (i=is-nghost; i<=ie+nghost; i++) {
00128     U1d[i] = Prim1D_to_Cons1D(&W1d[i],&Bxc[i]);
00129   }
00130 
00131   for (i=il; i<=ie+nghost; i++) {
00132     Wl[i] = W1d[i-1];
00133     Wr[i] = W1d[i  ];
00134 
00135     Ul[i] = U1d[i-1];
00136     Ur[i] = U1d[i  ];
00137   }
00138 
00139 /*--- Step 1c ------------------------------------------------------------------
00140  * No source terms needed */
00141 
00142 /*--- Step 1d ------------------------------------------------------------------
00143  * Compute flux in x1-direction */
00144 
00145   for (i=il; i<=ie+nghost; i++) {
00146     fluxes(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1Flux[i]);
00147   }
00148 
00149 /*=== STEPS 2-4: Not needed in 1D ===*/
00150 
00151 /*=== STEP 5: Update cell-centered variables to half-timestep ================*/
00152 
00153 /*--- Step 5a ------------------------------------------------------------------
00154  * Update cell-centered variables (including B2c and B3c) to half-timestep
00155  */
00156 
00157   for (i=il; i<=iu; i++) {
00158     Uhalf[i].d   -= hdtodx1*(x1Flux[i+1].d  - x1Flux[i].d );
00159     Uhalf[i].M1  -= hdtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00160     Uhalf[i].M2  -= hdtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00161     Uhalf[i].M3  -= hdtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00162     Uhalf[i].E   -= hdtodx1*(x1Flux[i+1].E  - x1Flux[i].E );
00163 #ifdef MHD
00164     Uhalf[i].B2c -= hdtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00165     Uhalf[i].B3c -= hdtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00166 #endif /* MHD */
00167 #if (NSCALARS > 0)
00168     for (n=0; n<NSCALARS; n++)
00169       Uhalf[i].s[n] -= hdtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00170 #endif
00171 #ifdef FIRST_ORDER_FLUX_CORRECTION
00172           x1FluxP[i] = x1Flux[i];
00173 #endif
00174   }
00175 
00176 /*=== STEP 6: Add source terms to predict values at half-timestep ============*/
00177 
00178 /*--- Step 6a ------------------------------------------------------------------
00179  * Add source terms from a static gravitational potential for 0.5*dt to predict
00180  * step.  To improve conservation of total energy, we average the energy
00181  * source term computed at cell faces.
00182  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
00183  */
00184 
00185   if (StaticGravPot != NULL){
00186     for (i=il; i<=iu; i++) {
00187       cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00188       phic = (*StaticGravPot)((x1            ),x2,x3);
00189       phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00190       phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00191 
00192       Uhalf[i].M1 -= hdtodx1*pG->U[ks][js][i].d*(phir-phil);
00193       Uhalf[i].E -= hdtodx1*(x1Flux[i  ].d*(phic - phil) +
00194                              x1Flux[i+1].d*(phir - phic));
00195     }
00196   }
00197 
00198 /*=== STEP 7: Conserved->Primitive variable inversion at t^{n+1/2} ===========*/
00199         
00200 /* Invert conserved variables at t^{n+1/2} to primitive variables. With FOFC, 
00201  * if cell-centered d < 0, P< 0, or v^2 > 1, correct by switching back to 
00202  * values at beginning of step, rendering update first order in time for that
00203  * cell.
00204  */
00205 
00206 #ifdef FIRST_ORDER_FLUX_CORRECTION        
00207   negd = 0;
00208   negP = 0;
00209   superl = 0;
00210   flag_cell = 0;
00211 #endif
00212   for (i=il; i<=iu; i++) {
00213     Whalf[i] = Cons_to_Prim(&Uhalf[i]);
00214 #ifdef FIRST_ORDER_FLUX_CORRECTION   
00215     if (Whalf[i].d < 0.0) {
00216       flag_cell = 1;
00217       negd++;
00218     }
00219     if (Whalf[i].P < 0.0) {
00220       flag_cell = 1;
00221       negP++;
00222     }
00223     Vsq = SQR(Whalf[i].V1) + SQR(Whalf[i].V2) + SQR(Whalf[i].V3);
00224     if (Vsq > 1.0) {
00225       flag_cell = 1;
00226       superl++;
00227     }
00228     if (flag_cell != 0) {
00229       Whalf[i].d = W[i].d;
00230       Whalf[i].V1 = W[i].V1;
00231       Whalf[i].V2 = W[i].V2;
00232       Whalf[i].V3 = W[i].V3;
00233       Whalf[i].P = W[i].P;
00234       flag_cell=0;
00235     }
00236 #endif
00237   }
00238 
00239 #ifdef FIRST_ORDER_FLUX_CORRECTION  
00240   if (negd > 0 || negP > 0 || superl > 0)
00241     printf("[Step7]: %i cells had d<0; %i cells had P<0; %i cells had v>1\n"
00242                                  ,negd,negP,superl);
00243 #endif
00244 
00245 
00246 /*=== STEP 8: Compute second-order L/R x1-interface states ===================*/
00247 
00248 /*--- Step 8a ------------------------------------------------------------------
00249  * Load 1D vector of primitive variables;
00250  * W = (d, V1, V2, V3, P, B2c, B3c, s[n])
00251  */
00252 
00253   for (i=il; i<=iu; i++) {
00254     W1d[i].d  = Whalf[i].d;
00255     W1d[i].Vx = Whalf[i].V1;
00256     W1d[i].Vy = Whalf[i].V2;
00257     W1d[i].Vz = Whalf[i].V3;
00258     W1d[i].P  = Whalf[i].P;
00259 #ifdef MHD
00260     W1d[i].By = Whalf[i].B2c;
00261     W1d[i].Bz = Whalf[i].B3c;
00262     Bxc[i]  = Whalf[i].B1c;
00263 #endif /* MHD */
00264 #if (NSCALARS > 0)
00265     for (n=0; n<NSCALARS; n++) W1d[i].s[n] = Whalf[i].s[n];
00266 #endif /* NSCALARS */
00267   }
00268 
00269 /*--- Step 8b ------------------------------------------------------------------
00270  * Compute L/R states on x1-interfaces, store into arrays
00271  */
00272 
00273   lr_states(pG,W1d,Bxc,pG->dt,pG->dx1,is,ie,Wl,Wr,cart_x1);
00274 
00275 #ifdef FIRST_ORDER_FLUX_CORRECTION
00276   for (i=il; i<=iu; i++) {
00277     Vsq = SQR(Wl[i].Vx) + SQR(Wl[i].Vy) + SQR(Wl[i].Vz);
00278     if (Vsq > 1.0){
00279       Wl[i] = W1d[i];
00280       Wr[i] = W1d[i];
00281     }
00282     Vsq = SQR(Wr[i].Vx) + SQR(Wr[i].Vy) + SQR(Wr[i].Vz);
00283     if (Vsq > 1.0){
00284       Wl[i] = W1d[i];
00285       Wr[i] = W1d[i];
00286     }
00287   }
00288 #endif
00289 
00290   for (i=is; i<=ie+1; i++) {
00291     Wl_x1Face[i] = Wl[i];
00292     Wr_x1Face[i] = Wr[i];
00293   }
00294 
00295 /*=== STEPS 9-10: Not needed in 1D ===*/
00296 
00297 /*=== STEP 11: Compute x1-Flux ===============================================*/
00298 
00299 /*--- Step 11b -----------------------------------------------------------------
00300  * Compute second-order fluxes in x1-direction
00301  */
00302 
00303   for (i=is; i<=ie+1; i++) {
00304     Ul[i] = Prim1D_to_Cons1D(&Wl_x1Face[i],&Bxi[i]);
00305     Ur[i] = Prim1D_to_Cons1D(&Wr_x1Face[i],&Bxi[i]);
00306 
00307     fluxes(Ul[i],Ur[i],Wl_x1Face[i],Wr_x1Face[i],Bxi[i],&x1Flux[i]);
00308   }
00309 
00310 /*=== STEP 12: Not needed in 1D ===*/
00311         
00312 /*=== STEP 13: Add source terms for a full timestep using n+1/2 states =======*/
00313        
00314 /*--- Step 13a -----------------------------------------------------------------
00315  * Add gravitational source terms due to a Static Potential
00316  * To improve conservation of total energy, we average the energy
00317  * source term computed at cell faces.
00318  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
00319  */
00320 
00321   if (StaticGravPot != NULL){
00322     for (i=is; i<=ie; i++) {
00323       cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00324       phic = (*StaticGravPot)((x1            ),x2,x3);
00325       phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00326       phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00327 
00328       pG->U[ks][js][i].M1 -= dtodx1*Uhalf[i].d*(phir-phil);
00329 #ifndef BAROTROPIC
00330       pG->U[ks][js][i].E -= dtodx1*(x1Flux[i  ].d*(phic - phil) +
00331                                     x1Flux[i+1].d*(phir - phic));
00332 #endif
00333     }
00334   }
00335 
00336 /*=== STEP 14: Update cell-centered values for a full timestep ===============*/
00337 
00338 /*--- Step 14a -----------------------------------------------------------------
00339  * Update cell-centered variables in pG (including B2c and B3c) using x1-Fluxes
00340  */
00341 
00342   for (i=is; i<=ie; i++) {
00343     pG->U[ks][js][i].d  -= dtodx1*(x1Flux[i+1].d  - x1Flux[i].d );
00344     pG->U[ks][js][i].M1 -= dtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00345     pG->U[ks][js][i].M2 -= dtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00346     pG->U[ks][js][i].M3 -= dtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00347 #ifndef BAROTROPIC
00348     pG->U[ks][js][i].E  -= dtodx1*(x1Flux[i+1].E  - x1Flux[i].E );
00349 #endif /* BAROTROPIC */
00350 #ifdef MHD
00351     pG->U[ks][js][i].B2c -= dtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00352     pG->U[ks][js][i].B3c -= dtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00353 /* For consistency, set B2i and B3i to cell-centered values.  */
00354     pG->B2i[ks][js][i] = pG->U[ks][js][i].B2c;
00355     pG->B3i[ks][js][i] = pG->U[ks][js][i].B3c;
00356 #endif /* MHD */
00357 #if (NSCALARS > 0)
00358     for (n=0; n<NSCALARS; n++)
00359       pG->U[ks][js][i].s[n] -= dtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00360 #endif
00361   }
00362 
00363 #ifdef FIRST_ORDER_FLUX_CORRECTION
00364 /*=== STEP 15: First-order flux correction ===================================*/
00365 
00366 /*--- Step 15a -----------------------------------------------------------------
00367  * If cell-centered d or P have gone negative, or if v^2 > 1, correct
00368  * by using 1st order predictor fluxes */
00369         
00370   for (i=is; i<=ie; i++) {
00371       Wcheck = check_Prim(&(pG->U[ks][js][i]));
00372       if (Wcheck.d < 0.0) {
00373         flag_cell = 1;
00374         BadCell.i = i;
00375         BadCell.j = js;
00376         BadCell.k = ks;
00377         negd++;
00378       }
00379       if (Wcheck.P < 0.0) {
00380         flag_cell = 1;
00381         BadCell.i = i;
00382         BadCell.j = js;
00383         BadCell.k = ks;
00384         negP++;
00385       }
00386       Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00387       if (Vsq > 1.0) {
00388         flag_cell = 1;
00389         BadCell.i = i;
00390         BadCell.j = js;
00391         BadCell.k = ks;
00392         superl++;
00393       }
00394       if (flag_cell != 0) {
00395         FixCell(pG, BadCell);
00396         flag_cell=0;
00397       }
00398 
00399   }
00400 
00401   if (negd > 0 || negP > 0 || superl > 0){
00402     printf("[Step15a]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00403     printf("[Step15a]: %i cells had v>1 at 1st correction\n",superl);
00404   }
00405 
00406 /*--- Step 15b -----------------------------------------------------------------
00407  * In SR the first-order flux correction can fail to fix an unphysical state.
00408  * We must fix these cells in order to avoid NaN's at the next timestep,
00409  * particuarly if v^2 > 1. We have 2 approaches; firstly, we use the entropy
00410  * equation (which we have not applied a 1st order flux correction to) to
00411  * calculate the pressure and the Lorentz factor of the gas. If this produces
00412  * and unphysical state, then we floor the pressure and iterate on v^2 until
00413  * v^2 < 1. Possibly could improved by averaging density and pressure from
00414  * adjacent cells and then calculating pressure.
00415  */
00416 
00417 #ifdef MHD
00418   fail = 0;
00419   negd = 0;
00420   negP = 0;
00421   final = 0;
00422   superl = 0;
00423   for (i=is; i<=ie; i++) {
00424     flag_cell=0;
00425     Wcheck = check_Prim(&(pG->U[ks][js][i]));
00426     if (Wcheck.d < 0.0) {
00427       flag_cell = 1;
00428       negd++;
00429     }
00430     if (Wcheck.P < 0.0) {
00431       flag_cell = 1;
00432       negP++;
00433     }
00434     Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00435     if (Vsq > 1.0) {
00436       flag_cell = 1;
00437       superl++;
00438     }
00439     if (flag_cell != 0) {
00440       final++;
00441       Wcheck = fix_vsq (&(pG->U[ks][js][i]));
00442       U = Prim_to_Cons(&Wcheck);
00443       pG->U[ks][js][i].d = U.d;
00444       pG->U[ks][js][i].M1 = U.M1;
00445       pG->U[ks][js][i].M2 = U.M2;
00446       pG->U[ks][js][i].M3 = U.M3;
00447       pG->U[ks][js][i].E = U.E;
00448       Wcheck = check_Prim(&(pG->U[ks][js][i]));
00449       Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00450       if (Wcheck.d < 0.0 || Wcheck.P < 0.0 || Vsq > 1.0){
00451         fail++;
00452       }
00453     }
00454   }
00455 
00456   if (negd > 0 || negP > 0 || superl > 0) {
00457     printf("[Step15b]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00458     printf("[Step15b]: %i cells had v>1 after 1st order correction\n",superl);
00459     printf("[Step15b]: %i cells required a  final fix\n",final);
00460     printf("[Step15b]: %i cells had an unphysical state\n",fail);
00461   }
00462 #endif /* MHD */
00463 #endif /* FIRST_ORDER_FLUX_CORRECTION */
00464 
00465 
00466 
00467 #ifdef STATIC_MESH_REFINEMENT
00468 /*--- Step 13d -----------------------------------------------------------------
00469  * With SMR, store fluxes at boundaries of child and parent grids.
00470  */
00471 
00472 /* x1-boundaries of child Grids (interior to THIS Grid) */
00473   for (ncg=0; ncg<pG->NCGrid; ncg++) {
00474     for (dim=0; dim<2; dim++){
00475       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
00476 
00477         if (dim==0) i = pG->CGrid[ncg].ijks[0];
00478         if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
00479 
00480         pG->CGrid[ncg].myFlx[dim][ks][js].d  = x1Flux[i].d;
00481         pG->CGrid[ncg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00482         pG->CGrid[ncg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00483         pG->CGrid[ncg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00484 #ifndef BAROTROPIC
00485         pG->CGrid[ncg].myFlx[dim][ks][js].E  = x1Flux[i].E;
00486 #endif /* BAROTROPIC */
00487 #ifdef MHD
00488         pG->CGrid[ncg].myFlx[dim][ks][js].B1c = 0.0;
00489         pG->CGrid[ncg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00490         pG->CGrid[ncg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00491 #endif /* MHD */
00492 #if (NSCALARS > 0)
00493         for (n=0; n<NSCALARS; n++)
00494           pG->CGrid[ncg].myFlx[dim][ks][js].s[n]  = x1Flux[i].s[n];
00495 #endif
00496       }
00497     }
00498   }
00499 
00500 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
00501   for (npg=0; npg<pG->NPGrid; npg++) {
00502     for (dim=0; dim<2; dim++){
00503       if (pG->PGrid[npg].myFlx[dim] != NULL) {
00504 
00505         if (dim==0) i = pG->PGrid[npg].ijks[0];
00506         if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
00507 
00508         pG->PGrid[npg].myFlx[dim][ks][js].d  = x1Flux[i].d;
00509         pG->PGrid[npg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00510         pG->PGrid[npg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00511         pG->PGrid[npg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00512 #ifndef BAROTROPIC
00513         pG->PGrid[npg].myFlx[dim][ks][js].E  = x1Flux[i].E;
00514 #endif /* BAROTROPIC */
00515 #ifdef MHD
00516         pG->PGrid[npg].myFlx[dim][ks][js].B1c = 0.0;
00517         pG->PGrid[npg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00518         pG->PGrid[npg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00519 #endif /* MHD */
00520 #if (NSCALARS > 0)
00521         for (n=0; n<NSCALARS; n++)
00522           pG->PGrid[npg].myFlx[dim][ks][js].s[n]  = x1Flux[i].s[n];
00523 #endif
00524       }
00525     }
00526   }
00527 
00528 #endif /* STATIC_MESH_REFINEMENT */
00529 
00530   return;
00531 }
00532 
00533 /*----------------------------------------------------------------------------*/
00534 /*! \fn void integrate_init_1d(MeshS *pM)
00535  *  \brief Allocate temporary integration arrays */
00536 void integrate_init_1d(MeshS *pM)
00537 {
00538   int size1=0,nl,nd;
00539 
00540 /* Cycle over all Grids on this processor to find maximum Nx1 */
00541   for (nl=0; nl<(pM->NLevels); nl++){
00542     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00543       if (pM->Domain[nl][nd].Grid != NULL) {
00544         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00545           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00546         }
00547       }
00548     }
00549   }
00550 
00551   size1 = size1 + 2*nghost;
00552 
00553   if ((Wl_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00554   if ((Wr_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00555   if ((x1Flux   =(Cons1DS*)malloc(size1*sizeof(Cons1DS))) ==NULL) goto on_error;
00556 #ifdef FIRST_ORDER_FLUX_CORRECTION
00557   if ((x1FluxP  =(Cons1DS*)malloc(size1*sizeof(Cons1DS))) ==NULL) goto on_error;
00558 #endif /* FIRST_ORDER_FLUX_CORRECTION */
00559 
00560   if ((Bxc = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00561   if ((Bxi = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00562 
00563   if ((U1d= (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00564   if ((Ul = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00565   if ((Ur = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00566   if ((W1d= (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00567   if ((Wl = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00568   if ((Wr = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00569 
00570   if ((Uhalf = (ConsS*)malloc(size1*sizeof(ConsS)))==NULL) goto on_error;
00571   if ((Whalf = (PrimS*)malloc(size1*sizeof(PrimS)))==NULL) goto on_error;
00572   if ((W     = (PrimS*)malloc(size1*sizeof(PrimS)))==NULL) goto on_error;
00573 
00574   return;
00575 
00576   on_error:
00577     integrate_destruct();
00578     ath_error("[integrate_init]: malloc returned a NULL pointer\n");
00579 }
00580 
00581 /*----------------------------------------------------------------------------*/
00582 /*! \fn void integrate_destruct_1d(void)
00583  *  \brief Free temporary integration arrays */
00584 void integrate_destruct_1d(void)
00585 {
00586   if (Wl_x1Face != NULL) free(Wl_x1Face);
00587   if (Wr_x1Face != NULL) free(Wr_x1Face);
00588   if (x1Flux    != NULL) free(x1Flux);
00589 #ifdef FIRST_ORDER_FLUX_CORRECTION
00590   if (x1FluxP   != NULL) free(x1FluxP);
00591 #endif
00592 
00593 
00594   if (Bxc != NULL) free(Bxc);
00595   if (Bxi != NULL) free(Bxi);
00596 
00597   if (U1d      != NULL) free(U1d);
00598   if (Ul       != NULL) free(Ul);
00599   if (Ur       != NULL) free(Ur);
00600   if (W1d      != NULL) free(W1d);
00601   if (Wl       != NULL) free(Wl);
00602   if (Wr       != NULL) free(Wr);
00603 
00604   if (Uhalf    != NULL) free(Uhalf);
00605   if (Whalf    != NULL) free(Whalf);
00606   if (W        != NULL) free(W);
00607 
00608   return;
00609 }
00610 
00611 #ifdef FIRST_ORDER_FLUX_CORRECTION
00612 /*----------------------------------------------------------------------------*/
00613 /*! \fn static void FixCell(GridS *pG, Int3Vect indx)
00614  *  \brief Uses first order fluxes to fix negative d,P or superluminal v
00615  */ 
00616 static void FixCell(GridS *pG, Int3Vect indx)
00617 {
00618   int ks=pG->ks,js=pG->js;
00619   Real dtodx1=pG->dt/pG->dx1;
00620   Cons1DS x1FD_i, x1FD_ip1, x2FD_j, x2FD_jp1;
00621   
00622   /* Compute difference of predictor and corrector fluxes at cell faces */
00623   
00624   x1FD_i.d    = x1Flux[indx.i  ].d - x1FluxP[indx.i  ].d;
00625   x1FD_ip1.d  = x1Flux[indx.i+1].d - x1FluxP[indx.i+1].d;
00626         
00627   x1FD_i.Mx    = x1Flux[indx.i  ].Mx - x1FluxP[indx.i  ].Mx;
00628   x1FD_ip1.Mx  = x1Flux[indx.i+1].Mx - x1FluxP[indx.i+1].Mx;
00629         
00630   x1FD_i.My    = x1Flux[indx.i  ].My - x1FluxP[indx.i  ].My;
00631   x1FD_ip1.My  = x1Flux[indx.i+1].My - x1FluxP[indx.i+1].My;
00632         
00633   x1FD_i.Mz    = x1Flux[indx.i  ].Mz - x1FluxP[indx.i  ].Mz;
00634   x1FD_ip1.Mz  = x1Flux[indx.i+1].Mz - x1FluxP[indx.i+1].Mz;
00635         
00636 #ifdef MHD
00637   x1FD_i.By    = x1Flux[indx.i  ].By - x1FluxP[indx.i  ].By;
00638   x1FD_ip1.By  = x1Flux[indx.i+1].By - x1FluxP[indx.i+1].By;
00639         
00640   x1FD_i.Bz    = x1Flux[indx.i  ].Bz - x1FluxP[indx.i  ].Bz;
00641   x1FD_ip1.Bz  = x1Flux[indx.i+1].Bz - x1FluxP[indx.i+1].Bz;
00642 #endif
00643         
00644 #ifndef BAROTROPIC
00645   x1FD_i.E    = x1Flux[indx.i  ].E - x1FluxP[indx.i  ].E;
00646   x1FD_ip1.E  = x1Flux[indx.i+1].E - x1FluxP[indx.i+1].E;
00647 #endif /* BAROTROPIC */
00648         
00649   /* Use flux differences to correct bad cell */
00650   pG->U[ks][js][indx.i].d  += dtodx1*(x1FD_ip1.d  - x1FD_i.d );
00651   pG->U[ks][js][indx.i].M1 += dtodx1*(x1FD_ip1.Mx - x1FD_i.Mx);
00652   pG->U[ks][js][indx.i].M2 += dtodx1*(x1FD_ip1.My - x1FD_i.My);
00653   pG->U[ks][js][indx.i].M3 += dtodx1*(x1FD_ip1.Mz - x1FD_i.Mz);
00654 #ifdef MHD
00655   pG->U[ks][js][indx.i].B2c += dtodx1*(x1FD_ip1.By - x1FD_i.By);
00656   pG->U[ks][js][indx.i].B3c += dtodx1*(x1FD_ip1.Bz - x1FD_i.Bz);
00657   /* For consistency, set B2i and B3i to cell-centered values.  */
00658   pG->B2i[ks][js][indx.i] = pG->U[ks][js][indx.i].B2c;
00659   pG->B3i[ks][js][indx.i] = pG->U[ks][js][indx.i].B3c;
00660 #endif
00661 #ifndef BAROTROPIC
00662   pG->U[ks][js][indx.i].E  += dtodx1*(x1FD_ip1.E  - x1FD_i.E );
00663 #endif /* BAROTROPIC */
00664         
00665         
00666   /* Use flux differences to correct bad cell neighbors at i-1 and i+1 */      
00667   if (indx.i > pG->is) {
00668     pG->U[ks][js][indx.i-1].d  += dtodx1*(x1FD_i.d );
00669     pG->U[ks][js][indx.i-1].M1 += dtodx1*(x1FD_i.Mx);
00670     pG->U[ks][js][indx.i-1].M2 += dtodx1*(x1FD_i.My);
00671     pG->U[ks][js][indx.i-1].M3 += dtodx1*(x1FD_i.Mz);
00672 #ifdef MHD
00673     pG->U[ks][js][indx.i-1].B2c += dtodx1*(x1FD_i.By);
00674     pG->U[ks][js][indx.i-1].B3c += dtodx1*(x1FD_i.Bz);
00675     /* For consistency, set B2i and B3i to cell-centered values.  */
00676     pG->B2i[ks][js][indx.i-1] = pG->U[ks][js][indx.i-1].B2c;
00677     pG->B3i[ks][js][indx.i-1] = pG->U[ks][js][indx.i-1].B3c;
00678 #endif
00679 #ifndef BAROTROPIC
00680     pG->U[ks][js][indx.i-1].E  += dtodx1*(x1FD_i.E );
00681 #endif /* BAROTROPIC */
00682   }
00683         
00684   if (indx.i < pG->ie) {
00685     pG->U[ks][js][indx.i+1].d  -= dtodx1*(x1FD_ip1.d );
00686     pG->U[ks][js][indx.i+1].M1 -= dtodx1*(x1FD_ip1.Mx);
00687     pG->U[ks][js][indx.i+1].M2 -= dtodx1*(x1FD_ip1.My);
00688     pG->U[ks][js][indx.i+1].M3 -= dtodx1*(x1FD_ip1.Mz);
00689 #ifdef MHD
00690     pG->U[ks][js][indx.i+1].B2c -= dtodx1*(x1FD_ip1.By);
00691     pG->U[ks][js][indx.i+1].B3c -= dtodx1*(x1FD_ip1.Bz);
00692     /* For consistency, set B2i and B3i to cell-centered values.  */
00693     pG->B2i[ks][js][indx.i+1] = pG->U[ks][js][indx.i+1].B2c;
00694     pG->B3i[ks][js][indx.i+1] = pG->U[ks][js][indx.i+1].B3c;
00695 #endif MHD
00696 #ifndef BAROTROPIC
00697     pG->U[ks][js][indx.i+1].E  -= dtodx1*(x1FD_ip1.E );
00698 #endif /* BAROTROPIC */
00699   }
00700         
00701 }
00702 #endif /* FIRST_ORDER_FLUX_CORRECTION */
00703 
00704 #endif /* VL_INTEGRATOR */
00705 
00706 #endif /* SPECIAL_RELATIVITY */
integrators/integrate_1d_vl_sr.c
