Merge pull request #459 from lanl/jmm/split-tables

Enable electron/ion subtables for SpinerEOS and EOSPAC backends

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR459]](https://github.com/lanl/singularity-eos/pull/459) Add electron and ion tables to EOSPAC and SpinerEOS backends
 - [[PR453]](https://github.com/lanl/singularity-eos/pull/453) A PT space PTE solver
 - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS
 

doc/sphinx/src/models.rst

@@ -218,9 +218,26 @@ density.
 ``singularity-eos`` assumes that, when 3T physics is active, electrons
 and ions are each described by a separate equation of state
 object. Several models are specifically designed to represent, e.g.,
-the electron equation of state or ion equation of state. The tabulated
-models may also support loading tables specifically for electron or
-ion equations of state.
+the electron equation of state or ion equation of state.
+
+The tabulated models may also support loading tables specifically for
+electron or ion equations of state. In these cases, an ``enum class``
+specifies which component of the material is being requst:
+
+.. code-block:: cpp
+
+   enum class TableSplit { Total, ElectronOnly, IonCold };
+
+where here ``ElectronOnly`` is the free electrons of the material,
+corresponding to Sesame 304 tables, ``IonCold`` is the ions plus cold
+curve (i.e., the ions if you don't know what a cold curve is),
+corresponding to the Sesame 303 tables, and ``Total`` is the sum of
+free electrons and ions, correspodning to the Sesame 301 tables.
+
+.. note::
+
+  The tabulated equations of state have an implicitly assumed
+  ionization fraction that depends on the electron temperature.
 
 Available EOS Information and Nomenclature
 ------------------------------------------
@@ -1665,15 +1682,19 @@ The constructor for ``SpinerEOSDependsRhoT`` is given by two overloads:
 .. code-block:: cpp
 
   SpinerEOSDependsRhoT(const std::string &filename, int matid,
+                       const singularity::TableSplit = singularity::TableSplit::Total,
                        bool reproduciblity_mode = false);
   SpinerEOSDependsRhoT(const std::string &filename, const std::string &materialName,
+                       const singularity::TableSplit = singularity::TableSplit::Total,
                        bool reproducibility_mode = false);
 
 where here ``filename`` is the input file, ``matid`` is the unique
 material ID in the database in the file, ``materialName`` is the name
 of the material in the file, and ``reproducability_mode`` is a boolean
 which slightly changes how initial guesses for root finds are
-computed. The constructor for ``SpinerEOSDependsRhoSie`` is identical.
+computed. The ``TableSplit`` option selects electron, ion, or total
+equation of state tables for partial ionization. The constructor for
+``SpinerEOSDependsRhoSie`` is identical.
 
 .. note::
 
@@ -1736,6 +1757,12 @@ and files to locate files in the `sesame`_ database, and
 database need not be provided by the command line. For how to specify
 `sesame`_ file locations, see the `eospac`_ manual.
 
+.. note::
+
+  To enable 3T subtables with ``SpinerEOS``, you must set
+  ``ionization=true`` in the input file for ``sesame2spiner`` for the
+  desired material.
+
 Piecewise Spiner Grids
 ````````````````````````
 
@@ -2134,27 +2161,40 @@ EOSPAC EOS
     Entropy is not yet available for this EOS
 
 This is a striaghtforward wrapper of the `EOSPAC`_ library for the
-`Sesame`_ database. The constructor for the ``EOSPAC`` model looks like
+`Sesame`_ database. The constructor for the ``EOSPAC`` model has several overloads
 
 .. code-block::
 
-  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)
+  EOSPAC(int matid, TableSplit split, 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);
+  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);
+  EOSPAC(int matid, const Options &opts);
 
 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, ``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. 
+``split`` is the ``TableSplit`` parameter for electron, ion, or total
+tables used in partial ionization. ``invert_at_setup`` specifies
+whether or not pre-compute tables of 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::
 

sesame2spiner/examples/unit_tests/steel.dat

@@ -2,7 +2,7 @@
 #  Example input deck for sesame2spiner,
 #  a tool for converting eospac to spiner
 #  Author: Jonah Miller ([email protected])
-#  © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+#  © 2021-2025. 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
@@ -17,3 +17,4 @@
 
 matid=4272
 name=steel
+ionization = true

sesame2spiner/generate_files.cpp

@@ -1,7 +1,7 @@
 //======================================================================
 // sesame2spiner tool for converting eospac to spiner
 // Author: Jonah Miller ([email protected])
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2025. 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
@@ -46,7 +46,8 @@ using namespace EospacWrapper;
 
 herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds,
                     const Bounds &lTBounds, const Bounds &leBounds,
-                    const std::string &name, Verbosity eospacWarn) {
+                    const std::string &name, const bool addSubtables,
+                    Verbosity eospacWarn) {
 
   const int matid = metadata.matid;
   std::string sMatid = std::to_string(matid);
@@ -94,49 +95,47 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh
   coldGroup =
       H5Gcreate(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
-  {
-    DataBox P, sie, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, mask;
-    eosDataOfRhoSie(matid, lRhoBounds, leBounds, P, T, bMod, dPdRho, dPdE, dTdRho, dTdE,
-                    dEdRho, mask, eospacWarn);
-    status += P.saveHDF(leGroup, SP5::Fields::P);
-    status += T.saveHDF(leGroup, SP5::Fields::T);
-    status += bMod.saveHDF(leGroup, SP5::Fields::bMod);
-    status += dPdRho.saveHDF(leGroup, SP5::Fields::dPdRho);
-    status += dPdE.saveHDF(leGroup, SP5::Fields::dPdE);
-    status += dTdRho.saveHDF(leGroup, SP5::Fields::dTdRho);
-    status += dTdE.saveHDF(leGroup, SP5::Fields::dTdE);
-    status += dEdRho.saveHDF(leGroup, SP5::Fields::dEdRho);
-    status += mask.saveHDF(leGroup, SP5::Fields::mask);
-  }
-
-  {
-    DataBox P, sie, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, dEdT, mask;
-    eosDataOfRhoT(matid, lRhoBounds, lTBounds, P, sie, bMod, dPdRho, dPdE, dTdRho, dTdE,
-                  dEdRho, dEdT, mask, eospacWarn);
-    status += P.saveHDF(lTGroup, SP5::Fields::P);
-    status += sie.saveHDF(lTGroup, SP5::Fields::sie);
-    status += bMod.saveHDF(lTGroup, SP5::Fields::bMod);
-    status += dPdRho.saveHDF(lTGroup, SP5::Fields::dPdRho);
-    status += dPdE.saveHDF(lTGroup, SP5::Fields::dPdE);
-    status += dTdRho.saveHDF(lTGroup, SP5::Fields::dTdRho);
-    status += dTdE.saveHDF(lTGroup, SP5::Fields::dTdE);
-    status += dEdRho.saveHDF(lTGroup, SP5::Fields::dEdRho);
-    status += dEdT.saveHDF(lTGroup, SP5::Fields::dEdT);
-    status += mask.saveHDF(lTGroup, SP5::Fields::mask);
-  }
+  status += saveTablesRhoSie(leGroup, matid, TableSplit::Total, lRhoBounds, leBounds,
+                             eospacWarn);
+  status +=
+      saveTablesRhoT(lTGroup, matid, TableSplit::Total, lRhoBounds, lTBounds, eospacWarn);
   {
     DataBox P, sie, dPdRho, dEdRho, bMod, mask, transitionMask;
     eosColdCurves(matid, lRhoBounds, P, sie, dPdRho, dEdRho, bMod, mask, eospacWarn);
-    eosColdCurveMask(matid, lRhoBounds, leBounds.grid.nPoints(), sie, transitionMask,
-                     eospacWarn);
+    // currently unused
+    // eosColdCurveMask(matid, lRhoBounds, leBounds.grid.nPoints(), sie, transitionMask,
+    //                  eospacWarn);
 
     status += P.saveHDF(coldGroup, SP5::Fields::P);
     status += sie.saveHDF(coldGroup, SP5::Fields::sie);
     status += bMod.saveHDF(coldGroup, SP5::Fields::bMod);
     status += dPdRho.saveHDF(coldGroup, SP5::Fields::dPdRho);
     status += dEdRho.saveHDF(coldGroup, SP5::Fields::dEdRho);
-    status += mask.saveHDF(coldGroup, SP5::Fields::mask);
-    status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask);
+    // currently unused
+    // status += mask.saveHDF(coldGroup, SP5::Fields::mask);
+    // status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask);
+  }
+
+  if (addSubtables) {
+    int i = 0;
+    std::vector<TableSplit> splits = {TableSplit::ElectronOnly, TableSplit::IonCold};
+    std::vector<std::string> grpnames = {SP5::SubTable::electronOnly,
+                                         SP5::SubTable::ionCold};
+    for (auto split : splits) {
+      std::string grpname = grpnames[i++];
+      {
+        hid_t grp =
+            H5Gcreate(leGroup, grpname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        status += saveTablesRhoSie(grp, matid, split, lRhoBounds, leBounds, eospacWarn);
+        status += H5Gclose(grp);
+      }
+      {
+        hid_t grp =
+            H5Gcreate(lTGroup, grpname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        status += saveTablesRhoT(grp, matid, split, lRhoBounds, lTBounds, eospacWarn);
+        status += H5Gclose(grp);
+      }
+    }
   }
 
   status += H5Gclose(leGroup);
@@ -219,8 +218,14 @@ herr_t saveAllMaterials(const std::string &savename,
                 << lRhoBounds << lTBounds << leBounds << std::endl;
     }
 
-    status +=
-        saveMaterial(file, metadata, lRhoBounds, lTBounds, leBounds, name, eospacWarn);
+    const bool add_subtables = params[i].Get("ionization", false);
+    if (eospacWarn == Verbosity::Debug) {
+      std::cout << "Adding subtables for partial ionization? " << add_subtables
+                << std::endl;
+    }
+
+    status += saveMaterial(file, metadata, lRhoBounds, lTBounds, leBounds, name,
+                           add_subtables, eospacWarn);
     if (status != H5_SUCCESS) {
       std::cerr << "WARNING: problem with HDf5" << std::endl;
     }
@@ -234,6 +239,45 @@ herr_t saveAllMaterials(const std::string &savename,
   return status;
 }
 
+herr_t saveTablesRhoSie(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds,
+                        const Bounds &leBounds, Verbosity eospacWarn) {
+  herr_t status = 0;
+  DataBox P, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, mask;
+  eosDataOfRhoSie(matid, split, lRhoBounds, leBounds, P, T, bMod, dPdRho, dPdE, dTdRho,
+                  dTdE, dEdRho, mask, eospacWarn);
+  status += P.saveHDF(loc, SP5::Fields::P);
+  status += T.saveHDF(loc, SP5::Fields::T);
+  status += bMod.saveHDF(loc, SP5::Fields::bMod);
+  status += dPdRho.saveHDF(loc, SP5::Fields::dPdRho);
+  status += dPdE.saveHDF(loc, SP5::Fields::dPdE);
+  status += dTdRho.saveHDF(loc, SP5::Fields::dTdRho);
+  status += dTdE.saveHDF(loc, SP5::Fields::dTdE);
+  status += dEdRho.saveHDF(loc, SP5::Fields::dEdRho);
+  // currently unused
+  // status += mask.saveHDF(loc, SP5::Fields::mask);
+  return status;
+}
+
+herr_t saveTablesRhoT(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds,
+                      const Bounds &lTBounds, Verbosity eospacWarn) {
+  herr_t status = 0;
+  DataBox P, sie, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, dEdT, mask;
+  eosDataOfRhoT(matid, split, lRhoBounds, lTBounds, P, sie, bMod, dPdRho, dPdE, dTdRho,
+                dTdE, dEdRho, dEdT, mask, eospacWarn);
+  status += P.saveHDF(loc, SP5::Fields::P);
+  status += sie.saveHDF(loc, SP5::Fields::sie);
+  status += bMod.saveHDF(loc, SP5::Fields::bMod);
+  status += dPdRho.saveHDF(loc, SP5::Fields::dPdRho);
+  status += dPdE.saveHDF(loc, SP5::Fields::dPdE);
+  status += dTdRho.saveHDF(loc, SP5::Fields::dTdRho);
+  status += dTdE.saveHDF(loc, SP5::Fields::dTdE);
+  status += dEdRho.saveHDF(loc, SP5::Fields::dEdRho);
+  status += dEdT.saveHDF(loc, SP5::Fields::dEdT);
+  // Currently unused
+  // status += mask.saveHDF(loc, SP5::Fields::mask);
+  return status;
+}
+
 void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params &params,
                   Bounds &lRhoBounds, Bounds &lTBounds, Bounds &leBounds) {
 

sesame2spiner/generate_files.hpp

@@ -1,7 +1,7 @@
 //======================================================================
 // sesame2spiner tool for converting eospac to spiner
 // Author: Jonah Miller ([email protected])
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2025. 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
@@ -42,12 +42,25 @@ constexpr Real T_SPLIT_POINT_DEFAULT = 1e4;
 
 herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds,
                     const Bounds &lTBounds, const Bounds &leBounds,
-                    const std::string &name, Verbosity eospacWarn = Verbosity::Quiet);
+                    const std::string &name, const bool addSubtables,
+                    Verbosity eospacWarn = Verbosity::Quiet);
+inline herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata,
+                           const Bounds &lRhoBounds, const Bounds &lTBounds,
+                           const Bounds &leBounds, const std::string &name,
+                           Verbosity eospacWarn = Verbosity::Quiet) {
+  return saveMaterial(loc, metadata, lRhoBounds, lTBounds, leBounds, name, false,
+                      eospacWarn);
+}
 
 herr_t saveAllMaterials(const std::string &savename,
                         const std::vector<std::string> &filenames, bool printMetadata,
                         Verbosity eospacWarn);
 
+herr_t saveTablesRhoSie(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds,
+                        const Bounds &leBounds, Verbosity eospacWarn = Verbosity::Quiet);
+herr_t saveTablesRhoT(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds,
+                      const Bounds &lTBounds, Verbosity eospacWarn = Verbosity::Quiet);
+
 void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params &params,
                   Bounds &lRhoBounds, Bounds &lTBounds, Bounds &leBounds);