Add mass fraction multi-phase table to sesame2spiner

eospac-wrapper/eospac_wrapper.cpp

@@ -335,6 +335,14 @@ void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacW
   eosCheckError(errorCode, "eos_GetTableCmnts", eospacWarn);
 }
 
+void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString,
+                          Verbosity eospacWarn) {
+  EOS_INTEGER errorCode = EOS_OK;
+  EOS_INTEGER infoTypes[1] = {infoItem};
+  eos_GetTableMetaData(table, infoTypes, infoString, &errorCode);
+  eosCheckError(errorCode, "eos_GetTableMetaData", eospacWarn);
+}
+
 void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn) {
   EOS_CHAR errorMessage[EOS_MaxErrMsgLen];
   if (errorCode != EOS_OK && eospacWarn != Verbosity::Quiet) {
@@ -343,6 +351,31 @@ void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eos
   }
 }
 
+int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn) {
+  // Note this opens and closes the table
+  int exists = 1;
+  int NT = 1;
+  EOS_INTEGER matid[1] = {matid_};
+  EOS_INTEGER tableType[1] = {tableType_};
+  EOS_INTEGER tableHandle[1];
+  EOS_INTEGER errorCode = EOS_OK;
+
+  eos_CreateTables(&NT, tableType, matid, tableHandle, &errorCode);
+  eosCheckError(errorCode, "eos_CreateTables", eospacWarn);
+  eos_LoadTables(&NT, &tableHandle[0], &errorCode);
+  if (errorCode != EOS_OK) {
+    if (eosErrorCodesEqual(EOS_NO_DATA_TABLE, errorCode)) {
+      exists = 0;
+    } else {
+      // something else happended
+      eosCheckError(errorCode, "eos_LoadTables", eospacWarn);
+      exists = -1;
+    }
+  }
+  eosSafeDestroy(NT, tableHandle, eospacWarn);
+  return exists;
+}
+
 std::string eosErrorString(EOS_INTEGER errorCode) {
   // Algorithmicallly generated by parsing the EOSPAC docs
   // I'm sorry. It's gross. ~JMM

eospac-wrapper/eospac_wrapper.hpp

@@ -122,13 +122,14 @@ void eosSafeTableInfo(EOS_INTEGER *table, EOS_INTEGER numInfoItems,
                       EOS_INTEGER infoItems[], EOS_REAL infoVals[], Verbosity eospacWarn);
 
 void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacWarn);
-
+void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString,
+                          Verbosity eospacWarn);
 void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn);
 std::string eosErrorString(EOS_INTEGER errorCode);
 void eosSafeDestroy(int ntables, EOS_INTEGER tableHandles[], Verbosity eospacWarn);
 std::string getName(std::string comment);
 bool eosErrorCodesEqual(EOS_INTEGER lhs, EOS_INTEGER rhs);
-
+int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn);
 } // namespace EospacWrapper
 
 #endif // _EOSPAC_WRAPPER_EOSPAC_WRAPPER_HPP_

sesame2spiner/generate_files.cpp

@@ -53,7 +53,7 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh
   std::string sMatid = std::to_string(matid);
 
   herr_t status = 0;
-  hid_t matGroup, lTGroup, leGroup, coldGroup;
+  hid_t matGroup, lTGroup, leGroup, coldGroup, mfGroup;
 
   matGroup = H5Gcreate(loc, sMatid.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   status += H5Lcreate_soft(sMatid.c_str(), loc, name.c_str(), H5P_DEFAULT, H5P_DEFAULT);
@@ -116,6 +116,25 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh
     // status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask);
   }
 
+  {
+    DataBox mf, mask;
+    Bounds nphBounds;
+    std::string phase_names;
+
+    if (eosMassFraction(matid, lRhoBounds, lTBounds, nphBounds, mf, mask, phase_names,
+                        eospacWarn)) {
+      mfGroup = H5Gcreate(matGroup, SP5::Depends::massFrac, H5P_DEFAULT, H5P_DEFAULT,
+                          H5P_DEFAULT);
+
+      status += mf.saveHDF(mfGroup, SP5::Fields::massFrac);
+      const int numphases = mf.dim(1);
+      status += H5LTset_attribute_int(mfGroup, ".", "numphases", &numphases, 1);
+      status +=
+          H5LTset_attribute_string(mfGroup, ".", "phase names", phase_names.c_str());
+      status += H5Gclose(mfGroup);
+    }
+  }
+
   if (addSubtables) {
     int i = 0;
     std::vector<TableSplit> splits = {TableSplit::ElectronOnly, TableSplit::IonCold};

sesame2spiner/io_eospac.cpp

@@ -328,6 +328,80 @@ void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie,
   eosSafeDestroy(NT, tableHandle, eospacWarn);
 }
 
+bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds,
+                     Bounds &nphBounds, DataBox &Ms, DataBox &mask,
+                     std::string &phase_names, Verbosity eospacWarn) {
+  using namespace EospacWrapper;
+
+  EOS_INTEGER errorCode = EOS_OK;
+  constexpr int NT = 2;
+  int ntables = NT;
+  EOS_INTEGER tableHandle[NT];
+  EOS_INTEGER tableType[NT] = {EOS_M_DT, EOS_Comment};
+  std::vector<EOS_INTEGER> matid_v(NT, matid);
+  std::vector<std::string> table_names = {"EOS_M_DT", "Eos_Material_Phases"};
+
+  const int exists = eosCheckTableExistence(EOS_M_DT, matid, eospacWarn);
+  if (exists > 0) {
+    eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_M_DT", "Eos_Material_Phases"},
+                eospacWarn);
+  } else {
+    phase_names = std::string("");
+    return false;
+  }
+  // indep vars
+  // Reuses the EOS log temp and log rho bounds
+  std::vector<EOS_REAL> rhos, Ts, phs;
+  makeInterpPoints(rhos, lRhoBounds);
+  makeInterpPoints(Ts, lTBounds);
+
+  EOS_REAL infoVals[1];
+  EOS_INTEGER infoItems[1] = {EOS_NUM_PHASES};
+  eosSafeTableInfo(&tableHandle[0], 1, infoItems, infoVals, eospacWarn);
+  const int nph = static_cast<int>(infoVals[0]);
+
+  EOS_CHAR infoString[EOS_META_DATA_STRLEN];
+  infoString[0] = '\0';
+  eosSafeTableMetaData(&tableHandle[1], EOS_Material_Phases, infoString, eospacWarn);
+  phase_names = std::string(infoString);
+
+  DataBox dMdt, dMdr;
+  Ms.resize(rhos.size(), Ts.size(), nph);
+  Ms.setRange(1, lTBounds.grid);
+  Ms.setRange(2, lRhoBounds.grid);
+
+  dMdr.copyMetadata(Ms);
+  dMdt.copyMetadata(Ms);
+  mask.copyMetadata(Ms);
+
+  // Interpolatable vars
+  EOS_INTEGER nXYPairs = rhos.size() * Ts.size();
+  std::vector<EOS_REAL> M_pack(nXYPairs * nph), DMDR_T(nXYPairs), DMDT_R(nXYPairs),
+      rho_flat(nXYPairs), T_flat(nXYPairs);
+
+  // prepare flat data structures
+  for (std::size_t j = 0, iflat = 0; j < rhos.size(); ++j) {
+    for (std::size_t i = 0; i < Ts.size(); ++i, iflat++) {
+      rho_flat[iflat] = densityToSesame(rhos[j]);
+      T_flat[iflat] = temperatureToSesame(Ts[i]);
+    }
+  }
+
+  const bool no_errors = eosSafeInterpolate(&tableHandle[0], nXYPairs, rho_flat.data(),
+                                            T_flat.data(), M_pack.data(), DMDR_T.data(),
+                                            DMDT_R.data(), "EOS_M_DT", eospacWarn);
+
+  for (size_t j = 0, iflat = 0; j < rhos.size(); j++) {
+    for (size_t i = 0; i < Ts.size(); i++, iflat++) {
+      for (size_t k = 0; k < nph; k++) {
+        Ms(j, i, k) = M_pack[k * nXYPairs + iflat];
+      }
+    }
+  }
+  eosSafeDestroy(NT, tableHandle, eospacWarn);
+  return true;
+}
+
 void makeInterpPoints(std::vector<EOS_REAL> &v, const Bounds &b) {
   v.resize(b.grid.nPoints());
   for (size_t i = 0; i < v.size(); i++) {

sesame2spiner/io_eospac.hpp

@@ -76,6 +76,9 @@ void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &si
 void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie,
                       const DataBox &sieColdCurve, DataBox &mask,
                       Verbosity eospacWarn = Verbosity::Quiet);
+bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds,
+                     Bounds &nphBounds, DataBox &mf, DataBox &mask,
+                     std::string &phase_names, Verbosity eospacWarn = Verbosity::Quiet);
 
 void makeInterpPoints(std::vector<EOS_REAL> &v, const Bounds &b);
 

singularity-eos/base/sp5/singularity_eos_sp5.hpp

@@ -26,6 +26,7 @@ namespace Depends {
 constexpr char logRhoLogSie[] = "dependsLogRhoLogSie";
 constexpr char logRhoLogT[] = "dependsLogRhoLogT";
 constexpr char coldCurve[] = "coldCurve";
+constexpr char massFrac[] = "massFrac";
 } // namespace Depends
 
 namespace SubTable {
@@ -68,6 +69,7 @@ constexpr char dEdRho[] = "dEdRho";
 constexpr char dEdT[] = "dEdT";
 constexpr char mask[] = "mask";
 constexpr char transitionMask[] = "transition mask";
+constexpr char massFrac[] = "mass fractions";
 } // namespace Fields
 
 } // namespace SP5