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/Cooling/test_cooling.cpp

@@ -37,6 +37,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<CoolingTest> {

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"
@@ -35,6 +36,7 @@ template <typename problem_t> struct EOS_Traits {
 	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 hydrogen_mass_code_units = hydrogen_mass_cgs;
 };
 
 template <typename problem_t> class EOS
@@ -84,8 +86,16 @@ 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);
+		// const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
+		// Tgas = Eint / (rho * c_v);
+		chem_eos_t chemstate;
+		chemstate.rho = rho;
+		chemstate.e = Eint / rho;
+		chemstate.mu = mean_molecular_weight_ / hydrogen_mass_code_units;
+		// amrex::Print() << "Tgas before " << chemstate.rho << " " << chemstate.e << " "<< chemstate.mu << " " << std::endl;
+		eos(eos_input_re, chemstate);
+		// amrex::Print() << "Tgas after "<< chemstate << std::endl;
+		Tgas = chemstate.T;
 	}
 #endif
 	return Tgas;

src/FCQuantities/test_fc_quantities.cpp

@@ -29,6 +29,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<FCQuantities> {

src/HydroBlast2D/test_hydro2d_blast.cpp

@@ -29,6 +29,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<BlastProblem> {

src/HydroBlast3D/test_hydro3d_blast.cpp

@@ -36,6 +36,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<SedovProblem> {

src/HydroContact/test_hydro_contact.cpp

@@ -25,6 +25,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ContactProblem> {

src/HydroHighMach/test_hydro_highmach.cpp

@@ -32,6 +32,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<HighMachProblem> {

src/HydroKelvinHelmholz/test_hydro2d_kh.cpp

@@ -26,6 +26,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<KelvinHelmholzProblem> {

src/HydroLeblanc/test_hydro_leblanc.cpp

@@ -32,6 +32,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroQuirk/test_quirk.cpp

@@ -45,6 +45,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<QuirkProblem> {

src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp

@@ -25,6 +25,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<RichtmeyerMeshkovProblem> {

src/HydroSMS/test_hydro_sms.cpp

@@ -25,6 +25,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShocktube/test_hydro_shocktube.cpp

@@ -30,6 +30,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp

@@ -35,6 +35,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroShuOsher/test_hydro_shuosher.cpp

@@ -27,6 +27,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroVacuum/test_hydro_vacuum.cpp

@@ -29,6 +29,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {

src/HydroWave/test_hydro_wave.cpp

@@ -25,6 +25,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<WaveProblem> {

src/ODEIntegration/test_ode.cpp

@@ -19,6 +19,7 @@ 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 gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 struct ODEUserData {

src/PassiveScalar/test_scalars.cpp

@@ -27,6 +27,7 @@ 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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ScalarProblem> {

src/PrimordialChem/test_primordial_chem.cpp

@@ -40,13 +40,6 @@ using amrex::Real;
 struct PrimordialChemTest {
 }; // dummy type to allow compile-type polymorphism via template specialization
 
-// Currently, microphysics uses its own EOS, and this one below is used by hydro. Need to only have one EOS at some point.
-template <> struct quokka::EOS_Traits<PrimordialChemTest> {
-	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;
-};
-
 template <> struct Physics_Traits<PrimordialChemTest> {
 	// cell-centred
 	static constexpr bool is_hydro_enabled = true;

src/RadBeam/test_radiation_beam.cpp

@@ -36,6 +36,7 @@ template <> struct quokka::EOS_Traits<BeamProblem> {
 	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<BeamProblem> {

src/RadForce/test_radiation_force.cpp

@@ -48,6 +48,7 @@ template <> struct quokka::EOS_Traits<TubeProblem> {
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = gamma_gas;
 	static constexpr double cs_isothermal = a0; // only used when gamma = 1
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<TubeProblem> {

src/RadMarshak/test_radiation_marshak.cpp

@@ -34,6 +34,7 @@ template <> struct quokka::EOS_Traits<SuOlsonProblem> {
 	static constexpr double mean_molecular_mass = 1.0;
 	static constexpr double boltzmann_constant = 1.0;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = 1.0;
 };
 
 template <> struct RadSystem_Traits<SuOlsonProblem> {

src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp

@@ -29,6 +29,7 @@ template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
 	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<SuOlsonProblemCgs> {

src/RadMarshakCGS/test_radiation_marshak_cgs.cpp

@@ -36,6 +36,7 @@ template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
 	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<SuOlsonProblemCgs> {

src/RadMatterCoupling/test_radiation_matter_coupling.cpp

@@ -34,6 +34,7 @@ template <> struct quokka::EOS_Traits<CouplingProblem> {
 	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<CouplingProblem> {

src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp

@@ -36,6 +36,7 @@ template <> struct quokka::EOS_Traits<CouplingProblem> {
 	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<CouplingProblem> {

src/RadPulse/test_radiation_pulse.cpp

@@ -30,6 +30,7 @@ template <> struct quokka::EOS_Traits<PulseProblem> {
 	static constexpr double mean_molecular_weight = 1.0;
 	static constexpr double boltzmann_constant = (2. / 3.);
 	static constexpr double gamma = 5. / 3.;
+	static constexpr hydrogen_mass_code_units = 1.0;
 };
 
 template <> struct RadSystem_Traits<PulseProblem> {

src/RadShadow/test_radiation_shadow.cpp

@@ -37,6 +37,7 @@ template <> struct quokka::EOS_Traits<ShadowProblem> {
 	static constexpr double mean_molecular_weight = 10. * quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<ShadowProblem> {

src/RadTophat/test_radiation_tophat.cpp

@@ -41,6 +41,7 @@ template <> struct quokka::EOS_Traits<TophatProblem> {
 	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = 5. / 3.;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<TophatProblem> {

src/RadTube/test_radiation_tube.cpp

@@ -40,6 +40,7 @@ template <> struct quokka::EOS_Traits<TubeProblem> {
 	static constexpr double mean_molecular_weight = mu;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
 	static constexpr double gamma = gamma_gas;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<TubeProblem> {

src/RadhydroShell/test_radhydro_shell.cpp

@@ -48,6 +48,7 @@ template <> struct quokka::EOS_Traits<ShellProblem> {
 	static constexpr double mean_molecular_weight = 2.2 * m_H;
 	static constexpr double boltzmann_constant = k_B;
 	static constexpr double gamma = gamma_gas;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct RadSystem_Traits<ShellProblem> {

src/RadhydroShock/test_radhydro_shock.cpp

@@ -65,6 +65,7 @@ template <> struct quokka::EOS_Traits<ShockProblem> {
 	static constexpr double mean_molecular_weight = mu;
 	static constexpr double boltzmann_constant = k_B;
 	static constexpr double gamma = gamma_gas;
+	static constexpr double hydrogen_mass_code_units = 1.0;
 };
 
 template <> struct Physics_Traits<ShockProblem> {

src/RadhydroShockCGS/test_radhydro_shock_cgs.cpp

@@ -68,6 +68,7 @@ template <> struct quokka::EOS_Traits<ShockProblem> {
 	static constexpr double mean_molecular_weight = m_H;
 	static constexpr double boltzmann_constant = k_B;
 	static constexpr double gamma = gamma_gas;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShockProblem> {

src/RadhydroSimulation.hpp

@@ -52,6 +52,8 @@
 #include "radiation_system.hpp"
 #include "simulation.hpp"
 
+#include "eos.H"
+
 // Simulation class should be initialized only once per program (i.e., is a singleton)
 template <typename problem_t> class RadhydroSimulation : public AMRSimulation<problem_t>
 {
@@ -127,13 +129,25 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 		defineComponentNames();
 		// read in runtime parameters
 		readParmParse();
+		// initialize Microphysics params
+		init_extern_parameters();
+		// initialize Microphysics EOS
+		amrex::Real small_temp = 1e-10;
+		amrex::Real small_dens = 1e-100;
+		eos_init(small_temp, small_dens);
 	}
 
 	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc) : AMRSimulation<problem_t>(BCs_cc)
 	{
 		defineComponentNames();
 		// read in runtime parameters
 		readParmParse();
+		// initialize Microphysics params
+		init_extern_parameters();
+		// initialize Microphysics EOS
+		amrex::Real small_temp = 1e-10;
+		amrex::Real small_dens = 1e-100;
+		eos_init(small_temp, small_dens);
 	}
 
 	[[nodiscard]] static auto getScalarVariableNames() -> std::vector<std::string>;

src/RayleighTaylor2D/test_hydro2d_rt.cpp

@@ -26,6 +26,7 @@ template <> struct quokka::EOS_Traits<RTProblem> {
 	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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<RTProblem> {

src/RayleighTaylor3D/test_hydro3d_rt.cpp

@@ -26,6 +26,7 @@ template <> struct quokka::EOS_Traits<RTProblem> {
 	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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<RTProblem> {

src/ShockCloud/cloud.cpp

@@ -42,6 +42,7 @@ template <> struct quokka::EOS_Traits<ShockCloud> {
 	static constexpr double gamma = 5. / 3.; // default value
 	static constexpr double mean_molecular_weight = NAN;
 	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct Physics_Traits<ShockCloud> {

src/SphericalCollapse/spherical_collapse.cpp

@@ -30,6 +30,7 @@ template <> struct quokka::EOS_Traits<CollapseProblem> {
 	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 hydrogen_mass_code_units = quokka::hydrogen_mass_cgs;
 };
 
 template <> struct HydroSystem_Traits<CollapseProblem> {