Athena: rsolvers/hllc.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file hllc.c
00004  *  \brief Computes 1D fluxes using the HLLC Riemann solver.
00005  *
00006  * PURPOSE: Computes 1D fluxes using the HLLC Riemann solver, an extension of
00007  *   the HLLE fluxes to include the contact wave.  Currently only works for
00008  *   hydrodynamics.  For an extension to MHD, see hlld.c
00009  *
00010  * REFERENCES:
00011  * - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics",
00012  *   2nd ed., Springer-Verlag, Berlin, (1999) chpt. 10.
00013  *
00014  * - P. Batten, N. Clarke, C. Lambert, and D. M. Causon,
00015  *   "On the Choice of Wavespeeds for the HLLC Riemann Solver", 
00016  *   SIAM J. Sci. & Stat. Comp. 18, 6, 1553-1570, (1997).
00017  *
00018  * CONTAINS PUBLIC FUNCTIONS: 
00019  * - fluxes() - all Riemann solvers in Athena must have this function name and
00020  *              use the same argument list as defined in rsolvers/prototypes.h*/
00021 /*============================================================================*/
00022 
00023 #include <math.h>
00024 #include <stdio.h>
00025 #include <stdlib.h>
00026 #include "../defs.h"
00027 #include "../athena.h"
00028 #include "../globals.h"
00029 #include "prototypes.h"
00030 #include "../prototypes.h"
00031 
00032 #ifdef HLLC_FLUX
00033 #ifndef SPECIAL_RELATIVITY
00034 
00035 #ifdef MHD
00036 #error : The HLLC flux only works for hydro.
00037 #endif /* MHD */
00038 
00039 /*----------------------------------------------------------------------------*/
00040 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00041  *            const Prim1DS Wl, const Prim1DS Wr,
00042  *            const Real Bxi, Cons1DS *pFlux)
00043  *  \brief Computes 1D fluxes
00044  *   Input Arguments:
00045  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface 
00046  *   Output Arguments:
00047  *   - pFlux = pointer to fluxes of CONSERVED variables at cell interface 
00048  */
00049 
00050 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00051             const Prim1DS Wl, const Prim1DS Wr,
00052             const Real Bxi, Cons1DS *pFlux)
00053 {
00054   Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe;
00055 #ifndef BAROTROPIC
00056   Real hroe;
00057 #endif
00058   Real ev[NWAVE];
00059   Real *pFl, *pFr, *pF;
00060   Cons1DS Fl,Fr;
00061   int n;
00062   Real cfl,cfr,bp,bm,tmp;
00063   Real al,ar; /* Min and Max wave speeds */
00064   Real am,cp; /* Contact wave speed and pressure */
00065   Real tl,tr,dl,dr,sl,sm,sr;
00066 
00067 /*--- Step 1. ------------------------------------------------------------------
00068  * Convert left- and right- states in conserved to primitive variables.
00069  */
00070 /*
00071   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00072   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00073 */
00074 
00075 /*--- Step 2. ------------------------------------------------------------------
00076  * Compute Roe-averaged data from left- and right-states
00077  */
00078 
00079   sqrtdl = sqrt((double)Wl.d);
00080   sqrtdr = sqrt((double)Wr.d);
00081   isdlpdr = 1.0/(sqrtdl + sqrtdr);
00082 
00083   droe  = sqrtdl*sqrtdr;
00084   v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00085   v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00086   v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00087 
00088 /*
00089  * Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
00090  * rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
00091  */
00092 
00093 #ifndef ISOTHERMAL
00094   hroe = ((Ul.E + Wl.P)/sqrtdl + (Ur.E + Wr.P)/sqrtdr)*isdlpdr;
00095 #endif
00096 
00097 /*--- Step 3. ------------------------------------------------------------------
00098  * Compute eigenvalues using Roe-averaged values
00099  */
00100 
00101 #ifdef ISOTHERMAL
00102   esys_roe_iso_hyd(v1roe, v2roe, v3roe, ev, NULL, NULL);
00103 #else
00104   esys_roe_adb_hyd(v1roe, v2roe, v3roe, hroe, ev, NULL, NULL);
00105 #endif /* ISOTHERMAL */
00106 
00107 /*--- Step 4. ------------------------------------------------------------------
00108  * Compute the max and min wave speeds
00109  */
00110 
00111 #ifdef ISOTHERMAL
00112   cfl = cfr = Iso_csound;
00113 #else
00114   cfl = sqrt((double)(Gamma*Wl.P/Wl.d));
00115   cfr = sqrt((double)(Gamma*Wr.P/Wr.d));
00116 #endif
00117 
00118   ar = MAX(ev[NWAVE-1],(Wr.Vx + cfr));
00119   al = MIN(ev[0]      ,(Wl.Vx - cfl));
00120 
00121   bp = ar > 0.0 ? ar : 0.0;
00122   bm = al < 0.0 ? al : 0.0;
00123 
00124 /*--- Step 5. ------------------------------------------------------------------
00125  * Compute the contact wave speed and Pressure
00126  */
00127 
00128 #ifdef ISOTHERMAL
00129   tl = Wl.d*Iso_csound2 + (Wl.Vx - al)*Ul.Mx;
00130   tr = Wr.d*Iso_csound2 + (Wr.Vx - ar)*Ur.Mx;
00131 #else
00132   tl = Wl.P + (Wl.Vx - al)*Ul.Mx;
00133   tr = Wr.P + (Wr.Vx - ar)*Ur.Mx;
00134 #endif
00135 
00136   dl =   Ul.Mx - Ul.d*al;
00137   dr = -(Ur.Mx - Ur.d*ar);
00138 
00139   tmp = 1.0/(dl + dr);
00140 /* Determine the contact wave speed... */
00141   am = (tl - tr)*tmp;
00142 /* ...and the pressure at the contact surface */
00143   cp = (dl*tr + dr*tl)*tmp;
00144   if(cp < 0.0) ath_perr(1,"[hllc flux]: Contact Pressure = %g\n",cp);
00145   cp = cp > 0.0 ? cp : 0.0;
00146 
00147 /*--- Step 6. ------------------------------------------------------------------
00148  * Compute L/R fluxes along the line bm, bp
00149  */
00150 
00151   Fl.d  = Ul.Mx - bm*Ul.d;
00152   Fr.d  = Ur.Mx - bp*Ur.d;
00153 
00154   Fl.Mx = Ul.Mx*(Wl.Vx - bm);
00155   Fr.Mx = Ur.Mx*(Wr.Vx - bp);
00156 
00157   Fl.My = Ul.My*(Wl.Vx - bm);
00158   Fr.My = Ur.My*(Wr.Vx - bp);
00159 
00160   Fl.Mz = Ul.Mz*(Wl.Vx - bm);
00161   Fr.Mz = Ur.Mz*(Wr.Vx - bp);
00162 
00163 #ifdef ISOTHERMAL
00164   Fl.Mx += Wl.d*Iso_csound2;
00165   Fr.Mx += Wr.d*Iso_csound2;
00166 #else
00167   Fl.Mx += Wl.P;
00168   Fr.Mx += Wr.P;
00169 
00170   Fl.E  = Ul.E*(Wl.Vx - bm) + Wl.P*Wl.Vx;
00171   Fr.E  = Ur.E*(Wr.Vx - bp) + Wr.P*Wr.Vx;
00172 #endif /* ISOTHERMAL */
00173 
00174 #if (NSCALARS > 0)
00175   for (n=0; n<NSCALARS; n++) {
00176     Fl.s[n] = Fl.d*Wl.r[n];
00177     Fr.s[n] = Fr.d*Wr.r[n];
00178   }
00179 #endif
00180 
00181 /*--- Step 7. ------------------------------------------------------------------
00182  * Compute flux weights or scales
00183  */
00184 
00185   if (am >= 0.0) {
00186     sl =  am/(am - bm);
00187     sr = 0.0;
00188     sm = -bm/(am - bm);
00189   }
00190   else {
00191     sl =  0.0;
00192     sr = -am/(bp - am);
00193     sm =  bp/(bp - am);
00194   }
00195 
00196 /*--- Step 8. ------------------------------------------------------------------
00197  * Compute the HLLC flux at interface
00198  */
00199   pFl = (Real *)&(Fl);
00200   pFr = (Real *)&(Fr);
00201   pF  = (Real *)pFlux;
00202   for (n=0; n<(NWAVE+NSCALARS); n++) pF[n] = sl*pFl[n] + sr*pFr[n];
00203 
00204 /* Add the weighted contribution of the flux along the contact */
00205   pFlux->Mx += sm*cp;
00206 #ifndef ISOTHERMAL
00207   pFlux->E  += sm*cp*am;
00208 #endif /* ISOTHERMAL */
00209 
00210 #ifdef CYLINDRICAL
00211   if (al > 0.0) {
00212 #ifndef ISOTHERMAL
00213     pFlux->Pflux = Wl.P;
00214 #else /* ISOTHERMAL */
00215     pFlux->Pflux = Wl.d*Iso_csound2;
00216 #endif /* ISOTHERMAL */
00217   }
00218   else if (ar < 0.0) {
00219 #ifndef ISOTHERMAL
00220     pFlux->Pflux = Wr.P;
00221 #else /* ISOTHERMAL */
00222     pFlux->Pflux = Wr.d*Iso_csound2;
00223 #endif /* ISOTHERMAL */
00224   }
00225   else {
00226 #ifndef ISOTHERMAL
00227     pFlux->Pflux = cp;
00228 #else /* ISOTHERMAL */
00229     if (am >= 0.0) {
00230       pFlux->Pflux = Wl.d*(al-Wl.Vx)/(al-am);
00231     }
00232     else {
00233       pFlux->Pflux = Wr.d*(ar-Wr.Vx)/(ar-am);
00234     }
00235 #endif /* ISOTHERMAL */
00236   }
00237 #endif /* CYLINDRICAL */
00238 
00239   return;
00240 }
00241 #endif /* SPECIAL_RELATIVITY */
00242 #endif /* HLLC_FLUX */