Make sure that DensityEnergyFromPressureTemperature works for all equations of state #449

Yurlungur · 2024-12-24T23:20:49Z

PR Summary

The easiest way to get the Jacobian for the PT-space PTE solver we are building is to ensure that the DensityEnergyFromPressureTemperature call works consistently for all equations of state. This MR makes sure thats the case and threads this call through all the tests.

To ease this process, I added a base class implementation. However, I think most classes will want their own implementations, and I added them when appropriate.

PR Checklist

Adds a test for any bugs fixed. Adds tests for new features.
Format your changes by using the make format command after configuring with cmake.
Document any new features, update documentation for changes made.
Make sure the copyright notice on any files you modified is up to date.
After creating a pull request, note it in the CHANGELOG.md file.
LANL employees: make sure tests pass both on the github CI and on the Darwin CI

If preparing for a new release, in addition please check the following:

Update the version in cmake.
Move the changes in the CHANGELOG.md file under a new header for the new release, and reset the categories.
Ensure that any when='@main' dependencies are updated to the release version in the package.py

… thats the issue

Yurlungur · 2024-12-25T00:14:59Z

.github/workflows/tests_minimal.yml

@@ -36,4 +36,4 @@ jobs:
                  ..
            make -j4
            make install
-            make test
+            ctest --output-on-failure


I think I'd like to leave this permanently. You still get a summary out at the end, but here we can see what went wrong with a test. Useful for debugging, especially when the CI machine and your laptop disagree...

Yurlungur · 2024-12-25T00:16:00Z

singularity-eos/eos/eos_base.hpp

+  // JMM: This method is often going to be overloaded for special cases.
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
+                                       Indexer_t &&lambda, Real &rho, Real &sie) const {


This default implementation is usually not optimal but does seem to work if a given EOS doesn't need special cases. In general, though, an EOS should overload and implement its own.

Yurlungur · 2024-12-25T00:17:28Z

singularity-eos/eos/eos_base.hpp

+    // JMM: This needs to not fail and instead return something sane
+    if (status != Status::SUCCESS) {
+      PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n");
+      rho = rhoguess;
+    }


This is a function that may be passed data and I see it as similar to a PTE solve. Even if it fails out, something reasonable must be done.

singularity-eos/eos/eos_helmholtz.hpp

Yurlungur · 2024-12-25T00:19:24Z

singularity-eos/eos/eos_mgusup.hpp

+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const {
+    if (_s > 1) {
+      return 0.99 * robust::ratio(_s * _rho0, _s - 1);
+    } else { // for s <= 1, no maximum, but we need to pick something.
+      return 1e3 * _rho0;
+    }
+  }


These are the conditions where the reference isotherm is a positive temperature.

Honestly the $s < 0$ case implies that the fundamental shock derivative is negative and the EOS is no longer convex, so it's not a very useful EOS. I think what you picked is as good as any for that case.

👍 I will add a comment to this effect.

Yurlungur · 2024-12-25T00:19:59Z

singularity-eos/eos/eos_mgusup.hpp

+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumDensity() const { return 10 * robust::EPS(); }


Since the reference isotherm scales as _rho0/rho, the EOS diverges at zero density. This bounds it at least to some degree.

Yurlungur · 2024-12-25T00:21:58Z

singularity-eos/eos/eos_stellar_collapse.hpp

+  auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); };
+  auto status = regula_falsi(lPofRT, lP, lrguess, lRhoMin_, lRhoMax_, ROOT_THRESH,
+                             ROOT_THRESH, lrguess);


Since the EOS is tabulated in logrho, use logrho for the root find.

Yurlungur · 2024-12-25T00:23:08Z

singularity-eos/eos/eos_vinet.hpp

@@ -131,6 +131,10 @@ class Vinet : public EosBase<Vinet> {
  ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
                         Real &bmod, Real &dpde, Real &dvdt,
                         Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; }


EOS depends on $\left(\frac{\rho_0}{\rho}\right)^{1/3}$ so a reasonable bound before the EOS will become ill behaved is this.

Yurlungur · 2024-12-25T00:23:31Z

singularity-eos/eos/eos_vinet.hpp

-  // Density/Energy from P/T not unique, if used will give error
-  template <typename Indexer_t = Real *>
-  PORTABLE_INLINE_FUNCTION void
-  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                       Indexer_t &&lambda, Real &rho, Real &sie) const;


use the default in the base class

Yurlungur · 2024-12-25T00:24:22Z

test/eos_unit_test_helpers.hpp

+template <typename EOS, typename Indexer_t = Real *>
+PORTABLE_INLINE_FUNCTION bool
+CheckRhoSieFromPT(EOS eos, Real rho, Real T,
+                  Indexer_t &&lambda = static_cast<Real *>(nullptr)) {


Convenience function to let me test the functionality more easily

Yurlungur · 2024-12-25T00:34:02Z

Ok that was finicky but tests on github now seem to be passing. re-git tests appear to be having some network issues. Ping @rbberger for after the holidays. :)

singularity-eos/eos/eos_base.hpp

singularity-eos/eos/eos_helmholtz.hpp

jhp-lanl · 2025-01-03T00:39:20Z

singularity-eos/eos/eos_mgusup.hpp

+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const {
+    if (_s > 1) {
+      return 0.99 * robust::ratio(_s * _rho0, _s - 1);
+    } else { // for s <= 1, no maximum, but we need to pick something.
+      return 1e3 * _rho0;
+    }
+  }


Honestly the $s < 0$ case implies that the fundamental shock derivative is negative and the EOS is no longer convex, so it's not a very useful EOS. I think what you picked is as good as any for that case.

singularity-eos/eos/eos_mgusup.hpp

… anyway and I think they might be interacting badly with GPUs

Yurlungur · 2025-01-04T00:31:53Z

OK... that took longer than I thought it would. All tests now passing once again and the machinery is more flexible and robust. @jhp-lanl please re-review when you get the chance. I believe I've addressed all your concerns.

I also noticed that the DensityEnergyFromPressureTemperature function wasn't documented. It's documented now.

Yurlungur · 2025-01-08T17:52:55Z

Ping @jhp-lanl

jhp-lanl · 2025-01-08T19:30:18Z

Ping @jhp-lanl

I'll try to re-review that by the end of today.

Yurlungur · 2025-01-08T21:58:45Z

Snuck in introspection for pressure into this MR as well.

jhp-lanl

Approving because none of my comments are really blocking. That said, I do think the name change to MaximumPressureAtTemperature() would be preferable. I also think the documentation should be changed to indicate that there could be a temperature dependence for the maximum pressure if an EOS supports it.

jhp-lanl · 2025-01-09T03:46:10Z

doc/sphinx/src/using-eos.rst

+.. cpp:function:: Real MaximumPressureFromTemperature() const;
+
+provides a maximum possible pressure an equation of state
+supports. (Most models are unbounded in pressure.) This is again
+useful for root finds.


If you're modeling this off of what I did in the Gruneisen EOS, this isn't correct. This gives the maximum pressure at a given temperature. If you increase the temperature, the maximum pressure will increase.

For a sesame table, this would simply be $P(\rho_\mathrm{max}, T)$ assuming that $\rho_\mathrm{max}$ is not inside the vapor dome for some strange reason.

I know our convention is to use "From" in our lookup functions, but I'd argue this is subtly different and should be "At" instead. I'd argue our lookups are saying "Return a quantity given the provided state" ("Given" would be more accurate in the naming, but I'm not proposing it here). By contrast, this function does not take a complete thermodynamic state, instead taking a single variable and returning the information to form a complete state. Temperature by itself is not a complete thermodynamic state.

I completely agree--I change the name. For the tabulated equations of state, I actually do pull it out and record it now (thinking ahead), but I let this float as the table interpolation does extrapolate off the top of the table. It's not particularly trustworthy, but it's not a hard limit. Only Gruneisen currently reports a maximum pressure bound. On the other hand, I do have each table report the minimum pressure, which is important.

doc/sphinx/src/using-eos.rst

singularity-eos/eos/eos_eospac.hpp

singularity-eos/eos/eos_variant.hpp

Yurlungur · 2025-01-09T16:48:07Z

Thanks for the review, @jhp-lanl ! Agree with your comments and have addressed them. Will merge once tests pass on re-git.

jonahm-LANL added 6 commits December 23, 2024 12:28

add RhoE of PT to base class and implement where relevant

4a210b5

formatting

50c08e7

Tests pass on macbook. Still need to add a few tests

acd9def

swap default to findRoot because it respects bounds

608a772

tests complete for everything except eospac and spiner

0bdead9

update zsplit

0a4b61a

Yurlungur added bug Something isn't working enhancement New feature or request Testing Additions/changes to the testing infrastruture Robustness Ensures that existing features work as intended labels Dec 24, 2024

Yurlungur requested review from dholladay00, chadmeyer and jhp-lanl December 24, 2024 23:20

Yurlungur self-assigned this Dec 24, 2024

changelog

f1d1fc4

Yurlungur mentioned this pull request Dec 24, 2024

regula_falsi solver does not respect regions of validity #450

Open

jonahm-LANL added 9 commits December 24, 2024 16:24

missing auto

c5f75a5

it works on my machine...

838b720

output on failure please

365dd65

just change tests to output-on-failure... and also loosen EPS in case…

376223f

… thats the issue

lets try this...

9b083fa

make default implementation more robust maybe I hope

c3182ff

put default rho gues close to STP

1034b38

put default rho gues close to STP

d6e80ee

re-tighten bounds after initial guess problem resolved

5f79307

Yurlungur commented Dec 25, 2024

View reviewed changes

oops missing boolean not

0e61f98

jhp-lanl approved these changes Jan 3, 2025

View reviewed changes

add several comments based on Jeff's feedback

cc8db8f

jonahm-LANL added 17 commits January 3, 2025 12:27

formatting

8cc2508

namespace

6dc8cae

tweak how to check this residual to work when P = 0

bc5f10e

oops need to do magnitudes

813e773

come up with sensible bounds for power series EOSs.

e0e587c

I give up. 1e-4 it is

cf286de

comment

2cdcc08

remove unused variable

1d1873f

tighter tolerances

ceec265

relax tolerance for chekc

8424efc

remove the maximum densities... we're never going to hit those bounds…

7d9c90a

… anyway and I think they might be interacting badly with GPUs

also add sanity check for rhomax vs rhomin

ee5e374

missed leading underscore

38bba46

stupid C++ deleting copy constructor for const correctness

a74bd26

choose worse, but more generic guess... does that help?

5f8e8e0

is it because rho is uninitialized?

74ee377

update docs

9e4f844

add introspection for pressure

1cd2f9d

Yurlungur mentioned this pull request Jan 8, 2025

[WIP] pt space pte solver #453

Open

13 tasks

jhp-lanl approved these changes Jan 9, 2025

View reviewed changes

jonahm-LANL added 2 commits January 9, 2025 08:58

add default nullptr argument

6967e55

change name

144b90b

Yurlungur mentioned this pull request Jan 9, 2025

More care for EOS extrapolation for tables #454

Open

Yurlungur merged commit 3f78b83 into main Jan 9, 2025
8 checks passed

Yurlungur deleted the jmm/pt-of-re-everywhere branch January 9, 2025 17:33

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Make sure that DensityEnergyFromPressureTemperature works for all equations of state #449

Make sure that DensityEnergyFromPressureTemperature works for all equations of state #449

Yurlungur commented Dec 24, 2024 •

edited

Loading

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

jhp-lanl Jan 3, 2025

Yurlungur Jan 3, 2025

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur Dec 25, 2024

Yurlungur commented Dec 25, 2024

jhp-lanl Jan 3, 2025

Yurlungur commented Jan 4, 2025

Yurlungur commented Jan 8, 2025

jhp-lanl commented Jan 8, 2025

Yurlungur commented Jan 8, 2025

jhp-lanl left a comment

jhp-lanl Jan 9, 2025

Yurlungur Jan 9, 2025

Yurlungur commented Jan 9, 2025

		PORTABLE_FORCEINLINE_FUNCTION
		Real MinimumDensity() const { return 10 * robust::EPS(); }

Make sure that DensityEnergyFromPressureTemperature works for all equations of state #449

Make sure that DensityEnergyFromPressureTemperature works for all equations of state #449

Conversation

Yurlungur commented Dec 24, 2024 • edited Loading

PR Summary

PR Checklist

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Yurlungur commented Dec 25, 2024

Choose a reason for hiding this comment

Yurlungur commented Jan 4, 2025

Yurlungur commented Jan 8, 2025

jhp-lanl commented Jan 8, 2025

Yurlungur commented Jan 8, 2025

jhp-lanl left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Yurlungur commented Jan 9, 2025

Yurlungur commented Dec 24, 2024 •

edited

Loading