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: diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index ffbef0d..41af71f 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); @@ -129,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); @@ -155,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!"); } 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")) { 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 diff --git a/src/utils/units.cpp b/src/utils/units.cpp index 1ac899d..da9809e 100644 --- a/src/utils/units.cpp +++ b/src/utils/units.cpp @@ -37,6 +37,7 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { length_ = 1.; time_ = 1.; mass_ = 1.; + temp_ = 1.; } else { std::string unit_conversion = pin->GetOrAddString("artemis", "unit_conversion", "base"); @@ -44,10 +45,12 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { 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); + temp_ = 1.0; } else { PARTHENON_FAIL("Unit conversion not recognized! Choices are [base, ppd]"); } @@ -99,16 +102,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_;