Athena: rsolvers/hlld.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file hlld.c
00004  *  \brief Computes 1D fluxes using the HLLD Riemann solver.
00005  *
00006  * PURPOSE: Computes 1D fluxes using the HLLD Riemann solver, an extension of
00007  *   the HLLE solver to MHD.  Only works for MHD problems.  SEPARATE code
00008  *   blocks for adiabatic and isothermal equations of state.
00009  *
00010  * REFERENCES:
00011  * - T. Miyoshi & K. Kusano, "A multi-state HLL approximate Riemann solver
00012  *   for ideal MHD", JCP, 208, 315 (2005)
00013  * - A. Mignone, "A simple and accurate Riemann solver for isothermal MHD",
00014  *   JPC, 225, 1427 (2007)
00015  *
00016  * HISTORY: Adiabatic version written by Brian Biskeborn, May 8, 2006,
00017  *             COS sophmore project.
00018  *          Isothermal version written by Nicole Lemaster, May 1, 2008.
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 "../defs.h"
00029 #include "../athena.h"
00030 #include "../globals.h"
00031 #include "prototypes.h"
00032 #include "../prototypes.h"
00033 
00034 #ifdef HLLD_FLUX
00035 #ifndef SPECIAL_RELATIVITY
00036 
00037 #define SMALL_NUMBER 1e-8
00038 
00039 #ifndef MHD
00040 #error : The HLLD flux only works for mhd.
00041 #endif /* MHD */
00042 
00043 #ifndef ISOTHERMAL
00044 
00045 /*----------------------------------------------------------------------------*/
00046 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00047  *           const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00048  *  \brief Compute 1D fluxes
00049  * Input Arguments:
00050  * - Bxi = B in direction of slice at cell interface
00051  * - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00052  *
00053  * Output Arguments:
00054  * - Flux = fluxes of CONSERVED variables at cell interface
00055  */
00056 
00057 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00058             const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00059 {
00060   Cons1DS Ulst,Uldst,Urdst,Urst;       /* Conserved variable for all states */
00061   Prim1DS Wlst,Wrst;                   /* Primitive variables for all states */
00062   Cons1DS Fl,Fr;                       /* Fluxes for left & right states */
00063   Real spd[5];                        /* signal speeds, left to right */
00064 /*  Real maxspd; */
00065   Real sdl,sdr,sdml,sdmr;             /* S_i-u_i, S_i-S_M (i=L or R) */
00066   Real pbl,pbr;                       /* Magnetic pressures */
00067   Real cfl,cfr,cfmax;                 /* Cf (left & right), max(cfl,cfr) */
00068   Real gpl,gpr,gpbl,gpbr;             /* gamma*P, gamma*P + B */
00069   Real sqrtdl,sqrtdr;                 /* sqrt of the L* & R* densities */
00070   Real invsumd;                       /* 1/(sqrtdl + sqrtdr) */
00071   Real ptl,ptr,ptst;                  /* total pressures */
00072   Real vbstl,vbstr;                   /* v_i* dot B_i* for i=L or R */
00073   Real Bxsig;                         /* sign(Bx) = 1 for Bx>0, -1 for Bx<0 */
00074   Real Bxsq = SQR(Bxi);               /* Bx^2 */
00075   Real tmp;                      /* Temp variable for repeated calculations */
00076 #if (NSCALARS > 0)
00077   int n;
00078 #endif
00079 
00080 
00081 
00082 /*--- Step 1. ------------------------------------------------------------------
00083  * Convert left- and right- states in conserved to primitive variables.
00084  */
00085 
00086 /*
00087   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00088   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00089 */
00090 
00091 /*--- Step 2. ------------------------------------------------------------------
00092  * Compute left & right wave speeds according to Miyoshi & Kusano, eqn. (67)
00093  */
00094 
00095   pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00096   pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00097   gpl  = Gamma * Wl.P;
00098   gpr  = Gamma * Wr.P;
00099   gpbl = gpl + 2.0*pbl;
00100   gpbr = gpr + 2.0*pbr;
00101 
00102   cfl = sqrt((gpbl + sqrt(SQR(gpbl)-4*gpl*Bxsq))/(2.0*Wl.d));
00103   cfr = sqrt((gpbr + sqrt(SQR(gpbr)-4*gpr*Bxsq))/(2.0*Wr.d));
00104   cfmax = MAX(cfl,cfr);
00105 
00106   if(Wl.Vx <= Wr.Vx) {
00107     spd[0] = Wl.Vx - cfmax;
00108     spd[4] = Wr.Vx + cfmax;
00109   }
00110   else {
00111     spd[0] = Wr.Vx - cfmax;
00112     spd[4] = Wl.Vx + cfmax;
00113   }
00114 
00115 /*  maxspd = MAX(fabs(spd[0]),fabs(spd[4])); */
00116 
00117 /*--- Step 3. ------------------------------------------------------------------
00118  * Compute L/R fluxes
00119  */
00120 
00121   /* total pressure */
00122   ptl = Wl.P + pbl;
00123   ptr = Wr.P + pbr;
00124 
00125   Fl.d  = Ul.Mx;
00126   Fl.Mx = Ul.Mx*Wl.Vx + ptl - Bxsq;
00127   Fl.My = Ul.d*Wl.Vx*Wl.Vy - Bxi*Ul.By;
00128   Fl.Mz = Ul.d*Wl.Vx*Wl.Vz - Bxi*Ul.Bz;
00129   Fl.E  = Wl.Vx*(Ul.E + ptl - Bxsq) - Bxi*(Wl.Vy*Ul.By + Wl.Vz*Ul.Bz);
00130   Fl.By = Ul.By*Wl.Vx - Bxi*Wl.Vy;
00131   Fl.Bz = Ul.Bz*Wl.Vx - Bxi*Wl.Vz;
00132 
00133   Fr.d  = Ur.Mx;
00134   Fr.Mx = Ur.Mx*Wr.Vx + ptr - Bxsq;
00135   Fr.My = Ur.d*Wr.Vx*Wr.Vy - Bxi*Ur.By;
00136   Fr.Mz = Ur.d*Wr.Vx*Wr.Vz - Bxi*Ur.Bz;
00137   Fr.E  = Wr.Vx*(Ur.E + ptr - Bxsq) - Bxi*(Wr.Vy*Ur.By + Wr.Vz*Ur.Bz);
00138   Fr.By = Ur.By*Wr.Vx - Bxi*Wr.Vy;
00139   Fr.Bz = Ur.Bz*Wr.Vx - Bxi*Wr.Vz;
00140 
00141 #if (NSCALARS > 0)
00142   for (n=0; n<NSCALARS; n++) {
00143     Fl.s[n] = Fl.d*Wl.r[n];
00144     Fr.s[n] = Fr.d*Wr.r[n];
00145   }
00146 #endif
00147 
00148 /*--- Step 4. ------------------------------------------------------------------
00149  * Return upwind flux if flow is supersonic
00150  */
00151 
00152   if(spd[0] >= 0.0){
00153     *pFlux = Fl;
00154 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00155     pFlux->Pflux = ptl;
00156 #endif 
00157     return;
00158   }
00159 
00160   if(spd[4] <= 0.0){
00161     *pFlux = Fr;
00162 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00163     pFlux->Pflux = ptr;
00164 #endif 
00165     return;
00166   }
00167 
00168 /*--- Step 5. ------------------------------------------------------------------
00169  * Compute middle and Alfven wave speeds
00170  */
00171 
00172   sdl = spd[0] - Wl.Vx;
00173   sdr = spd[4] - Wr.Vx;
00174 
00175   /* S_M: eqn (38) of Miyoshi & Kusano */
00176   spd[2] = (sdr*Wr.d*Wr.Vx - sdl*Wl.d*Wl.Vx - ptr + ptl) /
00177            (sdr*Wr.d-sdl*Wl.d);
00178 
00179   sdml   = spd[0] - spd[2];
00180   sdmr   = spd[4] - spd[2];
00181   /* eqn (43) of Miyoshi & Kusano */
00182   Ulst.d = Ul.d * sdl/sdml;
00183   Urst.d = Ur.d * sdr/sdmr;
00184   sqrtdl = sqrt(Ulst.d);
00185   sqrtdr = sqrt(Urst.d);
00186 
00187   /* eqn (51) of Miyoshi & Kusano */
00188   spd[1] = spd[2] - fabs(Bxi)/sqrtdl;
00189   spd[3] = spd[2] + fabs(Bxi)/sqrtdr;
00190 
00191 /*--- Step 6. ------------------------------------------------------------------
00192  * Compute intermediate states
00193  */
00194 
00195   ptst = ptl + Ul.d*sdl*(sdl-sdml);
00196  
00197 /* Ul* */
00198   /* eqn (39) of M&K */
00199   Ulst.Mx = Ulst.d * spd[2];
00200 //   if((fabs(spd[2]/Wl.Vx-1.0)<SMALL_NUMBER) ||
00201 //      (fabs(spd[2])/fabs(spd[0]) <= SMALL_NUMBER &&
00202 //       fabs(Wl.Vx)/fabs(spd[0]) <= SMALL_NUMBER)) {
00203 //     Ulst.My = Ulst.d * Wl.Vy;
00204 //     Ulst.Mz = Ulst.d * Wl.Vz;
00205 // 
00206 //     Ulst.By = Ul.By;
00207 //     Ulst.Bz = Ul.Bz;
00208 //   }
00209   if (fabs(Ul.d*sdl*sdml/Bxsq-1.0) < SMALL_NUMBER) {
00210     /* Degenerate case */
00211     Ulst.My = Ulst.d * Wl.Vy;
00212     Ulst.Mz = Ulst.d * Wl.Vz;
00213 
00214     Ulst.By = Ul.By;
00215     Ulst.Bz = Ul.Bz;
00216   }
00217   else {
00218     /* eqns (44) and (46) of M&K */
00219     tmp = Bxi*(sdl-sdml)/(Ul.d*sdl*sdml-Bxsq);
00220     Ulst.My = Ulst.d * (Wl.Vy - Ul.By*tmp);
00221     Ulst.Mz = Ulst.d * (Wl.Vz - Ul.Bz*tmp);
00222 //     if(Ul.By == 0.0 && Ul.Bz == 0.0) {
00223 //       Ulst.By = 0.0;
00224 //       Ulst.Bz = 0.0;
00225 //     }
00226 //     else {
00227 //       /* eqns (45) and (47) of M&K */
00228 //       tmp = (Ul.d*SQR(sdl)-Bxsq)/(Ul.d*sdl*sdml - Bxsq);
00229 //       Ulst.By = Ul.By * tmp;
00230 //       Ulst.Bz = Ul.Bz * tmp;
00231 //     }
00232 
00233     /* eqns (45) and (47) of M&K */
00234     tmp = (Ul.d*SQR(sdl)-Bxsq)/(Ul.d*sdl*sdml - Bxsq);
00235     Ulst.By = Ul.By * tmp;
00236     Ulst.Bz = Ul.Bz * tmp;
00237   }
00238   vbstl = (Ulst.Mx*Bxi+Ulst.My*Ulst.By+Ulst.Mz*Ulst.Bz)/Ulst.d;
00239   /* eqn (48) of M&K */
00240   Ulst.E = (sdl*Ul.E - ptl*Wl.Vx + ptst*spd[2] +
00241             Bxi*(Wl.Vx*Bxi+Wl.Vy*Ul.By+Wl.Vz*Ul.Bz - vbstl))/sdml;
00242   Wlst = Cons1D_to_Prim1D(&Ulst,&Bxi);
00243 
00244 
00245 /* Ur* */
00246   /* eqn (39) of M&K */
00247   Urst.Mx = Urst.d * spd[2];
00248 //   if((fabs(spd[2]/Wr.Vx-1.0)<SMALL_NUMBER) ||
00249 //      (fabs(spd[2])/fabs(spd[4]) <= SMALL_NUMBER &&
00250 //       fabs(Wr.Vx)/fabs(spd[4]) <= SMALL_NUMBER)) {
00251 //     Urst.My = Urst.d * Wr.Vy;
00252 //     Urst.Mz = Urst.d * Wr.Vz;
00253 // 
00254 //     Urst.By = Ur.By;
00255 //     Urst.Bz = Ur.Bz;
00256 //   }
00257   if (fabs(Ur.d*sdr*sdmr/Bxsq-1.0) < SMALL_NUMBER) {
00258     /* Degenerate case */
00259     Urst.My = Urst.d * Wr.Vy;
00260     Urst.Mz = Urst.d * Wr.Vz;
00261 
00262     Urst.By = Ur.By;
00263     Urst.Bz = Ur.Bz;
00264   }
00265   else {
00266     /* eqns (44) and (46) of M&K */
00267     tmp = Bxi*(sdr-sdmr)/(Ur.d*sdr*sdmr-Bxsq);
00268     Urst.My = Urst.d * (Wr.Vy - Ur.By*tmp);
00269     Urst.Mz = Urst.d * (Wr.Vz - Ur.Bz*tmp);
00270 
00271 //     if(Ur.By == 0.0 && Ur.Bz == 0.0) {
00272 //       Urst.By = 0.0;
00273 //       Urst.Bz = 0.0;
00274 //     }
00275 //     else {
00276 //       /* eqns (45) and (47) of M&K */
00277 //       tmp = (Ur.d*SQR(sdr)-Bxsq)/(Ur.d*sdr*sdmr - Bxsq);
00278 //       Urst.By = Ur.By * tmp;
00279 //       Urst.Bz = Ur.Bz * tmp;
00280 //     }
00281 
00282     /* eqns (45) and (47) of M&K */
00283     tmp = (Ur.d*SQR(sdr)-Bxsq)/(Ur.d*sdr*sdmr - Bxsq);
00284     Urst.By = Ur.By * tmp;
00285     Urst.Bz = Ur.Bz * tmp;
00286   }
00287   vbstr = (Urst.Mx*Bxi+Urst.My*Urst.By+Urst.Mz*Urst.Bz)/Urst.d;
00288   /* eqn (48) of M&K */
00289   Urst.E = (sdr*Ur.E - ptr*Wr.Vx + ptst*spd[2] +
00290             Bxi*(Wr.Vx*Bxi+Wr.Vy*Ur.By+Wr.Vz*Ur.Bz - vbstr))/sdmr;
00291   Wrst = Cons1D_to_Prim1D(&Urst,&Bxi);
00292 
00293 
00294 /* Ul** and Ur** - if Bx is zero, same as *-states */
00295 //   if(Bxi == 0.0) {
00296   if(0.5*Bxsq/MIN(pbl,pbr) < SQR(SMALL_NUMBER)) {
00297     Uldst = Ulst;
00298     Urdst = Urst;
00299   }
00300   else {
00301     invsumd = 1.0/(sqrtdl + sqrtdr);
00302     if(Bxi > 0) Bxsig =  1;
00303     else        Bxsig = -1;
00304 
00305     Uldst.d = Ulst.d;
00306     Urdst.d = Urst.d;
00307 
00308     Uldst.Mx = Ulst.Mx;
00309     Urdst.Mx = Urst.Mx;
00310 
00311     /* eqn (59) of M&K */
00312     tmp = invsumd*(sqrtdl*Wlst.Vy + sqrtdr*Wrst.Vy + Bxsig*(Urst.By-Ulst.By));
00313     Uldst.My = Uldst.d * tmp;
00314     Urdst.My = Urdst.d * tmp;
00315 
00316     /* eqn (60) of M&K */
00317     tmp = invsumd*(sqrtdl*Wlst.Vz + sqrtdr*Wrst.Vz + Bxsig*(Urst.Bz-Ulst.Bz));
00318     Uldst.Mz = Uldst.d * tmp;
00319     Urdst.Mz = Urdst.d * tmp;
00320 
00321     /* eqn (61) of M&K */
00322     tmp = invsumd*(sqrtdl*Urst.By + sqrtdr*Ulst.By +
00323                    Bxsig*sqrtdl*sqrtdr*(Wrst.Vy-Wlst.Vy));
00324     Uldst.By = Urdst.By = tmp;
00325 
00326     /* eqn (62) of M&K */
00327     tmp = invsumd*(sqrtdl*Urst.Bz + sqrtdr*Ulst.Bz +
00328                    Bxsig*sqrtdl*sqrtdr*(Wrst.Vz-Wlst.Vz));
00329     Uldst.Bz = Urdst.Bz = tmp;
00330 
00331     /* eqn (63) of M&K */
00332     tmp = spd[2]*Bxi + (Uldst.My*Uldst.By + Uldst.Mz*Uldst.Bz)/Uldst.d;
00333     Uldst.E = Ulst.E - sqrtdl*Bxsig*(vbstl - tmp);
00334     Urdst.E = Urst.E + sqrtdr*Bxsig*(vbstr - tmp);
00335   }
00336 
00337 /*--- Step 7. ------------------------------------------------------------------
00338  * Compute flux
00339  */
00340 
00341   if(spd[1] >= 0) {
00342 /* return Fl* */
00343     pFlux->d  = Fl.d  + spd[0]*(Ulst.d  - Ul.d);
00344     pFlux->Mx = Fl.Mx + spd[0]*(Ulst.Mx - Ul.Mx);
00345     pFlux->My = Fl.My + spd[0]*(Ulst.My - Ul.My);
00346     pFlux->Mz = Fl.Mz + spd[0]*(Ulst.Mz - Ul.Mz);
00347     pFlux->E  = Fl.E  + spd[0]*(Ulst.E  - Ul.E);
00348     pFlux->By = Fl.By + spd[0]*(Ulst.By - Ul.By);
00349     pFlux->Bz = Fl.Bz + spd[0]*(Ulst.Bz - Ul.Bz);
00350   }
00351   else if(spd[2] >= 0) {
00352 /* return Fl** */
00353     tmp = spd[1] - spd[0];
00354     pFlux->d  = Fl.d  - spd[0]*Ul.d  - tmp*Ulst.d  + spd[1]*Uldst.d;
00355     pFlux->Mx = Fl.Mx - spd[0]*Ul.Mx - tmp*Ulst.Mx + spd[1]*Uldst.Mx;
00356     pFlux->My = Fl.My - spd[0]*Ul.My - tmp*Ulst.My + spd[1]*Uldst.My;
00357     pFlux->Mz = Fl.Mz - spd[0]*Ul.Mz - tmp*Ulst.Mz + spd[1]*Uldst.Mz;
00358     pFlux->E  = Fl.E  - spd[0]*Ul.E  - tmp*Ulst.E  + spd[1]*Uldst.E;
00359     pFlux->By = Fl.By - spd[0]*Ul.By - tmp*Ulst.By + spd[1]*Uldst.By;
00360     pFlux->Bz = Fl.Bz - spd[0]*Ul.Bz - tmp*Ulst.Bz + spd[1]*Uldst.Bz;
00361   }
00362   else if(spd[3] > 0) {
00363 /* return Fr** */
00364     tmp = spd[3] - spd[4];
00365     pFlux->d  = Fr.d  - spd[4]*Ur.d  - tmp*Urst.d  + spd[3]*Urdst.d;
00366     pFlux->Mx = Fr.Mx - spd[4]*Ur.Mx - tmp*Urst.Mx + spd[3]*Urdst.Mx;
00367     pFlux->My = Fr.My - spd[4]*Ur.My - tmp*Urst.My + spd[3]*Urdst.My;
00368     pFlux->Mz = Fr.Mz - spd[4]*Ur.Mz - tmp*Urst.Mz + spd[3]*Urdst.Mz;
00369     pFlux->E  = Fr.E  - spd[4]*Ur.E  - tmp*Urst.E  + spd[3]*Urdst.E;
00370     pFlux->By = Fr.By - spd[4]*Ur.By - tmp*Urst.By + spd[3]*Urdst.By;
00371     pFlux->Bz = Fr.Bz - spd[4]*Ur.Bz - tmp*Urst.Bz + spd[3]*Urdst.Bz;
00372   }
00373   else {
00374 /* return Fr* */
00375     pFlux->d  = Fr.d  + spd[4]*(Urst.d  - Ur.d);
00376     pFlux->Mx = Fr.Mx + spd[4]*(Urst.Mx - Ur.Mx);
00377     pFlux->My = Fr.My + spd[4]*(Urst.My - Ur.My);
00378     pFlux->Mz = Fr.Mz + spd[4]*(Urst.Mz - Ur.Mz);
00379     pFlux->E  = Fr.E  + spd[4]*(Urst.E  - Ur.E);
00380     pFlux->By = Fr.By + spd[4]*(Urst.By - Ur.By);
00381     pFlux->Bz = Fr.Bz + spd[4]*(Urst.Bz - Ur.Bz);
00382   }
00383 
00384 /* Fluxes of passively advected scalars, computed from density flux */
00385 #if (NSCALARS > 0)
00386   if (pFlux->d >= 0.0) {
00387     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00388   } else {
00389     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00390   }
00391 #endif
00392 
00393 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00394   pFlux->Pflux = ptst;
00395 #endif
00396   return;
00397 }
00398 
00399 #else /* ISOTHERMAL */
00400 
00401 /*----------------------------------------------------------------------------*/
00402 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00403  *          const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00404  *  \brief Compute 1D fluxes
00405  * Input Arguments:
00406  * - Bxi = B in direction of slice at cell interface
00407  * - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00408  *
00409  * Output Arguments:
00410  * - Flux = fluxes of CONSERVED variables at cell interface
00411  */
00412 
00413 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00414             const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00415 {
00416   Cons1DS Ulst,Ucst,Urst;              /* Conserved variable for all states */
00417   Cons1DS Fl,Fr;                       /* Fluxes for left & right states */
00418   Real spd[5],idspd;                  /* signal speeds, left to right */
00419   Real pbl,pbr;                       /* Magnetic pressures */
00420   Real cfl,cfr;                       /* Cf (left & right) */
00421   Real gpl,gpr,gpbl,gpbr;             /* gamma*P, gamma*P + B */
00422   Real mxhll,ustar,dhll,sqrtdhll;     /* Mx, vel, and density in star states */
00423   Real fdhll,fmxhll;                  /* HLL fluxes (for star states) */
00424   Real ptl,ptr;                       /* total pressures */
00425   Real Bxsq = SQR(Bxi);               /* Bx^2 */
00426   Real tmp,mfact,bfact,X;             /* temporary variables */
00427   int n;
00428 
00429 //   /* FROM ROE */
00430 //   Real sqrtdl,sqrtdr,isdlpdr,idroe,v1roe,b2roe,b3roe,x,y;
00431 //   Real bt_starsq,vaxsq,twid_csq,ct2,tsum,tdif,cf2_cs2,cfsq,cf;
00432 
00433 /*--- Step 1. ------------------------------------------------------------------
00434  * Convert left- and right- states in conserved to primitive variables.
00435  */
00436 
00437 /*
00438   pbl = Cons1D_to_Prim1D(&Ul,&Wl,&Bxi);
00439   pbr = Cons1D_to_Prim1D(&Ur,&Wr,&Bxi);
00440 */
00441 
00442 /*--- Step 2. ------------------------------------------------------------------
00443  * Compute left & right wave speeds according to Mignone, eqns. (39), (8a),
00444  * and (9)
00445  */
00446 
00447   pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00448   pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00449   gpl  = Wl.d*Iso_csound2;
00450   gpr  = Wr.d*Iso_csound2;
00451   gpbl = gpl + 2.0*pbl;
00452   gpbr = gpr + 2.0*pbr;
00453 
00454   cfl = sqrt((gpbl + sqrt(SQR(gpbl)-4*gpl*Bxsq))/(2.0*Wl.d));
00455   cfr = sqrt((gpbr + sqrt(SQR(gpbr)-4*gpr*Bxsq))/(2.0*Wr.d));
00456 
00457   spd[0] = MIN(Wl.Vx-cfl,Wr.Vx-cfr);
00458   spd[4] = MAX(Wl.Vx+cfl,Wr.Vx+cfr);
00459 
00460 //   /* COMPUTE ROE AVERAGES */
00461 //   sqrtdl = sqrt((double)Wl.d);
00462 //   sqrtdr = sqrt((double)Wr.d);
00463 //   isdlpdr = 1.0/(sqrtdl + sqrtdr);
00464 //   idroe  = 1.0/(sqrtdl*sqrtdr);
00465 //   v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00466 //   b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00467 //   b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00468 //   x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))*SQR(isdlpdr);
00469 //   y = 0.5*(Wl.d + Wr.d)*idroe;
00470 // 
00471 //   /* COMPUTE FAST MAGNETOSONIC SPEED */
00472 //   bt_starsq = (SQR(b2roe) + SQR(b3roe))*y;
00473 //   vaxsq = Bxsq*idroe;
00474 //   twid_csq = Iso_csound2 + x;
00475 //   ct2 = bt_starsq*idroe;
00476 //   tsum = vaxsq + ct2 + twid_csq;
00477 //   cfsq = 0.5*(tsum + sqrt((double)(SQR(tsum) - 4.0*twid_csq*vaxsq)));
00478 //   cf = sqrt((double)cfsq);
00479 // 
00480 //   spd[0] = MIN(Wl.Vx-cfl,v1roe-cf);
00481 //   spd[4] = MAX(v1roe+cf,Wr.Vx+cfr);
00482 
00483 /*--- Step 3. ------------------------------------------------------------------
00484  * Compute L/R fluxes
00485  */
00486 
00487   /* total pressure */
00488   ptl = gpl + pbl;
00489   ptr = gpr + pbr;
00490 
00491   Fl.d  = Ul.Mx;
00492   Fl.Mx = Ul.Mx*Wl.Vx + ptl - Bxsq;
00493   Fl.My = Ul.d*Wl.Vx*Wl.Vy - Bxi*Ul.By;
00494   Fl.Mz = Ul.d*Wl.Vx*Wl.Vz - Bxi*Ul.Bz;
00495   Fl.By = Ul.By*Wl.Vx - Bxi*Wl.Vy;
00496   Fl.Bz = Ul.Bz*Wl.Vx - Bxi*Wl.Vz;
00497 
00498   Fr.d  = Ur.Mx;
00499   Fr.Mx = Ur.Mx*Wr.Vx + ptr - Bxsq;
00500   Fr.My = Ur.d*Wr.Vx*Wr.Vy - Bxi*Ur.By;
00501   Fr.Mz = Ur.d*Wr.Vx*Wr.Vz - Bxi*Ur.Bz;
00502   Fr.By = Ur.By*Wr.Vx - Bxi*Wr.Vy;
00503   Fr.Bz = Ur.Bz*Wr.Vx - Bxi*Wr.Vz;
00504 
00505 #if (NSCALARS > 0)
00506   for (n=0; n<NSCALARS; n++) {
00507     Fl.s[n] = Fl.d*Wl.r[n];
00508     Fr.s[n] = Fr.d*Wr.r[n];
00509   }
00510 #endif
00511 
00512 /*--- Step 4. ------------------------------------------------------------------
00513  * Return upwind flux if flow is supersonic
00514  */
00515 
00516   /* eqn. (38a) of Mignone */
00517   if(spd[0] >= 0.0){
00518     *pFlux = Fl;
00519     return;
00520   }
00521 
00522   /* eqn. (38e) of Mignone */
00523   if(spd[4] <= 0.0){
00524     *pFlux = Fr;
00525     return;
00526   }
00527 
00528 /*--- Step 5. ------------------------------------------------------------------
00529  * Compute hll averages and Alfven wave speed
00530  */
00531 
00532   /* inverse of difference between right and left signal speeds */
00533   idspd = 1.0/(spd[4]-spd[0]);
00534 
00535   /* rho component of U^{hll} from Mignone eqn. (15);
00536    * uses F_L and F_R from eqn. (6) */
00537   dhll = (spd[4]*Ur.d-spd[0]*Ul.d-Fr.d+Fl.d)*idspd;
00538   sqrtdhll = sqrt(dhll);
00539 
00540   /* rho and mx components of F^{hll} from Mignone eqn. (17) */
00541   fdhll = (spd[4]*Fl.d-spd[0]*Fr.d+spd[4]*spd[0]*(Ur.d-Ul.d))*idspd;
00542   fmxhll = (spd[4]*Fl.Mx-spd[0]*Fr.Mx+spd[4]*spd[0]*(Ur.Mx-Ul.Mx))*idspd;
00543 
00544   /* ustar from paragraph between eqns. (23) and (24) */
00545   ustar = fdhll/dhll;
00546 
00547   /* mx component of U^{hll} from Mignone eqn. (15); paragraph referenced
00548    * above states that mxhll should NOT be used to compute ustar */
00549   mxhll = (spd[4]*Ur.Mx-spd[0]*Ul.Mx-Fr.Mx+Fl.Mx)*idspd;
00550 
00551   /* S*_L and S*_R from Mignone eqn. (29) */
00552   spd[1] = ustar - fabs(Bxi)/sqrtdhll;
00553   spd[3] = ustar + fabs(Bxi)/sqrtdhll;
00554 
00555 /*--- Step 6. ------------------------------------------------------------------
00556  * Compute intermediate states
00557  */
00558 
00559 /* Ul* */
00560   /* eqn. (20) of Mignone */
00561   Ulst.d = dhll;
00562   /* eqn. (24) of Mignone */
00563   Ulst.Mx = mxhll;
00564 
00565   tmp = (spd[0]-spd[1])*(spd[0]-spd[3]);
00566 //   if (tmp == 0) {
00567   if ((fabs(spd[0]/spd[1]-1.0) < SMALL_NUMBER) 
00568         || (fabs(spd[0]/spd[3]-1.0) < SMALL_NUMBER)) {
00569     /* degenerate case described below eqn. (39) */
00570     Ulst.My = Ul.My;
00571     Ulst.Mz = Ul.Mz;
00572     Ulst.By = Ul.By;
00573     Ulst.Bz = Ul.Bz;
00574   } else {
00575     mfact = Bxi*(ustar-Wl.Vx)/tmp;
00576     bfact = (Ul.d*SQR(spd[0]-Wl.Vx)-Bxsq)/(dhll*tmp);
00577 
00578     /* eqn. (30) of Mignone */
00579     Ulst.My = dhll*Wl.Vy-Ul.By*mfact;
00580     /* eqn. (31) of Mignone */
00581     Ulst.Mz = dhll*Wl.Vz-Ul.Bz*mfact;
00582     /* eqn. (32) of Mignone */
00583     Ulst.By = Ul.By*bfact;
00584     /* eqn. (33) of Mignone */
00585     Ulst.Bz = Ul.Bz*bfact;
00586   }
00587 
00588 /* Ur* */
00589   /* eqn. (20) of Mignone */
00590   Urst.d = dhll;
00591   /* eqn. (24) of Mignone */
00592   Urst.Mx = mxhll;
00593 
00594   tmp = (spd[4]-spd[1])*(spd[4]-spd[3]);
00595 //   if (tmp == 0) {
00596   if ((fabs(spd[4]/spd[1]-1.0) < SMALL_NUMBER) 
00597         || (fabs(spd[4]/spd[3]-1.0) < SMALL_NUMBER)) {
00598     /* degenerate case described below eqn. (39) */
00599     Urst.My = Ur.My;
00600     Urst.Mz = Ur.Mz;
00601     Urst.By = Ur.By;
00602     Urst.Bz = Ur.Bz;
00603   } else {
00604     mfact = Bxi*(ustar-Wr.Vx)/tmp;
00605     bfact = (Ur.d*SQR(spd[4]-Wr.Vx)-Bxsq)/(dhll*tmp);
00606 
00607     /* eqn. (30) of Mignone */
00608     Urst.My = dhll*Wr.Vy-Ur.By*mfact;
00609     /* eqn. (31) of Mignone */
00610     Urst.Mz = dhll*Wr.Vz-Ur.Bz*mfact;
00611     /* eqn. (32) of Mignone */
00612     Urst.By = Ur.By*bfact;
00613     /* eqn. (33) of Mignone */
00614     Urst.Bz = Ur.Bz*bfact;
00615   }
00616 
00617 /* Uc* */
00618   /* from below eqn. (37) of Mignone */
00619   X = sqrtdhll*SIGN(Bxi);
00620   /* eqn. (20) of Mignone */
00621   Ucst.d = dhll;
00622   /* eqn. (24) of Mignone */
00623   Ucst.Mx = mxhll;
00624   /* eqn. (34) of Mignone */
00625   Ucst.My = 0.5*(Ulst.My+Urst.My+X*(Urst.By-Ulst.By));
00626   /* eqn. (35) of Mignone */
00627   Ucst.Mz = 0.5*(Ulst.Mz+Urst.Mz+X*(Urst.Bz-Ulst.Bz));
00628   /* eqn. (36) of Mignone */
00629   Ucst.By = 0.5*(Ulst.By+Urst.By+(Urst.My-Ulst.My)/X);
00630   /* eqn. (37) of Mignone */
00631   Ucst.Bz = 0.5*(Ulst.Bz+Urst.Bz+(Urst.Mz-Ulst.Mz)/X);
00632 
00633 /*--- Step 7. ------------------------------------------------------------------
00634  * Compute flux
00635  */
00636 
00637   if(spd[1] >= 0) {
00638 /* return (Fl+Sl*(Ulst-Ul)), eqn. (38b) of Mignone */
00639     pFlux->d  = Fl.d  + spd[0]*(Ulst.d  - Ul.d);
00640     pFlux->Mx = Fl.Mx + spd[0]*(Ulst.Mx - Ul.Mx);
00641     pFlux->My = Fl.My + spd[0]*(Ulst.My - Ul.My);
00642     pFlux->Mz = Fl.Mz + spd[0]*(Ulst.Mz - Ul.Mz);
00643     pFlux->By = Fl.By + spd[0]*(Ulst.By - Ul.By);
00644     pFlux->Bz = Fl.Bz + spd[0]*(Ulst.Bz - Ul.Bz);
00645   }
00646   else if (spd[3] <= 0) {
00647 /* return (Fr+Sr*(Urst-Ur)), eqn. (38d) of Mignone */
00648     pFlux->d  = Fr.d  + spd[4]*(Urst.d  - Ur.d);
00649     pFlux->Mx = Fr.Mx + spd[4]*(Urst.Mx - Ur.Mx);
00650     pFlux->My = Fr.My + spd[4]*(Urst.My - Ur.My);
00651     pFlux->Mz = Fr.Mz + spd[4]*(Urst.Mz - Ur.Mz);
00652     pFlux->By = Fr.By + spd[4]*(Urst.By - Ur.By);
00653     pFlux->Bz = Fr.Bz + spd[4]*(Urst.Bz - Ur.Bz);
00654   }
00655   else {
00656 /* return Fcst, eqn. (38c) of Mignone, using eqn. (24) */
00657     pFlux->d = dhll*ustar;
00658     pFlux->Mx = fmxhll;
00659     pFlux->My = Ucst.My*ustar - Bxi*Ucst.By;
00660     pFlux->Mz = Ucst.Mz*ustar - Bxi*Ucst.Bz;
00661     pFlux->By = Ucst.By*ustar - Bxi*Ucst.My/Ucst.d;
00662     pFlux->Bz = Ucst.Bz*ustar - Bxi*Ucst.Mz/Ucst.d;
00663   }
00664 
00665 /* Fluxes of passively advected scalars, computed from density flux */
00666 #if (NSCALARS > 0)
00667   if (pFlux->d >= 0.0) {
00668     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00669   } else {
00670     for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00671   }
00672 #endif
00673 
00674   return;
00675 }
00676 
00677 #endif /* ISOTHERMAL */
00678 #endif /* SPECIAL_RELATIVITY */
00679 #endif /* HLLD_FLUX */