Merge pull request #101 from lanl/dholladay00/sap_ramp

Dholladay00/sap ramp

.github/workflows/formatting.yml

@@ -18,6 +18,6 @@ jobs:
           run: export DEBIAN_FRONTEND=noninteractive
         - name: install dependencies
           run: |
-            sudo apt-get install -y --force-yes -qq git clang-format
+            sudo apt-get install -y --force-yes -qq git clang-format-12
         - name: check formatting
-          run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format -style=file -i && git diff --exit-code --ignore-submodules
+          run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format-12 -style=file -i && git diff --exit-code --ignore-submodules

CMakeLists.txt

@@ -211,7 +211,7 @@ if(SINGULARITY_USE_KOKKOS)
   ####
   # TODO: Once it's oppurtune, we will
   #       not support Kokkos < 3.3
-  if(Kokkos_VERSION_MAJOR VERSION_LESS "3.3")
+  if(${Kokkos_VERSION_MAJOR}.${Kokkos_VERSION_MINOR} VERSION_LESS "3.3")
     message(WARNING "`Kokkos` version [${Kokkos_VERSION_MAJOR}.${Kokkos_VERSION_MINOR}] is DEPRECATED, and `singularity-eos` will not support versions < '3.3' in the very near future.")
     if(SINGULARITY_USE_CUDA)
       get_filename_component(_compiler_exe ${CMAKE_CXX_COMPILER} NAME)

doc/sphinx/index.rst

@@ -19,6 +19,7 @@ Documentation approved for unlimited release. LA-UR-21-31131.
    src/building
    src/using-eos
    src/models
+   src/modifiers
    src/using-closures
    src/contributing
    src/sphinx-doc

doc/sphinx/src/modifiers.rst

@@ -0,0 +1,9 @@
+.. _modifiers:
+
+EOS Modifiers
+==============
+
+An  equation of  state  modifier (perhaps  intuitively) *modifies*  an
+equation of state in some pre-prescribed way. For example, one can use
+a modifier to  shift the assumed zero-point energy of  the equation of
+state, add a unit system, or account for porosity.

example/CMakeLists.txt

@@ -15,13 +15,11 @@
 add_executable(get_sound_speed_press
   get_sound_speed_press.cpp)
 target_link_libraries(get_sound_speed_press PRIVATE
-  eos
-  singularity-eos::libs
-  singularity-eos::flags)
+  singularity-eos::singularity-eos)
 
-add_executable(get_sesame_state
-  get_sesame_state.cpp)
-target_link_libraries(get_sesame_state PRIVATE
-  eos
-  singularity-eos::libs
-  singularity-eos::flags)
+if(SINGULARITY_USE_EOSPAC)
+  add_executable(get_sesame_state
+    get_sesame_state.cpp)
+  target_link_libraries(get_sesame_state PRIVATE
+    singularity-eos::singularity-eos)
+endif()

singularity-eos/CMakeLists.txt

@@ -16,11 +16,13 @@ set(EOS_HEADERS
     base/constants.hpp
     base/eos_error.hpp
     base/fast-math/logs.hpp
+    base/robust_utils.hpp
     base/root-finding-1d/root_finding.hpp
     eos/eos.hpp
     eos/eos_builder.hpp
     eos/eos_variant.hpp
     eos/modifiers/eos_unitsystem.hpp
+    eos/modifiers/ramps_eos.hpp
     eos/modifiers/scaled_eos.hpp
     eos/modifiers/shifted_eos.hpp
     eos/modifiers/relativistic_eos.hpp

singularity-eos/base/robust_utils.hpp

@@ -0,0 +1,57 @@
+//------------------------------------------------------------------------------
+// © 2021-2022. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_
+#define SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_
+
+#include <limits>
+#include <ports-of-call/portability.hpp>
+
+namespace singularity {
+namespace robust {
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() {
+  return 10 * std::numeric_limits<T>::min();
+}
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() {
+  return 10 * std::numeric_limits<T>::epsilon();
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) {
+  return std::max(val, EPS<T>());
+}
+
+PORTABLE_FORCEINLINE_FUNCTION
+Real make_bounded(const Real val, const Real vmin, const Real vmax) {
+  return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS()));
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) {
+  return (T(0) <= val) - (val < T(0));
+}
+
+template <typename A, typename B>
+PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) {
+  return a / (b + sgn(b) * SMALL<B>());
+}
+
+} // namespace robust
+} // namespace singularity
+
+#endif // SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_

singularity-eos/eos/eos.hpp

@@ -40,6 +40,7 @@
 #include <singularity-eos/base/eos_error.hpp>
 #include <singularity-eos/base/variadic_utils.hpp>
 #include <singularity-eos/eos/modifiers/eos_unitsystem.hpp>
+#include <singularity-eos/eos/modifiers/ramps_eos.hpp>
 #include <singularity-eos/eos/modifiers/relativistic_eos.hpp>
 #include <singularity-eos/eos/modifiers/scaled_eos.hpp>
 #include <singularity-eos/eos/modifiers/shifted_eos.hpp>
@@ -528,9 +529,12 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); }
+  PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return _n_lambda; }
-  inline RootFinding1D::Status rootStatus() const { return status_; }
-  inline TableStatus tableStatus() const { return whereAmI_; }
+  PORTABLE_INLINE_FUNCTION
+  RootFinding1D::Status rootStatus() const { return status_; }
+  PORTABLE_INLINE_FUNCTION
+  TableStatus tableStatus() const { return whereAmI_; }
   RootFinding1D::RootCounts counts;
   void Finalize();
   static std::string EosType() { return std::string("SpinerEOSDependsRhoT"); }
@@ -744,6 +748,7 @@ class SpinerEOSDependsRhoSie : public EosBase<SpinerEOSDependsRhoSie> {
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); }
 
+  PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return _n_lambda; }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     static constexpr char s1[]{"SpinerEOS Parameters:"};
@@ -752,6 +757,7 @@ class SpinerEOSDependsRhoSie : public EosBase<SpinerEOSDependsRhoSie> {
     printf("%s\n\t%s\n\t%s%i\n", s1, s2, s3, matid_);
     return;
   }
+  PORTABLE_INLINE_FUNCTION
   RootFinding1D::Status rootStatus() const { return status_; }
   RootFinding1D::RootCounts counts;
   static std::string EosType() { return std::string("SpinerEOSDependsRhoSie"); }
@@ -1296,8 +1302,14 @@ static constexpr const auto scaled_of_shifted =
 static constexpr const auto unit_or_rel =
     transform_variadic_list(partial_eos_list, apply_to_partial);
 // create combined list
-static constexpr const auto combined_list = singularity::detail::concat(
+static constexpr const auto combined_list_1 = singularity::detail::concat(
     full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel);
+// make a ramped eos of everything
+static constexpr const auto ramped_all =
+    transform_variadic_list(combined_list_1, al<BilinearRampEOS>{});
+// final combined list
+static constexpr const auto combined_list = singularity::detail::concat(
+    full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel, ramped_all);
 // a function that returns a Variant from a typelist
 template <typename... Ts>
 struct tl_to_Variant_struct {

singularity-eos/eos/eos_builder.cpp

@@ -27,7 +27,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
   bool shifted = (modifiers.count(EOSModifier::Shifted) > 0);
   bool relativistic = (modifiers.count(EOSModifier::Relativistic) > 0);
   bool units = (modifiers.count(EOSModifier::UnitSystem) > 0);
-  if ((shifted || scaled || relativistic || units) && !(isModifiable(type))) {
+  bool ramped = (modifiers.count(EOSModifier::BilinearRamp) > 0);
+  if ((shifted || scaled || relativistic || units || ramped) && !(isModifiable(type))) {
     EOS_ERROR("Modifiers not supported for this EOS");
   }
   if (relativistic && (shifted || scaled || units)) {
@@ -53,6 +54,16 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
     use_length_time =
         mpark::get<bool>(modifiers[EOSModifier::UnitSystem]["use_length_time"]);
   }
+  Real r0 = 0;
+  Real a = 1;
+  Real b = 0;
+  Real c = 0;
+  if (ramped) {
+    r0 = mpark::get<Real>(modifiers[EOSModifier::BilinearRamp]["r0"]);
+    a = mpark::get<Real>(modifiers[EOSModifier::BilinearRamp]["a"]);
+    b = mpark::get<Real>(modifiers[EOSModifier::BilinearRamp]["b"]);
+    c = mpark::get<Real>(modifiers[EOSModifier::BilinearRamp]["c"]);
+  }
   Real rho_unit = 1;
   Real sie_unit = 1;
   Real temp_unit = 1;
@@ -93,7 +104,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
     if (relativistic) {
       return makeRelativistic(std::move(g), cl);
     }
-    return applyShiftAndScale(std::move(g), scaled, shifted, scale, shift);
+    return applyShiftAndScaleAndBilinearRamp(std::move(g), scaled, shifted, ramped, scale,
+                                             shift, r0, a, b, c);
   }
 #ifdef SPINER_USE_HDF
   if (type == EOSType::SpinerEOSDependsRhoT || type == EOSType::SpinerEOSDependsRhoSie) {
@@ -120,7 +132,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
         if (relativistic) {
           return makeRelativistic(std::move(s), cl);
         }
-        return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift);
+        return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped,
+                                                 scale, shift, r0, a, b, c);
       }
     } else {
       string materialName = mpark::get<string>(base_params["materialName"]);
@@ -143,7 +156,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
         if (relativistic) {
           return makeRelativistic(std::move(s), cl);
         }
-        return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift);
+        return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped,
+                                                 scale, shift, r0, a, b, c);
       }
     }
   }
@@ -165,7 +179,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
     if (relativistic) {
       return makeRelativistic(std::move(s), cl);
     }
-    return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift);
+    return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale,
+                                             shift, r0, a, b, c);
   }
 #endif
   if (type == EOSType::Gruneisen) {

singularity-eos/eos/eos_builder.hpp

@@ -43,14 +43,15 @@ enum class EOSType {
   StellarCollapse
 #endif
 };
-enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem };
+enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem, BilinearRamp };
 
 // evil type erasure
 using param_t = mpark::variant<bool, int, Real, std::string>;
 using params_t = std::map<std::string, param_t>;
 using modifiers_t = std::map<EOSModifier, params_t>;
 const params_t NO_PARAMS;
 
+// TODO: Extend if needed
 const std::unordered_set<EOSType> modifiable({EOSType::IdealGas
 #ifdef SPINER_USE_HDF
                                               ,
@@ -83,6 +84,11 @@ EOS makeRelativistic(T &&eos, Real cl) {
   return RelativisticEOS<T>(std::move(eos), cl);
 }
 
+template <typename T>
+EOS makeBilinearRamp(T &&eos, Real r0, Real a, Real b, Real c) {
+  return BilinearRampEOS<T>(std::forward<T>(eos), r0, a, b, c);
+}
+
 template <typename T>
 EOS applyShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shift) {
   if (shifted && scaled) {
@@ -98,6 +104,38 @@ EOS applyShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shif
   }
   return eos;
 }
+
+template <typename T, template <typename> class W, typename... ARGS>
+EOS applyWrappedShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shift,
+                              ARGS... args) {
+  if (shifted && scaled) {
+    ShiftedEOS<T> a(std::forward<T>(eos), shift);
+    ScaledEOS<ShiftedEOS<T>> b(std::move(a), scale);
+    W<ScaledEOS<ShiftedEOS<T>>> c(std::move(b), args...);
+    return c;
+  }
+  if (shifted) {
+    ShiftedEOS<T> sh_eos(std::forward<T>(eos), shift);
+    return W<ShiftedEOS<T>>(std::move(sh_eos), args...);
+  }
+  if (scaled) {
+    ScaledEOS<T> sc_eos(std::forward<T>(eos), scale);
+    return W<ScaledEOS<T>>(std::move(sc_eos), args...);
+  }
+  return W<T>(std::forward<T>(eos), args...);
+}
+
+template <typename T>
+EOS applyShiftAndScaleAndBilinearRamp(T &&eos, bool scaled, bool shifted, bool ramped,
+                                      Real scale, Real shift, Real r0, Real a, Real b,
+                                      Real c) {
+  if (ramped) {
+    return applyWrappedShiftAndScale<T, BilinearRampEOS>(
+        std::forward<T>(eos), scaled, shifted, scale, shift, r0, a, b, c);
+  } else {
+    return applyShiftAndScale(std::forward<T>(eos), scaled, shifted, scale, shift);
+  }
+}
 } // namespace EOSBuilder
 
 } // namespace singularity