Add support for unit systems: CGS, CONSTANTS, CUSTOM (#782)

### Description

In this PR we add the ability to allow the user to choose from one of
the three unit systems: `CGS`, `CUSTOM`, and `CONSTANTS`.

**CGS:**  
By setting `unit_system = UnitSystem::CGS`, the user opt to use the CGS
units and no constants definition are required. If `is_radiation_enabled
= true`, the user need to also specify `c_hat_over_c`. See examples in
HydroShuOsher and RadTube tests.
```C++
template <> struct Physics_Traits<TubeProblem> {
	static constexpr UnitSystem unit_system = UnitSystem::CGS;
	...
};

template <> struct RadSystem_Traits<TubeProblem> {
	static constexpr double c_hat_over_c = 0.1;
	...
};
```

**CONSTANTS:**  
By setting `unit_system = UnitSystem::CONSTANTS`, the user need to
define `boltzmann_constant` and `gravitational_constant` for hydro
simulations. If `is_radiation_enabled = true`, the user need to also
define `c_light`, `c_hat_over_c`, and `radiation_constant`. See examples
in the `RadhydroShock` test.
```C++
template <> struct Physics_Traits<problem_t> {
	static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS;
	static constexpr double gravitational_constant = 1.0;
	static constexpr double boltzmann_constant = 1.0;
	static constexpr double c_light = 1.0;
	static constexpr double radiation_constant = 1.0;
	...
};

template <> struct RadSystem_Traits<problem_t> {
	static constexpr double c_hat_over_c = 1.0;
	...
};
```

**CUSTOM:**
By setting `unit_system = UnitSystem::CUSTOM`, the user opt to use a
custom unit system. The values of the following variables need to be
given in cgs units: `unit_length, unit_mass, unit_time,
unit_temperature`. If `is_radiation_enabled = true`, the user need to
also define `c_hat` in code units. See examples in the `RadLineCooling`
and `RadhydroShockCGS` tests.
```C++
template <> struct Physics_Traits<problem_t> {
	// A custom unit system is used here to replicate a dimentionless unit system (c = k_B = a_rad = G = 1), for testing units conversion. This is similar to, but not exact, the Planck unit system.
	static constexpr UnitSystem unit_system = UnitSystem::CUSTOM;
	static constexpr double unit_length = 1.733039549e-33;
	static constexpr double unit_mass = 2.333695323e-05;
	static constexpr double unit_time = 5.780797690e-44;
	static constexpr double unit_temperature = 1.519155670e+32;
	...
};

template <> struct RadSystem_Traits<problem_t> {
	static constexpr double c_hat_over_c = 0.1;
	...
};
```

Other changes include:
1. The following variables are available in `AMRSimulation` class:
`unit_length, unit_mass, unit_time, unit_temperature`. They give the
unit length, mass, time, and temperature in CGS units (cm, g, s, and K).
2. The code will write the following variables into `metadata.yaml`
located inside every plotfile and checkpoint outputs: `unit_length,
unit_mass, unit_time, unit_temperature, c, c_hat, k_B, a_rad, G`. The
values of `unit_xxx` are given in CGS units and the values of `c, c_hat,
k_B, a_rad, G` are given in code units.
3. The gravitational constant is no longer a runtime parameter. It's
defined by the unit system the user specifies or by
`Physics_Traits::gravitational_constant`.

A typical `metadata.yaml` file looks like this:

```
a_rad: 7.5657313567241239e-15
c_hat: 402955198.55200708
c: 29979245800
G: 6.6742800000000007e-08
unit_length: 1
unit_time: 1
k_B: 1.3806488e-16
unit_temperature: 1
unit_mass: 1
```

Footnote: When `CONSTANTS` is chosen, I set `unit_xxx` to `NAN` for two
reasons. First, when radiation is turned off, the units are
unconstrained because it's not possible to derive `unit_xxx` from only
two constants, `boltzmann_constant` and `gravitational_constant`.
Second, the `CONSTANTS` unit system is only used for testing purpose and
we don't care about the physical dimensions of the variables.

### Related issues
Resolves #780 

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above). ✅ 2024-10-27
- [x] I have added a link to any related issues see (see above). ✅
2024-10-27
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
✅ 2024-10-27
- [x] I have added tests for any new physics that this PR adds to the
code. ✅ 2024-10-27
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR. ✅ 2024-10-27



### Appendix: Units conversion

Here I show how to convert between CGS, CONSTANTS, and CUSTOM unit
system.

Let the code units of length, mass, time, and temperature expressed in
CGS units be $u_l, u_m, u_t, u_T$. Let $G$, $k_B$, $c$, and $a_R$ be the
values of the gravitational constant, Boltzmann constant, speed of
light, and radiation constant in CGS units, respectively. Assume the
values of those constants in code units are $\hat{G}, \hat{k}_B,
\hat{c}, \hat{a}_R$. The following relations apply:

$$
\begin{aligned}
\bar{G} & \equiv G / \hat{G} = u_l^3 u_m^{-1} u_t^{-2}  \\
\bar{k}_B & \equiv k_B / \hat{k}_B = u_l^2 u_m u_t^{-2} u_T^{-1} \\
\bar{c} & \equiv c / \hat{c} = u_l  u_t^{-1} \\
\bar{a}_R & \equiv a_R/\hat{a}_R = u_l^{-1} u_m u_t^{-2} u_T^{-4}
\end{aligned}
$$

To convert from $\hat{G}$ etc into $u_l$ etc, we want to solve for $u_l$
etc from this set of four equations. By taking logarithm on both sides
of the equations, we can turn them into a system of linear equations:

$$
\begin{aligned}
& \left(\begin{array}{cccc}
3 & -1 & -2 & 0 \\
2 & 1 & -2 & -1 \\
1 & 0 & -1 & 0 \\
-1 & 1 & -2 & -4
\end{array}\right)\left(\begin{array}{l}
\log u_l \\
\log u_m \\
\log u_t \\
\log u_T
\end{array}\right)=\left(\begin{array}{l}
\log \bar{G} \\
\log \bar{k}_B \\
\log \bar{c} \\
\log \bar{a}_R
\end{array}\right)
\end{aligned}
$$

A simple matrix inversion gives

$$
\begin{aligned}
\left(\begin{array}{l}
\log u_l \\
\log u_m \\
\log u_t \\
\log u_T
\end{array}\right)=\left(\begin{array}{cccc}
1 / 2 & 2 / 3 & -2 & -1 / 6 \\
-1 / 2 & 2 / 3 & 0 & -1 / 6 \\
1 / 2 & 2 / 3 & -3 & -1 / 6 \\
-1 / 2 & -1 / 3 & 2 & -1 / 6
\end{array}\right)\left(\begin{array}{c}
\log \bar{G} \\
\log \bar{k}_B \\
\log \bar{c} \\
\log \bar{a}_R
\end{array}\right)
\end{aligned}
$$

or, expressing explicitly,

$$
\begin{aligned}
u_l&=\bar{G}^{1/2} \bar{c}^{-2} \bar{d} \\
u_r&=\bar{G}^{-1/2} \bar{d} \\
u_t&=\bar{G}^{1/2} \bar{c}^{-3} \bar{d} \\
u_T&=\bar{G}^{-1/2} \bar{c}^2\left(\bar{k}^2 \bar{a}_R\right)^{-1 / 6}
\end{aligned}
$$

where $\bar{d} \equiv \left(\bar{k}^4 / \bar{a}_R\right)^{1 / 6}$. These
equations are used to derive `unit_length` etc in the `RadLineCooling`
test.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

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