Compute pressure from microphysics eos

.gitmodules

@@ -10,4 +10,4 @@
 [submodule "extern/Microphysics"]
 	path = extern/Microphysics
 	url = https://github.com/psharda/Microphysics
-	branch = development
+	branch = compute-cs

src/EOS.hpp

@@ -46,8 +46,15 @@ template <typename problem_t> class EOS
 	[[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
+	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;
 
       private:
 	static constexpr amrex::Real gamma_ = EOS_Traits<problem_t>::gamma;
@@ -60,11 +67,11 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeTgasFromEin
 										  const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
     -> amrex::Real
 {
-	// return temperature for an ideal gas
+	// return temperature for an ideal gas given density and internal energy
 	amrex::Real Tgas = NAN;
 
 #ifdef PRIMORDIAL_CHEM
-	burn_t chemstate;
+	eos_t chemstate;
 	chemstate.rho = rho;
 	chemstate.e = Eint / rho;
 	// initialize array of number densities
@@ -100,11 +107,11 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromTga
 										  const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
     -> amrex::Real
 {
-	// return internal energy density for a gamma-law ideal gas
+	// return internal energy density given density and temperature
 	amrex::Real Eint = NAN;
 
 #ifdef PRIMORDIAL_CHEM
-	burn_t chemstate;
+	eos_t chemstate;
 	chemstate.rho = rho;
 	// Define and initialize Tgas here
 	amrex::Real const Tgas_value = Tgas;
@@ -136,6 +143,45 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromTga
 	return Eint;
 }
 
+template <typename problem_t>
+AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure,
+										  const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
+    -> amrex::Real
+{
+	// return internal energy density given density and pressure
+	amrex::Real Eint = NAN;
+
+#ifdef PRIMORDIAL_CHEM
+	eos_t chemstate;
+	chemstate.rho = rho;
+	chemstate.p = Pressure;
+	// 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);
+	Eint = chemstate.e * chemstate.rho;
+#else
+	if constexpr (gamma_ != 1.0) {
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.p = Pressure;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rp, estate);
+		Eint = estate.e * rho;
+	}
+#endif
+	return Eint;
+}
+
 template <typename problem_t>
 AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
 EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real Tgas,
@@ -145,7 +191,7 @@ EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re
 	amrex::Real dEint_dT = NAN;
 
 #ifdef PRIMORDIAL_CHEM
-	burn_t chemstate;
+	eos_t chemstate;
 	chemstate.rho = rho;
 	// we don't need Tgas to find chemstate.dedT, but we still need to initialize chemstate.T because we are using the 'rt' EOS mode
 	chemstate.T = NAN;
@@ -176,6 +222,88 @@ 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>::ComputePressure(amrex::Real rho, amrex::Real Eint,
+									      const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
+    -> amrex::Real
+{
+	// return pressure for an ideal gas
+	amrex::Real P = NAN;
+#ifdef PRIMORDIAL_CHEM
+	eos_t chemstate;
+	chemstate.rho = rho;
+	chemstate.e = Eint / rho;
+	// 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_re, chemstate);
+	P = chemstate.p;
+#else
+	if constexpr (gamma_ != 1.0) {
+		chem_eos_t estate;
+		estate.rho = rho;
+		// if rho is 0, pass 0 to state.e
+		if (rho == 0.0) {
+			estate.e = 0;
+		} else {
+			estate.e = Eint / rho;
+		}
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_re, estate);
+		P = estate.p;
+	}
+#endif
+	return P;
+}
+
+template <typename problem_t>
+AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure,
+										const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
+    -> amrex::Real
+{
+	// return sound speed for an ideal gas
+	amrex::Real cs = NAN;
+
+#ifdef PRIMORDIAL_CHEM
+	eos_t chemstate;
+	chemstate.rho = rho;
+	chemstate.p = Pressure;
+	// 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);
+	cs = chemstate.cs;
+#else
+	if constexpr (gamma_ != 1.0) {
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.p = Pressure;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rp, estate);
+		cs = estate.cs;
+	}
+#endif
+	return cs;
+}
+
 } // namespace quokka
 
 #endif // EOS_HPP_

src/HydroBlast2D/test_hydro2d_blast.cpp

@@ -27,7 +27,7 @@ struct BlastProblem {
 
 template <> struct quokka::EOS_Traits<BlastProblem> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroBlast3D/test_hydro3d_blast.cpp

@@ -34,7 +34,7 @@ bool test_passes = false; // if one of the energy checks fails, set to false
 
 template <> struct quokka::EOS_Traits<SedovProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroContact/test_hydro_contact.cpp

@@ -70,16 +70,15 @@ template <> void RadhydroSimulation<ContactProblem>::setInitialConditionsOnGrid(
 		AMREX_ASSERT(!std::isnan(rho));
 		AMREX_ASSERT(!std::isnan(P));
 
-		const auto gamma = quokka::EOS_Traits<ContactProblem>::gamma;
 		for (int n = 0; n < ncomp_cc; ++n) {
 			state_cc(i, j, k, n) = 0.;
 		}
 		state_cc(i, j, k, HydroSystem<ContactProblem>::density_index) = rho;
 		state_cc(i, j, k, HydroSystem<ContactProblem>::x1Momentum_index) = rho * vx;
 		state_cc(i, j, k, HydroSystem<ContactProblem>::x2Momentum_index) = 0.;
 		state_cc(i, j, k, HydroSystem<ContactProblem>::x3Momentum_index) = 0.;
-		state_cc(i, j, k, HydroSystem<ContactProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx);
-		state_cc(i, j, k, HydroSystem<ContactProblem>::internalEnergy_index) = P / (gamma - 1.);
+		state_cc(i, j, k, HydroSystem<ContactProblem>::energy_index) = quokka::EOS<ContactProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx);
+		state_cc(i, j, k, HydroSystem<ContactProblem>::internalEnergy_index) = quokka::EOS<ContactProblem>::ComputeEintFromPres(rho, P);
 	});
 }
 
@@ -112,13 +111,13 @@ void RadhydroSimulation<ContactProblem>::computeReferenceSolution(amrex::MultiFa
 				stateExact(i, j, k, n) = 0.;
 			}
 
-			const auto gamma = quokka::EOS_Traits<ContactProblem>::gamma;
 			stateExact(i, j, k, HydroSystem<ContactProblem>::density_index) = rho;
 			stateExact(i, j, k, HydroSystem<ContactProblem>::x1Momentum_index) = rho * vx;
 			stateExact(i, j, k, HydroSystem<ContactProblem>::x2Momentum_index) = 0.;
 			stateExact(i, j, k, HydroSystem<ContactProblem>::x3Momentum_index) = 0.;
-			stateExact(i, j, k, HydroSystem<ContactProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx);
-			stateExact(i, j, k, HydroSystem<ContactProblem>::internalEnergy_index) = P / (gamma - 1.);
+			stateExact(i, j, k, HydroSystem<ContactProblem>::energy_index) =
+			    quokka::EOS<ContactProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx);
+			stateExact(i, j, k, HydroSystem<ContactProblem>::internalEnergy_index) = quokka::EOS<ContactProblem>::ComputeEintFromPres(rho, P);
 		});
 	}
 

src/HydroHighMach/test_hydro_highmach.cpp

@@ -30,7 +30,7 @@ struct HighMachProblem {
 
 template <> struct quokka::EOS_Traits<HighMachProblem> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroKelvinHelmholz/test_hydro2d_kh.cpp

@@ -24,7 +24,7 @@ struct KelvinHelmholzProblem {
 
 template <> struct quokka::EOS_Traits<KelvinHelmholzProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroQuirk/test_quirk.cpp

@@ -43,7 +43,7 @@ struct QuirkProblem {
 
 template <> struct quokka::EOS_Traits<QuirkProblem> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 
@@ -125,8 +125,8 @@ template <> void RadhydroSimulation<QuirkProblem>::setInitialConditionsOnGrid(qu
 		state_cc(i, j, k, HydroSystem<QuirkProblem>::x1Momentum_index) = rho * vx;
 		state_cc(i, j, k, HydroSystem<QuirkProblem>::x2Momentum_index) = rho * vy;
 		state_cc(i, j, k, HydroSystem<QuirkProblem>::x3Momentum_index) = rho * vz;
-		state_cc(i, j, k, HydroSystem<QuirkProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * v_sq;
-		state_cc(i, j, k, HydroSystem<QuirkProblem>::internalEnergy_index) = P / (gamma - 1.);
+		state_cc(i, j, k, HydroSystem<QuirkProblem>::energy_index) = quokka::EOS<QuirkProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * v_sq;
+		state_cc(i, j, k, HydroSystem<QuirkProblem>::internalEnergy_index) = quokka::EOS<QuirkProblem>::ComputeEintFromPres(rho, P);
 	});
 }
 

src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp

@@ -23,7 +23,7 @@ struct RichtmeyerMeshkovProblem {
 
 template <> struct quokka::EOS_Traits<RichtmeyerMeshkovProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroShocktube/test_hydro_shocktube.cpp

@@ -28,7 +28,7 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp

@@ -33,7 +33,7 @@ template <> struct SimulationData<ShocktubeProblem> {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/HydroVacuum/test_hydro_vacuum.cpp

@@ -27,7 +27,7 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = NAN;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = C::k_B;
 };
 

src/PassiveScalar/test_scalars.cpp

@@ -70,16 +70,15 @@ template <> void RadhydroSimulation<ScalarProblem>::setInitialConditionsOnGrid(q
 			scalar = 0.0;
 		}
 
-		const auto gamma = quokka::EOS_Traits<ScalarProblem>::gamma;
 		for (int n = 0; n < state_cc.nComp(); ++n) {
 			state_cc(i, j, k, n) = 0.;
 		}
 		state_cc(i, j, k, HydroSystem<ScalarProblem>::density_index) = rho;
 		state_cc(i, j, k, HydroSystem<ScalarProblem>::x1Momentum_index) = rho * vx;
 		state_cc(i, j, k, HydroSystem<ScalarProblem>::x2Momentum_index) = 0.;
 		state_cc(i, j, k, HydroSystem<ScalarProblem>::x3Momentum_index) = 0.;
-		state_cc(i, j, k, HydroSystem<ScalarProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx);
-		state_cc(i, j, k, HydroSystem<ScalarProblem>::internalEnergy_index) = P / (gamma - 1.);
+		state_cc(i, j, k, HydroSystem<ScalarProblem>::energy_index) = quokka::EOS<ScalarProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx);
+		state_cc(i, j, k, HydroSystem<ScalarProblem>::internalEnergy_index) = quokka::EOS<ScalarProblem>::ComputeEintFromPres(rho, P);
 		state_cc(i, j, k, HydroSystem<ScalarProblem>::scalar0_index) = scalar;
 	});
 }
@@ -112,16 +111,16 @@ void RadhydroSimulation<ScalarProblem>::computeReferenceSolution(amrex::MultiFab
 				scalar = 0.0;
 			}
 
-			const auto gamma = quokka::EOS_Traits<ScalarProblem>::gamma;
 			for (int n = 0; n < ncomp; ++n) {
 				stateExact(i, j, k, n) = 0.;
 			}
 			stateExact(i, j, k, HydroSystem<ScalarProblem>::density_index) = rho;
 			stateExact(i, j, k, HydroSystem<ScalarProblem>::x1Momentum_index) = rho * vx;
 			stateExact(i, j, k, HydroSystem<ScalarProblem>::x2Momentum_index) = 0.;
 			stateExact(i, j, k, HydroSystem<ScalarProblem>::x3Momentum_index) = 0.;
-			stateExact(i, j, k, HydroSystem<ScalarProblem>::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx);
-			stateExact(i, j, k, HydroSystem<ScalarProblem>::internalEnergy_index) = P / (gamma - 1.);
+			stateExact(i, j, k, HydroSystem<ScalarProblem>::energy_index) =
+			    quokka::EOS<ScalarProblem>::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx);
+			stateExact(i, j, k, HydroSystem<ScalarProblem>::internalEnergy_index) = quokka::EOS<ScalarProblem>::ComputeEintFromPres(rho, P);
 			stateExact(i, j, k, HydroSystem<ScalarProblem>::scalar0_index) = scalar;
 		});
 	}
@@ -166,7 +165,7 @@ void RadhydroSimulation<ScalarProblem>::computeReferenceSolution(amrex::MultiFab
 				const auto s = val_exact.at(HydroSystem<ScalarProblem>::scalar0_index)[i];
 				const auto vx = xmom / rho;
 				const auto Eint = E - 0.5 * rho * (vx * vx);
-				const auto P = (quokka::EOS_Traits<ScalarProblem>::gamma - 1.) * Eint;
+				const auto P = quokka::EOS<ScalarProblem>::ComputePressure(rho, Eint);
 				d_exact.push_back(rho);
 				vx_exact.push_back(vx);
 				P_exact.push_back(P);
@@ -180,7 +179,7 @@ void RadhydroSimulation<ScalarProblem>::computeReferenceSolution(amrex::MultiFab
 				const auto fs = values.at(HydroSystem<ScalarProblem>::scalar0_index)[i];
 				const auto fvx = fxmom / frho;
 				const auto fEint = fE - 0.5 * frho * (fvx * fvx);
-				const auto fP = (quokka::EOS_Traits<ScalarProblem>::gamma - 1.) * fEint;
+				const auto fP = quokka::EOS<ScalarProblem>::ComputePressure(frho, fEint);
 				d_final.push_back(frho);
 				vx_final.push_back(fvx);
 				P_final.push_back(fP);

src/RadMarshak/test_radiation_marshak.cpp

@@ -31,7 +31,7 @@ constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson;
 constexpr double T_initial = 1.0e-2;
 
 template <> struct quokka::EOS_Traits<SuOlsonProblem> {
-	static constexpr double mean_molecular_mass = 1.0;
+	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = 1.0;
 	static constexpr double gamma = 5. / 3.;
 };

src/RadMarshakCGS/test_radiation_marshak_cgs.cpp

@@ -33,7 +33,7 @@ constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson;
 constexpr double T_initial = 1.0e4; // K
 
 template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
-	static constexpr double mean_molecular_mass = C::m_u;
+	static constexpr double mean_molecular_weight = C::m_u;
 	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };