diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 0000000000..ceeed6452a --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,32 @@ +name: Generate JOSS paper + +on: + push: + branches: [ main ] + pull_request: + branches: [ main ] + release: + types: + - created + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: joss-paper/paper.md + - name: Upload + uses: actions/upload-artifact@v1 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: joss-paper/paper.pdf diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ebdf8f869..6f58bf01ef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,8 +3,12 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR#362]](https://github.com/lanl/singularity-eos/pull/362) Add lambda to thermalqs - [[PR#339]](https://github.com/lanl/singularity-eos/pull/339) Added COMPONENTS to singularity-eos CMake install, allowing to select a minimal subset needed e.g. for Fortran bindings only - [[PR#336]](https://github.com/lanl/singularity-eos/pull/336) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a pressure power law expansion in eta = 1-V/V0. PowerMG. +- [[PR334]](https://github.com/lanl/singularity-eos/pull/334) Include plugins infrastructure +- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. MGUsup. +- [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release - [[PR#357]](https://github.com/lanl/singularity-eos/pull/357) Added support for C++17 (e.g., needed when using newer Kokkos). ### Fixed (Repair bugs, etc) @@ -13,11 +17,8 @@ - [[PR341]](https://github.com/lanl/singularity-eos/pull/341) Short-circuit HDF5 machinery when cray-wrappers used in-tree - [[PR340]](https://github.com/lanl/singularity-eos/pull/335) Fix in-tree builds with plugin infrastructure - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - -### Added (new features/APIs/variables/...) -- [[PR334]](https://github.com/lanl/singularity-eos/pull/334) Include plugins infrastructure -- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. MGUsup. -- [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release +- [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver +- [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d2349059c..1a7ec465ba 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ # ------------------------------------------------------------------------------# -# © 2021-2023. Triad National Security, LLC. All rights reserved. This program +# © 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. @@ -340,10 +340,12 @@ if(SINGULARITY_USE_SPINER) # singularity_enable_hdf5(singularity-eos_Interface) endif() endif() +# Both the interface (headers) and the library (compiled) need to link to Kokkos +# see get_sg_eos.cpp if(SINGULARITY_USE_KOKKOS) - singularity_enable_kokkos(singularity-eos_Interface) + singularity_enable_kokkos(singularity-eos_Common) if(SINGULARITY_USE_KOKKOSKERNELS) - singularity_enable_kokkoskernels(singularity-eos_Interface) + singularity_enable_kokkoskernels(singularity-eos_Common) endif() endif() if(SINGULARITY_USE_EIGEN) @@ -408,7 +410,9 @@ if(SINGULARITY_BUILD_CLOSURE) add_library(singularity-eos_Library) set_target_properties(singularity-eos_Library PROPERTIES OUTPUT_NAME singularity-eos) add_library(singularity-eos::singularity-eos_Library ALIAS singularity-eos_Library) - target_link_libraries(singularity-eos_Library INTERFACE singularity-eos_Common) + # Public scope ensures explicit Kokkos dependency for library and anything + # that links to it + target_link_libraries(singularity-eos_Library PUBLIC singularity-eos_Common) target_link_libraries(singularity-eos INTERFACE singularity-eos_Library) target_sources(singularity-eos_Library PRIVATE ${eos_srcs}) message(VERBOSE "EOS Library Sources:\n\t${eos_srcs}") diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 38e24d7150..2e31c25f6e 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -500,11 +500,18 @@ currently: * ``thermalqs::temperature`` * ``thermalqs::specific_heat`` * ``thermalqs::bulk_modulus`` +* ``thermalqs::do_lambda`` * ``thermalqs::all_values`` however, most EOS models only specify that they prefer density and temperature or density and specific internal energy. +.. note:: + + The ``thermalqs::do_lambda`` flag is a bit special. It specifies that + eos-specific operations are to be performed on the additional + quantities passed in through the ``lambda`` variable. + .. _eos builder section: EOS Builder diff --git a/joss-paper/.gitignore b/joss-paper/.gitignore new file mode 100644 index 0000000000..865faf10cb --- /dev/null +++ b/joss-paper/.gitignore @@ -0,0 +1 @@ +auto diff --git a/joss-paper/paper.md b/joss-paper/paper.md new file mode 100644 index 0000000000..bb68faccfe --- /dev/null +++ b/joss-paper/paper.md @@ -0,0 +1,269 @@ +--- +title: 'Singularity-EOS: Performance Portable Equations of State and Mixed Cell Closures' +tags: + - C++ + - Fortran + - Python + - equations of state + - thermodynamics + - performance portability +authors: + - name: Jonah M. Miller + orcid: 0000-0001-6432-7860 + affiliation: "1, 2" + corresponding: true + - name: Daniel A. Holladay + - orcid: 0000-0002-0673-9741 + affiliation: "2, 3" + - name: Jeffrey H. Peterson + orcid: 0000-0001-9425-4674 + affiliation: 4 + - name: Christopher M. Mauney + affiliation: "2, 5" + - name: Richard Berger + orcid: 0000-0002-3044-8266 + affiliation: 3 + - name: Anna Pietarila Graham + affiliation: 5 + - name: Karen Tsai + affiliation: 3 + - name: Brandon Barker + orcid: 0000-0002-8825-0893 + affiliation: "1, 2, 6, 7" + - name: Alexander Holas + affiliation: "2, 3, 8" + - name: Ann E. Mattsson + affiliation: 9 + - name: Mariam Gogilashvili + affiliation: "1, 2, 7, 10" + - name: Joshua C. Dolence + affiliation: "1, 2" + - name: Chad D. Meyer + affiliation: 11 + - name: Sriram Swaminarayan + affiliation: 3 + - name: Christoph Junghans + orcid: 0000-0003-0925-1458 + affiliation: 3 +affiliations: + - name: CCS-2, Computational Physics and Methods, Los Alamos National Laboratory, USA + index: 1 + - name: Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos, NM + index: 2 + - name: CCS-7, Applied Computer Scienc, Los Alamos National Laboratory, USA + index: 3 + - name: XCP-2, Eulerian Codes, Los Alamos National Laboratory, USA + index: 4 + - name: HPC-ENV, HPC Environments, Los Alamos National Laboratory, USA + index: 5 + - name: Department of Physics and Astronomy, Michigan State University, USA + index: 6 + - name: Center for Nonlinear Studies, Los Alamos National Laboratory, USA + index: 7 + - name: Heidelberg Institute for Theoretical Studies, Germany + index: 8 + - name: XCP-5, Materials and Physical Data, Los Alamos National Laboratory, USA + index: 9 + - name: Department of Physics, Florida State University, USA + index: 10 + - name: XCP-4, Continuum Models and Numerical Methods, Los Alamos National Laboratory, USA +date: 13 April 2024 +bibliography: refs.bib + +--- + +# Summary + +We present Singularity-EOS, a new performance-portable library for +equations of state and related capabilities. Singularity-EOS provides +a large set of analytic equations of state, such as the Gruneisen +equation of state, and tabulated equation of state data under a +unified interface. It also provides support capabilities around these +equations of state, such as Python wrappers, solvers for finding +pressure-temperature equilibrium between multiple equations of state, +and a unique *modifier* framework, allowing the user to transform a +base equation of state, for example by shifting or scaling the +specific internal energy. All capabilities are performance portable, +meaning they compile and run on both CPU and GPU for a wide variety of +architectures. + +# Statement of need + +When expressed mathematically for continuous materials, the laws of +conservation of mass, energy, and momentum form the Navier-Stokes +equations of fluid dynamics. In the limit of zero molecular viscosity, +they become the Euler equations. These laws have been used to describe +phenomena as disparate as flow of air over an airplane wing, bacterial +motion in fluids, and the cataclysmic deaths of stars. However, the +fluid equations are not complete, and the system must be *closed* by a +description of the material at a sub-continuum (e.g., molecular or +atomic) scale. This closure is commonly called the *equation of state* +(EOS). + +Equations of state vary from the simple ideal gas law, to +sophisticated descriptions multi-phase descriptions of the lattice +structure of ice or wood, to models of quark-gluon plasma and nuclear +pasta at ultra high densities. A common form to write an equation of +state is as a pair of relations: + +$$P = P(\rho, T, \vec{\lambda}) \text{ and } \varepsilon = \varepsilon(\rho, T, \vec{\lambda}),$$ + +which relate the pressure $P$ and specific internal energy +$\varepsilon$ to density $\rho$, temperature $T$, and potentially some +unknown set of additional quantities $\vec{\lambda}$. However, other +representations are possible, and in common parlance an EOS is the +collection of knowledge needed to reconstruct some intrinsic +thermodynamic quantities from others. For example, the speed of sound +through a material or the specific heat capacity, which are +thermodynamic derivatives of the pressure and the specific internal +energy, are both determined by the EOS. + +In multi-material fluid dynamics simulations, one often will end up +with a so-called *mixed cell*, where two materials exist within +the same simulation zone. This can be an artifact of the numerical +representation; for example a steel bar and the surrounding air may +end up sharing a finite volume cell if the boundaries of the cell do +not align exactly with the surface of the steel bar. Or it may +represent physical reality; for example, air is a mixture of nitrogen +and oxygen gases, as well as water vapor. Regardless of the nature of +the mixed cell, one must somehow provide to the fluid code what the +material properties of the cell are as a whole. This is called a +*mixed cell closure.* One such closure is +*pressure-temperature equilibrium* (PTE), where all materials +in the cell are assumed to be at the same pressure and temperature. + +# State of the Field + +Typically fluid dynamics codes each develop an EOS package +individually to meet a given problem's needs. Databases of tabulated +equations of state, such as the Sesame [@sesame] and Stellar +Collapse [@stellarcollapseweb] databases often come with tabulated +data readers, for example, the EOSPAC library +[@PimentelDavidA2021EUMV] and Stellar Collapse library +[@stellarcollapsetables]. However, these libraries typically do +not include analytic equations of state or provide a unified API. They +also don't provide extra equation-of-state capabilities, such as +equilibrium solvers or production hardening. With a few exceptions, +these libraries are also typically not GPU-capable. + +We present Singularity-EOS, which aims to be a "one stop shop" for +EOS models for fluid and continuum dynamics codes. It provides a +unified interface for both analytic and tabulated equations of +state. It also provides useful surrounding capabilities, such as +Python wrappers, *modifiers*, which allow the user to transform a +given EOS, and solvers which can find the state in which multiple +EOS's are in PTE. To support usability, the library is extensively +documented and tested and supports builds through both ``cmake`` +and ``Spack`` [@spack]. + +Singularity-EOS leverages the ``Kokkos`` [@Kokkos] library for +performance portability, meaning the code can run on both CPUs and +GPUs, as well as other accelerators. This fills an important need, as +modern super computing capabilities increasingly rely on GPUs for +performance. Singularity-EOS is now used in the ongoing open-source +[Phoebus](https://github.com/lanl/phoebus) project which has a +separate code paper in-prep. + +# Design Principles and Feature Highlights + +Here we enumerate several design principles underlying +Singularity-EOS, and highlight a few feature of the library. + +## Flexibility in loop patterns + +Singularity-EOS provides both scalar and vector APIs, allowing the +user to make EOS calls on both single points in thermodynamic space, +and on collections of points. The vector calls may be more performant +(as they may vectorize), however care is made to ensure both APIs +operate at acceptable performance, to accommodate different code +structures downstream. + +## Flexibility in memory layout + +The vector calls in Singularity-EOS use an *accessor* API and (with a +few exceptions) accept any C++ object that has a ``operator[]`` +function defined. This allows users to lay out their memory as they +see fit and use Singularity-EOS even on strided or sparsely allocated +memory. + +## Expose APIs to aid performance + +Many equations of state are most naturally represented as functions of +density and temperature. However, fluid codes require pressure as a +function of density and internal energy. Extracting this often +requires computing a root find to invert the relation + +$$\varepsilon = \varepsilon(\rho, T).$$ + +In these cases, we expose an initial guess for temperature, which +helps the solution rapidly converge. Similarly, the performance of a +sequence of EOS calls may depend on the ordering of the calls. For +example, if both temperature and pressure are required from an +equation of state that requires inversion, requesting pressure first +will be less performant than requesting temperature first, as the +former requires two root finds, and the latter requires only one. To +enable this, we expose a function ``FillEos``, in which the user may +request multiple quantities at once, and the code uses ordering +knowledge to compute them as performantly as possible. + +## Performance-portable polymorphism + +Accelerators provide new challenges to standard object-oriented +programming. In particular, not all compiler stacks (such as Sycl +[@SYCL] or OpenMP Target Offload [@chandra2001parallel]) +support relocatable device code, which is required for standard C++ +polymorphism. Even in programming models, such as CUDA [@cuda], +which do support relocatable device code, polymorphism can be slower +than naively expected, and the user-level API can be cumbersome, +requiring operations such as ``placement new``. To sidestep these +issues, we use the C++ language feature ``std::variant`` to +implement a polymorphism mechanism that works on device. + +## Modifiers + +A given code may need to modify an EOS model to make it suitable for a +given application. For example, the zero-point of the energy may need +to be shifted, a porosity model may need to be added, or the unit +system may need to be changed. We implement this with a system of +*modifiers*, which can be applied on top of an EOS in a generic +way. Modifiers may also be chained. + +## Fast log-lookups + +To span the required orders of magnitude, tabulated equations of state +are often tabulated on log-spaced grids. Logarithms and exponentials +are, however, expensive operations and the performance of lookups can +suffer. We instead use the not-quite-transcendental lookups described +in [@NQT] to significantly enhance performance of log-like lookups. + +## Extensibility via modular parts and plugins + +Singularity-EOS is designed to be extensible. The +``std::variant``-based polymorphism, combined with modifiers, as +described above, already provides significant flexibility. However, +downstream codes may wish to add functionality to the library. This +may be implemented in several ways. **First**, as Singularity-EOS is +open source, contributions from downstream developers are +welcome. **Second**, a C++ code that depends on Singularity-EOS may +implement their own models and include them in a local variant +object. Singularity-EOS provides tooling to build variants up +iteratively. **Finally**, Singularity-EOS provides a flexible plugin +infrastructure that allows downstream users to add capability to the +core library locally by telling the build system to include a locally +downloaded plugin. This final capability allows downstream users to +share code with each other, even when committing that code to +Singularity-EOS proper is not possible due to, e.g., licensing issues. + +# Acknowledgements + +This work was supported through the Laboratory Directed Research and +Development program, the Center for Space and Earth Sciences, and the +center for Nonlinear Studies under project numbers 20240477CR-SES and +20220564ECR at Los Alamos National Laboratory (LANL). LANL is operated +by Triad National Security, LLC, for the National Nuclear Security +Administration of U.S. Department of Energy (Contract +No. 89233218CNA000001). This research used resources provided by the +Darwin testbed at LANL which is funded by the Computational Systems +and Software Environments subprogram of LANL's Advanced Simulation and +Computing program (NNSA/DOE). This work is approved for unlimited +release with report number LA-UR-24-23364. diff --git a/joss-paper/refs.bib b/joss-paper/refs.bib new file mode 100644 index 0000000000..3c4086b8a6 --- /dev/null +++ b/joss-paper/refs.bib @@ -0,0 +1,177 @@ +@ARTICLE{fornax, + author = {{Skinner}, M. Aaron and {Dolence}, Joshua C. and {Burrows}, Adam and {Radice}, David and {Vartanyan}, David}, + title = "{FORNAX: A Flexible Code for Multiphysics Astrophysical Simulations}", + journal = {\apjs}, + keywords = {methods: numerical, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - High Energy Astrophysical Phenomena, Astrophysics - Solar and Stellar Astrophysics, Physics - Computational Physics}, + year = 2019, + month = mar, + volume = {241}, + number = {1}, + eid = {7}, + pages = {7}, + doi = {10.3847/1538-4365/ab007f}, +archivePrefix = {arXiv}, + eprint = {1806.07390}, + primaryClass = {astro-ph.IM}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2019ApJS..241....7S}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{nubhlight, + author = {{Miller}, Jonah M. and {Ryan}, Ben. R. and {Dolence}, Joshua C.}, + title = "{{\ensuremath{\nu}} bhlight: Radiation GRMHD for Neutrino-driven Accretion Flows}", + journal = {\apjs}, + keywords = {accretion, accretion disks, black hole physics, magnetohydrodynamics: MHD, methods: numerical, neutrinos, radiative transfer, Astrophysics - Instrumentation and Methods for Astrophysics, General Relativity and Quantum Cosmology, Physics - Computational Physics}, + year = 2019, + month = apr, + volume = {241}, + number = {2}, + eid = {30}, + pages = {30}, + doi = {10.3847/1538-4365/ab09fc}, +archivePrefix = {arXiv}, + eprint = {1903.09273}, + primaryClass = {astro-ph.IM}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2019ApJS..241...30M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@techreport{sesame, +author = "Lyon, Stanford P. and Johnson, James D.", +title = "Sesame: The Los Alamos National Laboratory Equation of State Database", +institution = "Los Alamos National Laboratory", +number="LA-UR-92-3407", +year = "1992"} + +@article{stellarcollapsetables, + author={Evan O'Connor and Christian D Ott}, + title={A new open-source code for spherically symmetric stellar collapse to neutron stars and black holes}, + journal={Classical and Quantum Gravity}, + volume={27}, + number={11}, + pages={114103}, + url={http://stacks.iop.org/0264-9381/27/i=11/a=114103}, + year={2010}, +} + +@Misc{stellarcollapseweb, + author = {Evan O'Connor and Christian D Ott}, + title = {Stellar Collapse: Microphysics}, + year = {2010}, + url = {https://stellarcollapse.org/equationofstate}, + note = {Online access. \url{https://stellarcollapse.org/equationofstate}}, +} + +@book{PimentelDavidA2021EUMV, +year = {2021}, +title = {EOSPAC User’s Manual: V.6.5}, +language = {eng}, +address = {United States}, +author = {Pimentel, David A}, +keywords = {EOSPAC SESAME Equation of State Interpolation Thermodynamics ; MATHEMATICS AND COMPUTING}, +organization = {Los Alamos National Lab. (LANL), Los Alamos, NM (United States)}, +} + +@ARTICLE{SullivanWeak, + author = {{Sullivan}, Chris and {O'Connor}, Evan and {Zegers}, Remco G.~T. and {Grubb}, Thomas and {Austin}, Sam M.}, + title = "{The Sensitivity of Core-collapse Supernovae to Nuclear Electron Capture}", + journal = {\apj}, + keywords = {neutrinos, nuclear reactions, nucleosynthesis, abundances, equation of state, supernovae: general, Astrophysics - High Energy Astrophysical Phenomena, Nuclear Experiment, Nuclear Theory}, + year = 2016, + month = jan, + volume = {816}, + number = {1}, + eid = {44}, + pages = {44}, + doi = {10.3847/0004-637X/816/1/44}, +archivePrefix = {arXiv}, + eprint = {1508.07348}, + primaryClass = {astro-ph.HE}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2016ApJ...816...44S}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@book{press2007numerical, + title={Numerical Recipes with Source Code CD-ROM 3rd Edition: The Art of Scientific Computing}, + author={Press, W.H. and Teukolsky, S.A. and Vetterling, W.T. and Flannery, B.P.}, + isbn={9780521884075}, + lccn={2007062003}, + series={Numerical Recipes: The Art of Scientific Computing}, + url={https://books.google.com/books?id=DyykEZo4fwUC}, + year={2007}, + publisher={Cambridge University Press} +} + +@ARTICLE{Kokkos, + author={Trott, Christian R. and Lebrun-Grandié, Damien and Arndt, Daniel and Ciesko, Jan and Dang, Vinh and Ellingwood, Nathan and Gayatri, Rahulkumar and Harvey, Evan and Hollman, Daisy S. and Ibanez, Dan and Liber, Nevin and Madsen, Jonathan and Miles, Jeff and Poliakoff, David and Powell, Amy and Rajamanickam, Sivasankaran and Simberg, Mikael and Sunderland, Dan and Turcksin, Bruno and Wilke, Jeremiah}, + journal={IEEE Transactions on Parallel and Distributed Systems}, + title={Kokkos 3: Programming Model Extensions for the Exascale Era}, + year={2022}, + volume={33}, + number={4}, + pages={805-817}, + doi={10.1109/TPDS.2021.3097283}} + + @misc{cuda, + author={NVIDIA and Vingelmann, Péter and Fitzek, Frank H.P.}, + title={CUDA, release: 10.2.89}, + year={2020}, + url={https://developer.nvidia.com/cuda-toolkit}, +} + + @book{chandra2001parallel, + title={Parallel programming in OpenMP}, + author={Chandra, Rohit and Dagum, Leo and Kohr, David and Menon, Ramesh and Maydan, Dror and McDonald, Jeff}, + year={2001}, + publisher={Morgan kaufmann} +} + +@inproceedings{SYCL, +author = {Reyes, Ruyman and Brown, Gordon and Burns, Rod and Wong, Michael}, +title = {SYCL 2020: More than Meets the Eye}, +year = {2020}, +isbn = {9781450375313}, +publisher = {Association for Computing Machinery}, +address = {New York, NY, USA}, +url = {https://doi.org/10.1145/3388333.3388649}, +doi = {10.1145/3388333.3388649}, +booktitle = {Proceedings of the International Workshop on OpenCL}, +articleno = {4}, +numpages = {1}, +keywords = {HPC, SYCL, CUDA, OpenCL}, +location = {Munich, Germany}, +series = {IWOCL '20} +} + +@ARTICLE{NQT, + author = {{Miller}, Jonah M. and {Dolence}, Joshua C. and {Holladay}, Daniel}, + title = "{Not-Quite Transcendental Functions and their Applications}", + journal = {arXiv e-prints}, + keywords = {Physics - Computational Physics, Computer Science - Performance, Mathematics - Numerical Analysis}, + year = 2022, + month = jun, + eid = {arXiv:2206.08957}, + pages = {arXiv:2206.08957}, + doi = {10.48550/arXiv.2206.08957}, +archivePrefix = {arXiv}, + eprint = {2206.08957}, + primaryClass = {physics.comp-ph}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2022arXiv220608957M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@inproceedings{spack, +author = {Gamblin, Todd and LeGendre, Matthew and Collette, Michael R. and Lee, Gregory L. and Moody, Adam and de Supinski, Bronis R. and Futral, Scott}, +title = {The Spack Package Manager: Bringing Order to HPC Software Chaos}, +year = {2015}, +isbn = {9781450337236}, +publisher = {Association for Computing Machinery}, +address = {New York, NY, USA}, +url = {https://doi.org/10.1145/2807591.2807623}, +doi = {10.1145/2807591.2807623}, +booktitle = {Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis}, +articleno = {40}, +numpages = {12}, +location = {Austin, Texas}, +series = {SC '15} +} diff --git a/python/module.cpp b/python/module.cpp index 81db516a72..a6b3ad1c79 100644 --- a/python/module.cpp +++ b/python/module.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 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 @@ -141,6 +141,7 @@ PYBIND11_MODULE(singularity_eos, m) { thermalqs.attr("temperature") = pybind11::int_(thermalqs::temperature); thermalqs.attr("specific_heat") = pybind11::int_(thermalqs::specific_heat); thermalqs.attr("bulk_modulus") = pybind11::int_(thermalqs::bulk_modulus); + thermalqs.attr("do_lambda") = pybind11::int_(thermalqs::do_lambda); thermalqs.attr("all_values") = pybind11::int_(thermalqs::all_values); py::module eos_units = m.def_submodule("eos_units"); diff --git a/singularity-eos/base/constants.hpp b/singularity-eos/base/constants.hpp index bdbf21ca6d..303361534a 100644 --- a/singularity-eos/base/constants.hpp +++ b/singularity-eos/base/constants.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 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 @@ -27,7 +27,8 @@ constexpr unsigned long pressure = (1 << 2); constexpr unsigned long temperature = (1 << 3); constexpr unsigned long specific_heat = (1 << 4); constexpr unsigned long bulk_modulus = (1 << 5); -constexpr unsigned long all_values = (1 << 6) - 1; +constexpr unsigned long do_lambda = (1 << 6); +constexpr unsigned long all_values = (1 << 7) - 1; } // namespace thermalqs constexpr size_t MAX_NUM_LAMBDAS = 3; diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index bc5504d505..ad0018b131 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 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 @@ -16,7 +16,9 @@ #define _SINGULARITY_EOS_CLOSURE_MIXED_CELL_MODELS_ #include +#include #include +#include #include #include @@ -52,6 +54,9 @@ constexpr Real line_search_fac = 0.5; constexpr Real vfrac_safety_fac = 0.95; constexpr Real minimum_temperature = 1.e-9; constexpr Real maximum_temperature = 1.e9; +constexpr Real temperature_limit = 1.0e15; +constexpr Real default_tguess = 300.; +constexpr Real min_dtde = 1.0e-16; } // namespace mix_params namespace mix_impl { @@ -167,8 +172,9 @@ class PTESolverBase { Real vsum = 0; for (int m = 0; m < nmat; ++m) vsum += vfrac[m]; + PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive"); for (int m = 0; m < nmat; ++m) - vfrac[m] *= vfrac_total / vsum; + vfrac[m] *= robust::ratio(vfrac_total, vsum); } // Finalize restores the temperatures, energies, and pressures to unscaled values from // the internally scaled quantities used by the solvers @@ -212,6 +218,8 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { + PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction"); + PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density"); rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -223,13 +231,13 @@ class PTESolverBase { // set volume fractions for (int m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(T); - const Real vmax = std::min(0.9 * rhobar[m] / rho_min, 1.0); + const Real vmax = std::min(0.9 * robust::ratio(rhobar[m], rho_min), 1.0); vfrac[m] = (vfrac[m] > 0.0 ? std::min(vmax, vfrac[m]) : vmax); vsum += vfrac[m]; } // Normalize vfrac for (int m = 0; m < nmat; ++m) { - vfrac[m] *= vfrac_total / vsum; + vfrac[m] *= robust::ratio(vfrac_total, vsum); } } @@ -241,12 +249,14 @@ class PTESolverBase { if (Tnorm > 0.0) { Tguess = Tnorm; } else { - Tguess = 300.0; + Tguess = mix_params::default_tguess; for (int m = 0; m < nmat; ++m) Tguess = std::max(Tguess, temp[m]); } + PORTABLE_REQUIRE(Tguess > 0., "Non-positive temperature guess for PTE"); // check for sanity. basically checks that the input temperatures weren't garbage - assert(Tguess < 1.0e15); + PORTABLE_REQUIRE(Tguess < mix_params::temperature_limit, + "Very large input temperature or temperature guess"); // iteratively increase temperature guess until all rho's are above rho_at_pmin const Real Tfactor = 10.0; bool rho_fail; @@ -256,7 +266,7 @@ class PTESolverBase { rho_fail = false; for (int m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(Tguess); - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); if (rho[m] < rho_min) { rho_fail = true; Tguess *= Tfactor; @@ -299,7 +309,7 @@ class PTESolverBase { // intialize rhobar array and final density InitRhoBarandRho(); - uscale = rho_total * sie_total; + uscale = rho_total * abs(sie_total); // guess some non-zero temperature to start const Real Tguess = GetTguess(); @@ -311,14 +321,14 @@ class PTESolverBase { temp[m] = 1.0; sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tguess, Cache[m]); // note the scaling of pressure - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tguess, sie[m], Cache[m], - false) / - uscale; + press[m] = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tguess, sie[m], Cache[m], false), + uscale); } // note the scaling of the material internal energy densities for (int m = 0; m < nmat; ++m) - u[m] = sie[m] * rhobar[m] / uscale; + u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } PORTABLE_INLINE_FUNCTION @@ -339,13 +349,16 @@ class PTESolverBase { Real rhoBsum = 0.0; Real Asum = 0.0; for (int m = 0; m < nmat; m++) { - Asum += vfrac[m] * press[m] / temp[m]; - rhoBsum += rho[m] * vfrac[m] * sie[m] / temp[m]; + Asum += vfrac[m] * robust::ratio(press[m], temp[m]); + rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } + PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess"); Asum *= uscale / Tnorm; rhoBsum /= Tnorm; + PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density"); Tideal = uscale / rhoBsum / Tnorm; - Pideal = Tnorm * Tideal * Asum / vfrac_total / uscale; + PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum"); + Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } // Compute the ideal EOS PTE solution and replace the initial guess if it has a lower @@ -374,10 +387,10 @@ class PTESolverBase { res_norm_old += res[m] * res[m]; } // check if the volume fractions are reasonable - const Real alpha = Pideal / Tideal; + const Real alpha = robust::ratio(Pideal, Tideal); for (int m = 0; m < nmat; ++m) { - vfrac[m] *= press[m] / (temp[m] * alpha); - if (rhobar[m] / vfrac[m] < eos[m].RhoPmin(Tnorm * Tideal)) { + vfrac[m] *= robust::ratio(press[m], (temp[m] * alpha)); + if (robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { // abort because this is putting this material into a bad state for (int n = m; n >= 0; n--) vfrac[n] = vtemp[n]; @@ -386,13 +399,15 @@ class PTESolverBase { } // fill in the rest of the state for (int m = 0; m < nmat; ++m) { - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); + const Real sie_m = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * Tideal, Cache[m]); - u[m] = rhobar[m] * sie_m / uscale; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tideal, sie_m, - Cache[m], false) / - uscale; + u[m] = rhobar[m] * robust::ratio(sie_m, uscale); + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tideal, + sie_m, Cache[m], false), + uscale); } // fill in the residual solver->Residual(); @@ -415,7 +430,7 @@ class PTESolverBase { // did work, fill in temp and energy density for (int m = 0; m < nmat; ++m) { temp[m] = Tideal; - sie[m] = uscale * u[m] / rhobar[m]; + sie[m] = uscale * robust::ratio(u[m], rhobar[m]); } } } @@ -493,7 +508,7 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU( iter++; } - const Real alpha = (u_tot - ulo) / (uhi - ulo); + const Real alpha = robust::ratio((u_tot - ulo), (uhi - ulo)); return FastMath::pow2((1.0 - alpha) * lTlo + alpha * lThi); } @@ -611,27 +626,29 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert{}; Real e_pert = eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]); - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, - Cache[m], false) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; - dedv[m] = (rhobar[m] * e_pert / uscale - u[m]) / dv; + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, + e_pert, Cache[m], false), + uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); + dedv[m] = robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dv); ////////////////////////////// // perturb temperature ////////////////////////////// Real dT = Tequil * derivative_eps; e_pert = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * (Tequil + dT), Cache[m]); - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), - e_pert, Cache[m], false) / - uscale; - dpdT[m] = (p_pert - press[m]) / dT; - dedT_sum += (rhobar[m] * e_pert / uscale - u[m]) / dT; + p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], + Tnorm * (Tequil + dT), e_pert, + Cache[m], false), + uscale); + dpdT[m] = robust::ratio((p_pert - press[m]), dT); + dedT_sum += robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dT); } // Fill in the Jacobian @@ -657,7 +674,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -665,14 +682,14 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; m++) { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = -0.95 * Tequil / dx[nmat]; + scale = robust::ratio(-0.95 * Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -692,15 +709,16 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { Tequil = Ttemp + scale * dx[nmat]; for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( rho[m], Tnorm * Tequil, Cache[m]); - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); temp[m] = Tequil; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, sie[m], - Cache[m], false) / - uscale; + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, + sie[m], Cache[m], false), + uscale); } Residual(); return ResidualNorm(); @@ -784,15 +802,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // volume fractions have been potentially reset to ensure densitites are // larger than rho(Pmin(Tequil)); set the physical density to reflect // this change in volume fraction - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tequil, Cache[m]); uscale += sie[m] * rho[m]; // note the scaling of pressure press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]); } for (int m = 0; m < nmat; ++m) { - press[m] /= uscale; - u[m] = sie[m] * rhobar[m] / uscale; + press[m] = robust::ratio(press[m], uscale); + u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } Residual(); return ResidualNorm(); @@ -838,11 +856,11 @@ class PTESolverFixedT : public mix_impl::PTESolverBase ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); - Real p_pert = - eos[m].PressureFromDensityTemperature(rho_pert, Tequil, Cache[m]) / uscale; - dpdv[m] = (p_pert - press[m]) / dv; + Real p_pert = robust::ratio( + eos[m].PressureFromDensityTemperature(rho_pert, Tequil, Cache[m]), uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); } // Fill in the Jacobian @@ -864,15 +882,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m]); } } // control how big of a step toward rho = rho(Pmin) is allowed for (int m = 0; m < nmat; m++) { const Real rho_min = eos[m].RhoPmin(Tequil); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } } // Now apply the overall scaling @@ -891,12 +909,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase } for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature(rho[m], Tequil, Cache[m]); - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]) / uscale; + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]), uscale); } Residual(); return ResidualNorm(); @@ -990,8 +1009,9 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // note the scaling of the material internal energy densities for (int m = 0; m < nmat; ++m) { - u[m] = sie[m] * rhobar[m] / uscale; - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]) / uscale; + u[m] = robust::ratio(sie[m] * rhobar[m], uscale); + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]), uscale); } Residual(); // Set the current guess for the equilibrium temperature. Note that this is already @@ -1006,7 +1026,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Real vsum = 0.0; for (int m = 0; m < nmat; ++m) { vsum += vfrac[m]; - residual[m] = Pequil / uscale - press[m]; + residual[m] = robust::ratio(Pequil, uscale) - press[m]; } residual[nmat] = vfrac_total - vsum; } @@ -1036,23 +1056,25 @@ class PTESolverFixedP : public mix_impl::PTESolverBase ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert{}; Real e_pert{}; - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, - Cache[m], true) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, + e_pert, Cache[m], true), + uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); ////////////////////////////// // perturb temperature ////////////////////////////// Real dT = Tequil * derivative_eps; - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), - e_pert, Cache[m], true) / - uscale; - dpdT[m] = (p_pert - press[m]) / dT; + p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], + Tnorm * (Tequil + dT), e_pert, + Cache[m], true), + uscale); + dpdT[m] = robust::ratio((p_pert - press[m]), dT); } // Fill in the Jacobian @@ -1074,7 +1096,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -1082,14 +1104,14 @@ class PTESolverFixedP : public mix_impl::PTESolverBase for (int m = 0; m < nmat; m++) { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = -0.95 * Tequil / dx[nmat]; + scale = -0.95 * robust::ratio(Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -1109,13 +1131,14 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Tequil = Ttemp + scale * dx[nmat]; for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( rho[m], Tnorm * Tequil, Cache[m]); - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tnorm * Tequil, Cache[m]) / - uscale; - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tnorm * Tequil, Cache[m]), + uscale); + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); temp[m] = Tequil; } Residual(); @@ -1228,7 +1251,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { error_p += residual[m + 1] * residual[m + 1]; error_t += residual[m + nmat] * residual[m + nmat]; } - mean_t /= rho_total; + mean_t = robust::ratio(mean_t, rho_total); error_p = std::sqrt(error_p); error_t = std::sqrt(error_t); // Check for convergence @@ -1248,31 +1271,33 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert; - Real t_pert = - eos[m].TemperatureFromDensityInternalEnergy(rho_pert, sie[m], Cache[m]) / Tnorm; - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, sie[m], - Cache[m], false) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; - dtdv[m] = (t_pert - temp[m]) / dv; + Real t_pert = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho_pert, sie[m], Cache[m]), Tnorm); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, + sie[m], Cache[m], false), + uscale); + dpdv[m] = robust::ratio(p_pert - press[m], dv); + dtdv[m] = robust::ratio(t_pert - temp[m], dv); ////////////////////////////// // perturb energies ////////////////////////////// const Real de = std::abs(u[m]) * derivative_eps; - Real e_pert = (u[m] + de) / rhobar[m]; - - t_pert = - eos[m].TemperatureFromDensityInternalEnergy(rho[m], uscale * e_pert, Cache[m]) / - Tnorm; - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, - uscale * e_pert, Cache[m], false) / - uscale; - dpde[m] = (p_pert - press[m]) / de; - dtde[m] = (t_pert - temp[m]) / de; - if (std::abs(dtde[m]) < 1.e-16) { // must be on the cold curve + Real e_pert = robust::ratio(u[m] + de, rhobar[m]); + + t_pert = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho[m], uscale * e_pert, Cache[m]), + Tnorm); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, + uscale * e_pert, Cache[m], false), + uscale); + dpde[m] = robust::ratio(p_pert - press[m], de); + dtde[m] = robust::ratio(t_pert - temp[m], de); + if (std::abs(dtde[m]) < mix_params::min_dtde) { // must be on the cold curve dtde[m] = derivative_eps; } } @@ -1310,20 +1335,20 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed if (scale * dx[m] < -0.1 * vfrac[m]) { - scale = -0.1 * vfrac[m] / dx[m]; + scale = -0.1 * robust::ratio(vfrac[m], dx[m]); } // try to control steps toward T = 0 const Real dt = (dtdv[m] * dx[m] + dtde[m] * dx[m + nmat]); if (scale * dt < -0.1 * temp[m]) { - scale = -0.1 * temp[m] / dt; + scale = -0.1 * robust::ratio(temp[m], dt); } const Real tt = temp[m] + scale * dt; const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * temp[m]), eos[m].RhoPmin(Tnorm * tt)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // Now apply the overall scaling @@ -1344,14 +1369,15 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { } for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = utemp[m] + scale * dx[nmat + m]; - sie[m] = uscale * u[m] / rhobar[m]; - temp[m] = - eos[m].TemperatureFromDensityInternalEnergy(rho[m], sie[m], Cache[m]) / Tnorm; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], sie[m], - Cache[m], false) / - uscale; + sie[m] = uscale * robust::ratio(u[m], rhobar[m]); + temp[m] = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho[m], sie[m], Cache[m]), Tnorm); + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], + sie[m], Cache[m], false), + uscale); } Residual(); return ResidualNorm(); @@ -1401,7 +1427,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // backtrack Real err_mid = s.TestUpdate(0.5); if (err_mid < err && err_mid < err_old) { - scale = 0.75 + 0.5 * (err_mid - err) / (err - 2.0 * err_mid + err_old); + scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old); } else { scale = line_search_fac; } diff --git a/test/python_bindings.py b/test/python_bindings.py index 5ae87f006d..a0ff479ad7 100644 --- a/test/python_bindings.py +++ b/test/python_bindings.py @@ -61,7 +61,8 @@ def testConstants(self): thermalqs.pressure | thermalqs.temperature | thermalqs.specific_heat | - thermalqs.bulk_modulus) + thermalqs.bulk_modulus | + thermalqs.do_lambda) def testIdealGas(self): eos = singularity_eos.IdealGas(1,1)