Simple MACAW EOS #560

bgclayton-lanl · 2025-09-08T21:51:34Z

PR Summary

This update implements the Simple MACAW EOS from Tariq and Eduardo: https://www.osti.gov/biblio/2479474

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

…emperature.

… more things to do.

Yurlungur · 2025-09-09T18:53:13Z

Thanks for adding this, @bgclayton-lanl ! (And nice to virtually meet you.) Sorry for the delay; @jhp-lanl and I will take a look at this. Just haven't gotten to it yet. This is a particularly busy week for me.

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

Thanks for adding this, @bgclayton-lanl ! (And nice to virtually meet you.) Sorry for the delay; @jhp-lanl and I will take a look at this. Just haven't gotten to it yet. This is a particularly busy week for me.

Hi, good to meet you. I'm sure there will be a lot of things to change or tweak. The main thing I wasn't too sure about was the unit test part (Catch2). I saw in some other unit tests that there are different things for "getting on device" and "Kokkos" and all of these things. But I just kept it basic for the time being. But I can update them to replicate the testing from other EOS.

jhp-lanl

Some small changes, but overall it looks great!

doc/sphinx/src/models.rst

jhp-lanl · 2025-09-10T22:59:41Z

doc/sphinx/src/models.rst

+In order to maintain thermodynamic consistency, a sufficient set of constraints
+is given by :math:`A > 0`, :math:`B > 0`, :math:`C^\infty_v > 0`, :math:`v_0 >
+0`, :math:`T_0 > 0`, and :math:`\Gamma_c > 0`. If one wants to maintain
+thermodynamic consistency, we also require :math:`\Gamma_c \in (0,1)`. 


Great addition! But why is there an upper bound on $\Gamma_c$?

If $\Gamma_C \geq 1$, then the Simple MACAW may not be thermodynamically consistent. Basically, the isothermal sound speed could be negative. The specific details are mentioned in the paper.

Two thoughts here after reading the paper:

The upper bound on $\Gamma$ is actually for thermodynamic stability not thermodynamic consistency. We routinely deal with EOS that are not fully stable (i.e. VdW loops).

The upper bound on Gamma is overly restrictive. It is a guarantee that the isothermal sound speed is ~~positive~~ not imaginary (the paper says $c_T$ > 0 in the text, but really it's $c_T^2 > 0$), but it's possible to achieve that over a domain of interest without this restriction or for a particular selection of $B$.

As a result, I'd like to provide our users with the ability to decide what their domain is and how to deal with this constraint. I think the best course of action would be to reproduce equation 22 from the paper and discuss that $\Gamma > 1$ could produce an imaginary thermal sound speed for large volumes and low temperatures. We can warn on $\Gamma > 1$, but ultimately let the user decide if their parameters are justified.

singularity-eos/eos/eos_simple_macaw.hpp

jhp-lanl · 2025-09-10T23:28:22Z

singularity-eos/eos/eos_simple_macaw.hpp

+    PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative");
+    PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)");


For redundancy, can you also explain the upper-bound here as well as in the documentation? It's just a bit surprising to me. A comment is fine.

If the upper bound is more of a desirable thing rather than a required thing, then we can use a PORTABLE_WARN() if the value is above 1.

Given the discussion above (and after skimming the paper), I think that while _Gc > 0 is a hard restriction, we should only warn if _Gc > 1. Please add a comment about thermodynamic stability being the impetus for the restriction.

Oh, and I should add that although they give the stability criterion as an open interval, I think the stability criterion is actually a closed interval on the high side.

singularity-eos/eos/eos_simple_macaw.hpp

jhp-lanl · 2025-09-10T23:35:04Z

singularity-eos/eos/eos_simple_macaw.hpp

+template <typename Indexer_t>
+PORTABLE_INLINE_FUNCTION void
+SimpleMACAW::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod,
+                  const unsigned long output, Indexer_t &&lambda) const {
+  if (output & thermalqs::density && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie);
+  }
+  if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::density || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
+  }
+  if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) {
+    sie = InternalEnergyFromDensityPressure(rho, press);
+  }
+  if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::temperature)
+    temp = TemperatureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::bulk_modulus)
+    bmod = BulkModulusFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::specific_heat)
+    cv = SpecificHeatFromDensityInternalEnergy(rho, sie);
+}


You can have this in-line if you'd like. I don't think there's a reason to define it separate from where you declare it in the header. @Yurlungur correct me if I'm wrong

you're correct. Only reason to separate it out is for clarity. Choice left to the author.

jhp-lanl · 2025-09-10T23:37:57Z

singularity-eos/eos/eos_simple_macaw.hpp

+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }


We have some cruft in our code unfortunately. I think we'd ideally like to update these to be PORTABLE_REQUIRE statements instead.

@Yurlungur thoughts on this versus PORTABLE_ALWAYS_REQUIRE?

FYI @bgclayton-lanl the ALWAYS essentially is for causing an abort regardless of whether it's a debug build or not. Without ALWAYS, the check is only performed in debug builds.

PORTABLE_ALWAYS_REQUIRE would work. You could also use PORTABLE_ALWAYS_THROW_OR_ABORT with a useful error message.

I'm not sure I understand here. What needs to be changed? Are you saying to do a PORTABLE_REQUIRE or PORTABLE_ALWAYS_REQUIRE on the specific thermo flag that gets passed in?

So the current structure if if (condition) {ERROR}. That's essentiallyPORTABLE_REQUIRE(!condition). The latter is I think easier to read and ensures consistent error generation/handling across platforms.

As for ALWAYS part, @Yurlungur or @jonahm-LANL can you advise on what you think is appropriate here? My philosophy is that assert statements (i.e. those that are only checked at debug) should be reserved for highly unlikely conditions the developer thinks shouldn't happen; if they trigger, then they tell the developer that they need to refactor their logic. Illegal that a user can easily set (either by mistake or on purpose) should be always checked in my opinion unless they can be shown to be a performance bottleneck. I would argue that this falls into the latter category, so I would propose changing these error conditions to PORTABLE_ALWAYS_REQUIRE statements (with the appropriate negation of course).

If @Yurlungur doesn't give his opinion though, just go with the ALWAYS version for now.

I think this should always fail if you hit it, so I think PORTABLE_ALWAYS.

singularity-eos/eos/eos_simple_macaw.hpp

jhp-lanl · 2025-09-11T00:27:01Z

test/test_eos_simple_macaw.cpp

+// Create the EOS
+auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc);
+auto eos = host_eos.GetOnDevice();
+
+// Only run exception checking when we aren't offloading to the GPUs
+SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") {
+  GIVEN("Parameters for a Simple MACAW EOS") {
+    Real rho = 0.5;
+    for (int i = 0; i < 10; i++) {
+      rho += rho + i; // cylce through a variety of densities
+      THEN("The temperature evaluated on the cold curve should produce zero") {
+        Real v = 1.0 / rho;
+        Real e = eos.SieColdCurve(v);
+        REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0);
+      }
+    }
+  }
+}


I'd just group all of this under one scenario to test the SimpleMACAW eos. Although we have our tests fairly isolated, I think it's still good practice to scope things inside the SCENARIO.

Also, there are specific ways to test in loops with Catch2. Essentially the scenario, etc. all form a leaf in a tree and each leaf needs to be unique. In essence it means that some concatenation of all the names of the sections (GIVEN, THEN, etc.) needs to be unique.

So this test would look like this:

SCENARIO("Testing the SimpleMACAW EOS", "[SimpleMACAWEOS]") { GIVEN("Parameters for a Simple MACAW EOS") { constexpr Real A = 7.3; constexpr Real B = 3.9; constexpr Real Cvinf = 0.000389; constexpr Real v0 = 1. / 8.952; constexpr Real T0 = 150.; constexpr Real Gc = 0.5; auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); auto eos = host_eos.GetOnDevice(); // Not needed since no dynamic memory WHEN("A set of densities is provided") { Real rho = 0.5; for (int i = 0; i < 10; i++) { rho += rho + i; // cycle through a variety of densities DYNAMIC_SECTION("For a given density " << rho << " and energy from the cold curve") { Real v = 1.0 / rho; Real e = eos.SieColdCurve(v); INFO("rho = " << rho << " v = " << v << " e = " << e); THEN("The temperature at this density and energy should be zero") { REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); } } } } } }

All of the GIVEN, THEN, WHEN, etc. macros all the same from a code perspective. It really just changes how the test prints out when it fails. So feel free to change out whatever you like. The DYNAMIC_SELECTION macro is just there for ease of changing the name for each iteration of the loop (needed when you declare a section inside a loop like that).

I also added the INFO macro to print out the inputs to the lookup if it fails. That can help diagnose issues.

jhp-lanl · 2025-09-11T03:54:49Z

Oh, and one more thing: can you make sure to include some test that looks something like this:

WHEN("We put the SimpleMACAW EOS into a Variant") {
  using EOS = Variant<SimpleMACAW>;
  EOS eos_in_variant = eos;
  THEN("The SimpleMACAW model conforms to the Variant API") {}
}

Alternatively, you could do this with your inversion test and just use the eos_in_variant object for the lookups.

The point here is that the test will fail to compile if you are missing any API functions in the Variant class (see eos_variant.hpp). But don't worry about the vector overloads since those are handled by the base class. You just need to worry about the scalar lookup functions (anything where the thermodynamic variables are Real or the return type is Real. If you followed a working EOS then it should work fine already.

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

doc/sphinx/src/models.rst

@@ -1320,6 +1322,58 @@ This constructor also optionally accepts `MeanAtomicProperties` for
 the atomic mass and number as a final optional parameter.


+.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256
+
+Simple MACAW


Good documentation! Nice work.

Yurlungur · 2025-09-11T13:33:37Z

singularity-eos/eos/eos_simple_macaw.hpp

+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
+      const Real rho, const Real temperature,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+      /* Equation (15) */


Nice including equation numbers

Yurlungur · 2025-09-11T13:34:37Z

singularity-eos/eos/eos_simple_macaw.hpp

+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real SieColdCurve(
+      const Real rho, Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {


Do you want these functions exposed in the public API for this EOS? If so, you can leave them here. If not, move them to private scope and add a trailing underscore to the name.

I think in general it should be. In my own research, I can imagine ways in which one may need to use the cold curve. The main point is that the cold curve dictates the admissible set. But for the thermal portions, I don't know when someone might need those. I think I could move those ones into the private scope.

Yurlungur · 2025-09-11T13:36:40Z

singularity-eos/eos/eos_simple_macaw.hpp

+  static std::string EosPyType() { return EosType(); }
+
+ private:
+  Real _A, _B, _Cvinf, _v0, _T0, _Gc;


Not your fault for this---we're not consistent throughout the code. But in our newer EOS's we are trying to use a trailing underscore for private cope, rather than a leading underscore. Do you mind changing the naming convention here?

Yep, I can change this one. When should a leading underscore be used in general?

I don't think we use leading underscores anywhere. It's really just a matter of taste/convention.

Yurlungur · 2025-09-11T13:37:37Z

singularity-eos/eos/eos_simple_macaw.hpp

+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }


PORTABLE_ALWAYS_REQUIRE would work. You could also use PORTABLE_ALWAYS_THROW_OR_ABORT with a useful error message.

Yurlungur · 2025-09-11T13:38:16Z

singularity-eos/eos/eos_simple_macaw.hpp

+template <typename Indexer_t>
+PORTABLE_INLINE_FUNCTION void
+SimpleMACAW::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod,
+                  const unsigned long output, Indexer_t &&lambda) const {
+  if (output & thermalqs::density && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie);
+  }
+  if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::density || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
+  }
+  if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) {
+    sie = InternalEnergyFromDensityPressure(rho, press);
+  }
+  if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::temperature)
+    temp = TemperatureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::bulk_modulus)
+    bmod = BulkModulusFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::specific_heat)
+    cv = SpecificHeatFromDensityInternalEnergy(rho, sie);
+}


you're correct. Only reason to separate it out is for clarity. Choice left to the author.

Yurlungur · 2025-09-11T13:39:38Z

test/test_eos_simple_macaw.cpp

+using singularity::SimpleMACAW;
+using EOS = singularity::Variant<SimpleMACAW>;
+
+constexpr Real REAL_TOL = std::numeric_limits<Real>::epsilon() * 1e4; // 1e-12 for double


In this case is there a reason not to just set this to 1e-12? Is the 1e4 times machine epsilon significant or is this just what worked?

I was just following what another EOS test was using. I've changed it to 1e-12 now.

Yurlungur · 2025-09-11T13:41:30Z

test/test_eos_simple_macaw.cpp

+// Create the EOS
+auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc);
+auto eos = host_eos.GetOnDevice();


We probably don't want this here. You should put this inside the setup for a given scenario. (e.g., in line 48). That's so that at the end of the scenario you can call

eos.Finalize(); host_eos.Finalize();

and Catch will set up and tear down the EOS as appropriate.

I know that for an analytic EOS this really doesn't matter. But it does matter for more complicated EOS's so we should try to consistently use the known "safe" and "correct" style.

Yurlungur · 2025-09-11T13:44:01Z

test/test_eos_simple_macaw.cpp

+  }
+}
+
+SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][PT-RhoSie Equivalence]") {


See comment above about set up and tear down of the EOS object.

Yurlungur · 2025-09-16T18:08:51Z

singularity-eos/eos/eos_simple_macaw.hpp

+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }


I think this should always fail if you hit it, so I think PORTABLE_ALWAYS.

bgclayton-lanl requested a review from jhp-lanl September 8, 2025 21:51

bgclayton-lanl added the enhancement New feature or request label Sep 8, 2025

bgclayton-lanl added 7 commits September 9, 2025 10:03

First attempt at implementing an EOS haven't tested compiling yet.

ff01bf7

Got it compiling, need to test and add more tests.

7155400

Working on iterative method for computing density from pressure and t…

9417ea7

…emperature.

Got the nonlinear solve working for density from pressure temp. A few…

16fdc83

… more things to do.

Cleaned up the code some more and added another test.

d9b861c

Added documentation for the simple MACAW some more comments.

8f77c6d

Reordered the parameters to fix a clang sanitizer issue.

b9e6b59

bgclayton-lanl force-pushed the bgclayton_simple_macaw branch from 4df10b1 to b9e6b59 Compare September 9, 2025 16:03

jhp-lanl approved these changes Sep 11, 2025

View reviewed changes

Made suggested changes. A few more thing need to be finalized.

49ffdf2

Yurlungur reviewed Sep 16, 2025

View reviewed changes

Merge branch 'main' into bgclayton_simple_macaw

8cc311a

Simple MACAW EOS #560

Are you sure you want to change the base?

Simple MACAW EOS #560

Conversation

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

Uh oh!

PR Summary

PR Checklist

Uh oh!

Yurlungur commented Sep 9, 2025

Uh oh!

bgclayton-lanl commented Sep 10, 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!

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!

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!

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!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jhp-lanl 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!

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!

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!

Uh oh!

bgclayton-lanl commented Sep 8, 2025 •

edited

Loading