Merge pull request #265 from lanl/rberger_unit_system_modifier

Add missing UnitSystem to variant and EOSPAC

CHANGELOG.md

@@ -21,6 +21,7 @@
 - [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
 
 ### Added (new features/APIs/variables/...)
+- [[PR265]](https://github.com/lanl/singularity-eos/pull/265) Add missing UnitSystem modifier combinations to variant and EOSPAC
 - [[PR279]](https://github.com/lanl/singularity-eos/pull/279) added noble-abel EoS
 - [[PR274]](https://github.com/lanl/singularity-eos/pull/274) added a stiffened gas EoS
 - [[PR254]](https://github.com/lanl/singularity-eos/pull/254) added ability to peel off modifiers as needed

python/modifier_scaled.cpp

@@ -32,7 +32,4 @@ void create_scaled(py::module_ &m, tl<Ts...>) {
   auto l = {(scaled_eos_class<Ts>(m), 0)...};
 }
 
-void create_scaled_eos_classes(py::module_ &m) {
-  create_scaled(m, singularity::scaled);
-  create_scaled(m, singularity::scaled_of_shifted);
-}
+void create_scaled_eos_classes(py::module_ &m) { create_scaled(m, singularity::scaled); }

python/module.cpp

@@ -119,11 +119,11 @@ PYBIND11_MODULE(singularity_eos, m) {
     .def(py::init<int, bool>(), py::arg("matid"), py::arg("invert_at_setup")=false);
 #endif
 
+  create_unit_system_eos_classes(m);
+  create_relativistic_eos_classes(m);
   create_shifted_eos_classes(m);
   create_scaled_eos_classes(m);
   create_bilinear_ramp_eos_classes(m);
-  create_relativistic_eos_classes(m);
-  create_unit_system_eos_classes(m);
 
   m.def("pAlpha2BilinearRampParams", [](const IdealGas &eos, const Real alpha0, const Real Pe, const Real Pc){
     Real r0, a, b, c;

singularity-eos/eos/eos.hpp

@@ -77,38 +77,60 @@ static constexpr const auto full_eos_list =
        EOSPAC
 #endif // SINGULARITY_USE_EOSPAC
        >{};
-// eos's that get relativistic and unit system modifiers
-static constexpr const auto partial_eos_list =
+// eos's that get relativistic modifier
+static constexpr const auto relativistic_eos_list =
     tl<IdealGas
 #ifdef SINGULARITY_USE_SPINER_WITH_HDF5
        ,
        SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse
 #endif // SINGULAIRTY_USE_SPINER_WITH_HDF5
        >{};
+// eos's that get unit system modifier
+static constexpr const auto unit_system_eos_list =
+    tl<IdealGas
+#ifdef SPINER_USE_HDF
+       ,
+       SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse
+#endif // SPINER_USE_HDF
+#ifdef SINGULARITY_USE_EOSPAC
+       ,
+       EOSPAC
+#endif // SINGULARITY_USE_EOSPAC
+       >{};
 // modifiers that get applied to all eos's
 static constexpr const auto apply_to_all = al<ScaledEOS, ShiftedEOS>{};
-// modifiers thet get applied to a subset of eos's
-static constexpr const auto apply_to_partial = al<UnitSystem, RelativisticEOS>{};
-// variadic list of eos's with shifted or scaled modifiers
-static constexpr const auto shifted =
-    transform_variadic_list(full_eos_list, al<ShiftedEOS>{});
-static constexpr const auto scaled =
-    transform_variadic_list(full_eos_list, al<ScaledEOS>{});
-// variadic list of Scaled<Shifted<T>>'s
-static constexpr const auto scaled_of_shifted =
-    transform_variadic_list(shifted, al<ScaledEOS>{});
 // variadic list of UnitSystem<T>'s
 static constexpr const auto unit_system =
-    transform_variadic_list(partial_eos_list, al<UnitSystem>{});
+    transform_variadic_list(unit_system_eos_list, al<UnitSystem>{});
 // variadic list of Relativistic<T>'s
 static constexpr const auto relativistic =
-    transform_variadic_list(partial_eos_list, al<RelativisticEOS>{});
+    transform_variadic_list(relativistic_eos_list, al<RelativisticEOS>{});
+// variadic list of eos's with shifted or scaled modifiers
+static constexpr const auto shifted_1 =
+    transform_variadic_list(full_eos_list, al<ShiftedEOS>{});
+static constexpr const auto scaled_1 =
+    transform_variadic_list(full_eos_list, al<ScaledEOS>{});
 // relativistic and unit system modifiers
 static constexpr const auto unit_or_rel =
-    transform_variadic_list(partial_eos_list, apply_to_partial);
+    singularity::detail::concat(unit_system, relativistic);
+// variadic list of eos with shifted, relativistic or unit system modifiers
+static constexpr const auto shifted_of_unit_or_rel =
+    transform_variadic_list(unit_or_rel, al<ShiftedEOS>{});
+// combined list of all shifted EOS
+static constexpr const auto shifted =
+    singularity::detail::concat(shifted_1, shifted_of_unit_or_rel);
+// variadic list of eos with scaled, relativistic or unit system modifiers
+static constexpr const auto scaled_of_unit_or_rel =
+    transform_variadic_list(unit_or_rel, al<ScaledEOS>{});
+// variadic list of Scaled<Shifted<T>>'s
+static constexpr const auto scaled_of_shifted =
+    transform_variadic_list(shifted, al<ScaledEOS>{});
+// combined list of all scaled EOS
+static constexpr const auto scaled =
+    singularity::detail::concat(scaled_1, scaled_of_unit_or_rel, scaled_of_shifted);
 // create combined list
-static constexpr const auto combined_list_1 = singularity::detail::concat(
-    full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel);
+static constexpr const auto combined_list_1 =
+    singularity::detail::concat(full_eos_list, shifted, scaled, unit_or_rel);
 // make a ramped eos of everything
 static constexpr const auto ramped_all =
     transform_variadic_list(combined_list_1, al<BilinearRampEOS>{});

singularity-eos/eos/eos_eospac.hpp

@@ -656,8 +656,8 @@ class EOSPAC : public EosBase<EOSPAC> {
     EOS_REAL *R = const_cast<EOS_REAL *>(&rhos[0]);
     EOS_REAL *T = const_cast<EOS_REAL *>(&temperatures[0]);
     EOS_REAL *E = scratch + 0 * num;
-    EOS_REAL *DEDT = &cvs[0];
     EOS_REAL *DEDR = scratch + 1 * num;
+    EOS_REAL *DEDT = &cvs[0];
 
     EOS_INTEGER table = EofRT_table_;
 
@@ -682,20 +682,17 @@ class EOSPAC : public EosBase<EOSPAC> {
       }
     }
 
-    options[nopts] = EOS_F_CONVERT;
-    values[nopts] = cvFromSesame(1.0);
-
-    if (transform.f.is_set()) {
-      values[nopts] *= transform.f.get();
-    }
-    ++nopts;
-
     eosSafeInterpolate(&table, num, R, T, E, DEDR, DEDT, "EofRT", Verbosity::Quiet,
                        options, values, nopts);
 
+    const Real y = transform.y.is_set() ? (1.0 / transform.y.get()) : 1.0;
+    const Real f = transform.f.is_set() ? transform.f.get() : 1.0;
+
     portableFor(
         cname, 0, num, PORTABLE_LAMBDA(const int i) {
-          cvs[i] = std::max(cvs[i], 0.0); // Here we do something to the data!
+          cvs[i] =
+              f * y *
+              cvFromSesame(std::max(DEDT[i], 0.0)); // Here we do something to the data!
         });
   }
 

singularity-eos/eos/modifiers/eos_unitsystem.hpp

@@ -252,6 +252,153 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
     return inv_temp_unit_ * t_.MinimumTemperature();
   }
 
+  // vector implementations
+  template <typename LambdaIndexer>
+  inline void TemperatureFromDensityInternalEnergy(
+      const Real *rhos, const Real *sies, Real *temperatures, Real *scratch,
+      const int num, LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    transform.f.apply(inv_temp_unit_);
+    t_.TemperatureFromDensityInternalEnergy(rhos, sies, temperatures, scratch, num,
+                                            std::forward<LambdaIndexer>(lambdas),
+                                            std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void PressureFromDensityTemperature(const Real *rhos, const Real *temperatures,
+                                             Real *pressures, Real *scratch,
+                                             const int num, LambdaIndexer &&lambdas,
+                                             Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    transform.f.apply(inv_press_unit_);
+    t_.PressureFromDensityTemperature(rhos, temperatures, pressures, scratch, num,
+                                      std::forward<LambdaIndexer>(lambdas),
+                                      std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void
+  PressureFromDensityInternalEnergy(const Real *rhos, const Real *sies, Real *pressures,
+                                    Real *scratch, const int num, LambdaIndexer &&lambdas,
+                                    Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    transform.f.apply(inv_press_unit_);
+    t_.PressureFromDensityInternalEnergy(rhos, sies, pressures, scratch, num,
+                                         std::forward<LambdaIndexer>(lambdas),
+                                         std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void SpecificHeatFromDensityTemperature(
+      const Real *rhos, const Real *temperatures, Real *cvs, Real *scratch, const int num,
+      LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    transform.f.apply(inv_cv_unit_);
+    t_.SpecificHeatFromDensityTemperature(rhos, temperatures, cvs, scratch, num,
+                                          std::forward<LambdaIndexer>(lambdas),
+                                          std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void SpecificHeatFromDensityInternalEnergy(
+      const Real *rhos, const Real *sies, Real *cvs, Real *scratch, const int num,
+      LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    transform.f.apply(inv_cv_unit_);
+    t_.SpecificHeatFromDensityInternalEnergy(rhos, sies, cvs, scratch, num,
+                                             std::forward<LambdaIndexer>(lambdas),
+                                             std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void BulkModulusFromDensityTemperature(
+      const Real *rhos, const Real *temperatures, Real *bmods, Real *scratch,
+      const int num, LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    transform.f.apply(inv_bmod_unit_);
+    t_.BulkModulusFromDensityTemperature(rhos, temperatures, bmods, scratch, num,
+                                         std::forward<LambdaIndexer>(lambdas),
+                                         std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void BulkModulusFromDensityInternalEnergy(
+      const Real *rhos, const Real *sies, Real *bmods, Real *scratch, const int num,
+      LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    transform.f.apply(inv_bmod_unit_);
+    t_.BulkModulusFromDensityInternalEnergy(rhos, sies, bmods, scratch, num,
+                                            std::forward<LambdaIndexer>(lambdas),
+                                            std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void GruneisenParamFromDensityTemperature(
+      const Real *rhos, const Real *temperatures, Real *gm1s, Real *scratch,
+      const int num, LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    t_.GruneisenParamFromDensityTemperature(rhos, temperatures, gm1s, scratch, num,
+                                            std::forward<LambdaIndexer>(lambdas),
+                                            std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void GruneisenParamFromDensityInternalEnergy(
+      const Real *rhos, const Real *sies, Real *gm1s, Real *scratch, const int num,
+      LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    t_.GruneisenParamFromDensityInternalEnergy(rhos, sies, gm1s, scratch, num,
+                                               std::forward<LambdaIndexer>(lambdas),
+                                               std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void InternalEnergyFromDensityTemperature(
+      const Real *rhos, const Real *temperatures, Real *sies, Real *scratch,
+      const int num, LambdaIndexer &&lambdas, Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    transform.f.apply(inv_sie_unit_);
+    t_.InternalEnergyFromDensityTemperature(rhos, temperatures, sies, scratch, num,
+                                            std::forward<LambdaIndexer>(lambdas),
+                                            std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void EntropyFromDensityTemperature(const Real *rhos, const Real *temperatures,
+                                            Real *entropies, Real *scratch, const int num,
+                                            LambdaIndexer &&lambdas,
+                                            Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(temp_unit_);
+    transform.f.apply(inv_entropy_unit_);
+    t_.EntropyFromDensityTemperature(rhos, temperatures, entropies, scratch, num,
+                                     std::forward<LambdaIndexer>(lambdas),
+                                     std::forward<Transform>(transform));
+  }
+
+  template <typename LambdaIndexer>
+  inline void
+  EntropyFromDensityInternalEnergy(const Real *rhos, const Real *sies, Real *entropies,
+                                   Real *scratch, const int num, LambdaIndexer &&lambdas,
+                                   Transform &&transform = Transform()) const {
+    transform.x.apply(rho_unit_);
+    transform.y.apply(sie_unit_);
+    transform.f.apply(inv_entropy_unit_);
+    t_.EntropyFromDensityInternalEnergy(rhos, sies, entropies, scratch, num,
+                                        std::forward<LambdaIndexer>(lambdas),
+                                        std::forward<Transform>(transform));
+  }
+
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return t_.nlambda(); }
 

test/profile_eospac.cpp

@@ -367,6 +367,16 @@ int main(int argc, char *argv[]) {
       success = get_timing(ncycles, ncells_1d, X_MIN, X_MAX, Y_MIN, Y_MAX, "EOSPAC", eos);
     }
 
+    {
+      constexpr Real rho_unit = 1e1;
+      constexpr Real sie_unit = 1e-1;
+      constexpr Real temp_unit = 0.75;
+      auto eos = UnitSystem<EOSPAC>(EOSPAC(airID), eos_units_init::ThermalUnitsInit(),
+                                    rho_unit, sie_unit, temp_unit);
+      success = get_timing(ncycles, ncells_1d, X_MIN, X_MAX, Y_MIN, Y_MAX,
+                           "UnitSystem<EOSPAC>", eos);
+    }
+
     {
       auto eos = ShiftedEOS<EOSPAC>(EOSPAC(airID), 1.0);
       success = get_timing(ncycles, ncells_1d, X_MIN, X_MAX, Y_MIN, Y_MAX,