From 861a48986260c20e51d0e5bfac8d8d19029ed10e Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Sat, 18 Jan 2025 10:49:35 -0700 Subject: [PATCH 1/6] Add temperature to units. Add temperature conversion factor to constants. Have base unit specifiers be built on top of preset base values --- src/utils/units.cpp | 37 +++++++++++++++++++++---------------- src/utils/units.hpp | 10 ++++++++-- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/src/utils/units.cpp b/src/utils/units.cpp index 1ac899d..4522701 100644 --- a/src/utils/units.cpp +++ b/src/utils/units.cpp @@ -32,25 +32,29 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { } else { PARTHENON_FAIL("Physical unit system not recognized! Choices are [scalefree, cgs]"); } - - if (physical_units_ == PhysicalUnits::scalefree) { - length_ = 1.; - time_ = 1.; - mass_ = 1.; - } else { + length_ = 1.; + time_ = 1.; + mass_ = 1.; + temp_ = 1.; + if (physical_units_ != PhysicalUnits::scalefree) { std::string unit_conversion = pin->GetOrAddString("artemis", "unit_conversion", "base"); - if (unit_conversion == "base") { - length_ = pin->GetOrAddReal("artemis", "length", 1.); - time_ = pin->GetOrAddReal("artemis", "time", 1.); - mass_ = pin->GetOrAddReal("artemis", "mass", 1.); - } else if (unit_conversion == "ppd") { + if (unit_conversion == "ppd") { length_ = AU; mass_ = Msolar; time_ = Year / (2. * M_PI); + } else if (unit_conversion == "base") { + // do nothing } else { PARTHENON_FAIL("Unit conversion not recognized! Choices are [base, ppd]"); } + // not that these multiplied by whatever values were previously set. + // For example, if unit_conversion=ppd and mass = 10.0, then + // that sets mass_ to 10 MSolar. + length_ *= pin->GetOrAddReal("artemis", "length", 1.); + time_ *= pin->GetOrAddReal("artemis", "time", 1.); + mass_ *= pin->GetOrAddReal("artemis", "mass", 1.); + temp_ *= pin->GetOrAddReal("artemis", "temperature", 1.); } // Remaining conversion factors @@ -99,16 +103,17 @@ Constants::Constants(Units &units) { const Real length = units.GetLengthCodeToPhysical(); const Real time = units.GetTimeCodeToPhysical(); const Real mass = units.GetMassCodeToPhysical(); + const Real temp = units.GetTemperatureCodeToPhysical(); + const Real energy = mass * std::pow(length / time, 2); // Convert constants to code units G_code_ = G_ * std::pow(length, -3) / mass * std::pow(time, 2); - kb_code_ = - kb_ * std::pow(time, 2) / mass * std::pow(length, -2); // 1 K = 1 code unit temp + kb_code_ = kb_ * temp / energy; c_code_ = c_ * time / length; - h_code_ = h_ * time / mass * std::pow(length, -2); - ar_code_ = ar_ * length * time * time / mass; + h_code_ = h_ / (energy * time); + ar_code_ = ar_ * std::pow(temp, 4) * std::pow(length, 3) / energy; amu_code_ = amu_ / mass; - eV_code_ = eV_ * std::pow(time, 2) / mass * std::pow(length, -2); + eV_code_ = eV_ / energy; Msolar_code_ = Msolar_ / mass; AU_code_ = AU_ / length; Rjup_code_ = Rjup_ / length; diff --git a/src/utils/units.hpp b/src/utils/units.hpp index 0aa05cb..60eb700 100644 --- a/src/utils/units.hpp +++ b/src/utils/units.hpp @@ -52,6 +52,11 @@ class Units { KOKKOS_INLINE_FUNCTION Real GetMassPhysicalToCode() const { return 1. / mass_; } + KOKKOS_INLINE_FUNCTION + Real GetTemperatureCodeToPhysical() const { return temp_; } + KOKKOS_INLINE_FUNCTION + Real GetTemperaturePhysicalToCode() const { return 1. / temp_; } + KOKKOS_INLINE_FUNCTION Real GetSpeedCodeToPhysical() const { return length_ / time_; } KOKKOS_INLINE_FUNCTION @@ -83,9 +88,9 @@ class Units { Real GetOpacityPhysicalToCode() const { return mass_ / (length_ * length_); } KOKKOS_INLINE_FUNCTION - Real GetSpecificHeatCodeToPhysical() const { return energy_ / mass_; } + Real GetSpecificHeatCodeToPhysical() const { return energy_ / (mass_ * temp_); } KOKKOS_INLINE_FUNCTION - Real GetSpecificHeatPhysicalToCode() const { return mass_ / energy_; } + Real GetSpecificHeatPhysicalToCode() const { return mass_ * temp_ / energy_; } inline std::string GetSystemName() const { return (physical_units_ == PhysicalUnits::scalefree) ? "Scale free" : "CGS"; @@ -98,6 +103,7 @@ class Units { Real length_; Real time_; Real mass_; + Real temp_; Real energy_; Real number_density_; From 1beaf57b7fb34d44b07a285cab611e2e68bcac11 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Sat, 18 Jan 2025 10:50:01 -0700 Subject: [PATCH 2/6] Unit system aware singularity eos --- src/gas/gas.cpp | 6 +++++- src/utils/eos/eos.hpp | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index ffbef0d..62ad27e 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -113,7 +113,11 @@ std::shared_ptr Initialize(ParameterInput *pin, PARTHENON_REQUIRE(mu > 0, "Only positive mean molecular weight allowed!"); cv = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * mu); } - EOS eos_host = singularity::IdealGas(gamma - 1., cv); + EOS eos_host = singularity::UnitSystem( + singularity::IdealGas(gamma - 1., cv), + singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimeCodeToPhysical(), + units.GetMassCodeToPhysical(), units.GetLengthCodeToPhysical(), + units.GetTemperatureCodeToPhysical()); EOS eos_device = eos_host.GetOnDevice(); params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); diff --git a/src/utils/eos/eos.hpp b/src/utils/eos/eos.hpp index 85f5b8e..d1703fa 100644 --- a/src/utils/eos/eos.hpp +++ b/src/utils/eos/eos.hpp @@ -22,7 +22,7 @@ namespace ArtemisUtils { static constexpr int lambda_max_vals = 1; // Variant containing all EOSs to be used in Artemis. -using EOS = singularity::Variant; +using EOS = singularity::Variant>; } // namespace ArtemisUtils From 633ab213260742a9b821c314eb19d319199840cf Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Sun, 19 Jan 2025 14:52:26 -0700 Subject: [PATCH 3/6] Add temperature conversion to opacities. --- src/gas/gas.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 62ad27e..41af71f 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -133,11 +133,12 @@ std::shared_ptr Initialize(ParameterInput *pin, const Real length = units.GetLengthCodeToPhysical(); const Real time = units.GetTimeCodeToPhysical(); const Real mass = units.GetMassCodeToPhysical(); + const Real temp = units.GetTemperatureCodeToPhysical(); if (opacity_model_name == "none") { - opacity = NonCGSUnits(Gray(0.0), time, mass, length, 1.); + opacity = NonCGSUnits(Gray(0.0), time, mass, length, temp); } else if (opacity_model_name == "constant") { const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0); - opacity = NonCGSUnits(Gray(kappa_a), time, mass, length, 1.); + opacity = NonCGSUnits(Gray(kappa_a), time, mass, length, temp); } else if (opacity_model_name == "shocktube_a") { const Real coef_kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0); @@ -159,10 +160,10 @@ std::shared_ptr Initialize(ParameterInput *pin, std::string scattering_model_name = pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none"); if (scattering_model_name == "none") { - scattering = NonCGSUnitsS(GrayS(0.0, 1.0), time, mass, length, 1.); + scattering = NonCGSUnitsS(GrayS(0.0, 1.0), time, mass, length, temp); } else if (scattering_model_name == "constant") { const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); - scattering = NonCGSUnitsS(GrayS(kappa_s, 1.0), time, mass, length, 1.); + scattering = NonCGSUnitsS(GrayS(kappa_s, 1.0), time, mass, length, temp); } else { PARTHENON_FAIL("Scattering model not recognized!"); } From cf6fce23648e810be36622c68becc0bfc8b6a233 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 21 Jan 2025 17:28:13 -0700 Subject: [PATCH 4/6] Add temperature unit to splash --- src/utils/artemis_utils.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/artemis_utils.cpp b/src/utils/artemis_utils.cpp index 573c2ed..7be05b8 100644 --- a/src/utils/artemis_utils.cpp +++ b/src/utils/artemis_utils.cpp @@ -53,6 +53,7 @@ void PrintArtemisConfiguration(Packages_t &packages) { printf(" [L] = %.2e\n", units.GetLengthCodeToPhysical()); printf(" [M] = %.2e\n", units.GetMassCodeToPhysical()); printf(" [T] = %.2e\n", units.GetTimeCodeToPhysical()); + printf(" [K] = %.2e\n", units.GetTemperatureCodeToPhysical()); printf(" Active physics: %s", msg.c_str()); if (params.Get("do_nbody")) { From 02d61f1a79aa26b5acd65622ee370944bebbe9ca Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 21 Jan 2025 17:31:08 -0700 Subject: [PATCH 5/6] Revert allowing setting of base scale factors in non base specifiers --- src/utils/units.cpp | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/utils/units.cpp b/src/utils/units.cpp index 4522701..da9809e 100644 --- a/src/utils/units.cpp +++ b/src/utils/units.cpp @@ -32,29 +32,28 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { } else { PARTHENON_FAIL("Physical unit system not recognized! Choices are [scalefree, cgs]"); } - length_ = 1.; - time_ = 1.; - mass_ = 1.; - temp_ = 1.; - if (physical_units_ != PhysicalUnits::scalefree) { + + if (physical_units_ == PhysicalUnits::scalefree) { + length_ = 1.; + time_ = 1.; + mass_ = 1.; + temp_ = 1.; + } else { std::string unit_conversion = pin->GetOrAddString("artemis", "unit_conversion", "base"); - if (unit_conversion == "ppd") { + if (unit_conversion == "base") { + length_ = pin->GetOrAddReal("artemis", "length", 1.); + time_ = pin->GetOrAddReal("artemis", "time", 1.); + mass_ = pin->GetOrAddReal("artemis", "mass", 1.); + temp_ = pin->GetOrAddReal("artemis", "temperature", 1.); + } else if (unit_conversion == "ppd") { length_ = AU; mass_ = Msolar; time_ = Year / (2. * M_PI); - } else if (unit_conversion == "base") { - // do nothing + temp_ = 1.0; } else { PARTHENON_FAIL("Unit conversion not recognized! Choices are [base, ppd]"); } - // not that these multiplied by whatever values were previously set. - // For example, if unit_conversion=ppd and mass = 10.0, then - // that sets mass_ to 10 MSolar. - length_ *= pin->GetOrAddReal("artemis", "length", 1.); - time_ *= pin->GetOrAddReal("artemis", "time", 1.); - mass_ *= pin->GetOrAddReal("artemis", "mass", 1.); - temp_ *= pin->GetOrAddReal("artemis", "temperature", 1.); } // Remaining conversion factors From 94f09a4a385303b19b8af9b6967991a781c9f685 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 21 Jan 2025 17:34:26 -0700 Subject: [PATCH 6/6] Update params for new temperature unit --- src/artemis_params.yaml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/artemis_params.yaml b/src/artemis_params.yaml index e61ddb0..73e9d2e 100644 --- a/src/artemis_params.yaml +++ b/src/artemis_params.yaml @@ -33,13 +33,13 @@ artemis: cgs: _type: opt _description: "CGS unit system" - unit_conversion: + unit_specifier: _type: "string" _default: "base" _description: "How to provide unit conversions between code and physical units" base: _type: opt - _description: "Provide base unit conversions (length, time, mass)" + _description: "Provide base unit conversions (length, time, mass, temperature)" ppd: _type: opt _description: "AU, Year/(2 pi), solar mass units for protoplanetary disks" @@ -55,6 +55,10 @@ artemis: _type: "Real" _default: "1.0" _description: "Physical units value of mass equal to 1 code unit" + temperature: + _type: "Real" + _default: "1.0" + _description: "Physical units value of temperature equal to 1 code unit" physics: