Athena: reconstruction/lr_states

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file lr_states_prim3.c
00004  *  \brief Third order (piecewise parabolic) spatial reconstruction in the
00005  *   primitive variables using the extremum-preserving limiters of
00006  *   Colella & Sekora.
00007  *
00008  * PURPOSE: Third order (piecewise parabolic) spatial reconstruction in the
00009  *   primitive variables using the extremum-preserving limiters of
00010  *   Colella & Sekora.  With the CTU integrator, a time-evolution
00011  *   (characteristic tracing) step is used to interpolate interface values
00012  *   to the half time level {n+1/2}.
00013  *
00014  *   Limiting is performed in the primitive (rather than characteristic)
00015  *   variables.  When used with the VL integrator, an eigenvalue decomposition
00016  *   is NOT needed.
00017  *
00018  * NOTATION:
00019  * - W_{L,i-1/2} is reconstructed value on the left-side of interface at i-1/2
00020  * - W_{R,i-1/2} is reconstructed value on the right-side of interface at i-1/2
00021  *
00022  *   The L- and R-states at the left-interface in each cell are indexed i.
00023  * - W_{L,i-1/2} is denoted by Wl[i  ];   W_{R,i-1/2} is denoted by Wr[i  ]
00024  * - W_{L,i+1/2} is denoted by Wl[i+1];   W_{R,i+1/2} is denoted by Wr[i+1]
00025  *
00026  *   Internally, in this routine, Wlv and Wrv are the reconstructed values on
00027  *   the left-and right-side of cell center.  Thus (see Step 19),
00028  * -   W_{L,i-1/2} = Wrv(i-1);  W_{R,i-1/2} = Wlv(i)
00029  *
00030  * REFERENCE:
00031  * - P. Colella & P. Woodward, "The piecewise parabolic method (PPM) for
00032  *     gas-dynamical simulations", JCP, 54, 174 (1984).
00033  * - P. Colella & M. Sekora, "A limiter for PPM that preserves accuracy at
00034  *     smooth extrema", JCP, submitted (2007)
00035  *
00036  * CONTAINS PUBLIC FUNCTIONS:
00037  * - lr_states()          - computes L/R states
00038  * - lr_states_init()     - initializes memory for static global arrays
00039  * - lr_states_destruct() - frees memory for static global arrays             */
00040 /*============================================================================*/
00041 
00042 #include <math.h>
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include "../defs.h"
00046 #include "../athena.h"
00047 #include "../globals.h"
00048 #include "prototypes.h"
00049 #include "../prototypes.h"
00050 
00051 #ifdef THIRD_ORDER_PRIM
00052 
00053 static Real **pW=NULL, **Whalf=NULL;
00054 
00055 /*----------------------------------------------------------------------------*/
00056 /*! \fn void lr_states(const GridS *pG, const Prim1DS W[], const Real Bxc[],
00057  *               const Real dt, const Real dx, const int il, const int iu,
00058  *               Prim1DS Wl[], Prim1DS Wr[], const int dir)
00059  *  \brief Computes L/R states
00060  *
00061  * Input Arguments:
00062  * - W = PRIMITIVE variables at cell centers along 1-D slice
00063  * - Bxc = B in direction of slice at cell center
00064  * - dtodx = dt/dx
00065  * - il,iu = lower and upper indices of zone centers in slice
00066  * W and Bxc must be initialized over [il-3:iu+3]
00067  *
00068  * Output Arguments:
00069  * - Wl,Wr = L/R-states of PRIMITIVE variables at interfaces over [il:iu+1]
00070  */
00071 
00072 void lr_states(const GridS *pG, const Prim1DS W[], const Real Bxc[],
00073                const Real dt, const Real dx, const int il, const int iu,
00074                Prim1DS Wl[], Prim1DS Wr[], const int dir)
00075 {
00076   int i,n,m;
00077   Real lim_slope,qa,qb,qc,qx;
00078   Real ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00079   Real d2Wc[NWAVE+NSCALARS],d2Wl[NWAVE+NSCALARS];
00080   Real d2Wr[NWAVE+NSCALARS],d2W [NWAVE+NSCALARS];
00081   Real d2Wlim[NWAVE+NSCALARS];
00082   Real Wlv[NWAVE+NSCALARS],Wrv[NWAVE+NSCALARS];
00083   Real dW[NWAVE+NSCALARS],W6[NWAVE+NSCALARS];
00084   Real *pWl, *pWr;
00085   Real dtodx = dt/dx;
00086 
00087 /* Set pointer to primitive variables */
00088   for (i=il-3; i<=iu+3; i++) pW[i] = (Real*)&(W[i]);
00089 
00090 #ifdef CTU_INTEGRATOR /* zero eigenmatrices if using CTU integrator */
00091   for (n=0; n<NWAVE; n++) {
00092     for (m=0; m<NWAVE; m++) {
00093       rem[n][m] = 0.0;
00094       lem[n][m] = 0.0;
00095     }
00096   }
00097 #endif /* CTU_INTEGRATOR */
00098 
00099 /*--- Step 1. ------------------------------------------------------------------
00100  * Compute interface states (CS eqns 12-15) over entire 1D pencil.  Using usual
00101  * Athena notation that index i for face-centered quantities denotes L-edge
00102  * (interface i-1/2), then Whalf[i] = W[i-1/2]. */
00103 
00104   for (i=il-1; i<=iu+2; i++) {
00105     for (n=0; n<(NWAVE+NSCALARS); n++) {
00106       Whalf[i][n]=(7.0*(pW[i-1][n]+pW[i][n]) - (pW[i-2][n]+pW[i+1][n]))/12.0;
00107     }
00108     for (n=0; n<(NWAVE+NSCALARS); n++) {
00109       d2Wc[n] = 3.0*(pW[i-1][n] - 2.0*Whalf[i][n] + pW[i][n]);
00110       d2Wl[n] = (pW[i-2][n] - 2.0*pW[i-1][n] + pW[i  ][n]);
00111       d2Wr[n] = (pW[i-1][n] - 2.0*pW[i  ][n] + pW[i+1][n]);
00112       d2Wlim[n] = 0.0;
00113       lim_slope = MIN(fabs(d2Wl[n]),fabs(d2Wr[n]));
00114       if (d2Wc[n] > 0.0 && d2Wl[n] > 0.0 && d2Wr[n] > 0.0) {
00115         d2Wlim[n] = SIGN(d2Wc[n])*MIN(1.25*lim_slope,fabs(d2Wc[n]));
00116       }
00117       if (d2Wc[n] < 0.0 && d2Wl[n] < 0.0 && d2Wr[n] < 0.0) {
00118         d2Wlim[n] = SIGN(d2Wc[n])*MIN(1.25*lim_slope,fabs(d2Wc[n]));
00119       }
00120     }
00121     for (n=0; n<(NWAVE+NSCALARS); n++) {
00122       Whalf[i][n] = 0.5*((pW[i-1][n]+pW[i][n]) - d2Wlim[n]/3.0);
00123     }
00124   }
00125 
00126 /*====================== START BIG LOOP OVER i ============================*/
00127   for (i=il-1; i<=iu+1; i++) {
00128 
00129 /*--- Step 2. ------------------------------------------------------------------
00130  * Compute L/R values
00131  * Wlv = W at left  side of cell-center = W[i-1/2] = a_{j,-} in CS
00132  * Wrv = W at right side of cell-center = W[i+1/2] = a_{j,+} in CS
00133  */
00134 
00135     for (n=0; n<(NWAVE+NSCALARS); n++) {
00136       Wlv[n] = Whalf[i  ][n];
00137       Wrv[n] = Whalf[i+1][n];
00138     }
00139 
00140 /*--- Step 3. ------------------------------------------------------------------
00141  * Construct parabolic interpolant (CS eqn 16-19) */
00142 
00143     for (n=0; n<(NWAVE+NSCALARS); n++) {
00144       qa = (Wrv[n]-pW[i][n])*(pW[i][n]-Wlv[n]);
00145       qb = (pW[i-1][n]-pW[i][n])*(pW[i][n]-pW[i+1][n]);
00146       if (qa <= 0.0 && qb <= 0.0) {
00147         qc = 6.0*(pW[i][n] - 0.5*(Wlv[n]+Wrv[n]));
00148         d2W [n] = -2.0*qc;
00149         d2Wc[n] = (pW[i-1][n] - 2.0*pW[i  ][n] + pW[i+1][n]);
00150         d2Wl[n] = (pW[i-2][n] - 2.0*pW[i-1][n] + pW[i  ][n]);
00151         d2Wr[n] = (pW[i  ][n] - 2.0*pW[i+1][n] + pW[i+2][n]);
00152         d2Wlim[n] = 0.0;
00153         lim_slope = MIN(fabs(d2Wl[n]),fabs(d2Wr[n]));
00154         lim_slope = MIN(fabs(d2Wc[n]),lim_slope);
00155         if (d2Wc[n] > 0.0 && d2Wl[n] > 0.0 && d2Wr[n] > 0.0 && d2W[n] > 0.0) {
00156           d2Wlim[n] = SIGN(d2W[n])*MIN(1.25*lim_slope,fabs(d2W[n]));
00157         }
00158         if (d2Wc[n] < 0.0 && d2Wl[n] < 0.0 && d2Wr[n] < 0.0 && d2W[n] < 0.0) {
00159           d2Wlim[n] = SIGN(d2W[n])*MIN(1.25*lim_slope,fabs(d2W[n]));
00160         }
00161         if (d2W[n] == 0.0) {
00162           Wlv[n] = pW[i][n];
00163           Wrv[n] = pW[i][n];
00164         } else {
00165           Wlv[n] = pW[i][n] + (Wlv[n] - pW[i][n])*d2Wlim[n]/d2W[n];
00166           Wrv[n] = pW[i][n] + (Wrv[n] - pW[i][n])*d2Wlim[n]/d2W[n];
00167         }
00168       }
00169     }
00170 
00171 /*--- Step 4. ------------------------------------------------------------------
00172  * Monotonize again (CS eqn 21-22) */
00173 
00174 /*  First tests showed following does not work.  Not sure why, so leave out.
00175     for (n=0; n<(NWAVE+NSCALARS); n++) {
00176       s = SIGN(pW[i+1][n] - pW[i-1][n]);
00177       alphal = Wlv[n] - pW[i][n];
00178       alphar = Wrv[n] - pW[i][n];
00179       if (fabs(alphar)>=2.0*fabs(alphal)) {
00180         dalpha = pW[i+1][n] - pW[i][n]; 
00181         dI = -0.25*alphar*alphar/(alphar + alphal);
00182         if (s*dI >= s*dalpha) {
00183           Wlv[n] = pW[i][n]-2.0*(dalpha + s*sqrt(dalpha*dalpha-dalpha*alphar));
00184           Wrv[n] = pW[i][n]-2.0*(dalpha + s*sqrt(dalpha*dalpha-dalpha*alphal));
00185         }
00186       }
00187       if (fabs(alphal)>=2.0*fabs(alphar)){
00188         dalpha = pW[i-1][n] - pW[i][n]; 
00189         dI = -0.25*alphal*alphal/(alphar + alphal);
00190         if (s*dI >= s*dalpha) {
00191           Wlv[n] = pW[i][n]-2.0*(dalpha + s*sqrt(dalpha*dalpha-dalpha*alphar));
00192           Wrv[n] = pW[i][n]-2.0*(dalpha + s*sqrt(dalpha*dalpha-dalpha*alphal));
00193         }
00194       }
00195     }
00196 */
00197 
00198 /* Monotonize again (CW eqn 1.10), ensure they lie between neighboring
00199  * cell-centered vals */
00200 
00201     for (n=0; n<(NWAVE+NSCALARS); n++) {
00202       qa = (Wrv[n]-pW[i][n])*(pW[i][n]-Wlv[n]);
00203       qb = Wrv[n]-Wlv[n];
00204       qc = 6.0*(pW[i][n] - 0.5*(Wlv[n]+Wrv[n]));
00205       if (qa <= 0.0) {
00206         Wlv[n] = pW[i][n];
00207         Wrv[n] = pW[i][n];
00208       } else if ((qb*qc) > (qb*qb)) {
00209         Wlv[n] = 3.0*pW[i][n] - 2.0*Wrv[n];
00210       } else if ((qb*qc) < -(qb*qb)) {
00211         Wrv[n] = 3.0*pW[i][n] - 2.0*Wlv[n];
00212       }
00213     }
00214 
00215     for (n=0; n<(NWAVE+NSCALARS); n++) {
00216       Wlv[n] = MAX(MIN(pW[i][n],pW[i-1][n]),Wlv[n]);
00217       Wlv[n] = MIN(MAX(pW[i][n],pW[i-1][n]),Wlv[n]);
00218       Wrv[n] = MAX(MIN(pW[i][n],pW[i+1][n]),Wrv[n]);
00219       Wrv[n] = MIN(MAX(pW[i][n],pW[i+1][n]),Wrv[n]);
00220     }
00221 
00222 /*--- Step 5. ------------------------------------------------------------------
00223  * Set L/R values */
00224 
00225     pWl = (Real *) &(Wl[i+1]);
00226     pWr = (Real *) &(Wr[i]);
00227 
00228     for (n=0; n<(NWAVE+NSCALARS); n++) {
00229       pWl[n] = Wrv[n];
00230       pWr[n] = Wlv[n];
00231     }
00232 
00233 #ifdef CTU_INTEGRATOR /* only include steps below if CTU integrator */
00234 
00235 /*--- Step 6. ------------------------------------------------------------------
00236  * Compute eigensystem in primitive variables.  */
00237 
00238 #ifdef HYDRO
00239 #ifdef ISOTHERMAL
00240     esys_prim_iso_hyd(W[i].d,W[i].Vx,             ev,rem,lem);
00241 #else
00242     esys_prim_adb_hyd(W[i].d,W[i].Vx,Gamma*W[i].P,ev,rem,lem);
00243 #endif /* ISOTHERMAL */
00244 #endif /* HYDRO */
00245 
00246 #ifdef MHD
00247 #ifdef ISOTHERMAL
00248     esys_prim_iso_mhd(
00249       W[i].d,W[i].Vx,             Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00250 #else
00251     esys_prim_adb_mhd(
00252       W[i].d,W[i].Vx,Gamma*W[i].P,Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00253 #endif /* ISOTHERMAL */
00254 #endif /* MHD */
00255 
00256 /*--- Step 7. ------------------------------------------------------------------
00257  * Compute coefficients of interpolation parabolae (CW eqn 1.5) */
00258 
00259     for (n=0; n<(NWAVE+NSCALARS); n++) {
00260       dW[n] = Wrv[n] - Wlv[n];
00261       W6[n] = 6.0*(pW[i][n] - 0.5*(Wlv[n] + Wrv[n]));
00262     }
00263 
00264 /*--- Step 8. ------------------------------------------------------------------
00265  * Integrate linear interpolation function over domain of dependence defined by
00266  * max(min) eigenvalue (CW eqn 1.12)
00267  */
00268 
00269     qx = TWO_3RDS*MAX(ev[NWAVE-1],0.0)*dtodx;
00270     for (n=0; n<(NWAVE+NSCALARS); n++) {
00271       pWl[n] -= 0.75*qx*(dW[n] - (1.0 - qx)*W6[n]);
00272     }
00273 
00274     qx = -TWO_3RDS*MIN(ev[0],0.0)*dtodx;
00275     for (n=0; n<(NWAVE+NSCALARS); n++) {
00276       pWr[n] += 0.75*qx*(dW[n] + (1.0 - qx)*W6[n]);
00277     }
00278 
00279 /*--- Step 9. ------------------------------------------------------------------
00280  * Then subtract amount of each wave m that does not reach the interface
00281  * during timestep (CW eqn 3.5ff).  For HLL fluxes, must subtract waves that
00282  * move in both directions, but only to 2nd order.
00283  */
00284 
00285     for (n=0; n<NWAVE; n++) {
00286       if (ev[n] > 0.) {
00287         qa  = 0.0;
00288         qb = 0.5*dtodx*(ev[NWAVE-1]-ev[n]);
00289         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[NWAVE-1]*ev[NWAVE-1] - ev[n]*ev[n]);
00290         for (m=0; m<NWAVE; m++) {
00291           qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00292         }
00293         for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00294 /* For HLL fluxes, subtract wave moving away from interface as well. */
00295 #if defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00296         qa  = 0.0;
00297         qb = 0.5*dtodx*(ev[0]-ev[n]);
00298         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[0]*ev[0] - ev[n]*ev[n]);
00299         for (m=0; m<NWAVE; m++) {
00300           qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00301         }
00302 
00303         for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00304 #endif /* HLL_FLUX */
00305       }
00306     }
00307 
00308     for (n=0; n<NWAVE; n++) {
00309       if (ev[n] < 0.) {
00310         qa = 0.0;
00311         qb = 0.5*dtodx*(ev[0]-ev[n]);
00312         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[0]*ev[0] - ev[n]*ev[n]);
00313         for (m=0; m<NWAVE; m++) {
00314           qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00315         }
00316         for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00317 /* For HLL fluxes, subtract wave moving away from interface as well. */
00318 #if defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00319         qa = 0.0;
00320         qb = 0.5*dtodx*(ev[NWAVE-1]-ev[n]);
00321         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[NWAVE-1]*ev[NWAVE-1] - ev[n]*ev[n]);
00322         for (m=0; m<NWAVE; m++) {
00323           qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00324         }
00325 
00326         for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00327 #endif /* HLL_FLUX */
00328       }
00329     }
00330 
00331 /* Wave subtraction for advected quantities */
00332     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00333       if (W[i].Vx > 0.) {
00334         qb = 0.5*dtodx*(ev[NWAVE-1]-W[i].Vx);
00335         qc = 0.5*dtodx*dtodx*TWO_3RDS*(SQR(ev[NWAVE-1]) - SQR(W[i].Vx));
00336         pWl[n] += qb*(dW[m]-W6[m]) + qc*W6[m];
00337       } else if (W[i].Vx < 0.) {
00338         qb = 0.5*dtodx*(ev[0]-W[i].Vx);
00339         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[0]*ev[0] - W[i].Vx*W[i].Vx);
00340         pWr[n] += qb*(dW[m]+W6[m]) + qc*W6[m];
00341       }
00342     }
00343 
00344 #endif /* CTU_INTEGRATOR */
00345 
00346   } /*======================= END BIG LOOP OVER i ========================*/
00347 
00348   return;
00349 }
00350 
00351 /*----------------------------------------------------------------------------*/
00352 /*! \fn void lr_states_init(MeshS *pM)
00353  *  \brief Allocate enough memory for work arrays */
00354 
00355 void lr_states_init(MeshS *pM)
00356 {
00357   int nmax,size1=0,size2=0,size3=0,nl,nd;
00358 
00359 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
00360   for (nl=0; nl<(pM->NLevels); nl++){
00361     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00362       if (pM->Domain[nl][nd].Grid != NULL) {
00363         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00364           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00365         }
00366         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00367           size2 = pM->Domain[nl][nd].Grid->Nx[1];
00368         }
00369         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00370           size3 = pM->Domain[nl][nd].Grid->Nx[2];
00371         }
00372       }
00373     }
00374   }
00375 
00376   size1 = size1 + 2*nghost;
00377   size2 = size2 + 2*nghost;
00378   size3 = size3 + 2*nghost;
00379   nmax = MAX((MAX(size1,size2)),size3);
00380 
00381   if ((pW = (Real**)malloc(nmax*sizeof(Real*))) == NULL) goto on_error;
00382 
00383   if ((Whalf = (Real**)calloc_2d_array(nmax, NWAVE, sizeof(Real))) == NULL)
00384     goto on_error;
00385 
00386   return;
00387   on_error:
00388     lr_states_destruct();
00389     ath_error("[lr_states_init]: malloc returned a NULL pointer\n");
00390 }
00391 
00392 /*----------------------------------------------------------------------------*/
00393 /*! \fn void lr_states_destruct(void)
00394  *  \brief Free memory used by work arrays */
00395 
00396 void lr_states_destruct(void)
00397 {
00398   if (pW != NULL) free(pW);
00399   if (Whalf != NULL) free_2d_array(Whalf);
00400   return;
00401 }
00402 
00403 #endif /* THIRD_ORDER_PRIM */
reconstruction/lr_states_prim3.c
