Piecewise grids for spiner EOS

CHANGELOG.md

@@ -3,8 +3,10 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Piecewise grids for Spiner EOS.
 
 ### Fixed (Repair bugs, etc)
+- [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Includes a fix for extrapolation of specific internal energy in SpinerEOS.
 - [[PR401]](https://github.com/lanl/singularity-eos/pull/401) Fix for internal energy scaling in PTE closure
 
 ### Changed (changing behavior/API/variables/...)

doc/sphinx/src/building.rst

@@ -643,11 +643,11 @@ packages. This can save lots of time! Note, however, that external
 packages are loosely constrained and may not be correctly configured for
 the requested package.
 
-*NB*: By default, Spack will try to download the package source from the
-repository associated with the package. This behavior can be overrided
-with Spack *mirrors* , but that is beyond the scope of this doc.
+.. note::
 
-.. code:: bash
+  By default, Spack will try to download the package source from the
+  repository associated with the package. This behavior can be overrided
+  with Spack *mirrors* , but that is beyond the scope of this doc.
 
 Now, we can use Spack similarly to ``module load``,
 

doc/sphinx/src/contributing.rst

@@ -256,19 +256,19 @@ instantiated when ``singularity-eos`` is compiled.  Therefore to
 exercise all code paths, it is best to create an ``EOS`` type
 instantiated as
 
-.. code-block:: c++
+.. code-block:: cpp
 
-    #include <singularity-eos/eos/eos.hpp>
-    using EOS = singularity::Variant<MyNewEOS>;``.
-    EOS my_eos = MyNewEOS(parameter1, parameter2, ...)
+  #include <singularity-eos/eos/eos.hpp>
+  using EOS = singularity::Variant<MyNewEOS, ...>;
+  EOS my_eos = MyNewEOS(parameter1, parameter2, ...)
 
 in order to properly test the functionality of a new EOS. Simply using the
 new class as the type such as
 
-.. code-block:: c++
+.. code-block:: cpp
 
-    #include <singularity-eos/eos/eos.hpp>
-    auto my_eos = my_new_eos(parameter1, parameter2, ...)
+  #include <singularity-eos/eos/eos.hpp>
+  auto my_eos = my_new_eos(parameter1, parameter2, ...)
 
 won't ensure that the new EOS is working correctly in singularity with the
 static polymorphism of the ``EOS`` type.

doc/sphinx/src/integration.rst

@@ -8,7 +8,7 @@ or setting ``CMAKE_PREFIX_PATH`` manually, a CMake project can integrate
 Singularity-EOS via ``find_package(singularity-eos)`` and using the provided
 targets to link to it.
 
-.. code:: cmake
+.. code::
 
    find_package(singularity-eos)
    ...
@@ -40,7 +40,7 @@ dependencies, but only the compiled Fortran bindings provided by
 Singularity-EOS. To avoid unnecessary dependency checks by CMake, a Fortran
 project would integrate Singularity-EOS as follows:
 
-.. code:: cmake
+.. code::
 
    find_package(singularity-eos COMPONENTS Library)
    ...

doc/sphinx/src/models.rst

@@ -1364,10 +1364,12 @@ the energy offset of the products EOS is given by
 Practically, this means :math:`e_0` should be positive for any energetic material.
 
 To provide the energy offset to the Davis Products EOS, `the energy shift
-modifier<modifiers shifted EOS>`_ should be used. Note that the convention there
+modifier`_ should be used. Note that the convention there
 is that the shift is positive, so :math:`-e_0` should be provided to the shift
 modifier.
 
+.. _the energy shift modifier: modifiers shifted EOS
+
 The constructor for the Davis Products EOS is
 
 .. code-block:: cpp
@@ -1471,8 +1473,7 @@ key-value pairs. For exampe the following input deck is for air:
 
   matid = 5030
   # These set the number of grid points per decade
-  # for each variable. The default is 50 points
-  # per decade.
+  # for each variable. 
   numrho/decade = 40
   numT/decade = 40
   numSie/decade = 40
@@ -1491,16 +1492,114 @@ key-value pairs. For exampe the following input deck is for air:
   shrinklTBounds = 0.15
   shrinkleBounds = 0.5
 
-The only required value in an input file is the matid, in this
-case 5030. All other values will be inferred from the original sesame
-database if possible and if no value in the input file is
-provided. Comments are prefixed with ``#``.
+Comments are prefixed with ``#``. `eospac`_ uses environment variables
+and files to locate files in the `sesame`_ database, and
+``sesame2spiner`` uses `eospac`_. So the location of the ``sesame``
+database need not be provided by the command line. For how to specify
+`sesame`_ file locations, see the `eospac`_ manual.
+
+Piecewise Spiner Grids
+````````````````````````
+
+``sesame2spiner`` also supports grids with different resolutions in
+different parts of the table. We call these **piecewise** grids. By
+default grids are now piecewise. Piecewise grids can be disabled with
+
+.. code-block::
+
+  # defaults are true
+  piecewiseRho = false
+  piecewiseT = false
+  piecewiseSie = false
+
+These options may be true or false. The default is true. When
+piecewise grids are active, the density-temperature (or
+density-energy) grid is built as a Cartesian product grid of grids of
+non-uniform resolutions. The density grid gets split into three
+pieces, a region ``[rhoMin, rhoFineMin]``, a region ``[rhoFineMin,
+rhoFineMax]``, and a region ``[rhoFineMin, rhoMax]``. The
+``numrho/decade`` parameter sets the number of points per decade in
+the central refined region. The regions at lower and higher density
+have ``rhoCoarseFactorLo`` and ``rhoCoarseFactorHi`` fewer points per
+decade respectively compared to the finer region.
+
+Typically the fine region should be roughly centered around the normal
+density for a material, which is usually a challenging region to
+capture. If you neglect to set ``rhoFineMin`` and ``rhoFineMax``,
+``sesame2spiner`` will set the central refined region to be a region
+of diameter ``rhoFineDiameterDecades`` (in log space) around the
+material's normal density.
+
+The temperature grid has two regions, a more finely spaced region at
+low temperatures and a less finely spaced region at high
+temperatures. The regions are spearated by a temperature
+``TSplitPoint``. The default is :math:`10^4` Kelvin. The energy grid
+follows the temperature grid, with the energy split point
+corresponding to the temperature split point. The coarser
+high-temperature temperature and energy grids are coarsened by a
+factor of ``TCoarseFactor`` and ``sieCoarseFactor`` respectively.
+
+A diagram of a density-temperature grid is shown below. The region
+with temperatures below ``TSplitPoint`` is refined in temperature. The
+region between ``rhoFineMin`` and ``rhoFineMax`` is refined in
+density.
+
+.. image:: phase_diagram.png
+  :width: 400
+  :alt: An example piecewise density-temperature grid.
+
+
+Thus the input block for piecewise grid might look like this:
+
+.. code-block::
+
+  # Below, all right-hand-sides are set to their default values.
+  piecewiseRho = true
+  piecewiseT = true
+  piecewiseSie = true
+
+  # the fine resolution for rho.
+  numrho/decade = 350
+  # width of the fine region for rho
+  rhoFineDiameterDecades = 1.5
+  # the lower density region is 3x less refined
+  rhoCoarseFactorLo = 3
+  # the higher density region is 5x less refined
+  rhoCoarseFactorHi = 5
+
+  # the fine resolution for T
+  numT/decade = 100
+  # the point demarking the coarse and fine regions in temperature
+  TSplitPoint = 1e4
+  # it's usually wise to to not let
+  # temperature get too small in log space if you do this
+  Tmin = 1
+  # The coarser region (above the split point) is 50 percent less refined
+  TCoarseFactor = 1.5
+
+  # energy has the split point sie(rhonormal, TSplitPoint)
+  # but we may still specify the resolution
+  numSie/decade = 100
+  sieCoarseFactor = 1.5
+
+.. note::
+
+  For all grid types, the only required value in an input file is the
+  matid. Table bounds and normal density will be inferred from the
+  sesame metadata if possible and if no value in the original input
+  file is provided. Table densities and positions and sizes of refined
+  regions are not inferred from the table, but are chosen with
+  the default values listed in the above code block.
+
+.. note::
 
-`eospac`_ uses environment variables and files to locate files in the
-`sesame`_ database, and ``sesame2spiner`` uses `eospac`_. So the
-location of the ``sesame`` database need not be provided by the
-command line. For how to specify `sesame`_ file locations, see the
-`eospac`_ manual.
+  Both the flat and hierarchical grids attempt to align their grids so
+  that there is a grid point in density and temperature exactly at
+  room temperature and normal density. This is because normal density
+  and room temperature is a particularly important point in phase
+  space, as it is the point in phase space a piece of material sitting
+  on your desk would be at. This is called an *anchor* point for the
+  mesh.
 
 SAP Polynomial EOS
 ``````````````````

doc/sphinx/src/phase_diagram.tex

@@ -0,0 +1,31 @@
+\documentclass{standalone}
+\usepackage{tikz}
+
+\begin{document}
+
+\begin{tikzpicture}
+  \draw[dashed] (4.5, 0) -- ++(0,8);
+  \draw (4.5, 0) node[below] {$\rho$ normal};
+  \draw[red, dashed] (2.5, 0) node[below] {rhoFineMin} -- ++(0,8);
+  \draw[red, dashed] (6.5, 0) node[below] {rhoFineMax} -- ++(0,8);
+  \draw[red, |-|] (2.5,8.5) -- (6.5,8.5);
+  \draw[red] (4.5, 8.5) node[above] {rhoFineDiameter};
+
+  \draw[blue, dashed] (0, 5) -- ++(8, 0) node[right] {TSplit};
+  \draw[dashed] (0, 2.5) -- ++(8, 0) node[right] {TAnchor $=298$ K};
+
+  % Draw axes
+  \draw[ultra thick,->] (0,0) -- (9,0) node[below] {\huge $\log_{10}\rho$};
+  \draw[ultra thick,->] (0,0) -- (0,9) node[left] {\huge $\log_{10}T$};
+
+  \draw (4.5, 2.5) node {\huge fine region};
+  \draw (4.5, 7) node {\huge coarse in T};
+
+  \draw (1.25, 1) node {\large coarse in $\rho$};
+  \draw (7.75, 1) node {\large coarse in $\rho$};
+
+  \draw (1.25, 7) node {\large coarse in all};
+  \draw (7.75, 7) node {\large coarse in all};
+\end{tikzpicture}
+
+\end{document}

doc/sphinx/src/sphinx-doc.rst

@@ -27,7 +27,7 @@ Building documentation locally
 
 While you can rely on the CI to build the documentation associated with your
 branch, you can also very easily build sphinx documentation locally through
-python. These instructions also _do not_ require admin access and are usable
+python. These instructions also **do not** require admin access and are usable
 with shared machines or python distributions.
 
 First, ensure that you are running a modern version of python (i.e. python 3 of

doc/sphinx/src/using-eos.rst

@@ -404,13 +404,17 @@ A modified equation of state can be built up iteratively. To check if
 the equation of state currently stored in the variant can modified,
 you may call
 
-.. cpp:function:: bool ModifiedInVariant<Mod>() const;
+.. code-block:: cpp
+
+  bool ModifiedInVariant<Mod>() const;
 
 where ``Mod`` is the type of the modifier you want to apply, for
 example ``ShiftedEOS``. If this function returns true, then you can
 apply a modifier with the function
 
-.. cpp:function:: Variant Modify<Mod>(Args &&..args) const;
+.. code-block:: cpp
+
+  Variant Modify<Mod>(Args &&..args) const;
 
 where again ``Mod`` is the modifier you wish to apply, and ``args``
 are the arguments to the constructor for that modifier, e.g., the

sesame2spiner/examples/aluminum.dat

@@ -0,0 +1,21 @@
+# ======================================================================
+#  Example input deck for sesame2spiner,
+#  a tool for converting eospac to spiner
+#  Author: Jonah Miller ([email protected])
+#  © 2024. 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.
+# ======================================================================
+
+matid=3720
+name=aluminum
+rhomin = 1e-5
+Tmin = 1

sesame2spiner/examples/helium.dat

@@ -0,0 +1,18 @@
+# ======================================================================
+#  Example input deck for sesame2spiner,
+#  a tool for converting eospac to spiner
+#  Author: Jonah Miller ([email protected])
+#  © 2024. 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.
+# ======================================================================
+matid = 5762
+name = helium

sesame2spiner/examples/steel.dat

@@ -17,10 +17,18 @@
 
 matid = 4272
 name = steel
-rhomin = 1e-2
+rhomin = 1e-5
 Tmin = 1
+
+#piecewiseRho = false
+#piecewiseT = false
+#piecewiseSie = false
+#numrho = 1334
+#numT = 738
+#numsie = 634
+
 # These shrink logarithm of bounds
 # by a fraction of the total interval <= 1
-shrinklRhoBounds = 0.15,
-shrinklTBounds = 0.15,
-shrinkleBounds = 0.5
+#shrinklRhoBounds = 0.15,
+#shrinklTBounds = 0.15,
+#shrinkleBounds = 0.5

sesame2spiner/examples/unit_tests/copper.dat

@@ -17,4 +17,6 @@
 
 matid=3337
 name=copper
-rhomin=0.
+
+piecewiseRho = false
+numrho/decade = 50