Athena: rsolvers/roe.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file roe.c
00004  *  \brief Computes 1D fluxes using Roe's linearization.
00005  *
00006  * PURPOSE: Computes 1D fluxes using Roe's linearization.  When Roe's method
00007  * fails because of negative density or pressure in the intermediate states,
00008  * the fluxes are computed with the HLLE solver instead.
00009  *
00010  * REFERENCES:
00011  * - P. Roe, "Approximate Riemann solvers, parameter vectors, and difference
00012  *   schemes", JCP, 43, 357 (1981).
00013  *
00014  * CONTAINS PUBLIC FUNCTIONS: 
00015  * - fluxes() - all Riemann solvers in Athena must have this function name and
00016  *              use the same argument list as defined in rsolvers/prototypes.h
00017  * - HLLE_FUNCTION - since the HLLE solver is requird in conjunction with the
00018  *     Roe solver, the HLLE solver is included from the file hlle.c at the
00019  *     end of this file, and the HLLE flux function is renamed to be different
00020  *     than "fluxes" (the name used in Athena for R-solvers).                 */
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 /* maximum wavespeed used by H-correction, value passed from integrator */
00033 Real etah=0.0;
00034 
00035 #ifdef ROE_FLUX
00036 /*! \fn void flux_hlle(const Cons1DS Ul, const Cons1DS Ur,
00037  *               const Prim1DS Wl, const Prim1DS Wr,
00038  *               const Real Bxi, Cons1DS *pFlux);
00039  *  \brief Function prototype for HLLE fluxes */
00040 void flux_hlle(const Cons1DS Ul, const Cons1DS Ur,
00041                const Prim1DS Wl, const Prim1DS Wr,
00042                const Real Bxi, Cons1DS *pFlux);
00043 
00044 /* Test the intermediate states in the approximate Riemann solution. */
00045 #define TEST_INTERMEDIATE_STATES
00046 
00047 /*----------------------------------------------------------------------------*/
00048 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00049  *            const Prim1DS Wl, const Prim1DS Wr,
00050  *            const Real Bxi, Cons1DS *pFlux)
00051  *  \brief Computes 1D fluxes
00052  *   Input Arguments:
00053  *   - Bxi = B in direction of slice at cell interface
00054  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00055  *   Output Arguments:
00056  *   - pFlux = pointer to fluxes of CONSERVED variables at cell interface
00057  */
00058 
00059 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00060             const Prim1DS Wl, const Prim1DS Wr,
00061             const Real Bxi, Cons1DS *pFlux)
00062 {
00063   Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe,pbl=0.0,pbr=0.0;
00064 #ifndef ISOTHERMAL
00065   Real hroe;
00066 #endif /* ISOTHERMAL */
00067 #ifdef MHD
00068   Real b2roe,b3roe,x,y;
00069 #endif
00070   Real coeff[NWAVE];
00071   Real ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00072   Real dU[NWAVE],a[NWAVE];
00073 #ifdef TEST_INTERMEDIATE_STATES
00074   Real u_inter[NWAVE];
00075 #ifdef ADIABATIC
00076   Real p_inter=0.0;
00077 #endif
00078 #endif /* TEST_INTERMEDIATE_STATES */
00079   Real *pUl, *pUr, *pFl, *pFr, *pF;
00080   Cons1DS Fl,Fr;
00081   int n,m,hlle_flag;
00082 #ifdef CYLINDRICAL
00083   Real Eint,Emag,Ekin,coeff2[NWAVE];
00084 #endif
00085 
00086   for (n=0; n<NWAVE; n++) {
00087     for (m=0; m<NWAVE; m++) {
00088       rem[n][m] = 0.0;
00089       lem[n][m] = 0.0;
00090     }
00091   }
00092 
00093 /*--- Step 1. ------------------------------------------------------------------
00094  * Convert left- and right- states in conserved to primitive variables.
00095  */
00096 
00097 /*
00098   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00099   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00100 */
00101 
00102 /*--- Step 2. ------------------------------------------------------------------
00103  * Compute Roe-averaged data from left- and right-states
00104  */
00105 
00106   sqrtdl = sqrt((double)Wl.d);
00107   sqrtdr = sqrt((double)Wr.d);
00108   isdlpdr = 1.0/(sqrtdl + sqrtdr);
00109 
00110   droe  = sqrtdl*sqrtdr;
00111   v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00112   v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00113   v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00114 
00115 /* The Roe average of the magnetic field is defined differently  */
00116 
00117 #ifdef MHD
00118   b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00119   b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00120   x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))/(SQR(sqrtdl + sqrtdr));
00121   y = 0.5*(Wl.d + Wr.d)/droe;
00122   pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00123   pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00124 #endif
00125 
00126 /*
00127  * Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
00128  * rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
00129  */
00130 
00131 #ifndef ISOTHERMAL
00132   hroe  = ((Ul.E + Wl.P + pbl)/sqrtdl + (Ur.E + Wr.P + pbr)/sqrtdr)*isdlpdr;
00133 #endif
00134 
00135 /*--- Step 3. ------------------------------------------------------------------
00136  * Compute eigenvalues and eigenmatrices using Roe-averaged values
00137  */
00138 
00139 #ifdef HYDRO
00140 #ifdef ISOTHERMAL
00141   esys_roe_iso_hyd(v1roe,v2roe,v3roe,     ev,rem,lem);
00142 #else
00143   esys_roe_adb_hyd(v1roe,v2roe,v3roe,hroe,ev,rem,lem);
00144 #endif /* ISOTHERMAL */
00145 #endif /* HYDRO */
00146 
00147 #ifdef MHD
00148 #ifdef ISOTHERMAL
00149   esys_roe_iso_mhd(droe,v1roe,v2roe,v3roe,     Bxi,b2roe,b3roe,x,y,ev,rem,lem);
00150 #else
00151   esys_roe_adb_mhd(droe,v1roe,v2roe,v3roe,hroe,Bxi,b2roe,b3roe,x,y,ev,rem,lem);
00152 #endif /* ISOTHERMAL */
00153 #endif /* MHD */
00154 
00155 /*--- Step 4. ------------------------------------------------------------------
00156  * Compute L/R fluxes 
00157  */
00158 
00159   Fl.d  = Ul.Mx;
00160   Fr.d  = Ur.Mx;
00161 
00162   Fl.Mx = Ul.Mx*Wl.Vx;
00163   Fr.Mx = Ur.Mx*Wr.Vx;
00164 
00165   Fl.My = Ul.Mx*Wl.Vy;
00166   Fr.My = Ur.Mx*Wr.Vy;
00167 
00168   Fl.Mz = Ul.Mx*Wl.Vz;
00169   Fr.Mz = Ur.Mx*Wr.Vz;
00170 
00171 #ifdef ISOTHERMAL
00172   Fl.Mx += Wl.d*Iso_csound2;
00173   Fr.Mx += Wr.d*Iso_csound2;
00174 #else
00175   Fl.Mx += Wl.P;
00176   Fr.Mx += Wr.P;
00177 
00178   Fl.E  = (Ul.E + Wl.P)*Wl.Vx;
00179   Fr.E  = (Ur.E + Wr.P)*Wr.Vx;
00180 #endif /* ISOTHERMAL */
00181 
00182 #ifdef MHD
00183   Fl.Mx -= 0.5*(Bxi*Bxi - SQR(Wl.By) - SQR(Wl.Bz));
00184   Fr.Mx -= 0.5*(Bxi*Bxi - SQR(Wr.By) - SQR(Wr.Bz));
00185 
00186   Fl.My -= Bxi*Wl.By;
00187   Fr.My -= Bxi*Wr.By;
00188     
00189   Fl.Mz -= Bxi*Wl.Bz;
00190   Fr.Mz -= Bxi*Wr.Bz;
00191 
00192 #ifndef ISOTHERMAL
00193   Fl.E += (pbl*Wl.Vx - Bxi*(Bxi*Wl.Vx + Wl.By*Wl.Vy + Wl.Bz*Wl.Vz));
00194   Fr.E += (pbr*Wr.Vx - Bxi*(Bxi*Wr.Vx + Wr.By*Wr.Vy + Wr.Bz*Wr.Vz));
00195 #endif /* ISOTHERMAL */
00196 
00197   Fl.By = Wl.By*Wl.Vx - Bxi*Wl.Vy;
00198   Fr.By = Wr.By*Wr.Vx - Bxi*Wr.Vy;
00199 
00200   Fl.Bz = Wl.Bz*Wl.Vx - Bxi*Wl.Vz;
00201   Fr.Bz = Wr.Bz*Wr.Vx - Bxi*Wr.Vz;
00202 #endif /* MHD */
00203 
00204 #if (NSCALARS > 0)
00205   for (n=0; n<NSCALARS; n++) {
00206     Fl.s[n] = Fl.d*Wl.r[n];
00207     Fr.s[n] = Fr.d*Wr.r[n];
00208   }
00209 #endif
00210 
00211 /*--- Step 5. ------------------------------------------------------------------
00212  * Return upwind flux if flow is supersonic
00213  */
00214 
00215   if(ev[0] >= 0.0){
00216     *pFlux = Fl;
00217 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00218     pFlux->Pflux = Wl.P;
00219 #ifdef MHD
00220     pFlux->Pflux += pbl;
00221 #endif /* MHD */
00222 #endif /* CYLINDRICAL && NOT BAROTROPIC */
00223     return;
00224   }
00225 
00226   if(ev[NWAVE-1] <= 0.0){
00227     *pFlux = Fr;
00228 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00229     pFlux->Pflux = Wl.P;
00230 #ifdef MHD
00231     pFlux->Pflux += pbr;
00232 #endif /* MHD */
00233 #endif /* CYLINDRICAL && NOT BAROTROPIC */
00234     return;
00235   }
00236 
00237 /*--- Step 6. ------------------------------------------------------------------
00238  * Compute projection of dU onto L eigenvectors ("vector A")
00239  */
00240 
00241   pUr = (Real *)&(Ur);
00242   pUl = (Real *)&(Ul);
00243 
00244   for (n=0; n<NWAVE; n++) dU[n] = pUr[n] - pUl[n];
00245   for (n=0; n<NWAVE; n++) {
00246     a[n] = 0.0;
00247     for (m=0; m<NWAVE; m++) a[n] += lem[n][m]*dU[m];
00248   }
00249 
00250 /*--- Step 7. ------------------------------------------------------------------
00251  * Check that the density and pressure in the intermediate states are positive.
00252  * If not, set hlle_flag=1 if d_inter<0; hlle_flag=2 if p_inter<0, get HLLE
00253  * fluxes, and return
00254  */
00255 
00256   hlle_flag = 0;
00257 #ifdef TEST_INTERMEDIATE_STATES
00258 
00259   for (n=0; n<NWAVE; n++) u_inter[n] = pUl[n];
00260   for (n=0; n<NWAVE-1; n++) {
00261     for (m=0; m<NWAVE; m++) u_inter[m] += a[n]*rem[m][n];
00262     if(ev[n+1] > ev[n]) {
00263       if (u_inter[0] <= 0.0) {
00264         hlle_flag=1;
00265         break;
00266       }
00267 #ifdef ADIABATIC
00268       p_inter = u_inter[4] - 0.5*
00269         (SQR(u_inter[1])+SQR(u_inter[2])+SQR(u_inter[3]))/u_inter[0];
00270 #ifdef MHD
00271       p_inter -= 0.5*(SQR(u_inter[NWAVE-2])+SQR(u_inter[NWAVE-1])+SQR(Bxi));
00272 #endif
00273       if (p_inter < 0.0) {
00274         hlle_flag=2;
00275         break;
00276       }
00277 #endif
00278     }
00279   }
00280 
00281   if (hlle_flag != 0) {
00282     flux_hlle(Ul,Ur,Wl,Wr,Bxi,pFlux);
00283     return;
00284   }
00285 
00286 #endif /* TEST_INTERMEDIATE_STATES */
00287 
00288 /*--- Step 8. ------------------------------------------------------------------
00289  * Compute Roe flux */
00290 
00291   pFl = (Real *)&(Fl);
00292   pFr = (Real *)&(Fr);
00293   pF  = (Real *)(pFlux); 
00294 
00295   for (m=0; m<NWAVE; m++) {
00296     coeff[m] = 0.5*MAX(fabs(ev[m]),etah)*a[m];
00297 #ifdef CYLINDRICAL
00298     coeff2[m] = 0.5*SIGN(ev[m])*a[m];
00299 #endif
00300   }
00301   for (n=0; n<NWAVE; n++) {
00302     pF[n] = 0.5*(pFl[n] + pFr[n]);
00303 #ifdef CYLINDRICAL
00304     u_inter[n] = 0.5*(pUl[n] + pUr[n]);
00305 #endif
00306     for (m=0; m<NWAVE; m++) {
00307       pF[n] -= coeff[m]*rem[n][m];
00308 #ifdef CYLINDRICAL
00309       u_inter[n] -= coeff2[m]*rem[n][m];
00310 #endif
00311     }
00312   }
00313 
00314 /* Fluxes of passively advected scalars, computed from density flux */
00315 #if (NSCALARS > 0)
00316   if (pFlux->d >= 0.0) {
00317     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00318   } else {
00319     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00320   }
00321 #endif
00322 
00323 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00324 /* "TOTAL PRESSURE FLUX" COMPUTED FROM AVERAGED STAR-STATES */
00325   Emag = 0.0;
00326 #ifdef MHD
00327   Emag = 0.5*(SQR(u_inter[NWAVE-2])+SQR(u_inter[NWAVE-1])+SQR(Bxi));
00328 #endif /* MHD */
00329   Ekin = 0.5*(SQR(u_inter[1])+SQR(u_inter[2])+SQR(u_inter[3]))/u_inter[0];
00330   Eint = u_inter[4] - Emag - Ekin;
00331 
00332   pFlux->Pflux = Eint*Gamma_1 + Emag;
00333 #endif /* CYLINDRICAL && NOT BAROTROPIC */
00334 
00335   return;
00336 }
00337 
00338 /* Include HLLE solver below, and rename HLLE flux function name */
00339 #define HLLE_FUNCTION flux_hlle
00340 #include "hlle.c"
00341 #undef HLLE_FUNCTION
00342 
00343 #endif /* ROE_FLUX */