Athena: rsolvers/force.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file force.c
00004  *  \brief Computes 1D fluxes using a Riemann solver similar, but not
00005  *   identical, to Toro's FORCE (First-ORder-CEntred) flux.
00006  *
00007  * PURPOSE: Computes 1D fluxes using a Riemann solver similar, but not
00008  *   identical, to Toro's FORCE (First-ORder-CEntred) flux.  It uses the
00009  *   average of a Lax-Wendroff and HLLE flux for sub-sonic flow, and the
00010  *   appropriate upwind flux otherwise.
00011  *
00012  * HISTORY: -- TAG -- 3/4/2005
00013  *
00014  * REFERENCES:
00015  * - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics",
00016  *   2nd ed., Springer-Verlag, Berlin, (1999), section 7.4.2.
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 FORCE_FLUX
00033 /*----------------------------------------------------------------------------*/
00034 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur, 
00035  *            const Prim1DS Wl, const Prim1DS Wr,
00036  *            const Real Bxi, Cons1DS *pFlux)
00037  *  \brief Computes 1D fluxes
00038  *   Input Arguments:
00039  *  -  Bxi = B in direction of 1D slice at cell interface
00040  *  -  Ul,Ur = L/R-states of CONSERVED variables at cell interface
00041  *   Output Arguments:
00042  *  -  pFlux = pointer to fluxes of CONSERVED variables at cell interface
00043  */
00044 
00045 void fluxes(const Cons1DS Ul, const Cons1DS Ur, 
00046             const Prim1DS Wl, const Prim1DS Wr,
00047             const Real Bxi, Cons1DS *pFlux)
00048 {
00049   Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe,pbl=0.0,pbr=0.0;
00050   Real asq,vaxsq=0.0,qsq,cfsq,cfl,cfr,bp,bm,ct2=0.0,tmp;
00051 #ifndef ISOTHERMAL
00052   Real hroe;
00053 #endif
00054 #ifdef MHD
00055   Real b2roe,b3roe,x,y;
00056 #endif
00057   Real ev[NWAVE],al,ar;
00058   Real *pFl, *pFc, *pFr, *pUc, *pF;
00059   Prim1DS Wc;
00060   Cons1DS Fl, Fc, Fr, Uc;
00061   int n;
00062 
00063 /* The first 5 steps are identical to those in hlle fluxes */
00064 /*--- Step 1. ------------------------------------------------------------------
00065  * Convert left- and right- states in conserved to primitive variables.
00066  */
00067 /*
00068   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00069   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00070 */
00071 
00072 /*--- Step 2. ------------------------------------------------------------------
00073  * Compute Roe-averaged data from left- and right-states
00074  */
00075 
00076   sqrtdl = sqrt((double)Wl.d);
00077   sqrtdr = sqrt((double)Wr.d);
00078   isdlpdr = 1.0/(sqrtdl + sqrtdr);
00079 
00080   droe  = sqrtdl*sqrtdr;
00081   v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00082   v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00083   v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00084 
00085 /*
00086  * The Roe average of the magnetic field is defined differently.
00087  */
00088 
00089 #ifdef MHD
00090   b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00091   b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00092   x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))/(SQR(sqrtdl + sqrtdr));
00093   y = 0.5*(Wl.d + Wr.d)/droe;
00094   pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00095   pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00096 #endif
00097 
00098 /*
00099  * Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
00100  * rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
00101  */
00102 
00103 #ifndef ISOTHERMAL
00104   hroe  = ((Ul.E + Wl.P + pbl)/sqrtdl + (Ur.E + Wr.P + pbr)/sqrtdr)*isdlpdr;
00105 #endif
00106 
00107 /*--- Step 3. ------------------------------------------------------------------
00108  * Compute eigenvalues using Roe-averaged values
00109  */
00110 
00111 #ifdef HYDRO
00112 #ifdef ISOTHERMAL
00113   esys_roe_iso_hyd(v1roe, v2roe, v3roe, ev, NULL, NULL);
00114 #else
00115   esys_roe_adb_hyd(v1roe, v2roe, v3roe, hroe, ev, NULL, NULL);
00116 #endif /* ISOTHERMAL */
00117 #endif /* HYDRO */
00118 
00119 #ifdef MHD
00120 #ifdef ISOTHERMAL
00121  esys_roe_iso_mhd(droe,v1roe,v2roe,v3roe,     Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00122 #else
00123  esys_roe_adb_mhd(droe,v1roe,v2roe,v3roe,hroe,Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00124 #endif /* ISOTHERMAL */
00125 #endif /* MHD */
00126 
00127 /*--- Step 4. ------------------------------------------------------------------
00128  * Compute the max and min wave speeds
00129  */
00130 
00131 /* left state */
00132 #ifdef ISOTHERMAL
00133   asq = Iso_csound2;
00134 #else
00135   asq = Gamma*Wl.P/Wl.d;
00136 #endif
00137 #ifdef MHD
00138   vaxsq = Bxi*Bxi/Wl.d;
00139   ct2 = (Ul.By*Ul.By + Ul.Bz*Ul.Bz)/Wl.d;
00140 #endif
00141   qsq = vaxsq + ct2 + asq;
00142   tmp = vaxsq + ct2 - asq;
00143   cfsq = 0.5*(qsq + sqrt((double)(tmp*tmp + 4.0*asq*ct2)));
00144   cfl = sqrt((double)cfsq);
00145 
00146 /* right state */
00147 #ifdef ISOTHERMAL
00148   asq = Iso_csound2;
00149 #else
00150   asq = Gamma*Wr.P/Wr.d; 
00151 #endif
00152 #ifdef MHD
00153   vaxsq = Bxi*Bxi/Wr.d;
00154   ct2 = (Ur.By*Ur.By + Ur.Bz*Ur.Bz)/Wr.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   cfr = sqrt((double)cfsq);
00160 
00161 /* take max/min of Roe eigenvalues and L/R state wave speeds */
00162   ar = MAX(ev[NWAVE-1],(Wr.Vx + cfr));
00163   al = MIN(ev[0]      ,(Wl.Vx - cfl));
00164 
00165   bp = MAX(ar, 0.0);
00166   bm = MIN(al, 0.0);
00167 
00168 /*--- Step 5. ------------------------------------------------------------------
00169  * Compute L/R fluxes along the line bm, bp
00170  */
00171 
00172   Fl.d  = Ul.Mx - bm*Ul.d;
00173   Fr.d  = Ur.Mx - bp*Ur.d;
00174 
00175   Fl.Mx = Ul.Mx*(Wl.Vx - bm);
00176   Fr.Mx = Ur.Mx*(Wr.Vx - bp);
00177 
00178   Fl.My = Ul.My*(Wl.Vx - bm);
00179   Fr.My = Ur.My*(Wr.Vx - bp);
00180 
00181   Fl.Mz = Ul.Mz*(Wl.Vx - bm);
00182   Fr.Mz = Ur.Mz*(Wr.Vx - bp);
00183 
00184 #ifdef ISOTHERMAL
00185   Fl.Mx += Wl.d*Iso_csound2;
00186   Fr.Mx += Wr.d*Iso_csound2;
00187 #else
00188   Fl.Mx += Wl.P;
00189   Fr.Mx += Wr.P;
00190 
00191   Fl.E  = Ul.E*(Wl.Vx - bm) + Wl.P*Wl.Vx;
00192   Fr.E  = Ur.E*(Wr.Vx - bp) + Wr.P*Wr.Vx;
00193 #endif /* ISOTHERMAL */
00194 
00195 #ifdef MHD
00196   Fl.Mx -= 0.5*(Bxi*Bxi - SQR(Wl.By) - SQR(Wl.Bz));
00197   Fr.Mx -= 0.5*(Bxi*Bxi - SQR(Wr.By) - SQR(Wr.Bz));
00198 
00199   Fl.My -= Bxi*Wl.By;
00200   Fr.My -= Bxi*Wr.By;
00201     
00202   Fl.Mz -= Bxi*Wl.Bz;
00203   Fr.Mz -= Bxi*Wr.Bz;
00204 
00205 #ifndef ISOTHERMAL
00206   Fl.E += (pbl*Wl.Vx - Bxi*(Bxi*Wl.Vx + Wl.By*Wl.Vy + Wl.Bz*Wl.Vz));
00207   Fr.E += (pbr*Wr.Vx - Bxi*(Bxi*Wr.Vx + Wr.By*Wr.Vy + Wr.Bz*Wr.Vz));
00208 #endif /* ISOTHERMAL */
00209 
00210   Fl.By = Wl.By*(Wl.Vx - bm) - Bxi*Wl.Vy;
00211   Fr.By = Wr.By*(Wr.Vx - bp) - Bxi*Wr.Vy;
00212 
00213   Fl.Bz = Wl.Bz*(Wl.Vx - bm) - Bxi*Wl.Vz;
00214   Fr.Bz = Wr.Bz*(Wr.Vx - bp) - Bxi*Wr.Vz;
00215 #endif /* MHD */
00216 
00217 #if (NSCALARS > 0)
00218   for (n=0; n<NSCALARS; n++) {
00219     Fl.s[n] = Fl.d*Wl.x[n];
00220     Fr.s[n] = Fr.d*Wr.x[n];
00221   }
00222 #endif
00223 
00224 /*--- Step 6. ------------------------------------------------------------------
00225  * For supersonic flow, return the upwind flux.
00226  */
00227 
00228   if(al >= 0.0){
00229     *pFlux = Fl;
00230     return;
00231   }
00232 
00233   if(ar <= 0.0){
00234     *pFlux = Fr;
00235     return;
00236   }
00237 
00238 /*--- Step 7. ------------------------------------------------------------------
00239  * Compute the LW flux, start with the HLL mean state */
00240 
00241   pFl = (Real *)&(Fl);
00242   pFr = (Real *)&(Fr);
00243   pUc = (Real *)&(Uc);
00244   tmp = 1.0/(ar - al);
00245   for (n=0; n<(NWAVE+NSCALARS); n++){
00246     pUc[n] = (pFl[n] - pFr[n])*tmp;
00247   }
00248 
00249 /* Convert the HLL mean state to primitive variables */
00250   Cons1D_to_Prim1D(&Uc,&Wc,&Bxi);
00251 
00252 /* Compute the LW flux along the line dx/dt = 0 */
00253 
00254   Fc.d  = Uc.Mx;
00255   Fc.Mx = Uc.Mx*Wc.Vx;
00256   Fc.My = Uc.My*Wc.Vx;
00257   Fc.Mz = Uc.Mz*Wc.Vx;
00258 
00259 #ifdef ISOTHERMAL
00260   Fc.Mx += Wc.d*Iso_csound2;
00261 #else
00262   Fc.Mx += Wc.P;
00263 
00264   Fc.E  = Uc.E*Wc.Vx + Wc.P*Wc.Vx;
00265 #endif /* ISOTHERMAL */
00266 
00267 #ifdef MHD
00268   Fc.Mx -= 0.5*(Bxi*Bxi - SQR(Wc.By) - SQR(Wc.Bz));
00269   Fc.My -= Bxi*Wc.By;
00270   Fc.Mz -= Bxi*Wc.Bz;
00271 
00272 #ifndef ISOTHERMAL
00273   Fc.E += (pbl*Wc.Vx - Bxi*(Bxi*Wc.Vx + Wc.By*Wc.Vy + Wc.Bz*Wc.Vz));
00274 #endif /* ISOTHERMAL */
00275 
00276   Fc.By = Wc.By*Wc.Vx - Bxi*Wc.Vy;
00277   Fc.Bz = Wc.Bz*Wc.Vx - Bxi*Wc.Vz;
00278 #endif /* MHD */
00279 
00280 #if (NSCALARS > 0)
00281   for (n=0; n<NSCALARS; n++) {
00282     Fc.s[n] = Fc.d*Wc.x[n];
00283   }
00284 #endif
00285 
00286 
00287 /*--- Step 8. ------------------------------------------------------------------
00288  * Compute the average of the Lax-Wendroff & HLLE flux
00289  */
00290   pFl = (Real *)&(Fl);
00291   pFc = (Real *)&(Fc);
00292   pFr = (Real *)&(Fr);
00293   pF  = (Real *)pFlux;
00294   tmp = 0.25*(bp + bm)/(bp - bm);
00295   for (n=0; n<(NWAVE+NSCALARS); n++){
00296     pF[n] = 0.5*pFc[n] + 0.25*(pFl[n] + pFr[n]) + (pFl[n] - pFr[n])*tmp;
00297   }
00298 
00299   return;
00300 }
00301 #endif /* FORCE_FLUX */