Athena: rsolvers/exact_sr.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file exact_sr.c
00004  *  \brief Computes 1D fluxes using exact special relativistic Riemann solver.
00005  *
00006  * PURPOSE: Computes 1D fluxes using exact special relativistic Riemann solver.
00007  *   currently works only for hydrodynamics
00008  *
00009  * REFERENCES:
00010  * - Rezzolla, Zanotti, and Pons. "An Improved Exact Riemann Solver for
00011  *   Multidimensional Relativistic Flows." 2002. 
00012  *
00013  * HISTORY:
00014  * - April-2010:  Written by Nick Hand. 
00015  *
00016  * CONTAINS PUBLIC FUNCTIONS:
00017  * - fluxes() - all Riemann solvers in Athena must have this function name and
00018  *              use the same argument list as defined in rsolvers/prototypes.h
00019  *============================================================================*/
00020 
00021 #include <math.h>
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <float.h>
00025 #include "../defs.h"
00026 #include "../athena.h"
00027 #include "../globals.h"
00028 #include "prototypes.h"
00029 #include "../prototypes.h"
00030 
00031 #ifdef EXACT_FLUX
00032 #ifdef SPECIAL_RELATIVITY
00033 
00034 #ifdef MHD
00035 #error : The exact flux for MHD has not been implemented.
00036 #endif /* MHD */
00037 
00038 #if (NSCALARS > 0)
00039 #error : Passive scalars have not been implemented in the exact flux.
00040 #endif /* NSCALARS */
00041 
00042 #define EPS 3.0e-11 /* EPS is the relative precision. */
00043 #define JMAX 40     /* Maximum number of allowed bisections */
00044 enum WaveType {Two_S, RS, SR, Two_R, None}; 
00045 
00046 static enum WaveType wave; /* the wave pattern of the Riemann problem */ 
00047 
00048 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00049             const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pF);
00050 
00051 void gauleg(double x1, double x2, Real x[], Real w[], int n);
00052 
00053 Real rtbis_vel(Real (*func)(const Prim1DS Wl, const Prim1DS Wr, Real p), const Prim1DS Wl, 
00054            const Prim1DS Wr, Real x1, Real x2, Real xacc);  
00055 
00056 Real getDelVRel(const Prim1DS Wl, const Prim1DS Wr, Real p);
00057 
00058 Real getVRel_2R(const Prim1DS Wl, const Prim1DS Wr, Real p); 
00059 
00060 Real getVRel_RS(const Prim1DS Wl, const Prim1DS Wr, Real p); 
00061 
00062 Real getVRel_2S(const Prim1DS Wl, const Prim1DS Wr, Real p);
00063 
00064 Real getVlim_2S(const Prim1DS Wl, const Prim1DS Wr);
00065 
00066 Real getVlim_RS(const Prim1DS Wl, const Prim1DS Wr); 
00067 
00068 Real getVlim_2R(const Prim1DS Wl, const Prim1DS Wr);
00069 
00070 Real integrateRaref(const Prim1DS W, Real a, Real b);
00071 
00072 Real getP(const Prim1DS Wl, const Prim1DS Wr); 
00073 
00074 void getShockVars(const Prim1DS Wa, Real Pb, char dir, Real *pJ, Real *pV_shock, Real *pd);
00075  
00076 Real getVb_Shock(const Prim1DS Wa, Real Pb, char dir);
00077 
00078 Real getVb_Raref(const Prim1DS Wa, Real Pb, char dir);
00079 
00080 Real getXi(const Prim1DS Wa, Real P, Real vxc, char dir);
00081 
00082 Real rtbis_xi(Real (*func)(const Prim1DS Wa, Real P, Real vx, char), const Prim1DS Wa, 
00083                  char dir, Real f_x1, Real f_x2, Real x1, Real x2, Real xacc);
00084 
00085 void getVelT_Raref(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz);
00086 
00087 void getVelT_Shock(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz);
00088 
00089 void setFluxes(Real vx, Real vy, Real vz, Real P, Real d,Cons1DS *pF );
00090 
00091 /*--------------------------------------------------------------------------*/
00092 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00093  *            const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pF)
00094  *  \brief Computes 1D fluxes using exact special relativistic Riemann solver.
00095  *
00096  *   Input Arguments:
00097  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00098  *   - Bx = B in direction of slice at cell interface
00099  *   Output Arguments:
00100  *   - pF = pointer to fluxes of CONSERVED variables at cell interface
00101  *
00102  *--------------------------------------------------------------------------*/
00103 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00104             const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pF)
00105 {
00106   Real pc, vxc, dcl, dcr; /* pressure, normal velocity, density in center region*/
00107   Real vl_shock, vr_shock;  /* left/right shock velocity */
00108   Real vr_hd, vr_tl, vl_hd, vl_tl;  /* raref head/tail velocities */
00109   Real eps, xacc; 
00110   Real vt, p, vx, vy, vz, C, d; 
00111   Real signY = 1.0; 
00112   Real signZ = 1.0;
00113   Real TOL = 1e-5;
00114  
00115   /* check if intial data states are the same */
00116   if ((fabs(Wl.P-Wr.P) <= TOL)  && (fabs(Wl.Vx-Wr.Vx) <= TOL)) {  
00117     if (vxc >= 0.0) {
00118       setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF);
00119       return; 
00120     }
00121     else {
00122       setFluxes(Wl.Vx, Wr.Vy, Wr.Vz, Wl.P, Wr.d, pF);
00123       return;
00124     }
00125   }
00126 
00127   /* get the pressure in the intermediate state */
00128   pc = getP(Wl, Wr);
00129 
00130   /* determine density and velocity behind waves */
00131   if (wave == RS) {
00132     vxc = getVb_Shock(Wr, pc, 'R');
00133     getShockVars(Wr, pc, 'R', NULL, &vr_shock, &dcr); 
00134     dcl = Wl.d*pow((double)(pc/Wl.P), (double) (1.0/Gamma));
00135   }
00136   else if (wave == SR) {
00137     vxc = getVb_Shock(Wl, pc, 'L'); 
00138     getShockVars(Wl, pc, 'L', NULL, &vl_shock, &dcl); 
00139     dcr = Wr.d*pow((double)(pc/Wr.P), (double) (1.0/Gamma));
00140   }
00141   else if (wave == Two_R) {
00142     vxc = getVb_Raref(Wl, pc, 'L'); 
00143     dcl = Wl.d*pow((double)(pc/Wl.P), (double) (1.0/Gamma));
00144     dcr = Wr.d*pow((double)(pc/Wr.P), (double) (1.0/Gamma));
00145   }
00146   else if (wave == Two_S) {
00147     vxc = getVb_Shock(Wl, pc, 'L');
00148     getShockVars(Wl, pc, 'L', NULL, &vl_shock, &dcl); 
00149     getShockVars(Wr, pc, 'R', NULL, &vr_shock, &dcr);    
00150   }
00151   
00152 /*-----------------------------------------------------------------
00153 * Calculate the interface flux if the wave speeds are such that we aren't
00154 * actually in the intermediate state */
00155   
00156   if (pc > Wl.P) {
00157     /* left shock wave */
00158 
00159     if (vl_shock >= 0.0) {
00160       /* to left of shock */   
00161       setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF);
00162       return; 
00163     }
00164   }
00165   else {
00166     /* left rarefaction */
00167     
00168    /* calculate velocity at head and tail of raref */
00169     vl_hd = getXi(Wl, Wl.P, Wl.Vx, 'L');  
00170     vl_tl = getXi(Wl, pc, vxc, 'L'); 
00171     
00172     if (vl_hd >= 0.0) {
00173       /* To left of rarefaction */
00174       setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF); 
00175       return;
00176     } 
00177     else if (vl_tl >= 0.0) {
00178       /* Inside rarefaction fan */
00179       
00180       /* machine precision */ 
00181       eps = 1.0f; 
00182       do eps /= 2.0f; 
00183       while ((float)(1.0 + (eps/2.0)) != 1.0);
00184       xacc = eps*0.5*(Wl.P + pc); 
00185 
00186       /* get vx, p, d in raref fan */
00187       p = rtbis_xi(getXi, Wl, 'L', vl_hd, vl_tl, Wl.P, pc, xacc); 
00188       vx = getVb_Raref(Wl, p, 'L'); 
00189       d = Wl.d*pow((double)(p/Wl.P), (double) (1.0/Gamma));
00190      
00191       getVelT_Raref(Wl, p, vx, &vy, &vz);     
00192       setFluxes(vx, vy, vz, p, d, pF); 
00193       return; 
00194     } 
00195   }
00196 
00197   if (pc > Wr.P) {
00198     /* right shock wave */
00199     
00200     /* right shock speed */
00201    
00202     if (vr_shock <= 0.0) {
00203       /* to right of shock */
00204 
00205       setFluxes(Wr.Vx, Wr.Vy, Wr.Vz, Wr.P, Wr.d, pF); 
00206       return; 
00207     }
00208   }
00209   else {
00210     /* right rarefaction */
00211    
00212     /* calculate velocity at head and tail of raref */
00213     vr_hd = getXi(Wr, Wr.P, Wr.Vx, 'R');  
00214     vr_tl = getXi(Wr, pc, vxc, 'R'); 
00215     
00216     if (vr_hd <= 0.0) {
00217       /* To right of rarefaction */
00218       
00219       setFluxes(Wr.Vx, Wr.Vy, Wr.Vz, Wr.P, Wr.d, pF); 
00220       return;
00221       
00222     }
00223     else if (vr_tl <= 0.0) {
00224       /* Inside rarefaction fan */
00225       
00226       /* machine precision */ 
00227       eps = 1.0f; 
00228       do eps /= 2.0f; 
00229       while ((float)(1.0 + (eps/2.0)) != 1.0); 
00230       xacc = eps*0.5*(Wr.P + pc); 
00231 
00232       /* get vx, p, d in raref fan */
00233       p = rtbis_xi(getXi, Wr, 'R', vr_hd, vr_tl, Wr.P, pc, xacc); 
00234       vx = getVb_Raref(Wr, p, 'R'); 
00235       d = Wr.d*pow((double)(p/Wr.P), (double) (1.0/Gamma));
00236      
00237       getVelT_Raref(Wr, p, vx, &vy, &vz);  
00238       setFluxes(vx, vy, vz, p, d, pF); 
00239       return;
00240     }
00241   }
00242 
00243 /* We are in the intermediate state */
00244 /*---------------------------------------------------------------------
00245  * Calculate the interface flux */
00246 
00247   if (vxc >= 0.0) {
00248     
00249     if (Wl.P > pc) {
00250       /* left rarefaction */
00251       getVelT_Raref(Wl, pc, vxc, &vy, &vz);
00252       setFluxes(vxc, vy, vz, pc, dcl, pF); 
00253       return;
00254     }
00255     else {
00256       /*left shock */
00257       getVelT_Shock(Wl, pc, vxc, &vy, &vz);
00258       setFluxes(vxc, vy, vz, pc, dcl, pF); 
00259       return;
00260     }
00261    
00262   }
00263   else {
00264     
00265     if (pc > Wr.P) {
00266       /* right shock */
00267       getVelT_Shock(Wr, pc, vxc, &vy, &vz);    
00268       setFluxes(vxc, vy, vz, pc, dcr, pF); 
00269       return;
00270     }
00271     else {
00272       /* right rarefaction */
00273       
00274       getVelT_Raref(Wr, pc, vxc, &vy, &vz);
00275       setFluxes(vxc, vy, vz, pc, dcr, pF); 
00276       return;
00277     }
00278   }
00279 }
00280 /*----------------------------------------------------------------------------*/
00281 /*! \fn Real integrateRaref(const Prim1DS Wa, Real a, Real b)
00282  *  \brief Numerically integrate the integral for rarefactions (eq 3.22 in  RZP)
00283  * using Gaussian-Legendre quadrature */
00284 Real integrateRaref(const Prim1DS Wa, Real a, Real b)
00285 {
00286   
00287   int i;
00288   Real integral, sum, diff, f, xx; 
00289   Real * x; 
00290   Real * w; 
00291   int n = 10; /* number of integration points */
00292   Real k, h, vt, G, A;
00293   Real dd, ccs2, hh;
00294 
00295 
00296   /* allocate memory for absicssas and weights */
00297   x = (Real *)malloc((size_t) ((n+1)*sizeof(Real)));
00298   w = (Real *)malloc((size_t) ((n+1)*sizeof(Real)));
00299  
00300   /* get the abscissas and weights for integration */
00301   gauleg(-1.0, 1.0, x, w, n);
00302 
00303   sum = 0.5 * (b + a); 
00304   diff = 0.5 * (b - a); 
00305   integral = 0.0; 
00306 
00307   for (i = 1; i <= n; i++) {
00308     xx = diff*x[i] + sum; 
00309 
00310     h = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d); 
00311     vt = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz)); 
00312   
00313     G = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00314     A = h*G*vt;
00315 
00316     dd = Wa.d*pow((double)(xx/Wa.P), (double) (1.0/Gamma)); 
00317     ccs2 = Gamma*(Gamma - 1.)*xx / (Gamma*xx + (Gamma - 1.)*dd); 
00318     hh = 1.0 + xx*Gamma/(dd*(Gamma - 1.0)); 
00319   
00320     f = sqrt((double) (hh*hh + A*A*(1.0-ccs2)))/(dd*sqrt((double)ccs2)*(hh*hh + A*A));
00321 
00322     integral = integral + diff*w[i]*f; 
00323   }
00324   
00325   free(x); 
00326   free(w); 
00327 
00328   return integral; 
00329 }
00330 
00331 /*----------------------------------------------------------------------------*/
00332 /*! \fn void getShockVars(const Prim1DS Wa, Real Pb, char dir, Real *pJ, 
00333  *                        Real *pV_s, Real *pD)
00334  *  \brief Determine the mass flux, shock velocity and density across a shock 
00335  *  wave */
00336 void getShockVars(const Prim1DS Wa, Real Pb, char dir, Real *pJ, Real *pV_s, 
00337                   Real *pD)
00338 {
00339   Real sign, ha, Ga;
00340   Real db, hb;  /* specific enthalpy and density behind wave */
00341   Real A, B, C,D; 
00342   Real d, J, v_s;
00343   Real TOL = 1e-5;
00344 
00345   if (dir == 'L') sign = -1.0; 
00346   if (dir == 'R') sign = 1.0;
00347 
00348   /* check if pressure is constant across wave */
00349   if (fabs(Wa.P - Pb) <= TOL) { 
00350     if (pJ != NULL)
00351       *pJ = 0.0; 
00352     if (pD != NULL)
00353       *pD = Wa.d;
00354     if (pV_s != NULL)
00355       *pV_s = Wa.Vx;
00356     
00357     return;
00358   }
00359   
00360   ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00361   Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00362 
00363   /* calculate specific enthalpy */ 
00364   A = 1.0 + (Gamma - 1.0)*(Wa.P-Pb)/(Gamma*Pb); 
00365   B = 1.0 - A;
00366   C = ha*(Wa.P - Pb)/Wa.d - ha*ha; 
00367 
00368   /*check for unphysical specific enthalpies*/
00369   if (C > (B*B/(4.0*A)))
00370     ath_error("[exact flux]: Unphysical specific enthalpy in intermediate state");
00371     
00372   hb = (-B + sqrt((double)(B*B - 4.0*A*C)))/(2.0*A); 
00373 
00374   d = Gamma*Pb/((Gamma - 1.0)*(hb - 1.0));  
00375   if (pD != NULL)
00376     *pD = d; 
00377   
00378   J =  sign*sqrt((double)((Pb - Wa.P)/(ha/Wa.d - hb/d)));
00379   if (pJ != NULL)
00380     *pJ = J;
00381  
00382   A = Wa.d*Wa.d*Ga*Ga;
00383   if (pV_s != NULL) 
00384     *pV_s = (A*Wa.Vx +sign*fabs((float) J)*
00385              sqrt((double)(J*J + A*(1.0 - Wa.Vx*Wa.Vx))))/(A + J*J);
00386   
00387   return;
00388 }
00389 
00390 /*----------------------------------------------------------------------------*/
00391 /*! \fn Real getVlim_2R(const Prim1DS Wl, const Prim1DS Wr)
00392  *  \brief Determine vlim for two rarefactions (eq. 4.12 in RZP) */
00393 Real getVlim_2R(const Prim1DS Wl, const Prim1DS Wr)
00394 {
00395   
00396   Real vxlc, vxrc;  
00397   Real integral; 
00398   Real vlb, vrb; 
00399 
00400   /* compute vxl,c from RZP eq 4.13 */
00401   vlb = getVb_Raref(Wl, 0.0, 'L');
00402   vxlc = (Wl.Vx - vlb)/(1.0 - Wl.Vx*vlb);
00403 
00404   /* compute vxr,c from RZP eq 4.14 */
00405   vrb = getVb_Raref(Wr, 0.0, 'R');
00406   vxrc = (Wr.Vx - vrb)/(1.0 - Wr.Vx*vrb);
00407 
00408   return (vxlc - vxrc) / (1.0 - vxlc*vxrc); 
00409   
00410 }
00411 
00412 /*----------------------------------------------------------------------------*/
00413 /*! \fn Real getVlim_RS(const Prim1DS Wl, const Prim1DS Wr)
00414  *  \brief Determine vlim for the case of rarefaction-shock (eq 4.10) */
00415 Real getVlim_RS(const Prim1DS Wl, const Prim1DS Wr)
00416 {
00417  
00418   Real a, b, integral; 
00419   Real vlb, vrb, vlc, vrc; 
00420   
00421   if (Wl.P < Wr.P) { 
00422     /* shock to left and rarefaction to right */
00423     /* limiting p in center is Wl.P */
00424     vlb = Wl.Vx;
00425     vrb = getVb_Raref(Wr, Wl.P, 'R');
00426   }
00427   else {
00428     /* shock to right and rarefaction to left */
00429     /* limiting p in center is Wr.P */
00430     vlb = getVb_Raref(Wl, Wr.P, 'L');
00431     vrb = Wr.Vx;
00432   }
00433      
00434   vlc = (Wl.Vx - vlb)/(1.0 - Wl.Vx*vlb);
00435   vrc = (Wr.Vx - vrb)/(1.0 - Wr.Vx*vrb);
00436     
00437   return (vlc - vrc)/(1.0 - vlc*vrc);
00438 } 
00439 
00440 /*--------------------------------------------------------------------------*/
00441 /*! \fn Real getVlim_2S(const Prim1DS Wl, const Prim1DS Wr)
00442  *  \brief Determine vlim for the case of two shocks (eq 4.5) */
00443 Real getVlim_2S(const Prim1DS Wl, const Prim1DS Wr)
00444 {  
00445   if (Wl.P > Wr.P) {
00446     Real vlc, vrc, vrb; 
00447 
00448     vrb = getVb_Shock(Wr, Wl.P, 'R');  
00449     vrc = (Wr.Vx - vrb)/(1.0 - vrb*Wr.Vx); 
00450     vlc = 0.0; 
00451     return (vlc - vrc)/(1.0 - vlc*vrc);  
00452   }
00453   else {
00454     Real vlc, vrc, vlb; 
00455 
00456     vlb = getVb_Shock(Wl, Wr.P, 'L'); 
00457     vlc = (Wl.Vx - vlb)/(1.0 - vlb*Wl.Vx); 
00458     vrc = 0.0; 
00459     return (vlc - vrc)/(1.0 - vlc*vrc);
00460     
00461   }
00462 }
00463 
00464 /*----------------------------------------------------------------------------*/       
00465 /*! \fn Real getDelVRel(const Prim1DS Wl, const Prim1DS Wr, Real p)
00466  *  \brief Determine vRel for the wave pattern indicated by wave */
00467 Real getDelVRel(const Prim1DS Wl, const Prim1DS Wr, Real p)
00468 {
00469   Real vRel; 
00470   Real vRel_0 = (Wl.Vx - Wr.Vx)/(1.0 - Wl.Vx*Wr.Vx); 
00471   
00472   if (wave == Two_R) {
00473     vRel = getVRel_2R(Wl, Wr, p); 
00474     return vRel - vRel_0; 
00475   }
00476   
00477   else if (wave == RS || wave == SR) {
00478     vRel = getVRel_RS(Wl, Wr, p);
00479     return vRel - vRel_0; 
00480   }
00481   else if (wave == Two_S) {
00482     vRel = getVRel_2S(Wl, Wr, p);
00483     return vRel - vRel_0; 
00484   }
00485 
00486   return -1.0; /* should never reach here*/
00487 }
00488 /*----------------------------------------------------------------------------*/
00489 /*! \fn Real getP(const Prim1DS Wl, const Prim1DS Wr)
00490  *  \brief Determine the pressure in the intermediate state using interval 
00491  *  bisection */
00492 Real getP(const Prim1DS Wl, const Prim1DS Wr)
00493 {
00494   Real vRel, vRS, vSS, vRR; 
00495   Real pmin, pmax;  /* brackets for pressure */ 
00496   Real xacc; 
00497   Real pc;  /* pressure in intermediate state*/
00498   float tol, eps; 
00499 
00500   /* 1. Compute the relative x-velocity between Wl and Wr */
00501   vRel = (Wl.Vx - Wr.Vx) / (1.0 - Wl.Vx*Wr.Vx); 
00502   
00503   /* 2. Determine the wave pattern by comparing vRel with the limiting values of
00504      the relative velocity and bracket pressure in intermed state */
00505   vSS = getVlim_2S(Wl, Wr); 
00506   vRS = getVlim_RS(Wl, Wr); 
00507 
00508   if (vRel <= vRS) {
00509     wave = Two_R;
00510     pmin = 0.0; 
00511     if (Wl.P < Wr.P) pmax = Wl.P; 
00512     else pmax = Wr.P;
00513   }
00514   else if (vRel > vRS && vRel <= vSS) {
00515     if (Wl.P > Wr.P) {
00516       wave = RS; 
00517       pmin = Wr.P;
00518       pmax = Wl.P;
00519     }
00520     else { 
00521       wave = SR; 
00522       pmin = Wl.P; 
00523       pmax = Wr.P;
00524     } 
00525   }
00526   else {
00527     wave = Two_S; 
00528     if (Wl.P > Wr.P) pmin = Wl.P;
00529     else pmin = Wr.P;
00530     pmax = 1000 * 0.5*(Wl.P + Wr.P); 
00531   }
00532 
00533   /* 3. compute the machine precision */
00534   eps = 1.0f; 
00535   do eps /= 2.0f; 
00536   while ((float)(1.0 + (eps/2.0)) != 1.0); 
00537   xacc = eps*0.5*(pmin + pmax); 
00538 
00539   /* 4. use interval bisection to find p in the intermediate state */
00540   pc = rtbis_vel(getDelVRel, Wl, Wr, pmin, pmax, xacc); 
00541   return pc; 
00542   
00543 }
00544 
00545 /*---------------------------------------------------------------------------*/
00546 /*! \fn Real getVRel_2R(const Prim1DS Wl, const Prim1DS Wr, Real p)
00547  *  \brief Determine vRel for the case of two rarefactions */
00548 Real getVRel_2R(const Prim1DS Wl, const Prim1DS Wr, Real p)
00549 {
00550 
00551   Real vlb, vrb; /* the normal velocity behind the left and right rarefactions */
00552   Real integral; 
00553   Real vlc, vrc; /* normal velocity of left/right waves in frame of contact */
00554 
00555   /* determine x-velocity behind left/right waves */
00556   vlb = getVb_Raref(Wl, p, 'L'); 
00557  
00558   vrb = getVb_Raref(Wr, p, 'R'); 
00559 
00560   /* determine velocities in frame of contact discontinuity */ 
00561   vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb); 
00562   vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb); 
00563 
00564   return (vlc - vrc)/(1.0 - vrc*vlc); 
00565 }
00566 
00567 /*---------------------------------------------------------------------------*/
00568 /*! \fn Real getVRel_RS(const Prim1DS Wl, const Prim1DS Wr, Real p)
00569  *  \brief Determine vRel for the case of rarefaction-shock wave */
00570 Real getVRel_RS(const Prim1DS Wl, const Prim1DS Wr, Real p)
00571 {
00572   Real vlb, vrb; /* normal velocity behind the waves */
00573   Real vlc, vrc; /* normal velocity in frame of contact */
00574   Real v_shock, G_shock, J, ha, Ga; 
00575   Real num, denom, integral; 
00576 
00577   if (wave == SR) {
00578     /* shock wave moving to left and rarefaction to right */
00579     
00580     /* calculate vlb for shock*/
00581     vlb = getVb_Shock(Wl, p, 'L'); 
00582 
00583     /*calculate vrb for rarefaction */
00584     vrb = getVb_Raref(Wr, p, 'R');
00585 
00586     vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb); 
00587     vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb); 
00588 
00589     return (vlc - vrc)/(1.0 - vrc*vlc); 
00590   }
00591   else {
00592     /* shock wave moving to right and rarefaction to left */
00593     
00594     /* calculate vrb for shock  */
00595     vrb = getVb_Shock(Wr, p, 'R'); 
00596 
00597     /*calculate vlb for raref*/
00598     vlb = getVb_Raref(Wl, p, 'L'); 
00599 
00600     vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb); 
00601     vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb); 
00602 
00603     return (vlc - vrc)/(1.0 - vrc*vlc); 
00604   } 
00605 }
00606 /*---------------------------------------------------------------------------*/
00607 /*! \fn Real getVRel_2S(const Prim1DS Wl, const Prim1DS Wr, Real p)
00608  *  \brief Determine vRel for the case of two shocks */
00609 Real getVRel_2S(const Prim1DS Wl, const Prim1DS Wr, Real p)
00610 {
00611   Real vlb, vrb; /* normal velocity behind the shock waves */
00612   Real vlc, vrc; /* normal velocity in frame of contact */
00613   char dir;
00614   Real v_shock, G_shock, J, ha, Ga; 
00615   Real num, denom; 
00616 
00617   /* calculate vlb */
00618   vlb = getVb_Shock(Wl, p, 'L');
00619 
00620   /* calculate vrb */ 
00621   vrb = getVb_Shock(Wr, p, 'R'); 
00622 
00623   /* determine velocities in frame of contact discontinuity */ 
00624   vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb); 
00625   vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb);
00626 
00627   return (vlc - vrc)/(1.0 - vlc*vrc); 
00628 }
00629 /*---------------------------------------------------------------------------*/
00630 /*! \fn Real getVb_Shock(const Prim1DS Wa, Real Pb, char dir) 
00631  *  \brief Determine the normal velocity behind a shock wave */
00632 Real getVb_Shock(const Prim1DS Wa, Real Pb, char dir)
00633 {
00634   Real num, denom;
00635   Real G_shock, ha, Ga;
00636   Real J, v_shock; 
00637   
00638   getShockVars(Wa, Pb, dir, &J, &v_shock, NULL); 
00639   if (fabs(J) == 0.0) return Wa.Vx;
00640 
00641   G_shock = 1.0 / sqrt((double) (1.0 - v_shock*v_shock));
00642   
00643   /* specific enthalpy and gamma factor for left state */
00644   ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00645   Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00646 
00647   num = ha*Ga*Wa.Vx + G_shock*(Pb - Wa.P)/J; 
00648   denom = ha*Ga + (Pb - Wa.P)*(G_shock*Wa.Vx/J + 1.0/(Wa.d*Ga));  
00649 
00650   return num/denom; 
00651 }
00652 /*----------------------------------------------------------------------------*/
00653 /*! \fn Real getVb_Raref(const Prim1DS Wa, Real Pb, char dir)
00654  *  \brief Determine the normal velocity behind a rarefaction wave */
00655 Real getVb_Raref(const Prim1DS Wa, Real Pb, char dir)
00656 {
00657   Real integral, sign, B; 
00658   
00659   if (dir == 'L') sign = -1.0; 
00660   if (dir == 'R') sign = 1.0; 
00661   
00662   integral = integrateRaref(Wa, Wa.P, Pb);
00663   B = 0.5*log((double)((1. + Wa.Vx)/(1. - Wa.Vx))) + sign*integral; 
00664               
00665   return tanh(B); 
00666 }
00667 /*----------------------------------------------------------------------------*/
00668 /*! \fn Real getXi(const Prim1DS Wa, Real P, Real vxc, char dir) 
00669  *  \brief Determine the self-similarity variable xi for rarefactions (eq 3.15) 
00670 */
00671 Real getXi(const Prim1DS Wa, Real P, Real vxc, char dir) 
00672 {
00673   Real vt, h, G, A; /* vars associated with state ahead of
00674                       rarefaction wave */
00675   Real dc, hc, vtc, cs, v2, num, denom; /* vars inside wave */
00676   Real sign;  
00677 
00678   if (dir == 'L') sign = -1.0; 
00679   if (dir == 'R') sign = 1.0; 
00680   
00681   vt = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz)); 
00682   h = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);  
00683   G = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00684   A = h*G*vt;
00685 
00686   dc = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00687   hc =  1.0 + Gamma*P/((Gamma-1.0)*dc);  
00688   vtc = A*sqrt((double)((1.0 - vxc*vxc)/(hc*hc + A*A)));
00689   cs = sqrt((double)(Gamma*(Gamma-1.0)*P/((Gamma -1.0)*dc + Gamma*P)));
00690 
00691   v2 = vxc*vxc + vtc*vtc; 
00692   num = vxc*(1-cs*cs) + sign*cs*sqrt((double)((1-v2)*(1.-v2*cs*cs-vxc*vxc*(1.-cs*cs))));
00693   denom = 1.0 - v2*cs*cs; 
00694   
00695   return num/denom; 
00696 }
00697 /*---------------------------------------------------------------------------*/ 
00698 /*! \fn void getVelT_Raref(const Prim1DS Wa, Real P, Real vxb, Real *pVy, 
00699  *                         Real *pVz) 
00700  *  \brief Determine the tangential velocities vy and vz behind a rarefaction 
00701  *  wave */
00702 void getVelT_Raref(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz) 
00703 {
00704   Real va_t, ha, Ga, A; 
00705   Real db, hb, vb_t, Gb, C;
00706   Real signY = 1.0; 
00707   Real signZ = 1.0;
00708 
00709   if (Wa.Vy == 0.0 && Wa.Vz == 0.0) { 
00710     *pVy = 0.0;  
00711     *pVz = 0.0;  
00712     return;
00713   }
00714   
00715   va_t = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz)); 
00716   ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);  
00717   Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00718   A = ha*Ga*va_t;
00719 
00720   db = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00721   hb =  1.0 + Gamma*P/((Gamma-1.0)*db);  
00722 
00723   vb_t = A*sqrt((double)((1.0 - vxb*vxb)/(hb*hb + A*A)));
00724   Gb = 1.0/sqrt((double)(1.0 -vxb*vxb - vb_t*vb_t));
00725   
00726  
00727   if (Wa.Vy == 0.0) { 
00728     if (Wa.Vz < 0.0) signZ = -1.0; 
00729     *pVy = 0.0;
00730     *pVz = signZ*vb_t;
00731   }
00732   else if (Wa.Vz == 0.0) {
00733     if (Wa.Vy < 0.0) signY = -1.0; 
00734       *pVz = 0.0;
00735       *pVy = signY*vb_t;
00736   }
00737   else {
00738     C = (Wa.Vy*Wa.Vy) / (Wa.Vz*Wa.Vz);
00739     if (Wa.Vy < 0.0) signY = -1.0;
00740     if (Wa.Vz < 0.0) signZ = -1.0; 
00741     
00742     *pVz = signZ*sqrt((double)(vb_t*vb_t / (1.0 + C))); 
00743     *pVy = signY*sqrt((double)(vb_t*vb_t / (1.0 + 1.0/C))); 
00744   }
00745 }
00746 /*----------------------------------------------------------------------------*/
00747 /*! \fn void getVelT_Shock(const Prim1DS Wa, Real P, Real vxb, Real *pVy, 
00748  *                         Real *pVz)
00749  *  \brief Determine the tangential velocitiies vy and vz behind a shock wave. 
00750  *
00751  *  Note that this function differs from Pons et al. (2000) equation 4.13, 
00752  *  which is only valid for one of vz or vy equal to zero. */
00753 void getVelT_Shock(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz)
00754 {
00755   Real va_t, ha, Ga, Ay, Az; 
00756   Real hb, db, vy, vz, Gb;
00757   Real Cy, Cz; 
00758   Real D; 
00759 
00760   if (Wa.Vy == 0.0 && Wa.Vz == 0.0) { 
00761     *pVy = 0.0;  
00762     *pVz = 0.0;  
00763     return;
00764   }
00765 
00766   ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);  
00767   Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz)); 
00768   
00769   Ay = ha*Ga*Wa.Vy;
00770   Az = ha*Ga*Wa.Vz; 
00771   db = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00772   hb =  1.0 + Gamma*P/((Gamma-1.0)*db);  
00773 
00774   Cz = Az*Az/(hb*hb+Az*Az);
00775   Cy = Ay*Ay/(hb*hb+Ay*Ay);
00776   D = (1.0 - Cz*Cy);
00777 
00778   vz = sqrt((double)(Cz*(1-vxb*vxb)*(1-Cy)/D));
00779   vy = sqrt((double)(Cy*(1-vxb*vxb)*(1-Cz)/D));
00780   
00781   Gb = 1.0/sqrt((double)(1.0 - vxb*vxb - vy*vy - vz*vz));
00782 
00783   if (Wa.Vy >= 0.0) 
00784     *pVy = vy; 
00785   else 
00786     *pVy = -vy;
00787   
00788   if (Wa.Vz >= 0.0) 
00789     *pVz = vz; 
00790   else
00791     *pVz = -vz;
00792 }
00793 /*---------------------------------------------------------------------------*/ 
00794 /*! \fn void setFluxes(Real vx, Real vy, Real vz, Real P, Real d, Cons1DS *pF)
00795  *  \brief Set the corresponding fluxes as the fields of pF */
00796 void setFluxes(Real vx, Real vy, Real vz, Real P, Real d, Cons1DS *pF)
00797 {
00798   Real G, h, D, Sx, Sy, Sz; 
00799   
00800   if ((vx*vx + vy*vy + vz*vz) >= 1.0)
00801     ath_error("[exact_flux]: Superluminal velocities vx = %f, vy = %f, vz = %f\n",
00802               vx, vy, vz);
00803 
00804   G = 1.0/sqrt((double)(1.0 - vx*vx - vy*vy - vz*vz));
00805   h = 1.0 + Gamma*P/((Gamma-1.0)*d); 
00806   D = d*G; 
00807   Sx = d*h*G*G*vx;
00808   Sy = d*h*G*G*vy;
00809   Sz = d*h*G*G*vz; 
00810 
00811   pF->d  = D*vx;
00812   pF->Mx = Sx*vx + P;
00813   pF->My = Sy*vx;
00814   pF->Mz = Sz*vx;
00815   pF->E = Sx;
00816 }
00817 /*----------------------------------------------------------------------------*/
00818 /*! \fn void gauleg(double x1, double x2, Real x[], Real w[], int n)
00819  *  \brief Given the lower and upper limits of integration x1 and x2, and 
00820  *   the abscissas and weights of the Gauss-Legendre n-point quadrature formula.
00821  *  
00822  *  Given the lower and upper limits of integration x1 and x2, and given n, 
00823  *   this routine returns arrays x[1..n] and w[1..n] of length n, containing 
00824  *   the abscissas and weights of the Gauss-Legendre n-point quadrature formula.
00825  *   See Numerical Recipes in C, Press et al. */
00826 void gauleg(double x1, double x2, Real x[], Real w[], int n)
00827 {
00828   int m,j,i;
00829   Real z1,z,xm,xl,pp,p3,p2,p1;
00830   m=(n+1)/2; /* The roots are symmetric, so we only find half of them. */
00831   xm=0.5*(x2+x1);
00832   xl=0.5*(x2-x1);
00833   
00834   for (i=1; i<=m; i++) { /* Loop over the desired roots. */
00835     z=cos(3.141592654*(i-0.25)/(n+0.5));
00836     /* Starting with the above approximation to the ith root, we enter */
00837     /* the main loop of refinement by Newton's method.                 */
00838     do {
00839       p1=1.0;
00840       p2=0.0;
00841       for (j=1;j<=n;j++) { /* Recurrence to get Legendre polynomial. */
00842         p3=p2;
00843         p2=p1;
00844         p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
00845       }
00846       /* p1 is now the desired Legendre polynomial. We next compute */
00847       /* pp, its derivative, by a standard relation involving also  */
00848       /* p2, the polynomial of one lower order.                     */
00849       pp=n*(z*p1-p2)/(z*z-1.0);
00850       z1=z;
00851       z=z1-p1/pp; /* Newton's method. */
00852     } while (fabs(z-z1) > EPS);
00853     x[i]=xm-xl*z;      /* Scale the root to the desired interval, */
00854     x[n+1-i]=xm+xl*z;  /* and put in its symmetric counterpart.   */
00855     w[i]=2.0*xl/((1.0-z*z)*pp*pp); /* Compute the weight             */
00856     w[n+1-i]=w[i];                 /* and its symmetric counterpart. */
00857   }
00858 }
00859 
00860 /*-------------------------------------------------------------------------*/
00861 
00862 /*! \fn Real rtbis_vel(Real (*func)(const Prim1DS Wl, const Prim1DS Wr, Real p),
00863  *             const Prim1DS Wl, const Prim1DS Wr, Real x1, Real x2, Real xacc)
00864  *  \brief Implements bisection root finding algorithm to solve the velocity 
00865  *  eqn. 
00866  *
00867  * Assumes func(x1) and func(x2) have opposite signs without a check */
00868 Real rtbis_vel(Real (*func)(const Prim1DS Wl, const Prim1DS Wr, Real p),
00869                const Prim1DS Wl, const Prim1DS Wr, Real x1, Real x2, Real xacc)
00870 {
00871   Real j; 
00872   Real dx, f, fmid, xmid, rtb; 
00873   Real vRel_0 = (Wl.Vx-Wr.Vx)/(1.0-Wl.Vx*Wr.Vx);
00874 
00875   if (wave == RS || wave == SR) {
00876     f = getVlim_RS(Wl, Wr) - vRel_0;
00877     fmid = getVlim_2S(Wl, Wr) - vRel_0; 
00878   }
00879  
00880   if (wave == Two_R) {
00881     f = getVlim_2R(Wl, Wr) - vRel_0; 
00882     fmid = getVlim_RS(Wl, Wr) - vRel_0; 
00883   }
00884   
00885   if (wave == Two_S) {
00886     f = getVlim_2S(Wl, Wr) - vRel_0; 
00887     fmid = (*func)(Wl, Wr, x2); 
00888   }
00889   if (f < 0.0) {
00890     dx = x2-x1;
00891     rtb =x1; 
00892   }
00893   else {
00894     dx = x1-x2; 
00895     rtb = x2; 
00896   }
00897   
00898   for (j = 1; j <= JMAX; j++) {
00899     dx = dx / 2.0; 
00900     xmid = rtb + dx;
00901     fmid = (*func)(Wl, Wr, xmid);  
00902     if (fmid <= 0.0) rtb = xmid;
00903     if (fabs(dx) < xacc || fmid == 0.0) {
00904       return rtb;
00905     }
00906   }   
00907 }
00908 /*----------------------------------------------------------------------------*/
00909 /*! \fn Real rtbis_xi(Real (*func)(const Prim1DS Wa, Real P, Real vx, char), 
00910  *                    const Prim1DS Wa, char dir, Real f_x1, Real f_x2, Real x1,
00911  *                    Real x2, Real xacc)
00912  *  \brief Implements bisection root finding algorithm to find the rarefaction 
00913  *  head and tail velocity. 
00914  *  Assumes func(x1) and func(x2) have opposite signs without a check */ 
00915    
00916 Real rtbis_xi(Real (*func)(const Prim1DS Wa, Real P, Real vx, char), 
00917               const Prim1DS Wa, char dir, Real f_x1, Real f_x2, Real x1, 
00918               Real x2, Real xacc) 
00919 {
00920   Real j, vx; 
00921   Real dx, f, fmid, xmid, rtb; 
00922 
00923   f = f_x1; 
00924   fmid = f_x2; 
00925 
00926   if (f < 0.0) {
00927     dx = x2-x1;
00928     rtb =x1; 
00929   }
00930   else {
00931     dx = x1-x2; 
00932     rtb = x2; 
00933   }
00934   
00935   for (j = 1; j <= JMAX; j++) {
00936     dx = dx / 2.0; 
00937     xmid = rtb + dx;
00938     vx = getVb_Raref(Wa, xmid, dir); 
00939     fmid = (*func)(Wa, xmid, vx, dir);
00940  
00941     
00942     if (fmid <= 0.0) rtb = xmid;
00943     if (fabs(dx) < xacc || fmid == 0.0) {
00944       return rtb;
00945     }
00946   }
00947   
00948  
00949 }
00950 #endif /* Special Relativity */
00951 #endif /* Exact flux */