From 4fa830c2efac48aa56fb2083f137cba5cf07cf92 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 4 Nov 2024 10:59:51 +1100 Subject: [PATCH] Add support for unit systems: CGS, CONSTANTS, CUSTOM (#782) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ### 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 { static constexpr UnitSystem unit_system = UnitSystem::CGS; ... }; template <> struct RadSystem_Traits { 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 { 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 { 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 { // 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 { 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 | 47 ++++++++++++++ src/hydro/EOS.hpp | 19 ++++-- src/hydro/NSCBC_inflow.hpp | 2 +- src/physics_info.hpp | 13 ++++ src/problems/Advection/test_advection.cpp | 1 + src/problems/Advection2D/test_advection2d.cpp | 8 ++- .../test_advection_semiellipse.cpp | 1 + src/problems/BinaryOrbitCIC/binary_orbit.cpp | 2 +- src/problems/Cooling/test_cooling.cpp | 2 +- .../FCQuantities/test_fc_quantities.cpp | 2 +- .../HydroBlast2D/test_hydro2d_blast.cpp | 2 +- .../HydroBlast3D/test_hydro3d_blast.cpp | 2 +- .../HydroContact/test_hydro_contact.cpp | 2 +- .../HydroHighMach/test_hydro_highmach.cpp | 2 +- .../HydroKelvinHelmholz/test_hydro2d_kh.cpp | 2 +- .../HydroLeblanc/test_hydro_leblanc.cpp | 3 +- src/problems/HydroQuirk/test_quirk.cpp | 2 +- .../test_hydro2d_rm.cpp | 2 +- src/problems/HydroSMS/test_hydro_sms.cpp | 2 +- .../HydroShocktube/test_hydro_shocktube.cpp | 2 +- .../test_hydro_shocktube_cma.cpp | 2 +- .../HydroShuOsher/test_hydro_shuosher.cpp | 2 +- .../HydroVacuum/test_hydro_vacuum.cpp | 2 +- src/problems/HydroWave/test_hydro_wave.cpp | 2 +- src/problems/NSCBC/channel.cpp | 2 +- src/problems/NSCBC/vortex.cpp | 4 +- src/problems/ODEIntegration/test_ode.cpp | 1 - src/problems/PassiveScalar/test_scalars.cpp | 2 +- src/problems/PopIII/popiii.cpp | 1 + .../PrimordialChem/test_primordial_chem.cpp | 1 + src/problems/RadBeam/test_radiation_beam.cpp | 6 +- src/problems/RadDust/test_rad_dust.cpp | 10 +-- src/problems/RadDustMG/test_rad_dust_MG.cpp | 12 ++-- .../RadForce/test_radiation_force.cpp | 6 +- .../RadLineCooling/test_rad_line_cooling.cpp | 19 ++++-- .../test_rad_line_cooling_MG.cpp | 16 ++--- .../RadMarshak/test_radiation_marshak.cpp | 9 ++- .../test_radiation_marshak_asymptotic.cpp | 14 +--- .../test_radiation_marshak_cgs.cpp | 8 +-- .../test_radiation_marshak_dust.cpp | 17 ++--- .../test_radiation_marshak_dust_and_PE.cpp | 13 ++-- .../test_radiation_marshak_Vaytet.cpp | 6 +- .../test_radiation_matter_coupling.cpp | 8 +-- .../test_radiation_matter_coupling_rsla.cpp | 11 ++-- .../RadPulse/test_radiation_pulse.cpp | 10 +-- .../RadShadow/test_radiation_shadow.cpp | 13 ++-- .../RadStreaming/test_radiation_streaming.cpp | 10 +-- .../test_radiation_streaming_y.cpp | 12 ++-- .../RadSuOlson/test_radiation_SuOlson.cpp | 11 ++-- .../RadTophat/test_radiation_tophat.cpp | 10 ++- src/problems/RadTube/test_radiation_tube.cpp | 6 +- src/problems/RadhydroBB/test_radhydro_bb.cpp | 10 +-- .../RadhydroPulse/test_radhydro_pulse.cpp | 12 ++-- .../test_radhydro_pulse_dyn.cpp | 15 ++--- .../test_radhydro_pulse_grey.cpp | 15 ++--- .../test_radhydro_pulse_MG_const_kappa.cpp | 12 ++-- .../test_radhydro_pulse_MG_int.cpp | 12 ++-- .../RadhydroShell/test_radhydro_shell.cpp | 18 +++--- .../RadhydroShock/test_radhydro_shock.cpp | 12 ++-- .../test_radhydro_shock_cgs.cpp | 11 ++-- .../test_radhydro_shock_multigroup.cpp | 6 +- .../test_radhydro_uniform_advecting.cpp | 12 ++-- src/problems/RandomBlast/blast.cpp | 6 +- .../RayleighTaylor2D/test_hydro2d_rt.cpp | 2 +- .../RayleighTaylor3D/test_hydro3d_rt.cpp | 2 +- src/problems/ShockCloud/cloud.cpp | 2 +- .../SphericalCollapse/spherical_collapse.cpp | 4 +- src/problems/StarCluster/star_cluster.cpp | 4 +- src/radiation/radiation_system.hpp | 48 ++++++++++++-- src/simulation.hpp | 64 +++++++++++++++++-- tests/RadStreamingY.in | 2 +- tests/RadTube.in | 4 +- tests/SphericalCollapse.in | 2 - tests/radshock.in | 4 +- tests/radshock_dimensionless.in | 3 +- 75 files changed, 400 insertions(+), 246 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 42c7cc462..f6fff3738 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -321,6 +321,53 @@ template void QuokkaSimulation::defineComponentN } } +// initialize metadata +template void AMRSimulation::initializeSimulationMetadata() +{ + if constexpr (Physics_Traits::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::boltzmann_constant; + simulationMetadata_["G"] = Physics_Traits::gravitational_constant; + if constexpr (Physics_Traits::is_radiation_enabled) { + simulationMetadata_["c"] = Physics_Traits::c_light; + simulationMetadata_["c_hat"] = Physics_Traits::c_light * RadSystem_Traits::c_hat_over_c; + simulationMetadata_["a_rad"] = Physics_Traits::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::unit_system == UnitSystem::CGS) { + k_B = C::k_B; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + // Have to do a conversion because EOS class is not accessible here + k_B = C::k_B / + (Physics_Traits::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_mass / + (Physics_Traits::unit_time * Physics_Traits::unit_time) / Physics_Traits::unit_temperature); + } + simulationMetadata_["k_B"] = k_B; + simulationMetadata_["G"] = Gconst_; + if constexpr (Physics_Traits::is_radiation_enabled) { + simulationMetadata_["c"] = RadSystem::c_light_; + simulationMetadata_["c_hat"] = RadSystem::c_hat_; + simulationMetadata_["a_rad"] = RadSystem::radiation_constant_; + } + } +} + template auto QuokkaSimulation::getScalarVariableNames() -> std::vector { // return vector of names for the passive scalars diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index a8c5a5e60..a54cd0162 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -38,6 +38,9 @@ template struct EOS_Traits { template class EOS { + private: + static constexpr amrex::Real gamma_ = EOS_Traits::gamma; + static constexpr amrex::Real mean_molecular_weight_ = EOS_Traits::mean_molecular_weight; public: static constexpr int nmscalars_ = Physics_Traits::numMassScalars; @@ -65,10 +68,18 @@ template class EOS ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure, std::optional> const &massScalars = {}) -> amrex::Real; - private: - static constexpr amrex::Real gamma_ = EOS_Traits::gamma; - static constexpr amrex::Real boltzmann_constant_ = EOS_Traits::boltzmann_constant; - static constexpr amrex::Real mean_molecular_weight_ = EOS_Traits::mean_molecular_weight; + static constexpr amrex::Real boltzmann_constant_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return C::k_B; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return Physics_Traits::boltzmann_constant; + } else if constexpr (Physics_Traits::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::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_mass / + (Physics_Traits::unit_time * Physics_Traits::unit_time) / Physics_Traits::unit_temperature); + } + }(); }; template diff --git a/src/hydro/NSCBC_inflow.hpp b/src/hydro/NSCBC_inflow.hpp index 19b308fae..a46e2ba13 100644 --- a/src/hydro/NSCBC_inflow.hpp +++ b/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::boltzmann_constant / quokka::EOS_Traits::mean_molecular_weight; + const Real R = quokka::EOS::boltzmann_constant_ / quokka::EOS_Traits::mean_molecular_weight; // see SymPy notebook for derivation quokka::valarray::nvar_> dQ_dx{}; diff --git a/src/physics_info.hpp b/src/physics_info.hpp index e41ff1fda..5597d0edb 100644 --- a/src/physics_info.hpp +++ b/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 +// 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 struct Physics_Traits { // cell-centred @@ -14,6 +18,15 @@ template 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 diff --git a/src/problems/Advection/test_advection.cpp b/src/problems/Advection/test_advection.cpp index 460c028d1..0a7b68b71 100644 --- a/src/problems/Advection/test_advection.cpp +++ b/src/problems/Advection/test_advection.cpp @@ -34,6 +34,7 @@ template <> struct Physics_Traits { 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 const &exact_arr, diff --git a/src/problems/Advection2D/test_advection2d.cpp b/src/problems/Advection2D/test_advection2d.cpp index 4b9f0b1d9..7810d1831 100644 --- a/src/problems/Advection2D/test_advection2d.cpp +++ b/src/problems/Advection2D/test_advection2d.cpp @@ -39,6 +39,7 @@ template <> struct Physics_Traits { 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 const &prob_lo, @@ -106,7 +107,12 @@ template <> void AdvectionSimulation::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; diff --git a/src/problems/AdvectionSemiellipse/test_advection_semiellipse.cpp b/src/problems/AdvectionSemiellipse/test_advection_semiellipse.cpp index e2bd8250f..0f7187994 100644 --- a/src/problems/AdvectionSemiellipse/test_advection_semiellipse.cpp +++ b/src/problems/AdvectionSemiellipse/test_advection_semiellipse.cpp @@ -33,6 +33,7 @@ template <> struct Physics_Traits { 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 const &exact_arr, diff --git a/src/problems/BinaryOrbitCIC/binary_orbit.cpp b/src/problems/BinaryOrbitCIC/binary_orbit.cpp index 7878d0d0f..5d4ce3c98 100644 --- a/src/problems/BinaryOrbitCIC/binary_orbit.cpp +++ b/src/problems/BinaryOrbitCIC/binary_orbit.cpp @@ -35,7 +35,6 @@ template <> struct quokka::EOS_Traits { 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 { @@ -49,6 +48,7 @@ template <> struct Physics_Traits { 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 { diff --git a/src/problems/Cooling/test_cooling.cpp b/src/problems/Cooling/test_cooling.cpp index 87e058aa7..0b4b3e5fd 100644 --- a/src/problems/Cooling/test_cooling.cpp +++ b/src/problems/Cooling/test_cooling.cpp @@ -27,7 +27,6 @@ constexpr double seconds_in_year = 3.154e7; template <> struct quokka::EOS_Traits { 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 { @@ -39,6 +38,7 @@ template <> 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; }; template <> struct SimulationData { diff --git a/src/problems/FCQuantities/test_fc_quantities.cpp b/src/problems/FCQuantities/test_fc_quantities.cpp index eb8820960..9fb7371f6 100644 --- a/src/problems/FCQuantities/test_fc_quantities.cpp +++ b/src/problems/FCQuantities/test_fc_quantities.cpp @@ -28,7 +28,6 @@ struct FCQuantities { template <> struct quokka::EOS_Traits { 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 { @@ -40,6 +39,7 @@ template <> struct Physics_Traits { // 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 diff --git a/src/problems/HydroBlast2D/test_hydro2d_blast.cpp b/src/problems/HydroBlast2D/test_hydro2d_blast.cpp index bc4548f1a..115fd3206 100644 --- a/src/problems/HydroBlast2D/test_hydro2d_blast.cpp +++ b/src/problems/HydroBlast2D/test_hydro2d_blast.cpp @@ -28,7 +28,6 @@ struct BlastProblem { template <> struct quokka::EOS_Traits { 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 { @@ -40,6 +39,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroBlast3D/test_hydro3d_blast.cpp b/src/problems/HydroBlast3D/test_hydro3d_blast.cpp index 9d601462b..49f742f27 100644 --- a/src/problems/HydroBlast3D/test_hydro3d_blast.cpp +++ b/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 { 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 { @@ -46,6 +45,7 @@ template <> 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; }; // declare global variables diff --git a/src/problems/HydroContact/test_hydro_contact.cpp b/src/problems/HydroContact/test_hydro_contact.cpp index 47bf32f20..b78b853dc 100644 --- a/src/problems/HydroContact/test_hydro_contact.cpp +++ b/src/problems/HydroContact/test_hydro_contact.cpp @@ -24,7 +24,6 @@ struct ContactProblem { template <> struct quokka::EOS_Traits { 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 { @@ -36,6 +35,7 @@ template <> 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; }; constexpr double v_contact = 0.0; // contact wave velocity diff --git a/src/problems/HydroHighMach/test_hydro_highmach.cpp b/src/problems/HydroHighMach/test_hydro_highmach.cpp index f70c05355..a8ecf6446 100644 --- a/src/problems/HydroHighMach/test_hydro_highmach.cpp +++ b/src/problems/HydroHighMach/test_hydro_highmach.cpp @@ -31,7 +31,6 @@ struct HighMachProblem { template <> struct quokka::EOS_Traits { 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 { @@ -43,6 +42,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroKelvinHelmholz/test_hydro2d_kh.cpp b/src/problems/HydroKelvinHelmholz/test_hydro2d_kh.cpp index e292bdce9..fa6cd4dac 100644 --- a/src/problems/HydroKelvinHelmholz/test_hydro2d_kh.cpp +++ b/src/problems/HydroKelvinHelmholz/test_hydro2d_kh.cpp @@ -25,7 +25,6 @@ struct KelvinHelmholzProblem { template <> struct quokka::EOS_Traits { 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 { @@ -41,6 +40,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroLeblanc/test_hydro_leblanc.cpp b/src/problems/HydroLeblanc/test_hydro_leblanc.cpp index e9c9199a4..c17ce9e44 100644 --- a/src/problems/HydroLeblanc/test_hydro_leblanc.cpp +++ b/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 { 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 { @@ -43,6 +43,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroQuirk/test_quirk.cpp b/src/problems/HydroQuirk/test_quirk.cpp index a60ada952..64ea765b4 100644 --- a/src/problems/HydroQuirk/test_quirk.cpp +++ b/src/problems/HydroQuirk/test_quirk.cpp @@ -38,7 +38,6 @@ struct QuirkProblem { template <> struct quokka::EOS_Traits { 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 { @@ -54,6 +53,7 @@ template <> 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; }; constexpr Real dl = 3.692; diff --git a/src/problems/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp b/src/problems/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp index 09b48ba60..1ec363919 100644 --- a/src/problems/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp +++ b/src/problems/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp @@ -22,7 +22,6 @@ struct RichtmeyerMeshkovProblem { template <> struct quokka::EOS_Traits { 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 { @@ -38,6 +37,7 @@ template <> 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; }; template <> void QuokkaSimulation::computeAfterTimestep() diff --git a/src/problems/HydroSMS/test_hydro_sms.cpp b/src/problems/HydroSMS/test_hydro_sms.cpp index 7d2baeafc..eb9bb32e1 100644 --- a/src/problems/HydroSMS/test_hydro_sms.cpp +++ b/src/problems/HydroSMS/test_hydro_sms.cpp @@ -24,7 +24,6 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { 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 { @@ -36,6 +35,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroShocktube/test_hydro_shocktube.cpp b/src/problems/HydroShocktube/test_hydro_shocktube.cpp index cc19860de..a9371df45 100644 --- a/src/problems/HydroShocktube/test_hydro_shocktube.cpp +++ b/src/problems/HydroShocktube/test_hydro_shocktube.cpp @@ -29,7 +29,6 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { 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 { @@ -41,6 +40,7 @@ template <> 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; }; // left- and right- side shock states diff --git a/src/problems/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp b/src/problems/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp index e5495455e..802497595 100644 --- a/src/problems/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp +++ b/src/problems/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp @@ -34,7 +34,6 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { 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 { @@ -46,6 +45,7 @@ template <> 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; }; // left- and right- side shock states diff --git a/src/problems/HydroShuOsher/test_hydro_shuosher.cpp b/src/problems/HydroShuOsher/test_hydro_shuosher.cpp index ef847a179..f558a7aa1 100644 --- a/src/problems/HydroShuOsher/test_hydro_shuosher.cpp +++ b/src/problems/HydroShuOsher/test_hydro_shuosher.cpp @@ -26,7 +26,6 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { 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 { @@ -38,6 +37,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroVacuum/test_hydro_vacuum.cpp b/src/problems/HydroVacuum/test_hydro_vacuum.cpp index 14774aa05..28bdcfda0 100644 --- a/src/problems/HydroVacuum/test_hydro_vacuum.cpp +++ b/src/problems/HydroVacuum/test_hydro_vacuum.cpp @@ -28,7 +28,6 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { 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 { @@ -40,6 +39,7 @@ template <> 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; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/HydroWave/test_hydro_wave.cpp b/src/problems/HydroWave/test_hydro_wave.cpp index 00d329e4b..b11eeb5fa 100644 --- a/src/problems/HydroWave/test_hydro_wave.cpp +++ b/src/problems/HydroWave/test_hydro_wave.cpp @@ -24,7 +24,6 @@ struct WaveProblem { template <> struct quokka::EOS_Traits { 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 { @@ -36,6 +35,7 @@ template <> 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; }; constexpr double rho0 = 1.0; // background density diff --git a/src/problems/NSCBC/channel.cpp b/src/problems/NSCBC/channel.cpp index 59b6d5e7d..21a36834d 100644 --- a/src/problems/NSCBC/channel.cpp +++ b/src/problems/NSCBC/channel.cpp @@ -51,7 +51,6 @@ struct Channel { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.1; static constexpr double mean_molecular_weight = 28.96 * C::m_u; // air - static constexpr double boltzmann_constant = C::k_B; }; template <> struct Physics_Traits { @@ -61,6 +60,7 @@ template <> struct Physics_Traits { static constexpr int numPassiveScalars = numMassScalars + 1; // number of passive scalars static constexpr bool is_radiation_enabled = false; static constexpr int nGroups = 1; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; // global variables needed for Dirichlet boundary condition and initial conditions diff --git a/src/problems/NSCBC/vortex.cpp b/src/problems/NSCBC/vortex.cpp index c5742f4a8..1fa9c79ba 100644 --- a/src/problems/NSCBC/vortex.cpp +++ b/src/problems/NSCBC/vortex.cpp @@ -48,7 +48,6 @@ struct Vortex { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; static constexpr double mean_molecular_weight = 28.96 * C::m_u; // air - static constexpr double boltzmann_constant = C::k_B; }; template <> struct Physics_Traits { @@ -58,6 +57,7 @@ template <> struct Physics_Traits { static constexpr int numPassiveScalars = numMassScalars + 1; // number of passive scalars static constexpr bool is_radiation_enabled = false; static constexpr int nGroups = 1; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; // global variables needed for Dirichlet boundary condition and initial conditions @@ -83,7 +83,7 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::gr const amrex::GpuArray prob_hi = grid_elem.prob_hi_; constexpr Real gamma = quokka::EOS_Traits::gamma; - constexpr Real R = quokka::EOS_Traits::boltzmann_constant / quokka::EOS_Traits::mean_molecular_weight; + constexpr Real R = C::k_B / quokka::EOS_Traits::mean_molecular_weight; const Real c = std::sqrt(gamma * R * T_ref); const Real G = ::G_vortex; diff --git a/src/problems/ODEIntegration/test_ode.cpp b/src/problems/ODEIntegration/test_ode.cpp index 947d57f98..5fa885101 100644 --- a/src/problems/ODEIntegration/test_ode.cpp +++ b/src/problems/ODEIntegration/test_ode.cpp @@ -19,7 +19,6 @@ constexpr double rho0 = 0.01 * (C::m_p + C::m_e); // g cm^-3 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; diff --git a/src/problems/PassiveScalar/test_scalars.cpp b/src/problems/PassiveScalar/test_scalars.cpp index 7861834a2..bdc7addd5 100644 --- a/src/problems/PassiveScalar/test_scalars.cpp +++ b/src/problems/PassiveScalar/test_scalars.cpp @@ -26,7 +26,6 @@ struct ScalarProblem { template <> struct quokka::EOS_Traits { 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 { @@ -38,6 +37,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; // number of radiation groups; has to be defined even though radiation is disabled + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; constexpr double v_contact = 2.0; // contact wave velocity diff --git a/src/problems/PopIII/popiii.cpp b/src/problems/PopIII/popiii.cpp index d1ff9d069..39e3cb90d 100644 --- a/src/problems/PopIII/popiii.cpp +++ b/src/problems/PopIII/popiii.cpp @@ -48,6 +48,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct SimulationData { diff --git a/src/problems/PrimordialChem/test_primordial_chem.cpp b/src/problems/PrimordialChem/test_primordial_chem.cpp index 9d9d83aef..b0f363c2f 100644 --- a/src/problems/PrimordialChem/test_primordial_chem.cpp +++ b/src/problems/PrimordialChem/test_primordial_chem.cpp @@ -49,6 +49,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct SimulationData { diff --git a/src/problems/RadBeam/test_radiation_beam.cpp b/src/problems/RadBeam/test_radiation_beam.cpp index df9cc02a8..f0f6d55e2 100644 --- a/src/problems/RadBeam/test_radiation_beam.cpp +++ b/src/problems/RadBeam/test_radiation_beam.cpp @@ -28,14 +28,11 @@ constexpr double c = 2.99792458e10; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -49,6 +46,7 @@ template <> 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; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index 52b460dfc..f736002a6 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -40,14 +40,11 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; }; @@ -67,6 +64,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = k_B; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index cb504ec2a..98c68480f 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -16,7 +16,6 @@ struct DustProblem { constexpr int beta_order_ = 1; // order of beta in the radiation four-force constexpr double c = 1.0e8; -constexpr double chat = c; constexpr double v0 = 0.0; constexpr double chi0 = 10000.0; @@ -24,7 +23,6 @@ constexpr double T0 = 1.0; constexpr double rho0 = 1.0; constexpr double a_rad = 1.0; constexpr double mu = 1.0; -constexpr double k_B = 1.0; constexpr double max_time = 1.0e-5; constexpr double delta_time = 1.0e-8; @@ -40,7 +38,6 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -53,12 +50,15 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 4; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = 1.0; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; static constexpr double energy_unit = 1.; diff --git a/src/problems/RadForce/test_radiation_force.cpp b/src/problems/RadForce/test_radiation_force.cpp index d5934ebb4..61fd659b1 100644 --- a/src/problems/RadForce/test_radiation_force.cpp +++ b/src/problems/RadForce/test_radiation_force.cpp @@ -47,7 +47,6 @@ constexpr double Lx = (a0 * a0) / g0; // cm template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = gamma_gas; static constexpr double cs_isothermal = a0; // only used when gamma = 1 }; @@ -62,12 +61,11 @@ template <> struct Physics_Traits { static constexpr bool is_mhd_enabled = false; // number of radiation groups static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = 10. * (Mach1 * a0); - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 10. * (Mach1 * a0) / C::c_light; static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; static constexpr int beta_order = 1; diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index c8f78e675..4f12b5012 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -44,7 +44,6 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -57,15 +56,25 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + // A custom unit system is used here to replicate a dimentionless unit system (c = k_B = a_rad = G = 1), for testing units conversion + 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; + // Equivalently, set + // static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + // static constexpr double boltzmann_constant = k_B; + // static constexpr double gravitational_constant = 1.0; + // static constexpr double c_light = c; + // static constexpr double radiation_constant = a_rad; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 0; static constexpr double energy_unit = nu_unit; + static constexpr double c_hat_over_c = chat / c; }; template <> struct ISM_Traits { @@ -202,7 +211,7 @@ auto problem_main() -> int sim.initDt_ = the_dt; sim.maxDt_ = the_dt; sim.maxTimesteps_ = max_timesteps; - sim.plotfileInterval_ = -1; + sim.plotfileInterval_ = 10000000; // initialize sim.setInitialConditions(); diff --git a/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp b/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp index cdd520ed9..629bfaeb8 100644 --- a/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp +++ b/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp @@ -19,8 +19,7 @@ struct CoolingProblemMG { constexpr int n_groups_ = 4; constexpr amrex::GpuArray rad_boundaries_ = {1.00000000e-03, 1.77827941e-02, 3.16227766e-01, 5.62341325e+00, 1.00000000e+02}; -constexpr double c = 1.0; -constexpr double chat = c; +constexpr double chat_over_c = 1.0; constexpr double v0 = 0.0; constexpr double kappa0 = 0.0; @@ -29,7 +28,6 @@ constexpr double rho0 = 1.0; // matter density constexpr double a_rad = 1.0; constexpr double mu = 1.5; // mean molecular weight; so that C_V = 1.0 constexpr double C_V = 1.0; -constexpr double k_B = 1.0; constexpr double nu_unit = 1.0; constexpr double Erad_bar = a_rad * T0 * T0 * T0 * T0; @@ -50,7 +48,6 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -63,12 +60,15 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_groups_; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = 1.0; + static constexpr double radiation_constant = a_rad; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat_over_c; static constexpr double Erad_floor = Erad_floor_; static constexpr int beta_order = 0; static constexpr double energy_unit = nu_unit; @@ -235,7 +235,7 @@ auto problem_main() -> int std::exp(-cooling_rate * t_exact) * (cooling_rate * T0 - heating_rate_ + heating_rate_ * std::exp(cooling_rate * t_exact)) / cooling_rate; const double T_exact_solution = Egas_exact_solution / C_V; Tgas_exact_vec.push_back(T_exact_solution); - const double Erad_line_exact_solution = -(Egas_exact_solution - C_V * T0 - heating_rate_ * t_exact) * (chat / c); + const double Erad_line_exact_solution = -(Egas_exact_solution - C_V * T0 - heating_rate_ * t_exact) * chat_over_c; Erad_line_exact_vec.push_back(Erad_line_exact_solution); t_exact_vec.push_back(t_exact); t_exact += t_end / N_dt; diff --git a/src/problems/RadMarshak/test_radiation_marshak.cpp b/src/problems/RadMarshak/test_radiation_marshak.cpp index 81e031bbf..89a0e4171 100644 --- a/src/problems/RadMarshak/test_radiation_marshak.cpp +++ b/src/problems/RadMarshak/test_radiation_marshak.cpp @@ -32,13 +32,11 @@ constexpr double T_initial = 1.0e-2; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 1.0; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = c; + static constexpr double c_hat_over_c = 1.0; static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; @@ -53,6 +51,11 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index e32e29abd..f39c27df9 100644 --- a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -29,14 +29,11 @@ constexpr double Erad_floor_ = a_rad * T_initial * T_initial * T_initial * T_ini template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = Erad_floor_; static constexpr int beta_order = 0; }; @@ -50,6 +47,7 @@ template <> 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; }; template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> amrex::Real @@ -206,14 +204,6 @@ auto problem_main() -> int // const int nx = 60; // [18 == matches resolution of McClarren & Lowrie (2008)] // const double Lx = 0.66; // cm - // Problem initialization - std::cout << "radiation constant (code units) = " << RadSystem_Traits::radiation_constant << "\n"; - std::cout << "c_light (code units) = " << RadSystem_Traits::c_light << "\n"; - std::cout << "rho * c_v = " << rho0 * c_v << "\n"; - std::cout << "initial_dt = " << initial_dt << "\n"; - std::cout << "max_dt = " << max_dt << "\n"; - std::cout << "max_time = " << max_time << std::endl; - constexpr int nvars = RadSystem::nvar_; amrex::Vector BCs_cc(nvars); for (int n = 0; n < nvars; ++n) { diff --git a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp index 818639850..a3edcdd76 100644 --- a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp +++ b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp @@ -27,21 +27,18 @@ constexpr double eps_SuOlson = 1.0; // dimensionless constexpr double kappa = 577.0; // g cm^-2 (opacity) constexpr double rho0 = 10.0; // g cm^-3 (matter density) constexpr double T_hohlraum = 3.481334e6; // K -constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4 +constexpr double a_rad = C::a_rad; // erg cm^-3 K^-4 // constexpr double c = 2.99792458e10; // cm s^-1 constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson; constexpr double T_initial = 1.0e4; // K template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -55,6 +52,7 @@ template <> 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; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp index 7aab2e713..e5f5d8b69 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -19,8 +19,7 @@ struct MarshakProblem { AMREX_GPU_MANAGED double kappa1 = NAN; // dust opacity at IR AMREX_GPU_MANAGED double kappa2 = NAN; // dust opacity at FUV -constexpr double c = 1.0; // speed of light -constexpr double chat = 1.0; // reduced speed of light +constexpr double c = 1.0; // speed of light constexpr double rho0 = 1.0; constexpr double CV = 1.0; constexpr double mu = 1.5 / CV; // mean molecular weight @@ -42,7 +41,6 @@ static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_ template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; @@ -55,12 +53,15 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_group_; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = 1.0; + static constexpr double radiation_constant = a_rad; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 0; static constexpr double energy_unit = 1.0; @@ -151,8 +152,8 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect // const auto Erads = RadSystem::ComputeThermalRadiation(T_rad_L, radBoundaries_); quokka::valarray const Erads = {erad_floor, EradL}; - const double c_light = c; - const auto Frads = Erads * c_light; + const double c_device = c; + const auto Frads = Erads * c_device; if (i < lo[0]) { // streaming inflow boundary diff --git a/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp index 991b05d2d..6d4096ac5 100644 --- a/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp +++ b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp @@ -24,8 +24,7 @@ constexpr bool dust_on = 1; constexpr bool PE_on = 1; constexpr double gas_dust_coupling_threshold_ = 1.0e-4; -constexpr double c = 1.0; // speed of light -constexpr double chat = 1.0; // reduced speed of light +constexpr double c = 1.0; // speed of light constexpr double rho0 = 1.0; constexpr double CV = 1.0; constexpr double mu = 1.5 / CV; // mean molecular weight @@ -41,7 +40,6 @@ static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_ template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; @@ -54,12 +52,15 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_group_; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; static constexpr double energy_unit = 1.0; diff --git a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp index 5e9a75173..e9d39fe24 100644 --- a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp +++ b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp @@ -100,7 +100,6 @@ constexpr double Erad_floor_ = a_rad * T_initial * T_initial * T_initial * T_ini template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; @@ -113,12 +112,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_groups_; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = Erad_floor_; static constexpr int beta_order = 0; static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz diff --git a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp index 11be6e3d3..5d09df71b 100644 --- a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp +++ b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp @@ -21,7 +21,7 @@ struct CouplingProblem { // Su & Olson (1997) test problem constexpr double eps_SuOlson = 1.0; -constexpr double a_rad = 7.5646e-15; // cgs +constexpr double a_rad = C::a_rad; // cgs constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson; template <> struct SimulationData { @@ -32,14 +32,11 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -53,6 +50,7 @@ template <> 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; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp index 97d5a27ac..39e4029c7 100644 --- a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp +++ b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp @@ -19,11 +19,12 @@ struct CouplingProblem { }; // dummy type to allow compile-type polymorphism via template specialization // constexpr double c = 2.99792458e10; // cgs -constexpr double c_rsla = 0.1 * c_light_cgs_; +constexpr double chat_over_c = 0.1; +constexpr double c_rsla = chat_over_c * C::c_light; // Su & Olson (1997) test problem constexpr double eps_SuOlson = 1.0; -constexpr double a_rad = 7.5646e-15; // cgs +constexpr double a_rad = C::a_rad; // cgs constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson; template <> struct SimulationData { @@ -34,14 +35,11 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_rsla; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = chat_over_c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -55,6 +53,7 @@ template <> 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; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadPulse/test_radiation_pulse.cpp b/src/problems/RadPulse/test_radiation_pulse.cpp index cde6a25cc..da1f7586d 100644 --- a/src/problems/RadPulse/test_radiation_pulse.cpp +++ b/src/problems/RadPulse/test_radiation_pulse.cpp @@ -30,14 +30,11 @@ constexpr double initial_time = 1.0e-8; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 1.0; - static constexpr double boltzmann_constant = (2. / 3.); static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 0; }; @@ -51,6 +48,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = (2. / 3.); + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; AMREX_GPU_HOST_DEVICE diff --git a/src/problems/RadShadow/test_radiation_shadow.cpp b/src/problems/RadShadow/test_radiation_shadow.cpp index f4e583ac2..d99e1d41e 100644 --- a/src/problems/RadShadow/test_radiation_shadow.cpp +++ b/src/problems/RadShadow/test_radiation_shadow.cpp @@ -26,19 +26,16 @@ constexpr double rho_clump = 1.0; // g cm^-3 (matter density) constexpr double T_hohlraum = 1740.; // K constexpr double T_initial = 290.; // K -constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4 -constexpr double c = 2.99792458e10; // cm s^-1 +constexpr double a_rad = C::a_rad; // erg cm^-3 K^-4 +constexpr double c = C::c_light; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 10. * C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = c; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -50,6 +47,7 @@ template <> struct Physics_Traits { static constexpr bool is_radiation_enabled = true; static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real @@ -228,8 +226,7 @@ auto problem_main() -> int // Print radiation epsilon ("stiffness parameter" from Su & Olson). // (if epsilon is smaller than tolerance, there can be unphysical radiation shocks.) const auto dt_CFL = CFL_number * std::min(Lx / nx, Ly / ny) / c; - const auto c_v = quokka::EOS_Traits::boltzmann_constant / - (quokka::EOS_Traits::mean_molecular_weight * (quokka::EOS_Traits::gamma - 1.0)); + const auto c_v = C::k_B / (quokka::EOS_Traits::mean_molecular_weight * (quokka::EOS_Traits::gamma - 1.0)); const auto epsilon = 4.0 * a_rad * (T_initial * T_initial * T_initial) * sigma0 * (c * dt_CFL) / c_v; amrex::Print() << "radiation epsilon (stiffness parameter) = " << epsilon << "\n"; diff --git a/src/problems/RadStreaming/test_radiation_streaming.cpp b/src/problems/RadStreaming/test_radiation_streaming.cpp index 1746e9e5b..0df0c7676 100644 --- a/src/problems/RadStreaming/test_radiation_streaming.cpp +++ b/src/problems/RadStreaming/test_radiation_streaming.cpp @@ -25,7 +25,6 @@ constexpr double rho = 1.0; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 1.0; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; @@ -38,12 +37,15 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = 1.0; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = 1.0; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; }; diff --git a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp index 2c6894291..94f43e9ad 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -22,13 +22,12 @@ constexpr int direction = 1; constexpr double initial_Erad = 1.0e-5; constexpr double initial_Egas = 1.0e-5; constexpr double c = 1.0; // speed of light -constexpr double chat = c; // reduced speed of light +constexpr double chat = 0.3; // reduced speed of light constexpr double kappa0 = 1.0e-10; // opacity constexpr double rho = 1.0; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 1.0; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; @@ -41,12 +40,15 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = 1.0; + static constexpr double radiation_constant = 1.0; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = 1.0; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; }; diff --git a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp index 7d64fa4db..b03ef0d93 100644 --- a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp +++ b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp @@ -42,11 +42,8 @@ constexpr double Q = (1.0 / (2.0 * x0)); // do NOT change this constexpr double S = Q * (a_rad * (T_hohlraum * T_hohlraum * T_hohlraum * T_hohlraum)); // erg cm^{-3} template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = c; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double mean_molecular_mass = 1.0; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; @@ -54,7 +51,6 @@ template <> struct RadSystem_Traits { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 1.0; - static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; @@ -68,6 +64,11 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = 1.0; + static constexpr double radiation_constant = 1.0; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadTophat/test_radiation_tophat.cpp b/src/problems/RadTophat/test_radiation_tophat.cpp index eff81e77d..b0f8adbf7 100644 --- a/src/problems/RadTophat/test_radiation_tophat.cpp +++ b/src/problems/RadTophat/test_radiation_tophat.cpp @@ -31,19 +31,16 @@ constexpr double T_hohlraum = 500. / kelvin_to_eV; // K [== 500 eV] constexpr double T_initial = 50. / kelvin_to_eV; // K [== 50 eV] constexpr double c_v = (1.0e15 * 1.0e-6 * kelvin_to_eV); // erg g^-1 K^-1 -constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4 -constexpr double c = 2.99792458e10; // cm s^-1 +constexpr double a_rad = C::a_rad; // erg cm^-3 K^-4 +constexpr double c = C::c_light; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; }; @@ -57,6 +54,7 @@ template <> 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; }; template <> AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index 046519773..b35ae1d9f 100644 --- a/src/problems/RadTube/test_radiation_tube.cpp +++ b/src/problems/RadTube/test_radiation_tube.cpp @@ -42,7 +42,6 @@ constexpr double a0 = 4.0295519855200705e7; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = gamma_gas; }; @@ -56,12 +55,11 @@ template <> struct Physics_Traits { static constexpr bool is_mhd_enabled = false; // number of radiation groups static constexpr int nGroups = 2; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = 10.0 * a0; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 10.0 * a0 / C::c_light; static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::k_B; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0.01 * T0, 3.3 * T0, 1000. * T0}; // Kelvin diff --git a/src/problems/RadhydroBB/test_radhydro_bb.cpp b/src/problems/RadhydroBB/test_radhydro_bb.cpp index f7f57f55b..aeb9bfc07 100644 --- a/src/problems/RadhydroBB/test_radhydro_bb.cpp +++ b/src/problems/RadhydroBB/test_radhydro_bb.cpp @@ -99,7 +99,6 @@ constexpr double erad_floor = a_rad * 1e-30; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -112,12 +111,15 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_groups_; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = k_B; + static constexpr double gravitational_constant = 1.0; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; static constexpr double energy_unit = nu_unit; diff --git a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp index 181ed8821..7285620b8 100644 --- a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp +++ b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp @@ -34,26 +34,20 @@ AMREX_GPU_MANAGED double v0_adv = 1.0e6; // NOLINT template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; }; @@ -67,6 +61,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct Physics_Traits { // cell-centred @@ -77,6 +72,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; AMREX_GPU_HOST_DEVICE diff --git a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp index 2415c3e84..13cdb2768 100644 --- a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -21,8 +21,7 @@ constexpr double T1 = 2.0e7; // K (temperature) constexpr double rho0 = 1.2; // g cm^-3 (matter density) constexpr double a_rad = C::a_rad; constexpr double c = C::c_light; // speed of light (cgs) -constexpr double chat = c; -constexpr double width = 24.0; // cm, width of the pulse +constexpr double width = 24.0; // cm, width of the pulse constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-10; constexpr double mu = 2.33 * C::m_u; constexpr double k_B = C::k_B; @@ -35,26 +34,20 @@ AMREX_GPU_MANAGED double v0_adv = 3.0e7; // NOLINT template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; }; @@ -68,6 +61,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct Physics_Traits { // cell-centred @@ -78,6 +72,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; AMREX_GPU_HOST_DEVICE diff --git a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp index 8c6536e0f..bcc366e36 100644 --- a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp +++ b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp @@ -18,8 +18,6 @@ constexpr double T0 = 1.0e7; // K (temperature) constexpr double T1 = 2.0e7; // K (temperature) constexpr double rho0 = 1.2; // g cm^-3 (matter density) constexpr double a_rad = C::a_rad; -constexpr double c = C::c_light; // speed of light (cgs) -constexpr double chat = c; constexpr double width = 24.0; // cm, width of the pulse constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-10; constexpr double mu = 2.33 * C::m_u; @@ -37,26 +35,20 @@ AMREX_GPU_MANAGED double v0_adv = 1.0e6; // NOLINT template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 2; }; @@ -70,6 +62,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct Physics_Traits { // cell-centred @@ -80,6 +73,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; AMREX_GPU_HOST_DEVICE @@ -171,7 +165,6 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(q const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); const double v0 = v0_adv; - // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; state_cc(i, j, k, RadSystem::radEnergy_index) = Erad; state_cc(i, j, k, RadSystem::x1RadFlux_index) = 4. / 3. * v0 * Erad; state_cc(i, j, k, RadSystem::x2RadFlux_index) = 0; diff --git a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp index 93bc13c9a..18ceebfc8 100644 --- a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp +++ b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp @@ -85,7 +85,6 @@ auto compute_exact_rho(const double x) -> double template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -98,12 +97,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; }; @@ -150,7 +148,6 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka: template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -163,12 +160,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_groups_; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr double energy_unit = h_planck; static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; diff --git a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp index befd25c4d..e0d3d354b 100644 --- a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp +++ b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp @@ -92,12 +92,10 @@ constexpr double coeff_ = h_planck * nu_ref / (k_B * T_ref); // = 4.799243073 = template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; @@ -110,6 +108,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = n_groups_; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct Physics_Traits { // cell-centred @@ -120,12 +119,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr bool compute_v_over_c_terms = true; static constexpr double energy_unit = h_planck; @@ -134,9 +132,7 @@ template <> struct RadSystem_Traits { static constexpr OpacityModel opacity_model = opacity_model_; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = erad_floor; static constexpr bool compute_v_over_c_terms = true; static constexpr int beta_order = 1; diff --git a/src/problems/RadhydroShell/test_radhydro_shell.cpp b/src/problems/RadhydroShell/test_radhydro_shell.cpp index 13bfc1e86..5e60a3592 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -36,24 +36,21 @@ struct ShellProblem { // if false, use octant symmetry constexpr bool simulate_full_box = true; -constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4 -constexpr double c = 2.99792458e10; // cm s^-1 -constexpr double a0 = 2.0e5; // ('reference' sound speed) [cm s^-1] -constexpr double chat = 860. * a0; // cm s^-1 -constexpr double k_B = C::k_B; // erg K^-1 -constexpr double m_H = C::m_u; // atomic mass unit +constexpr double a_rad = C::a_rad; +constexpr double c = C::c_light; +constexpr double a0 = 2.0e5; // ('reference' sound speed) [cm s^-1] +constexpr double chat = 860. * a0; // cm s^-1 +constexpr double k_B = C::k_B; // erg K^-1 +constexpr double m_H = C::m_u; // atomic mass unit constexpr double gamma_gas = 5. / 3.; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 2.2 * m_H; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = gamma_gas; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -71,6 +68,7 @@ template <> 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; }; constexpr amrex::Real Msun = 2.0e33; // g diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index dd04e87a9..42e7c88d2 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -43,6 +43,8 @@ constexpr double v1 = (Mach0 * c_s0) * (rho0 / rho1); constexpr double chat = 10.0 * (v0 + c_s0); // reduced speed of light +constexpr double Ggrav = 1.0; // dimensionless gravitational constant; arbitrary + constexpr double Erad0 = a_rad * (T0 * T0 * T0 * T0); constexpr double Egas0 = rho0 * c_v * T0; constexpr double Erad1 = a_rad * (T1 * T1 * T1 * T1); @@ -54,16 +56,13 @@ constexpr double shock_position = 0.0130; // 0.0132; // cm // we initialize slightly to the left...) template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = gamma_gas; }; @@ -76,6 +75,11 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = k_B; + static constexpr double gravitational_constant = Ggrav; + static constexpr double c_light = c; + static constexpr double radiation_constant = a_rad; }; template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 940054e7b..9555be7ae 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -55,16 +55,13 @@ constexpr double shock_position = 0.01305; // 0.0132; // cm (shock position drif constexpr double Lx = 0.01575; // cm template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_p + C::m_e; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = gamma_gas; }; @@ -77,6 +74,12 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + // A custom unit system is used here to replicate the CGS units, for testing units conversion + static constexpr UnitSystem unit_system = UnitSystem::CUSTOM; + static constexpr double unit_length = 1.0; // cm + static constexpr double unit_mass = 1.0; // g + static constexpr double unit_time = 1.0; // s + static constexpr double unit_temperature = 1.0; // K }; template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index 40552c359..641472c88 100644 --- a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -52,12 +52,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 5; + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = Erad_floor_; static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.00000000e+15, 1.00000000e+16, 1.00000000e+17, @@ -70,7 +69,6 @@ template <> struct RadSystem_Traits { template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_p + C::m_e; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = gamma_gas; }; diff --git a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp index c1085a384..8e6375659 100644 --- a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp +++ b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp @@ -32,7 +32,7 @@ constexpr double c = 1.0e8; constexpr int beta_order_ = 2; // order of beta in the radiation four-force constexpr double v0 = 1e-2 * c; constexpr double kappa0 = 1.0e5; -constexpr double chat = 1.0e8; +constexpr double chat_over_c = 1.0; constexpr double T0 = 1.0; // temperature constexpr double rho0 = 1.0; // matter density @@ -52,14 +52,11 @@ constexpr double Erad_beta2 = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad0; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; + static constexpr double c_hat_over_c = chat_over_c; static constexpr double Erad_floor = 0.0; static constexpr int beta_order = beta_order_; }; @@ -73,6 +70,11 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = k_B; + static constexpr double c_light = c; + static constexpr double gravitational_constant = 1.0; + static constexpr double radiation_constant = a_rad; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RandomBlast/blast.cpp b/src/problems/RandomBlast/blast.cpp index bcf749315..81eddd3a0 100644 --- a/src/problems/RandomBlast/blast.cpp +++ b/src/problems/RandomBlast/blast.cpp @@ -28,6 +28,7 @@ #include "fundamental_constants.H" #include "hydro/hydro_system.hpp" #include "math/quadrature.hpp" +#include "physics_info.hpp" using amrex::Real; @@ -47,12 +48,12 @@ template <> struct Physics_Traits { static constexpr int numMassScalars = 0; static constexpr int numPassiveScalars = numMassScalars + 1; static constexpr int nGroups = 1; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; }; constexpr Real Tgas0 = 1.0e4; // K @@ -287,6 +288,9 @@ template <> void QuokkaSimulation::ErrorEst(int lev, amrex::TagBoxA auto problem_main() -> int { + // This problem is only implemented in CGS units because the Grackle cooling tables are in CGS units. + static_assert(Physics_Traits::unit_system == UnitSystem::CGS); + // read parameters amrex::ParmParse const pp; diff --git a/src/problems/RayleighTaylor2D/test_hydro2d_rt.cpp b/src/problems/RayleighTaylor2D/test_hydro2d_rt.cpp index 231742451..164d4c875 100644 --- a/src/problems/RayleighTaylor2D/test_hydro2d_rt.cpp +++ b/src/problems/RayleighTaylor2D/test_hydro2d_rt.cpp @@ -21,7 +21,6 @@ struct RTProblem { template <> struct quokka::EOS_Traits { 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 { @@ -37,6 +36,7 @@ template <> 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; }; amrex::Real constexpr g_x = 0; diff --git a/src/problems/RayleighTaylor3D/test_hydro3d_rt.cpp b/src/problems/RayleighTaylor3D/test_hydro3d_rt.cpp index f4ab35faa..60c2a9962 100644 --- a/src/problems/RayleighTaylor3D/test_hydro3d_rt.cpp +++ b/src/problems/RayleighTaylor3D/test_hydro3d_rt.cpp @@ -22,7 +22,6 @@ struct RTProblem { template <> struct quokka::EOS_Traits { 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 { @@ -38,6 +37,7 @@ template <> 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; }; amrex::Real constexpr g_x = 0; diff --git a/src/problems/ShockCloud/cloud.cpp b/src/problems/ShockCloud/cloud.cpp index 5a2bf7f69..aec244207 100644 --- a/src/problems/ShockCloud/cloud.cpp +++ b/src/problems/ShockCloud/cloud.cpp @@ -50,12 +50,12 @@ template <> struct Physics_Traits { static constexpr int numMassScalars = 0; static constexpr int numPassiveScalars = numMassScalars + 3; static constexpr int nGroups = 1; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; }; template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; }; /// global variables diff --git a/src/problems/SphericalCollapse/spherical_collapse.cpp b/src/problems/SphericalCollapse/spherical_collapse.cpp index 5e87b9ef6..064ffe735 100644 --- a/src/problems/SphericalCollapse/spherical_collapse.cpp +++ b/src/problems/SphericalCollapse/spherical_collapse.cpp @@ -29,7 +29,6 @@ struct CollapseProblem { template <> struct quokka::EOS_Traits { 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 { @@ -45,6 +44,9 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = C::k_B; + static constexpr double gravitational_constant = 1.0; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) diff --git a/src/problems/StarCluster/star_cluster.cpp b/src/problems/StarCluster/star_cluster.cpp index 84127f519..f76830935 100644 --- a/src/problems/StarCluster/star_cluster.cpp +++ b/src/problems/StarCluster/star_cluster.cpp @@ -38,7 +38,6 @@ template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.0; static constexpr double cs_isothermal = 1.0; // dimensionless static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; }; template <> struct HydroSystem_Traits { @@ -54,6 +53,9 @@ template <> 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::CONSTANTS; + static constexpr double boltzmann_constant = C::k_B; + static constexpr amrex::Real gravitational_constant = 1.0; }; template <> struct SimulationData { diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 6c10fa13e..cdf778919 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -22,6 +22,7 @@ #include "AMReX_REAL.H" // internal headers +#include "fundamental_constants.H" #include "hydro/EOS.hpp" #include "hyperbolic_system.hpp" #include "math/math_impl.hpp" @@ -71,9 +72,7 @@ enum class OpacityModel { // this struct is specialized by the user application code // template struct RadSystem_Traits { - static constexpr double c_light = c_light_cgs_; - static constexpr double c_hat = c_light_cgs_; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double c_hat_over_c = 1.0; static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries = {0., inf}; @@ -189,9 +188,32 @@ template class RadSystem : public HyperbolicSystem::c_light; - static constexpr double c_hat_ = RadSystem_Traits::c_hat; - static constexpr double radiation_constant_ = RadSystem_Traits::radiation_constant; + + static constexpr amrex::Real c_light_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return c_light_cgs_; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return Physics_Traits::c_light; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + // c / c_bar = u_l / u_t + return c_light_cgs_ / (Physics_Traits::unit_length / Physics_Traits::unit_time); + } + }(); + static constexpr double c_hat_ = c_light_ * RadSystem_Traits::c_hat_over_c; + + static constexpr double radiation_constant_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return C::a_rad; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return Physics_Traits::radiation_constant; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + // a_rad / a_rad_bar = 1 / u_l * u_m / u_t^2 / u_T^4 + return C::a_rad / (1.0 / Physics_Traits::unit_length * Physics_Traits::unit_mass / + (Physics_Traits::unit_time * Physics_Traits::unit_time) / + (Physics_Traits::unit_temperature * Physics_Traits::unit_temperature * + Physics_Traits::unit_temperature * Physics_Traits::unit_temperature)); + } + }(); static constexpr int beta_order_ = RadSystem_Traits::beta_order; @@ -227,9 +249,21 @@ template class RadSystem : public HyperbolicSystem::mean_molecular_weight; - static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; static constexpr double gamma_ = quokka::EOS_Traits::gamma; + static constexpr amrex::Real boltzmann_constant_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return C::k_B; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return Physics_Traits::boltzmann_constant; + } else if constexpr (Physics_Traits::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::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_mass / + (Physics_Traits::unit_time * Physics_Traits::unit_time) / Physics_Traits::unit_temperature); + } + }(); + // static functions static void ComputeMaxSignalSpeed(amrex::Array4 const &cons, array_t &maxSignal, amrex::Box const &indexRange); diff --git a/src/simulation.hpp b/src/simulation.hpp index 3b5876cf6..08d64debf 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -335,6 +335,9 @@ template class AMRSimulation : public amrex::AmrCore void kickParticlesAllLevels(amrex::Real dt); void driftParticlesAllLevels(amrex::Real dt); + // simulation metadata + void initializeSimulationMetadata(); + #ifdef AMREX_USE_ASCENT void AscentCustomActions(conduit::Node const &blueprintMesh); void RenderAscent(); @@ -389,7 +392,56 @@ template class AMRSimulation : public amrex::AmrCore amrex::Vector cellUpdatesEachLevel_; // gravity - amrex::Real Gconst_ = C::Gconst; // gravitational constant G + static constexpr amrex::Real Gconst_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return C::Gconst; // gravitational constant G, CGS units + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return Physics_Traits::gravitational_constant; // gravitational constant G, user defined + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + // G / G_bar = u_l^3 / u_m / u_t^2 + return C::Gconst / + (Physics_Traits::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_length / + Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time)); + } + }(); + + // unit length, mass, time, temperature + static constexpr amrex::Real unit_length = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_length; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return 1.0; + } else { + return NAN; + } + }(); + static constexpr amrex::Real unit_mass = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_mass; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return 1.0; + } else { + return NAN; + } + }(); + static constexpr amrex::Real unit_time = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_time; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return 1.0; + } else { + return NAN; + } + }(); + static constexpr amrex::Real unit_temperature = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_temperature; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return 1.0; + } else { + return NAN; + } + }(); // tracer particles #ifdef AMREX_PARTICLES @@ -475,6 +527,10 @@ template void AMRSimulation::initialize() } } + if constexpr (Physics_Traits::is_hydro_enabled || Physics_Traits::is_radiation_enabled) { + initializeSimulationMetadata(); + } + #ifdef AMREX_USE_ASCENT // initialize Ascent conduit::Node ascent_options; @@ -627,12 +683,6 @@ template void AMRSimulation::readParameters() maxWalltime_ = 3600 * hours + 60 * minutes + seconds; amrex::Print() << fmt::format("Setting walltime limit to {} hours, {} minutes, {} seconds.\n", hours, minutes, seconds); } - - // set gravity runtime parameters - { - const amrex::ParmParse hpp("gravity"); - hpp.query("Gconst", Gconst_); - } } template void AMRSimulation::setInitialConditions() diff --git a/tests/RadStreamingY.in b/tests/RadStreamingY.in index c30fcd360..5fcb9407b 100644 --- a/tests/RadStreamingY.in +++ b/tests/RadStreamingY.in @@ -1,4 +1,4 @@ -max_time = 0.2 +max_time = 1.0 # ***************************************************************** # Problem size and geometry diff --git a/tests/RadTube.in b/tests/RadTube.in index 82ecf3c01..3c6084ebe 100644 --- a/tests/RadTube.in +++ b/tests/RadTube.in @@ -19,4 +19,6 @@ amr.blocking_factor = 4 # grid size must be divisible by this do_reflux = 0 do_subcycle = 0 -suppress_output = 1 \ No newline at end of file +suppress_output = 1 + +plotfile_interval = 10000000 \ No newline at end of file diff --git a/tests/SphericalCollapse.in b/tests/SphericalCollapse.in index cb5c328f9..cb51ee1d8 100644 --- a/tests/SphericalCollapse.in +++ b/tests/SphericalCollapse.in @@ -27,8 +27,6 @@ stop_time = 0.05 derived_vars = gpot -gravity.Gconst = 1.0 # gravitational constant - do_reflux = 1 do_subcycle = 0 do_tracers = 0 # turn on tracer particles diff --git a/tests/radshock.in b/tests/radshock.in index 36dc68345..3fba25884 100644 --- a/tests/radshock.in +++ b/tests/radshock.in @@ -13,9 +13,11 @@ amr.v = 0 # verbosity in Amr # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 512 4 4 +amr.n_cell = 128 4 4 amr.max_level = 0 # number of levels = max_level + 1 amr.blocking_factor = 4 # grid size must be divisible by this do_reflux = 0 do_subcycle = 0 + +plotfile_interval = 10000000 \ No newline at end of file diff --git a/tests/radshock_dimensionless.in b/tests/radshock_dimensionless.in index 063f20bf1..c41c593f2 100644 --- a/tests/radshock_dimensionless.in +++ b/tests/radshock_dimensionless.in @@ -13,12 +13,13 @@ amr.v = 0 # verbosity in Amr # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 256 4 4 +amr.n_cell = 128 4 4 amr.max_level = 0 # number of levels = max_level + 1 amr.blocking_factor = 4 # grid size must be divisible by this amr.max_grid_size_x = 256 do_reflux = 0 do_subcycle = 0 +plotfile_interval = 10000000 radiation.print_iteration_counts = 1 \ No newline at end of file