Implement general EOS equations in Riemann Solver (#338)

* add nonlinear wavespeed correction

* first step: for ideal gas EOS

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* rename to match other variables

* store mass scalars in HydroState

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* add massScalars

* template HLLC with problem_t

* find dedr from Microphysics EOS

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* add num mass scalars in template

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* use microphysics to find eos quantities

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix bug

* added thermodynamic derivatives

* added expression for C_tilde_P

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* cosmetic change

* club derivative functions into one to increase speed

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* dont need to declare dedr dpdr dedp

* remove comments

* calculate Gamma for fundamental derivative

* use second order wavespeed correction

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* cleanup

* updated microp

* implement full equation for G

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* simplified expression for G

* updated microp: correct expression for d2pdr2_s

* define G directly in microp

* updated microp: add primordial chem thermo derivatives

---------

Co-authored-by: Ben Wibking <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/EOS.hpp

@@ -9,6 +9,7 @@
 
 #include <cmath>
 #include <optional>
+#include <tuple>
 
 #include "AMReX_Array.H"
 #include "AMReX_GpuQualifiers.H"
@@ -43,16 +44,24 @@ template <typename problem_t> class EOS
 	static constexpr int nmscalars_ = Physics_Traits<problem_t>::numMassScalars;
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputeTgasFromEint(amrex::Real rho, amrex::Real Eint, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {}) -> amrex::Real;
+
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputeEintFromTgas(amrex::Real rho, amrex::Real Tgas, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {}) -> amrex::Real;
+
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {})
 	    -> amrex::Real;
+
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputeEintTempDerivative(amrex::Real rho, amrex::Real Tgas, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {})
 	    -> amrex::Real;
+
+	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
+	ComputeOtherDerivatives(amrex::Real rho, amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {});
+
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputePressure(amrex::Real rho, amrex::Real Eint, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {}) -> amrex::Real;
+
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {}) -> amrex::Real;
 
@@ -187,7 +196,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
 EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real Tgas,
 					  const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars) -> amrex::Real
 {
-	// compute derivative of internal energy w/r/t temperature
+	// compute derivative of internal energy w/r/t temperature, given density and temperature
 	amrex::Real dEint_dT = NAN;
 
 #ifdef PRIMORDIAL_CHEM
@@ -222,6 +231,61 @@ EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re
 	return dEint_dT;
 }
 
+template <typename problem_t>
+AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDerivatives(const amrex::Real rho, const amrex::Real P,
+										      const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
+{
+	// compute derivative of specific internal energy w/r/t density, given density and pressure
+	amrex::Real deint_dRho = NAN;
+	// compute derivative of specific internal energy w/r/t density, given density and pressure
+	amrex::Real deint_dP = NAN;
+	// compute derivative of density w/r/t pressure, given density and pressure
+	amrex::Real dRho_dP = NAN;
+	// compute derivative of pressure w/r/t density at constant entropy, given density and pressure (needed for the fundamental derivative G)
+	amrex::Real dP_dRho_s = NAN;
+	// fundamental derivative
+	amrex::Real G = NAN;
+
+#ifdef PRIMORDIAL_CHEM
+	eos_t chemstate;
+	chemstate.rho = rho;
+	chemstate.p = P;
+	// initialize array of number densities
+	for (int ii = 0; ii < NumSpec; ++ii) {
+		chemstate.xn[ii] = -1.0;
+	}
+
+	if (massScalars) {
+		const auto &massArray = *massScalars;
+		for (int nn = 0; nn < nmscalars_; ++nn) {
+			chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho)
+		}
+	}
+
+	eos(eos_input_rp, chemstate);
+	deint_dRho = chemstate.dedr;
+	deint_dP = 1.0 / chemstate.dpde;
+	dRho_dP = 1.0 / (chemstate.dpdr * C::k_B / boltzmann_constant_);
+	dP_dRho_s = chemstate.cs * chemstate.cs;
+	G = chemstate.G;
+
+#else
+	if constexpr (gamma_ != 1.0) {
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.p = P;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rp, estate);
+		deint_dRho = estate.dedr;
+		deint_dP = 1.0 / estate.dpde;
+		dRho_dP = 1.0 / (estate.dpdr * C::k_B / boltzmann_constant_);
+		dP_dRho_s = estate.cs * estate.cs;
+		G = estate.G;
+	}
+#endif
+	return std::make_tuple(deint_dRho, deint_dP, dRho_dP, dP_dRho_s, G);
+}
+
 template <typename problem_t>
 AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputePressure(amrex::Real rho, amrex::Real Eint,
 									      const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)

src/HLLC.hpp

@@ -7,6 +7,7 @@
 #include <AMReX_REAL.H>
 
 #include "ArrayView.hpp"
+#include "EOS.hpp"
 #include "HydroState.hpp"
 #include "valarray.hpp"
 
@@ -17,9 +18,9 @@ namespace quokka::Riemann
 //  Minoshima & Miyoshi, "A low-dissipation HLLD approximate Riemann solver
 //  	for a very wide range of Mach numbers," JCP (2021).]
 //
-template <int N_scalars, int fluxdim>
-AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars> const &sL, quokka::HydroState<N_scalars> const &sR, const double gamma,
-					      const double du, const double dw) -> quokka::valarray<double, fluxdim>
+template <typename problem_t, int N_scalars, int N_mscalars, int fluxdim>
+AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars, N_mscalars> const &sL, quokka::HydroState<N_scalars, N_mscalars> const &sR,
+					      const double gamma, const double du, const double dw) -> quokka::valarray<double, fluxdim>
 {
 	// compute Roe averages
 
@@ -34,17 +35,49 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars> cons
 	const double H_R = (sR.E + sR.P) / sR.rho; // sR specific enthalpy
 	const double H_tilde = (wl * H_L + wr * H_R) * norm;
 	double cs_tilde = NAN;
+
+	// compute nonlinear wavespeed correction [Rider, Computers & Fluids 28 (1999) 741-777]
+	// (Only applicable in compressions. Significantly reduces slow-moving shock oscillations.)
+
+	const double dU = sL.u - sR.u;
+	double S_L = NAN;
+	double S_R = NAN;
 	if (gamma != 1.0) {
-		// TODO(ben): implement Roe average for general EOS
-		cs_tilde = std::sqrt((gamma - 1.) * (H_tilde - 0.5 * vsq_tilde));
+
+		auto [dedr_L, dedp_L, drdp_L, dpdr_s_L, G_L] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sL.rho, sL.P, sL.massScalar);
+		auto [dedr_R, dedp_R, drdp_R, dpdr_s_R, G_R] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sR.rho, sR.P, sR.massScalar);
+
+		// equation A.5a of Kershaw+1998
+		// need specific internal energy here
+		const double C_tilde_rho = 0.5 * ((sL.Eint / sL.rho) + (sR.Eint / sR.rho) + sL.rho * dedr_L + sR.rho * dedr_R);
+
+		// equation A.5b of Kershaw+1998
+		const double C_tilde_P = 0.5 * ((sL.Eint / sL.rho) * drdp_L + (sR.Eint / sR.rho) * drdp_R + sL.rho * dedp_L + sR.rho * dedp_R);
+
+		// equation 4.12 of Kershaw+1998
+		cs_tilde = std::sqrt((1.0 / C_tilde_P) * (H_tilde - 0.5 * vsq_tilde - C_tilde_rho));
+
+		const double s_NL = 0.5 * G_L * std::max(dU, 0.); // second-order wavespeed correction
+		const double s_NR = 0.5 * G_R * std::max(dU, 0.);
+
+		// compute wave speeds following Batten et al. (1997)
+		S_L = std::min(sL.u - (sL.cs + s_NL), u_tilde - (cs_tilde + s_NL));
+		S_R = std::max(sR.u + (sR.cs + s_NR), u_tilde + (cs_tilde + s_NR));
+
 	} else {
 		cs_tilde = 0.5 * (sL.cs + sR.cs);
-	}
+		const double G_gamma_L = 1.0;
+		const double G_gamma_R = 1.0;
+
+		const double G_L = 0.5 * (G_gamma_L + 1.);
+		const double G_R = 0.5 * (G_gamma_R + 1.);
 
-	// compute wave speeds following Batten et al. (1997)
+		const double s_NL = 0.5 * G_L * std::max(dU, 0.);
+		const double s_NR = 0.5 * G_R * std::max(dU, 0.);
 
-	const double S_L = std::min(sL.u - sL.cs, u_tilde - cs_tilde);
-	const double S_R = std::max(sR.u + sR.cs, u_tilde + cs_tilde);
+		S_L = std::min(sL.u - (sL.cs + s_NL), u_tilde - (cs_tilde + s_NL));
+		S_R = std::max(sR.u + (sR.cs + s_NR), u_tilde + (cs_tilde + s_NR));
+	}
 
 	// carbuncle correction [Eq. 10 of Minoshima & Miyoshi (2021)]
 
@@ -113,4 +146,4 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars> cons
 }
 } // namespace quokka::Riemann
 
-#endif // HLLC_HPP_
+#endif // HLLC_HPP_

src/HydroContact/test_hydro_contact.cpp

@@ -146,7 +146,7 @@ void RadhydroSimulation<ContactProblem>::computeReferenceSolution(amrex::MultiFa
 				const auto E = val_exact.at(HydroSystem<ContactProblem>::energy_index)[i];
 				const auto vx = xmom / rho;
 				const auto Eint = E - 0.5 * rho * (vx * vx);
-				const auto P = (quokka::EOS_Traits<ContactProblem>::gamma - 1.) * Eint;
+				const auto P = quokka::EOS<ContactProblem>::ComputePressure(rho, Eint);
 				d_exact.push_back(rho);
 				vx_exact.push_back(vx);
 				P_exact.push_back(P);

src/HydroSMS/test_hydro_sms.cpp

@@ -194,13 +194,13 @@ void RadhydroSimulation<ShocktubeProblem>::computeReferenceSolution(amrex::Multi
 			amrex::Real vx = vx_arr[i];
 			amrex::Real P = P_arr[i];
 
-			const auto gamma = quokka::EOS_Traits<ShocktubeProblem>::gamma;
 			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::density_index) = rho;
 			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::x1Momentum_index) = rho * vx;
 			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::x2Momentum_index) = 0.;
 			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::x3Momentum_index) = 0.;
-			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx);
-			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::internalEnergy_index) = P / (gamma - 1.);
+			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::energy_index) =
+			    quokka::EOS<ShocktubeProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx);
+			stateExact(i, j, k, HydroSystem<ShocktubeProblem>::internalEnergy_index) = quokka::EOS<ShocktubeProblem>::ComputeEintFromPres(rho, P);
 		});
 	}
 
@@ -221,7 +221,7 @@ void RadhydroSimulation<ShocktubeProblem>::computeReferenceSolution(amrex::Multi
 
 			amrex::Real xvel = xmom / rho;
 			amrex::Real Eint = Egas - xmom * xmom / (2.0 * rho);
-			amrex::Real pressure = (quokka::EOS_Traits<ShocktubeProblem>::gamma - 1.) * Eint;
+			amrex::Real pressure = quokka::EOS<ShocktubeProblem>::ComputePressure(rho, Eint);
 
 			d.at(i) = rho;
 			vx.at(i) = xvel;

src/HydroState.hpp

@@ -5,18 +5,19 @@
 
 namespace quokka
 {
-template <int N> struct HydroState {
-	double rho;		      // density
-	double u;		      // normal velocity component
-	double v;		      // transverse velocity component
-	double w;		      // 2nd transverse velocity component
-	double P;		      // pressure
-	double cs;		      // adiabatic sound speed
-	double E;		      // total energy density
-	double Eint;		      // internal energy density
-	std::array<double, N> scalar; // passive scalars
+template <int Nall, int Nmass> struct HydroState {
+	double rho;				   // density
+	double u;				   // normal velocity component
+	double v;				   // transverse velocity component
+	double w;				   // 2nd transverse velocity component
+	double P;				   // pressure
+	double cs;				   // adiabatic sound speed
+	double E;				   // total energy density
+	double Eint;				   // internal energy density
+	std::array<double, Nall> scalar;	   // passive scalars
+	amrex::GpuArray<double, Nmass> massScalar; // mass scalars
 };
 
 } // namespace quokka
 
-#endif // HYDROSTATE_HPP_
+#endif // HYDROSTATE_HPP_

src/hydro_system.hpp

@@ -917,7 +917,7 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 			velW_index = x2Velocity_index;
 		}
 
-		quokka::HydroState<nscalars_> sL{};
+		quokka::HydroState<nscalars_, nmscalars_> sL{};
 		sL.rho = rho_L;
 		sL.u = x1LeftState(i, j, k, velN_index);
 		sL.v = x1LeftState(i, j, k, velV_index);
@@ -927,7 +927,7 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 		sL.E = E_L;
 		sL.Eint = Eint_L;
 
-		quokka::HydroState<nscalars_> sR{};
+		quokka::HydroState<nscalars_, nmscalars_> sR{};
 		sR.rho = rho_R;
 		sR.u = x1RightState(i, j, k, velN_index);
 		sR.v = x1RightState(i, j, k, velV_index);
@@ -937,12 +937,17 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 		sR.E = E_R;
 		sR.Eint = Eint_R;
 
-		// The remaining components are passive scalars, so just copy them from
+		// The remaining components are mass scalars and passive scalars, so just copy them from
 		// x1LeftState and x1RightState into the (left, right) state vectors U_L and
 		// U_R
 		for (int n = 0; n < nscalars_; ++n) {
 			sL.scalar[n] = x1LeftState(i, j, k, scalar0_index + n);
 			sR.scalar[n] = x1RightState(i, j, k, scalar0_index + n);
+			// also store mass scalars separately
+			if (n < nmscalars_) {
+				sL.massScalar[n] = x1LeftState(i, j, k, scalar0_index + n);
+				sR.massScalar[n] = x1RightState(i, j, k, scalar0_index + n);
+			}
 		}
 
 		// difference in normal velocity along normal axis
@@ -964,7 +969,7 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 #endif
 
 		// solve the Riemann problem in canonical form
-		quokka::valarray<double, nvar_> F_canonical = quokka::Riemann::HLLC<nscalars_, nvar_>(sL, sR, gamma_, du, dw);
+		quokka::valarray<double, nvar_> F_canonical = quokka::Riemann::HLLC<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR, gamma_, du, dw);
 		quokka::valarray<double, nvar_> F = F_canonical;
 
 		// add artificial viscosity