Temperature unit

CHANGELOG.md

@@ -9,7 +9,8 @@
 - [[PR638]](https://github.com/lanl/singularity-eos/pull/638) Fix Bugs in Davis Reactants and Products.  Bulk modulus is now isentropic.
 
 ### Fixed (Repair bugs, etc)
-
+- [[PR63]](https://github.com/lanl/singularity-eos/pull/639) Fixed UnitSystem temperature bug. Temperature is now treated in the same as the time, mass, and length unit factors.
+
 ### Changed (changing behavior/API/variables/...)
 
 ### Infrastructure (changes irrelevant to downstream codes)

singularity-eos/eos/modifiers/eos_unitsystem.hpp

@@ -87,7 +87,8 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
       // as the conversion factors of the base units.
       : UnitSystem(std::forward<T>(t), eos_units_init::thermal_units_init_tag,
                    (length_unit * length_unit * length_unit) / mass_unit,
-                   (time_unit * time_unit) / (length_unit * length_unit), temp_unit) {}
+                   (time_unit * time_unit) / (length_unit * length_unit),
+                   1.0 / temp_unit) {}
   UnitSystem(T &&t, const Real rho_unit, const Real sie_unit, const Real temp_unit)
       : UnitSystem(std::forward<T>(t), eos_units_init::thermal_units_init_tag, rho_unit,
                    sie_unit, temp_unit) {}

test/test_eos_modifiers.cpp

@@ -381,3 +381,52 @@ SCENARIO("Modifiers propagate introspection bounds correctly", "[Modifiers]") {
     }
   }
 }
+SCENARIO("UnitSystem  modifier converts units correctly",
+         "[Modifiers][Units][IdealGas]") {
+  GIVEN("An IdealGas EOS in non-cgs units") {
+    constexpr Real mu = 2.3;
+    constexpr Real gm1 = 0.4;
+    constexpr Real amu = 1.6605e-24;       // g
+    constexpr Real kb = 1.3807e-16;        // cm^2 g / s^2 / K
+    constexpr Real kbmu = kb / (mu * amu); // cm^2 / s^2 /K
+    constexpr Real cv = kbmu / gm1;
+
+    // conversions from code to cgs
+    constexpr Real length = 1.495978707e14;     // 10 AU
+    constexpr Real mass = 1.988416e33;          // Msun
+    constexpr Real time = 1.5871820329271033e8; // year/(2*pi)*((10 AU)**3/msun)**0.5
+    constexpr Real temperature = 2.4753531e4;   // Such that  T = vkep**2/(kb*mu*mp)
+
+    constexpr Real kbmu_code =
+        kbmu * (time / length) * (time / length) * temperature; // L^2/[t^2 T]
+    constexpr Real cv_code = kbmu_code / gm1;
+
+    auto eos =
+        UnitSystem<IdealGas>(IdealGas(gm1, cv), eos_units_init::LengthTimeUnitsInit(),
+                             1. / time, 1. / mass, 1. / length, 1. / temperature);
+
+    // P = rho kbmu T
+    // e = cv * T
+    // cs^2 = gamma*P/rho
+    constexpr Real rho_code = 1.0;
+    constexpr Real temp_code = 1.0;
+    constexpr Real pres_code = rho_code * temp_code * kbmu_code;
+    constexpr Real e_code = cv_code * temp_code;
+
+    THEN("Pressure is correct") {
+
+      REQUIRE(isClose(eos.PressureFromDensityTemperature(rho_code, temp_code), pres_code,
+                      1.e-12));
+    }
+    THEN("Specific Heat is correct") {
+
+      REQUIRE(isClose(eos.SpecificHeatFromDensityTemperature(rho_code, temp_code),
+                      cv_code, 1.e-12));
+    }
+    THEN("SIE is correct") {
+
+      REQUIRE(isClose(eos.InternalEnergyFromDensityTemperature(rho_code, temp_code),
+                      e_code, 1e-12));
+    }
+  }
+}