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-pressure

src/EOS.hpp

@@ -48,6 +48,8 @@ template <typename problem_t> class EOS
 	[[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;
 
       private:
 	static constexpr amrex::Real gamma_ = EOS_Traits<problem_t>::gamma;
@@ -176,6 +178,44 @@ 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; // cannot use burn_t because it does not have pressure
+	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 * rho;
+#else
+	if constexpr (gamma_ != 1.0) {
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.e = Eint / rho;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_re, estate);
+		P = estate.p * rho;
+	}
+#endif
+	return P;
+}
+
 } // namespace quokka
 
 #endif // EOS_HPP_

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/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/hydro_system.hpp

@@ -150,7 +150,7 @@ template <typename problem_t> void HydroSystem<problem_t>::ConservedToPrimitive(
 		const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz);
 		const auto Eint_cons = E - kinetic_energy;
 
-		const amrex::Real Pgas = Eint_cons * (HydroSystem<problem_t>::gamma_ - 1.0);
+		const amrex::Real Pgas = ComputePressure(cons[bx], i, j, k);
 		const amrex::Real eint_cons = Eint_cons / rho;
 		const amrex::Real eint_aux = Eint_aux / rho;
 
@@ -202,9 +202,7 @@ template <typename problem_t> auto HydroSystem<problem_t>::maxSignalSpeedLocal(a
 					if constexpr (is_eos_isothermal()) {
 						cs = cs_iso_;
 					} else {
-						const auto Etot = cons[bx](i, j, k, HydroSystem<problem_t>::energy_index);
-						const auto Eint = Etot - kinetic_energy;
-						const auto P = Eint * (HydroSystem<problem_t>::gamma_ - 1.0);
+						const auto P = ComputePressure(cons[bx], i, j, k);
 						cs = std::sqrt(HydroSystem<problem_t>::gamma_ * P / rho);
 					}
 					return {cs + abs_vel};
@@ -233,11 +231,8 @@ void HydroSystem<problem_t>::ComputeMaxSignalSpeed(amrex::Array4<const amrex::Re
 		if constexpr (is_eos_isothermal()) {
 			cs = cs_iso_;
 		} else {
-			const auto E = cons(i, j, k, energy_index); // *total* gas energy per unit volume
-			AMREX_ASSERT(!std::isnan(E));
-			const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz);
-			const auto thermal_energy = E - kinetic_energy;
-			const auto P = thermal_energy * (HydroSystem<problem_t>::gamma_ - 1.0);
+			const auto P = ComputePressure(cons, i, j, k);
+			AMREX_ASSERT(!std::isnan(P));
 			cs = std::sqrt(HydroSystem<problem_t>::gamma_ * P / rho);
 		}
 		AMREX_ASSERT(cs > 0.);
@@ -265,7 +260,7 @@ template <typename problem_t> auto HydroSystem<problem_t>::CheckStatesValid(amre
 					const auto vz = pz / rho;
 					const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz);
 					const auto thermal_energy = E - kinetic_energy;
-					const auto P = thermal_energy * (HydroSystem<problem_t>::gamma_ - 1.0);
+					const auto P = ComputePressure(cons[bx], i, j, k);
 
 					bool negativeDensity = (rho <= 0.);
 					bool negativePressure = (P <= 0.);
@@ -300,7 +295,10 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem<problem_t>::ComputePressure
 	const auto vz = pz / rho;
 	const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz);
 	const auto thermal_energy = E - kinetic_energy;
-	const auto P = thermal_energy * (HydroSystem<problem_t>::gamma_ - 1.0);
+
+	amrex::GpuArray<Real, nmscalars_> massScalars = RadSystem<problem_t>::ComputeMassScalars(cons, i, j, k);
+	amrex::Real P = quokka::EOS<problem_t>::ComputePressure(rho, thermal_energy, massScalars);
+
 	return P;
 }