Athena: reconstruction/lr_states

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