Line data    Source code

       1              : #ifndef SOLVER_LOCAL_RIEMANN_HLLE_H
       2              : #define SOLVER_LOCAL_RIEMANN_HLLE_H
       3              : 
       4              : #include "IO/ParmParse.H"
       5              : #include "Solver/Local/Riemann/Riemann.H"
       6              : 
       7              : /// A bunch of solvers
       8              : namespace Solver
       9              : {
      10              : /// Local solvers
      11              : namespace Local
      12              : {
      13              : 
      14              : namespace Riemann
      15              : {
      16              : 
      17              : class HLLE : public Riemann
      18              : {
      19              : public:
      20              : 
      21              : 
      22              :     static constexpr const char* name = "hlle";
      23            0 :     HLLE (IO::ParmParse &pp, std::string name) 
      24            0 :     {pp_queryclass(name,*this);}
      25              :     HLLE (IO::ParmParse &pp) 
      26              :     {pp_queryclass(*this);}
      27              :     HLLE () 
      28              :     {
      29              :         IO::ParmParse pp;
      30              :         pp_queryclass(*this);
      31              :     }
      32              : 
      33              :     int lowmach = 0;
      34              :     Set::Scalar cutoff = NAN;
      35              : 
      36            0 :     static void Parse(HLLE & value, IO::ParmParse & pp)
      37              :     {
      38              :         // toggle to apply a low mach correction approximation
      39            0 :         pp.query_default("lowmach",value.lowmach,false);
      40              :         // cutoff value if using the low mach approximation
      41            0 :         pp.query_default("cutoff",value.cutoff,0.1);
      42            0 :     }
      43              : 
      44            0 :     virtual Flux Solve(State lo, State hi, Set::Scalar gamma, Set::Scalar p_ref, Set::Scalar small) override
      45              :     {
      46              :         // Grab and label the input states
      47            0 :         Set::Scalar rho_L = lo.rho       ,  rho_R = hi.rho;
      48            0 :         Set::Scalar Mn_L  = lo.M_normal  ,  Mn_R  = hi.M_normal  ;
      49              :         //Set::Scalar Mt_L  = lo.M_tangent ,  Mt_R  = hi.M_tangent ;
      50            0 :         Set::Scalar E_L   = lo.E         ,  E_R   = hi.E         ;
      51              : 
      52              :         // Ensure no negative densities
      53            0 :         rho_L = std::max(0.0,rho_L);
      54            0 :         rho_R = std::max(0.0,rho_R);
      55              : 
      56              : 
      57              :         // Calculate primitives
      58            0 :         Set::Scalar u_L   = Mn_L/(rho_L+small),  u_R   = Mn_R/(rho_R+small);
      59            0 :         Set::Scalar v_L   = Mn_L/(rho_L+small),  v_R   = Mn_R/(rho_R+small);
      60            0 :         Set::Scalar p_L  = (gamma - 1.0) * ( E_L - 0.5*rho_L*(u_L*u_L + v_L*v_L) ) + p_ref;
      61            0 :         Set::Scalar p_R  = (gamma - 1.0) * ( E_R - 0.5*rho_R*(u_R*u_R + v_R*v_R) ) + p_ref;
      62            0 :         Set::Scalar c_L = std::sqrt(gamma * p_L / (rho_L + small));
      63            0 :         Set::Scalar c_R = std::sqrt(gamma * p_R / (rho_R + small));
      64              : 
      65              :         // Max wave speeds
      66            0 :         Set::Scalar S_L = std::min(u_L-c_L, u_R - c_R);
      67            0 :         Set::Scalar S_R = std::max(u_L+c_L, u_R + c_R);
      68              : 
      69              :         
      70            0 :         if (lowmach)
      71              :         {
      72            0 :             Set::Scalar u_RL = 0.5 * (u_L + u_R);
      73            0 :             Set::Scalar a_RL = 0.5 * (c_L + c_R);
      74            0 :             Set::Scalar Mach = std::abs(u_RL) / (a_RL + small);
      75            0 :             Set::Scalar psi = std::max(Mach, cutoff); // Mach cutoff
      76            0 :             S_L *= psi;
      77            0 :             S_R *= psi;
      78              :         }
      79              : 
      80              :         Flux F_L (
      81              :             lo.M_normal,
      82            0 :             lo.M_normal * u_L + p_L,
      83            0 :             lo.M_tangent * u_L,
      84            0 :             u_L * (lo.E + p_L)
      85            0 :             );
      86              :         Flux F_R (
      87              :             hi.M_normal,
      88            0 :             hi.M_normal * u_R + p_R,
      89            0 :             hi.M_tangent * u_R,
      90            0 :             u_R * (hi.E + p_R)
      91            0 :             );
      92              : 
      93            0 :         if (S_L >= 0.0)
      94              :         {
      95            0 :             return F_L;
      96              :         }
      97            0 :         else if (S_R <= 0.0)
      98              :         {
      99            0 :             return F_R;
     100              :         }
     101              :         else
     102              :         {
     103            0 :             State dState = hi - lo;
     104            0 :             Set::Scalar dS = S_R - S_L;
     105            0 :             return (S_R*F_L - S_L*F_R + S_L*S_R*Flux(dState)) / dS;
     106              :         }
     107              :     }
     108              : };
     109              : }
     110              : }
     111              : }
     112              : 
     113              : #endif