Athena: rsolvers/hllc_sr.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file hllc_sr.c
00004  *  \brief Computes 1D fluxes using the relativistic HLLC Riemann solver. 
00005  *
00006  * PURPOSE: Computes 1D fluxes using the relativistic HLLC Riemann solver, 
00007  *   an extension of the HLLE fluxes to include the contact wave.  Currently 
00008  *   only works for hydrodynamics.  For an extension to MHD, see hlld_sr.c
00009  *
00010  * REFERENCES:
00011  * - A. Mignone and G. Bodo, "An HLLC Riemann solver for relativistic flows",
00012  *   Mon. Not. R. Astron. Soc. 364, 126-136 (2005)
00013  *
00014  * - A. Mignone and G. Bodo, "An HLLC Solver for Relativistic Flows - II:
00015  *   Magnetohydrodynamics", arxiv:astro-ph/0601640v1 (2006)
00016  *
00017  * HISTORY: Written by Jonathan FUlton, February 2009
00018  *          Extended to MHD by Kris Beckwith, Spring 2010
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 HLLC_FLUX
00035 #ifdef SPECIAL_RELATIVITY
00036 
00037 #if (NSCALARS > 0)
00038 #error : The SR HLLC flux does not work with passive scalars.
00039 #endif
00040 
00041 #ifdef MHD
00042 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p);
00043 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00044                         const Real Bx, Real* low, Real* high);
00045 void getMaxSignalSpeeds_echo(const Prim1DS Wl, const Prim1DS Wr,
00046                          const Real Bx, Real* low, Real* high);
00047 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00048 void getVChar_echo (const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00049 /* solves quartic equation defined by a and returns roots in root
00050  * returns the number of Real roots 
00051  * error specifies an accuracy
00052  * currently force four Real solutions b/c it's physical */
00053 int QUARTIC (Real b, Real c, Real d, Real e, Real z[]);
00054 
00055 /* solves cubic equation defined by a and stores roots in root
00056  * returns number of Real roots */
00057 int CUBIC(Real b, Real c, Real d, Real z[]);
00058 #endif
00059 
00060 #ifdef HYDRO
00061 /*----------------------------------------------------------------------------*/
00062 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00063  *      const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00064  *  \brief Computes 1D fluxes
00065  *   Input Arguments:
00066  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface 
00067  *   - Wl,Wr = L/R-states of PRIMITIVE variables at cell interface 
00068  *   Output Arguments:
00069  *   - pFlux = pointer to fluxes of CONSERVED variables at cell interface 
00070  */
00071 
00072 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00073             const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00074 {
00075   Cons1DS Fl,Fr,Fhll,Uhll,Usl,Usr;
00076   Real rhl, rhr, csl, csr, cslsq, csrsq, vsql, vsqr, gammasql, gammasqr;
00077   Real ssl, ssr, radl, radr, lmdapl, lmdapr, lmdaml, lmdamr, lmdatlmda;
00078   Real lmdal,lmdar; /* Left and Right wave speeds */
00079   Real lmdas; /* Contact wave speed */
00080   Real ovlrmll;
00081   Real a,b,c,quad,rad;
00082   Real den,ps; /* PressUre in inner region */
00083 
00084 /*--- Step 1. ------------------------------------------------------------------
00085  * Compute the max and min wave speeds used in Mignone 
00086  */
00087 
00088   rhl = Wl.d + Wl.P * Gamma / Gamma_1; /* Mignone Eq 3.5 */
00089   rhr = Wr.d + Wr.P * Gamma / Gamma_1;
00090 
00091   csl = sqrt(Gamma * Wl.P / rhl); /* Mignone Eq 4 */
00092   csr = sqrt(Gamma * Wr.P / rhr);
00093 
00094   cslsq = SQR(csl);
00095   csrsq = SQR(csr);
00096 
00097   vsql = SQR(Wl.Vx) + SQR(Wl.Vy) + SQR(Wl.Vz);
00098   vsqr = SQR(Wr.Vx) + SQR(Wr.Vy) + SQR(Wr.Vz);
00099 
00100   gammasql = 1.0 / (1.0 - vsql);
00101   gammasqr = 1.0 / (1.0 - vsqr);
00102 
00103   ssl = cslsq / ( gammasql * (1.0 - cslsq) ); /* Mignone Eq 22.5 */
00104   ssr = csrsq / ( gammasqr * (1.0 - csrsq) );
00105 
00106   radl = sqrt( ssl*(1.0-SQR(Wl.Vx)+ssl) ); /* Mignone Eq 23 (radical part) */
00107   radr = sqrt( ssr*(1.0-SQR(Wr.Vx)+ssr) );
00108 
00109   lmdapl = (Wl.Vx + radl) / (1.0 + ssl); /* Mignone Eq 23 */
00110   lmdapr = (Wr.Vx + radr) / (1.0 + ssr);
00111   lmdaml = (Wl.Vx - radl) / (1.0 + ssl);
00112   lmdamr = (Wr.Vx - radr) / (1.0 + ssr);
00113 
00114   lmdal = MIN(lmdaml, lmdamr); /* Mignone Eq 21 */
00115   lmdar = MAX(lmdapl, lmdapr);
00116   
00117 
00118 /*--- Step 2. ------------------------------------------------------------------
00119  * Compute L/R fluxes according to Mignone 2
00120  */
00121 
00122   Fl.d  = Ul.d * Wl.Vx;
00123   Fl.Mx = Ul.Mx * Wl.Vx + Wl.P;
00124   Fl.My = Ul.My * Wl.Vx;
00125   Fl.Mz = Ul.Mz * Wl.Vx;
00126   Fl.E  = Ul.Mx;
00127 
00128   Fr.d  = Ur.d * Wr.Vx;
00129   Fr.Mx = Ur.Mx * Wr.Vx + Wr.P;
00130   Fr.My = Ur.My * Wr.Vx;
00131   Fr.Mz = Ur.Mz * Wr.Vx;
00132   Fr.E  = Ur.Mx;
00133 
00134 /*--- Step 3. ------------------------------------------------------------------
00135  * Compute HLL flux using Mignone Eq 11 (necessary for computing lmdas (Eq 18)
00136  * Compute HLL conserved quantities using Mignone eq 9
00137  */
00138 
00139   ovlrmll = 1.0 / ( lmdar - lmdal );
00140   lmdatlmda = lmdal*lmdar;
00141 
00142   Fhll.d  = (lmdar*Fl.d  - lmdal*Fr.d  + lmdatlmda * (Ur.d  - Ul.d) ) * ovlrmll;
00143   Fhll.Mx = (lmdar*Fl.Mx - lmdal*Fr.Mx + lmdatlmda * (Ur.Mx - Ul.Mx)) * ovlrmll;
00144   Fhll.My = (lmdar*Fl.My - lmdal*Fr.My + lmdatlmda * (Ur.My - Ul.My)) * ovlrmll;
00145   Fhll.Mz = (lmdar*Fl.Mz - lmdal*Fr.Mz + lmdatlmda * (Ur.Mz - Ul.Mz)) * ovlrmll;
00146   Fhll.E  = (lmdar*Fl.E  - lmdal*Fr.E  + lmdatlmda * (Ur.E  - Ul.E )) * ovlrmll;
00147 
00148   Uhll.d  = (lmdar * Ur.d  - lmdal * Ul.d  + Fl.d  - Fr.d ) * ovlrmll;
00149   Uhll.Mx = (lmdar * Ur.Mx - lmdal * Ul.Mx + Fl.Mx - Fr.Mx) * ovlrmll;
00150   Uhll.My = (lmdar * Ur.My - lmdal * Ul.My + Fl.My - Fr.My) * ovlrmll;
00151   Uhll.Mz = (lmdar * Ur.Mz - lmdal * Ul.Mz + Fl.Mz - Fr.Mz) * ovlrmll;
00152   Uhll.E  = (lmdar * Ur.E  - lmdal * Ul.E  + Fl.E  - Fr.E ) * ovlrmll;
00153 
00154 /*--- Step 4. ------------------------------------------------------------------
00155  * Compute contact wave speed using larger root from Mignone Eq 18
00156  * Physical root is the root with the minus sign
00157  */
00158 
00159   /* quadratic formUla calcUlation */
00160 
00161   a = Fhll.E;
00162   b = -(Uhll.E + Fhll.Mx);
00163   c = Uhll.Mx;
00164 
00165 
00166   quad = -0.5*(b + SIGN(b)*sqrt(b*b - 4.0*a*c));
00167   lmdas = c/quad;
00168 
00169 /*--- Step 5. ------------------------------------------------------------------
00170  * Determine intercell flux according to Mignone 13
00171  */
00172 
00173   if( lmdal >= 0.0){ /* Fl */
00174     /* intercell flux is left flux */
00175     pFlux->d  = Fl.d;
00176     pFlux->Mx = Fl.Mx;
00177     pFlux->My = Fl.My;
00178     pFlux->Mz = Fl.Mz;
00179     pFlux->E  = Fl.E;
00180 
00181    return;
00182   }
00183   else if( lmdas >= 0.0){ /* Fls */
00184 
00185     /* Mignone 2006 Eq 48 */
00186     ps = -Fhll.E*lmdas + Fhll.Mx;
00187 
00188     /* now calcUlate Usl with Mignone Eq 16 */
00189     den = 1.0 / (lmdal - lmdas);
00190 
00191     Usl.d  = Ul.d * (lmdal - Wl.Vx) * den;
00192     Usl.Mx = (Ul.Mx * (lmdal - Wl.Vx) + ps - Wl.P) * den;
00193     Usl.My = Ul.My * (lmdal - Wl.Vx) * den;
00194     Usl.Mz = Ul.Mz * (lmdal - Wl.Vx) * den;
00195     Usl.E  = (Ul.E * (lmdal - Wl.Vx) + ps * lmdas - Wl.P * Wl.Vx) * den;
00196 
00197     /* now calcUlate Fsr using Mignone Eq 14 */
00198 
00199     pFlux->d  = lmdal*(Usl.d  - Ul.d ) + Fl.d;
00200     pFlux->Mx = lmdal*(Usl.Mx - Ul.Mx) + Fl.Mx;
00201     pFlux->My = lmdal*(Usl.My - Ul.My) + Fl.My;
00202     pFlux->Mz = lmdal*(Usl.Mz - Ul.Mz) + Fl.Mz;
00203     pFlux->E  = lmdal*(Usl.E  - Ul.E ) + Fl.E;
00204 
00205     return;
00206   }
00207   else if( lmdar >= 0.0){ /* Frs */
00208 
00209     /* Mignone 2006 Eq 48 */
00210     ps = -Fhll.E*lmdas + Fhll.Mx;
00211 
00212     /* now calcUlate Usr with Mignone Eq 16 */
00213     den = 1.0 / (lmdar - lmdas);
00214 
00215     Usr.d  = Ur.d * (lmdar - Wr.Vx) * den;
00216     Usr.Mx = (Ur.Mx * (lmdar - Wr.Vx) + ps - Wr.P) * den;
00217     Usr.My = Ur.My * (lmdar - Wr.Vx) * den;
00218     Usr.Mz = Ur.Mz * (lmdar - Wr.Vx) * den;
00219     Usr.E  = (Ur.E * (lmdar - Wr.Vx) + ps * lmdas - Wr.P * Wr.Vx) * den;
00220 
00221     /* now calcUlate Fsr using Mignone Eq 14 */
00222 
00223     pFlux->d  = lmdar*(Usr.d  - Ur.d ) + Fr.d;
00224     pFlux->Mx = lmdar*(Usr.Mx - Ur.Mx) + Fr.Mx;
00225     pFlux->My = lmdar*(Usr.My - Ur.My) + Fr.My;
00226     pFlux->Mz = lmdar*(Usr.Mz - Ur.Mz) + Fr.Mz;
00227     pFlux->E  = lmdar*(Usr.E  - Ur.E ) + Fr.E;
00228 
00229     return;
00230   }
00231   else{ /* Fr */
00232     /* intercell flux is right flux */
00233     pFlux->d  = Fr.d;
00234     pFlux->Mx = Fr.Mx;
00235     pFlux->My = Fr.My;
00236     pFlux->Mz = Fr.Mz;
00237     pFlux->E  = Fr.E;
00238 
00239     return;
00240   }
00241   
00242   /* need to deal with scalar fluxes */
00243 }
00244 
00245 #endif /* HYDRO */
00246 
00247 #ifdef MHD
00248 /*----------------------------------------------------------------------------*/
00249 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00250  *          const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00251  *  \brief Computes 1D fluxes
00252  *   Input Arguments:
00253  *  -  Ul,Ur = L/R-states of CONSERVED variables at cell interface 
00254  *  -  Wl,Wr = L/R-states of PRIMITIVE variables at cell interface 
00255  *   Output Arguments:
00256  *  -  pFlux = pointer to fluxes of CONSERVED variables at cell interface 
00257  */
00258 
00259 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00260             const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00261 {
00262   Cons1DS Fl,Fr,Fhll,Uhll,Usl,Usr;
00263   Prim1DS Wsl,Wsr,Whll;
00264   Real Sl, Sr, Pl, Pr;
00265   Real Sla, Sra;
00266   Real dS_1, scrh;
00267   Real Bx, Bys, Bzs;
00268   Real BtFBt, Bt2, FBt2;
00269   Real a, b, c;
00270   Real ps, vxs, vys, vzs, gammas_2, vBs, V2l, V2r;
00271   Real vxl, vxr, alpha_l, alpha_r;
00272   int switch_to_hll,wave_speed_fail;
00273 
00274   wave_speed_fail = 0;
00275   switch_to_hll = 0;
00276         
00277 /*--- Step 1. ------------------------------------------------------------------
00278  * Compute the max and min wave speeds used in Mignone 
00279  */
00280   getMaxSignalSpeeds_pluto(Wl,Wr,Bxi,&Sl,&Sr);
00281         
00282   if (Sl != Sl) {
00283     wave_speed_fail = 1;
00284     printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00285     Sl = -1.0;
00286     Sr =  1.0;
00287   }
00288         
00289   if (Sr != Sr) {
00290     wave_speed_fail = 1;
00291     printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00292     Sl = -1.0;
00293     Sr = 1.0;
00294   }
00295         
00296   if (Sl < -1.0) {
00297     wave_speed_fail = 1;
00298     printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00299     Sl = -1.0;
00300     Sr = 1.0;
00301   }
00302   if (Sr > 1.0) {
00303     wave_speed_fail = 1;
00304     printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00305     Sl = -1.0;
00306     Sr = 1.0;
00307   }
00308 
00309 /*--- Step 1a. -----------------------------------------------------------------
00310  * If PLUTO wavespeeds are bad, fall back to the estimate used in ECHO
00311  */
00312   if (wave_speed_fail){
00313     getMaxSignalSpeeds_echo (Wl,Wr,Bxi,&Sla,&Sra);
00314         
00315     if (Sla != Sla) {
00316       switch_to_hll = 1;
00317       printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00318       Sla = -1.0;
00319       Sra =  1.0;
00320     }
00321         
00322     if (Sra != Sra) {
00323       switch_to_hll = 1;
00324       printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00325       Sla = -1.0;
00326       Sra = 1.0;
00327     }
00328         
00329     if (Sla < -1.0) {
00330       switch_to_hll = 1;
00331       printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00332       Sla = -1.0;
00333       Sra = 1.0;
00334     }
00335     if (Sra > 1.0) {
00336       switch_to_hll = 1;
00337       printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00338       Sla = -1.0;
00339       Sra = 1.0;
00340     }
00341 
00342     Sl = Sla;
00343     Sr = Sra;
00344 
00345   }
00346 
00347   /* compute L/R fluxes */
00348   flux_LR(Ul,Wl,&Fl,Bxi,&Pl);
00349   flux_LR(Ur,Wr,&Fr,Bxi,&Pr);
00350         
00351 /*--- Step 2. ------------------------------------------------------------------
00352  * Construct HLL fluxes & average state
00353  */
00354   dS_1 = 1.0/(Sr - Sl);
00355 
00356   Uhll.d  = (Sr*Ur.d  - Sl*Ul.d  + Fl.d  - Fr.d ) * dS_1;
00357   Uhll.Mx = (Sr*Ur.Mx - Sl*Ul.Mx + Fl.Mx - Fr.Mx) * dS_1;
00358   Uhll.My = (Sr*Ur.My - Sl*Ul.My + Fl.My - Fr.My) * dS_1;
00359   Uhll.Mz = (Sr*Ur.Mz - Sl*Ul.Mz + Fl.Mz - Fr.Mz) * dS_1;
00360   Uhll.E  = (Sr*Ur.E  - Sl*Ul.E  + Fl.E  - Fr.E ) * dS_1;
00361   Uhll.By = (Sr*Ur.By - Sl*Ul.By + Fl.By - Fr.By) * dS_1;
00362   Uhll.Bz = (Sr*Ur.Bz - Sl*Ul.Bz + Fl.Bz - Fr.Bz) * dS_1;
00363                 
00364   Fhll.d  = (Sr*Fl.d  - Sl*Fr.d  + Sl*Sr*(Ur.d  - Ul.d )) * dS_1;
00365   Fhll.Mx = (Sr*Fl.Mx - Sl*Fr.Mx + Sl*Sr*(Ur.Mx - Ul.Mx)) * dS_1;
00366   Fhll.My = (Sr*Fl.My - Sl*Fr.My + Sl*Sr*(Ur.My - Ul.My)) * dS_1;
00367   Fhll.Mz = (Sr*Fl.Mz - Sl*Fr.Mz + Sl*Sr*(Ur.Mz - Ul.Mz)) * dS_1;
00368   Fhll.E  = (Sr*Fl.E  - Sl*Fr.E  + Sl*Sr*(Ur.E  - Ul.E )) * dS_1;
00369   Fhll.By = (Sr*Fl.By - Sl*Fr.By + Sl*Sr*(Ur.By - Ul.By)) * dS_1;
00370   Fhll.Bz = (Sr*Fl.Bz - Sl*Fr.Bz + Sl*Sr*(Ur.Bz - Ul.Bz)) * dS_1;
00371 
00372 
00373   if (switch_to_hll) {
00374     pFlux->d = Fhll.d;
00375     pFlux->Mx = Fhll.Mx;
00376     pFlux->My = Fhll.My;
00377     pFlux->Mz = Fhll.Mz;
00378     pFlux->E = Fhll.E;
00379     pFlux->By = Fhll.By;
00380     pFlux->Bz = Fhll.Bz;
00381                 
00382     return;
00383   }
00384 
00385         
00386 /*--- Step 3. ------------------------------------------------------------------
00387  * Compute fluxes based on wave speeds (Mignone et al. eqn 26)
00388  */
00389   if(Sl >= 0.0){
00390     pFlux->d  = Fl.d;
00391     pFlux->Mx = Fl.Mx;
00392     pFlux->My = Fl.My;
00393     pFlux->Mz = Fl.Mz;
00394     pFlux->E  = Fl.E;
00395     pFlux->By = Fl.By;
00396     pFlux->Bz = Fl.Bz;
00397                 
00398     return;
00399    
00400   }
00401   else if(Sr <= 0.0){
00402     pFlux->d  = Fr.d;
00403     pFlux->Mx = Fr.Mx;
00404     pFlux->My = Fr.My;
00405     pFlux->Mz = Fr.Mz;
00406     pFlux->E  = Fr.E;
00407     pFlux->By = Fr.By;
00408     pFlux->Bz = Fr.Bz;
00409                 
00410     return;
00411   }
00412   else {
00413 
00414     switch_to_hll = 0;
00415                 
00416     /* Construct HLLC fluxes */
00417     vxl = Wl.Vx;
00418     vxr = Wr.Vx;
00419                 
00420     Bx  = Bxi;
00421     Bys = Uhll.By;
00422     Bzs = Uhll.Bz;
00423                 
00424     if (fabs(Bx) < 1.0e-12) {
00425       a  = Fhll.E;
00426       b  = - (Fhll.Mx + Uhll.E);
00427       c  = Uhll.Mx;
00428     } else {
00429       BtFBt = Uhll.By*Fhll.By + Uhll.Bz*Fhll.Bz;
00430       Bt2 = Uhll.By*Uhll.By + Uhll.Bz*Uhll.Bz;
00431       FBt2 = Fhll.By*Fhll.By + Fhll.Bz*Fhll.Bz;                
00432                         
00433       a  = Fhll.E - BtFBt;
00434       b  = Bt2 + FBt2 - (Fhll.Mx + Uhll.E);
00435       c  = Uhll.Mx - BtFBt;
00436     }
00437 
00438     if (fabs(a) > 1.e-12){
00439       scrh = 1.0 + sqrt(1.0 - 4.0*a*c/(b*b));
00440       if (scrh != scrh) {
00441         switch_to_hll = 1;
00442       }
00443       vxs  = - 2.0*c/(b*scrh);
00444     } else {
00445       vxs = -c/b;
00446     }
00447     if ((vxs != vxs || vxs > 1.0) && (switch_to_hll == 0)) {
00448       switch_to_hll = 1;
00449     }
00450                 
00451     if (switch_to_hll) {
00452       pFlux->d = Fhll.d;
00453       pFlux->Mx = Fhll.Mx;
00454       pFlux->My = Fhll.My;
00455       pFlux->Mz = Fhll.Mz;
00456       pFlux->E = Fhll.E;
00457       pFlux->By = Fhll.By;
00458       pFlux->Bz = Fhll.Bz;
00459 
00460       if (pFlux->d != pFlux->d) {
00461       }
00462                 
00463       return;
00464     }
00465                 
00466     if (fabs(Bx) < 1.0e-12) {
00467                         
00468       /* -------------------------------
00469          the value of vy and vz
00470          is irrelevant in this case  
00471          ------------------------------- */
00472                         
00473       ps  = Fhll.Mx - Fhll.E*vxs;
00474 
00475       if (ps < 0) {
00476         switch_to_hll = 1;
00477       } else {                  
00478         alpha_l = (Sl - vxl)/(Sl - vxs);
00479         alpha_r = (Sr - vxr)/(Sr - vxs);
00480                         
00481         Usl.d = Ul.d*alpha_l;
00482         Usr.d = Ur.d*alpha_r;
00483                         
00484         Usl.E = (Sl*Ul.E - Fl.E + ps*vxs)/(Sl - vxs);
00485         Usr.E = (Sr*Ur.E - Fr.E + ps*vxs)/(Sr - vxs);
00486                         
00487         Usl.Mx = (Usl.E + ps)*vxs; 
00488         Usr.Mx = (Usr.E + ps)*vxs;
00489         Usl.My = Ul.My*alpha_l; 
00490         Usr.My = Ur.My*alpha_r; 
00491         Usl.Mz = Ul.Mz*alpha_l; 
00492         Usr.Mz = Ur.Mz*alpha_r;
00493                         
00494         Usl.By = Ul.By*alpha_l;
00495         Usr.By = Ur.By*alpha_r;
00496         Usl.Bz = Ul.Bz*alpha_l;
00497         Usr.Bz = Ur.Bz*alpha_r;
00498       }
00499                 
00500     } else {
00501 
00502       vys = (Bys*vxs - Fhll.By)/Bx;
00503       vzs = (Bzs*vxs - Fhll.Bz)/Bx;
00504 
00505       gammas_2 = vxs*vxs + vys*vys + vzs*vzs;
00506       gammas_2 = 1.0 - gammas_2;
00507       vBs = vxs*Bx + vys*Bys + vzs*Bzs;
00508                         
00509       ps = (Bx*vBs - Fhll.E)*vxs + (Bx*Bx*gammas_2) + Fhll.Mx;
00510       if (ps < 0) {
00511         switch_to_hll = 1;
00512       } else {
00513                         
00514         alpha_l = (Sl - vxl)/(Sl - vxs);
00515         alpha_r = (Sr - vxr)/(Sr - vxs);
00516                         
00517         if (alpha_l != alpha_l) {
00518           switch_to_hll = 1;
00519         }
00520                         
00521         if (alpha_r != alpha_r) {
00522           switch_to_hll = 1;
00523 
00524         }
00525 
00526         if (switch_to_hll == 0){
00527                         
00528           Usl.d = Ul.d*alpha_l;
00529           Usr.d = Ur.d*alpha_r;
00530                         
00531           Usl.E = (Sl*Ul.E - Fl.E + ps*vxs - vBs*Bx)/(Sl - vxs);
00532           Usr.E = (Sr*Ur.E - Fr.E + ps*vxs - vBs*Bx)/(Sr - vxs);
00533                         
00534           Usl.Mx = (Usl.E + ps)*vxs - vBs*Bx; 
00535           Usr.Mx = (Usr.E + ps)*vxs - vBs*Bx;
00536           Usl.My = (Sl*Ul.My - Fl.My - Bx*(Bys*gammas_2 + vBs*vys))/(Sl - vxs); 
00537           Usr.My = (Sr*Ur.My - Fr.My - Bx*(Bys*gammas_2 + vBs*vys))/(Sr - vxs);  
00538           Usl.Mz = (Sl*Ul.Mz - Fl.Mz - Bx*(Bzs*gammas_2 + vBs*vzs))/(Sl - vxs); 
00539           Usr.Mz = (Sr*Ur.Mz - Fr.Mz - Bx*(Bzs*gammas_2 + vBs*vzs))/(Sr - vxs);
00540                         
00541           Usl.By = Usr.By = Bys;
00542           Usl.Bz = Usr.Bz = Bzs;
00543         }
00544       }
00545     }
00546 
00547     if (switch_to_hll) {
00548       pFlux->d = Fhll.d;
00549       pFlux->Mx = Fhll.Mx;
00550       pFlux->My = Fhll.My;
00551       pFlux->Mz = Fhll.Mz;
00552       pFlux->E = Fhll.E;
00553       pFlux->By = Fhll.By;
00554       pFlux->Bz = Fhll.Bz;
00555                 
00556       return;
00557     }
00558                 
00559     /*  ----  Compute HLLC flux  ----  */
00560     if (vxs > 0.0) {
00561       pFlux->d  = Fl.d  + Sl*(Usl.d  - Ul.d );
00562       pFlux->Mx = Fl.Mx + Sl*(Usl.Mx - Ul.Mx);
00563       pFlux->My = Fl.My + Sl*(Usl.My - Ul.My);
00564       pFlux->Mz = Fl.Mz + Sl*(Usl.Mz - Ul.Mz);
00565       pFlux->E  = Fl.E  + Sl*(Usl.E  - Ul.E );
00566       pFlux->By = Fl.By + Sl*(Usl.By - Ul.By);
00567       pFlux->Bz = Fl.Bz + Sl*(Usl.Bz - Ul.Bz);
00568                         
00569       if (pFlux->d != pFlux->d) {
00570 
00571         pFlux->d = Fhll.d;
00572         pFlux->Mx = Fhll.Mx;
00573         pFlux->My = Fhll.My;
00574         pFlux->Mz = Fhll.Mz;
00575         pFlux->E = Fhll.E;
00576         pFlux->By = Fhll.By;
00577         pFlux->Bz = Fhll.Bz;
00578 
00579       }
00580                         
00581       return;
00582                         
00583     } else {
00584         pFlux->d  = Fr.d  + Sr*(Usr.d  - Ur.d );
00585         pFlux->Mx = Fr.Mx + Sr*(Usr.Mx - Ur.Mx);
00586         pFlux->My = Fr.My + Sr*(Usr.My - Ur.My);
00587         pFlux->Mz = Fr.Mz + Sr*(Usr.Mz - Ur.Mz);
00588         pFlux->E  = Fr.E  + Sr*(Usr.E  - Ur.E );
00589         pFlux->By = Fr.By + Sr*(Usr.By - Ur.By);
00590         pFlux->Bz = Fr.Bz + Sr*(Usr.Bz - Ur.Bz);
00591                         
00592       if (pFlux->d != pFlux->d) {
00593 
00594         pFlux->d = Fhll.d;
00595         pFlux->Mx = Fhll.Mx;
00596         pFlux->My = Fhll.My;
00597         pFlux->Mz = Fhll.Mz;
00598         pFlux->E = Fhll.E;
00599         pFlux->By = Fhll.By;
00600         pFlux->Bz = Fhll.Bz;
00601 
00602       }
00603                         
00604       return;
00605     }
00606 
00607   }
00608 
00609 }
00610 
00611 /*! \fn void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00612  *            const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00613  *  \brief Calculate entropy flux */
00614 void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00615             const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00616 {
00617   Real Fl, Fr;
00618   Real USl, USr;
00619   Real WSl, WSr;
00620   Real Uhll, Fhll;
00621   Real Pl, Pr;
00622   Real Sl, Sr;
00623   Real Sla, Sra;
00624   Real dS_1;
00625   int wave_speed_fail;
00626 
00627   wave_speed_fail = 0;
00628         
00629 /*--- Step 1. ------------------------------------------------------------------
00630  * Compute the max and min wave speeds used in Mignone 
00631  */
00632   getMaxSignalSpeeds_pluto(Wl,Wr,Bx,&Sl,&Sr);
00633         
00634   if (Sl != Sl) {
00635     wave_speed_fail = 1;
00636     printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00637     Sl = -1.0;
00638     Sr =  1.0;
00639   }
00640         
00641   if (Sr != Sr) {
00642     wave_speed_fail = 1;
00643     printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00644     Sl = -1.0;
00645     Sr = 1.0;
00646   }
00647         
00648   if (Sl < -1.0) {
00649     wave_speed_fail = 1;
00650     printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00651     Sl = -1.0;
00652     Sr = 1.0;
00653   }
00654   if (Sr > 1.0) {
00655     wave_speed_fail = 1;
00656     printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00657     Sl = -1.0;
00658     Sr = 1.0;
00659   }
00660 
00661 /*--- Step 1a. -----------------------------------------------------------------
00662  * If PLUTO wavespeeds are bad, fall back to the estimate used in ECHO
00663  */
00664   if (wave_speed_fail){
00665     getMaxSignalSpeeds_echo (Wl,Wr,Bx,&Sla,&Sra);
00666         
00667     if (Sla != Sla) {
00668       printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00669       Sla = -1.0;
00670       Sra =  1.0;
00671     }
00672         
00673     if (Sra != Sra) {
00674       printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00675       Sla = -1.0;
00676       Sra = 1.0;
00677     }
00678         
00679     if (Sla < -1.0) {
00680       printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00681       Sla = -1.0;
00682       Sra = 1.0;
00683     }
00684     if (Sra > 1.0) {
00685       printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00686       Sla = -1.0;
00687       Sra = 1.0;
00688     }
00689 
00690     Sl = Sla;
00691     Sr = Sra;
00692 
00693   }
00694 
00695   /* compute L/R fluxes */
00696   WSl = Wl.P*pow(Wl.d,1.0-Gamma);
00697   WSr = Wr.P*pow(Wr.d,1.0-Gamma);
00698   USl = WSl * Ul.d/Wl.d;
00699   USr = WSr * Ur.d/Wr.d;
00700   Fl = USl * Wl.Vx;
00701   Fr = USr * Wr.Vx;
00702 
00703   if(Sl >= 0.0){
00704     *pFlux = Fl;
00705     return;
00706   }
00707   else if(Sr <= 0.0){
00708     *pFlux = Fr;
00709     return;
00710   }
00711   else{
00712     /* Compute HLL average state */
00713     dS_1 = 1.0/(Sr - Sl);
00714     *pFlux = (Sr*Fl  - Sl*Fr  + Sl*Sr*(USr  - USl)) * dS_1;
00715     return;
00716   }
00717 }
00718 
00719 /*! \fn void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p)
00720  *  \brief Calculate LR flux. */
00721 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p)
00722 {
00723   Real wtg2, pt, g, g2, g_1,g_2, h, gmmr, theta;
00724   Real bx, by, bz, vB, b2, Bmag2;
00725         
00726   /* calcUlate enthalpy */
00727         
00728   theta = W.P/W.d;
00729   gmmr = Gamma / Gamma_1;
00730         
00731   h = 1.0 + gmmr*theta;
00732         
00733   /* calcUlate gamma */
00734   g   = U.d/W.d;
00735   g2  = SQR(g);
00736   g_1 = 1.0/g;
00737   g_2 = 1.0/g2;
00738         
00739   pt = W.P;
00740   wtg2 = W.d*h*g2;
00741         
00742   vB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00743   Bmag2 = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00744         
00745   bx = g*(  Bx*g_2 + vB*W.Vx);
00746   by = g*(W.By*g_2 + vB*W.Vy);
00747   bz = g*(W.Bz*g_2 + vB*W.Vz);
00748         
00749   b2 = Bmag2*g_2 + vB*vB;
00750         
00751   pt += 0.5*b2;
00752   wtg2 += b2*g2;
00753         
00754   flux->d  = U.d*W.Vx;
00755   flux->Mx = wtg2*W.Vx*W.Vx + pt;
00756   flux->My = wtg2*W.Vy*W.Vx;
00757   flux->Mz = wtg2*W.Vz*W.Vx;
00758   flux->E  = U.Mx;
00759 
00760   flux->Mx -= bx*bx;
00761   flux->My -= by*bx;
00762   flux->Mz -= bz*bx;
00763   flux->By = W.Vx*W.By - Bx*W.Vy;
00764   flux->Bz = W.Vx*W.Bz - Bx*W.Vz;
00765         
00766   *p = pt;
00767 }
00768 
00769 /*! \fn void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00770  *                            const Real Bx, Real* low, Real* high)
00771  *  \brief
00772  */
00773 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00774                               const Real Bx, Real* low, Real* high)
00775 {
00776         
00777   Real lml,lmr;        /* smallest roots, Mignone Eq 55 */
00778   Real lpl,lpr;        /* largest roots, Mignone Eq 55 */
00779   Real al,ar;
00780         
00781   getVChar_pluto(Wl,Bx,&lml,&lpl);
00782   getVChar_pluto(Wr,Bx,&lmr,&lpr);
00783         
00784   *low =  MIN(lml, lmr);
00785   *high = MAX(lpl, lpr);
00786 }
00787 
00788 /*! \fn void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00789  *  \brief
00790  */
00791 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00792 {
00793   Real rhoh,vsq,bsq;
00794   Real cssq,vasq,asq;
00795   Real Vx2,gamma,gamma2;
00796   Real Bx2,Bsq,vDotB,vDotBsq,b0,bx;
00797   Real w_1,a0,a1,a2,a3,a4,Q;
00798   Real scrh,scrh1,scrh2,eps2;
00799   Real a2_w,one_m_eps2,lambda[5];
00800   int iflag;
00801         
00802   rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00803         
00804   Vx2 = SQR(W.Vx);
00805   vsq = Vx2 + SQR(W.Vy) + SQR(W.Vz);
00806   if (vsq > 1.0){
00807     /*printf("[getVChar]: |v|= %f > 1\n",vsq);*/        
00808     *lm = -1.0;
00809     *lp = 1.0;  
00810     return;             
00811   }
00812   gamma = 1.0 / sqrt(1 - vsq);    
00813   gamma2 = SQR(gamma);
00814         
00815   Bx2 = SQR(Bx);
00816   Bsq = Bx2 + SQR(W.By) + SQR(W.Bz);
00817   vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00818   vDotBsq = SQR(vDotB);
00819   b0 = gamma * vDotB;
00820   bx = Bx/gamma2 + W.Vx*vDotB;
00821   bsq = Bsq / gamma2 + SQR(vDotB);
00822         
00823   cssq = (Gamma * W.P) / (rhoh);
00824   vasq = bsq / (rhoh + bsq);
00825         
00826   if (bsq < 0.0) bsq = 0.0;
00827   if (cssq < 0.0) cssq = 0.0;
00828   if (cssq > 1.0) cssq = 1.0;
00829   if (vasq > 1.0) bsq = rhoh + bsq;
00830         
00831   if (vsq < 1.0e-12) {
00832     w_1  = 1.0/(rhoh + bsq);   
00833     eps2 = cssq + bsq*w_1*(1.0 - cssq);
00834     a0   = cssq*Bx*Bx*w_1;
00835     a1   = - a0 - eps2;
00836     scrh = a1*a1 - 4.0*a0;
00837     if (scrh < 0.0) scrh = 0.0;
00838                 
00839     scrh = sqrt(0.5*(-a1 + sqrt(scrh)));
00840     *lp =  scrh;
00841     *lm = -scrh;
00842     return;
00843   }
00844         
00845   w_1 = 1.0/(rhoh + bsq);   
00846         
00847   if (Bx < 1.0e-14) {
00848                 
00849     eps2  = cssq + bsq*w_1*(1.0 - cssq);
00850                 
00851     scrh1 = (1.0 - eps2)*gamma2;
00852     scrh2 = cssq*vDotBsq*w_1 - eps2;
00853                 
00854     a2  = scrh1 - scrh2;
00855     a1  = -2.0*W.Vx*scrh1;
00856     a0  = Vx2*scrh1 + scrh2;
00857                 
00858     *lp = 0.5*(-a1 + sqrt(a1*a1 - 4.0*a2*a0))/a2;
00859     *lm = 0.5*(-a1 - sqrt(a1*a1 - 4.0*a2*a0))/a2;
00860                 
00861     return;
00862   }
00863         
00864   scrh1 = bx;  /* -- this is bx/u0 -- */
00865   scrh2 = scrh1*scrh1;  
00866         
00867   a2_w       = cssq*w_1;
00868   eps2       = (cssq*rhoh + bsq)*w_1;
00869   one_m_eps2 = gamma2*rhoh*(1.0 - cssq)*w_1;
00870         
00871   /* ---------------------------------------
00872      Define coefficients for the quartic  
00873      --------------------------------------- */
00874         
00875   scrh = 2.0*(a2_w*vDotB*scrh1 - eps2*W.Vx);
00876   a4 = one_m_eps2 - a2_w*vDotBsq + eps2;
00877   a3 = - 4.0*W.Vx*one_m_eps2 + scrh;
00878   a2 =   6.0*Vx2*one_m_eps2 + a2_w*(vDotBsq - scrh2) + eps2*(Vx2 - 1.0);
00879   a1 = - 4.0*W.Vx*Vx2*one_m_eps2 - scrh;
00880   a0 = Vx2*Vx2*one_m_eps2 + a2_w*scrh2 - eps2*Vx2;
00881         
00882   if (a4 < 1.e-12){
00883     /*printPrim1D(W);*/
00884     printf("[MAX_CH_SPEED]: Can not divide by a4 in MAX_CH_SPEED\n");
00885                 
00886     *lm = -1.0;
00887     *lp = 1.0;
00888                 
00889     return;
00890   }
00891         
00892   scrh = 1.0/a4;
00893         
00894   a3 *= scrh;
00895   a2 *= scrh;
00896   a1 *= scrh;
00897   a0 *= scrh;
00898   iflag = QUARTIC(a3, a2, a1, a0, lambda);
00899         
00900   if (iflag){
00901     printf ("Can not find max speed:\n");
00902     /*SHOW(uprim,i);*/
00903     printf("QUARTIC: f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
00904            a0*a4, a1*a4, a2*a4);
00905     printf("+ x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
00906     printf("[MAX_CH_SPEED]: Failed to find wave speeds");
00907                 
00908     *lm = -1.0;
00909     *lp = 1.0;
00910                 
00911     return;
00912   }
00913         
00914   *lp = MIN(1.0,MAX(lambda[3], lambda[2]));
00915   *lp = MIN(1.0,MAX(*lp, lambda[1]));
00916   *lp = MIN(1.0,MAX(*lp, lambda[0]));
00917 
00918   *lm = MAX(-1.0,MIN(lambda[3], lambda[2]));
00919   *lm = MAX(-1.0,MIN(*lm, lambda[1]));
00920   *lm = MAX(-1.0,MIN(*lm, lambda[0]));
00921         
00922   return;
00923         
00924 }
00925 
00926 /*! \fn void getMaxSignalSpeeds_echo (const Prim1DS Wl, const Prim1DS Wr,
00927  *                            const Real Bx, Real* low, Real* high)
00928  *  \brief
00929  */
00930 void getMaxSignalSpeeds_echo (const Prim1DS Wl, const Prim1DS Wr,
00931                               const Real Bx, Real* low, Real* high)
00932 {
00933         
00934   Real lml,lmr;        /* smallest roots, Mignone Eq 55 */
00935   Real lpl,lpr;        /* largest roots, Mignone Eq 55 */
00936   Real al,ar;
00937         
00938   getVChar_echo(Wl,Bx,&lml,&lpl);
00939   getVChar_echo(Wr,Bx,&lmr,&lpr);
00940         
00941   *low =  MIN(lml, lmr);
00942   *high = MAX(lpl, lpr);
00943 }
00944 
00945 /*! \fn void getVChar_echo(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00946  *  \brief
00947  */
00948 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00949 {
00950   Real rhoh,vsq,bsq;
00951   Real cssq,vasq,asq;
00952   Real gamma,gamma2;
00953   Real Bsq,vDotB,b0,bx;
00954   Real tmp1,tmp2,tmp3,tmp4,tmp5;
00955   Real vm,vp;
00956         
00957   rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00958         
00959   vsq = SQR(W.Vx) + SQR(W.Vy) + SQR(W.Vz);
00960   gamma = 1.0 / sqrt(1 - vsq);    
00961   gamma2 = SQR(gamma);
00962         
00963   Bsq = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00964   vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00965   b0 = gamma * vDotB;
00966   bx = Bx/gamma2 + W.Vx*vDotB;
00967   bsq = Bsq / gamma2 + SQR(vDotB);
00968         
00969   cssq = (Gamma * W.P) / (rhoh);
00970   vasq = bsq / (rhoh + bsq);
00971   asq = cssq + vasq - (cssq*vasq);
00972         
00973   if (cssq < 0.0) cssq = 0.0;
00974   if (vasq > 0.0) vasq = 0.0;
00975   if (asq < 0.0) asq = 0.0;
00976   if (cssq > 1.0) cssq = 1.0;
00977   if (vasq > 1.0) vasq = 1.0;
00978   if (asq > 1.0) asq = 1.0;
00979         
00980   tmp1 = (1.0 - asq);
00981   tmp2 = (1.0 - vsq);
00982   tmp3 = (1.0 - vsq*asq);
00983   tmp4 = SQR(W.Vx);
00984   tmp5 = 1.0 / tmp3;
00985         
00986   vm = tmp1*W.Vx - sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00987   vp = tmp1*W.Vx + sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00988   vm *=tmp5;
00989   vp *=tmp5;
00990         
00991   if (vp > vm) {
00992     *lm = vm;
00993     *lp = vp;
00994   } else {
00995     *lm = vp;
00996     *lp = vm;
00997   }
00998 }
00999 
01000 /* ******************************************** */
01001 /*! \fn int QUARTIC (Real b, Real c, Real d, Real e, Real z[])
01002  *  \brief Solve a quartic equation */
01003 /* 
01004  *
01005  * PURPOSE:
01006  *
01007  *   Solve a quartic equation in the form 
01008  *
01009  *    - z^4 + bz^3 + cz^2 + dz + e = 0
01010  *
01011  *   For its purpose, it is assumed that ALL 
01012  *   roots are real. This makes things faster.
01013  *
01014  *
01015  * ARGUMENTS
01016  *
01017  * - b, c,
01018  * - d, e  (IN)  = coefficient of the quartic
01019  *                 z^4 + bz^3 + cz^2 + dz + e = 0
01020  *
01021  * - z[]   (OUT) = a vector containing the 
01022  *                 (real) roots of the quartic
01023  *   
01024  *
01025  * REFERENCE:
01026  *
01027  * - http://www.1728.com/quartic2.htm 
01028  * 
01029  *
01030  */
01031 /********************************************** */
01032 int QUARTIC (Real b, Real c, Real d, 
01033              Real e, Real z[])
01034 {
01035   int    n, ifail;
01036   Real b2, f, g, h;
01037   Real a2, a1, a0, u[4];
01038   Real p, q, r, s;
01039   static Real three_256 = 3.0/256.0;
01040   static Real one_64 = 1.0/64.0;
01041         
01042   b2 = b*b;
01043         
01044   f = c - b2*0.375;
01045   g = d + b2*b*0.125 - b*c*0.5;
01046   h = e - b2*b2*three_256 + 0.0625*b2*c - 0.25*b*d;
01047         
01048   a2 = 0.5*f;
01049   a1 = (f*f - 4.0*h)*0.0625;
01050   a0 = -g*g*one_64;
01051         
01052   ifail = CUBIC(a2, a1, a0, u);
01053         
01054   if (ifail)return(1);
01055         
01056   if (u[1] < 1.e-14){
01057                 
01058     p = sqrt(u[2]);
01059     s = 0.25*b;
01060     z[0] = z[2] = - p - s;
01061     z[1] = z[3] = + p - s;
01062                 
01063   }else{
01064                 
01065     p = sqrt(u[1]);
01066     q = sqrt(u[2]);
01067                 
01068     r = -0.125*g/(p*q);
01069     s =  0.25*b;
01070                 
01071     z[0] = - p - q + r - s;
01072     z[1] =   p - q - r - s;
01073     z[2] = - p + q - r - s;
01074     z[3] =   p + q + r - s;
01075                 
01076   }  
01077         
01078   /* ----------------------------------------------
01079      verify that cmax and cmin satisfy original 
01080      equation
01081      ---------------------------------------------- */  
01082         
01083   for (n = 0; n < 4; n++){
01084     s = e + z[n]*(d + z[n]*(c + z[n]*(b + z[n])));
01085     if (s != s) {
01086       printf ("Nan found in QUARTIC \n");
01087       return(1);
01088     }
01089     if (fabs(s) > 1.e-6) {
01090       printf ("Solution does not satisfy f(z) = 0; f(z) = %12.6e\n",s);
01091       return(1);
01092     }
01093   }
01094         
01095   return(0);
01096   /*  
01097       printf (" z: %f ; %f ; %f ; %f\n",z[0], z[1], z[2], z[3]);
01098   */
01099 }
01100 /* *************************************************** */
01101 /*! \fn int CUBIC(Real b, Real c, Real d, Real z[])
01102  *  \brief Solve a cubic equation.
01103  *
01104  * PURPOSE:
01105  *
01106  *   Solve a cubic equation in the form 
01107  *
01108  *   -  z^3 + bz^2 + cz + d = 0
01109  *
01110  *   For its purpose, it is assumed that ALL 
01111  *   roots are real. This makes things faster.
01112  *
01113  *
01114  * ARGUMENTS
01115  *
01116  * - b, c, d (IN)  = coefficient of the cubic
01117  *                    z^3 + bz^2 + cz + d = 0
01118  *
01119  * - z[]   (OUT)   = a vector containing the 
01120  *                   (real) roots of the cubic.
01121  *                   Roots should be sorted
01122  *                   in increasing order.
01123  *   
01124  *
01125  * REFERENCE:
01126  *
01127  * - http://www.1728.com/cubic2.htm 
01128  *
01129  *
01130  *
01131  ***************************************************** */
01132 int CUBIC(Real b, Real c, Real d, Real z[])
01133 {
01134   Real b2, g2;
01135   Real f, g, h;
01136   Real i, i2, j, k, m, n, p;
01137   static Real one_3 = 1.0/3.0, one_27=1.0/27.0;
01138         
01139   b2 = b*b;
01140         
01141   /*  ----------------------------------------------
01142       the expression for f should be 
01143       f = c - b*b/3.0; however, to avoid negative
01144       round-off making h > 0.0 or g^2/4 - h < 0.0
01145       we let c --> c(1- 1.1e-16)
01146       ---------------------------------------------- */
01147         
01148   f  = c*(1.0 - 1.e-16) - b2*one_3;
01149   g  = b*(2.0*b2 - 9.0*c)*one_27 + d; 
01150   g2 = g*g;
01151   i2 = -f*f*f*one_27;
01152   h  = g2*0.25 - i2;
01153         
01154   /* --------------------------------------------
01155      Real roots are possible only when 
01156          
01157      h <= 0 
01158      -------------------------------------------- */
01159         
01160   if (h > 1.e-12){
01161     printf ("Only one real root (%12.6e)!\n", h);
01162   }
01163   if (i2 < 0.0){
01164     /*
01165       printf ("i2 < 0.0 %12.6e\n",i2);
01166       return(1);
01167     */
01168     i2 = 0.0;
01169   }
01170         
01171   /* --------------------------------------
01172      i^2 must be >= g2*0.25
01173      -------------------------------------- */
01174         
01175   i = sqrt(i2);       /*  > 0   */
01176   j = pow(i, one_3);  /*  > 0   */
01177   k = -0.5*g/i;
01178         
01179   /*  this is to prevent unpleseant situation 
01180       where both g and i are close to zero       */
01181         
01182   k = (k < -1.0 ? -1.0:k);
01183   k = (k >  1.0 ?  1.0:k);
01184         
01185   k = acos(k)*one_3;       /*  pi/3 < k < 0 */
01186         
01187   m = cos(k);              /*   > 0   */
01188   n = sqrt(3.0)*sin(k);    /*   > 0   */
01189   p = -b*one_3;
01190         
01191   z[0] = -j*(m + n) + p;
01192   z[1] = -j*(m - n) + p;
01193   z[2] =  2.0*j*m + p;
01194         
01195   /* ------------------------------------------------------
01196      Since j, m, n > 0, it should follow that from
01197     
01198      z0 = -jm - jn + p
01199      z1 = -jm + jn + p
01200      z2 = 2jm + p
01201          
01202      z2 is the greatest of the roots, while z0 is the 
01203      smallest one.
01204      ------------------------------------------------------ */
01205         
01206   return(0);
01207 }
01208 
01209 #endif /* MHD */
01210 
01211 #endif /* HLLC_FLUX */
01212 #endif /* SPECIAL_RELATIVITY */