Merge pull request #278 from lanl/apg/sap_integration

Apg/sap integration
lanl · Nov 22, 2023 · b08627a · b08627a
2 parents a92cccc + 384be5a
commit b08627a
Show file tree

Hide file tree

Showing 28 changed files with 771 additions and 257 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,7 @@
 - [[PR269]](https://github.com/lanl/singularity-eos/pull/269) Add SAP Polynomial EoS
 
 ### Fixed (Repair bugs, etc)
+- [[PR278]](https://github.com/lanl/singularity-eos/pull/278) Fixed EOSPAC unit conversion errors for scalar lookups
 - [[PR316]](https://github.com/lanl/singularity-eos/pull/316) removed `fmax-errors=3` from `singularity-eos` compile flags
 - [[PR296]](https://github.com/lanl/singularity-eos/pull/296) changed `CMAKE_SOURCE_DIR` to `PROJECT_SOURCE_DIR` to fix downstream submodule build
 - [[PR291]](https://github.com/lanl/singularity-eos/pull/291) package.py updates to reflect new CMake options 
@@ -24,6 +25,9 @@
 - [[PR308]](https://github.com/lanl/singularity-eos/pull/308) spack builds +fortran now compile via correct blocking out of interfaces via preprocessor ifdef
 
 ### Added (new features/APIs/variables/...)
+- [[PR278]](https://github.com/lanl/singularity-eos/pull/278) Added EOSPAC option functionality in class constructor
+- [[PR278]](https://github.com/lanl/singularity-eos/pull/278) Added a new function for returning the minimum energy as a function of density for an EOS (only EOSPAC at the moment)
+- [[PR278]](https://github.com/lanl/singularity-eos/pull/278) Added a new Fortran API for simple pressure and bulk moduli lookups
 - [[PR306]](https://github.com/lanl/singularity-eos/pull/306) Added generic Evaluate method
 - [[PR304]](https://github.com/lanl/singularity-eos/pull/304) added a Newton-Raphson root find for use with the Helmholtz EoS
 - [[PR265]](https://github.com/lanl/singularity-eos/pull/265) Add missing UnitSystem modifier combinations to variant and EOSPAC

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -382,8 +382,10 @@ endforeach()
 
 target_sources(singularity-eos PRIVATE ${_srcs} ${_headers})
 
-# make sure .mods are placed in build path, and installed along with includes
 if(SINGULARITY_USE_FORTRAN)
+  # Turn on preprocessor for fortran files
+  set_target_properties(singularity-eos PROPERTIES Fortran_PREPROCESS ON)
+  # make sure .mods are placed in build path, and installed along with includes
   set_target_properties(singularity-eos PROPERTIES Fortran_MODULE_DIRECTORY
                                                    ${CMAKE_CURRENT_BINARY_DIR}/fortran)
   target_include_directories(

diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst
@@ -1423,11 +1423,23 @@ This is a striaghtforward wrapper of the `EOSPAC`_ library for the
 
 .. code-block::
 
-  EOSPAC(int matid, bool invert_at_setup = false)
+  EOSPAC(int matid, bool invert_at_setup = false, Real insert_data = 0.0, eospacMonotonicity monotonicity = eospacMonotonicity::none, bool apply_smoothing = false, eospacSplit apply_splitting = eospacSplit::none, bool linear_interp = false)
 
-where ``matid`` is the unique material number in the database and
+where ``matid`` is the unique material number in the database,
 ``invert_at_setup`` specifies whether or not pre-compute tables of
-temperature as a function of density and energy.
+temperature as a function of density and energy, ``insert_data`` 
+inserts specified number of grid points between original grid points 
+in the `Sesame`_ table, ``monotonicity` enforces monotonicity in x, 
+y or both (:math:`monotonicityX/Y/XY`), ``apply_smoothing`` enables 
+data table smoothing that imposes a linear floor on temperature dependence, 
+forces linear temperature dependence for low temperature, and forces 
+linear density dependence for low and high density, ``apply_splitting`` 
+has the following options for ion data tables not found in the `Sesame`_ 
+database :. :math:`splitNumProp` uses the cold curve plus number-proportional 
+model, :math:`splitIdealGas` uses the cold curve plus ideal gas model 
+and :math:`splitCowan` uses the cold curve plus Cowan-nuclear model 
+for ions and the final option ``linear_interp`` uses linear instead of 
+bilinear interpolation. 
 
 Note for performance reasons this EOS uses a slightly different vector API.
 See :ref:`EOSPAC Vector Functions <eospac_vector>` for more details.

diff --git a/eospac-wrapper/eospac_wrapper.cpp b/eospac-wrapper/eospac_wrapper.cpp
@@ -14,15 +14,14 @@
 // publicly and display publicly, and to permit others to do so.
 //======================================================================
 
+#include "eospac_wrapper.hpp"
 #include <array>
+#include <eos_Interface.h>
 #include <iostream>
 #include <regex>
 #include <string>
 #include <vector>
 
-#include "eospac_wrapper.hpp"
-#include <eos_Interface.h>
-
 namespace EospacWrapper {
 
 void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
@@ -87,7 +86,9 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
 
     comments.resize(static_cast<int>(commentLen));
     metadata.comments.resize(comments.size());
-    eosSafeTableCmnts(&eospacComments, comments.data(), eospacWarn);
+
+    if (comments.size() > 0)
+      eosSafeTableCmnts(&eospacComments, comments.data(), eospacWarn);
     for (size_t i = 0; i < comments.size(); i++) {
       metadata.comments[i] = comments[i];
     }
@@ -107,16 +108,21 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
 
 EOS_INTEGER eosSafeLoad(int ntables, int matid, EOS_INTEGER tableType[],
                         EOS_INTEGER tableHandle[], Verbosity eospacWarn,
-                        bool invert_at_setup) {
+                        bool invert_at_setup, double insert_data,
+                        eospacMonotonicity monotonicity, bool apply_smoothing,
+                        eospacSplit apply_splitting, bool linear_interp) {
   std::vector<std::string> empty;
   return eosSafeLoad(ntables, matid, tableType, tableHandle, empty, eospacWarn,
-                     invert_at_setup);
+                     invert_at_setup, insert_data, monotonicity, apply_smoothing,
+                     apply_splitting, linear_interp);
 }
 
 EOS_INTEGER eosSafeLoad(int ntables, int matid, EOS_INTEGER tableType[],
                         EOS_INTEGER tableHandle[],
                         const std::vector<std::string> &table_names, Verbosity eospacWarn,
-                        bool invert_at_setup) {
+                        bool invert_at_setup, double insert_data,
+                        eospacMonotonicity monotonicity, bool apply_smoothing,
+                        eospacSplit apply_splitting, bool linear_interp) {
   EOS_INTEGER NTABLES[] = {ntables};
   std::vector<EOS_INTEGER> MATID(ntables, matid);
 
@@ -127,19 +133,83 @@ EOS_INTEGER eosSafeLoad(int ntables, int matid, EOS_INTEGER tableType[],
   eos_CreateTables(NTABLES, tableType, MATID.data(), tableHandle, &errorCode);
 
   if (invert_at_setup) {
-    EOS_INTEGER options[] = {EOS_INVERT_AT_SETUP, EOS_INSERT_DATA};
-    EOS_REAL values[] = {1., 4.};
+    EOS_REAL values[] = {1.};
     for (int i = 0; i < ntables; i++) {
-      if (tableType[i] == EOS_T_DUt) {
-        eos_SetOption(&(tableHandle[i]), &(options[0]), &(values[0]), &errorCode);
-        eos_SetOption(&(tableHandle[i]), &(options[1]), &(values[1]), &errorCode);
+      if (tableType[i] != EOS_Uc_D and tableType[i] != EOS_Pc_D and
+          tableType[i] != EOS_Pt_DT and tableType[i] != EOS_Ut_DT) {
+        eos_SetOption(&(tableHandle[i]), &EOS_INVERT_AT_SETUP, &(values[0]), &errorCode);
+        eosCheckError(errorCode, "eospac options: eos_invert_at_setup", eospacWarn);
       }
     }
   }
 
+  if (insert_data > 0) {
+    EOS_REAL values[] = {insert_data};
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_INSERT_DATA, &(values[0]), &errorCode);
+      eosCheckError(errorCode, "eospac options: eos_insert_data", eospacWarn);
+    }
+  }
+
+  // choice of which table types setOption is called on mimics SAP. Some table types are
+  // incmopatible whiel others lead to numerical issues.
+  if (monotonicity == eospacMonotonicity::monotonicityX ||
+      monotonicity == eospacMonotonicity::monotonicityXY) {
+    for (int i = 0; i < ntables; i++) {
+      if (tableType[i] != EOS_Uc_D and tableType[i] != EOS_Ut_DT and
+          tableType[i] != EOS_Ut_DPt) {
+        eos_SetOption(&(tableHandle[i]), &EOS_MONOTONIC_IN_X, &(EOS_NullVal), &errorCode);
+        eosCheckError(errorCode, "eospac options: eos_monotonic_in_x", eospacWarn);
+      }
+    }
+  }
+
+  if (monotonicity == eospacMonotonicity::monotonicityXY ||
+      monotonicity == eospacMonotonicity::monotonicityY) {
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_MONOTONIC_IN_Y, &(EOS_NullVal), &errorCode);
+      eosCheckError(errorCode, "eospac options: eos_monotonic_in_y", eospacWarn);
+    }
+  }
+
+  if (apply_smoothing) {
+    for (int i = 0; i < ntables; i++) {
+      if (tableType[i] == EOS_Pt_DUt || tableType[i] == EOS_T_DUt ||
+          tableType[i] == EOS_Ut_DPt) {
+        eos_SetOption(&(tableHandle[i]), &EOS_SMOOTH, &(EOS_NullVal), &errorCode);
+        eosCheckError(errorCode, "eospac options: eos_smooth", eospacWarn);
+      }
+    }
+  }
+
+  if (apply_splitting == eospacSplit::splitNumProp) {
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_SPLIT_NUM_PROP, &(EOS_NullVal), &errorCode);
+      eosCheckError(errorCode, "eospac options: eos_split_num_prop", eospacWarn);
+    }
+  } else if (apply_splitting == eospacSplit::splitIdealGas) {
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_SPLIT_IDEAL_GAS, &(EOS_NullVal), &errorCode);
+      eosCheckError(errorCode, "eospac options: eos_split_ideal_gas", eospacWarn);
+    }
+  } else if (apply_splitting == eospacSplit::splitCowan) {
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_SPLIT_COWAN, &(EOS_NullVal), &errorCode);
+      eosCheckError(errorCode, "eospac options: eos_split_cowan", eospacWarn);
+    }
+  }
+
+  if (linear_interp) {
+    for (int i = 0; i < ntables; i++) {
+      eos_SetOption(&(tableHandle[i]), &EOS_LINEAR, &(EOS_NullVal), &errorCode);
+      eosCheckError(errorCode, "eospac options: linear_interp", eospacWarn);
+    }
+  }
+
 #ifdef SINGULARITY_EOSPAC_SKIP_EXTRAP
   for (int i = 0; i < ntables; i++) {
-    eos_SetOption(&(tableHandle[i]), &EOS_SKIP_EXTRAP_CHECK, NULL, &errorCode);
+    eos_SetOption(&(tableHandle[i]), &EOS_SKIP_EXTRAP_CHECK, &(EOS_NullVal), &errorCode);
+    eosCheckError(errorCode, "eos_skip_extrap_check", eospacWarn);
   }
 #endif // SINGULARITY_EOSPAC_SKIP_EXTRAP
 
@@ -181,6 +251,7 @@ bool eosSafeInterpolate(EOS_INTEGER *table, EOS_INTEGER nxypairs, EOS_REAL xVals
   }
 
   eos_Interpolate(table, &nxypairs, xVals, yVals, var, dx, dy, &errorCode);
+
 #ifndef SINGULARITY_EOSPAC_SKIP_EXTRAP
   if (errorCode != EOS_OK && eospacWarn == Verbosity::Debug) {
     eos_GetErrorMessage(&errorCode, errorMessage);

diff --git a/eospac-wrapper/eospac_wrapper.hpp b/eospac-wrapper/eospac_wrapper.hpp
@@ -19,11 +19,18 @@
 
 #include <iostream>
 #include <string>
+#include <vector>
 
 #include <eos_Interface.h> // eospac API
 
 namespace EospacWrapper {
-
+enum class eospacSplit { none = 0, splitNumProp = 1, splitIdealGas = 2, splitCowan = 3 };
+enum class eospacMonotonicity {
+  none = 0,
+  monotonicityX = 1,
+  monotonicityXY = 2,
+  monotonicityY = 3
+};
 inline double densityToSesame(const double CodeTemp) { return CodeTemp; }
 inline double densityFromSesame(const double sesTemp) { return sesTemp; }
 inline double temperatureToSesame(const double CodeTemp) { return CodeTemp; }
@@ -87,11 +94,19 @@ void eosGetMetadata(int matid, SesameMetadata &metadata,
 
 EOS_INTEGER eosSafeLoad(int ntables, int matid, EOS_INTEGER tableType[],
                         EOS_INTEGER tableHandle[], Verbosity eospacWarn,
-                        bool invert_at_setup = false);
+                        bool invert_at_setup = false, EOS_REAL insert_data = 0.0,
+                        eospacMonotonicity monotonicity = eospacMonotonicity::none,
+                        bool apply_smoothing = false,
+                        eospacSplit apply_splitting = eospacSplit::none,
+                        bool linear_interp = false);
 EOS_INTEGER eosSafeLoad(int ntables, int matid, EOS_INTEGER tableType[],
                         EOS_INTEGER tableHandle[],
                         const std::vector<std::string> &table_names, Verbosity eospacWarn,
-                        bool invert_at_setup = false);
+                        bool invert_at_setup = false, EOS_REAL insert_data = 0.0,
+                        eospacMonotonicity monotonicity = eospacMonotonicity::none,
+                        bool apply_smoothing = false,
+                        eospacSplit apply_splitting = eospacSplit::none,
+                        bool linear_interp = false);
 
 // output is boolean mask. 1 for no errors. 0 for errors.
 bool eosSafeInterpolate(EOS_INTEGER *table, EOS_INTEGER nxypairs, EOS_REAL xVals[],

diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
@@ -35,7 +35,7 @@ namespace singularity {
 namespace eos_base {
 
 namespace impl {
-constexpr std::size_t MAX_NUM_CHARS = 81;
+constexpr std::size_t MAX_NUM_CHARS = 121;
 // Cuda doesn't have strcat, so we implement it ourselves
 PORTABLE_FORCEINLINE_FUNCTION
 char *StrCat(char *destination, const char *source) {
@@ -68,6 +68,7 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::InternalEnergyFromDensityTemperature;                       \
   using EosBase<EOSDERIVED>::PressureFromDensityTemperature;                             \
   using EosBase<EOSDERIVED>::PressureFromDensityInternalEnergy;                          \
+  using EosBase<EOSDERIVED>::MinInternalEnergyFromDensity;                               \
   using EosBase<EOSDERIVED>::SpecificHeatFromDensityTemperature;                         \
   using EosBase<EOSDERIVED>::SpecificHeatFromDensityInternalEnergy;                      \
   using EosBase<EOSDERIVED>::BulkModulusFromDensityTemperature;                          \
@@ -80,6 +81,7 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::EntropyFromDensityTemperature;                              \
   using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;                           \
   using EosBase<EOSDERIVED>::EntropyIsNotEnabled;                                        \
+  using EosBase<EOSDERIVED>::MinInternalEnergyIsNotEnabled;                              \
   using EosBase<EOSDERIVED>::IsModified;                                                 \
   using EosBase<EOSDERIVED>::UnmodifyOnce;                                               \
   using EosBase<EOSDERIVED>::GetUnmodifiedObject;
@@ -271,6 +273,35 @@ class EosBase {
     PressureFromDensityInternalEnergy(rhos, sies, pressures, num,
                                       std::forward<LambdaIndexer>(lambdas));
   }
+  ///
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void MinInternalEnergyFromDensity(ConstRealIndexer &&rhos, RealIndexer &&sies,
+                                           const int num, LambdaIndexer &&lambdas) const {
+    static auto const name = SG_MEMBER_FUNC_NAME();
+    static auto const cname = name.c_str();
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) {
+          sies[i] = copy.MinInternalEnergyFromDensity(rhos[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer,
+            typename = std::enable_if_t<!is_raw_pointer<RealIndexer, Real>::value>>
+  inline void MinInternalEnergyFromDensity(ConstRealIndexer &&rhos, RealIndexer &&sies,
+                                           Real * /*scratch*/, const int num,
+                                           LambdaIndexer &&lambdas) const {
+    MinInternalEnergyFromDensity(std::forward<ConstRealIndexer>(rhos),
+                                 std::forward<RealIndexer>(sies), num,
+                                 std::forward<LambdaIndexer>(lambdas));
+  }
+  template <typename LambdaIndexer>
+  inline void MinInternalEnergyFromDensity(const Real *rhos, Real *sies,
+                                           Real * /*scratch*/, const int num,
+                                           LambdaIndexer &&lambdas,
+                                           Transform && = Transform()) const {
+    MinInternalEnergyFromDensity(rhos, num, std::forward<LambdaIndexer>(lambdas));
+  }
+  ///
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos,
                                             ConstRealIndexer &&temperatures,
@@ -563,15 +594,30 @@ class EosBase {
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const { return 0.0; }
 
-  // Default entropy behavior is to return an error
+  // Default entropy behavior is to cause an error
   PORTABLE_FORCEINLINE_FUNCTION
   void EntropyIsNotEnabled(const char *eosname) const {
     // Construct the error message using char* so it works on device
     // WARNING: This needs to be updated if EOS names get longer
     // base msg length 32 + 5 chars = 37 chars
     // + 1 char for null terminator
     // maximum allowed EOS length = 44 chars
-    char msg[impl::MAX_NUM_CHARS] = "Entropy is not enabled for the '";
+    char msg[impl::MAX_NUM_CHARS] = "Singularity-EOS: Entropy is not enabled for the '";
+    impl::StrCat(msg, eosname);
+    impl::StrCat(msg, "' EOS");
+    PORTABLE_ALWAYS_THROW_OR_ABORT(msg);
+  }
+
+  // Default MinInternalEnergyFromDensity behavior is to cause an error
+  PORTABLE_FORCEINLINE_FUNCTION
+  void MinInternalEnergyIsNotEnabled(const char *eosname) const {
+    // Construct the error message using char* so it works on device
+    // WARNING: This needs to be updated if EOS names get longer
+    // base msg length 32 + 5 chars = 37 chars
+    // + 1 char for null terminator
+    // maximum allowed EOS length = 44 chars
+    char msg[impl::MAX_NUM_CHARS] =
+        "Singularity-EOS: MinInternalEnergyFromDensity() is not enabled for the '";
     impl::StrCat(msg, eosname);
     impl::StrCat(msg, "' EOS");
     PORTABLE_ALWAYS_THROW_OR_ABORT(msg);

diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp
@@ -60,6 +60,12 @@ class DavisReactants : public EosBase<DavisReactants> {
       const Real rho, const Real sie, Real *lambda = nullptr) const {
     return Ps(rho) + Gamma(rho) * rho * (sie - Es(rho));
   }
+
+  PORTABLE_INLINE_FUNCTION Real
+  MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const {
+    MinInternalEnergyIsNotEnabled("DavisReactants");
+    return 0.0;
+  }
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const {
     EntropyIsNotEnabled("DavisReactants");
@@ -162,6 +168,11 @@ class DavisProducts : public EosBase<DavisProducts> {
       const Real rho, const Real sie, Real *lambda = nullptr) const {
     return Ps(rho) + rho * Gamma(rho) * (sie - Es(rho));
   }
+  PORTABLE_INLINE_FUNCTION Real
+  MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const {
+    MinInternalEnergyIsNotEnabled("DavisProducts");
+    return 0.0;
+  }
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const {
     EntropyIsNotEnabled("DavisProducts");