Athena: rsolvers/hlld_sr.c Source File

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