Athena: reconstruction/lr_states

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file lr_states_ppm.c
00004  *  \brief Third order (piecewise parabolic) spatial reconstruction using
00005  *   characteristic interpolation in the primitive variables.  
00006  *
00007  * PURPOSE: Third order (piecewise parabolic) 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 18),
00022  * -   W_{L,i-1/2} = Wrv(i-1);  W_{R,i-1/2} = Wlv(i)
00023  *
00024  * REFERENCE:
00025  *   P. Colella & P. Woodward, "The piecewise parabolic method (PPM) for
00026  *   gas-dynamical simulations", JCP, 54, 174 (1984).
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 THIRD_ORDER_CHAR
00044 #ifdef SPECIAL_RELATIVITY
00045 #error : PPM reconstruction (order=3) cannot be used for special relativity.
00046 #endif /* SPECIAL_RELATIVITY */
00047 
00048 static Real **pW=NULL, **dWm=NULL, **Wim1h=NULL;
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  *
00056  * Input Arguments:
00057  * - W = PRIMITIVE variables at cell centers along 1-D slice
00058  * - Bxc = B in direction of slice at cell center
00059  * - dtodx = dt/dx
00060  * - il,iu = lower and upper indices of zone centers in slice
00061  * W and Bxc must be initialized over [il-3:iu+3]
00062  *
00063  * Output Arguments:
00064  * - Wl,Wr = L/R-states of PRIMITIVE variables at interfaces over [il:iu+1]
00065  */
00066 
00067 void lr_states(const GridS* pG, const Prim1DS W[], const Real Bxc[],
00068                const Real dt, const Real dx, const int il, const int iu,
00069                Prim1DS Wl[], Prim1DS Wr[], const int dir)
00070 {
00071   int i,n,m;
00072   Real lim_slope1,lim_slope2,qa,qb,qc,qx;
00073   Real ev    [NWAVE],rem    [NWAVE][NWAVE],lem    [NWAVE][NWAVE];
00074   Real ev_ip1[NWAVE],rem_ip1[NWAVE][NWAVE],lem_ip1[NWAVE][NWAVE];
00075   Real dWc[NWAVE+NSCALARS],dWl[NWAVE+NSCALARS];
00076   Real dWr[NWAVE+NSCALARS],dWg[NWAVE+NSCALARS];
00077   Real dac[NWAVE+NSCALARS],dal[NWAVE+NSCALARS];
00078   Real dar[NWAVE+NSCALARS],dag[NWAVE+NSCALARS],da[NWAVE+NSCALARS];
00079   Real Wlv[NWAVE+NSCALARS],Wrv[NWAVE+NSCALARS];
00080   Real dW[NWAVE+NSCALARS],W6[NWAVE+NSCALARS];
00081   Real *pWl, *pWr;
00082   Real qx1,qx2;
00083 
00084   /* ADDITIONAL VARIABLES REQUIRED FOR CYLINDRICAL COORDINATES */
00085   Real ql,qr,qxx1,qxx2,zc,zr,zl,q1,q2,gamma_curv;
00086   int hllallwave_flag = 0;
00087   const Real dtodx = dt/dx;
00088 #ifdef CYLINDRICAL
00089   const Real *r=pG->r, *ri=pG->ri;
00090 #endif
00091 
00092 /* Zero eigenmatrices, set pointer to primitive variables */
00093   for (n=0; n<NWAVE; n++) {
00094     for (m=0; m<NWAVE; m++) {
00095       rem[n][m] = 0.0;
00096       lem[n][m] = 0.0;
00097       rem_ip1[n][m] = 0.0;
00098       lem_ip1[n][m] = 0.0;
00099     }
00100   }
00101   for (i=il-3; i<=iu+3; i++) pW[i] = (Real*)&(W[i]);
00102 
00103 /*====================== START LOOP OVER il-2:il-1 ==========================*/
00104 
00105 /*--- Step 1. ------------------------------------------------------------------
00106  * Compute eigensystem in primitive variables.  */
00107 
00108   for (i=il-2; i<=il-1; i++) {
00109 
00110 #ifdef HYDRO
00111 #ifdef ISOTHERMAL
00112     esys_prim_iso_hyd(W[i].d,W[i].Vx,             ev,rem,lem);
00113 #else
00114     esys_prim_adb_hyd(W[i].d,W[i].Vx,Gamma*W[i].P,ev,rem,lem);
00115 #endif /* ISOTHERMAL */
00116 #endif /* HYDRO */
00117 
00118 #ifdef MHD
00119 #ifdef ISOTHERMAL
00120     esys_prim_iso_mhd(
00121       W[i].d,W[i].Vx,             Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00122 #else
00123     esys_prim_adb_mhd(
00124       W[i].d,W[i].Vx,Gamma*W[i].P,Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00125 #endif /* ISOTHERMAL */
00126 #endif /* MHD */
00127 
00128 /*--- Step 2. ------------------------------------------------------------------
00129  * Compute centered, L/R, and van Leer differences of primitive variables
00130  * Note we access contiguous array elements by indexing pointers for speed */
00131 
00132     for (n=0; n<(NWAVE+NSCALARS); n++) {
00133 #ifdef CYLINDRICAL
00134       if (dir==1) {
00135         dWc[n] = (pW[i+1][n]*r[i+1] - pW[i-1][n]*r[i-1])/r[i];
00136         dWl[n] = (pW[i  ][n]*r[i  ] - pW[i-1][n]*r[i-1])/ri[i];
00137         dWr[n] = (pW[i+1][n]*r[i+1] - pW[i  ][n]*r[i  ])/ri[i+1];
00138       }
00139       else
00140 #endif
00141       {
00142         dWc[n] = pW[i+1][n] - pW[i-1][n];
00143         dWl[n] = pW[i  ][n] - pW[i-1][n];
00144         dWr[n] = pW[i+1][n] - pW[i  ][n];
00145       }
00146       if (dWl[n]*dWr[n] > 0.0) {
00147         dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00148       } else {
00149         dWg[n] = 0.0;
00150       }
00151     }
00152 
00153 /*--- Step 3. ------------------------------------------------------------------
00154  * Project differences in primitive variables along characteristics */
00155 
00156     for (n=0; n<NWAVE; n++) {
00157       dac[n] = lem[n][0]*dWc[0];
00158       dal[n] = lem[n][0]*dWl[0];
00159       dar[n] = lem[n][0]*dWr[0];
00160       dag[n] = lem[n][0]*dWg[0];
00161       for (m=1; m<NWAVE; m++) {
00162         dac[n] += lem[n][m]*dWc[m];
00163         dal[n] += lem[n][m]*dWl[m];
00164         dar[n] += lem[n][m]*dWr[m];
00165         dag[n] += lem[n][m]*dWg[m];
00166       }
00167     }
00168 
00169 /* Advected variables are treated differently; for them the right and left
00170  * eigenmatrices are simply the identitiy matrix.
00171  */
00172 #if (NSCALARS > 0)
00173     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00174       dac[n] = dWc[n];
00175       dal[n] = dWl[n];
00176       dar[n] = dWr[n];
00177       dag[n] = dWg[n];
00178     }
00179 #endif
00180 
00181 /*--- Step 4. ------------------------------------------------------------------
00182  * Apply monotonicity constraints to characteristic projections */
00183 
00184     for (n=0; n<(NWAVE+NSCALARS); n++) {
00185       da[n] = 0.0;
00186       if (dal[n]*dar[n] > 0.0) {
00187         lim_slope1 = MIN(    fabs(dal[n]),fabs(dar[n]));
00188         lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00189         da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00190       }
00191     }
00192 
00193 /*--- Step 5. ------------------------------------------------------------------
00194  * Project monotonic slopes in characteristic back to primitive variables  */
00195 
00196     for (n=0; n<NWAVE; n++) {
00197       dWm[i][n] = da[0]*rem[n][0];
00198       for (m=1; m<NWAVE; m++) {
00199         dWm[i][n] += da[m]*rem[n][m];
00200       }
00201     }
00202 
00203 #if (NSCALARS > 0)
00204     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00205       dWm[i][n] = da[n];
00206     }
00207 #endif
00208 
00209 /*--- Step 6. ------------------------------------------------------------------
00210  * Limit velocity difference to sound speed
00211  * Limit velocity so momentum is always TVD (using only minmod limiter)
00212  * CURRENTLY NOT USED.  Was added to make code more robust for turbulence
00213  * simulations, but found it added noise to Noh shocktube.
00214  */ 
00215 
00216 #ifdef H_CORRECTION
00217 /* 
00218 #ifdef ISOTHERMAL
00219     qa = Iso_csound;
00220 #else
00221     qa = sqrt(Gamma*W[i].P/W[i].d);
00222 #endif
00223     dWm[i][1] = SIGN(dWm[i][1])*MIN(fabs(dWm[i][1]),qa);
00224 */
00225 #endif /* H_CORRECTION */
00226 /* 
00227     qa = W[i  ].Vx*W[i  ].d - W[i-1].Vx*W[i-1].d;
00228     qb = W[i+1].Vx*W[i+1].d - W[i  ].Vx*W[i  ].d;
00229     qc = W[i+1].Vx*W[i+1].d - W[i-1].Vx*W[i-1].d;
00230     qx = SIGN(qc)*MIN(2.0*MIN(fabs(qa),fabs(qb)), 0.5*fabs(qc));
00231 
00232     if ((-W[i].Vx*dWm[i][0]) > 0.0) {
00233       qa = 0.0;
00234       qb = -W[i].Vx*dWm[i][0];
00235     } else {
00236       qa = -W[i].Vx*dWm[i][0];
00237       qb = 0.0;
00238     }
00239     if (qx > 0.0) {
00240       qb += qx;
00241     } else {
00242       qa += qx;
00243     }
00244     qa = qa/W[i].d;
00245     qb = qb/W[i].d;
00246 
00247     dWm[i][1] = MIN(dWm[i][1],qb);
00248     dWm[i][1] = MAX(dWm[i][1],qa);
00249 */
00250   }
00251 /*====================== END LOOP OVER il-2:il-1 =========================*/
00252 
00253 
00254 /*--- Step 7. ------------------------------------------------------------------
00255  * Construct parabolic interpolant in primitive variables at left-interface
00256  * of cell il-1 ("W[il-3/2]", CW eqn 1.6) using linear TVD slopes at il-2 and
00257  * il-1 computed in Steps 2-7.
00258  */
00259 
00260   for (n=0; n<(NWAVE+NSCALARS); n++) {
00261 #ifdef CYLINDRICAL
00262     if (dir==1) {
00263       Wim1h[il-1][n] = (0.5*(pW[il-1][n]*r[il-1] + pW[il-2][n]*r[il-2])
00264                      - (dWm[il-1][n]*r[il-1] - dWm[il-2][n]*r[il-2])/6.0)/ri[il-1];
00265     }
00266     else
00267 #endif
00268     {
00269       Wim1h[il-1][n] =.5*(pW[il-1][n]+pW[il-2][n])-(dWm[il-1][n]-dWm[il-2][n])/6.;
00270     }
00271   }
00272 
00273 /* Loop over il-2:il-1 in Steps 2-7 was necessary to bootstrap method by
00274  * computing Wim1h[il-1].  Now repeat these steps for rest of grid.  Splitting
00275  * into two loops avoids calculating eigensystems twice per cell.
00276  *
00277  * At the start of the loop, rem and lem still store values at i=il-1 computed
00278  * at the end of Step 2.  For each i, the eigensystem at i+1 is stored in
00279  * rem_ip1 and lem_ip1.  At the end of the loop rem[lem] is then set to
00280  * rem_ip1[lem_ip1] in preparation for the next iteration.
00281  */
00282 
00283 /*========================= START BIG LOOP OVER i =========================*/
00284 /* Steps 8-14 below are identical to steps 1-7 above */
00285   for (i=il-1; i<=iu+1; i++) {
00286 
00287 /*--- Step 8. ------------------------------------------------------------------
00288  * Compute eigensystem in primitive variables.  */
00289 
00290 #ifdef HYDRO
00291 #ifdef ISOTHERMAL
00292     esys_prim_iso_hyd(W[i+1].d,W[i+1].Vx,               ev_ip1,rem_ip1,lem_ip1);
00293 #else
00294     esys_prim_adb_hyd(W[i+1].d,W[i+1].Vx,Gamma*W[i+1].P,ev_ip1,rem_ip1,lem_ip1);
00295 #endif /* ISOTHERMAL */
00296 #endif /* HYDRO */
00297 
00298 #ifdef MHD
00299 #ifdef ISOTHERMAL
00300     esys_prim_iso_mhd(W[i+1].d,W[i+1].Vx,
00301       Bxc[i+1],W[i+1].By,W[i+1].Bz,ev_ip1,rem_ip1,lem_ip1);
00302 #else
00303     esys_prim_adb_mhd(W[i+1].d,W[i+1].Vx,Gamma*W[i+1].P,
00304       Bxc[i+1],W[i+1].By,W[i+1].Bz,ev_ip1,rem_ip1,lem_ip1);
00305 #endif /* ISOTHERMAL */
00306 #endif /* MHD */
00307 
00308 /*--- Step 9. ------------------------------------------------------------------
00309  * Compute centered, L/R, and van Leer differences of primitive variables */
00310 
00311     for (n=0; n<(NWAVE+NSCALARS); n++) {
00312 #ifdef CYLINDRICAL
00313       if (dir==1) {
00314         dWc[n] = (pW[i+2][n]*r[i+2] - pW[i  ][n]*r[i  ])/r[i+1];
00315         dWl[n] = (pW[i+1][n]*r[i+1] - pW[i  ][n]*r[i  ])/ri[i+1];
00316         dWr[n] = (pW[i+2][n]*r[i+2] - pW[i+1][n]*r[i+1])/ri[i+2];
00317       }
00318       else
00319 #endif
00320       {
00321         dWc[n] = pW[i+2][n] - pW[i  ][n];
00322         dWl[n] = pW[i+1][n] - pW[i  ][n];
00323         dWr[n] = pW[i+2][n] - pW[i+1][n];
00324       }
00325       if (dWl[n]*dWr[n] > 0.0) {
00326         dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00327       } else {
00328         dWg[n] = 0.0;
00329       }
00330     }
00331 
00332 /*--- Step 10. -----------------------------------------------------------------
00333  * Project differences in primitive variables along characteristics */
00334 
00335     for (n=0; n<NWAVE; n++) {
00336       dac[n] = lem_ip1[n][0]*dWc[0];
00337       dal[n] = lem_ip1[n][0]*dWl[0];
00338       dar[n] = lem_ip1[n][0]*dWr[0];
00339       dag[n] = lem_ip1[n][0]*dWg[0];
00340       for (m=1; m<NWAVE; m++) {
00341         dac[n] += lem_ip1[n][m]*dWc[m];
00342         dal[n] += lem_ip1[n][m]*dWl[m];
00343         dar[n] += lem_ip1[n][m]*dWr[m];
00344         dag[n] += lem_ip1[n][m]*dWg[m];
00345       }
00346     }
00347 
00348 /* Advected variables are treated differently; for them the right and left
00349  * eigenmatrices are simply the identitiy matrix.
00350  */
00351 #if (NSCALARS > 0)
00352     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00353       dac[n] = dWc[n];
00354       dal[n] = dWl[n];
00355       dar[n] = dWr[n];
00356       dag[n] = dWg[n];
00357     }
00358 #endif
00359 
00360 /*--- Step 11. -----------------------------------------------------------------
00361  * Apply monotonicity constraints to characteristic projections */
00362 
00363     for (n=0; n<(NWAVE+NSCALARS); n++) {
00364       da[n] = 0.0;
00365       if (dal[n]*dar[n] > 0.0) {
00366         lim_slope1 = MIN(    fabs(dal[n]),fabs(dar[n]));
00367         lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00368         da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00369       }
00370     }
00371 
00372 /*--- Step 12. -----------------------------------------------------------------
00373  * Project monotonic slopes in characteristic back to primitive variables  */
00374 
00375     for (n=0; n<NWAVE; n++) {
00376       dWm[i+1][n] = da[0]*rem_ip1[n][0];
00377       for (m=1; m<NWAVE; m++) {
00378         dWm[i+1][n] += da[m]*rem_ip1[n][m];
00379       }
00380     }
00381 
00382 #if (NSCALARS > 0)
00383     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00384       dWm[i+1][n] = da[n];
00385     }
00386 #endif
00387 
00388 /*--- Step 13. -----------------------------------------------------------------
00389  * When H-correction defined, limit velocity difference to sound speed
00390  * Limit velocity so momentum is always TVD (using only minmod limiter)
00391  * CURRENTLY NOT USED.  Was added to make code more robust for turbulence
00392  * simulations, but found it added noise to Noh shocktube.
00393  */
00394 
00395 #ifdef H_CORRECTION
00396 /* 
00397 #ifdef ISOTHERMAL
00398     qa = Iso_csound;
00399 #else
00400     qa = sqrt(Gamma*W[i+1].P/W[i+1].d);
00401 #endif
00402     dWm[i+1][1] = SIGN(dWm[i+1][1])*MIN(fabs(dWm[i+1][1]),qa);
00403 */
00404 #endif /* H_CORRECTION */
00405 /* 
00406     qa = W[i+1].Vx*W[i+1].d - W[i  ].Vx*W[i  ].d;
00407     qb = W[i+2].Vx*W[i+2].d - W[i+1].Vx*W[i+1].d;
00408     qc = W[i+2].Vx*W[i+2].d - W[i  ].Vx*W[i  ].d;
00409     qx = SIGN(qc)*MIN(2.0*MIN(fabs(qa),fabs(qb)), 0.5*fabs(qc));
00410 
00411     if ((-W[i+1].Vx*dWm[i+1][0]) > 0.0) {
00412       qa = 0.0;
00413       qb = -W[i+1].Vx*dWm[i+1][0];
00414     } else {
00415       qa = -W[i+1].Vx*dWm[i+1][0];
00416       qb = 0.0;
00417     }
00418     if (qx > 0.0) {
00419       qb += qx;
00420     } else {
00421       qa += qx;
00422     }
00423     qa = qa/W[i+1].d;
00424     qb = qb/W[i+1].d;
00425 
00426     dWm[i+1][1] = MIN(dWm[i+1][1],qb);
00427     dWm[i+1][1] = MAX(dWm[i+1][1],qa);
00428 */
00429 /*--- Step 14. -----------------------------------------------------------------
00430  * Construct parabolic interpolant in primitive variables at left-interface
00431  * of cell i+1 ("W[i+1/2]", CW eqn 1.6) using linear TVD slopes at i and
00432  * i+1 computed in Steps 2-7.
00433  */
00434 
00435     for (n=0; n<(NWAVE+NSCALARS); n++) {
00436 #ifdef CYLINDRICAL
00437       if (dir==1) {
00438         Wim1h[i+1][n] = (0.5*(pW[i+1][n]*r[i+1] + pW[i][n]*r[i]) 
00439                       - (dWm[i+1][n]*r[i+1] - dWm[i][n]*r[i])/6.0)/ri[i+1];
00440       }
00441       else
00442 #endif
00443       {
00444         Wim1h[i+1][n] = 0.5*(pW[i+1][n]+pW[i][n]) - (dWm[i+1][n]-dWm[i][n])/6.0;
00445       }
00446     }
00447 
00448 /*--- Step 15. -----------------------------------------------------------------
00449  * Compute L/R values */
00450 
00451     for (n=0; n<(NWAVE+NSCALARS); n++) {
00452       Wlv[n] = Wim1h[i  ][n];
00453       Wrv[n] = Wim1h[i+1][n];
00454     }
00455 
00456 /*--- Step 16. -----------------------------------------------------------------
00457  * Monotonize again (CW eqn 1.10), ensure they lie between neighboring
00458  * cell-centered vals */
00459 
00460     gamma_curv = 0.0;
00461 #ifdef CYLINDRICAL
00462     if (dir==1) gamma_curv = dx/(6.0*r[i]);
00463 #endif
00464 
00465     for (n=0; n<(NWAVE+NSCALARS); n++) {
00466       qa = (Wrv[n]-pW[i][n])*(pW[i][n]-Wlv[n]);
00467       qb = Wrv[n]-Wlv[n];
00468       qc = 6.0*(pW[i][n] - 0.5*(Wlv[n]*(1.0-gamma_curv) + Wrv[n]*(1.0+gamma_curv)));
00469       if (qa <= 0.0) {
00470         Wlv[n] = pW[i][n];
00471         Wrv[n] = pW[i][n];
00472       } else if ((qb*qc) > (qb*qb)) {
00473         Wlv[n] = (6.0*pW[i][n] - Wrv[n]*(4.0+3.0*gamma_curv))/(2.0-3.0*gamma_curv);
00474       } else if ((qb*qc) < -(qb*qb)) {
00475         Wrv[n] = (6.0*pW[i][n] - Wlv[n]*(4.0-3.0*gamma_curv))/(2.0+3.0*gamma_curv);
00476       }
00477     }
00478 
00479     for (n=0; n<(NWAVE+NSCALARS); n++) {
00480       Wlv[n] = MAX(MIN(pW[i][n],pW[i-1][n]),Wlv[n]);
00481       Wlv[n] = MIN(MAX(pW[i][n],pW[i-1][n]),Wlv[n]);
00482       Wrv[n] = MAX(MIN(pW[i][n],pW[i+1][n]),Wrv[n]);
00483       Wrv[n] = MIN(MAX(pW[i][n],pW[i+1][n]),Wrv[n]);
00484     }
00485 
00486 /*--- Step 17. -----------------------------------------------------------------
00487  * Compute coefficients of interpolation parabolae (CW eqn 1.5) */
00488 
00489     for (n=0; n<(NWAVE+NSCALARS); n++) {
00490       dW[n] = Wrv[n] - Wlv[n];
00491       W6[n] = 6.0*(pW[i][n] - 0.5*(Wlv[n]*(1.0-gamma_curv) + Wrv[n]*(1.0+gamma_curv)));
00492     }
00493 
00494 /*--- Step 18. -----------------------------------------------------------------
00495  * Integrate linear interpolation function over domain of dependence defined by
00496  * max(min) eigenvalue (CW eqn 1.12)
00497  */
00498 
00499     pWl = (Real *) &(Wl[i+1]);
00500     pWr = (Real *) &(Wr[i]);
00501 
00502 #ifndef CTU_INTEGRATOR 
00503 
00504     for (n=0; n<(NWAVE+NSCALARS); n++) {
00505       pWl[n] = Wrv[n];
00506       pWr[n] = Wlv[n];
00507     }
00508 
00509 
00510 #elif defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00511     for (n=0; n<NWAVE; n++) {
00512       pWl[n] = Wrv[n];
00513       pWr[n] = Wlv[n];
00514     }
00515 
00516 #ifdef HLL_ALL_WAVE
00517     hllallwave_flag = 1;
00518 #endif
00519 
00520     for (n=0; n<NWAVE; n++) {
00521       qa = 0.0;
00522       if (hllallwave_flag || ev[n] > 0.0) {
00523         qx1 = 0.5*dtodx*ev[n];
00524         qb  = qx1;
00525         qc  = FOUR_3RDS*SQR(qx1);
00526 #ifdef CYLINDRICAL
00527         if (dir==1) {
00528           qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00529           qb  -= qxx1;
00530           qc  -= 2.0*qx1*qxx1;
00531         }
00532 #endif
00533         for (m=0; m<NWAVE; m++) {
00534           qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00535         }
00536         for (m=0; m<NWAVE; m++) pWl[m] -= qa*rem[m][n];
00537       }
00538 
00539       qa = 0.0;
00540       if (hllallwave_flag || ev[n] < 0.0) {
00541         qx2 = 0.5*dtodx*ev[n];
00542         qb  = qx2;
00543         qc  = FOUR_3RDS*SQR(qx2);
00544 #ifdef CYLINDRICAL
00545         if (dir==1) {
00546           qxx2 = SQR(qx2)*dx/(3.0*(ri[i]-dx*qx2));
00547           qb  -= qxx2;
00548           qc  -= 2.0*qx2*qxx2;
00549         }
00550 #endif
00551         for (m=0; m<NWAVE; m++) {
00552           qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00553         }
00554         for (m=0; m<NWAVE; m++) pWr[m] -= qa*rem[m][n];
00555       }
00556     }
00557 
00558 
00559 #else /* include steps 18-19 only if using CTU integrator */
00560 
00561     qx1  = 0.5*MAX(ev[NWAVE-1],0.0)*dtodx;
00562     qxx1 = 0.0;
00563 #ifdef CYLINDRICAL
00564     if (dir==1) 
00565       qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00566 #endif
00567     for (n=0; n<(NWAVE+NSCALARS); n++) {
00568       pWl[n] = Wrv[n] - qx1 *(dW[n] - (1.0-FOUR_3RDS*qx1)*W6[n])
00569                       + qxx1*(dW[n] - (1.0-      2.0*qx1)*W6[n]);
00570     }
00571 
00572     qx2  = -0.5*MIN(ev[0],0.0)*dtodx;
00573     qxx2 = 0.0;
00574 #ifdef CYLINDRICAL
00575     if (dir==1) 
00576       qxx2 = SQR(qx2)*dx/(3.0*(ri[i]+dx*qx2));
00577 #endif
00578     for (n=0; n<(NWAVE+NSCALARS); n++) {
00579       pWr[n] = Wlv[n] + qx2 *(dW[n] + (1.0-FOUR_3RDS*qx2)*W6[n])
00580                       + qxx2*(dW[n] + (1.0-      2.0*qx2)*W6[n]);
00581     }
00582 
00583 
00584 /*--- Step 19. -----------------------------------------------------------------
00585  * Then subtract amount of each wave m that does not reach the interface
00586  * during timestep (CW eqn 3.5ff).  For HLL fluxes, must subtract waves that
00587  * move in both directions, but only to 2nd order.
00588  */
00589 
00590     for (n=0; n<NWAVE; n++) {
00591       if (ev[n] >= 0.0) {
00592         qa  = 0.0;
00593         qx1 = 0.5*dtodx*ev[NWAVE-1];
00594         qx2 = 0.5*dtodx*ev[n];
00595         qb  = qx1 - qx2;
00596         qc  = FOUR_3RDS*(SQR(qx1) - SQR(qx2));
00597 #ifdef CYLINDRICAL
00598         if (dir==1) {
00599           qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00600           qxx2 = SQR(qx2)*dx/(3.0*(ri[i+1]-dx*qx2));
00601           qb -= qxx1 - qxx2;
00602           qc -= 2.0*(qx1*qxx1 - qx2*qxx2);
00603         }
00604 #endif
00605         for (m=0; m<NWAVE; m++) {
00606           qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00607         }
00608         for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00609       }
00610     }
00611 
00612     for (n=0; n<NWAVE; n++) {
00613       if (ev[n] <= 0.0) {
00614         qa = 0.0;
00615         qx1 = 0.5*dtodx*ev[0];
00616         qx2 = 0.5*dtodx*ev[n];
00617         qb  = qx1 - qx2;
00618         qc  = FOUR_3RDS*(SQR(qx1) - SQR(qx2));
00619 #ifdef CYLINDRICAL
00620         if (dir==1) {
00621           qxx1 = SQR(qx1)*dx/(3.0*(r[i]-dx*qx1));
00622           qxx2 = SQR(qx2)*dx/(3.0*(r[i]-dx*qx2));
00623           qb -= qxx1 - qxx2;
00624           qc -= 2.0*(qx1*qxx1 - qx2*qxx2);
00625         }
00626 #endif
00627         for (m=0; m<NWAVE; m++) {
00628           qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00629         }
00630         for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00631       }
00632     }
00633 
00634 /* Wave subtraction for advected quantities */
00635     for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00636       if (W[i].Vx > 0.) {
00637         qb = 0.5*dtodx*(ev[NWAVE-1]-W[i].Vx);
00638         qc = 0.5*dtodx*dtodx*TWO_3RDS*(SQR(ev[NWAVE-1]) - SQR(W[i].Vx));
00639         pWl[n] += qb*(dW[m]-W6[m]) + qc*W6[m];
00640       } else if (W[i].Vx < 0.) {
00641         qb = 0.5*dtodx*(ev[0]-W[i].Vx);
00642         qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[0]*ev[0] - W[i].Vx*W[i].Vx);
00643         pWr[n] += qb*(dW[m]+W6[m]) + qc*W6[m];
00644       }
00645     }
00646 
00647 #endif /* CTU_INTEGRATOR */
00648 
00649 /*--- Step 20. -----------------------------------------------------------------
00650  * Save eigenvalues and eigenmatrices at i+1 for use in next iteration */
00651 
00652     for (m=0; m<NWAVE; m++) {
00653       ev[m] = ev_ip1[m];
00654       for (n=0; n<NWAVE; n++) {
00655         rem[m][n] = rem_ip1[m][n];
00656         lem[m][n] = lem_ip1[m][n];
00657       }
00658     }
00659 
00660   } /*====================== END BIG LOOP OVER i =========================*/
00661 
00662   return;
00663 }
00664 
00665 
00666 /*----------------------------------------------------------------------------*/
00667 /*! \fn void lr_states_init(MeshS *pM)
00668  *  \brief Allocate enough memory for work arrays */
00669 
00670 void lr_states_init(MeshS *pM)
00671 {
00672   int nmax,size1=0,size2=0,size3=0,nl,nd;
00673 
00674 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
00675   for (nl=0; nl<(pM->NLevels); nl++){
00676     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00677       if (pM->Domain[nl][nd].Grid != NULL) {
00678         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00679           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00680         }
00681         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00682           size2 = pM->Domain[nl][nd].Grid->Nx[1];
00683         }
00684         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00685           size3 = pM->Domain[nl][nd].Grid->Nx[2];
00686         }
00687       }
00688     }
00689   }
00690 
00691   size1 = size1 + 2*nghost;
00692   size2 = size2 + 2*nghost;
00693   size3 = size3 + 2*nghost;
00694   nmax = MAX((MAX(size1,size2)),size3);
00695 
00696   if ((pW = (Real**)malloc(nmax*sizeof(Real*))) == NULL) goto on_error;
00697 
00698   if ((dWm = (Real**)calloc_2d_array(nmax, NWAVE, sizeof(Real))) == NULL)
00699     goto on_error;
00700 
00701   if ((Wim1h = (Real**)calloc_2d_array(nmax, NWAVE, sizeof(Real))) == NULL)
00702     goto on_error;
00703 
00704   return;
00705   on_error:
00706     lr_states_destruct();
00707     ath_error("[lr_states_init]: malloc returned a NULL pointer\n");
00708 }
00709 
00710 /*----------------------------------------------------------------------------*/
00711 /*! \fn void lr_states_destruct(void)
00712  *  \brief Free memory used by work arrays */
00713 
00714 void lr_states_destruct(void)
00715 {
00716   if (pW != NULL) free(pW);
00717   if (dWm != NULL) free_2d_array(dWm);
00718   if (Wim1h != NULL) free_2d_array(Wim1h);
00719   return;
00720 }
00721 
00722 #endif /* THIRD_ORDER_CHAR */
reconstruction/lr_states_ppm.c
