From fce46d7d3d00fd322bfe18fda9daa8d32214f297 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 26 Oct 2024 17:50:27 +1100 Subject: [PATCH 01/34] reduce nx in shock tests --- tests/radshock.in | 2 +- tests/radshock_dimensionless.in | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/radshock.in b/tests/radshock.in index 36dc68345..0826fa41b 100644 --- a/tests/radshock.in +++ b/tests/radshock.in @@ -13,7 +13,7 @@ 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 diff --git a/tests/radshock_dimensionless.in b/tests/radshock_dimensionless.in index 063f20bf1..01c23ddaa 100644 --- a/tests/radshock_dimensionless.in +++ b/tests/radshock_dimensionless.in @@ -13,7 +13,7 @@ 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 From b7276fcdc612c9cab7796935e351d819ec282530 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 26 Oct 2024 17:51:38 +1100 Subject: [PATCH 02/34] define cgs and k_B: passing shuosher and shock_cgs tests --- src/hydro/EOS.hpp | 9 ++++++++- src/physics_info.hpp | 14 ++++++++++++++ src/problems/HydroShuOsher/test_hydro_shuosher.cpp | 2 +- .../RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 2 +- 4 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index a8c5a5e60..cc948ea37 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -67,8 +67,15 @@ template class EOS 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 { + return Physics_Traits::boltzmann_constant; + } + }(); }; template diff --git a/src/physics_info.hpp b/src/physics_info.hpp index e41ff1fda..d3ddbf189 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -2,8 +2,16 @@ #define PHYSICS_INFO_HPP_ #include "physics_numVars.hpp" +#include "fundamental_constants.H" #include +// enum for unit system +enum class UnitSystem { + CGS = 0, + CONSTANTS = 1, + CUSTOM = 2 +}; + // this struct is specialized by the user application code. template struct Physics_Traits { // cell-centred @@ -14,6 +22,12 @@ 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 c_hat = C::c_light; // Radiation + static constexpr double radiation_constant = C::a_rad; // Radiation }; // this struct stores the indices at which quantities start 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/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 940054e7b..b6738c08a 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -64,7 +64,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; }; @@ -77,6 +76,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 <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real From 6fedb5582cc22007f9e7e73e5265daad9799d09b Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 26 Oct 2024 23:28:15 +1100 Subject: [PATCH 03/34] cgs and k_B: shock test --- src/hydro/EOS.hpp | 3 ++- src/physics_info.hpp | 5 ++--- src/problems/RadhydroShock/test_radhydro_shock.cpp | 1 - 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index cc948ea37..544de17a3 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -69,10 +69,11 @@ template class EOS static constexpr amrex::Real gamma_ = EOS_Traits::gamma; 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 { + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { return Physics_Traits::boltzmann_constant; } }(); diff --git a/src/physics_info.hpp b/src/physics_info.hpp index d3ddbf189..af53f4d84 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -25,9 +25,8 @@ template struct Physics_Traits { 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 c_hat = C::c_light; // Radiation - static constexpr double radiation_constant = C::a_rad; // Radiation + // static constexpr double c_light = C::c_light; // Radiation + // static constexpr double radiation_constant = C::a_rad; // Radiation }; // this struct stores the indices at which quantities start diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index dd04e87a9..0d8d5d175 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -63,7 +63,6 @@ template <> struct RadSystem_Traits { 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; }; From 9a9dbead8dadcfefc300bc339c630519344cdb16 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 26 Oct 2024 23:34:27 +1100 Subject: [PATCH 04/34] finish CGS and CONSTANTS definition, passing shock and shockcgs --- .../RadhydroShock/test_radhydro_shock.cpp | 2 ++ .../test_radhydro_shock_cgs.cpp | 2 -- src/radiation/radiation_system.hpp | 19 +++++++++++++++++-- 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index 0d8d5d175..adba7387a 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -75,6 +75,8 @@ 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; }; 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 b6738c08a..bf6b8ce76 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -55,9 +55,7 @@ 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 Erad_floor = 0.; static constexpr int beta_order = 1; }; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 6c10fa13e..c0687a33e 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" @@ -189,9 +190,23 @@ template class RadSystem : public HyperbolicSystem::c_light; + + static constexpr amrex::Real c_light_ = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { + return C::c_light; + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + return RadSystem_Traits::c_light; + } + }(); static constexpr double c_hat_ = RadSystem_Traits::c_hat; - static constexpr double radiation_constant_ = RadSystem_Traits::radiation_constant; + + 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 RadSystem_Traits::radiation_constant; + } + }(); static constexpr int beta_order_ = RadSystem_Traits::beta_order; From acf7f06d3722874b554150b4d5ce4e388e66df31 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 00:05:23 +1100 Subject: [PATCH 05/34] add G into CGS and CONSTANTS; passing star_cluster test --- src/problems/StarCluster/star_cluster.cpp | 2 ++ src/simulation.hpp | 14 +++++++------- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/problems/StarCluster/star_cluster.cpp b/src/problems/StarCluster/star_cluster.cpp index 84127f519..869c6c5c7 100644 --- a/src/problems/StarCluster/star_cluster.cpp +++ b/src/problems/StarCluster/star_cluster.cpp @@ -54,6 +54,8 @@ 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 amrex::Real gravitational_constant = 1.0; }; template <> struct SimulationData { diff --git a/src/simulation.hpp b/src/simulation.hpp index 3b5876cf6..778bdb001 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -389,7 +389,13 @@ 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 + } + }(); // tracer particles #ifdef AMREX_PARTICLES @@ -627,12 +633,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() From 432ef8a235e3d72cb7cea9ad1e9f6ed501037eec Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 01:03:37 +1100 Subject: [PATCH 06/34] compute constants from unit_length ... when using UnitSystem::CUSTUM --- src/hydro/EOS.hpp | 6 +++++- src/problems/RadhydroShock/test_radhydro_shock.cpp | 3 +++ src/radiation/radiation_system.hpp | 12 +++++++++++- src/simulation.hpp | 5 +++++ 4 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index 544de17a3..753661637 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -69,12 +69,16 @@ template class EOS static constexpr amrex::Real gamma_ = EOS_Traits::gamma; 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 Physics_Traits::boltzmann_constant / (Physics_Traits::unit_length * Physics_Traits::unit_length / Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time) * Physics_Traits::unit_temperature); + } else { + static_assert(false, "Invalid unit system"); } }(); }; diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index adba7387a..b1b001311 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); @@ -77,6 +79,7 @@ template <> struct Physics_Traits { 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; }; template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index c0687a33e..f374014b2 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -193,9 +193,14 @@ template class RadSystem : public HyperbolicSystem::unit_system == UnitSystem::CGS) { - return C::c_light; + return c_light_cgs_; } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { return RadSystem_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); + } else { + static_assert(false, "Invalid unit system"); } }(); static constexpr double c_hat_ = RadSystem_Traits::c_hat; @@ -205,6 +210,11 @@ template class RadSystem : public HyperbolicSystem::unit_system == UnitSystem::CONSTANTS) { return RadSystem_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)); + } else { + static_assert(false, "Invalid unit system"); } }(); diff --git a/src/simulation.hpp b/src/simulation.hpp index 778bdb001..bb349db3a 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -394,6 +394,11 @@ template class AMRSimulation : public amrex::AmrCore 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 Physics_Traits::gravitational_constant / (Physics_Traits::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_length / Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time)); + } else { + static_assert(false, "Invalid unit system"); } }(); From 2688902b25f587c32041f512d3a9b354c5e1eca7 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 01:19:08 +1100 Subject: [PATCH 07/34] fix: use C::Gconst --- src/hydro/EOS.hpp | 2 +- src/simulation.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index 753661637..e450d360f 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -76,7 +76,7 @@ template class EOS 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 Physics_Traits::boltzmann_constant / (Physics_Traits::unit_length * Physics_Traits::unit_length / Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time) * Physics_Traits::unit_temperature); + 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); } else { static_assert(false, "Invalid unit system"); } diff --git a/src/simulation.hpp b/src/simulation.hpp index bb349db3a..9d5c4ce10 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -396,7 +396,7 @@ template class AMRSimulation : public amrex::AmrCore 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 Physics_Traits::gravitational_constant / (Physics_Traits::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_length / Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time)); + 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)); } else { static_assert(false, "Invalid unit system"); } From 92a7e5f44ff1c262b5c638724c3caf228c8de270 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 10:37:38 +1100 Subject: [PATCH 08/34] mod shock_cgs and line_cooling problems to test units conversion --- src/hydro/EOS.hpp | 4 ++-- .../RadLineCooling/test_rad_line_cooling.cpp | 16 +++++++++++++--- .../RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 7 ++++++- src/radiation/radiation_system.hpp | 2 +- 4 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index e450d360f..a4129a4d4 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -75,8 +75,8 @@ template class EOS } 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); + // 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); } else { static_assert(false, "Invalid unit system"); } diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index c8f78e675..f882a148e 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,26 @@ 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 dimentionless units (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; }; 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; + // Equivalently, set + // static constexpr double c_light = c; + // static constexpr double radiation_constant = a_rad; }; template <> struct ISM_Traits { diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index bf6b8ce76..24f0504d8 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -74,7 +74,12 @@ 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; + // 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/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index f374014b2..dab8bb97a 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -252,7 +252,7 @@ template class RadSystem : public HyperbolicSystem::mean_molecular_weight; - static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; + static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant_; static constexpr double gamma_ = quokka::EOS_Traits::gamma; // static functions From e28f2436fced47295ae4465ab865e4477d522002 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 11:04:40 +1100 Subject: [PATCH 09/34] add unit_length to Physics_Traits --- src/physics_info.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/physics_info.hpp b/src/physics_info.hpp index af53f4d84..9538bbcfa 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -25,8 +25,10 @@ template struct Physics_Traits { 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 From 6f7509c32d1a1e6700a387beddac8a9117d8a134 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 15:54:12 +1100 Subject: [PATCH 10/34] finishing initializeSimulationMetadata and writing units to metadata --- src/physics_info.hpp | 6 +-- .../RadLineCooling/test_rad_line_cooling.cpp | 2 +- src/problems/RadTube/test_radiation_tube.cpp | 4 +- src/radiation/radiation_system.hpp | 14 ++++- src/simulation.hpp | 52 +++++++++++++++++++ tests/RadLineCooling.in | 4 +- tests/RadTube.in | 4 +- 7 files changed, 74 insertions(+), 12 deletions(-) diff --git a/src/physics_info.hpp b/src/physics_info.hpp index 9538bbcfa..8ae474cd9 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -6,11 +6,7 @@ #include // enum for unit system -enum class UnitSystem { - CGS = 0, - CONSTANTS = 1, - CUSTOM = 2 -}; +enum class UnitSystem { CGS, CONSTANTS, CUSTOM }; // this struct is specialized by the user application code. template struct Physics_Traits { diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index f882a148e..1a8a598d6 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -56,7 +56,7 @@ 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 dimentionless units (c = k_B = a_rad = G = 1), for testing units conversion + // 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; diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index 046519773..c64e2499c 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 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/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index dab8bb97a..567ba39cc 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -252,9 +252,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); + } else { + static_assert(false, "Invalid unit system"); + } + }(); + // 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 9d5c4ce10..de00084f1 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(); @@ -402,6 +405,36 @@ template class AMRSimulation : public amrex::AmrCore } }(); + // 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 { + return 1.0; + } + }(); + static constexpr amrex::Real unit_mass = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_mass; + } else { + return 1.0; + } + }(); + static constexpr amrex::Real unit_time = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_time; + } else { + return 1.0; + } + }(); + static constexpr amrex::Real unit_temperature = []() constexpr { + if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { + return Physics_Traits::unit_temperature; + } else { + return 1.0; + } + }(); + // tracer particles #ifdef AMREX_PARTICLES void InitParticles(); // create tracer particles @@ -430,6 +463,23 @@ template auto AMRSimulation::getNewMF_fc() const return state_new_fc_; } +// initialize metadata +template void AMRSimulation::initializeSimulationMetadata() +{ + if constexpr (Physics_Traits::unit_system != UnitSystem::CONSTANTS) { + simulationMetadata_["unit_length"] = unit_length; + simulationMetadata_["unit_mass"] = unit_mass; + simulationMetadata_["unit_time"] = unit_time; + simulationMetadata_["unit_temperature"] = unit_temperature; + } else { + // 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"] = "undefined"; + simulationMetadata_["unit_mass"] = "undefined"; + simulationMetadata_["unit_time"] = "undefined"; + simulationMetadata_["unit_temperature"] = "undefined"; + } +} + template void AMRSimulation::initialize() { BL_PROFILE("AMRSimulation::initialize()"); @@ -486,6 +536,8 @@ template void AMRSimulation::initialize() } } + initializeSimulationMetadata(); + #ifdef AMREX_USE_ASCENT // initialize Ascent conduit::Node ascent_options; diff --git a/tests/RadLineCooling.in b/tests/RadLineCooling.in index e0f16b7ac..6272cda90 100644 --- a/tests/RadLineCooling.in +++ b/tests/RadLineCooling.in @@ -23,4 +23,6 @@ radiation.dust_gas_interaction_coeff = 1e-20 do_reflux = 0 do_subcycle = 0 -suppress_output = 1 \ No newline at end of file +suppress_output = 1 + +checkpoint_interval = 10000000 \ No newline at end of file diff --git a/tests/RadTube.in b/tests/RadTube.in index 82ecf3c01..9fcdf3c71 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 + +checkpoint_interval = 10000000 \ No newline at end of file From 5bd9e13190f318906dcc8a80d411cff176923dff Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 16:59:05 +1100 Subject: [PATCH 11/34] write both unit_xxx and constants into metadata for all unit systems --- src/QuokkaSimulation.hpp | 44 +++++++++++++++++++++++++++++++++ src/simulation.hpp | 17 ------------- tests/radshock_dimensionless.in | 1 + 3 files changed, 45 insertions(+), 17 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 42c7cc462..8f231c345 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -321,6 +321,50 @@ 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. + // units + simulationMetadata_["unit_length"] = "undefined"; + simulationMetadata_["unit_mass"] = "undefined"; + simulationMetadata_["unit_time"] = "undefined"; + simulationMetadata_["unit_temperature"] = "undefined"; + + // constants + simulationMetadata_["k_B"] = Physics_Traits::boltzmann_constant; + simulationMetadata_["G"] = Physics_Traits::gravitational_constant; + if constexpr (Physics_Traits::is_radiation_enabled) { + simulationMetadata_["c"] = RadSystem_Traits::c_light; + simulationMetadata_["c_hat"] = RadSystem_Traits::c_hat; + simulationMetadata_["a_rad"] = RadSystem_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/simulation.hpp b/src/simulation.hpp index de00084f1..3a091484a 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -463,23 +463,6 @@ template auto AMRSimulation::getNewMF_fc() const return state_new_fc_; } -// initialize metadata -template void AMRSimulation::initializeSimulationMetadata() -{ - if constexpr (Physics_Traits::unit_system != UnitSystem::CONSTANTS) { - simulationMetadata_["unit_length"] = unit_length; - simulationMetadata_["unit_mass"] = unit_mass; - simulationMetadata_["unit_time"] = unit_time; - simulationMetadata_["unit_temperature"] = unit_temperature; - } else { - // 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"] = "undefined"; - simulationMetadata_["unit_mass"] = "undefined"; - simulationMetadata_["unit_time"] = "undefined"; - simulationMetadata_["unit_temperature"] = "undefined"; - } -} - template void AMRSimulation::initialize() { BL_PROFILE("AMRSimulation::initialize()"); diff --git a/tests/radshock_dimensionless.in b/tests/radshock_dimensionless.in index 01c23ddaa..c7295b289 100644 --- a/tests/radshock_dimensionless.in +++ b/tests/radshock_dimensionless.in @@ -20,5 +20,6 @@ amr.max_grid_size_x = 256 do_reflux = 0 do_subcycle = 0 +checkpoint_interval = 10000000 radiation.print_iteration_counts = 1 \ No newline at end of file From 621fe4ecb0b554aec45b14c80fee0891fa3dd419 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 17:01:42 +1100 Subject: [PATCH 12/34] update in files --- tests/RadLineCooling.in | 2 +- tests/RadTube.in | 2 +- tests/radshock.in | 2 ++ tests/radshock_dimensionless.in | 2 +- 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/RadLineCooling.in b/tests/RadLineCooling.in index 6272cda90..b77624b49 100644 --- a/tests/RadLineCooling.in +++ b/tests/RadLineCooling.in @@ -25,4 +25,4 @@ do_reflux = 0 do_subcycle = 0 suppress_output = 1 -checkpoint_interval = 10000000 \ No newline at end of file +plotfile_interval = 10000000 \ No newline at end of file diff --git a/tests/RadTube.in b/tests/RadTube.in index 9fcdf3c71..3c6084ebe 100644 --- a/tests/RadTube.in +++ b/tests/RadTube.in @@ -21,4 +21,4 @@ do_reflux = 0 do_subcycle = 0 suppress_output = 1 -checkpoint_interval = 10000000 \ No newline at end of file +plotfile_interval = 10000000 \ No newline at end of file diff --git a/tests/radshock.in b/tests/radshock.in index 0826fa41b..3fba25884 100644 --- a/tests/radshock.in +++ b/tests/radshock.in @@ -19,3 +19,5 @@ 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 c7295b289..c41c593f2 100644 --- a/tests/radshock_dimensionless.in +++ b/tests/radshock_dimensionless.in @@ -20,6 +20,6 @@ amr.max_grid_size_x = 256 do_reflux = 0 do_subcycle = 0 -checkpoint_interval = 10000000 +plotfile_interval = 10000000 radiation.print_iteration_counts = 1 \ No newline at end of file From babfcefe321c3d4b7cdc3f0cb81f751ba4492626 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 17:50:11 +1100 Subject: [PATCH 13/34] add comments --- src/physics_info.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics_info.hpp b/src/physics_info.hpp index 8ae474cd9..c35964f96 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -5,7 +5,7 @@ #include "fundamental_constants.H" #include -// enum for unit system +// enum for unit system, one of CGS, CONSTANTS, CUSTOM enum class UnitSystem { CGS, CONSTANTS, CUSTOM }; // this struct is specialized by the user application code. From 7316bb8198d4d9759e5b74ca5d566adc27c30703 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sun, 27 Oct 2024 17:59:36 +1100 Subject: [PATCH 14/34] set unit_xxx to NAN when units_system == CONSTANTS --- src/simulation.hpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/simulation.hpp b/src/simulation.hpp index 3a091484a..2f708b514 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -409,29 +409,37 @@ template class AMRSimulation : public amrex::AmrCore static constexpr amrex::Real unit_length = []() constexpr { if constexpr (Physics_Traits::unit_system == UnitSystem::CUSTOM) { return Physics_Traits::unit_length; - } else { + } 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 { + } 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 { + } 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 { + } else if constexpr (Physics_Traits::unit_system == UnitSystem::CGS) { return 1.0; + } else { + return NAN; } }(); From 25d6c4457e5fb91a6c19a1b052dd50af946214b6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 27 Oct 2024 07:11:22 +0000 Subject: [PATCH 15/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/QuokkaSimulation.hpp | 9 ++++++--- src/hydro/EOS.hpp | 6 ++++-- src/physics_info.hpp | 4 ++-- .../RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 6 +++--- src/radiation/radiation_system.hpp | 11 ++++++++--- src/simulation.hpp | 4 +++- 6 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 8f231c345..3e8018e2e 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -325,8 +325,9 @@ template void QuokkaSimulation::defineComponentN 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. - // units + // 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. units simulationMetadata_["unit_length"] = "undefined"; simulationMetadata_["unit_mass"] = "undefined"; simulationMetadata_["unit_time"] = "undefined"; @@ -353,7 +354,9 @@ template void AMRSimulation::initializeSimulatio 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); + 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_; diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index a4129a4d4..b629a2b22 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -72,11 +72,13 @@ template class EOS 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) { + } 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); + 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); } else { static_assert(false, "Invalid unit system"); } diff --git a/src/physics_info.hpp b/src/physics_info.hpp index c35964f96..2f880fdb0 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -1,8 +1,8 @@ #ifndef PHYSICS_INFO_HPP_ // NOLINT #define PHYSICS_INFO_HPP_ -#include "physics_numVars.hpp" #include "fundamental_constants.H" +#include "physics_numVars.hpp" #include // enum for unit system, one of CGS, CONSTANTS, CUSTOM @@ -19,7 +19,7 @@ template struct Physics_Traits { 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 boltzmann_constant = C::k_B; // Hydro, EOS static constexpr double gravitational_constant = C::Gconst; // gravity static constexpr double unit_length = 1.0; static constexpr double unit_mass = 1.0; diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 24f0504d8..763b7263b 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -76,9 +76,9 @@ template <> struct Physics_Traits { 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_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 }; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 567ba39cc..56fa3f5fd 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -212,7 +212,10 @@ template class RadSystem : public HyperbolicSystem::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)); + 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)); } else { static_assert(false, "Invalid unit system"); } @@ -257,11 +260,13 @@ template class RadSystem : public HyperbolicSystem::unit_system == UnitSystem::CGS) { return C::k_B; - } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { + } 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); + 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); } else { static_assert(false, "Invalid unit system"); } diff --git a/src/simulation.hpp b/src/simulation.hpp index 2f708b514..b39e8ea04 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -399,7 +399,9 @@ template class AMRSimulation : public amrex::AmrCore 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)); + 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)); } else { static_assert(false, "Invalid unit system"); } From 2010aa7af14721ac7764f8e569d8787e845790be Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 11:21:03 +1100 Subject: [PATCH 16/34] set unit_mass to NAN if CONSTNATS --- src/QuokkaSimulation.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 3e8018e2e..1ceb134ec 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -327,11 +327,11 @@ template void AMRSimulation::initializeSimulatio 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. units - simulationMetadata_["unit_length"] = "undefined"; - simulationMetadata_["unit_mass"] = "undefined"; - simulationMetadata_["unit_time"] = "undefined"; - simulationMetadata_["unit_temperature"] = "undefined"; + // 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; From 1c2e29868cb21fac96e5014eef6b94303647cd26 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 11:35:10 +1100 Subject: [PATCH 17/34] move c_light and c_hat to PhysicsTraits --- src/QuokkaSimulation.hpp | 6 +++--- src/physics_info.hpp | 3 +++ src/problems/RadLineCooling/test_rad_line_cooling.cpp | 2 +- src/problems/RadTube/test_radiation_tube.cpp | 2 +- src/problems/RadhydroShock/test_radhydro_shock.cpp | 6 +++--- src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 2 +- src/radiation/radiation_system.hpp | 8 +++----- 7 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 1ceb134ec..cab6821d5 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -337,9 +337,9 @@ template void AMRSimulation::initializeSimulatio simulationMetadata_["k_B"] = Physics_Traits::boltzmann_constant; simulationMetadata_["G"] = Physics_Traits::gravitational_constant; if constexpr (Physics_Traits::is_radiation_enabled) { - simulationMetadata_["c"] = RadSystem_Traits::c_light; - simulationMetadata_["c_hat"] = RadSystem_Traits::c_hat; - simulationMetadata_["a_rad"] = RadSystem_Traits::radiation_constant; + simulationMetadata_["c"] = Physics_Traits::c_light; + simulationMetadata_["c_hat"] = Physics_Traits::c_hat; + simulationMetadata_["a_rad"] = Physics_Traits::radiation_constant; } } else { // units diff --git a/src/physics_info.hpp b/src/physics_info.hpp index 2f880fdb0..bfa408957 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -21,6 +21,9 @@ template struct Physics_Traits { 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 c_hat = 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; diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index 1a8a598d6..1234c26fd 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -58,6 +58,7 @@ template <> struct Physics_Traits { 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 c_hat = chat; static constexpr double unit_length = 1.733039549e-33; static constexpr double unit_mass = 2.333695323e-05; static constexpr double unit_time = 5.780797690e-44; @@ -69,7 +70,6 @@ template <> struct Physics_Traits { }; template <> struct RadSystem_Traits { - static constexpr double c_hat = chat; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 0; static constexpr double energy_unit = nu_unit; diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index c64e2499c..de2087700 100644 --- a/src/problems/RadTube/test_radiation_tube.cpp +++ b/src/problems/RadTube/test_radiation_tube.cpp @@ -56,10 +56,10 @@ template <> struct Physics_Traits { // number of radiation groups static constexpr int nGroups = 2; static constexpr UnitSystem unit_system = UnitSystem::CGS; + static constexpr double c_hat = 10.0 * a0; }; template <> struct RadSystem_Traits { - static constexpr double c_hat = 10.0 * a0; 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/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index b1b001311..d0eab6af5 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -56,9 +56,6 @@ 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 Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -80,6 +77,9 @@ template <> struct Physics_Traits { 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 c_hat = chat; + 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 763b7263b..59578ab50 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -55,7 +55,6 @@ 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_hat = chat; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -76,6 +75,7 @@ template <> struct Physics_Traits { 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 c_hat = chat; static constexpr double unit_length = 1.0; // cm static constexpr double unit_mass = 1.0; // g static constexpr double unit_time = 1.0; // s diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 56fa3f5fd..b98b7250a 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -72,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 Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries = {0., inf}; @@ -195,7 +193,7 @@ template class RadSystem : public HyperbolicSystem::unit_system == UnitSystem::CGS) { return c_light_cgs_; } else if constexpr (Physics_Traits::unit_system == UnitSystem::CONSTANTS) { - return RadSystem_Traits::c_light; + 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); @@ -203,13 +201,13 @@ template class RadSystem : public HyperbolicSystem::c_hat; + static constexpr double c_hat_ = Physics_Traits::c_hat; 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 RadSystem_Traits::radiation_constant; + 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 / From 55e4dd61cc7cb158157f96059c530de92df7ce1c Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 11:42:10 +1100 Subject: [PATCH 18/34] change c_hat to c_hat_over_c --- src/QuokkaSimulation.hpp | 2 +- src/physics_info.hpp | 1 - src/problems/RadLineCooling/test_rad_line_cooling.cpp | 2 +- src/problems/RadTube/test_radiation_tube.cpp | 2 +- src/problems/RadhydroShock/test_radhydro_shock.cpp | 2 +- src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 2 +- src/radiation/radiation_system.hpp | 4 ++-- 7 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index cab6821d5..f6fff3738 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -338,7 +338,7 @@ template void AMRSimulation::initializeSimulatio 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_hat; + simulationMetadata_["c_hat"] = Physics_Traits::c_light * RadSystem_Traits::c_hat_over_c; simulationMetadata_["a_rad"] = Physics_Traits::radiation_constant; } } else { diff --git a/src/physics_info.hpp b/src/physics_info.hpp index bfa408957..8618b9136 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -22,7 +22,6 @@ template struct Physics_Traits { 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 c_hat = 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; diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index 1234c26fd..0d7de1b1f 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -58,7 +58,6 @@ template <> struct Physics_Traits { 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 c_hat = chat; static constexpr double unit_length = 1.733039549e-33; static constexpr double unit_mass = 2.333695323e-05; static constexpr double unit_time = 5.780797690e-44; @@ -73,6 +72,7 @@ template <> struct RadSystem_Traits { 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; // Equivalently, set // static constexpr double c_light = c; // static constexpr double radiation_constant = a_rad; diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index de2087700..b35ae1d9f 100644 --- a/src/problems/RadTube/test_radiation_tube.cpp +++ b/src/problems/RadTube/test_radiation_tube.cpp @@ -56,10 +56,10 @@ template <> struct Physics_Traits { // number of radiation groups static constexpr int nGroups = 2; static constexpr UnitSystem unit_system = UnitSystem::CGS; - static constexpr double c_hat = 10.0 * a0; }; template <> struct RadSystem_Traits { + 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/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index d0eab6af5..42e7c88d2 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -56,6 +56,7 @@ constexpr double shock_position = 0.0130; // 0.0132; // cm // we initialize slightly to the left...) template <> struct RadSystem_Traits { + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -78,7 +79,6 @@ template <> struct Physics_Traits { static constexpr double boltzmann_constant = k_B; static constexpr double gravitational_constant = Ggrav; static constexpr double c_light = c; - static constexpr double c_hat = chat; static constexpr double radiation_constant = a_rad; }; diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 59578ab50..9555be7ae 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -55,6 +55,7 @@ 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_hat_over_c = chat / c; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; }; @@ -75,7 +76,6 @@ template <> struct Physics_Traits { 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 c_hat = chat; static constexpr double unit_length = 1.0; // cm static constexpr double unit_mass = 1.0; // g static constexpr double unit_time = 1.0; // s diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index b98b7250a..28ce70489 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -72,7 +72,7 @@ enum class OpacityModel { // this struct is specialized by the user application code // template struct RadSystem_Traits { - static constexpr double c_hat = c_light_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}; @@ -201,7 +201,7 @@ template class RadSystem : public HyperbolicSystem::c_hat; + 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) { From 2add5798c7b4664a3fee7c0cc3d198534288e161 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 11:47:39 +1100 Subject: [PATCH 19/34] cleanup --- src/problems/RadLineCooling/test_rad_line_cooling.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index 0d7de1b1f..f40c02d3f 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -66,6 +66,8 @@ template <> struct Physics_Traits { // 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 { @@ -73,9 +75,6 @@ template <> struct RadSystem_Traits { static constexpr int beta_order = 0; static constexpr double energy_unit = nu_unit; static constexpr double c_hat_over_c = chat / c; - // Equivalently, set - // static constexpr double c_light = c; - // static constexpr double radiation_constant = a_rad; }; template <> struct ISM_Traits { From ea0c0ee803d5882719a527e5fa4507ec0c3f758e Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 12:47:06 +1100 Subject: [PATCH 20/34] remove print messages in MarshakAsymptotic --- .../test_radiation_marshak_asymptotic.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index e32e29abd..f28a6f1ef 100644 --- a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -206,14 +206,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) { From ff75400b3735dfb3cf3834f1e3852545b0f42ef8 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 13:48:14 +1100 Subject: [PATCH 21/34] minor fixes --- src/hydro/EOS.hpp | 7 +++---- src/hydro/NSCBC_inflow.hpp | 2 +- src/simulation.hpp | 4 +++- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index b629a2b22..47c923871 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,6 @@ 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 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; 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/simulation.hpp b/src/simulation.hpp index b39e8ea04..64da1fb6e 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -529,7 +529,9 @@ template void AMRSimulation::initialize() } } - initializeSimulationMetadata(); + if constexpr (Physics_Traits::is_hydro_enabled || Physics_Traits::is_radiation_enabled) { + initializeSimulationMetadata(); + } #ifdef AMREX_USE_ASCENT // initialize Ascent From 6f84a04a136d86fe8d6073a58c955ffdc09c5df6 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 13:48:31 +1100 Subject: [PATCH 22/34] update all tests to adapt new unit systems --- src/problems/Advection/test_advection.cpp | 1 + .../test_advection_semiellipse.cpp | 1 + src/problems/BinaryOrbitCIC/binary_orbit.cpp | 2 +- src/problems/Cooling/test_cooling.cpp | 2 +- src/problems/FCQuantities/test_fc_quantities.cpp | 2 +- src/problems/HydroBlast2D/test_hydro2d_blast.cpp | 2 +- src/problems/HydroBlast3D/test_hydro3d_blast.cpp | 2 +- src/problems/HydroContact/test_hydro_contact.cpp | 2 +- .../HydroHighMach/test_hydro_highmach.cpp | 2 +- .../HydroKelvinHelmholz/test_hydro2d_kh.cpp | 2 +- src/problems/HydroLeblanc/test_hydro_leblanc.cpp | 3 ++- src/problems/HydroQuirk/test_quirk.cpp | 2 +- .../HydroRichtmeyerMeshkov/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 +- src/problems/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 | 2 +- src/problems/ODEIntegration/test_ode.cpp | 1 - src/problems/PassiveScalar/test_scalars.cpp | 2 +- .../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 ++++++------ src/problems/RadForce/test_radiation_force.cpp | 6 ++---- .../test_rad_line_cooling_MG.cpp | 16 ++++++++-------- .../RadMarshak/test_radiation_marshak.cpp | 9 ++++++--- .../test_radiation_marshak_asymptotic.cpp | 6 ++---- .../RadMarshakCGS/test_radiation_marshak_cgs.cpp | 8 +++----- .../test_radiation_marshak_dust.cpp | 14 +++++++------- .../test_radiation_marshak_dust_and_PE.cpp | 11 ++++++----- .../test_radiation_marshak_Vaytet.cpp | 6 ++---- .../test_radiation_matter_coupling.cpp | 8 +++----- .../test_radiation_matter_coupling_rsla.cpp | 11 +++++------ src/problems/RadPulse/test_radiation_pulse.cpp | 10 ++++++---- src/problems/RadShadow/test_radiation_shadow.cpp | 10 ++++------ .../RadStreaming/test_radiation_streaming.cpp | 10 ++++++---- .../RadStreamingY/test_radiation_streaming_y.cpp | 11 ++++++----- .../RadSuOlson/test_radiation_SuOlson.cpp | 11 ++++++----- src/problems/RadTophat/test_radiation_tophat.cpp | 10 ++++------ src/problems/RadhydroBB/test_radhydro_bb.cpp | 10 ++++++---- .../RadhydroPulse/test_radhydro_pulse.cpp | 12 ++++-------- .../RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp | 13 ++++--------- .../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 | 10 ++++------ .../test_radhydro_shock_multigroup.cpp | 6 ++---- .../test_radhydro_uniform_advecting.cpp | 12 +++++++----- src/problems/RandomBlast/blast.cpp | 2 +- .../RayleighTaylor2D/test_hydro2d_rt.cpp | 2 +- .../RayleighTaylor3D/test_hydro3d_rt.cpp | 2 +- src/problems/ShockCloud/cloud.cpp | 2 +- .../SphericalCollapse/spherical_collapse.cpp | 2 +- src/problems/StarCluster/star_cluster.cpp | 2 +- 57 files changed, 158 insertions(+), 184 deletions(-) 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/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/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..2e8b9f9a3 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 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/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/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 f28a6f1ef..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 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..1d1cc3419 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -20,7 +20,6 @@ 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 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,7 @@ 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 auto Frads = Erads * c; 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..1f6abd822 100644 --- a/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp +++ b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp @@ -25,7 +25,6 @@ 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 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..2afaa904f 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 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..4d12064dc 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -22,13 +22,11 @@ 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 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 +39,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 = 1.0; 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..e57bc0fca 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/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..a6f758542 100644 --- a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -21,7 +21,6 @@ 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; @@ -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..a1774c3f8 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -36,8 +36,8 @@ 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 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 @@ -46,14 +46,11 @@ 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/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..732014f98 100644 --- a/src/problems/RandomBlast/blast.cpp +++ b/src/problems/RandomBlast/blast.cpp @@ -47,12 +47,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 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..343268994 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,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/StarCluster/star_cluster.cpp b/src/problems/StarCluster/star_cluster.cpp index 869c6c5c7..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 { @@ -55,6 +54,7 @@ template <> struct Physics_Traits { 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; }; From b21d2d982d4e154125075ae358b20fa5564eddb6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 28 Oct 2024 02:49:03 +0000 Subject: [PATCH 23/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/physics_info.hpp | 4 ++-- .../RadMarshakDust/test_radiation_marshak_dust.cpp | 2 +- .../test_radiation_marshak_dust_and_PE.cpp | 2 +- src/problems/RadShadow/test_radiation_shadow.cpp | 2 +- src/problems/RadTophat/test_radiation_tophat.cpp | 2 +- src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp | 2 +- src/problems/RadhydroShell/test_radhydro_shell.cpp | 8 ++++---- 7 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/physics_info.hpp b/src/physics_info.hpp index 8618b9136..5597d0edb 100644 --- a/src/physics_info.hpp +++ b/src/physics_info.hpp @@ -21,8 +21,8 @@ template struct Physics_Traits { 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 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; diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp index 1d1cc3419..0ab782b82 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -19,7 +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 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 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 1f6abd822..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,7 +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 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 diff --git a/src/problems/RadShadow/test_radiation_shadow.cpp b/src/problems/RadShadow/test_radiation_shadow.cpp index 2afaa904f..f476ebd8f 100644 --- a/src/problems/RadShadow/test_radiation_shadow.cpp +++ b/src/problems/RadShadow/test_radiation_shadow.cpp @@ -27,7 +27,7 @@ constexpr double T_hohlraum = 1740.; // K constexpr double T_initial = 290.; // K constexpr double a_rad = C::a_rad; // erg cm^-3 K^-4 -constexpr double c = C::c_light; // cm s^-1 +constexpr double c = C::c_light; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = 10. * C::m_u; diff --git a/src/problems/RadTophat/test_radiation_tophat.cpp b/src/problems/RadTophat/test_radiation_tophat.cpp index e57bc0fca..b0f8adbf7 100644 --- a/src/problems/RadTophat/test_radiation_tophat.cpp +++ b/src/problems/RadTophat/test_radiation_tophat.cpp @@ -32,7 +32,7 @@ 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 = C::a_rad; // erg cm^-3 K^-4 -constexpr double c = C::c_light; // cm s^-1 +constexpr double c = C::c_light; // cm s^-1 template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = C::m_u; diff --git a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp index a6f758542..13cdb2768 100644 --- a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -21,7 +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 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; diff --git a/src/problems/RadhydroShell/test_radhydro_shell.cpp b/src/problems/RadhydroShell/test_radhydro_shell.cpp index a1774c3f8..5e60a3592 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -38,10 +38,10 @@ constexpr bool simulate_full_box = true; 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 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 { From 9446336f59260f8be143aafd2cc4d96e3ac782c6 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 14:31:24 +1100 Subject: [PATCH 24/34] remove some static_assert --- src/hydro/EOS.hpp | 2 -- src/radiation/radiation_system.hpp | 6 ------ src/simulation.hpp | 2 -- 3 files changed, 10 deletions(-) diff --git a/src/hydro/EOS.hpp b/src/hydro/EOS.hpp index 47c923871..a54cd0162 100644 --- a/src/hydro/EOS.hpp +++ b/src/hydro/EOS.hpp @@ -78,8 +78,6 @@ template class EOS 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); - } else { - static_assert(false, "Invalid unit system"); } }(); }; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 28ce70489..cdf778919 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -197,8 +197,6 @@ template class RadSystem : public HyperbolicSystem::unit_system == UnitSystem::CUSTOM) { // c / c_bar = u_l / u_t return c_light_cgs_ / (Physics_Traits::unit_length / Physics_Traits::unit_time); - } else { - static_assert(false, "Invalid unit system"); } }(); static constexpr double c_hat_ = c_light_ * RadSystem_Traits::c_hat_over_c; @@ -214,8 +212,6 @@ template class RadSystem : public HyperbolicSystem::unit_time * Physics_Traits::unit_time) / (Physics_Traits::unit_temperature * Physics_Traits::unit_temperature * Physics_Traits::unit_temperature * Physics_Traits::unit_temperature)); - } else { - static_assert(false, "Invalid unit system"); } }(); @@ -265,8 +261,6 @@ template class RadSystem : public HyperbolicSystem::unit_length * Physics_Traits::unit_length * Physics_Traits::unit_mass / (Physics_Traits::unit_time * Physics_Traits::unit_time) / Physics_Traits::unit_temperature); - } else { - static_assert(false, "Invalid unit system"); } }(); diff --git a/src/simulation.hpp b/src/simulation.hpp index 64da1fb6e..08d64debf 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -402,8 +402,6 @@ template class AMRSimulation : public amrex::AmrCore 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)); - } else { - static_assert(false, "Invalid unit system"); } }(); From 57efbf2ef8b553c20a127df576d406fed84de7b1 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 14:47:57 +1100 Subject: [PATCH 25/34] fix advection2d --- src/problems/Advection2D/test_advection2d.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/problems/Advection2D/test_advection2d.cpp b/src/problems/Advection2D/test_advection2d.cpp index 4b9f0b1d9..51c45f8fe 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, From 0ad4f466cdf70a06663ffe74c0b94025fd1db3a9 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 15:10:49 +1100 Subject: [PATCH 26/34] fix boltzmann_constant errors --- src/problems/NSCBC/vortex.cpp | 2 +- src/problems/RadShadow/test_radiation_shadow.cpp | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/problems/NSCBC/vortex.cpp b/src/problems/NSCBC/vortex.cpp index 2e8b9f9a3..1fa9c79ba 100644 --- a/src/problems/NSCBC/vortex.cpp +++ b/src/problems/NSCBC/vortex.cpp @@ -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/RadShadow/test_radiation_shadow.cpp b/src/problems/RadShadow/test_radiation_shadow.cpp index f476ebd8f..d99e1d41e 100644 --- a/src/problems/RadShadow/test_radiation_shadow.cpp +++ b/src/problems/RadShadow/test_radiation_shadow.cpp @@ -226,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"; From 7b4b8bc08c882cd7f269c63dad4b2803b825a4c4 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 15:30:22 +1100 Subject: [PATCH 27/34] fix chat error --- src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp | 3 ++- src/problems/RadStreamingY/test_radiation_streaming_y.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp index 0ab782b82..e5f5d8b69 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -152,7 +152,8 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect // const auto Erads = RadSystem::ComputeThermalRadiation(T_rad_L, radBoundaries_); quokka::valarray const Erads = {erad_floor, EradL}; - const auto Frads = Erads * c; + const double c_device = c; + const auto Frads = Erads * c_device; if (i < lo[0]) { // streaming inflow boundary diff --git a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp index 4d12064dc..8c83362fe 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -22,6 +22,7 @@ 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 = 0.3; // reduced speed of light constexpr double kappa0 = 1.0e-10; // opacity constexpr double rho = 1.0; @@ -47,7 +48,7 @@ template <> struct Physics_Traits { }; template <> struct RadSystem_Traits { - static constexpr double c_hat_over_c = 1.0; + static constexpr double c_hat_over_c = chat / c; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; }; From 36f0f27779927ee06f62605218539e2f65b4d386 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 28 Oct 2024 04:30:58 +0000 Subject: [PATCH 28/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/problems/RadStreamingY/test_radiation_streaming_y.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp index 8c83362fe..94f43e9ad 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -22,7 +22,7 @@ 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 = 0.3; // 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; From 8d3a70ca2d988fc83e1b7f673f65791c02055942 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Mon, 28 Oct 2024 16:07:49 +1100 Subject: [PATCH 29/34] fix bugs --- src/problems/Advection2D/test_advection2d.cpp | 7 ++++++- src/problems/PopIII/popiii.cpp | 1 + 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/problems/Advection2D/test_advection2d.cpp b/src/problems/Advection2D/test_advection2d.cpp index 51c45f8fe..7810d1831 100644 --- a/src/problems/Advection2D/test_advection2d.cpp +++ b/src/problems/Advection2D/test_advection2d.cpp @@ -107,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/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 { From 27ae3cd2159faded4979010c66b3b65327a0662a Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Tue, 29 Oct 2024 17:08:38 +1100 Subject: [PATCH 30/34] add static_assert unit_system == CGS in RandomBlast test --- src/problems/RandomBlast/blast.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/problems/RandomBlast/blast.cpp b/src/problems/RandomBlast/blast.cpp index 732014f98..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; @@ -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; From 6f7a4bc0c84bc472b432f75c3abb199f6955b820 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Wed, 30 Oct 2024 12:31:47 +1100 Subject: [PATCH 31/34] fix bug --- tests/RadStreamingY.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From f92913348b7c76cab7d9def2afa3fd8fab435d53 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 31 Oct 2024 11:15:35 +1100 Subject: [PATCH 32/34] set G = 1 in SphericalCollapse --- src/problems/SphericalCollapse/spherical_collapse.cpp | 4 +++- tests/SphericalCollapse.in | 2 -- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/problems/SphericalCollapse/spherical_collapse.cpp b/src/problems/SphericalCollapse/spherical_collapse.cpp index 343268994..c8085bf36 100644 --- a/src/problems/SphericalCollapse/spherical_collapse.cpp +++ b/src/problems/SphericalCollapse/spherical_collapse.cpp @@ -44,7 +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::CGS; + static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS; + static constexpr double boltzmann_constant = C::k_B; + static constexpr double gravity_constant = 1.0; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) 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 From 226a40121ad6e1f53b51e34af71609738bf7fce6 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 1 Nov 2024 13:24:49 +1100 Subject: [PATCH 33/34] fix minor bug --- src/problems/SphericalCollapse/spherical_collapse.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/SphericalCollapse/spherical_collapse.cpp b/src/problems/SphericalCollapse/spherical_collapse.cpp index c8085bf36..064ffe735 100644 --- a/src/problems/SphericalCollapse/spherical_collapse.cpp +++ b/src/problems/SphericalCollapse/spherical_collapse.cpp @@ -46,7 +46,7 @@ template <> struct Physics_Traits { 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 gravity_constant = 1.0; + static constexpr double gravitational_constant = 1.0; }; template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) From 51cf4960143e5c34f0ab663997aa8ee22114a35a Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 2 Nov 2024 11:20:03 +1100 Subject: [PATCH 34/34] pltfile from line_cooling --- src/problems/RadLineCooling/test_rad_line_cooling.cpp | 2 +- tests/RadLineCooling.in | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index f40c02d3f..4f12b5012 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -211,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/tests/RadLineCooling.in b/tests/RadLineCooling.in index b77624b49..e0f16b7ac 100644 --- a/tests/RadLineCooling.in +++ b/tests/RadLineCooling.in @@ -23,6 +23,4 @@ radiation.dust_gas_interaction_coeff = 1e-20 do_reflux = 0 do_subcycle = 0 -suppress_output = 1 - -plotfile_interval = 10000000 \ No newline at end of file +suppress_output = 1 \ No newline at end of file