Merge pull request #47 from lanl/dempsey/units

Add temperature to units and have units aware EOS

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:

src/gas/gas.cpp

@@ -113,7 +113,11 @@ std::shared_ptr<StateDescriptor> 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>(
+        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<StateDescriptor> 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>(Gray(0.0), time, mass, length, 1.);
+    opacity = NonCGSUnits<Gray>(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>(Gray(kappa_a), time, mass, length, 1.);
+    opacity = NonCGSUnits<Gray>(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<StateDescriptor> Initialize(ParameterInput *pin,
   std::string scattering_model_name =
       pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none");
   if (scattering_model_name == "none") {
-    scattering = NonCGSUnitsS<GrayS>(GrayS(0.0, 1.0), time, mass, length, 1.);
+    scattering = NonCGSUnitsS<GrayS>(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>(GrayS(kappa_s, 1.0), time, mass, length, 1.);
+    scattering = NonCGSUnitsS<GrayS>(GrayS(kappa_s, 1.0), time, mass, length, temp);
   } else {
     PARTHENON_FAIL("Scattering model not recognized!");
   }

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<bool>("do_nbody")) {

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<singularity::IdealGas>;
+using EOS = singularity::Variant<singularity::UnitSystem<singularity::IdealGas>>;
 
 } // namespace ArtemisUtils
 

src/utils/units.cpp

@@ -37,17 +37,20 @@ Units::Units(ParameterInput *pin, std::shared_ptr<StateDescriptor> pkg) {
     length_ = 1.;
     time_ = 1.;
     mass_ = 1.;
+    temp_ = 1.;
   } else {
     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.);
+      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;

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_;