Athena: reconstruction/lr_states

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file lr_states_plm.c
00004  *  \brief Second order (piecewise linear) spatial reconstruction using
00005  *   characteristic interpolation in the primitive variables.  
00006  *
00007  * PURPOSE: Second order (piecewise linear) spatial reconstruction using
00008  *   characteristic interpolation in the primitive variables.  With the CTU
00009  *   integrator, a time-evolution (characteristic tracing) step is used to
00010  *   interpolate interface values to the half time level {n+1/2}. 
00011  *
00012  * NOTATION: 
00013  * - W_{L,i-1/2} is reconstructed value on the left-side of interface at i-1/2
00014  * - W_{R,i-1/2} is reconstructed value on the right-side of interface at i-1/2
00015  *
00016  *   The L- and R-states at the left-interface in each cell are indexed i.
00017  * - W_{L,i-1/2} is denoted by Wl[i  ];   W_{R,i-1/2} is denoted by Wr[i  ]
00018  * - W_{L,i+1/2} is denoted by Wl[i+1];   W_{R,i+1/2} is denoted by Wr[i+1]
00019  *
00020  *   Internally, in this routine, Wlv and Wrv are the reconstructed values on
00021  *   the left-and right-side of cell center.  Thus (see Step 8),
00022  * -   W_{L,i-1/2} = Wrv(i-1);  W_{R,i-1/2} = Wlv(i)
00023  *
00024  * CONTAINS PUBLIC FUNCTIONS:
00025  * - lr_states()          - computes L/R states
00026  * - lr_states_init()     - initializes memory for static global arrays
00027  * - lr_states_destruct() - frees memory for static global arrays
00028  *============================================================================*/
00029 
00030 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include "../defs.h"
00034 #include "../athena.h"
00035 #include "../globals.h"
00036 #include "prototypes.h"
00037 #include "../prototypes.h"
00038 
00039 #ifdef SECOND_ORDER_CHAR
00040 #ifdef SPECIAL_RELATIVITY
00041 #error : PLM reconstruction (order=2) cannot be used for special relativity.
00042 #endif /* SPECIAL_RELATIVITY */
00043 
00044 static Real **pW=NULL;
00045 
00046 /*----------------------------------------------------------------------------*/
00047 /*! \fn void lr_states(const GridS *pG, const Prim1DS W[], const Real Bxc[], 
00048  *               const Real dt, const Real dx, const int il, const int iu, 
00049  *               Prim1DS Wl[], Prim1DS Wr[], const int dir)
00050  *  \brief Computes L/R states
00051  * Input Arguments:
00052  *   W = PRIMITIVE variables at cell centers along 1-D slice
00053  *   Bxc = B in direction of slice at cell center
00054  *   dtodx = dt/dx
00055  *   il,iu = lower and upper indices of zone centers in slice
00056  * W and Bxc must be initialized over [il-2:iu+2]
00057  *
00058  * Output Arguments:
00059  *   Wl,Wr = L/R-states of PRIMITIVE variables at interfaces over [il:iu+1]
00060  */
00061 
00062 void lr_states(const GridS *pG, const Prim1DS W[], const Real Bxc[], 
00063                const Real dt, const Real dx, const int il, const int iu, 
00064                Prim1DS Wl[], Prim1DS Wr[], const int dir)
00065 {
00066   int i,n,m;
00067   Real lim_slope1,lim_slope2,qa,qx;
00068   Real ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00069   Real dWc[NWAVE+NSCALARS],dWl[NWAVE+NSCALARS];
00070   Real dWr[NWAVE+NSCALARS],dWg[NWAVE+NSCALARS];
00071   Real dac[NWAVE+NSCALARS],dal[NWAVE+NSCALARS];
00072   Real dar[NWAVE+NSCALARS],dag[NWAVE+NSCALARS],da[NWAVE+NSCALARS];
00073   Real Wlv[NWAVE+NSCALARS],Wrv[NWAVE+NSCALARS];
00074   Real dW[NWAVE+NSCALARS],dWm[NWAVE+NSCALARS];
00075   Real *pWl, *pWr;
00076   Real qx1,qx2,C;
00077 
00078   /* ADDITIONAL VARIABLES REQUIRED FOR CYLINDRICAL COORDINATES */
00079   Real zl,zr,zc,gamma_curv,opg,omg,beta,betai;
00080   int hllallwave_flag = 0;
00081   const Real dtodx = dt/dx;
00082 #ifdef CYLINDRICAL
00083   const Real *r=pG->r, *ri=pG->ri;
00084 #endif /* CYLINDRICAL */
00085 
00086 /* Zero eigenmatrices, set pointer to primitive variables */
00087   for (n=0; n<NWAVE; n++) {
00088     for (m=0; m<NWAVE; m++) {
00089       rem[n][m] = 0.0;
00090       lem[n][m] = 0.0;
00091     }
00092   }
00093   for (i=il-2; i<=iu+2; i++) pW[i] = (Real*)&(W[i]);
00094 
00095 /*========================== START BIG LOOP OVER i =======================*/
00096   for (i=il-1; i<=iu+1; i++) {
00097 
00098 /*--- Step 1. ------------------------------------------------------------------
00099  * Compute eigensystem in primitive variables.  */
00100 
00101 #ifdef HYDRO
00102 #ifdef ISOTHERMAL
00103     esys_prim_iso_hyd(W[i].d,W[i].Vx,             ev,rem,lem);
00104 #else
00105     esys_prim_adb_hyd(W[i].d,W[i].Vx,Gamma*W[i].P,ev,rem,lem);
00106 #endif /* ISOTHERMAL */
00107 #endif /* HYDRO */
00108 
00109 #ifdef MHD
00110 #ifdef ISOTHERMAL
00111     esys_prim_iso_mhd(
00112       W[i].d,W[i].Vx,             Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00113 #else
00114     esys_prim_adb_mhd(
00115       W[i].d,W[i].Vx,Gamma*W[i].P,Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00116 #endif /* ISOTHERMAL */
00117 #endif /* MHD */
00118 
00119 /*--- Step 2. ------------------------------------------------------------------
00120  * Compute centered, L/R, and van Leer differences of primitive variables
00121  * Note we access contiguous array elements by indexing pointers for speed */
00122 
00123 #ifdef CYLINDRICAL
00124     if (dir==1) {
00125       /* compute cylindrical weighting factors */
00126      zc = 1.0/(1.0 - SQR(dx)/(12.0*r[i+1]*r[i-1]));
00127      zl = 1.0/(1.0 - SQR(dx)/(12.0*r[i  ]*r[i-1]));
00128      zr = 1.0/(1.0 - SQR(dx)/(12.0*r[i+1]*r[i  ]));
00129     }
00130 #endif
00131     for (n=0; n<(NWAVE+NSCALARS); n++) {
00132       dWc[n] = pW[i+1][n] - pW[i-1][n];
00133       dWl[n] = pW[i][n]   - pW[i-1][n];
00134       dWr[n] = pW[i+1][n] - pW[i][n];
00135 #ifdef CYLINDRICAL
00136       if (dir==1) {
00137         dWc[n] *= zc;
00138         dWl[n] *= zl;
00139         dWr[n] *= zr;
00140       }
00141 #endif
00142       if (dWl[n]*dWr[n] > 0.0) {
00143         dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00144       } else {
00145         dWg[n] = 0.0;
00146       }
00147     }
00148 
00149 /*--- Step 3. ------------------------------------------------------------------
00150  * Project differences in primitive variables along characteristics */
00151 
00152     for (n=0; n<NWAVE; n++) {
00153       dac[n] = lem[n][0]*dWc[0];
00154       dal[n] = lem[n][0]*dWl[0];
00155       dar[n] = lem[n][0]*dWr[0];
00156       dag[n] = lem[n][0]*dWg[0];
00157       for (m=1; m<NWAVE; m++) {
00158         dac[n] += lem[n][m]*dWc[m];
00159         dal[n] += lem[n][m]*dWl[m];
00160         dar[n] += lem[n][m]*dWr[m];
00161         dag[n] += lem[n][m]*dWg[m];
00162       }
00163     }
00164 
00165 /* Advected variables are treated differently; for them the right and left
00166  * eigenmatrices are simply the identitiy matrix.
00167  */
00168 #if (NSCALARS > 0)
00169     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00170       dac[n] = dWc[n];
00171       dal[n] = dWl[n];
00172       dar[n] = dWr[n];
00173       dag[n] = dWg[n];
00174     }
00175 #endif
00176 
00177 /*--- Step 4. ------------------------------------------------------------------
00178  * Apply monotonicity constraints to characteristic projections */
00179 
00180     for (n=0; n<(NWAVE+NSCALARS); n++) {
00181       da[n] = 0.0;
00182       if (dal[n]*dar[n] > 0.0) {
00183         lim_slope1 = MIN(    fabs(dal[n]),fabs(dar[n]));
00184         lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00185         da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00186       }
00187     }
00188 
00189 /*--- Step 5. ------------------------------------------------------------------
00190  * Project monotonic slopes in characteristic back to primitive variables  */
00191 
00192     for (n=0; n<NWAVE; n++) {
00193       dWm[n] = da[0]*rem[n][0];
00194       for (m=1; m<NWAVE; m++) {
00195         dWm[n] += da[m]*rem[n][m];
00196       }
00197     }
00198 
00199 #if (NSCALARS > 0)
00200     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00201       dWm[n] = da[n];
00202     }
00203 #endif
00204 
00205 /*--- Step 6. ------------------------------------------------------------------
00206  * Limit velocity difference to sound speed (deleted).  Was added to make
00207  * turbulence runs more robust, but found added noise to Noh shocktube and 
00208  * other tests.  See r995 and earlier for this step */
00209 
00210 /*--- Step 7. ------------------------------------------------------------------
00211  * Compute L/R values, ensure they lie between neighboring cell-centered vals */
00212 
00213     gamma_curv = 0.0;
00214 #ifdef CYLINDRICAL
00215     if (dir==1) gamma_curv = dx/(6.0*r[i]);
00216 #endif
00217     opg = 1.0 + gamma_curv;
00218     omg = 1.0 - gamma_curv;
00219     beta  = omg/opg;
00220     betai = opg/omg;
00221 
00222     for (n=0; n<(NWAVE+NSCALARS); n++) {
00223       Wlv[n] = pW[i][n] - 0.5*dWm[n]*opg;
00224       Wrv[n] = pW[i][n] + 0.5*dWm[n]*omg;
00225     }
00226 
00227     for (n=0; n<(NWAVE+NSCALARS); n++) {
00228       C = Wrv[n] + beta*Wlv[n];
00229       Wlv[n] = MAX(MIN(pW[i][n],pW[i-1][n]),Wlv[n]);
00230       Wlv[n] = MIN(MAX(pW[i][n],pW[i-1][n]),Wlv[n]);
00231       Wrv[n] = C - beta*Wlv[n];
00232 
00233       Wrv[n] = MAX(MIN(pW[i][n],pW[i+1][n]),Wrv[n]);
00234       Wrv[n] = MIN(MAX(pW[i][n],pW[i+1][n]),Wrv[n]);
00235       Wlv[n] = (C - Wrv[n])*betai;
00236     }
00237 
00238     for (n=0; n<(NWAVE+NSCALARS); n++) {
00239       dW[n] = Wrv[n] - Wlv[n];
00240     }
00241 
00242 /*--- Step 8. ------------------------------------------------------------------
00243  * Integrate linear interpolation function over domain of dependence defined by
00244  * max(min) eigenvalue
00245  */
00246 
00247     pWl = (Real *) &(Wl[i+1]);
00248     pWr = (Real *) &(Wr[i]);
00249 
00250 #ifndef CTU_INTEGRATOR
00251 
00252     for (n=0; n<(NWAVE+NSCALARS); n++) {
00253       pWl[n] = Wrv[n];
00254       pWr[n] = Wlv[n];
00255     }
00256 
00257 #elif defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00258 
00259     for (n=0; n<NWAVE; n++) {
00260       pWl[n] = Wrv[n];
00261       pWr[n] = Wlv[n];
00262     }
00263 
00264 #ifdef HLL_ALL_WAVE
00265     hllallwave_flag = 1;
00266 #endif
00267 
00268     for (n=0; n<NWAVE; n++) {
00269       qa = 0.0;
00270       if (hllallwave_flag || ev[n] > 0.) {
00271         qx = 0.5*dtodx*ev[n];
00272 #ifdef CYLINDRICAL
00273         if (dir==1) 
00274           qx *= 1.0 - dx*qx/(3.0*(ri[i+1]-dx*qx));
00275 #endif
00276         for (m=0; m<NWAVE; m++) {
00277           qa += lem[n][m]*qx*dW[m];
00278         }
00279         for (m=0; m<NWAVE; m++) pWl[m] -= qa*rem[m][n];
00280       }
00281 
00282       qa = 0.0;
00283       if (hllallwave_flag || ev[n] < 0.) {
00284         qx = 0.5*dtodx*ev[n];
00285 #ifdef CYLINDRICAL
00286         if (dir==1)
00287           qx *= 1.0 - dx*qx/(3.0*(ri[i]-dx*qx));
00288 #endif
00289         for (m=0; m<NWAVE; m++) {
00290           qa += lem[n][m]*qx*dW[m];
00291         }
00292         for (m=0; m<NWAVE; m++) pWr[m] -= qa*rem[m][n];
00293       }
00294     }
00295 
00296 #else  /* include steps 8-9 only if using CTU integrator (AND NOT HLL) */   
00297     qx = 0.5*MAX(ev[NWAVE-1],0.0)*dtodx;
00298 #ifdef CYLINDRICAL
00299     if (dir==1) 
00300       qx *= 1.0 - dx*qx/(3.0*(ri[i+1]-dx*qx));
00301 #endif
00302     for (n=0; n<(NWAVE+NSCALARS); n++) {
00303       pWl[n] = Wrv[n] - qx*dW[n];
00304     }
00305 
00306     qx = -0.5*MIN(ev[0],0.0)*dtodx;
00307 #ifdef CYLINDRICAL
00308     if (dir==1) 
00309       qx *= 1.0 + dx*qx/(3.0*(ri[i]+dx*qx));
00310 #endif
00311     for (n=0; n<(NWAVE+NSCALARS); n++) {
00312       pWr[n] = Wlv[n] + qx*dW[n];
00313     }
00314 
00315 
00316 /*--- Step 9. ------------------------------------------------------------------
00317  * Then subtract amount of each wave n that does not reach the interface
00318  * during timestep (CW eqn 3.5ff).  For HLL fluxes, must subtract waves that
00319  * move in both directions.
00320  */
00321 
00322     for (n=0; n<NWAVE; n++) {
00323       if (ev[n] >= 0.0) {
00324         qa  = 0.0;
00325         qx1 = 0.5*dtodx*ev[NWAVE-1];
00326         qx2 = 0.5*dtodx*ev[n];
00327 #ifdef CYLINDRICAL
00328         if (dir==1) {
00329           qx1 *= 1.0 - dx*qx1/(3.0*(ri[i+1]-dx*qx1));
00330           qx2 *= 1.0 - dx*qx2/(3.0*(ri[i+1]-dx*qx2));
00331         }
00332 #endif
00333         qx = qx1 - qx2;
00334         for (m=0; m<NWAVE; m++) {
00335           qa += lem[n][m]*qx*dW[m];
00336         }
00337         for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00338       }
00339     }
00340 
00341     for (n=0; n<NWAVE; n++) {
00342       if (ev[n] <= 0.0) {
00343         qa = 0.0;
00344         qx1 = -0.5*dtodx*ev[0];
00345         qx2 = -0.5*dtodx*ev[n];
00346 #ifdef CYLINDRICAL
00347         if (dir==1) {
00348           qx1 *= 1.0 + dx*qx1/(3.0*(ri[i]+dx*qx1));
00349           qx2 *= 1.0 + dx*qx2/(3.0*(ri[i]+dx*qx2));
00350         }
00351 #endif
00352         qx = -qx1 + qx2;
00353         for (m=0; m<NWAVE; m++) {
00354           qa += lem[n][m]*qx*dW[m];
00355         }
00356         for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00357       }
00358     }
00359 
00360 /* Wave subtraction for advected quantities */
00361     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00362       if (W[i].Vx > 0.) {
00363         pWl[n] += 0.5*dtodx*(ev[NWAVE-1]-W[i].Vx)*dW[n];
00364       } else if (W[i].Vx < 0.) {
00365         pWr[n] += 0.5*dtodx*(ev[0]-W[i].Vx)*dW[n];
00366       }
00367     }
00368 
00369 #endif /* CTU_INTEGRATOR */
00370 
00371   } /*===================== END BIG LOOP OVER i ===========================*/
00372 
00373   return;
00374 }
00375 
00376 
00377 /*----------------------------------------------------------------------------*/
00378 /*! \fn void lr_states_init(MeshS *pM)
00379  *  \brief Allocate enough memory for work arrays */
00380 
00381 void lr_states_init(MeshS *pM)
00382 {
00383   int nmax,size1=0,size2=0,size3=0,nl,nd;
00384 
00385 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
00386   for (nl=0; nl<(pM->NLevels); nl++){
00387     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00388       if (pM->Domain[nl][nd].Grid != NULL) {
00389         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00390           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00391         }
00392         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00393           size2 = pM->Domain[nl][nd].Grid->Nx[1];
00394         }
00395         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00396           size3 = pM->Domain[nl][nd].Grid->Nx[2];
00397         }
00398       }
00399     }
00400   }
00401 
00402   size1 = size1 + 2*nghost;
00403   size2 = size2 + 2*nghost;
00404   size3 = size3 + 2*nghost;
00405   nmax = MAX((MAX(size1,size2)),size3);
00406 
00407   if ((pW = (Real**)malloc(nmax*sizeof(Real*))) == NULL) goto on_error;
00408 
00409   return;
00410   on_error:
00411     lr_states_destruct();
00412     ath_error("[lr_states_init]: malloc returned a NULL pointer\n");
00413 }
00414 
00415 /*----------------------------------------------------------------------------*/
00416 /*! \fn void lr_states_destruct(void)
00417  *  \brief Free memory used by work arrays */
00418 
00419 void lr_states_destruct(void)
00420 {
00421   if (pW != NULL) free(pW);
00422   return;
00423 }
00424 
00425 #endif /* SECOND_ORDER_CHAR */
reconstruction/lr_states_plm.c
