diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b8436346d..17b168d37f 100644 --- a/CHANGELOG.md +++ b/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/...) diff --git a/doc/sphinx/src/building.rst b/doc/sphinx/src/building.rst index a25a60087a..a716f931ef 100644 --- a/doc/sphinx/src/building.rst +++ b/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``, diff --git a/doc/sphinx/src/contributing.rst b/doc/sphinx/src/contributing.rst index 1344cffe37..3dbd81cac8 100644 --- a/doc/sphinx/src/contributing.rst +++ b/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 - using EOS = singularity::Variant;``. - EOS my_eos = MyNewEOS(parameter1, parameter2, ...) + #include + using EOS = singularity::Variant; + 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 - auto my_eos = my_new_eos(parameter1, parameter2, ...) + #include + 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. diff --git a/doc/sphinx/src/integration.rst b/doc/sphinx/src/integration.rst index 495b85a07d..cd7c3a7b09 100644 --- a/doc/sphinx/src/integration.rst +++ b/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) ... diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 205aa421c6..e2db07f41b 100644 --- a/doc/sphinx/src/models.rst +++ b/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`_ 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 `````````````````` diff --git a/doc/sphinx/src/phase_diagram.png b/doc/sphinx/src/phase_diagram.png new file mode 100644 index 0000000000..469e91c5aa Binary files /dev/null and b/doc/sphinx/src/phase_diagram.png differ diff --git a/doc/sphinx/src/phase_diagram.tex b/doc/sphinx/src/phase_diagram.tex new file mode 100644 index 0000000000..b89408c184 --- /dev/null +++ b/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} \ No newline at end of file diff --git a/doc/sphinx/src/sphinx-doc.rst b/doc/sphinx/src/sphinx-doc.rst index d0c347139b..552e4c3c2a 100644 --- a/doc/sphinx/src/sphinx-doc.rst +++ b/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 diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 83bdde1d35..adabdbeb33 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/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() const; +.. code-block:: cpp + + bool ModifiedInVariant() 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(Args &&..args) const; +.. code-block:: cpp + + Variant Modify(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 diff --git a/sesame2spiner/examples/aluminum.dat b/sesame2spiner/examples/aluminum.dat new file mode 100644 index 0000000000..2d26089c5f --- /dev/null +++ b/sesame2spiner/examples/aluminum.dat @@ -0,0 +1,21 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 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 diff --git a/sesame2spiner/examples/helium.dat b/sesame2spiner/examples/helium.dat new file mode 100644 index 0000000000..2668c9eceb --- /dev/null +++ b/sesame2spiner/examples/helium.dat @@ -0,0 +1,18 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 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 diff --git a/sesame2spiner/examples/steel.dat b/sesame2spiner/examples/steel.dat index ea8d8f0c2b..77aa6301de 100644 --- a/sesame2spiner/examples/steel.dat +++ b/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 diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat index 16ba63d428..8f125df4ba 100644 --- a/sesame2spiner/examples/unit_tests/copper.dat +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -17,4 +17,6 @@ matid=3337 name=copper -rhomin=0. + +piecewiseRho = false +numrho/decade = 50 diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 7fea0e706b..7c608be956 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -252,9 +253,9 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params checkValInMatBounds(matid, "sieMin", sieMin, metadata.sieMin, metadata.sieMax); checkValInMatBounds(matid, "sieMax", sieMax, metadata.sieMin, metadata.sieMax); - Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.0); - Real shrinklTBounds = params.Get("shrinklTBounds", 0.0); - Real shrinkleBounds = params.Get("shrinkleBounds", 0.0); + Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.); + Real shrinklTBounds = params.Get("shrinklTBounds", 0.); + Real shrinkleBounds = params.Get("shrinkleBounds", 0.); shrinklRhoBounds = std::min(1., std::max(shrinklRhoBounds, 0.)); shrinklTBounds = std::min(1., std::max(shrinklTBounds, 0.)); @@ -273,21 +274,53 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params << "shrinkleBounds > 0 and siemin or siemax set" << std::endl; } - int ppdRho = params.Get("numrho/decade", PPD_DEFAULT); - int numRhoDefault = getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); + int ppdRho = params.Get("numrho/decade", PPD_DEFAULT_RHO); + int numRhoDefault = Bounds::getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); - int ppdT = params.Get("numT/decade", PPD_DEFAULT); - int numTDefault = getNumPointsFromPPD(TMin, TMax, ppdT); + int ppdT = params.Get("numT/decade", PPD_DEFAULT_T); + int numTDefault = Bounds::getNumPointsFromPPD(TMin, TMax, ppdT); - int ppdSie = params.Get("numSie/decade", PPD_DEFAULT); - int numSieDefault = getNumPointsFromPPD(sieMin, sieMax, ppdSie); + int ppdSie = params.Get("numSie/decade", PPD_DEFAULT_T); + int numSieDefault = Bounds::getNumPointsFromPPD(sieMin, sieMax, ppdSie); int numRho = params.Get("numrho", numRhoDefault); int numT = params.Get("numT", numTDefault); int numSie = params.Get("numsie", numSieDefault); - Real rhoAnchor = metadata.normalDensity; - Real TAnchor = 298.15; + constexpr Real TAnchor = 298.15; + Real rhoAnchor = params.Get("rho_fine_center", metadata.normalDensity); + if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { + std::cerr << "WARNING [" << matid << "] " + << "normal density ill defined. Setting it to a sensible default." + << std::endl; + rhoAnchor = 1; + } + + // Piecewise grids stuff + const bool piecewiseRho = params.Get("piecewiseRho", true); + const bool piecewiseT = params.Get("piecewiseT", true); + const bool piecewiseSie = params.Get("piecewiseSie", true); + + const Real ppd_factor_rho_lo = + params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); + const Real ppd_factor_rho_hi = + params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); + const Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real rho_fine_diameter = + params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); + const Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); + const Real rho_fine_center = rhoAnchor; + + // These override the rho center/diameter settings + Real rho_fine_min = params.Get("rhoFineMin", -1); + Real rho_fine_max = params.Get("rhoFineMax", -1); + if (rho_fine_min * rho_fine_max < 0) { + std::cerr << "WARNING [" << matid << "]: " + << "Either rhoFineMin or rhoFineMax is set while the other is still unset." + << " Both must be set to be sensible. Ignoring." << std::endl; + rho_fine_min = rho_fine_max = -1; + } // Forces density and temperature to be in a region where an offset // is not needed. This improves resolution at low densities and @@ -295,12 +328,57 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params // Extrapolation and other resolution tricks will be explored in the // future. - if (rhoMin < STRICTLY_POS_MIN) rhoMin = STRICTLY_POS_MIN; - if (TMin < STRICTLY_POS_MIN) TMin = STRICTLY_POS_MIN; + if (rhoMin < STRICTLY_POS_MIN_RHO) rhoMin = STRICTLY_POS_MIN_RHO; + if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; + + if (piecewiseRho) { + if (rho_fine_min > 0) { + lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, + rho_fine_min, rho_fine_max, ppdRho, ppd_factor_rho_lo, + ppd_factor_rho_hi, true, shrinklRhoBounds); + } else { + lRhoBounds = + Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, + ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, true, shrinklRhoBounds); + } + } else { + lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); + } + if (piecewiseT) { + lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, + ppd_factor_T, true, shrinklTBounds); + } else { + lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); + } + if (piecewiseSie) { + // compute temperature as a reasonable anchor point + constexpr int NT = 1; + constexpr EOS_INTEGER nXYPairs = 2; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER tableType[NT] = {EOS_Ut_DT}; + EOS_REAL rho[2], T[2], sie[2], dx[2], dy[2]; + { + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Ut_DT"}, Verbosity::Quiet); + EOS_INTEGER eospacEofRT = tableHandle[0]; + rho[0] = rho[1] = densityToSesame(rhoAnchor); + T[0] = temperatureToSesame(TAnchor); + T[1] = temperatureToSesame(TSplitPoint); + eosSafeInterpolate(&eospacEofRT, nXYPairs, rho, T, sie, dx, dy, "EofRT", + Verbosity::Quiet); + eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); + } + const Real sieAnchor = sie[0]; + const Real sieSplitPoint = sie[1]; + leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, + ppdSie, ppd_factor_sie, true, shrinkleBounds); + } else { + leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); + } - lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); - lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); - leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); + std::cout << "lRho bounds are\n" + << lRhoBounds << "lT bounds are\n" + << lTBounds << "lSie bounds are \n" + << leBounds << std::endl; return; } @@ -316,9 +394,3 @@ bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin } return true; } - -int getNumPointsFromPPD(Real min, Real max, int ppd) { - Bounds b(min, max, 3, true); - Real ndecades = b.grid.max() - b.grid.min(); - return static_cast(std::ceil(ppd * ndecades)); -} diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 90fb7d2276..914ec62ed8 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -30,8 +30,15 @@ using namespace EospacWrapper; -constexpr int PPD_DEFAULT = 50; -constexpr Real STRICTLY_POS_MIN = 1e-9; +constexpr int PPD_DEFAULT_RHO = 350; +constexpr int PPD_DEFAULT_T = 100; +constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; +constexpr Real STRICTLY_POS_MIN_T = 1e-2; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; +constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; +constexpr Real RHO_FINE_DIAMETER_DEFAULT = 1.5; +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, diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 86b780b415..dcbfce18eb 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -17,10 +17,6 @@ #ifndef _SESAME2SPINER_IO_EOSPAC_HPP_ #define _SESAME2SPINER_IO_EOSPAC_HPP_ -#include -#include -#include -#include #include #include // eospac API @@ -31,83 +27,16 @@ #include #include +#include #include -#include #include using EospacWrapper::Verbosity; -using DataBox = Spiner::DataBox; -using RegularGrid1D = Spiner::RegularGrid1D; - -// For logarithmic interpolation, quantities may be negative. -// If they are, use offset to ensure negative values make sense. -class Bounds { - public: - Bounds() {} - - Bounds(Real min, Real max, int N, Real offset) - : grid(RegularGrid1D(min, max, N)), offset(offset) {} - - Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, - Real anchor_point = std::numeric_limits::signaling_NaN()) - : offset(0) { - if (convertToLog) { - // Log scales can't handle negative numbers or exactly zero. To - // deal with that, we offset. - constexpr Real epsilon = std::numeric_limits::epsilon(); - const Real min_offset = 10 * std::abs(epsilon); - // 1.1 so that the y-intercept isn't identically zero - // when min < 0. - // min_offset to handle the case where min=0 - if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - - min += offset; - max += offset; - - min = singularity::FastMath::log10(std::abs(min)); - max = singularity::FastMath::log10(std::abs(max)); - Real delta = max - min; - min += 0.5 * shrinkRange * delta; - max -= 0.5 * shrinkRange * delta; - - if (!(std::isnan(anchor_point))) { - anchor_point += offset; - anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - } - } - - if (!(std::isnan(anchor_point))) { - if (min < anchor_point && anchor_point < max) { - Real dxguess = (max - min) / ((Real)N - 1); - int Nmax = static_cast((max - min) / dxguess); - int Nanchor = static_cast((anchor_point - min) / dxguess); - Real dx = (anchor_point - min) / static_cast(Nanchor + 1); - int Nmax_new = static_cast((max - min) / dx); - Nmax_new = std::max(Nmax, Nmax_new); - max = dx * (Nmax_new) + min; - N = Nmax_new; - } - } - - grid = RegularGrid1D(min, max, N); - } - - inline Real log2lin(Real xl) const { return singularity::FastMath::pow10(xl) - offset; } - inline Real i2lin(int i) const { return log2lin(grid.x(i)); } - - friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { - os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" - << " + " << b.offset << ", " - << "[N,dx] = [" << b.grid.nPoints() << ", " << b.grid.dx() << "]" - << "\n"; - return os; - } - - public: - RegularGrid1D grid; - Real offset; -}; +constexpr int NGRIDS = 3; +using Bounds = singularity::table_utils::Bounds; +using Grid_t = Spiner::PiecewiseGrid1D; +using DataBox = Spiner::DataBox; void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds, DataBox &P, DataBox &T, DataBox &bMods, DataBox &dPdRho, diff --git a/sesame2spiner/test.cpp b/sesame2spiner/test.cpp index d3f6f8f1ee..c7bd1bd766 100644 --- a/sesame2spiner/test.cpp +++ b/sesame2spiner/test.cpp @@ -39,8 +39,8 @@ #include "io_eospac.hpp" #include "parser.hpp" -using Spiner::DataBox; -using Spiner::RegularGrid1D; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; constexpr int matid = 5030; // dry air constexpr int matid2 = 4272; // stainless steel diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 399702f989..0096902e09 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -27,6 +27,7 @@ register_headers( base/fast-math/logs.hpp base/robust_utils.hpp base/root-finding-1d/root_finding.hpp + base/spiner_table_bounds.hpp base/variadic_utils.hpp base/math_utils.hpp base/constants.hpp diff --git a/singularity-eos/base/spiner_table_bounds.hpp b/singularity-eos/base/spiner_table_bounds.hpp new file mode 100644 index 0000000000..80d3dbd080 --- /dev/null +++ b/singularity-eos/base/spiner_table_bounds.hpp @@ -0,0 +1,286 @@ +//------------------------------------------------------------------------------ +// © 2021-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. +//------------------------------------------------------------------------------ + +#ifndef SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ +#define SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ +#ifdef SINGULARITY_USE_SPINER + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace singularity { +namespace table_utils { +// For logarithmic interpolation, quantities may be negative. +// If they are, use offset to ensure negative values make sense. +template +class Bounds { + public: + using RegularGrid1D = Spiner::RegularGrid1D; + using Grid_t = Spiner::PiecewiseGrid1D; + // for tag dispatch in constructors + using OneGrid = std::integral_constant; + using TwoGrids = std::integral_constant; + using ThreeGrids = std::integral_constant; + + Bounds() : offset(0), piecewise(false), linmin_(0), linmax_(0) {} + + Bounds(const Real min, const Real max, const int N, const Real offset) + : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), + offset(offset), piecewise(false), linmin_(min), linmax_(max) {} + Bounds(OneGrid, const Real min, const Real max, const int N, const Real offset) + : Bounds(min, max, N, offset) {} + + Bounds(Real min, Real max, int N, const bool convertToLog = false, + const Real shrinkRange = 0, + Real anchor_point = std::numeric_limits::signaling_NaN()) + : offset(0), piecewise(true), linmin_(min), linmax_(max) { + if (convertToLog) { + convertBoundsToLog_(min, max, shrinkRange); + if (!(std::isnan(anchor_point))) { + anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); + } + } + if (!(std::isnan(anchor_point))) { + adjustForAnchor_(min, max, N, anchor_point); + } + grid = Grid_t(std::vector{RegularGrid1D(min, max, N)}); + } + Bounds(OneGrid, Real min, Real max, int N, const bool convertToLog = false, + const Real shrinkRange = 0, + Real anchor_point = std::numeric_limits::signaling_NaN()) + : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} + + Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, + const Real ppd_fine, const Real ppd_factor, const bool convertToLog, + const Real shrinkRange = 0) + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; + + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + anchor_point = toLog_(anchor_point, offset); + splitPoint = toLog_(splitPoint, offset); + } + + checkInterval_(splitPoint, global_min, global_max, "Split point"); + checkInterval_(anchor_point, global_min, splitPoint, "Anchor point"); + + // add a point just to make sure we have enough points after adjusting for anchor + int N_fine = getNumPointsFromDensity(global_min, splitPoint, ppd_fine) + 1; + adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); + RegularGrid1D grid_lower(global_min, splitPoint, N_fine); + + const int N_upper = getNumPointsFromDensity(splitPoint, global_max, ppd_coarse); + RegularGrid1D grid_upper(splitPoint, global_max, N_upper); + + grid = Grid_t(std::vector{grid_lower, grid_upper}); + } + + Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real fine_min, + Real fine_max, const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + anchor_point = toLog_(anchor_point, offset); + fine_min = toLog_(fine_min, offset); + fine_max = toLog_(fine_max, offset); + } + checkInterval_(anchor_point, global_min, global_max, "Anchor point"); + + grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, + fine_min, fine_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); + } + + Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, + Real log_fine_diameter, const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + anchor_point = toLog_(anchor_point, offset); + } + + checkInterval_(anchor_point, global_min, global_max, "Anchor point"); + Real mid_min = anchor_point - 0.5 * log_fine_diameter; + Real mid_max = anchor_point + 0.5 * log_fine_diameter; + + grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, mid_min, + mid_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); + } + + inline Real log2lin(const Real xl) const { + // JMM: Need to guard this with the linear bounds passed in. The + // reason is that our fast math routines, while completely + // invertible at the level of machine epsilon, do introduce error + // at machine epsilon, which can bring us out of the interpolation + // range of eospac. + return std::min(linmax_, + std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); + } + inline Real i2lin(const int i) const { return log2lin(grid.x(i)); } + + friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { + os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" + << " + " << b.offset << "\n" + << "\tN = " << b.grid.nPoints() << "\n"; + for (int ig = 0; ig < b.grid.nGrids(); ++ig) { + os << "\t[ig,dx] = [" << ig << ", " << b.grid.dx(ig) << "]" + << "\n"; + } + return os; + } + + // This uses real logs + template + static int getNumPointsFromPPD(Real min, Real max, const T ppd) { + constexpr Real epsilon = std::numeric_limits::epsilon(); + const Real min_offset = 10 * std::abs(epsilon); + // 1.1 so that the y-intercept isn't identically zero + // when min < 0. + // min_offset to handle the case where min=0 + Real offset = 0; + if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; + min += offset; + max += offset; + + Real lmin = std::log10(min); + Real lmax = std::log10(max); + int N = getNumPointsFromDensity(lmin, lmax, ppd); + return N; + } + template + static int getNumPointsFromDensity(const Real min, const Real max, const T density) { + Real delta = max - min; + int N = std::max(2, static_cast(std::ceil(density * delta))); + return N; + } + + private: + Real toLog_(Real val, const Real offset) { + val += offset; + return singularity::FastMath::log10(std::abs(val)); + } + + void convertBoundsToLog_(Real &min, Real &max, const Real shrinkRange = 0) { + // Log scales can't handle negative numbers or exactly zero. To + // deal with that, we offset. + constexpr Real epsilon = std::numeric_limits::epsilon(); + const Real min_offset = 10 * std::abs(epsilon); + // 1.1 so that the y-intercept isn't identically zero + // when min < 0. + // min_offset to handle the case where min=0 + if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; + + min = toLog_(min, offset); + max = toLog_(max, offset); + + Real delta = max - min; + min += 0.5 * shrinkRange * delta; + max -= 0.5 * shrinkRange * delta; + } + + static void adjustForAnchor_(const Real min, Real &max, int &N, + const Real anchor_point) { + if (min < anchor_point && anchor_point < max) { + Real Nfrac = (anchor_point - min) / (max - min); + PORTABLE_REQUIRE((0 < Nfrac && Nfrac < 1), "anchor in bounds"); + + int Nanchor = static_cast(std::ceil(N * Nfrac)); + if ((Nanchor < 2) || (Nanchor >= N)) return; // not possible to shift this safely + + Real dx = (anchor_point - min) / static_cast(Nanchor - 1); + int Nmax_new = static_cast((max - min) / dx); + Real max_new = dx * (Nmax_new - 1) + min; + PORTABLE_REQUIRE(max_new <= max, "must not exceed table bounds"); + + N = Nmax_new; + max = max_new; + } else { + PORTABLE_ALWAYS_WARN("Anchor point out of bounds. Ignoring it."); + } + } + + static void checkInterval_(Real &p, const Real min, const Real max, + const std::string &name) { + if (p <= min) { + PORTABLE_ALWAYS_WARN(name + " less than minimum. Adjusting."); + Real eps = 0.1 * std::abs(min); + p = min + eps; + } + if (p >= max) { + PORTABLE_ALWAYS_WARN(name + " greater than maximum. Adjusting."); + Real eps = 0.1 * std::abs(max); + p = max - eps; + } + } + + static Grid_t gridFromIntervals_(ThreeGrids, Real global_min, Real global_max, + Real anchor_point, Real mid_min, Real mid_max, + const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi) { + const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; + const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; + + if (mid_min <= global_min) { + PORTABLE_ALWAYS_WARN( + "Table bounds refined minimum lower than global minimum. Adjusting."); + Real delta = std::abs(anchor_point - global_min); + mid_min = anchor_point - 0.9 * delta; + } + if (mid_max >= global_max) { + PORTABLE_ALWAYS_WARN( + "Table bounds refined maximum greater than global maximum. Adjusting."); + Real delta = std::abs(global_max - anchor_point); + mid_max = anchor_point + 0.9 * delta; + } + + // add a point just to make sure we have enough points after adjusting for anchor + int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; + adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); + RegularGrid1D grid_middle(mid_min, mid_max, N_fine); + + const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_lo); + RegularGrid1D grid_lower(global_min, mid_min, N_lower); + + const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_hi); + RegularGrid1D grid_upper(mid_max, global_max, N_upper); + + return Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); + } + + public: + Grid_t grid; + Real offset = 0; + bool piecewise = false; + + private: + Real linmin_, linmax_; +}; +} // namespace table_utils +} // namespace singularity + +#endif // SINGULARITY_USE_SPINER +#endif // SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 1649441faa..8b1d57343c 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -66,9 +66,11 @@ using namespace eos_base; we use log-linear extrapolation. */ class SpinerEOSDependsRhoT : public EosBase { - using DataBox = Spiner::DataBox; - public: + static constexpr int NGRIDS = 3; + using Grid_t = Spiner::PiecewiseGrid1D; + using DataBox = Spiner::DataBox; + // A weakly typed index map for lambdas struct Lambda { enum Index { lRho = 0, lT = 1 }; @@ -321,9 +323,9 @@ class SpinerEOSDependsRhoT : public EosBase { mitigated by Ye and (1-Ye) to control how important each term is. */ class SpinerEOSDependsRhoSie : public EosBase { - using DataBox = Spiner::DataBox; - public: + using Grid_t = SpinerEOSDependsRhoT::Grid_t; + using DataBox = SpinerEOSDependsRhoT::DataBox; struct SP5Tables { DataBox P, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho; }; @@ -528,8 +530,7 @@ class SpinerEOSDependsRhoSie : public EosBase { // replace lambdas with callable namespace callable_interp { - -using DataBox = Spiner::DataBox; +using DataBox = SpinerEOSDependsRhoT::DataBox; class l_interp { private: const DataBox &field; @@ -780,8 +781,7 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, // fill in minimum pressure as a function of temperature rho_at_pmin_.resize(numT_); - // rho_at_pmin_.setRange(0, P_.range(0)); - rho_at_pmin_.setRange(0, lTMin_, lTMax_, numT_); + rho_at_pmin_.setRange(0, P_.range(0)); for (int i = 0; i < numT_; i++) { Real pmin = std::numeric_limits::max(); int jmax = -1; diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index e3d42f6363..1bbf1dfc9a 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -220,6 +220,20 @@ class StellarCollapse : public EosBase { bool dependent_var_log); private: + class LogT { + public: + PORTABLE_INLINE_FUNCTION + LogT(const DataBox &field, const Real Ye, const Real lRho) + : field_(field), Ye_(Ye), lRho_(lRho) {} + PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { + return field_.interpToReal(Ye_, lT, lRho_); + } + + private: + const DataBox &field_; + const Real Ye_, lRho_; + }; + inline void LoadFromSP5File_(const std::string &filename); inline void LoadFromStellarCollapseFile_(const std::string &filename, bool filter_bmod); inline int readSCInt_(const hid_t &file_id, const std::string &name); @@ -379,25 +393,6 @@ class StellarCollapse : public EosBase { // Implementation details below // ====================================================================== -namespace callable_interp { - -class LogT { - public: - using DataBox = Spiner::DataBox; - PORTABLE_INLINE_FUNCTION - LogT(const DataBox &field, const Real Ye, const Real lRho) - : field_(field), Ye_(Ye), lRho_(lRho) {} - PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { - return field_.interpToReal(Ye_, lT, lRho_); - } - - private: - const DataBox &field_; - const Real Ye_, lRho_; -}; - -} // namespace callable_interp - // For some reason, the linker doesn't like this being a member field // of StellarCollapse. So we'll make it a global variable. constexpr char METADATA_NAME[] = "Metadata"; @@ -1164,7 +1159,7 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::lTFromlRhoSie_( } // Get log(sie) Real lE = e2le_(sie); - const callable_interp::LogT lEFunc(lE_, Ye, lRho); + const LogT lEFunc(lE_, Ye, lRho); status = regula_falsi(lEFunc, lE, lTGuess, lTMin_, lTMax_, ROOT_THRESH, ROOT_THRESH, lT, pcounts); if (status != RootFinding1D::Status::SUCCESS) { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4c654975ec..c0a127feb6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -34,6 +34,7 @@ add_executable( test_eos_vector.cpp test_math_utils.cpp test_variadic_utils.cpp + test_bounds.cpp ) add_executable( diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 9c19ddef33..5a4bfeb51c 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -47,10 +47,12 @@ #include -using namespace Spiner; using namespace singularity; using namespace EospacWrapper; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; + using duration = std::chrono::microseconds; constexpr char diffFileName[] = "diffs.sp5"; constexpr int NTIMES = 10; @@ -61,7 +63,7 @@ class Bounds { Bounds() {} Bounds(Real min, Real max, int N, Real offset) - : grid(RegularGrid1D(min, max, N)), offset(offset) {} + : grid(RegularGrid1D(min, max, N)), offset(offset) {} Bounds(Real min, Real max, int N) : offset(0) { constexpr Real epsilon = std::numeric_limits::epsilon(); @@ -71,7 +73,7 @@ class Bounds { max += offset; min = std::log10(std::abs(min)); max = std::log10(std::abs(max)); - grid = RegularGrid1D(min, max, N); + grid = RegularGrid1D(min, max, N); } PORTABLE_INLINE_FUNCTION Real log2lin(Real xl) const { return pow(10., xl) - offset; } @@ -86,7 +88,7 @@ class Bounds { } public: - RegularGrid1D grid; + RegularGrid1D grid; Real offset; }; @@ -146,59 +148,59 @@ int main(int argc, char *argv[]) { constexpr int NT = 3; EOS_INTEGER nXYPairs = (nFineRho) * (nFineT); - DataBox xVals(nFineRho, nFineT); - DataBox yVals(nFineRho, nFineT); - DataBox vars(nFineRho, nFineT); - DataBox dx(nFineRho, nFineT); - DataBox dy(nFineRho, nFineT); - DataBox pressEOSPAC(nFineRho, nFineT); - DataBox pressDiff_h(nFineRho, nFineT); - DataBox pressDiff_d(nFineRho, nFineT); + DataBox xVals(nFineRho, nFineT); + DataBox yVals(nFineRho, nFineT); + DataBox vars(nFineRho, nFineT); + DataBox dx(nFineRho, nFineT); + DataBox dy(nFineRho, nFineT); + DataBox pressEOSPAC(nFineRho, nFineT); + DataBox pressDiff_h(nFineRho, nFineT); + DataBox pressDiff_d(nFineRho, nFineT); #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; RView pressSpinerDView("pressSpinerDevice", nFineRho * nFineT); typename RView::HostMirror pressSpinerHView = Kokkos::create_mirror_view(pressSpinerDView); - DataBox pressSpiner_d(pressSpinerDView.data(), nFineRho, nFineT); - DataBox pressSpiner_hm(pressSpinerHView.data(), nFineRho, nFineT); + DataBox pressSpiner_d(pressSpinerDView.data(), nFineRho, nFineT); + DataBox pressSpiner_hm(pressSpinerHView.data(), nFineRho, nFineT); #else - DataBox pressSpiner_d(nFineRho, nFineT); - DataBox pressSpiner_hm = pressSpiner_d.slice(2, 0, nFineRho); + DataBox pressSpiner_d(nFineRho, nFineT); + DataBox pressSpiner_hm = pressSpiner_d.slice(2, 0, nFineRho); #endif - DataBox pressSpiner_h(nFineRho, nFineT); - DataBox tempSpiner_h(nFineRho, nFineT); - DataBox tempSpinerE_h(nFineRho, nFineT); - DataBox tempEOSPAC(nFineRho, nFineT); - DataBox tempDiff_h(nFineRho, nFineT); - DataBox tempDiff_d(nFineRho, nFineT); - DataBox tempDiffE_h(nFineRho, nFineT); - DataBox tempDiffE_d(nFineRho, nFineT); + DataBox pressSpiner_h(nFineRho, nFineT); + DataBox tempSpiner_h(nFineRho, nFineT); + DataBox tempSpinerE_h(nFineRho, nFineT); + DataBox tempEOSPAC(nFineRho, nFineT); + DataBox tempDiff_h(nFineRho, nFineT); + DataBox tempDiff_d(nFineRho, nFineT); + DataBox tempDiffE_h(nFineRho, nFineT); + DataBox tempDiffE_d(nFineRho, nFineT); #ifdef PORTABILITY_STRATEGY_KOKKOS RView tempSpinerDView("tempSpinerDevice", nFineRho * nFineT); RView tempSpinerDViewE("tempSpinerSieDevice", nFineRho * nFineT); RView::HostMirror tempSpinerHView = Kokkos::create_mirror_view(tempSpinerDView); RView::HostMirror tempSpinerHViewE = Kokkos::create_mirror_view(tempSpinerDViewE); - DataBox tempSpiner_d(tempSpinerDView.data(), nFineRho, nFineT); - DataBox tempSpiner_hm(tempSpinerHView.data(), nFineRho, nFineT); - DataBox tempSpinerE_d(tempSpinerDViewE.data(), nFineRho, nFineT); - DataBox tempSpinerE_hm(tempSpinerHViewE.data(), nFineRho, nFineT); + DataBox tempSpiner_d(tempSpinerDView.data(), nFineRho, nFineT); + DataBox tempSpiner_hm(tempSpinerHView.data(), nFineRho, nFineT); + DataBox tempSpinerE_d(tempSpinerDViewE.data(), nFineRho, nFineT); + DataBox tempSpinerE_hm(tempSpinerHViewE.data(), nFineRho, nFineT); RView rhos_v("rhos", nFineRho); RView Ts_v("Ts", nFineT); RView sies_v("sies", nFineT); - DataBox rhos(rhos_v.data(), nFineRho); - DataBox Ts(Ts_v.data(), nFineT); - DataBox sies(sies_v.data(), nFineT); + DataBox rhos(rhos_v.data(), nFineRho); + DataBox Ts(Ts_v.data(), nFineT); + DataBox sies(sies_v.data(), nFineT); #else - DataBox tempSpiner_d(nFineRho, nFineT); - DataBox tempSpiner_hm = tempSpiner_d.slice(2, 0, nFineRho); - DataBox tempSpinerE_d(nFineRho, nFineT); - DataBox tempSpinerE_hm = tempSpiner_h.slice(2, 0, nFineRho); - DataBox rhos(nFineRho); - DataBox Ts(nFineT); - DataBox sies(nFineT); + DataBox tempSpiner_d(nFineRho, nFineT); + DataBox tempSpiner_hm = tempSpiner_d.slice(2, 0, nFineRho); + DataBox tempSpinerE_d(nFineRho, nFineT); + DataBox tempSpinerE_hm = tempSpiner_h.slice(2, 0, nFineRho); + DataBox rhos(nFineRho); + DataBox Ts(nFineT); + DataBox sies(nFineT); #endif std::cout << "Comparing to EOSPAC for materials..." << std::endl; @@ -222,6 +224,7 @@ int main(int argc, char *argv[]) { SesameMetadata metadata; eosGetMetadata(matid, metadata, Verbosity::Debug); + std::cout << metadata << std::endl; hid_t idGroup = H5Gcreate(file, std::to_string(matid).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); @@ -239,15 +242,15 @@ int main(int argc, char *argv[]) { (Real *)PORTABLE_MALLOC(sizeof(Real) * nFineRho * nFineT * eos_host.nlambda()); Real *lambda_hp = (Real *)malloc(sizeof(Real) * nFineRho * nFineT * eos_host.nlambda()); - DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); - DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); + DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); + DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); - Real rhoMin = eosE_host.rhoMin(); - Real rhoMax = eosE_host.rhoMax(); - Real TMin = eosE_host.TMin(); - Real TMax = eosE_host.TMax(); - Real sieMin = eosE_host.sieMin(); - Real sieMax = eosE_host.sieMax(); + Real rhoMin = 1.1 * std::max(metadata.rhoMin, 1e-5); + Real rhoMax = 0.9 * metadata.rhoMax; + Real TMin = 1.1 * std::max(metadata.TMin, 1.0); + Real TMax = 0.9 * metadata.TMax; + Real sieMin = metadata.sieMin + 0.1 * std::abs(metadata.sieMin); + Real sieMax = 0.9 * metadata.sieMax; Bounds lRhoBounds(rhoMin, rhoMax, nFineRho); Bounds lTBounds(TMin, TMax, nFineT); @@ -352,10 +355,11 @@ int main(int argc, char *argv[]) { const Real p_true = pressureFromSesame(vars(j, i)); pressEOSPAC(j, i) = p_true; const Real diff = pressSpiner_h(j, i) - p_true; - pressDiff_h(j, i) = diff; - diffL2 += diff * diff; const Real mean_p = 0.5 * diff + p_true; - L2 += mean_p * mean_p; + pressDiff_h(j, i) = diff; + const Real reldiff = diff / (1e-10 + std::abs(mean_p)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -379,10 +383,11 @@ int main(int argc, char *argv[]) { // pressSpiner_d has been copied into pressSpiner_h const Real p_true = pressEOSPAC(j, i); const Real diff = pressSpiner_hm(j, i) - p_true; - pressDiff_d(j, i) = diff; const Real mean_p = 0.5 * diff + p_true; - diffL2 += diff * diff; - L2 += mean_p * mean_p; + pressDiff_d(j, i) = diff; + const Real reldiff = diff / (1e-10 + std::abs(mean_p)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -420,7 +425,7 @@ int main(int argc, char *argv[]) { Real rho = lRhoBounds.i2lin(j); for (int i = 0; i < nFineT; i++) { Real sie = leBounds.i2lin(i); - xVals(j, i) = rho; + xVals(j, i) = densityToSesame(rho); yVals(j, i) = sieToSesame(sie); } } @@ -461,10 +466,10 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTIMES; n++) { for (int j = 0; j < nFineRho; j++) { Real rho = lRhoBounds.i2lin(j); - DataBox lambda_hj = lambda_h.slice(j); + DataBox lambda_hj = lambda_h.slice(j); for (int i = 0; i < nFineT; i++) { Real sie = leBounds.i2lin(i); - DataBox lambda = lambda_hj.slice(i); + DataBox lambda = lambda_hj.slice(i); tempSpiner_h(j, i) = eos_host.TemperatureFromDensityInternalEnergy(rho, sie, lambda.data()); } @@ -484,7 +489,7 @@ int main(int argc, char *argv[]) { PORTABLE_LAMBDA(const int &j, const int &i) { Real rho = rhos(j); Real sie = sies(i); - DataBox lambda = lambda_d.slice(j).slice(i); + DataBox lambda = lambda_d.slice(j).slice(i); tempSpiner_d(j, i) = eos_device.TemperatureFromDensityInternalEnergy( rho, sie, lambda.data()); }); @@ -555,8 +560,9 @@ int main(int argc, char *argv[]) { diff = 0; } #endif - diffL2 += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -569,8 +575,9 @@ int main(int argc, char *argv[]) { const Real diff = tempSpinerE_h(j, i) - t_true; tempDiffE_h(j, i) = diff; const Real mean_t = 0.5 * diff + t_true; - diffL2E += diff * diff; - L2 += mean_t * mean_t; + const Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2E += reldiff * reldiff; + L2 += 1; } } diffL2E /= L2; @@ -609,8 +616,9 @@ int main(int argc, char *argv[]) { diff = 0; } #endif - diffL2_d += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2_d += reldiff * reldiff; + L2 += 1; } } diffL2_d /= L2; @@ -623,8 +631,9 @@ int main(int argc, char *argv[]) { const Real diff = tempSpinerE_hm(j, i) - t_true; tempDiffE_d(j, i) = diff; const Real mean_t = 0.5 * diff + t_true; - diffL2E_d += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2E_d += reldiff * reldiff; + L2 += 1; } } diffL2E_d /= L2; diff --git a/test/profile_eos.cpp b/test/profile_eos.cpp index dfdff7d6d1..96980e872c 100644 --- a/test/profile_eos.cpp +++ b/test/profile_eos.cpp @@ -40,8 +40,10 @@ using namespace singularity; using duration = std::chrono::microseconds; using dvec = std::vector; using ivec = std::vector; -using Spiner::RegularGrid1D; + using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; + #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; using RMirror = typename RView::HostMirror; diff --git a/test/profile_eospac.cpp b/test/profile_eospac.cpp index 71bcb66c9d..52daeb01e9 100644 --- a/test/profile_eospac.cpp +++ b/test/profile_eospac.cpp @@ -42,8 +42,9 @@ using namespace singularity; using duration = std::chrono::microseconds; using dvec = std::vector; using ivec = std::vector; -using Spiner::RegularGrid1D; + using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp new file mode 100644 index 0000000000..34d2f7da40 --- /dev/null +++ b/test/test_bounds.cpp @@ -0,0 +1,134 @@ +//------------------------------------------------------------------------------ +// © 2021-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. +//------------------------------------------------------------------------------ + +#include +#include + +#include +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +constexpr int NGRIDS = 3; +using Bounds = singularity::table_utils::Bounds; +using RegularGrid1D = Spiner::RegularGrid1D; +using Grid_t = Bounds::Grid_t; + +constexpr Real rho_normal = 3; // g/cm^3 +constexpr Real T_normal = 300; // Kelvin +constexpr Real rho_min = 1e-6; +constexpr Real rho_max = 1e3; +constexpr Real T_min = 1; +constexpr Real T_max = 1e9; +constexpr Real T_split = 1e4; +constexpr int N_per_decade_fine = 200; +constexpr Real N_factor = 5; + +constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e3; + +SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { + WHEN("We compute the number of points from points per decade") { + int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); + constexpr int NDECADES = 9; + THEN("We get the right number") { REQUIRE(np == NDECADES * N_per_decade_fine); } + } +} + +SCENARIO("Linear bounds in the bounds object", "[Bounds]") { + WHEN("We compute linear bounds ") { + constexpr Real min = -2; + constexpr Real max = 5; + constexpr int N = 21; + Bounds lin_bounds(min, max, N); + THEN("The min and max and point number correct") { + REQUIRE(lin_bounds.grid.min() == min); + REQUIRE(lin_bounds.grid.max() == max); + REQUIRE(lin_bounds.grid.nPoints() == N); + } + } +} + +SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { + WHEN("We compute logarithmic single-grid bounds") { + int np = Bounds::getNumPointsFromPPD(rho_min, rho_max, N_per_decade_fine); + Bounds lRhoBounds(rho_min, rho_max, np, true, 0.0, rho_normal); + THEN("The lower and upper bounds are right") { + REQUIRE(std::abs(lRhoBounds.grid.min() - singularity::FastMath::log10(rho_min)) <= + REAL_TOL); + REQUIRE(lRhoBounds.grid.max() <= + singularity::FastMath::log10(rho_max)); // shifted due to anchor + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(rho_normal); + int ianchor; + Spiner::weights_t w; + lRhoBounds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); + } + } + } +} + +SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { + WHEN("We compute a piecewise bounds object with three grids") { + Bounds bnds(Bounds::ThreeGrids(), rho_min, rho_max, rho_normal, 0.5, + N_per_decade_fine, N_factor, N_factor, true); + THEN("The bounds are right") { + Real lrmin = singularity::FastMath::log10(rho_min); + Real lrmax = singularity::FastMath::log10(rho_max); + REQUIRE(std::abs(bnds.grid.min() - lrmin) <= REAL_TOL); + REQUIRE(std::abs(bnds.grid.max() - lrmax) <= REAL_TOL); + REQUIRE(bnds.grid.nGrids() == 3); + AND_THEN( + "The total number of points is less than a uniform fine spacing would imply") { + REQUIRE(bnds.grid.nPoints() < N_per_decade_fine * (lrmax - lrmin)); + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(rho_normal); + int ianchor; + Spiner::weights_t w; + bnds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); + } + } + } + } + WHEN("We compute a piecewise bounds object with two grids") { + Bounds bnds(Bounds::TwoGrids(), T_min, T_max, T_normal, T_split, N_per_decade_fine, + N_factor, true); + THEN("The bounds are right") { + Real ltmin = singularity::FastMath::log10(T_min); + Real ltmax = singularity::FastMath::log10(T_max); + REQUIRE(std::abs(bnds.grid.min() - ltmin) <= REAL_TOL); + REQUIRE(std::abs(bnds.grid.max() - ltmax) <= REAL_TOL); + REQUIRE(bnds.grid.nGrids() == 2); + AND_THEN( + "The total number of points is less than a uniform fine spacing would imply") { + REQUIRE(bnds.grid.nPoints() < N_per_decade_fine * (ltmax - ltmin)); + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(T_normal); + int ianchor; + Spiner::weights_t w; + bnds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); + } + } + } + } +} diff --git a/utils/spiner b/utils/spiner index cc089ccb3a..0d0d7474f3 160000 --- a/utils/spiner +++ b/utils/spiner @@ -1 +1 @@ -Subproject commit cc089ccb3a1f57afdc21becaaab286217589cf7c +Subproject commit 0d0d7474f34cdc168c141185cc692b65ee8ad0c0