Apg/sap integration

eospac-wrapper/eospac_wrapper.cpp

@@ -19,7 +19,6 @@
 #include <regex>
 #include <string>
 #include <vector>
-
 #include "eospac_wrapper.hpp"
 #include <eos_Interface.h>
 
@@ -87,7 +86,8 @@ 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 (commentLen>0)    eosSafeTableCmnts(&eospacComments, comments.data(), eospacWarn);
     for (size_t i = 0; i < comments.size(); i++) {
       metadata.comments[i] = comments[i];
     }
@@ -107,16 +107,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 +132,81 @@ 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_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_T_DUt) {
-        eos_SetOption(&(tableHandle[i]), &(options[0]), &(values[0]), &errorCode);
-        eos_SetOption(&(tableHandle[i]), &(options[1]), &(values[1]), &errorCode);
+      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 +248,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);

eospac-wrapper/eospac_wrapper.hpp

@@ -22,8 +22,10 @@
 
 #include <eos_Interface.h> // eospac API
 
-namespace EospacWrapper {
 
+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 +89,15 @@ 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[],

singularity-eos/eos/eos_base.hpp

@@ -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;                          \
@@ -262,6 +263,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,

singularity-eos/eos/eos_davis.hpp

@@ -60,6 +60,15 @@ 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 {
+#ifndef NDEBUG
+      printf(
+          "WARNING: MinInternalEnergtyFromDensity is not defined for Davis EOS.");
+#endif
+    return 0.0;
+  }
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const {
     EntropyIsNotEnabled("DavisReactants");
@@ -162,6 +171,10 @@ 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 {
+    return 0.0;
+  }
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const {
     EntropyIsNotEnabled("DavisProducts");