Add support for unit systems: CGS, CONSTANTS, CUSTOM

src/QuokkaSimulation.hpp

@@ -321,6 +321,53 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::defineComponentN
 	}
 }
 
+// initialize metadata
+template <typename problem_t> void AMRSimulation<problem_t>::initializeSimulationMetadata()
+{
+	if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CONSTANTS) {
+		// if unit system is CONSTANTS, the units are not well defined unless all four constants, G, k_B, c, and a_rad, are defined. However, in a hydro
+		// simulation, only k_B is defined. In a radiation-hydrodynamics simulation, only k_B, c, and a_rad are defined. Besides, CONSTANTS is only used
+		// for testing purposes, so we don't care about the units in that case.
+		simulationMetadata_["unit_length"] = NAN;
+		simulationMetadata_["unit_mass"] = NAN;
+		simulationMetadata_["unit_time"] = NAN;
+		simulationMetadata_["unit_temperature"] = NAN;
+
+		// constants
+		simulationMetadata_["k_B"] = Physics_Traits<problem_t>::boltzmann_constant;
+		simulationMetadata_["G"] = Physics_Traits<problem_t>::gravitational_constant;
+		if constexpr (Physics_Traits<problem_t>::is_radiation_enabled) {
+			simulationMetadata_["c"] = Physics_Traits<problem_t>::c_light;
+			simulationMetadata_["c_hat"] = Physics_Traits<problem_t>::c_light * RadSystem_Traits<problem_t>::c_hat_over_c;
+			simulationMetadata_["a_rad"] = Physics_Traits<problem_t>::radiation_constant;
+		}
+	} else {
+		// units
+		simulationMetadata_["unit_length"] = unit_length;
+		simulationMetadata_["unit_mass"] = unit_mass;
+		simulationMetadata_["unit_time"] = unit_time;
+		simulationMetadata_["unit_temperature"] = unit_temperature;
+
+		// constants
+		double k_B = NAN;
+		if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CGS) {
+			k_B = C::k_B;
+		} else if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CUSTOM) {
+			// Have to do a conversion because EOS class is not accessible here
+			k_B = C::k_B /
+			      (Physics_Traits<problem_t>::unit_length * Physics_Traits<problem_t>::unit_length * Physics_Traits<problem_t>::unit_mass /
+			       (Physics_Traits<problem_t>::unit_time * Physics_Traits<problem_t>::unit_time) / Physics_Traits<problem_t>::unit_temperature);
+		}
+		simulationMetadata_["k_B"] = k_B;
+		simulationMetadata_["G"] = Gconst_;
+		if constexpr (Physics_Traits<problem_t>::is_radiation_enabled) {
+			simulationMetadata_["c"] = RadSystem<problem_t>::c_light_;
+			simulationMetadata_["c_hat"] = RadSystem<problem_t>::c_hat_;
+			simulationMetadata_["a_rad"] = RadSystem<problem_t>::radiation_constant_;
+		}
+	}
+}
+
 template <typename problem_t> auto QuokkaSimulation<problem_t>::getScalarVariableNames() -> std::vector<std::string>
 {
 	// return vector of names for the passive scalars

src/hydro/EOS.hpp

@@ -38,6 +38,9 @@ template <typename problem_t> struct EOS_Traits {
 
 template <typename problem_t> class EOS
 {
+      private:
+	static constexpr amrex::Real gamma_ = EOS_Traits<problem_t>::gamma;
+	static constexpr amrex::Real mean_molecular_weight_ = EOS_Traits<problem_t>::mean_molecular_weight;
 
       public:
 	static constexpr int nmscalars_ = Physics_Traits<problem_t>::numMassScalars;
@@ -65,10 +68,18 @@ template <typename problem_t> class EOS
 	ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure, std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> const &massScalars = {})
 	    -> amrex::Real;
 
-      private:
-	static constexpr amrex::Real gamma_ = EOS_Traits<problem_t>::gamma;
-	static constexpr amrex::Real boltzmann_constant_ = EOS_Traits<problem_t>::boltzmann_constant;
-	static constexpr amrex::Real mean_molecular_weight_ = EOS_Traits<problem_t>::mean_molecular_weight;
+	static constexpr amrex::Real boltzmann_constant_ = []() constexpr {
+		if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CGS) {
+			return C::k_B;
+		} else if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CONSTANTS) {
+			return Physics_Traits<problem_t>::boltzmann_constant;
+		} else if constexpr (Physics_Traits<problem_t>::unit_system == UnitSystem::CUSTOM) {
+			// k_B / k_B_bar = u_l^2 * u_m / u_t^2 / u_T
+			return C::k_B /
+			       (Physics_Traits<problem_t>::unit_length * Physics_Traits<problem_t>::unit_length * Physics_Traits<problem_t>::unit_mass /
+				(Physics_Traits<problem_t>::unit_time * Physics_Traits<problem_t>::unit_time) / Physics_Traits<problem_t>::unit_temperature);
+		}
+	}();
 };
 
 template <typename problem_t>

src/hydro/NSCBC_inflow.hpp

@@ -67,7 +67,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_inflow_x1_lower(quokka::valarray<
 	const Real eta_5 = 2.;
 	const Real eta_6 = 2.;
 
-	const Real R = quokka::EOS_Traits<problem_t>::boltzmann_constant / quokka::EOS_Traits<problem_t>::mean_molecular_weight;
+	const Real R = quokka::EOS<problem_t>::boltzmann_constant_ / quokka::EOS_Traits<problem_t>::mean_molecular_weight;
 
 	// see SymPy notebook for derivation
 	quokka::valarray<Real, HydroSystem<problem_t>::nvar_> dQ_dx{};

src/physics_info.hpp

@@ -1,9 +1,13 @@
 #ifndef PHYSICS_INFO_HPP_ // NOLINT
 #define PHYSICS_INFO_HPP_
 
+#include "fundamental_constants.H"
 #include "physics_numVars.hpp"
 #include <AMReX.H>
 
+// enum for unit system, one of CGS, CONSTANTS, CUSTOM
+enum class UnitSystem { CGS, CONSTANTS, CUSTOM };
+
 // this struct is specialized by the user application code.
 template <typename problem_t> struct Physics_Traits {
 	// cell-centred
@@ -14,6 +18,15 @@ template <typename problem_t> struct Physics_Traits {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
+	static constexpr double boltzmann_constant = C::k_B;	    // Hydro, EOS
+	static constexpr double gravitational_constant = C::Gconst; // gravity
+	static constexpr double c_light = C::c_light;		    // radiation
+	static constexpr double radiation_constant = C::a_rad;	    // radiation
+	static constexpr double unit_length = 1.0;
+	static constexpr double unit_mass = 1.0;
+	static constexpr double unit_time = 1.0;
+	static constexpr double unit_temperature = 1.0;
 };
 
 // this struct stores the indices at which quantities start

src/problems/Advection/test_advection.cpp

@@ -34,6 +34,7 @@ template <> struct Physics_Traits<SawtoothProblem> {
 	static constexpr bool is_radiation_enabled = false;
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 AMREX_GPU_DEVICE void ComputeExactSolution(int i, int j, int k, int n, amrex::Array4<amrex::Real> const &exact_arr,

src/problems/Advection2D/test_advection2d.cpp

@@ -39,6 +39,7 @@ template <> struct Physics_Traits<SquareProblem> {
 	static constexpr bool is_radiation_enabled = false;
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto exactSolutionAtIndex(int i, int j, amrex::GpuArray<Real, AMREX_SPACEDIM> const &prob_lo,
@@ -106,7 +107,12 @@ template <> void AdvectionSimulation<SquareProblem>::ErrorEst(int lev, amrex::Ta
 
 			Real const del_x = (state(i + 1, j, k, n) - state(i - 1, j, k, n)) / (2.0 * dx[0]);
 			Real const del_y = (state(i, j + 1, k, n) - state(i, j - 1, k, n)) / (2.0 * dx[1]);
-			Real const gradient_indicator = std::sqrt(del_x * del_x + del_y * del_y) / rho;
+			Real gradient_indicator = NAN;
+			if (rho > 0) {
+				gradient_indicator = std::sqrt(del_x * del_x + del_y * del_y) / rho;
+			} else {
+				gradient_indicator = 1.0e100;
+			}
 
 			if (gradient_indicator > eta_threshold && rho >= rho_min) {
 				tag(i, j, k) = amrex::TagBox::SET;

src/problems/AdvectionSemiellipse/test_advection_semiellipse.cpp

@@ -33,6 +33,7 @@ template <> struct Physics_Traits<SemiellipseProblem> {
 	static constexpr bool is_radiation_enabled = false;
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 AMREX_GPU_DEVICE void ComputeExactSolution(int i, int j, int k, int n, amrex::Array4<amrex::Real> const &exact_arr,

src/problems/BinaryOrbitCIC/binary_orbit.cpp

@@ -35,7 +35,6 @@ template <> struct quokka::EOS_Traits<BinaryOrbit> {
 	static constexpr double gamma = 1.0;	       // isothermal
 	static constexpr double cs_isothermal = 1.3e7; // cm s^{-1}
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<BinaryOrbit> {
@@ -49,6 +48,7 @@ template <> struct Physics_Traits<BinaryOrbit> {
 	static constexpr int numMassScalars = 0;		     // number of mass scalars
 	static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars
 	static constexpr int nGroups = 1;			     // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> struct SimulationData<BinaryOrbit> {

src/problems/Cooling/test_cooling.cpp

@@ -27,7 +27,6 @@ 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 = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<CoolingTest> {
@@ -39,6 +38,7 @@ template <> struct Physics_Traits<CoolingTest> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> struct SimulationData<CoolingTest> {

src/problems/FCQuantities/test_fc_quantities.cpp

@@ -28,7 +28,6 @@ struct FCQuantities {
 template <> struct quokka::EOS_Traits<FCQuantities> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<FCQuantities> {
@@ -40,6 +39,7 @@ template <> struct Physics_Traits<FCQuantities> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = true;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 constexpr double rho0 = 1.0;					     // background density

src/problems/HydroBlast2D/test_hydro2d_blast.cpp

@@ -28,7 +28,6 @@ struct BlastProblem {
 template <> struct quokka::EOS_Traits<BlastProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<BlastProblem> {
@@ -40,6 +39,7 @@ template <> struct Physics_Traits<BlastProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<BlastProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/HydroBlast3D/test_hydro3d_blast.cpp

@@ -30,7 +30,6 @@ 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 = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<SedovProblem> {
@@ -46,6 +45,7 @@ template <> struct Physics_Traits<SedovProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 // declare global variables

src/problems/HydroContact/test_hydro_contact.cpp

@@ -24,7 +24,6 @@ struct ContactProblem {
 template <> struct quokka::EOS_Traits<ContactProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ContactProblem> {
@@ -36,6 +35,7 @@ template <> struct Physics_Traits<ContactProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 constexpr double v_contact = 0.0; // contact wave velocity

src/problems/HydroHighMach/test_hydro_highmach.cpp

@@ -31,7 +31,6 @@ struct HighMachProblem {
 template <> struct quokka::EOS_Traits<HighMachProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<HighMachProblem> {
@@ -43,6 +42,7 @@ template <> struct Physics_Traits<HighMachProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<HighMachProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/HydroKelvinHelmholz/test_hydro2d_kh.cpp

@@ -25,7 +25,6 @@ struct KelvinHelmholzProblem {
 template <> struct quokka::EOS_Traits<KelvinHelmholzProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<KelvinHelmholzProblem> {
@@ -41,6 +40,7 @@ template <> struct Physics_Traits<KelvinHelmholzProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<KelvinHelmholzProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/HydroLeblanc/test_hydro_leblanc.cpp

@@ -17,6 +17,7 @@
 
 #include "QuokkaSimulation.hpp"
 #include "hydro/hydro_system.hpp"
+#include "physics_info.hpp"
 #include "radiation/radiation_system.hpp"
 #include "test_hydro_leblanc.hpp"
 #include "util/ArrayUtil.hpp"
@@ -31,7 +32,6 @@ struct ShocktubeProblem {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = (5. / 3.);
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
@@ -43,6 +43,7 @@ template <> struct Physics_Traits<ShocktubeProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<ShocktubeProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/HydroQuirk/test_quirk.cpp

@@ -38,7 +38,6 @@ struct QuirkProblem {
 template <> struct quokka::EOS_Traits<QuirkProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<QuirkProblem> {
@@ -54,6 +53,7 @@ template <> struct Physics_Traits<QuirkProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 constexpr Real dl = 3.692;

src/problems/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp

@@ -22,7 +22,6 @@ struct RichtmeyerMeshkovProblem {
 template <> struct quokka::EOS_Traits<RichtmeyerMeshkovProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<RichtmeyerMeshkovProblem> {
@@ -38,6 +37,7 @@ template <> struct Physics_Traits<RichtmeyerMeshkovProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<RichtmeyerMeshkovProblem>::computeAfterTimestep()

src/problems/HydroSMS/test_hydro_sms.cpp

@@ -24,7 +24,6 @@ struct ShocktubeProblem {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
@@ -36,6 +35,7 @@ template <> struct Physics_Traits<ShocktubeProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 template <> void QuokkaSimulation<ShocktubeProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/HydroShocktube/test_hydro_shocktube.cpp

@@ -29,7 +29,6 @@ struct ShocktubeProblem {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
@@ -41,6 +40,7 @@ template <> struct Physics_Traits<ShocktubeProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 // left- and right- side shock states

src/problems/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp

@@ -34,7 +34,6 @@ template <> struct SimulationData<ShocktubeProblem> {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = C::m_u;
-	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
@@ -46,6 +45,7 @@ template <> struct Physics_Traits<ShocktubeProblem> {
 	// face-centred
 	static constexpr bool is_mhd_enabled = false;
 	static constexpr int nGroups = 1; // number of radiation groups
+	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 // left- and right- side shock states