diff --git a/CHANGELOG.md b/CHANGELOG.md index 607c5a92104..f10e8526e04 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,15 +3,18 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR658]](https://github.com/lanl/singularity-eos/pull/XXX) Add `MinimumInternalEnergy`/`MaximumInternalEnergy` to the EOS introspection API, so energy bounds are reachable through modifiers and the `singularity::EOS` variant ### Fixed (Repair bugs, etc) ### Changed (changing behavior/API/variables/...) +- [[PR658]](https://github.com/lanl/singularity-eos/pull/XXX) `ScaledEOS::CheckParams` now requires a strictly positive scale factor, where it previously accepted any nonzero value. ### Infrastructure (changes irrelevant to downstream codes) - [[PR653]](https://github.com/lanl/singularity-eos/pull/653) Move pybind11 to a submodule rather than fetching it via cmake fetchcontent ### Deprecated (soon to be removed behavior/API/variables/...) +- [[PR658]](https://github.com/lanl/singularity-eos/pull/XXX) The table bounds accessors `sieMin()`/`sieMax()` on `SpinerEOSDependsRhoSie` and `StellarCollapse` are superseded by `MinimumInternalEnergy()`/`MaximumInternalEnergy()`, which every EOS provides and which work through modifiers and the `singularity::EOS` variant. The old names still work but will be removed. ## Release 1.12.1 Date: 08/10/2026 diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index e270de3839d..194bff89966 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -2277,6 +2277,16 @@ functions ``rhoMin()``, ``rhoMax()``, ``TMin()``, ``TMax()``, ``YeMin()``, ``YeMax()``, ``sieMin()``, and ``sieMax()``, which all return a ``Real`` number. +.. note:: + + ``sieMin()`` and ``sieMax()`` are superseded by + ``MinimumInternalEnergy()`` and ``MaximumInternalEnergy()``, described + in the :ref:`EOS API section`. Prefer the latter: they are + part of the general EOS API, so unlike ``sieMin``/``sieMax`` they are + available on every model, are transformed correctly by modifiers, and + are reachable through the ``singularity::EOS`` variant. The old names + still work but will be removed in a future release. + .. warning:: As with the SpinerEOS models, the stellar collapse models use fast logs. You can switch the logs to true logs with the @@ -2484,3 +2494,5 @@ See :ref:`EOSPAC Vector Functions ` for more details. .. _EOSPAC: https://laws.lanl.gov/projects/data/eos/eospacReleases.php + +This file was made in part with generative AI. diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst index d3281be3ce7..f0055df6b5b 100644 --- a/doc/sphinx/src/modifiers.rst +++ b/doc/sphinx/src/modifiers.rst @@ -121,6 +121,11 @@ where the first two parameters are the Gruneisen parameter and specific heat required by the ideal gas constructor and the latter is the scale. +The scale factor must be strictly positive. A negative scale is not +physically meaningful: it would invert the sign of energy and entropy, +and would turn every minimum bound reported by the bounds introspection +API into a maximum. ``CheckParams`` enforces this. + The Relativistic EOS --------------------- @@ -189,6 +194,33 @@ internal energy, and temperature. On the other hand, specifies the unit system by specifying units for time, mass, length, and temperature. +All of the bounds introspection methods described in the :ref:`EOS API +section`, including ``MinimumInternalEnergy`` and +``MaximumInternalEnergy``, are converted to the new unit system. Because +they are part of the general EOS API, they are also reachable through +the ``singularity::EOS`` variant, so one no longer needs to call +``GetUnmodifiedObject`` to retrieve energy bounds, which would return +them in the unmodified (cgs) unit system: + +.. code-block:: cpp + + using namespace singularity; + EOS my_eos = UnitSystem(/* ... */); + // bounds in the new unit system + Real sie_min = my_eos.MinimumInternalEnergy(); + Real sie_max = my_eos.MaximumInternalEnergy(); + +See :ref:`How Modifiers Transform the Energy +Bounds` below. + +.. note:: + + ``UnitSystem`` does *not* forward the model-specific table bounds + accessors ``sieMin()``/``sieMax()``. Those are not part of the general + EOS API, so they cannot be converted for an arbitrary underlying + model. Use ``MinimumInternalEnergy``/``MaximumInternalEnergy`` + instead, which every model provides. + Z-Split EOS ------------- @@ -373,3 +405,29 @@ Modifiers can also be undone, extracting the underlying EOS. Continuing the exam auto unmodified = my_eos.GetUnmodifiedObject(); will extract the underlying ``IdealGas`` EOS model out from the scale and shift. + +.. _modifier-energy-bounds: + +How Modifiers Transform the Energy Bounds +------------------------------------------ + +Specific internal energy is transformed by several modifiers, so the +energy bounds reported by ``MinimumInternalEnergy`` and +``MaximumInternalEnergy`` (described in the :ref:`EOS API +section`) are transformed to match: + +* ``UnitSystem`` divides both bounds by the energy unit. +* ``ShiftedEOS`` adds the shift to both bounds. +* ``ScaledEOS`` multiplies both bounds by the scale factor. Since the + scale factor is required to be positive, this preserves their ordering. +* ``RelativisticEOS`` and ``BilinearRampEOS`` leave energy alone, so the + bounds pass through unchanged. +* ``FlooredEnergy`` clamps energy to the per-density cold curve, which + lies inside the underlying energy range, so the global bounds pass + through unchanged. +* ``ZSplit`` scales energy by a factor that depends on the ionization + state passed through the ``lambda``, but the bounds introspection API + takes no ``lambda``. The bounds it reports are therefore the + *un-split* bounds of the underlying EOS. + +This file was made in part with generative AI. diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 85e2213227f..99eb80e87af 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1504,14 +1504,35 @@ and the function .. cpp:function:: Real MaximumDensity() const; provide bounds for valid inputs into a table, which can be used by a -root finder to meaningful bound the root search. +root finder to meaningful bound the root search. The corresponding +functions for specific internal energy are + +.. cpp:function:: Real MinimumInternalEnergy() const; + +and + +.. cpp:function:: Real MaximumInternalEnergy() const; + +For tabulated models these report the extent of the tabulated energies, +which is the natural way to seed a root find or a bounds array in +energy. Note that for a table in density and temperature, energy is a +dependent variable, so these are the extrema of the tabulated energy +field over the whole grid rather than the endpoints of an axis. + +.. note:: + + ``MinimumInternalEnergy`` is a property of the whole table. It is + *not* the same quantity as ``MinInternalEnergyFromDensity``, which is + the cold curve: the minimum energy at a *given* density. .. warning:: For unbounded equations of state, ``MinimumDensity`` and - ``MinimumTemperature`` will return zero, while ``MaximumDensity`` - will return a very large finite number. Which number you get, - however, is not guaranteed. You may wish to apply more sensible + ``MinimumTemperature`` will return zero, while ``MaximumDensity`` and + ``MaximumInternalEnergy`` will return very large finite numbers, and + ``MinimumInternalEnergy`` a very large negative one. (Energies are + legitimately negative, so zero is not a safe floor.) Which number you + get, however, is not guaranteed. You may wish to apply more sensible bounds in your own code. Similarly, diff --git a/plan_histories/MR658-2026-09-21-modifier-energy-bounds.md b/plan_histories/MR658-2026-09-21-modifier-energy-bounds.md new file mode 100644 index 00000000000..711adfb8ce9 --- /dev/null +++ b/plan_histories/MR658-2026-09-21-modifier-energy-bounds.md @@ -0,0 +1,350 @@ +# Plan: Expose Internal-Energy Bounds Through Modifiers (starting with `UnitSystem`) + +**MR**: #XXX (rename this file once the MR number is assigned) +**Date**: 2026-09-19 +**Status**: ✅ Phase 1 implemented — ✅ Phase 2 implemented (2026-09-21, at reviewer request) + +## Motivation + +A host code uses + +```cpp +using baseEOS = singularity::SpinerEOSDependsRhoSie; +using EOS = singularity::UnitSystem; +``` + +and needs the table's energy bounds to seed a root-finder / bounds array. Today the only +way to get them is to strip the modifier: + +```cpp +auto unmodified = eosh.GetUnmodifiedObject(); +sie_bounds(1) = unmodified.sieMin(); +sie_bounds(nbounds+1) = unmodified.sieMax(); +``` + +which returns values in the *base* (cgs) unit system, not in the unit system the host +code is actually working in. Every use site therefore has to remember to multiply by the +energy unit by hand, which is exactly the kind of conversion `UnitSystem` exists to +absorb. + +`UnitSystem` already converts the other introspection quantities +(`eos_unitsystem.hpp:275-293`: `MinimumDensity`, `MinimumTemperature`, `MaximumDensity`, +`MinimumPressure`, `MaximumPressureAtTemperature`, `RhoPmin`). Energy bounds are simply +missing from that set — for `UnitSystem` and for the introspection API in general. + +## Current Situation + +### Where bounds introspection lives + +- `EosBase` supplies permissive defaults for every model: + `eos_base.hpp:535` (`MinimumDensity() -> 0`), `:537` (`MinimumTemperature() -> 0`), + `:546` (`MaximumDensity() -> 1e100`), `:551` (`MinimumPressure() -> 0`), + `:554` (`MaximumPressureAtTemperature() -> 1e100`), `:558` (`RhoPmin() -> 0`). + There is **no** energy-bound member in this set. +- Modifiers either forward verbatim via `SG_ADD_MODIFIER_INTROSPECTION_METHODS` + (`eos_base.hpp:133`; used by `shifted_eos.hpp:429`, `floored_energy.hpp:413`, + `relativistic_eos.hpp:210`, `zsplit_eos.hpp:283`) or hand-write transformed versions + (`scaled_eos.hpp:253-269`, `eos_unitsystem.hpp:275-293`, `ramps_eos.hpp:253-269`). +- `Variant` re-exposes each one through a `PortsOfCall::visit` (`eos_variant.hpp:475-503`). + +### Where energy bounds exist today + +- `SpinerEOSDependsRhoSie::sieMin()/sieMax()` — `eos_spiner_rho_sie.hpp:265-270` +- `StellarCollapse::sieMin()/sieMax()` — `eos_stellar_collapse.hpp:225-226` +- EOSPAC knows them (`eospac_wrapper.cpp:70-71` fills `metadata.sieMin/sieMax`) but + `EOSPAC` itself only stores `rho_min_`/`temp_min_` (`eos_eospac.hpp:1212-1213`). +- `SpinerEOSDependsRhoT` and `Helmholtz` expose density/temperature bounds only + (`eos_spiner_rho_temp.hpp:237-255`, `eos_helmholtz.hpp:269-273`). + +These are ad-hoc, model-specific accessors — not part of the `EosBase` contract — which +is why no modifier and no variant can see them. + +## Two Scopes + +### Phase 1 — `UnitSystem` pass-through (small, unblocks the host code) + +Add to `eos_unitsystem.hpp`, adjacent to the existing introspection block (~line 293): + +```cpp + PORTABLE_FORCEINLINE_FUNCTION Real sieMin() const { return inv_sie_unit_ * t_.sieMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real sieMax() const { return inv_sie_unit_ * t_.sieMax(); } +``` + +`inv_sie_unit_` is the correct factor: the file's convention is "multiply by a unit to +convert to cgs", which is why `MinInternalEnergyFromDensity` (`eos_unitsystem.hpp:132`) +returns `inv_sie_unit_ * S`. + +**Why this does not require touching any other EOS.** `UnitSystem` is a class +template, so member function *bodies* are only instantiated when called. Implicit +instantiation of `UnitSystem` — including as a variant alternative — instantiates +member *declarations* only. The code above therefore compiles for every `T` in the +variant, and only fails if someone actually calls `sieMin()` on a `UnitSystem` whose +`T` lacks it. That failure is a clear compile error at the call site, not a silent wrong +answer. + +Optional, same cost, for symmetry with the base tables: `rhoMin()`, `rhoMax()`, +`TMin()`, `TMax()` pass-throughs (scaled by `inv_rho_unit_` / `inv_temp_unit_`). These +partly duplicate `MinimumDensity`/`MaximumDensity`/`MinimumTemperature`, so include them +only if the host code wants the base-table naming to survive the modifier. + +**Explicit limitation of Phase 1**: this works on the concrete +`UnitSystem` type only. It does *not* work through +`singularity::EOS`: + +- `Variant::GetUnmodifiedObject()` returns a `Variant`, not a concrete EOS + (`eos_variant.hpp:797`), so `.sieMin()` on the result does not compile. +- Inside a `visit` / `EvaluateHost` / `EvaluateDevice` lambda the body must compile for + *every* alternative, so a bare `eos.sieMin()` fails on `IdealGas`. A host-side + `if constexpr` + detection trait works but pushes a sentinel-value convention into + downstream code — i.e. a worse version of Phase 2. + +#### Phase 1 Implementation Record (2026-09-19) + +Implemented as described above, energy bounds only; the optional +`rhoMin`/`rhoMax`/`TMin`/`TMax` pass-throughs were **not** added, since +`MinimumDensity`/`MaximumDensity`/`MinimumTemperature` already cover them through the +modifier. + +Files changed: + +- `singularity-eos/eos/modifiers/eos_unitsystem.hpp`: `sieMin()`/`sieMax()` after the + `RhoPmin` introspection block, with a comment recording the lazy-instantiation + requirement. +- `test/test_eos_modifiers.cpp`: `sieMin()/sieMax()` added to the `BoundedGas` helper; + new `THEN` block in the "Modifiers propagate introspection bounds correctly" scenario + asserting the `sie_unit` conversion and the round trip back to base units; comment on + the `UnitSystem` scenario noting that it guards usability for a `T` without + the accessors. +- `doc/sphinx/src/modifiers.rst`: documented both methods in the unit system section, + including the compile-time-error caveat and the replaced `GetUnmodifiedObject` idiom. +- `CHANGELOG.md`: entry under `## Current develop` → `### Added` (PR number is a + placeholder). + +Verification: + +- `test/test_eos_modifiers.cpp` built and run standalone against the system Catch2: + all tests pass, 26 assertions in 3 test cases. +- Separate syntax-only compile of `UnitSystem::sieMin()/sieMax()` + (the host-code pattern that motivated this) and of `UnitSystem` with only + `MinimumDensity()` called, confirming the lazy-instantiation property in both + directions. +- `clang-format` (v21.1.4) applied; no reflow of surrounding code. + +Not verified: the full CMake test suite. The pre-existing `build/` directory cannot +reconfigure because its spack-installed `ports-of-call`/`spiner` prefixes under `/tmp` +have been removed, so the standalone compiles above were used instead. + +### Phase 2 — first-class energy bounds in the introspection API (optional, larger) + +Only needed if the bounds must be reachable from a runtime-polymorphic +`singularity::EOS`. + +Proposed names: `MinimumInternalEnergy()` / `MaximumInternalEnergy()`. Deliberately not +`MinimumEnergy`, to avoid confusion with the existing density-dependent cold curve +`MinInternalEnergyFromDensity`, which is a different quantity (a curve, not a table +bound). Naming is a review decision — flag it early. + +1. **`EosBase` defaults** (`eos_base.hpp`, next to `:535-551`): + `MinimumInternalEnergy() -> -1e100`, `MaximumInternalEnergy() -> 1e100`. Note the + default minimum must be *negative*: energies are legitimately negative for cold curves + and shifted EOS, so `0` is not a safe floor. Follow the existing `MaximumDensity` + comment (`eos_base.hpp:539-545`) on big-finite-number vs. infinity. +2. **`SG_ADD_MODIFIER_INTROSPECTION_METHODS`** (`eos_base.hpp:133`): add verbatim + forwarding of both methods. This covers `RelativisticEOS` (energy passes through + untouched, `relativistic_eos.hpp:72-73`) and `RampEOS` (modifies pressure only) + correctly for free. +3. **Modifiers needing a real transform** — these must *not* use the plain macro: + - `UnitSystem`: `inv_sie_unit_ * t_.M*InternalEnergy()`. + - `ScaledEOS`: modified energy is `scale_ * base` (`scaled_eos.hpp:73,80`), so + `scale_ * t_.M*InternalEnergy()`. + - `ShiftedEOS`: modified energy is `base + shift_` (`shifted_eos.hpp:75,85`). It + currently uses the plain macro (`shifted_eos.hpp:429`), so it needs an explicit + override; forwarding unchanged would be wrong by exactly `shift_`. +4. **Modifiers with a documented caveat rather than a transform**: + - `ZSplit` scales energy by a lambda-dependent factor (`zsplit_eos.hpp:65,88`), and the + introspection API takes no lambda. Forward unchanged and document that the bound is + the un-split bound. + - `FlooredEnergy` clamps energy to the cold curve per-density; the global bound is + unchanged, so forwarding is correct. +5. **Concrete models**: override in `SpinerEOSDependsRhoSie` and `StellarCollapse` + (trivial — delegate to existing `sieMin()/sieMax()`). `EOSPAC` can store the values + already computed in `eospac_wrapper.cpp:70-71`. `SpinerEOSDependsRhoT` and `Helmholtz` + are open questions (see below) and can keep the base defaults initially. +6. **`Variant`**: add two `visit`-based accessors alongside `eos_variant.hpp:475-503`. +7. **Python bindings**: add to the variant bindings; `SpinerEOSDependsRhoSie` and + `StellarCollapse` already expose `sieMin`/`sieMax` (`python/module.cpp:143-144,157-158`). + +#### Phase 2 Implementation Record (2026-09-21) + +Carried out at a reviewer's request. Decisions taken on the open questions: + +- **Naming**: `MinimumInternalEnergy()` / `MaximumInternalEnergy()`, matching the + spelled-out style of the rest of the introspection API. +- **Concrete model scope**: all four models that can know their bounds — + `SpinerEOSDependsRhoSie`, `StellarCollapse`, `EOSPAC`, and `SpinerEOSDependsRhoT`. + `Helmholtz` keeps the permissive base defaults. +- **Negative `ScaledEOS` scale**: swap min/max when `scale_ < 0`, so the reported + bounds stay ordered. The pre-existing `MinimumDensity`/`MaximumDensity` issue was + deliberately *not* touched, to keep the diff scoped. + +Deviation from the plan: step 2 proposed adding the two methods to +`SG_ADD_MODIFIER_INTROSPECTION_METHODS`, but `ShiftedEOS` uses that macro and needs a +real transform, so the macro would collide with its override. Instead the energy bounds +live in a second macro, `SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS`, applied only by the +modifiers that genuinely forward verbatim. + +Files changed: + +- `singularity-eos/eos/eos_base.hpp`: `MinimumInternalEnergy()`/`MaximumInternalEnergy()` + defaults (`-1e100` / `1e100`); new `SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS` macro; fixed + `SG_ADD_MODIFIER_INTROSPECTION_METHODS`'s parameter name (`t` → `t_`, open question 5). +- Verbatim forwarding via the new macro: `floored_energy.hpp`, `relativistic_eos.hpp`, + `zsplit_eos.hpp` (each with a comment on why forwarding is right for that modifier). +- `ramps_eos.hpp`: hand-written forwarding, since it hand-writes its whole introspection + block rather than using the macro. +- `shifted_eos.hpp`: explicit override adding `shift_` to both bounds. +- `scaled_eos.hpp`: explicit override multiplying by `scale_`, with the negative-scale swap. +- `eos_unitsystem.hpp`: explicit override dividing by the energy unit; the Phase 1 + `sieMin`/`sieMax` pass-throughs are kept, with a comment steering callers to the new + methods. +- `eos_spiner_rho_sie.hpp`, `eos_stellar_collapse.hpp`: one-line delegation to the + existing `sieMin()`/`sieMax()`. +- `eos_eospac.hpp`: new `sie_min_`/`sie_max_` members, filled from the `SesameMetadata` + the constructor already fetches (answering open question 4: yes); ordering check added + to `CheckParams`. +- `eos_spiner_rho_temp.hpp`: new `sie_min_`/`sie_max_` members plus a `setEnergyBounds_()` + helper called from both construction paths (`loadDataboxes_` and the from-EOS + constructor). Answering open question 3: the bounds are the extrema of the tabulated + `sie_` field, unioned with `sieCold_`, *not* corner evaluations — corners would assume + monotonicity in both arguments, which the table does not guarantee. Cached at load time + because a `DataBox::min()` scan per call would be O(numRho*numT). `sieMin()`/`sieMax()` + accessors added here too, for parity with the rho-sie table. +- `eos_variant.hpp`: two `visit`-based accessors. +- `python/module.hpp`: bound on the generic `eos_class` template, so every bound type + including the variant gets them. +- `test/test_eos_modifiers.cpp`: `BoundedGas` gains the two methods; new `THEN` blocks for + the shifted+scaled composition, the negative-scale swap, relativistic pass-through, the + unit-system conversion, and the permissive analytic defaults; new scenario exercising the + motivating case — energy bounds read off a `singularity::Variant` holding a + `UnitSystem`. +- `doc/sphinx/src/using-eos.rst`: documented both methods next to the density bounds, with + a note distinguishing them from `MinInternalEnergyFromDensity` and an extended warning + covering the negative default minimum. +- `doc/sphinx/src/modifiers.rst`: new "How Modifiers Transform the Energy Bounds" section + enumerating each modifier's transform; the `UnitSystem` section now leads with the new + API. +- `CHANGELOG.md`: second entry under `## Current develop` → `### Added`. +- Copyright years bumped and generative-AI notices added where missing. + +Verification: `clang-format` (v21.1.4) applied to all changed C++ files. **The build and +test run were handed off to the user** — not verified by me. + +#### Phase 1 Reversal: the `UnitSystem` pass-throughs were removed (2026-09-21) + +A reviewer objected to the Phase 1 `sieMin()`/`sieMax()` pass-throughs in +`eos_unitsystem.hpp` — specifically to the lazy-instantiation construct, on the grounds +that the compile error a new developer hits is confusing even with the explanatory +comment. Their suggestion was `if constexpr` + a `static_assert`. + +Resolution: **delete the pass-throughs instead.** Phase 2 makes them redundant — +`UnitSystem::MinimumInternalEnergy()` returns the same values, applies the same +`inv_sie_unit_` conversion, is unconditionally well formed for every `T` because `EosBase` +supplies a default, and works through the variant. The lazy-instantiation trick only +existed because `sieMin` was not part of the contract, which is exactly what Phase 2 +fixed. This removes the construct rather than improving its diagnostic. Safe to do because +Phase 1 was never merged — `c4af3a91` lives only on `buechler/table_bounds`, and the only +in-tree caller of the pass-through was the Phase 1 test. + +Also removed: the `sieMin()`/`sieMax()` aliases Phase 2 had added to +`SpinerEOSDependsRhoT` "for parity". Adding new instances of a name being signposted as +superseded works against the deprecation, so that table exposes only +`MinimumInternalEnergy`/`MaximumInternalEnergy`. + +Deprecation of the *base-table* `sieMin`/`sieMax` (on `SpinerEOSDependsRhoSie` and +`StellarCollapse`) was scoped to **docs + CHANGELOG only** — a note in `models.rst` and an +entry under `### Deprecated`. No `[[deprecated]]` attribute. Two reasons: + +1. **CI.** `SINGULARITY_STRICT_WARNINGS=ON` sets `-Wall -Werror` (`CMakeLists.txt:668`) and + both `.github/workflows/warnings.yml` and `sanitizer.yml` enable it. + `-Wdeprecated-declarations` is in `-Wall`, so the attribute turns every surviving + in-tree call into a hard error: `python/module.cpp:143-144,157-158` (member pointers + warn at bind time), `test/profile_stellar_collapse.cpp:121,146`, and the Phase 2 + delegations in `eos_spiner_rho_sie.hpp`/`eos_stellar_collapse.hpp`. There is also no + C++ deprecation macro in the repo; the only precedent is the hand-rolled per-compiler + `DEPRECATED_MODULE` in `singularity_eos.f90` (PR644), and these are + `PORTABLE_FORCEINLINE_FUNCTION`, so an `ALLOW_DEPRECATED`-style escape hatch would be + needed for device compilers. +2. **Family consistency.** `sieMin`/`sieMax` are one pair among eight table accessors. + `rhoMin`/`rhoMax`/`TMin` already duplicate `MinimumDensity`/`MaximumDensity`/ + `MinimumTemperature` and have never been deprecated; `TMax`, `YeMin` and `YeMax` have + **no** contract equivalent at all — there is no `MaximumTemperature()` in the library. + `models.rst:2277` presents all eight as one set, so attributing only the energy pair + would look arbitrary. Deprecating the family properly requires first adding + `MaximumTemperature()` and deciding about `YeMin`/`YeMax` — a follow-on MR. + +Python note: the new methods were bound on the generic `eos_class` template, so Python +users already have the replacement. Giving them a real runtime `DeprecationWarning` would +mean replacing `def_property_readonly("sieMin", &T::sieMin)` with a lambda calling +`PyErr_WarnEx`; deferred with the rest of the attribute work. + +## Sign Caveat (applies to both phases, worth noting in review) + +`ScaledEOS::CheckParams` only requires `|scale| > 0` (`scaled_eos.hpp:59-65`), so a +negative scale is legal. Multiplying a min bound by a negative number turns it into a max +bound. `ScaledEOS::MinimumDensity` (`scaled_eos.hpp:253`) already has this latent issue, +so Phase 2 should either `min`/`max`-swap when `scale_ < 0` or explicitly document that +negative scales are unsupported for bounds introspection. `UnitSystem` is safe — its +`CheckParams` requires strictly positive units (`eos_unitsystem.hpp:100-105`). + +## Testing + +- `test/test_eos_modifiers.cpp` already has a `BoundedGas` helper with hand-written bounds + (lines 64-82) and unit-system bound assertions (lines 328-375). Extend `BoundedGas` with + energy bounds and add: + - `UnitSystem`: energy bounds divided by `sie_unit`. + - Phase 2 only: `ShiftedEOS` shifts by `shift_`, `ScaledEOS` scales by `scale_`, + unmodified analytic EOS returns the permissive defaults, and the same values come back + through `singularity::EOS`. +- Phase 1 needs at least one test that instantiates `UnitSystem` for a `T` *without* + `sieMin` (e.g. `IdealGas`) and exercises the rest of its API, to lock in the + lazy-instantiation guarantee the design relies on. +- Tabulated coverage (`test/test_eos_tabulated.cpp`) for the Spiner round trip: + `UnitSystem(...).sieMin() * sie_unit == base.sieMin()`. + +## Documentation + +- `doc/sphinx/src/using-eos.rst`, "Methods Used for Mixed Cell Closures" (~lines + 1490-1535): document the new methods next to `MinimumDensity`/`MaximumDensity`, and + extend the existing `warning` block about unbounded EOS to cover the energy defaults. +- `doc/sphinx/src/modifiers.rst`: note how each modifier transforms energy bounds. + +## Bookkeeping (repo checklist) + +- `CHANGELOG.md` under `## Current develop` → `### Added`, with the `[[PRxxx]]` link form. +- Update copyright year on each modified file; `eos_unitsystem.hpp` already carries the + generative-AI notice (line 15), other touched files need one added if AI-assisted. +- `make format` after configuring. +- Rename this file to match the MR number. + +## Open Questions + +1. Ship Phase 1 alone, or Phase 1 + Phase 2 together? Phase 1 is ~6 lines and unblocks the + host code immediately; Phase 2 is the API-consistent answer but touches ~10 files and + needs review on naming and defaults. +2. `MinimumInternalEnergy` vs. `MinimumEnergy` vs. `SieMin` for the Phase 2 names. +3. `SpinerEOSDependsRhoT` energy bounds: derive from the tabulated `sie` field's min/max, + or from `InternalEnergyFromDensityTemperature` at the table corners? The corner + evaluation is only valid if `sie` is monotone in both arguments over the table. +4. Should `EOSPAC` store the `sieMin`/`sieMax` it already computes at load time? +5. Aside: `SG_ADD_MODIFIER_INTROSPECTION_METHODS(t)` (`eos_base.hpp:133`) names its + parameter `t` but its body uses `t_`. Harmless today because every caller passes `t_`, + but worth fixing while in the neighborhood. + +## Risk Assessment + +- **Phase 1**: very low. Additive, header-only, no existing call site changes behavior. The + only failure mode is a compile error if `sieMin` is requested from a `UnitSystem` over a + base EOS that lacks it. +- **Phase 2**: low but broad. Additive to the public API, but it adds two methods every + model inherits, and the `ShiftedEOS` transform is a behavior *correction* relative to + naive forwarding — it must land with the tests that pin it down. diff --git a/python/module.hpp b/python/module.hpp index 06e93174ded..89d1a6a4ff4 100644 --- a/python/module.hpp +++ b/python/module.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// © 2021-2026. 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 @@ -520,6 +520,8 @@ py::class_ eos_class(py::module_ & m, std::string name) { .def("MinimumDensity", &T::MinimumDensity) .def("MinimumTemperature", &T::MinimumTemperature) + .def("MinimumInternalEnergy", &T::MinimumInternalEnergy) + .def("MaximumInternalEnergy", &T::MaximumInternalEnergy) .def_property_readonly_static("nlambda", [](py::object) { return T::nlambda(); }) .def_property_readonly_static("PreferredInput", [](py::object) { return T::PreferredInput(); }) .def("PrintParams", &T::PrintParams) diff --git a/singularity-eos/base/constants.hpp b/singularity-eos/base/constants.hpp index b7edc9b3d9a..fc8442f6c37 100644 --- a/singularity-eos/base/constants.hpp +++ b/singularity-eos/base/constants.hpp @@ -41,6 +41,12 @@ enum class TableStatus { OnTable = 0, OffBottom = 1, OffTop = 2 }; constexpr Real ROOM_TEMPERATURE = 293; // K constexpr Real ATMOSPHERIC_PRESSURE = 1e6; +// The "unbounded" value reported by the bounds introspection API when a +// model has no real bound in a given variable: a big finite number +// rather than an actual infinity. See the comment on +// EosBase::MaximumDensity for why. Negate it for lower bounds. +constexpr Real BIG_FINITE_BOUND = 1e100; + struct SharedMemSettings { SharedMemSettings() = default; SharedMemSettings(char *data_, bool is_domain_root_) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 60823b3f0e1..17c63e224c0 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -11,7 +11,6 @@ // prepare derivative works, distribute copies to the public, perform // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ - #ifndef _SINGULARITY_EOS_EOS_EOS_BASE_ #define _SINGULARITY_EOS_EOS_EOS_BASE_ @@ -131,7 +130,7 @@ char *StrCat(char *destination, const char *source) { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } // for bounds introspection -#define SG_ADD_MODIFIER_INTROSPECTION_METHODS(t) \ +#define SG_ADD_MODIFIER_INTROSPECTION_METHODS(t_) \ PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { \ return t_.MinimumDensity(); \ } \ @@ -151,6 +150,18 @@ char *StrCat(char *destination, const char *source) { PORTABLE_FORCEINLINE_FUNCTION \ Real RhoPmin(const Real temp) const { return t_.RhoPmin(temp); } +// The energy bounds are kept in their own macro because several +// modifiers transform energy and therefore must supply their own +// versions. Use this macro only for modifiers that leave the energy +// scale untouched. +#define SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS(t_) \ + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { \ + return t_.MinimumInternalEnergy(); \ + } \ + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { \ + return t_.MaximumInternalEnergy(); \ + } + // These macros are to reduce boilerplate in vector API. They declare // all the different "default" vector methods that we define in the // base class. For an example of what this might concretize to (albeit @@ -543,7 +554,19 @@ class EosBase { // the max. On the other hand, it's more fraught if someone tries to // put it into a formula without guarding against it. PORTABLE_FORCEINLINE_FUNCTION - Real MaximumDensity() const { return 1e100; } + Real MaximumDensity() const { return BIG_FINITE_BOUND; } + + // Report the range of specific internal energies an EOS supports. + // Tabulated models report the extent of their energy axis; analytic + // models are unbounded, so the defaults are very large finite numbers + // (see the MaximumDensity comment above on the tradeoffs there). + // JMM/MB: The default minimum must be negative, not zero. Energies + // are legitimately negative for cold curves and for shifted EOS, so + // zero is not a safe floor. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumInternalEnergy() const { return -BIG_FINITE_BOUND; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumInternalEnergy() const { return BIG_FINITE_BOUND; } // These are for the PT space PTE solver to bound the iterations in // a safe range. @@ -551,7 +574,9 @@ class EosBase { Real MinimumPressure() const { return 0; } // Gruneisen EOS's often have a maximum density, which implies a maximum pressure. PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { + return BIG_FINITE_BOUND; + } PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 6ada2f49055..1853e7d1574 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -238,6 +238,7 @@ class EOSPAC : public EosBase { // TODO(JMM): More validation checks? PORTABLE_ALWAYS_REQUIRE(rho_min_ >= 0, "Non-negative minimum density"); PORTABLE_ALWAYS_REQUIRE(temp_min_ >= 0, "Non-negative minimum temperature"); + PORTABLE_ALWAYS_REQUIRE(sie_max_ > sie_min_, "Energy bounds must be ordered"); AZbar_.CheckParams(); } inline EOSPAC GetOnDevice() { return *this; } @@ -1212,6 +1213,8 @@ class EOSPAC : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return press_min_; } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { return sie_min_; } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { return sie_max_; } private: static constexpr const unsigned long _preferred_input = @@ -1240,6 +1243,10 @@ class EOSPAC : public EosBase { Real rho_min_ = 0; Real temp_min_ = 0; Real press_min_ = 0; + // Defaults match the permissive EosBase bounds, in case the table + // metadata is unavailable. + Real sie_min_ = -BIG_FINITE_BOUND; + Real sie_max_ = BIG_FINITE_BOUND; // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a // problem? Could it ever realistically overflow? EOS_INTEGER shared_size_, packed_size_; @@ -1356,6 +1363,8 @@ inline EOSPAC::EOSPAC(const int matid, TableSplit split, bool invert_at_setup, rho_min_ = m.rhoMin; temp_min_ = m.TMin; press_min_ = m.PMin; + sie_min_ = m.sieMin; + sie_max_ = m.sieMax; // use std::max to hydrogen, in case of bad table AZbar_.Abar = std::max(1.0, m.meanAtomicMass); diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index ba6dc20d36e..0d653be7422 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -139,9 +139,11 @@ class MGUsup : public EosBase { // Hugoniot pressure ill defined at reference density. On one side, // negative. On the other positive. PORTABLE_FORCEINLINE_FUNCTION - Real MinimumPressure() const { return -1e100; } + Real MinimumPressure() const { return -BIG_FINITE_BOUND; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { + return BIG_FINITE_BOUND; + } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index e6a9a69905f..822f35bab7a 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -166,9 +166,11 @@ class PowerMG : public EosBase { } // Essentially unbounded... I think. PORTABLE_FORCEINLINE_FUNCTION - Real MinimumPressure() const { return -1e100; } + Real MinimumPressure() const { return -BIG_FINITE_BOUND; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { + return BIG_FINITE_BOUND; + } inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index e36cbe8b837..72b7aeb6785 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -157,9 +157,11 @@ class SAP_Polynomial : public EosBase { // Essentially unbounded... I think. PORTABLE_FORCEINLINE_FUNCTION - Real MinimumPressure() const { return -1e100; } + Real MinimumPressure() const { return -BIG_FINITE_BOUND; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { + return BIG_FINITE_BOUND; + } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp index e83c05a3009..c919c1299c6 100644 --- a/singularity-eos/eos/eos_spiner_rho_sie.hpp +++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp @@ -11,7 +11,6 @@ // prepare derivative works, distribute copies to the public, perform // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ - #ifndef _SINGULARITY_EOS_EOS_EOS_SPINER_RHO_SIE_HPP_ #define _SINGULARITY_EOS_EOS_EOS_SPINER_RHO_SIE_HPP_ @@ -275,6 +274,8 @@ class SpinerEOSDependsRhoSieTransformable PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { return sieMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { return sieMax(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return PMin_; } PORTABLE_INLINE_FUNCTION diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp index 2d6e9a4f96c..b3906ce44ab 100644 --- a/singularity-eos/eos/eos_spiner_rho_temp.hpp +++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp @@ -12,6 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// This file was created in part with generative AI + #ifndef _SINGULARITY_EOS_EOS_EOS_SPINER_RHO_TEMP_HPP_ #define _SINGULARITY_EOS_EOS_EOS_SPINER_RHO_TEMP_HPP_ @@ -253,6 +255,15 @@ class SpinerEOSDependsRhoT : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + // The extent of the tabulated specific internal energy. Unlike + // density and temperature, energy is not an independent variable + // here, so these are the min and max over the whole (rho, T) grid + // rather than the endpoints of an axis. They are cached at load time + // because scanning the table on each call would be O(numRho*numT). + // No sieMin()/sieMax() aliases are provided here; those names are + // superseded on the tables that already have them. + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { return sie_min_; } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { return sie_max_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return PMin_; } @@ -278,6 +289,7 @@ class SpinerEOSDependsRhoT : public EosBase { hid_t coldGroup, hid_t mfGroup); inline void fixBulkModulus_(); inline void setlTColdCrit_(); + inline void setEnergyBounds_(); PORTABLE_FORCEINLINE_FUNCTION Real lRho_(const Real rho) const noexcept { @@ -356,6 +368,8 @@ class SpinerEOSDependsRhoT : public EosBase { int numRho_, numT_; Real lRhoMin_, lRhoMax_, rhoMax_; Real lTMin_, lTMax_, TMax_; + // Defaults match the permissive EosBase bounds until a table is loaded + Real sie_min_ = -BIG_FINITE_BOUND, sie_max_ = BIG_FINITE_BOUND; Real PMin_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; @@ -694,9 +708,45 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, Real dPdR = dPdRho_.interpToReal(lRhoNormal, lTNormal); dVdTNormal_ = dPdENormal_ * CvNormal_ / (rhoNormal_ * rhoNormal_ * dPdR); + setEnergyBounds_(); + return status; } +// Energy is a dependent variable for this table, so its bounds are the +// extrema of the tabulated sie field rather than the endpoints of an +// axis. +// +// The cold curve is unioned in because it is stored separately from the +// main (rho, T) grid, so a scan of sie_ alone does not see it. In +// practice the two nearly coincide: EOSPAC's cold curve tables +// (EOS_Pc_D/EOS_Uc_D, see io_eospac.cpp:eosColdCurves) use the lowest +// temperature isotherm, and sesame2spiner defaults the grid's TMin to +// that same isotherm nudged up by an epsilon (TinyShift in +// generate_files.cpp, needed to avoid EOSPAC extrapolation errors at the +// table edge). Likewise, the from-EOS constructor evaluates the cold +// curve at exactly the TMin used for row i == 0. So this union is +// normally a tie or an epsilon-scale correction, not a large one. +// +// It is kept because the near-coincidence is a property of how tables +// happen to be generated, not a guarantee: an input deck may set Tmin +// explicitly, or shrinklTBounds may pull the grid's bottom row up away +// from the cold curve. And the cold curve is a reachable *output*, not +// merely stored data -- TableStatus::OffBottom returns sieCold_ from +// sieFromlRhoTlT_, and MinInternalEnergyFromDensity returns it directly. +// Taking the min costs one O(numRho) scan at load time and guarantees +// MinimumInternalEnergy() never reports a floor above an energy this EOS +// will actually hand back, which would break the root-find bracketing +// these bounds exist to provide. +inline void SpinerEOSDependsRhoT::setEnergyBounds_() { + sie_min_ = sie_.min(); + sie_max_ = sie_.max(); + if (sieCold_.size() > 0) { + sie_min_ = std::min(sie_min_, sieCold_.min()); + sie_max_ = std::max(sie_max_, sieCold_.max()); + } +} + inline void SpinerEOSDependsRhoT::fixBulkModulus_() { // assumes all databoxes are the same size // TODO: do we need to smooth this data with a median filter @@ -1651,6 +1701,8 @@ inline SpinerEOSDependsRhoT::SpinerEOSDependsRhoT(const EOS &source_eos, dPdENormal_ = dPdE_.interpToReal(lRhoNormal, lTNormal); Real dPdR = dPdRho_.interpToReal(lRhoNormal, lTNormal); dVdTNormal_ = robust::ratio(dPdENormal_ * CvNormal_, rhoNormal_ * rhoNormal_ * dPdR); + + setEnergyBounds_(); } } // namespace singularity diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index b7c2e846483..1ccf7925bec 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// © 2021-2026. 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 @@ -12,6 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// This file was created in part with generative AI + #ifndef _SINGULARITY_EOS_EOS_EOS_STELLAR_COLLAPSE_HPP_ #define _SINGULARITY_EOS_EOS_EOS_STELLAR_COLLAPSE_HPP_ #include @@ -236,6 +238,8 @@ class StellarCollapse : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { return sieMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { return sieMax(); } constexpr static inline int nlambda() noexcept { return _n_lambda; } template static inline constexpr bool NeedsLambda() { diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index f966dc4d92d..a9c58436b6e 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -488,6 +488,18 @@ class Variant { return PortsOfCall::visit([](const auto &eos) { return eos.MaximumDensity(); }, eos_); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumInternalEnergy() const { + return PortsOfCall::visit([](const auto &eos) { return eos.MinimumInternalEnergy(); }, + eos_); + } + + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumInternalEnergy() const { + return PortsOfCall::visit([](const auto &eos) { return eos.MaximumInternalEnergy(); }, + eos_); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return PortsOfCall::visit([](const auto &eos) { return eos.MinimumPressure(); }, diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index ecda94839cd..07d5cdf807e 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -136,9 +136,11 @@ class Vinet : public EosBase { // Essentially unbounded... I think. PORTABLE_FORCEINLINE_FUNCTION - Real MinimumPressure() const { return -1e100; } + Real MinimumPressure() const { return -BIG_FINITE_BOUND; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { + return BIG_FINITE_BOUND; + } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 08a0b3045b7..0897bbac47b 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -292,6 +292,13 @@ class UnitSystem : public EosBase> { return inv_rho_unit_ * t_.RhoPmin(temp * temp_unit_); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { + return inv_sie_unit_ * t_.MinimumInternalEnergy(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { + return inv_sie_unit_ * t_.MaximumInternalEnergy(); + } + template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp index 81e06ddb2df..9e2066e90b0 100644 --- a/singularity-eos/eos/modifiers/floored_energy.hpp +++ b/singularity-eos/eos/modifiers/floored_energy.hpp @@ -411,6 +411,10 @@ class FlooredEnergy : public EosBase> { SG_ADD_MODIFIER_METHODS(T, t_); SG_ADD_MODIFIER_MEAN_METHODS(t_); SG_ADD_MODIFIER_INTROSPECTION_METHODS(t_); + // The floor clamps energy to the per-density cold curve, which lies + // inside the underlying energy range, so the global bounds are + // unchanged. + SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS(t_); private: T t_; diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 78c0648a2e3..9dad18dfbe4 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -268,6 +268,13 @@ class BilinearRampEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real RhoPmin(const Real temp) const { return t_.RhoPmin(temp); } + // The ramp modifies pressure only, so the energy bounds pass through. + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { + return t_.MinimumInternalEnergy(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { + return t_.MaximumInternalEnergy(); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 329162c9029..282bed3e5e2 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -12,6 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// This file was created in part with generative AI + #ifndef _SINGULARITY_EOS_EOS_RELATIVISTIC_EOS_ #define _SINGULARITY_EOS_EOS_RELATIVISTIC_EOS_ @@ -208,6 +210,9 @@ class RelativisticEOS : public EosBase> { SG_ADD_MODIFIER_METHODS(T, t_); SG_ADD_MODIFIER_MEAN_METHODS(t_); SG_ADD_MODIFIER_INTROSPECTION_METHODS(t_); + // Specific internal energy passes through unmodified, so the energy + // bounds forward verbatim. + SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS(t_); private: T t_; diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index c08ddcee616..50d59f8d79b 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -57,11 +57,15 @@ class ScaledEOS : public EosBase> { ScaledEOS() = default; PORTABLE_INLINE_FUNCTION void CheckParams() const { - PORTABLE_ALWAYS_REQUIRE(std::abs(scale_) > 0, "Scale must not be zero."); - PORTABLE_ALWAYS_REQUIRE(std::abs(inv_scale_) > 0, "Inverse scale must not be zero."); PORTABLE_ALWAYS_REQUIRE(!std::isnan(scale_), "Scale must be well defined."); PORTABLE_ALWAYS_REQUIRE(!std::isnan(inv_scale_), "Inverse scale must be well defined."); + // A negative scale is not physically meaningful: it would invert the + // sign of energy and entropy, and would turn every minimum bound into + // a maximum. Require strict positivity, as UnitSystem does for its + // unit factors. This also subsumes the zero check. + PORTABLE_ALWAYS_REQUIRE(scale_ > 0, "Scale must be positive."); + PORTABLE_ALWAYS_REQUIRE(inv_scale_ > 0, "Inverse scale must be positive."); t_.CheckParams(); } auto GetOnDevice() { return ScaledEOS(t_.GetOnDevice(), scale_); } @@ -268,6 +272,15 @@ class ScaledEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real RhoPmin(const Real temp) const { return inv_scale_ * t_.RhoPmin(temp); } + // The modified energy is scale_ times the base energy. CheckParams + // requires scale_ > 0, so this preserves the ordering of the bounds. + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { + return scale_ * t_.MinimumInternalEnergy(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { + return scale_ * t_.MaximumInternalEnergy(); + } + PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { return inv_scale_ * t_.MeanAtomicMass(); } PORTABLE_INLINE_FUNCTION diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index f87e44f080e..20071eae733 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -428,6 +428,16 @@ class ShiftedEOS : public EosBase> { SG_ADD_MODIFIER_MEAN_METHODS(t_); SG_ADD_MODIFIER_INTROSPECTION_METHODS(t_); + // The modified energy is the base energy plus the shift, so the + // energy bounds shift with it. These cannot use + // SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS, which forwards verbatim. + PORTABLE_FORCEINLINE_FUNCTION Real MinimumInternalEnergy() const { + return t_.MinimumInternalEnergy() + shift_; + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumInternalEnergy() const { + return t_.MaximumInternalEnergy() + shift_; + } + private: T t_; double shift_; diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index 827f3d9145f..7114f9bf2c4 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -12,6 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// This file was created in part with generative AI + #ifndef _SINGULARITY_EOS_EOS_ZSPLIT_EOS_ #define _SINGULARITY_EOS_EOS_ZSPLIT_EOS_ @@ -281,6 +283,10 @@ class ZSplit : public EosBase> { SG_ADD_MODIFIER_METHODS(T, t_); SG_ADD_MODIFIER_MEAN_METHODS(t_); SG_ADD_MODIFIER_INTROSPECTION_METHODS(t_); + // ZSplit scales energy by a lambda-dependent factor, but the bounds + // introspection API takes no lambda. The bounds reported here are + // therefore the *un-split* bounds of the underlying EOS. + SG_ADD_MODIFIER_ENERGY_BOUNDS_METHODS(t_); private: template diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp index 6eae711e4b3..a703905251d 100644 --- a/test/test_eos_modifiers.cpp +++ b/test/test_eos_modifiers.cpp @@ -53,6 +53,19 @@ using singularity::UnitSystem; /* A toy version of ideal gas where bounds have been placed on it so we can check these are passed through modifiers properly. + + The bound values below are arbitrary sentinels, not physically + derived: the tests that use them only check that modifiers transform + a bound correctly (scale it, shift it, convert its units), and never + evaluate the EOS at one. They are deliberately given distinct + magnitudes so that a modifier forwarding the *wrong* bound -- say, an + energy bound where it meant a temperature bound -- shows up as a test + failure rather than a coincidental pass. + + Note in particular that the energy bounds are not consistent with + sie = Cv * T for the Cv the tests construct this with. They do not + need to be, and tying them to the temperature bounds would defeat the + distinctness described above. */ class BoundedGas : public IdealGas { public: @@ -80,6 +93,12 @@ class BoundedGas : public IdealGas { PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real /*temp*/) const { return MinimumDensity(); } + + // Sentinel energy bounds; see the note on the class above. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumInternalEnergy() const { return 1e-4; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumInternalEnergy() const { return 1e12; } }; #ifndef SINGULARITY_BUILD_CLOSURE @@ -331,6 +350,8 @@ SCENARIO("Modifiers propagate introspection bounds correctly", "[Modifiers]") { const Real base_min_pres = bg.MinimumPressure(); const Real base_max_pres = bg.MaximumPressureAtTemperature(0.0); const Real base_rho_pmin = bg.RhoPmin(0.0); + const Real base_min_sie = bg.MinimumInternalEnergy(); + const Real base_max_sie = bg.MaximumInternalEnergy(); AND_GIVEN("A shifted, scaled EOS") { auto eos = ScaledEOS>( @@ -359,6 +380,34 @@ SCENARIO("Modifiers propagate introspection bounds correctly", "[Modifiers]") { THEN("RhoPmin returns the scaled minimum density") { REQUIRE(isClose(eos.RhoPmin(0.0), base_rho_pmin / scale, 1.e-12)); } + + THEN("The energy bounds are shifted and then scaled") { + REQUIRE( + isClose(eos.MinimumInternalEnergy(), scale * (base_min_sie + shift), 1.e-12)); + REQUIRE( + isClose(eos.MaximumInternalEnergy(), scale * (base_max_sie + shift), 1.e-12)); + } + } + + AND_GIVEN("A scaled EOS") { + auto eos = ScaledEOS(BoundedGas(gm1, Cv), scale); + + // ScaledEOS::CheckParams requires a positive scale, so the bounds + // cannot be inverted by a sign flip and stay ordered by construction. + THEN("The energy bounds are scaled and remain ordered") { + REQUIRE(isClose(eos.MinimumInternalEnergy(), scale * base_min_sie, 1.e-12)); + REQUIRE(isClose(eos.MaximumInternalEnergy(), scale * base_max_sie, 1.e-12)); + REQUIRE(eos.MinimumInternalEnergy() < eos.MaximumInternalEnergy()); + } + } + + AND_GIVEN("A relativistic EOS") { + auto eos = RelativisticEOS(BoundedGas(gm1, Cv), 1.0); + + THEN("The energy bounds pass through untouched") { + REQUIRE(isClose(eos.MinimumInternalEnergy(), base_min_sie, 1.e-12)); + REQUIRE(isClose(eos.MaximumInternalEnergy(), base_max_sie, 1.e-12)); + } } AND_GIVEN("A UnitSystem") { @@ -378,6 +427,53 @@ SCENARIO("Modifiers propagate introspection bounds correctly", "[Modifiers]") { EPS)); REQUIRE(isClose(us.RhoPmin(0.0), base_rho_pmin / rho_unit, EPS)); } + + THEN("The unit system converts the energy bounds") { + REQUIRE(isClose(us.MinimumInternalEnergy(), base_min_sie / sie_unit, EPS)); + REQUIRE(isClose(us.MaximumInternalEnergy(), base_max_sie / sie_unit, EPS)); + AND_THEN("Multiplying by the energy unit recovers the base bounds") { + REQUIRE(isClose(us.MinimumInternalEnergy() * sie_unit, base_min_sie, EPS)); + REQUIRE(isClose(us.MaximumInternalEnergy() * sie_unit, base_max_sie, EPS)); + } + } + } + } + + GIVEN("An unmodified analytic EOS") { + IdealGas ig(0.5, 2.0); + + THEN("The energy bounds are the permissive defaults") { + // The minimum must be negative. Energies are legitimately negative + // for cold curves and for shifted EOS, so zero is not a safe floor. + REQUIRE(ig.MinimumInternalEnergy() < 0); + REQUIRE(ig.MaximumInternalEnergy() > 0); + REQUIRE(ig.MinimumInternalEnergy() < ig.MaximumInternalEnergy()); + } + } +} + +SCENARIO("Energy bounds are reachable through the EOS variant", "[Modifiers][Variant]") { + GIVEN("A BoundedGas in a unit system, held in a variant") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + constexpr Real rho_unit = 2.0; + constexpr Real sie_unit = 3.0; + constexpr Real temp_unit = 4.0; + constexpr Real EPS = 10 * singularity::robust::EPS(); + + BoundedGas bg(gm1, Cv); + const Real base_min_sie = bg.MinimumInternalEnergy(); + const Real base_max_sie = bg.MaximumInternalEnergy(); + + // This is the motivating use case: the host code holds a + // runtime-polymorphic EOS and needs the energy bounds in its own + // unit system, without stripping the modifier. + using EOS = singularity::Variant>; + EOS eos = UnitSystem(BoundedGas(gm1, Cv), rho_unit, sie_unit, temp_unit); + + THEN("The variant reports the converted energy bounds") { + REQUIRE(isClose(eos.MinimumInternalEnergy(), base_min_sie / sie_unit, EPS)); + REQUIRE(isClose(eos.MaximumInternalEnergy(), base_max_sie / sie_unit, EPS)); } } }