Add support for unit systems: CGS, CONSTANTS, CUSTOM

chongchonghe commented 1 month ago

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.

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.

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.

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:

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).
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.
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. ✅ 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.

BenWibking commented 4 weeks ago

I think it would make sense to move all of the constants to Physics_Traits, so that the unit system is always defined by Physics_Traits and nothing else is needed.