Athena: rsolvers/exact.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file exact.c
00004  *  \brief Computes 1D fluxes using exact nonlinear Riemann solver.
00005  *
00006  * PURPOSE: Computes 1D fluxes using exact nonlinear Riemann solver.
00007  *   Currently only isothermal hydrodynamics has been implemented.  
00008  *
00009  * REFERENCES:
00010  * - R.J. LeVeque, "Numerical Methods for Conservation Laws", 2nd ed.,
00011  *   Birkhauser Verlag, Basel, (1992).
00012  *
00013  * - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics",
00014  *   2nd ed., Springer-Verlag, Berlin, (1999).
00015  *
00016  * HISTORY:
00017  * - dec-2006  Isothermal hydro version written by Nicole Lemaster
00018  * - mar-2010  Adiabatic hydro version written by Nick Hand
00019  *
00020  * CONTAINS PUBLIC FUNCTIONS:
00021  * - fluxes() - all Riemann solvers in Athena must have this function name and
00022  *              use the same argument list as defined in rsolvers/prototypes.h*/
00023 /*============================================================================*/
00024 
00025 #include <math.h>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <float.h>
00029 #include "../defs.h"
00030 #include "../athena.h"
00031 #include "../globals.h"
00032 #include "prototypes.h"
00033 #include "../prototypes.h"
00034 
00035 #ifdef EXACT_FLUX
00036 #ifndef SPECIAL_RELATIVITY
00037 
00038 #ifdef MHD
00039 #error : The exact flux for MHD has not been implemented.
00040 #endif /* MHD */
00041 
00042 #if (NSCALARS > 0)
00043 #error : Passive scalars have not been implemented in the exact flux.
00044 #endif /* NSCALARS */
00045 
00046 #ifdef ISOTHERMAL
00047 
00048 static void srder(double dm, double vl, double vr, double dmin, double dmax, 
00049                   double *y, double *dydx);
00050 static double rtsafe(void (*funcd)(double, double, double, double, double,
00051                                    double *, double *), double x1, double x2, double xacc, 
00052                      double vl, double vr, double dmin, double dmax);
00053 
00054 /*----------------------------------------------------------------------------*/
00055 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00056  *            const Prim1DS Wl, const Prim1DS Wr,
00057  *            const Real Bxi, Cons1DS *pF)
00058  *  \brief Calculates fluxes of CONSERVED variables
00059  *
00060  *   Input Arguments:
00061  *   - Bxi = B in direction of slice at cell interface
00062  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00063  *   Output Arguments:
00064  *   - pFlux = pointer to fluxes of CONSERVED variables at cell interface
00065  */
00066 
00067 
00068 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00069             const Prim1DS Wl, const Prim1DS Wr,
00070             const Real Bxi, Cons1DS *pF)
00071 {
00072   Real zl, zr, zm, dm, Vxm, Mxm, tmp, dmin, dmax;
00073   Real sl, sr;    /* Left and right going shock velocity */
00074   Real hdl, hdr;  /* Left and right going rarefaction head velocity */
00075   Real tll, tlr;  /* Left and right going rarefaction tail velocity */
00076   char soln;      /* two bits: 0=shock, 1=raref */
00077 
00078   if(!(Ul.d > 0.0)||!(Ur.d > 0.0))
00079     ath_error("[exact flux]: Non-positive densities: dl = %e  dr = %e\n",
00080               Ul.d, Ur.d);
00081 
00082 /*--- Step 1. ------------------------------------------------------------------
00083  * Compute the density and momentum of the intermediate state
00084  */
00085 
00086   zl = sqrt((double)Wl.d);
00087   zr = sqrt((double)Wr.d);
00088 
00089   /* --- 1-shock and 2-shock --- */
00090   soln = 0;
00091 
00092   /* Start by finding density if shocks on both left and right.
00093    * This will only be the case if dm > Wl.d and dm > Wr.d */
00094   tmp = zl*zr*(Wl.Vx - Wr.Vx)/(2.0*Iso_csound*(zl + zr));
00095   zm = tmp + sqrt((double)(tmp*tmp + zl*zr));
00096   dm = zm*zm;
00097 
00098   /* Get velocity from 1-shock formula */
00099   Vxm = Wl.Vx - Iso_csound*(dm-Wl.d)/(zm*zl);
00100 
00101   /* If left or right density is greater than intermediate density,
00102    * then at least one side has rarefaction instead of shock */
00103   dmin = MIN(Wl.d, Wr.d);
00104   dmax = MAX(Wl.d, Wr.d);
00105   if (dm < dmax) {
00106     /* --- 1-rarefaction and 2-rarefaction --- */
00107     soln = 3;
00108 
00109     /* Try rarefactions on both left and right, since it's a quicker
00110      * calculation than 1-shock+2-raref or 1-raref+2-shock */
00111     dm = zl*zr*exp((Wl.Vx-Wr.Vx)/(2.0*Iso_csound));
00112 
00113     /* Get velocity from 1-rarefaction formula */
00114     Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00115 
00116     /* If left or right density is smaller than intermediate density,
00117      * then we must instead have a combination of shock and rarefaction */
00118     if (dm > dmin) {
00119       /* --- EITHER 1-rarefaction and 2-shock ---
00120        * --- OR     1-shock and 2-rarefaction --- */
00121 
00122       /* Solve iteratively equation for shock and rarefaction
00123        * If Wl.d > Wr.d ==> 1-rarefaction and 2-shock
00124        * If Wr.d > Wl.d ==> 1-shock and 2-rarefaction */
00125       if (Wl.d > Wr.d) soln = 2; else soln = 1;
00126 
00127       dm = rtsafe(&srder,dmin,dmax, 2.0*DBL_EPSILON, Wl.Vx, Wr.Vx, dmin, dmax);
00128 
00129       /* Don't be foolish enough to take ln of zero */
00130       if ((dm > dmin) && (dm <= dmax)) {
00131         if (Wl.d > Wr.d) {
00132           /* Get velocity from 1-rarefaction formula */
00133           Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00134         } else {
00135           /* Get velocity from 2-rarefaction formula */
00136           Vxm = Wr.Vx + Iso_csound*log(dm/Wr.d);
00137         }
00138       } else {
00139         /* --- DEFAULT 1-rarefaction and 2-rarefaction --- */
00140         soln = 3;
00141 
00142         /* In the event that the intermediate density fails to fall between
00143          * the left and right densities (should only happen when left and
00144          * right densities differ only slightly and intermediate density
00145          * calculated in any step has significant truncation and/or roundoff
00146          * errors), default to rarefactions on both left and right */
00147         dm = zl*zr*exp((Wl.Vx-Wr.Vx)/(2.0*Iso_csound));
00148 
00149         /* Get velocity from 1-rarefaction formula */
00150         Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00151       }
00152     }
00153   }
00154 
00155   if (dm < 0.0)
00156     ath_error("[exact flux]: Solver finds negative density %5.4e\n", dm);
00157 
00158 /*--- Step 2. ------------------------------------------------------------------
00159  * Calculate the Interface Flux if the wave speeds are such that we aren't
00160  * actually in the intermediate state
00161  */
00162 
00163   if (soln & 2) { /* left rarefaction */
00164     /* The L-going rarefaction head/tail velocity */
00165     hdl = Wl.Vx - Iso_csound;
00166     tll = Vxm - Iso_csound;
00167 
00168     if (hdl >= 0.0) {
00169       /* To left of rarefaction */
00170       pF->d  = Ul.Mx;
00171       pF->Mx = Ul.Mx*(Wl.Vx) + Wl.d*Iso_csound2;
00172       pF->My = Ul.My*(Wl.Vx);
00173       pF->Mz = Ul.Mz*(Wl.Vx);
00174       return;
00175     } else if (tll >= 0.0) {
00176       /* Inside rarefaction fan */
00177       dm = Ul.d*exp(hdl/Iso_csound);
00178       Mxm = Ul.d*Iso_csound*exp(hdl/Iso_csound);
00179       Vxm = (dm == 0.0 ? 0.0 : Mxm / dm);
00180 
00181       pF->d  = Mxm;
00182       pF->Mx = Mxm*Vxm + dm*Iso_csound2;
00183       pF->My = Mxm*Wl.Vy;
00184       pF->Mz = Mxm*Wl.Vz;
00185       return;
00186     }
00187   } else { /* left shock */
00188     /* The L-going shock velocity */
00189     sl = Wl.Vx - Iso_csound*sqrt(dm)/zl;
00190 
00191     if(sl >= 0.0) {
00192       /* To left of shock */
00193       pF->d  = Ul.Mx;
00194       pF->Mx = Ul.Mx*(Wl.Vx) + Wl.d*Iso_csound2;
00195       pF->My = Ul.My*(Wl.Vx);
00196       pF->Mz = Ul.Mz*(Wl.Vx);
00197       return;
00198     }
00199   }
00200 
00201   if (soln & 1) { /* right rarefaction */
00202     /* The R-going rarefaction head/tail velocity */
00203     hdr = Wr.Vx + Iso_csound;
00204     tlr = Vxm + Iso_csound;
00205 
00206     if (hdr <= 0.0) {
00207       /* To right of rarefaction */
00208       pF->d  = Ur.Mx;
00209       pF->Mx = Ur.Mx*(Wr.Vx) + Wr.d*Iso_csound2;
00210       pF->My = Ur.My*(Wr.Vx);
00211       pF->Mz = Ur.Mz*(Wr.Vx);
00212       return;
00213     } else if (tlr <= 0.0) {
00214       /* Inside rarefaction fan */
00215       tmp = dm;
00216       dm = tmp*exp(-tlr/Iso_csound);
00217       Mxm = -tmp*Iso_csound*exp(-tlr/Iso_csound);
00218       Vxm = (dm == 0.0 ? 0.0 : Mxm / dm);
00219 
00220       pF->d  = Mxm;
00221       pF->Mx = Mxm*Vxm + dm*Iso_csound2;
00222       pF->My = Mxm*Wr.Vy;
00223       pF->Mz = Mxm*Wr.Vz;
00224       return;
00225     }
00226   } else { /* right shock */
00227     /* The R-going shock velocity */
00228     sr = Wr.Vx + Iso_csound*sqrt(dm)/zr;
00229 
00230     if(sr <= 0.0) {
00231       /* To right of shock */
00232       pF->d  = Ur.Mx;
00233       pF->Mx = Ur.Mx*(Wr.Vx) + Wr.d*Iso_csound2;
00234       pF->My = Ur.My*(Wr.Vx);
00235       pF->Mz = Ur.Mz*(Wr.Vx);
00236       return;
00237     }
00238   }
00239 
00240 /* If we make it this far, then we're in the intermediate state */
00241 
00242 /*--- Step 3. ------------------------------------------------------------------
00243  * Calculate the Interface Flux */
00244 
00245   if(Vxm >= 0.0){
00246     pF->d  = dm*Vxm;
00247     pF->Mx = dm*Vxm*Vxm + dm*Iso_csound2;
00248     pF->My = dm*Vxm*Wl.Vy;
00249     pF->Mz = dm*Vxm*Wl.Vz;
00250   }
00251   else{
00252     pF->d  = dm*Vxm;
00253     pF->Mx = dm*Vxm*Vxm + dm*Iso_csound2;
00254     pF->My = dm*Vxm*Wr.Vy;
00255     pF->Mz = dm*Vxm*Wr.Vz;
00256   }
00257 
00258   return;
00259 }
00260 
00261 /*! \fn static void srder(double dm, double vl, double vr, double dmin, 
00262  *                        double dmax, double *y, double *dydx)
00263  *  \brief Equation to solve iteratively for shock and rarefaction as well
00264  * as its derivative.  Used by rtsafe() */
00265 
00266 static void srder(double dm, double vl, double vr, double dmin, double dmax, 
00267                         double *y, double *dydx)
00268 {
00269   *y = (vr - vl) + Iso_csound*(log(dm/dmax) + (dm-dmin)/sqrt(dm*dmin));
00270   *dydx = Iso_csound/dm*(1.0 + 0.5*(dm+dmin)/sqrt(dm*dmin));
00271 
00272   return;
00273 }
00274 
00275 /*! \fn static double rtsafe(void (*funcd)(double, double, double, double, 
00276  *                          double, double *, double *), double x1, double x2, 
00277  *                          double xacc, double vl, double vr, double dmin, 
00278  *                          double dmax)
00279  *  \brief Numerical Recipes function rtsafe modified to use doubles and take
00280  * extra parameters to pass to the derivative function above */
00281 
00282 #define MAXIT 100
00283 
00284 static double rtsafe(void (*funcd)(double, double, double, double, double,
00285                         double *, double *), double x1, double x2, double xacc, 
00286                         double vl, double vr, double dmin, double dmax)
00287 {
00288         int j;
00289         double df,dx,dxold,f,fh,fl;
00290         double temp,xh,xl,rts;
00291 
00292         (*funcd)(x1,vl,vr,dmin,dmax,&fl,&df);
00293         (*funcd)(x2,vl,vr,dmin,dmax,&fh,&df);
00294         if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) {
00295                 return 0.0;
00296         }
00297         if (fl == 0.0) return x1;
00298         if (fh == 0.0) return x2;
00299         if (fl < 0.0) {
00300                 xl=x1;
00301                 xh=x2;
00302         } else {
00303                 xh=x1;
00304                 xl=x2;
00305         }
00306         rts=0.5*(x1+x2);
00307         dxold=fabs(x2-x1);
00308         dx=dxold;
00309         (*funcd)(rts,vl,vr,dmin,dmax,&f,&df);
00310         for (j=1;j<=MAXIT;j++) {
00311                 if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
00312                         || (fabs(2.0*f) > fabs(dxold*df))) {
00313                         dxold=dx;
00314                         dx=0.5*(xh-xl);
00315                         rts=xl+dx;
00316                                 if (xl == rts) return rts;
00317                 } else {
00318                         dxold=dx;
00319                         dx=f/df;
00320                         temp=rts;
00321                         rts -= dx;
00322                         if (temp == rts) return rts;
00323                 }
00324                 if (fabs(dx) < xacc) return rts;
00325                 (*funcd)(rts,vl,vr,dmin,dmax,&f,&df);
00326                 if (f < 0.0)
00327                         xl=rts;
00328                 else
00329                         xh=rts;
00330         }
00331         /* Cap the number of iterations but don't fail */
00332         return rts;
00333 }
00334 
00335 #undef MAXIT 
00336 
00337 #elif defined ADIABATIC
00338 
00339 /*------------------------------------------------------------------*/
00340 /*! \fn static Real guessP(const Prim1DS Wl, const Prim1DS Wr)
00341  *  \brief  Provides a guess value for the pressure in the center region. 
00342  *
00343  * The choice is made according to the Two-Shock approximate
00344  * Riemman solver algorithm. Returns pressure guess value p_0, or 
00345  *   TOL if p_0 is less than zero */
00346  
00347 static Real guessP(const Prim1DS Wl, const Prim1DS Wr)
00348 {
00349  
00350   Real al = sqrt((double) Gamma*Wl.P/Wl.d);   /* left sound speed */
00351   Real ar = sqrt((double) Gamma*Wr.P/Wr.d);   /* right sound speed */
00352   Real ppv;  /* inital guess for pressure */
00353   Real gl, gr, p_0;
00354   Real TOL = 1.0e-6; 
00355 
00356   ppv = 0.5*(Wl.P + Wr.P)-0.125*(Wr.Vx-Wl.Vx)*(Wl.d+Wr.d)*(al+ar);
00357 
00358   if (ppv < 0.0)
00359     ppv = 0.0;
00360 
00361   gl = sqrt((double) (2.0/(Wl.d*(Gamma+1)))/
00362             ((Gamma-1)*Wl.P/(Gamma+1) + ppv)); 
00363   gr = sqrt((double) (2.0/(Wr.d*(Gamma+1)))/
00364             ((Gamma-1)*Wr.P/(Gamma+1) + ppv));
00365   p_0 = (gl*Wl.P + gr*Wr.P - (Wr.Vx-Wl.Vx))/(gr + gl); 
00366 
00367   if (p_0 < 0.0)
00368     p_0 = TOL; 
00369 
00370   return p_0; 
00371 }  
00372 /*------------------------------------------------------------------*/
00373 
00374 /*! \fn static Real PFunc(const Prim1DS W, Real POld)
00375  *  \brief Evaluates the pressure function (see Toro) for the exact 
00376    riemann solver, given the pressure POld */
00377 static Real PFunc(const Prim1DS W, Real POld)
00378 {
00379   Real f; 
00380   Real a = sqrt((double) Gamma*W.P/W.d);
00381   Real Ak, Bk; 
00382 
00383   /* rarefaction wave */
00384   if (POld <= W.P)  
00385     f = 2*a/(Gamma-1) * (pow((double) POld/W.P, (double) 
00386                              (Gamma-1)/(2*Gamma)) - 1); 
00387 
00388   /* shock wave */
00389   else {
00390     Ak = 2/(W.d*(Gamma+1)); 
00391     Bk = W.P*(Gamma-1)/(Gamma+1); 
00392 
00393     f = (POld - W.P) * sqrt((double) Ak/(POld + Bk)); 
00394   }
00395 
00396   return f; 
00397 }
00398    
00399 /*------------------------------------------------------------------*/
00400 
00401 /*! \fn static Real PFuncDeriv(const Prim1DS W, Real POld)
00402  *  \brief  evaluate the derivative of the pressure function (see Toro) 
00403  * at the pressure value POld */
00404 static Real PFuncDeriv(const Prim1DS W, Real POld)
00405 {
00406   Real fDer; 
00407   Real a = sqrt((double) Gamma*W.P/W.d);
00408   Real Ak, Bk; 
00409 
00410   /* rarefaction wave */
00411   if (POld <= W.P)  
00412     fDer = 1/(a * W.d) * pow((double) POld/W.P, (double) -
00413                              (Gamma+1)/(2*Gamma)); 
00414 
00415   /* shock wave */
00416   else {
00417     Ak = 2/(W.d*(Gamma+1)); 
00418     Bk = W.P*(Gamma-1)/(Gamma+1); 
00419 
00420     fDer = sqrt((double) Ak/(POld + Bk))
00421       *(1.0 - 0.5*(POld - W.P)/(Bk + POld)); 
00422   }
00423 
00424   return fDer; 
00425 }
00426    
00427 /*------------------------------------------------------------------*/ 
00428     
00429 /*! \fn static Real getPC(const Prim1DS Wl, const Prim1DS Wr)
00430  *  \brief Uses Newton-Raphson iteration to find the alebraic root of the 
00431  * pressure function and determine the pressure solution in the 
00432  * center region, Pc. Fails if iteration diverges */
00433 static Real getPC(const Prim1DS Wl, const Prim1DS Wr)
00434 {
00435   Real POld; 
00436   Real VxDiff = Wr.Vx - Wl.Vx;
00437   int i = 0;   
00438   Real TOL = 1.0e-6;
00439   Real change, p, fr, fl, frder, flder; 
00440   Real MAX_ITER = 100; 
00441 
00442   POld = guessP(Wl, Wr);
00443   
00444   while (i < MAX_ITER) {
00445 
00446     fl = PFunc(Wl, POld); 
00447     fr = PFunc(Wr, POld); 
00448     flder = PFuncDeriv(Wl, POld); 
00449     frder = PFuncDeriv(Wr, POld); 
00450     
00451     p = POld - (fl + fr + VxDiff)/(flder + frder);
00452 
00453     change = 2.0*fabs((p-POld)/(p+POld)); 
00454     if (change <= TOL)
00455       return p;
00456     
00457     if (p < 0.0)
00458       p = TOL; 
00459 
00460     POld = p; 
00461     i++; 
00462   }
00463    
00464   /* exit failure due to divergence */
00465   ath_error("[exact_flux]: Divergence in Newton-Raphson \
00466 iteration: p = %e\n", p); 
00467   
00468 }
00469 /*------------------------------------------------------------------*/
00470 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur, 
00471  *          const Prim1DS Wl, const Prim1DS Wr, 
00472  *          const Real Bxi, Cons1DS *pF)
00473  *  \brief Computes fluxes of CONSERVES variables
00474  *
00475  *   Input Arguments:
00476  *   - Bxi = B in direction of slice at cell interface
00477  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00478  *   Output Arguments:
00479  *   - pF = pointer to fluxes of CONSERVED variables at cell interface
00480  */
00481 
00482 void fluxes(const Cons1DS Ul, const Cons1DS Ur, 
00483             const Prim1DS Wl, const Prim1DS Wr, 
00484             const Real Bxi, Cons1DS *pF)
00485 {
00486   Real pc, fr, fl; 
00487   Real dc, dcl, dcr, Vxc;
00488   Real sl, sr;    /* Left and right going shock velocity */
00489   Real hdl, hdr;  /* Left and right going rarefaction head velocity */
00490   Real tll, tlr;  /* Left and right going rarefaction tail velocity */
00491   Real al = sqrt((double) Gamma*Wl.P/Wl.d); /* left sound speed */
00492   Real ar = sqrt((double) Gamma*Wr.P/Wr.d); /* right sound speed */
00493   Real tmp1, tmp2;
00494   Real e, V, E; 
00495  
00496   if(!(Ul.d > 0.0)||!(Ur.d > 0.0))
00497     ath_error("[exact flux]: Non-positive densities: dl = %e  dr = %e\n", 
00498               Ul.d, Ur.d);
00499 
00500   pc = getPC(Wl, Wr); 
00501   
00502   /*----------------------------------------------------------------*/
00503   /* calculate Vxc */
00504   fr = PFunc(Wr, pc); 
00505   fl = PFunc(Wl, pc); 
00506   Vxc = 0.5*(Wl.Vx + Wr.Vx) + 0.5*(fr - fl); 
00507 
00508   /*-----------------------------------------------------------------*/
00509   /* calucate density to left of contact (dcl) */
00510   if (pc > Wl.P) {
00511 
00512     /* left shock wave */
00513     Real tmp = (Gamma - 1.0) / (Gamma + 1.0);
00514     dcl = Wl.d*(pc/Wl.P + tmp)/(tmp*pc/Wl.P + 1); 
00515   }
00516   else {
00517 
00518     /* left rarefaction wave */
00519     dcl = Wl.d*pow((double) (pc/Wl.P), (double) (1/Gamma)); 
00520   }
00521 
00522   if (dcl < 0.0)
00523      ath_error("[exact flux]: Solver finds negative density %5.4e\n", dcl);
00524   
00525   /*-----------------------------------------------------------------*/
00526   /* calculate density to the right of contact (dcr) */
00527   if (pc > Wr.P) {
00528 
00529     /* right shock wave */
00530     Real tmp = (Gamma - 1)/(Gamma + 1);
00531     dcr = Wr.d*(pc/Wr.P + tmp)/(tmp*pc/Wr.P + 1); 
00532   }
00533   else {
00534 
00535     /* right rarefaction wave */
00536     dcr = Wr.d*pow((double) (pc/Wr.P), (double) (1/Gamma)); 
00537   }
00538 
00539   if (dcr < 0.0)
00540     ath_error("[exact flux]: Solver finds negative density %5.4e\n", dcr);
00541  /*-----------------------------------------------------------------
00542   * Calculate the Interface Flux if the wave speeds are such that we aren't
00543   * actually in the intermediate state */
00544   
00545   if (pc > Wl.P) {
00546     /* left shock wave */
00547     
00548     /* left shock speed */
00549     sl = Wl.Vx - al*sqrt((double)(pc*(Gamma+1)/(2*Gamma*Wl.P) + 
00550                                   (Gamma-1)/(2*Gamma))); 
00551     if (sl >= 0.0) {
00552       /* to left of shock */
00553      
00554       e = Wl.P/(Wl.d*(Gamma-1));
00555       V = Wl.Vx*Wl.Vx + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00556       E = Wl.d*(0.5*V + e); 
00557       
00558       pF->E = Wl.Vx*(E + Wl.P);
00559       pF->d  = Ul.Mx;
00560       pF->Mx = Ul.Mx*(Wl.Vx) + Wl.P;
00561       pF->My = Ul.My*(Wl.Vx);
00562       pF->Mz = Ul.Mz*(Wl.Vx);
00563 
00564       return; 
00565     }
00566   }
00567   else {
00568     /* left rarefaction */
00569     
00570     Real alc = al*pow((double)(pc/Wl.P), (double)(Gamma-1)/(2*Gamma));
00571     
00572     hdl = Wl.Vx - al; 
00573     tll = Vxc - alc; 
00574     
00575     if (hdl >= 0.0) {
00576       /* To left of rarefaction */
00577         
00578       e = Wl.P/(Wl.d*(Gamma-1));
00579       V = Wl.Vx*Wl.Vx + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00580       E = Wl.d*(0.5*V + e); 
00581       
00582       pF->E =  Wl.Vx*(E + Wl.P);
00583       pF->d  = Ul.Mx;
00584       pF->Mx = Ul.Mx*(Wl.Vx) + Wl.P;
00585       pF->My = Ul.My*(Wl.Vx);
00586       pF->Mz = Ul.Mz*(Wl.Vx);
00587       return;
00588       
00589     } 
00590     else if (tll >= 0.0) {
00591       /* Inside rarefaction fan */
00592        
00593       
00594       tmp1 = 2/(Gamma + 1); 
00595       tmp2 = (Gamma - 1)/(al*(Gamma+1)); 
00596       
00597       dc = Wl.d*pow((double)(tmp1 + tmp2*Wl.Vx), (double)(2/(Gamma-1))); 
00598       Vxc = tmp1*(al + Wl.Vx*(Gamma-1)/2);
00599       pc = Wl.P*pow((double)(tmp1+tmp2*Wl.Vx),(double)(2*Gamma/(Gamma-1))); 
00600       
00601       e = pc/(dc*(Gamma-1));
00602       V = Vxc*Vxc + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00603       E = dc*(0.5*V + e);       
00604       
00605       pF->E = Vxc*(E + pc); 
00606       pF->d  = dc*Vxc;
00607       pF->Mx = dc*Vxc*Vxc + pc;
00608       pF->My = dc*Vxc*Wl.Vy;
00609       pF->Mz = dc*Vxc*Wl.Vz;
00610       return;
00611     } 
00612   }
00613 
00614   if (pc > Wr.P) {
00615     /* right shock wave */
00616     
00617     /* right shock speed */
00618     sr = Wr.Vx + ar*sqrt((double)(pc*(Gamma+1)/(2*Gamma*Wr.P) + 
00619                                   (Gamma-1)/(2*Gamma))); 
00620     if (sr <= 0.0) {
00621       /* to right of shock */
00622 
00623       e = Wr.P/(Wr.d*(Gamma-1));
00624       V = Wr.Vx*Wr.Vx + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00625       E = Wr.d*(0.5*V + e);     
00626       
00627       pF->E = Wr.Vx*(E + Wr.P);
00628       pF->d  = Ur.Mx;
00629       pF->Mx = Ur.Mx*(Wr.Vx) + Wr.P;
00630       pF->My = Ur.My*(Wr.Vx);
00631       pF->Mz = Ur.Mz*(Wr.Vx);
00632       return; 
00633     }
00634   }
00635   else {
00636     /* right rarefaction */
00637     
00638     Real arc = ar*pow((double)(pc/Wr.P), (double)(Gamma-1)/(2*Gamma)); 
00639       
00640     hdr = Wr.Vx + ar; 
00641     tlr = Vxc + arc; 
00642     
00643     if (hdr <= 0.0) {
00644       /* To right of rarefaction */
00645 
00646       e = Wr.P/(Wr.d*(Gamma-1));
00647       V = Wr.Vx*Wr.Vx + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00648       E = Wr.d*(0.5*V + e);     
00649       
00650       pF->E = Wr.Vx*(E + Wr.P);
00651       pF->d  = Ur.Mx;
00652       pF->Mx = Ur.Mx*(Wr.Vx) + Wr.P;
00653       pF->My = Ur.My*(Wr.Vx);
00654       pF->Mz = Ur.Mz*(Wr.Vx);
00655       return;
00656       
00657     } else if (tlr <= 0.0) {
00658       /* Inside rarefaction fan */
00659       
00660       tmp1 = 2/(Gamma + 1); 
00661       tmp2 = (Gamma - 1)/(ar*(Gamma+1)); 
00662       
00663       dc = Wr.d*pow((double)(tmp1 - tmp2*Wr.Vx), (double)(2/(Gamma-1))); 
00664       Vxc = tmp1*(-ar + Wr.Vx*(Gamma-1)/2);
00665       pc = Wr.P*pow((double)(tmp1-tmp2*Wr.Vx), (double)(2*Gamma/(Gamma-1))); 
00666       
00667       e = pc/(dc*(Gamma-1));
00668       V = Vxc*Vxc + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00669       E = dc*(0.5*V + e); 
00670       
00671       pF->E = Vxc*(E + pc);
00672       pF->d  = dc*Vxc;
00673       pF->Mx = dc*Vxc*Vxc + pc;
00674       pF->My = dc*Vxc*Wr.Vy;
00675       pF->Mz = dc*Vxc*Wr.Vz;
00676       return;
00677     }
00678   }
00679     
00680 /* We are in the intermediate state */
00681 
00682 /*---------------------------------------------------------------------
00683  * Calculate the Interface Flux */
00684   if (Vxc >= 0.0) {
00685 
00686     e = pc/(dcl*(Gamma-1));
00687     V = Vxc*Vxc + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00688     E = dcl*(0.5*V + e); 
00689     
00690     pF->E = Vxc*(E + pc); 
00691     pF->d  = dcl*Vxc;
00692     pF->Mx = dcl*Vxc*Vxc + pc;
00693     pF->My = dcl*Vxc*Wl.Vy;
00694     pF->Mz = dcl*Vxc*Wl.Vz;
00695   }
00696   else {
00697 
00698     e = pc/(dcr*(Gamma-1));
00699     V = Vxc*Vxc + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00700     E = dcr*(0.5*V + e); 
00701     
00702     pF->E = Vxc*(E + pc); 
00703     pF->d  = dcr*Vxc;
00704     pF->Mx = dcr*Vxc*Vxc + pc;
00705     pF->My = dcr*Vxc*Wr.Vy;
00706     pF->Mz = dcr*Vxc*Wr.Vz;
00707   }
00708   
00709   return;
00710 }
00711 /*------------------------------------------------------------------*/
00712 
00713 #endif /*ISOTHERMAL*/
00714 
00715 #endif /* SPECIAL_RELATIVITY */
00716 #endif /* EXACT_FLUX */