Athena: rsolvers/hlle.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file hlle.c
00004  *  \brief Computes 1D fluxes using the Harten-Lax-van Leer (HLLE) Riemann
00005  *  solver.
00006  *
00007  * PURPOSE: Computes 1D fluxes using the Harten-Lax-van Leer (HLLE) Riemann
00008  *   solver.  This flux is very diffusive, especially for contacts, and so it
00009  *   is not recommended for use in applications.  However, as shown by Einfeldt
00010  *   et al.(1991), it is positively conservative, so it is a useful option when
00011  *   other approximate solvers fail and/or when extra dissipation is needed.
00012  *
00013  * REFERENCES:
00014  * - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics",
00015  *   2nd ed., Springer-Verlag, Berlin, (1999) chpt. 10.
00016  *
00017  * - Einfeldt et al., "On Godunov-type methods near low densities",
00018  *   JCP, 92, 273 (1991)
00019  *
00020  * - A. Harten, P. D. Lax and B. van Leer, "On upstream differencing and
00021  *   Godunov-type schemes for hyperbolic conservation laws",
00022  *   SIAM Review 25, 35-61 (1983).
00023  *
00024  * - B. Einfeldt, "On Godunov-type methods for gas dynamics",
00025  *   SIAM J. Numerical Analysis 25, 294-318 (1988).
00026  *
00027  * CONTAINS PUBLIC FUNCTIONS:
00028  * - fluxes() - all Riemann solvers in Athena must have this function name and
00029  *              use the same argument list as defined in rsolvers/prototypes.h*/
00030 /*============================================================================*/
00031 
00032 #include <math.h>
00033 #include <stdio.h>
00034 #include <stdlib.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040 
00041 #ifndef SPECIAL_RELATIVITY
00042 #ifdef HLLE_FLUX
00043 #define HLLE_FUNCTION fluxes
00044 #endif
00045 
00046 #if defined(HLLE_FLUX) || defined (ROE_FLUX)
00047 /*----------------------------------------------------------------------------*/
00048 /*! \fn void HLLE_FUNCTION(const Cons1DS Ul, const Cons1DS Ur,
00049  *                     const Prim1DS Wl, const Prim1DS Wr,
00050  *                   const Real Bxi, Cons1DS *pFlux)
00051  *  \brief HLLE flux function wrapper.
00052  *
00053  * - HLLE_FUNCTION=fluxes if HLLE_FLUX defined
00054  * - HLLE_FUNCTION=hlle_flux if ROE_FLUX defined
00055  *   Input Arguments:
00056  *  -  Bxi = B in direction of 1D slice at cell interface
00057  *  -  Ul,Ur = L/R-states of CONSERVED variables at cell interface
00058  *   Output Arguments:
00059  *  -  pFlux = pointer to fluxes of CONSERVED variables at cell interface
00060  */
00061 
00062 void HLLE_FUNCTION(const Cons1DS Ul, const Cons1DS Ur,
00063                    const Prim1DS Wl, const Prim1DS Wr,
00064                    const Real Bxi, Cons1DS *pFlux)
00065 {
00066   Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe,pbl=0.0,pbr=0.0;
00067   Real asq,vaxsq=0.0,qsq,cfsq,cfl,cfr,bp,bm,ct2=0.0,tmp;
00068 #ifndef ISOTHERMAL
00069   Real hroe;
00070 #endif
00071 #ifdef MHD
00072   Real b2roe,b3roe,x,y;
00073 #endif
00074   Real ev[NWAVE],al,ar;
00075   Real *pFl, *pFr, *pF;
00076   Cons1DS Fl,Fr;
00077   int n;
00078 
00079 
00080 /*--- Step 1. ------------------------------------------------------------------
00081  * Convert left- and right- states in conserved to primitive variables.
00082  */
00083 
00084 /*
00085   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00086   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00087 */
00088 
00089 /*--- Step 2. ------------------------------------------------------------------
00090  * Compute Roe-averaged data from left- and right-states
00091  */
00092 
00093   sqrtdl = sqrt((double)Wl.d);
00094   sqrtdr = sqrt((double)Wr.d);
00095   isdlpdr = 1.0/(sqrtdl + sqrtdr);
00096 
00097   droe  = sqrtdl*sqrtdr;
00098   v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00099   v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00100   v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00101 
00102 /* The Roe average of the magnetic field is defined differently.  */
00103 
00104 #ifdef MHD
00105   b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00106   b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00107   x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))/(SQR(sqrtdl + sqrtdr));
00108   y = 0.5*(Wl.d + Wr.d)/droe;
00109   pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00110   pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00111 #endif
00112 
00113 /*
00114  * Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
00115  * rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
00116  */
00117 
00118 #ifndef ISOTHERMAL
00119   hroe  = ((Ul.E + Wl.P + pbl)/sqrtdl + (Ur.E + Wr.P + pbr)/sqrtdr)*isdlpdr;
00120 #endif
00121 
00122 /*--- Step 3. ------------------------------------------------------------------
00123  * Compute eigenvalues using Roe-averaged values, needed in step 4.
00124  */
00125 
00126 #ifdef HYDRO
00127 #ifdef ISOTHERMAL
00128   esys_roe_iso_hyd(v1roe, v2roe, v3roe,       ev, NULL, NULL);
00129 #else
00130   esys_roe_adb_hyd(v1roe, v2roe, v3roe, hroe, ev, NULL, NULL);
00131 #endif /* ISOTHERMAL */
00132 #endif /* HYDRO */
00133 
00134 #ifdef MHD
00135 #ifdef ISOTHERMAL
00136  esys_roe_iso_mhd(droe,v1roe,v2roe,v3roe,     Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00137 #else
00138  esys_roe_adb_mhd(droe,v1roe,v2roe,v3roe,hroe,Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00139 #endif /* ISOTHERMAL */
00140 #endif /* MHD */
00141 
00142 /*--- Step 4. ------------------------------------------------------------------
00143  * Compute the max and min wave speeds
00144  */
00145 
00146 /* left state */
00147 #ifdef ISOTHERMAL
00148   asq = Iso_csound2;
00149 #else
00150   asq = Gamma*Wl.P/Wl.d;
00151 #endif
00152 #ifdef MHD
00153   vaxsq = Bxi*Bxi/Wl.d;
00154   ct2 = (Ul.By*Ul.By + Ul.Bz*Ul.Bz)/Wl.d;
00155 #endif
00156   qsq = vaxsq + ct2 + asq;
00157   tmp = vaxsq + ct2 - asq;
00158   cfsq = 0.5*(qsq + sqrt((double)(tmp*tmp + 4.0*asq*ct2)));
00159   cfl = sqrt((double)cfsq);
00160 
00161 /* right state */
00162 #ifdef ISOTHERMAL
00163   asq = Iso_csound2;
00164 #else
00165   asq = Gamma*Wr.P/Wr.d; 
00166 #endif
00167 #ifdef MHD
00168   vaxsq = Bxi*Bxi/Wr.d;
00169   ct2 = (Ur.By*Ur.By + Ur.Bz*Ur.Bz)/Wr.d;
00170 #endif
00171   qsq = vaxsq + ct2 + asq;
00172   tmp = vaxsq + ct2 - asq;
00173   cfsq = 0.5*(qsq + sqrt((double)(tmp*tmp + 4.0*asq*ct2)));
00174   cfr = sqrt((double)cfsq);
00175 
00176 /* take max/min of Roe eigenvalues and L/R state wave speeds */
00177   ar = MAX(ev[NWAVE-1],(Wr.Vx + cfr));
00178   al = MIN(ev[0]      ,(Wl.Vx - cfl));
00179 
00180   bp = MAX(ar, 0.0);
00181   bm = MIN(al, 0.0);
00182 
00183 /*--- Step 5. ------------------------------------------------------------------
00184  * Compute L/R fluxes along the lines bm/bp: F_{L}-S_{L}U_{L}; F_{R}-S_{R}U_{R}
00185  */
00186 
00187   Fl.d  = Ul.Mx - bm*Ul.d;
00188   Fr.d  = Ur.Mx - bp*Ur.d;
00189 
00190   Fl.Mx = Ul.Mx*(Wl.Vx - bm);
00191   Fr.Mx = Ur.Mx*(Wr.Vx - bp);
00192 
00193   Fl.My = Ul.My*(Wl.Vx - bm);
00194   Fr.My = Ur.My*(Wr.Vx - bp);
00195 
00196   Fl.Mz = Ul.Mz*(Wl.Vx - bm);
00197   Fr.Mz = Ur.Mz*(Wr.Vx - bp);
00198 
00199 #ifdef ISOTHERMAL
00200   Fl.Mx += Wl.d*Iso_csound2;
00201   Fr.Mx += Wr.d*Iso_csound2;
00202 #else
00203   Fl.Mx += Wl.P;
00204   Fr.Mx += Wr.P;
00205 
00206   Fl.E  = Ul.E*(Wl.Vx - bm) + Wl.P*Wl.Vx;
00207   Fr.E  = Ur.E*(Wr.Vx - bp) + Wr.P*Wr.Vx;
00208 #endif /* ISOTHERMAL */
00209 
00210 #ifdef MHD
00211   Fl.Mx -= 0.5*(Bxi*Bxi - SQR(Wl.By) - SQR(Wl.Bz));
00212   Fr.Mx -= 0.5*(Bxi*Bxi - SQR(Wr.By) - SQR(Wr.Bz));
00213 
00214   Fl.My -= Bxi*Wl.By;
00215   Fr.My -= Bxi*Wr.By;
00216     
00217   Fl.Mz -= Bxi*Wl.Bz;
00218   Fr.Mz -= Bxi*Wr.Bz;
00219 
00220 #ifndef ISOTHERMAL
00221   Fl.E += (pbl*Wl.Vx - Bxi*(Bxi*Wl.Vx + Wl.By*Wl.Vy + Wl.Bz*Wl.Vz));
00222   Fr.E += (pbr*Wr.Vx - Bxi*(Bxi*Wr.Vx + Wr.By*Wr.Vy + Wr.Bz*Wr.Vz));
00223 #endif /* ISOTHERMAL */
00224 
00225   Fl.By = Wl.By*(Wl.Vx - bm) - Bxi*Wl.Vy;
00226   Fr.By = Wr.By*(Wr.Vx - bp) - Bxi*Wr.Vy;
00227 
00228   Fl.Bz = Wl.Bz*(Wl.Vx - bm) - Bxi*Wl.Vz;
00229   Fr.Bz = Wr.Bz*(Wr.Vx - bp) - Bxi*Wr.Vz;
00230 #endif /* MHD */
00231 
00232 #if (NSCALARS > 0)
00233   for (n=0; n<NSCALARS; n++) {
00234     Fl.s[n] = Fl.d*Wl.r[n];
00235     Fr.s[n] = Fr.d*Wr.r[n];
00236   }
00237 #endif
00238 
00239 #ifdef CYLINDRICAL
00240 #ifndef ISOTHERMAL
00241   Fl.Pflux = Wl.P;
00242   Fr.Pflux = Wr.P;
00243 #ifdef MHD
00244   Fl.Pflux += pbl;
00245   Fr.Pflux += pbr;
00246 #endif /* MHD */
00247 #endif /* ISOTHERMAL */
00248 #endif /* CYLINDRICAL */
00249 
00250 /*--- Step 6. ------------------------------------------------------------------
00251  * Compute the HLLE flux at interface.
00252  */
00253 
00254   pFl = (Real *)&(Fl);
00255   pFr = (Real *)&(Fr);
00256   pF  = (Real *)pFlux;
00257   tmp = 0.5*(bp + bm)/(bp - bm);
00258   for (n=0; n<(NWAVE+NSCALARS); n++){
00259     pF[n] = 0.5*(pFl[n] + pFr[n]) + (pFl[n] - pFr[n])*tmp;
00260   }
00261 
00262 #ifdef CYLINDRICAL
00263   n = NWAVE+NSCALARS;
00264   pF[n] = 0.5*(pFl[n] + pFr[n]) + (pFl[n] - pFr[n])*tmp;
00265 #endif 
00266 
00267   return;
00268 }
00269 #endif /* HLLE_FLUX */
00270 #endif