Fix for Carnahan-Starling EOS with DensityFromPressureTemperature #563

bgclayton-lanl · 2025-09-10T04:53:51Z

PR Summary

This fix corrects some of the fails with the analytic EOS tests regarding the Carnahan-Starling. Mainly it swaps our the root solver for the regular_falsi method. I haven't gone through the analysis to determine if the root is unique though. Maybe someone else has done this analysis before?

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

Yurlungur

Thanks for the fix! This looks like this is introducing test failures, but that might just be formatting changes. Can you run

make format_singularity-eos

in the build directory? Hopefully that's why tests are failing. But we might need to make other changes to see. After tests pass on github, we'll also need to make sure they pass on LANL systems by triggering tests on re-git.

You should also add a changelog entry.

singularity-eos/eos/eos_carnahan_starling.hpp

bgclayton-lanl · 2025-09-10T05:33:38Z

Thanks for the fix! This looks like this is introducing test failures, but that might just be formatting changes. Can you run
make format_singularity-eos
in the build directory? Hopefully that's why tests are failing. But we might need to make other changes to see. After tests pass on github, we'll also need to make sure they pass on LANL systems by triggering tests on re-git.

You should also add a changelog entry.

I currently don't have a way to run the clang formatting thing. I need to talk with Jeff more about this and see the best way to set this up. I also don't have access to run the re-git process either. But I'll go ahead and update the Changelog. Thank you.

Yurlungur · 2025-09-10T05:36:38Z

Ok Jeff or I can run the formatting and trigger the tests

jhp-lanl

Also, the gold standards for the test need to be updated (see test). Ideally we'd do something where we didn't need gold standards, but that's probably beyond the scope of this PR.

-------------------------------------------------------------------------------
Scenario: CarnahanStarling6
      Given: Parameters for a CarnahanStarling EOS
      Given: Pressure and temperature
       When: A rho(P, T) lookup is performed
       Then: The returned rho(P, T) should be equal to the true value
-------------------------------------------------------------------------------
/home/runner/work/singularity-eos/singularity-eos/test/test_eos_carnahan_starling.cpp:686
...............................................................................

/home/runner/work/singularity-eos/singularity-eos/test/eos_unit_test_helpers.hpp:82: FAILED:
  CHECK( isClose(z[i], ztrue[i], tol) )
with expansion:
  false
with message:
  i: 1, Pressure: 3.44505000000000000e+07, Temperature: 1.53209999999999991e+
  03, Value: 7.82907368899728694e-03, True Value: 7.82907368903815008e-03

/home/runner/work/singularity-eos/singularity-eos/test/eos_unit_test_helpers.hpp:82: FAILED:
  CHECK( isClose(z[i], ztrue[i], tol) )
with expansion:
  false
with message:
  i: 2, Pressure: 6.78877500000000000e+07, Temperature: 2.76605000000000018e+
  03, Value: 8.54539433273353993e-03, True Value: 8.54539433278823403e-03

===============================================================================

singularity-eos/eos/eos_carnahan_starling.hpp

…d the CHANGELOG

bgclayton-lanl · 2025-09-10T19:30:02Z

@Yurlungur @jhp-lanl I think I've gotten everything fixed. Let me know if I need to change anything else.

Yurlungur · 2025-09-11T13:19:55Z

Also, the gold standards for the test need to be updated (see test). Ideally we'd do something where we didn't need gold standards, but that's probably beyond the scope of this PR.

-------------------------------------------------------------------------------
Scenario: CarnahanStarling6
      Given: Parameters for a CarnahanStarling EOS
      Given: Pressure and temperature
       When: A rho(P, T) lookup is performed
       Then: The returned rho(P, T) should be equal to the true value
-------------------------------------------------------------------------------
/home/runner/work/singularity-eos/singularity-eos/test/test_eos_carnahan_starling.cpp:686
...............................................................................

/home/runner/work/singularity-eos/singularity-eos/test/eos_unit_test_helpers.hpp:82: FAILED:
  CHECK( isClose(z[i], ztrue[i], tol) )
with expansion:
  false
with message:
  i: 1, Pressure: 3.44505000000000000e+07, Temperature: 1.53209999999999991e+
  03, Value: 7.82907368899728694e-03, True Value: 7.82907368903815008e-03

/home/runner/work/singularity-eos/singularity-eos/test/eos_unit_test_helpers.hpp:82: FAILED:
  CHECK( isClose(z[i], ztrue[i], tol) )
with expansion:
  false
with message:
  i: 2, Pressure: 6.78877500000000000e+07, Temperature: 2.76605000000000018e+
  03, Value: 8.54539433273353993e-03, True Value: 8.54539433278823403e-03

===============================================================================

It doesn't look like this issue has been resolved. You can either change the tolerance in the unit test or change the value it finds.

Yurlungur · 2025-09-11T13:25:59Z

I tightened the root finder tolerance on regula falsi and ran formatting. Lets see if that resolves the failed tests. I also triggered the CI on re-git.

Yurlungur

Looks good. Thanks for the fix! Will only merge after @jhp-lanl also approves and tests pass.

Yurlungur · 2025-09-11T13:29:31Z

Ok those too-tight tolerances cause the solver to fail. I just switched to a lower tolerance on the test, as clearly there's some jitter in the root that the finder finds across methods and other choices.

Yurlungur · 2025-09-11T14:19:17Z

Ah I think I found an underlying problem, which may have been the issue the hole time. If the upper bound of the root find depends on 1/bb, we cannot allow bb to be exactly zero. I modified the code to enforce strictly bb > 0 and modified the tests appropriately.

jhp-lanl · 2025-09-11T14:30:23Z

Ah I think I found an underlying problem, which may have been the issue the hole time. If the upper bound of the root find depends on 1/bb, we cannot allow bb to be exactly zero. I modified the code to enforce strictly bb > 0 and modified the tests appropriately.

Well, but bb = 0 corresponds to an ideal gas, which I think is what the test was checking. I'd prefer to be able to allow that condition and just choose a better upper bound on the root find. I'll make suggestions in the code

jhp-lanl

I would add an early return for _bb == 0. Well, technically I suggested also taking into account very small values of _bb, but either way the spirit is the same.

singularity-eos/eos/eos_carnahan_starling.hpp

jhp-lanl · 2025-09-11T14:39:14Z

singularity-eos/eos/eos_carnahan_starling.hpp

-    auto status = findRoot(poly, press, guess, lo_bound, hi_bound, xtol, ytol, real_root);
-    if (status != Status::SUCCESS) {
+    static constexpr Real rho_low = 0.;
+    const Real rho_high = robust::ratio(1.0, _bb);


If _bb == 0 then we can just use the ideal gas inversion

Suggested change

const Real rho_high = robust::ratio(1.0, _bb);

// Use ideal gas value when _bb is zero

if (_bb <= std::numeric_limits<Real>::min()) {

return press / (_Cv * _gm1 * temperature);

}

const Real rho_high = 1.0 / _bb;

Note that I got rid of the robust::ratio now that we have an explicit guarantee that _bb > 0 at this point in the code.

I don't think we should rely on min. I think this should be a number slightly larger than the root finder tolerance.

Ah okay, yeah that's a good suggestion

Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com>

Yurlungur · 2025-09-11T15:22:19Z

I am fine with that change but am not going to be able to debug this further today. You guys can push changes. Something is concerning with this system still IMO. Tests are passing locally for me but not in the CI. This suggests some low-level floating point issues that have not been tracked down.

bgclayton-lanl · 2025-09-11T15:26:13Z

test/test_eos_carnahan_starling.cpp

          }
          array_compare(num, density, energy, h_pressure, pressure_true, "Density",
-                        "Energy");
+                        "Energy", 1e-5);


Is the 1e-5 tolerance fine? It seems kind of big.

It is kind of large, but it is about as large as I was seeing differences. If you run the tests locally yourself you can experiment and maybe find a tighter bound... or some other underlying reason why the tests are failing.

You can also update the gold values to the new values and keep a tighter tolerance

Yurlungur · 2025-09-11T15:48:52Z

singularity-eos/eos/eos_carnahan_starling.hpp

+      const Real term1 = x * (1. + eta + eta * eta - eta * eta * eta);
+      const Real term2 =
+          press / (_Cv * _gm1 * temperature) * math_utils::pow<3>(1. - eta);


Is this numerically safe? We presumably want this to return something sensible for zero temperature but I think right now it will give NaN or infinity. You can and should replace the / by robust::ratio but you may also want to special case the low pressure limit.

That's a good point. So if we allow $T = 0$, this actually implies that $e = q$ and hence $P = 0$. So I think it's reasonable to assume that the $\rho = 0$ (vacuum state). It looks like in general the density can be any value at zero temperature. So we may want to check that $T > 0$ and if not return $\rho = 0$.

bgclayton-lanl · 2025-09-13T04:31:25Z

Okay, sorry for all the issues. It started opening of the door to other bigger questions. The more fundamental question is how to define convergence criteria. But in any case, all the tests are passing, and I updated the tests to include a unit scaling (based on the reference state) on the epsilon for checking errors. Let me know if there are any other changes to do.

@jhp-lanl @Yurlungur

Yurlungur

This seems entirely sensible. Thanks for the careful digging and nice work @bgclayton-lanl !

@jhp-lanl or I will run formatting and trigger tests on re-git.

jhp-lanl

One of the jobs failed, but I'm rerunning it

jhp-lanl · 2025-09-15T17:22:04Z

All tests pass now

bgclayton-lanl added bug Something isn't working Robustness Ensures that existing features work as intended labels Sep 10, 2025

Yurlungur reviewed Sep 10, 2025

View reviewed changes

singularity-eos/eos/eos_carnahan_starling.hpp Show resolved Hide resolved

jhp-lanl requested changes Sep 10, 2025

View reviewed changes

singularity-eos/eos/eos_carnahan_starling.hpp Outdated Show resolved Hide resolved

singularity-eos/eos/eos_carnahan_starling.hpp Outdated Show resolved Hide resolved

singularity-eos/eos/eos_carnahan_starling.hpp Outdated Show resolved Hide resolved

bgclayton-lanl added 2 commits September 10, 2025 12:22

Updated the DensityFromPressureTemperature function for Carn-Star EOS

0334fd8

Cleaned up DensityFromPressureTemperature, added comments, and update…

9c04acc

…d the CHANGELOG

bgclayton-lanl force-pushed the bgclayton_carnahan_star_fix branch from 7d00fca to 9c04acc Compare September 10, 2025 18:24

Didn't mean to touch the CMakeLists

0fbf049

Yurlungur linked an issue Sep 11, 2025 that may be closed by this pull request

[rho,e](P, T) failing for Carnahan Starling on some systems #452

Closed

formatting

ac05bed

Tighten tolerance on root finder

a9e28ee

Yurlungur approved these changes Sep 11, 2025

View reviewed changes

loosen tolerances on test. return to 1e-12 on regula falsi

18c1075

jonahm-LANL added 2 commits September 11, 2025 10:08

Enforce bb > 0 strictly

fb19bf9

Do not set bb to too small of a number

3ab2ceb

jhp-lanl requested changes Sep 11, 2025

View reviewed changes

Apply suggestions from code review

f0938a7

Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com>

bgclayton-lanl commented Sep 11, 2025

View reviewed changes

Yurlungur reviewed Sep 11, 2025

View reviewed changes

bgclayton-lanl added 2 commits September 11, 2025 13:47

Updated the error tolerences with a reference value.

65df56d

Cleaned up the testing some more.

afb51ec

bgclayton-lanl added 5 commits September 11, 2025 15:53

Reduced a tolerance for one test.

23269a3

Reduced the tolerances some.

bf33260

Fixed a ill-defined recursive definition of

46fdd21

Fixed another bug with the use of .

cf6182a

Reduced the tolerance since the CI was failing.

b1ebf8f

Yurlungur approved these changes Sep 13, 2025

View reviewed changes

format tests

c41fb1a

jhp-lanl approved these changes Sep 15, 2025

View reviewed changes

jhp-lanl merged commit 019405e into main Sep 15, 2025
9 of 10 checks passed

jhp-lanl deleted the bgclayton_carnahan_star_fix branch September 15, 2025 17:22

-    const Real rho_high = robust::ratio(1.0, _bb);
+    // Use ideal gas value when _bb is zero
+    if (_bb <= std::numeric_limits<Real>::min()) {
+      return press / (_Cv * _gm1 * temperature);
+    }
+    const Real rho_high = 1.0 / _bb;

Fix for Carnahan-Starling EOS with DensityFromPressureTemperature #563

Fix for Carnahan-Starling EOS with DensityFromPressureTemperature #563

Uh oh!

Conversation

bgclayton-lanl commented Sep 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR Summary

PR Checklist

Uh oh!

Yurlungur left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

bgclayton-lanl commented Sep 10, 2025

Uh oh!

Yurlungur commented Sep 10, 2025

Uh oh!

jhp-lanl left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

bgclayton-lanl commented Sep 10, 2025

Uh oh!

Yurlungur commented Sep 11, 2025

Uh oh!

Yurlungur commented Sep 11, 2025

Uh oh!

Yurlungur left a comment

Choose a reason for hiding this comment

Uh oh!

Yurlungur commented Sep 11, 2025

Uh oh!

Yurlungur commented Sep 11, 2025

Uh oh!

jhp-lanl commented Sep 11, 2025

Uh oh!

jhp-lanl left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Yurlungur commented Sep 11, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bgclayton-lanl commented Sep 13, 2025

Uh oh!

Yurlungur left a comment

Choose a reason for hiding this comment

Uh oh!

jhp-lanl left a comment

Choose a reason for hiding this comment

Uh oh!

jhp-lanl commented Sep 15, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

bgclayton-lanl commented Sep 10, 2025 •

edited

Loading