[WIP] Hook the exact Riemann solver into Castro

Util/exact_riemann/exact_riemann_sample.H → Source/hydro/exact_riemann_sample.H

@@ -14,7 +14,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real
  const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r,
  const amrex::Real ustar, const amrex::Real pstar,
  const amrex::Real W_l, const amrex::Real W_r,
- const amrex::Real x, const amrex::Real xjump, const amrex::Real time,
+ const amrex::Real xi,
  amrex::Real& rho, amrex::Real& u, amrex::Real& p, amrex::Real* xn) {
 
  // get the initial sound speeds
@@ -47,11 +47,6 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real
 
  // This follows from the discussion around C&G Eq. 15
 
- // compute xi = x/t -- this is the similarity variable for the
- // solution
-
- amrex::Real xi = (x - xjump) / time;
-
  // check which side of the contact we need to worry about
 
  amrex::Real chi = std::copysign(1.0_rt, xi - ustar);

Source/hydro/riemann.cpp

@@ -71,7 +71,9 @@ Castro::cmpflx_plus_godunov(const Box& bx,
  {
 
 
- if (riemann_solver == 0 || riemann_solver == 1) {
+ if (riemann_solver == riemann_constants::TWO_SHOCK_CGF ||
+ riemann_solver == riemann_constants::TWO_SHOCK_CG ||
+ riemann_solver == riemann_constants::EXACT) {
  // approximate state Riemann solvers
 
  // first find the interface state on the current interface
@@ -120,7 +122,7 @@ Castro::cmpflx_plus_godunov(const Box& bx,
  }
  }
 
- } else if (riemann_solver == 2) {
+ } else if (riemann_solver == riemann_constants::HLLC) {
  // HLLC
  HLL::HLLC(i, j, k, idir,
  qm, qp,

Source/hydro/riemann_constants.H

@@ -6,6 +6,14 @@
 using namespace amrex::literals;
 
 namespace riemann_constants {
+
+ enum RiemannSolver : std::uint8_t {
+ TWO_SHOCK_CGF=0,
+ TWO_SHOCK_CG,
+ HLLC,
+ EXACT
+ };
+
  constexpr amrex::Real smlp1 = 1.e-10_rt;
  constexpr amrex::Real small = 1.e-8_rt;
  constexpr amrex::Real smallu = 1.e-12_rt;

Source/hydro/riemann_solvers.H

@@ -4,6 +4,8 @@
 #include <Castro_util.H>
 #include <riemann_type.H>
 #include <riemann_2shock_solvers.H>
+#include <exact_riemann_star_state.H>
+#include <exact_riemann_sample.H>
 #ifdef HYBRID_MOMENTUM
 #include <hybrid.H>
 #endif
@@ -318,18 +320,87 @@ riemann_state(const int i, const int j, const int k, const int idir,
 
 
  // Solve Riemann problem
- if (riemann_solver == 0) {
+ if (riemann_solver == riemann_constants::TWO_SHOCK_CGF) {
  // Colella, Glaz, & Ferguson solver
 
  TwoShock::riemannus(ql, qr, raux, qint);
 
- } else if (riemann_solver == 1) {
+ } else if (riemann_solver == riemann_constants::TWO_SHOCK_CG) {
  // Colella & Glaz solver
 
 #ifndef RADIATION
  TwoShock::riemanncg(ql, qr, raux, qint);
 #endif
 
+ } else if (riemann_solver == riemann_constants::EXACT) {
+ // Exact Riemann solver
+
+ // solve for the star state
+
+ amrex::Real ustar{};
+ amrex::Real pstar{};
+ amrex::Real W_l{};
+ amrex::Real W_r{};
+
+ amrex::Real xn_l[NumSpec];
+ for (int n = 0; n < NumSpec; ++n) {
+ xn_l[n] = qm(i,j,k,QFS+n);
+ }
+ amrex::Real xn_r[NumSpec];
+ for (int n = 0; n < NumSpec; ++n) {
+ xn_r[n] = qp(i,j,k,QFS+n);
+ }
+
+ riemann_star_state(ql.rho, ql.un, ql.p, xn_l,
+ qr.rho, qr.un, qr.p, xn_r,
+ ustar, pstar, W_l, W_r);
+
+ // sample the solution at xi = 0
+
+ amrex::Real xi{};
+
+ amrex::Real rho{};
+ amrex::Real u{};
+ amrex::Real p{};
+ amrex::Real xn_s[NumSpec]{};
+
+ riemann_sample(ql.rho, ql.un, ql.p, xn_l,
+ qr.rho, qr.un, qr.p, xn_r,
+ ustar, pstar, W_l, W_r,
+ xi,
+ rho, u, p, xn_s);
+
+ // pack the RiemannState interface state
+
+ qint.rho = rho;
+ qint.un = u;
+ qint.p = p;
+
+ // get the energy
+ eos_rep_t eos_state;
+ eos_state.rho = qint.rho;
+ eos_state.T = castro::T_guess;
+ eos_state.p = qint.p;
+ for (int n = 0; n < NumSpec; ++n) {
+ eos_state.xn[n] = xn_s[n];
+ }
+
+ eos(eos_input_rp, eos_state);
+
+ qint.rhoe = qint.rho * eos_state.e;
+ qint.gamc = eos_state.gam1;
+
+ // we need to manually do the transverse sampling based on the
+ // contact
+
+ if (ustar > 0) {
+ qint.ut = ql.ut;
+ qint.utt = ql.utt;
+ } else {
+ qint.ut = qr.ut;
+ qint.utt = qr.utt;
+ }
+
 #ifndef AMREX_USE_GPU
  } else {
  amrex::Error("ERROR: invalid value of riemann_solver");

Util/exact_riemann/README.md

@@ -1,9 +1,12 @@
-This is an exact Riemann solver for a general equation of state. It
-follows the outline for Colella & Glaz 1985, section 1. This is too
-slow to use in an actual hydro run, but instead is intended to
-generate exact solutions to the Riemann problem for comparison with
-Castro shocktube output. Several inputs files for Helmholtz EOS-based
-shocktubes are provided.
+# Exact Riemann solver
+
+This is a driver for the exact Riemann solver for a general equation
+of state. The main implementation is in Source/hydro.
+
+The exact Riemann solver follows the outline for Colella & Glaz 1985,
+section 1 and this driver is intended to generate exact solutions to
+the Riemann problem for comparison with Castro shocktube output.
+Several inputs files for Helmholtz EOS-based shocktubes are provided.
 
 This solver is used in Zingale & Katz (2015):
 

Util/exact_riemann/exact_riemann.H

@@ -78,10 +78,15 @@ exact_riemann() {
 
  amrex::Real rho, u, p, xn_s[NumSpec];
 
+ // compute xi = x/t -- this is the similarity variable for the
+ // solution
+
+ amrex::Real xi = (x - problem::xjump) / problem::t;
+
  riemann_sample(problem::rho_l, problem::u_l, problem::p_l, xn,
  problem::rho_r, problem::u_r, problem::p_r, xn,
  ustar, pstar, W_l, W_r,
- x, problem::xjump, problem::t,
+ xi,
  rho, u, p, xn_s);
 
  // get the thermodynamics for this state for output