Use microphysics EOS in Quokka's EOS.hpp for all quokka tests + use amu instead of hydrogen mass throughout Quokka

.gitmodules

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

src/CMakeLists.txt

@@ -95,7 +95,9 @@ link_libraries(hdf5::hdf5)
 include(CTest)
 
 #create an object library for files that should be compiled with all test problems
-set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp")
+set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" 
+                      "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp"
+                      "${gamma_law_sources}")
 #we don't use it anymore because it gives issues on Cray systems
 #add_library(QuokkaObj OBJECT ${QuokkaObjSources} ${gamma_law_sources})
 #if(AMReX_GPU_BACKEND MATCHES "CUDA")

src/CloudyCooling.hpp

@@ -81,7 +81,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 {
 	// interpolate cooling rates from Cloudy tables
 	const Real rhoH = rho * cloudy_H_mass_fraction; // mass density of H species
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 	const Real log_nH = std::log10(nH);
 	const Real log_T = std::log10(T);
 
@@ -102,7 +102,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 
 	// compute electron density
 	// N.B. it is absolutely critical to include the metal contribution here!
-	double n_e = (rho / quokka::hydrogen_mass_cgs) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / quokka::hydrogen_mass_cgs));
+	double n_e = (rho / (C::m_p + C::m_e)) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / (C::m_p + C::m_e)));
 	// the approximation for the metals contribution to e- fails at high densities (~1e3 or higher)
 	n_e = std::max(n_e, 1.0e-4 * nH);
 
@@ -118,7 +118,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 	// Compton term (CMB photons)
 	// [e.g., Hirata 2018: doi:10.1093/mnras/stx2854]
 	constexpr double Gamma_C = (8. * sigma_T * E_cmb) / (3. * electron_mass_cgs * c_light_cgs_);
-	constexpr double C_n = Gamma_C * quokka::boltzmann_constant_cgs / (5. / 3. - 1.0);
+	constexpr double C_n = Gamma_C * C::k_B / (5. / 3. - 1.0);
 	const double compton_CMB = -C_n * (T - T_cmb) * n_e;
 	Edot += compton_CMB;
 
@@ -129,13 +129,13 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto ComputeEgasFromTgas(double rho, do
 {
 	// convert Egas (internal gas energy) to temperature
 	const Real rhoH = rho * cloudy_H_mass_fraction;
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 
 	// compute mu from mu(T) table
 	const Real mu = interpolate2d(std::log10(nH), std::log10(Tgas), tables.log_nH, tables.log_Tgas, tables.meanMolWeight);
 
 	// compute thermal gas energy
-	const Real Egas = (rho / (quokka::hydrogen_mass_cgs * mu)) * quokka::boltzmann_constant_cgs * Tgas / (gamma - 1.);
+	const Real Egas = (rho / ((C::m_p + C::m_e) * mu)) * C::k_B * Tgas / (gamma - 1.);
 	return Egas;
 }
 
@@ -157,12 +157,12 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto ComputeTgasFromEgas(double rho, do
 
 	// solve for temperature given Eint (with fixed adiabatic index gamma)
 	const Real rhoH = rho * cloudy_H_mass_fraction;
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 	const Real log_nH = std::log10(nH);
 
 	// mean molecular weight (in Grackle tables) is defined w/r/t
 	// hydrogen_mass_cgs_
-	const Real C = (gamma - 1.) * Egas / (quokka::boltzmann_constant_cgs * (rho / quokka::hydrogen_mass_cgs));
+	const Real C = (gamma - 1.) * Egas / (C::k_B * (rho / (C::m_p + C::m_e)));
 
 	// solve for mu(T)*C == T.
 	// (Grackle does this with a fixed-point iteration. We use a more robust

src/Cooling/test_cooling.cpp

@@ -30,13 +30,13 @@ using amrex::Real;
 struct CoolingTest {
 }; // dummy type to allow compile-type polymorphism via template specialization
 
-constexpr double m_H = quokka::hydrogen_mass_cgs;
+constexpr double m_H = C::m_u;
 constexpr double seconds_in_year = 3.154e7;
 
 template <> struct quokka::EOS_Traits<CoolingTest> {
 	static constexpr double gamma = 5. / 3.; // default value
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<CoolingTest> {
@@ -131,7 +131,7 @@ template <> void RadhydroSimulation<CoolingTest>::setInitialConditionsOnGrid(quo
 		Real xmom = 0;
 		Real ymom = 0;
 		Real zmom = 0;
-		Real const P = 4.0e4 * quokka::boltzmann_constant_cgs; // erg cm^-3
+		Real const P = 4.0e4 * C::k_B; // erg cm^-3
 		Real Eint = (quokka::EOS_Traits<CoolingTest>::gamma - 1.) * P;
 
 		Real const Egas = RadSystem<CoolingTest>::ComputeEgasFromEint(rho, xmom, ymom, zmom, Eint);

src/EOS.hpp

@@ -14,6 +14,7 @@
 #include "AMReX_GpuQualifiers.H"
 #include "AMReX_REAL.H"
 #include "physics_info.hpp"
+#include <AMReX_Print.H>
 
 #include "burn_type.H"
 #include "eos.H"
@@ -25,16 +26,14 @@
 
 namespace quokka
 {
-static constexpr double boltzmann_constant_cgs = 1.380658e-16; // cgs
-static constexpr double hydrogen_mass_cgs = 1.6726231e-24;     // cgs
 
 // specify default values for ideal gamma-law EOS
 //
 template <typename problem_t> struct EOS_Traits {
 	static constexpr double gamma = 5. / 3.;     // default value
 	static constexpr double cs_isothermal = NAN; // only used when gamma = 1
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <typename problem_t> class EOS
@@ -84,8 +83,13 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeTgasFromEin
 	Tgas = chemstate.T;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		Tgas = Eint / (rho * c_v);
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.e = Eint / rho;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_re, estate);
+		// scale returned temperature in case boltzmann constant is dimensionless
+		Tgas = estate.T * C::k_B / boltzmann_constant_;
 	}
 #endif
 	return Tgas;
@@ -121,16 +125,20 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromTga
 	Eint = chemstate.e * chemstate.rho;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		Eint = rho * c_v * Tgas;
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.T = Tgas;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rt, estate);
+		Eint = estate.e * rho * boltzmann_constant_ / C::k_B;
 	}
 #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*/,
+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
@@ -157,8 +165,12 @@ EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re
 	dEint_dT = chemstate.dedT * chemstate.rho;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		dEint_dT = rho * c_v;
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.T = Tgas;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rt, estate);
+		dEint_dT = estate.dedT * rho * boltzmann_constant_ / C::k_B;
 	}
 #endif
 	return dEint_dT;

src/FCQuantities/test_fc_quantities.cpp

@@ -27,8 +27,8 @@ struct FCQuantities {
 
 template <> struct quokka::EOS_Traits<FCQuantities> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<FCQuantities> {

src/HydroBlast2D/test_hydro2d_blast.cpp

@@ -28,7 +28,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<BlastProblem> {

src/HydroBlast3D/test_hydro3d_blast.cpp

@@ -35,7 +35,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<SedovProblem> {

src/HydroContact/test_hydro_contact.cpp

@@ -23,8 +23,8 @@ struct ContactProblem {
 
 template <> struct quokka::EOS_Traits<ContactProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ContactProblem> {

src/HydroHighMach/test_hydro_highmach.cpp

@@ -31,7 +31,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<HighMachProblem> {

src/HydroKelvinHelmholz/test_hydro2d_kh.cpp

@@ -25,7 +25,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<KelvinHelmholzProblem> {

src/HydroLeblanc/test_hydro_leblanc.cpp

@@ -30,8 +30,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = (5. / 3.);
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroQuirk/test_quirk.cpp

@@ -44,7 +44,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<QuirkProblem> {

src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp

@@ -24,7 +24,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<RichtmeyerMeshkovProblem> {

src/HydroSMS/test_hydro_sms.cpp

@@ -23,8 +23,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShocktube/test_hydro_shocktube.cpp

@@ -29,7 +29,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp

@@ -34,7 +34,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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShuOsher/test_hydro_shuosher.cpp

@@ -25,8 +25,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroVacuum/test_hydro_vacuum.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 boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroWave/test_hydro_wave.cpp

@@ -23,8 +23,8 @@ struct WaveProblem {
 
 template <> struct quokka::EOS_Traits<WaveProblem> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<WaveProblem> {

src/ODEIntegration/test_ode.cpp

@@ -8,16 +8,18 @@
 ///
 
 #include "test_ode.hpp"
+#include "eos.H"
+#include "extern_parameters.H"
 #include "radiation_system.hpp"
 
 using amrex::Real;
 
-constexpr double Tgas0 = 6000.;				  // K
-constexpr double rho0 = 0.01 * quokka::hydrogen_mass_cgs; // g cm^-3
+constexpr double Tgas0 = 6000.;			  // K
+constexpr double rho0 = 0.01 * (C::m_p + C::m_e); // g cm^-3
 
 template <> struct quokka::EOS_Traits<ODETest> {
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
@@ -30,7 +32,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cooling_function(Real const rho, R
 	// use fitting function from Koyama & Inutsuka (2002)
 	Real gamma_heat = 2.0e-26;
 	Real lambda_cool = gamma_heat * (1.0e7 * std::exp(-114800. / (T + 1000.)) + 14. * std::sqrt(T) * std::exp(-92. / T));
-	Real rho_over_mh = rho / quokka::hydrogen_mass_cgs;
+	Real rho_over_mh = rho / (C::m_p + C::m_e);
 	Real cooling_source_term = rho_over_mh * gamma_heat - (rho_over_mh * rho_over_mh) * lambda_cool;
 	return cooling_source_term;
 }
@@ -52,6 +54,12 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto user_rhs(Real /*t*/, quokka::valar
 
 auto problem_main() -> int
 {
+	// initialize EOS
+	init_extern_parameters();
+	Real small_temp = 1e-10;
+	Real small_dens = 1e-100;
+	eos_init(small_temp, small_dens);
+
 	// set up initial conditions
 	const Real Eint0 = quokka::EOS<ODETest>::ComputeEintFromTgas(rho0, Tgas0);
 	const Real Edot0 = cooling_function(rho0, Tgas0);

src/PassiveScalar/test_scalars.cpp

@@ -25,8 +25,8 @@ struct ScalarProblem {
 
 template <> struct quokka::EOS_Traits<ScalarProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ScalarProblem> {