From 778a20b9c842c4cda63487778735ad40d2cf6431 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Tue, 17 Mar 2026 16:46:04 -0400
Subject: [PATCH 01/31] At least a simple implementation of
 PTDerivativesFromPerferred

---
 doc/sphinx/src/using-eos.rst        | 25 +++++++++++++++++++++++++
 singularity-eos/eos/eos_base.hpp    | 26 ++++++++++++++++++++++++++
 singularity-eos/eos/eos_variant.hpp | 14 ++++++++++++++
 3 files changed, 65 insertions(+)

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index faae9c7fbe7..bd33a235236 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -1334,6 +1334,31 @@ specific internal energy in :math:`erg/g`.
 
 The function
 
+.. code-block:: cpp
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
+
+computes the partial derivatives of density and specific internal
+energy with respect to pressure and temperature, with either pressure
+or temperature fixed. Density, energy, pressure, and temperature, must
+all be passed into this function, to maximize the ability for the
+underlying EOS to perform this calculation performantly. The intended
+use of this method is to compute the cell-averaged thermodynamic
+derivatives in a mixed cell in pressure-temperature equilibrium.
+
+.. warning::
+
+  Not all equations of state provide implementations of
+  ``PTDerivativesFromPreferred``, when an implementation is not
+  available, ``singularity-eos`` uses a finite differences
+  approximation.
+
+The function
+
 .. code-block:: cpp
 
    template<typename Indexer_t Real*>
diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 11cc5d261eb..8ce32a7782c 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -266,6 +266,32 @@ class EosBase {
     return sie + (P / rho) - T * S;
   }
 
+  // The base class computes these differences by finite
+  // differences. When available, this should absolutely be overloaded
+  // by an individual EOS, but this provides a default implementation
+  // at the very least.
+  //
+  // TODO(JMM): Should we have a vectorized version of this function?
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    Real r_pert, e_pert;
+
+    constexpr Real derivative_eps = 3.0e-6;
+    const Real dP = -P * derivative_eps;
+    const Real dT = T * derivative_eps;
+
+    DensityEnergyFromPressureTemperature(P + dP, T, lambda, r_pert, e_pert);
+    drdP_T = robust::ratio(r_pert - rho, dP);
+    dedP_T = robust::ratio(e_pert - sie, dP);
+
+    DensityEnergyFromPressureTemperature(P, T + dT, lambda, r_pert, e_pert);
+    drdT_P = robust::ratio(r_pert - rho, dT);
+    dedT_P = robust::ratio(e_pert - sie, dT);
+  }
+
   // Vector member functions
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index 691b6ad91ca..cb36ed8f5f5 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -280,6 +280,20 @@ class Variant {
         eos_);
   }
 
+  // TODO(JMM): Do we need a vectorized version of this call?
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    return PortsOfCall::visit(
+        [&](const auto &eos) {
+          return eos.PTDerivativesFromPreferred(rho, sie, P, T, lambda, dedP_T, drdP_T,
+                                                dedT_P, drdT_P);
+        },
+        eos_);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,

From 6c4e65f1a799976cac3581687430198e79d06952 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 10:44:29 -0400
Subject: [PATCH 02/31] factor out finite differences implementation

---
 singularity-eos/eos/eos_base.hpp | 39 ++++++++++++++++++++++----------
 1 file changed, 27 insertions(+), 12 deletions(-)

diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 8ce32a7782c..00aa7f9f1d7 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -184,6 +184,30 @@ struct MeanAtomicProperties {
   }
 };
 
+/* A version of PT Derivatives that uses only finite
+   differences. Factored out of the base class so it can be used to
+   cross-validate per-class implementations.
+*/
+template <typename EOS_t, typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void
+PTDerivativesByFiniteDifferences(const EOS_t &eos, const Real rho, const Real sie,
+                                 const Real P, const Real T, const Real derivative_eps,
+                                 Lambda_t &&lambda, Real &dedP_T, Real &drdP_T,
+                                 Real &dedT_P, Real &drdT_P) {
+  Real r_pert, e_pert;
+
+  const Real dP = -P * derivative_eps;
+  const Real dT = T * derivative_eps;
+
+  eos.DensityEnergyFromPressureTemperature(P + dP, T, lambda, r_pert, e_pert);
+  drdP_T = robust::ratio(r_pert - rho, dP);
+  dedP_T = robust::ratio(e_pert - sie, dP);
+
+  eos.DensityEnergyFromPressureTemperature(P, T + dT, lambda, r_pert, e_pert);
+  drdT_P = robust::ratio(r_pert - rho, dT);
+  dedT_P = robust::ratio(e_pert - sie, dT);
+}
+
 /*
 This is a CRTP that allows for static inheritance so that default behavior for
 various member functions can be defined.
@@ -277,19 +301,10 @@ class EosBase {
   PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
                              Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
                              Real &drdT_P) const {
-    Real r_pert, e_pert;
-
     constexpr Real derivative_eps = 3.0e-6;
-    const Real dP = -P * derivative_eps;
-    const Real dT = T * derivative_eps;
-
-    DensityEnergyFromPressureTemperature(P + dP, T, lambda, r_pert, e_pert);
-    drdP_T = robust::ratio(r_pert - rho, dP);
-    dedP_T = robust::ratio(e_pert - sie, dP);
-
-    DensityEnergyFromPressureTemperature(P, T + dT, lambda, r_pert, e_pert);
-    drdT_P = robust::ratio(r_pert - rho, dT);
-    dedT_P = robust::ratio(e_pert - sie, dT);
+    const CRTP &copy = *(static_cast<CRTP const *>(this));
+    PTDerivativesByFiniteDifferences(copy, rho, sie, P, T, derivative_eps, lambda, dedP_T,
+                                     drdP_T, dedT_P, drdT_P);
   }
 
   // Vector member functions

From 463ce742e19ca89b6eec9e969a1342d9efde31de Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 10:45:08 -0400
Subject: [PATCH 03/31] ideal gas derivatives

---
 singularity-eos/eos/eos_ideal.hpp | 19 +++++++++++++
 test/test_eos_ideal.cpp           | 47 +++++++++++++++++++++++++++++++
 2 files changed, 66 insertions(+)

diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
index 146a2b66a05..d057f2f564e 100644
--- a/singularity-eos/eos/eos_ideal.hpp
+++ b/singularity-eos/eos/eos_ideal.hpp
@@ -148,6 +148,25 @@ class IdealGas : public EosBase<IdealGas> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     return _gm1;
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    // P = (gm1) rho cv T = (gm1) rho sie
+
+    dedP_T = 0;
+    dedT_P = _Cv;
+
+    // => rho(P, T) = P / (gm1 cv T) = P / (gm1 sie)
+    drdP_T = robust::ratio(1.0, _gm1 * sie);
+
+    // rho(P, T) = (P / (gm1 cv)) T^{-1}
+    // => drdT_P = -1 (P / (gm1 cv)) T^{-2}
+    drdT_P = - robust::ratio(P, _gm1 * _Cv * T *T);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index dd2f5ab8d0d..c86abf0c2ad 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -198,6 +198,53 @@ SCENARIO("Ideal gas density energy from prssure temperature",
           nwrong);
       THEN("There are no errors") { REQUIRE(nwrong == 0); }
     }
+    WHEN("We check PT derivatives from preferred") {
+      constexpr int N = 100;
+      constexpr Real derivative_eps = 3e-6;
+      int nwrong = 0;
+      portableReduce(
+          "Check PT derivatives from preferred", 1, N,
+          PORTABLE_LAMBDA(const int i, int &nw) {
+            Real rho = i;
+            Real T = 100 * i;
+            Real P = device_eos.PressureFromDensityTemperature(rho, T);
+            Real sie = device_eos.InternalEnergyFromDensityTemperature(rho, T);
+            Real dedP_T, drdP_T, dedT_P, drdT_P;
+            Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+            device_eos.PTDerivativesFromPreferred(rho, sie, P, T,
+                                                  static_cast<Real *>(nullptr), dedP_T,
+                                                  drdP_T, dedT_P, drdT_P);
+            singularity::eos_base::PTDerivativesByFiniteDifferences(
+                device_eos, rho, sie, P, T, derivative_eps, static_cast<Real *>(nullptr),
+                dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
+            if (!isClose(dedP_T, dedP_T_fd)) {
+              printf("rho, sie, P, T, dedP_T expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, dedP_T, dedP_T_fd, dedP_T - dedP_T_fd);
+              nw += 1;
+            }
+            if (!isClose(drdP_T, drdP_T_fd)) {
+              printf("rho, sie, P, T, drdP_T expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, drdP_T, drdP_T_fd, drdP_T - drdP_T_fd);
+              nw += 1;
+            }
+            if (!isClose(dedT_P, dedT_P_fd)) {
+              printf("rho, sie, P, T, dedT_P expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, dedT_P, dedT_P_fd, dedT_P - dedT_P_fd);
+              nw += 1;
+            }
+            if (!isClose(drdT_P, drdT_P_fd)) {
+              printf("rho, sie, P, T, drdT_P expected, true, diff:\n\t"
+                     "%.14e %.14e %.14e %.14e %.14e %.14e %.14e\n",
+                     rho, sie, P, T, drdT_P, drdT_P_fd, drdT_P - drdT_P_fd);
+              nw += 1;
+            }
+          },
+          nwrong);
+      THEN("There are no errors") { REQUIRE(nwrong == 0); }
+    }
   }
 }
 

From 1f2bcb94f7f11f986df5619c53a9f369bad1444e Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 12:46:25 -0400
Subject: [PATCH 04/31] PT derivatives for spiner

---
 singularity-eos/eos/eos_spiner_rho_temp.hpp | 37 +++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 1ffb45ab8a4..3953daa1677 100644
--- a/singularity-eos/eos/eos_spiner_rho_temp.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp
@@ -196,6 +196,11 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -1005,6 +1010,38 @@ PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoT::DensityEnergyFromPressureTem
   sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
 }
 
+template <typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void
+PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                           Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                           Real &drdT_P) const {
+  const Real lT = lT_(temp);
+  const Real lRho = lRho_(rho);
+  const TableStatus whereAmI = getLocDependsRhoT_(lRho, lT);
+
+  if (whereAmI == TableStatus::OffBottom) {
+    dedP_T = robust::ratio(1., dPdECold_.interpToReal(lRho));
+    drdP_T = robust::ratio(1., dPdRhoCold_.interpToReal(lRho));
+    dedT_P = dedTCold_.interpToReal(lRho);
+    drdT_P = 0; // TODO(JMM) is this right?
+  } else if (whereAmI == TableStatus::OffTop) { // idela gas
+    const Real gm1 = gm1Max_.interpToReal(lRho);
+    const Real Cv = dEdTMax_.interpToReal(lRho);
+    dedP_T = 0;
+    dedT_P = Cv;
+    drdP_T = robust::ratio(1.0, gm1 * sie);
+    drdT_P = - robust::ratio(P, gm1 * Cv * T *T);
+  } else {
+    dedP_T = robust::ratio(1., dPdE_.interpToReal(lRho, lT));
+    drdP_T = robust::ratio(1., dPdRho_.interpToReal(lRho, lT));
+    dedT_P = dedT_.interpToReal(lRho, lT);
+    // dP = (dP/dr)_T dr + (dP/dT)_r dT
+    // dP = 0
+    // => 0 = (dP/dr)_T dr + (dP/dT)_r dT
+    drdT_P = -robust::ratio(dPdT_.interpToReal(lRho, lT), dPdRho_.interpToReal(lRho, lT));
+  }
+}
+
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 SpinerEOSDependsRhoT::FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv,

From 8ad0b203cdf4e4673015435580d4eb88a4647a46 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 13:15:19 -0400
Subject: [PATCH 05/31] thread through modifiers. Make a note about tables.

---
 singularity-eos/eos/eos_spiner_rho_sie.hpp    |  8 ++++++
 singularity-eos/eos/eos_spiner_rho_temp.hpp   | 25 ++++++++++++-------
 .../eos/modifiers/eos_unitsystem.hpp          | 15 +++++++++++
 .../eos/modifiers/floored_energy.hpp          | 11 ++++++++
 singularity-eos/eos/modifiers/ramps_eos.hpp   |  3 +++
 .../eos/modifiers/relativistic_eos.hpp        |  8 ++++++
 singularity-eos/eos/modifiers/scaled_eos.hpp  | 14 +++++++++++
 singularity-eos/eos/modifiers/shifted_eos.hpp | 13 ++++++++--
 singularity-eos/eos/modifiers/zsplit_eos.hpp  | 15 +++++++++++
 9 files changed, 101 insertions(+), 11 deletions(-)

diff --git a/singularity-eos/eos/eos_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp
index f0c2986fa4d..617ca8ea732 100644
--- a/singularity-eos/eos/eos_spiner_rho_sie.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp
@@ -215,6 +215,14 @@ class SpinerEOSDependsRhoSieTransformable
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  /*
+  // TODO(JMM): For now using FD. Fix this.
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
+  */
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 3953daa1677..0f8b49e75a0 100644
--- a/singularity-eos/eos/eos_spiner_rho_temp.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp
@@ -196,11 +196,14 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const;
+  /*
+  // TODO(JMM): For now using FD. Fix this.
   template <typename Lambda_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
                              Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
                              Real &drdT_P) const;
+  */
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -1010,37 +1013,41 @@ PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoT::DensityEnergyFromPressureTem
   sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
 }
 
+/*
+// TODO(JMM): I would like to use this, but shockingly we don't have
+// dPdT stored. So need to add that to sesame2spiner first. Possibly
+// should be a different MR.
 template <typename Lambda_t>
-PORTABLE_INLINE_FUNCTION void
-PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
-                           Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
-                           Real &drdT_P) const {
-  const Real lT = lT_(temp);
+PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoT::PTDerivativesFromPreferred(
+    const Real rho, const Real sie, const Real P, const Real T, Lambda_t &&lambda,
+    Real &dedP_T, Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+  const Real lT = lT_(T);
   const Real lRho = lRho_(rho);
   const TableStatus whereAmI = getLocDependsRhoT_(lRho, lT);
 
   if (whereAmI == TableStatus::OffBottom) {
     dedP_T = robust::ratio(1., dPdECold_.interpToReal(lRho));
     drdP_T = robust::ratio(1., dPdRhoCold_.interpToReal(lRho));
-    dedT_P = dedTCold_.interpToReal(lRho);
-    drdT_P = 0; // TODO(JMM) is this right?
+    dedT_P = dEdTCold_.interpToReal(lRho);
+    drdT_P = 0;                                 // TODO(JMM) is this right?
   } else if (whereAmI == TableStatus::OffTop) { // idela gas
     const Real gm1 = gm1Max_.interpToReal(lRho);
     const Real Cv = dEdTMax_.interpToReal(lRho);
     dedP_T = 0;
     dedT_P = Cv;
     drdP_T = robust::ratio(1.0, gm1 * sie);
-    drdT_P = - robust::ratio(P, gm1 * Cv * T *T);
+    drdT_P = -robust::ratio(P, gm1 * Cv * T * T);
   } else {
     dedP_T = robust::ratio(1., dPdE_.interpToReal(lRho, lT));
     drdP_T = robust::ratio(1., dPdRho_.interpToReal(lRho, lT));
-    dedT_P = dedT_.interpToReal(lRho, lT);
+    dedT_P = dEdT_.interpToReal(lRho, lT);
     // dP = (dP/dr)_T dr + (dP/dT)_r dT
     // dP = 0
     // => 0 = (dP/dr)_T dr + (dP/dT)_r dT
     drdT_P = -robust::ratio(dPdT_.interpToReal(lRho, lT), dPdRho_.interpToReal(lRho, lT));
   }
 }
+*/
 
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 8acbf19025b..9baf657c242 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -214,6 +214,21 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
     sie *= inv_sie_unit_;
   }
 
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho * rho_unit_, sie * sie_unit_, P * press_unit_,
+                                  temp * temp_unit_, lambda, dedP_T, drdP_T, dedT_P,
+                                  drdT_P);
+    // scale outputs
+    dedP_T *= robust::ratio(press_unit_, temp_unit_);
+    drdP_T *= robust::ratio(press_unit_, rho_unit_);
+    dedT_P *= robust::ratio(temp_unit_, sie_unit_);
+    drdT_P *= robust::ratio(temp_unit_, rho_unit_);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
index 314cd33a432..b2700d4b7ca 100644
--- a/singularity-eos/eos/modifiers/floored_energy.hpp
+++ b/singularity-eos/eos/modifiers/floored_energy.hpp
@@ -126,6 +126,17 @@ class FlooredEnergy : public EosBase<FlooredEnergy<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real temp,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const {
+    const Real min_sie = t_.MinInternalEnergyFromDensity(rho);
+    t_.PTDerivativesFromPreferred(rho, std::max(sie, min_sie), P, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp
index 1891b1d81ca..4bec46dd247 100644
--- a/singularity-eos/eos/modifiers/ramps_eos.hpp
+++ b/singularity-eos/eos/modifiers/ramps_eos.hpp
@@ -213,6 +213,9 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     // TODO(JMM): This doesn't seem right. dpdrho relevant.
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
+
+  /* JMM: Use the base class PTDerivativesFromPreferred for now */
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index 8863a7295a9..a63b15ba0b2 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -145,6 +145,14 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
   }
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho, sie, P, temp, lambda, dedP_T, drdP_T, dedT_P,
+                                  drdT_P);
+  }
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index 63851ded5fb..11cce2cb99d 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -141,6 +141,20 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
                                     Indexer_t &&lambda = nullptr) const {
     t_.InternalEnergyFromDensityPressure(scale_ * rho, P, sie, lambda);
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(scale_ * rho, sie, P, temp, lambda, dedP_T, drdP_T,
+                                  dedT_P, drdT_P);
+    dedP_T *= scale_;
+    drdP_T *= inv_scale_;
+    dedT_P *= scale_;
+    drdT_P *= inv_scale_;
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index cd0bad539b6..40a7614203e 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -134,6 +134,16 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
     sie += shift_;
   }
+
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    t_.PTDerivativesFromPreferred(rho, sie - shift_, P, temp, lambda, dedP_T, drdP_T,
+                                  dedT_P, drdT_P);
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
@@ -169,8 +179,7 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
         singularity::mfuncname::member_func_name(typeid(ShiftedEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
     const auto shift_val = shift_;
-    portableFor(
-        cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_val; });
+    portableFor(cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_val; });
   }
 
   template <typename LambdaIndexer>
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index ee9f260e400..9abc2cf1a2f 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -180,6 +180,21 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     sie *= scale;
   }
 
+  template <typename Lambda_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
+    const Real scale = GetScale_(lambda);
+    const Real iscale = GetInvScale_(lambda);
+    t_.PTDerivativesFromPreferred(rho, sie * iscale, P * iscale, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
+    // dedP scale cancels
+    drdP_T *= iscale;
+    dedT_P *= scale;
+    // drdT_P does not need scaling
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,

From ed3dd8036ba0696c4492250de5c7fff1596b0b67 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 13:28:44 -0400
Subject: [PATCH 06/31] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 63e5dbf91fb..83201622159 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR618]](https://github.com/lanl/singularity-eos/pull/618) Add PTDerivativesFromPreferred for computing derivatives of a mixture in a cell
 - [[PR600]](https://github.com/lanl/singularity-eos/pull/600) Added mass fractions for multiphase tables to sesame2spiner and added MassFractionsFromDensityTemperture/InternalEnergy to the SpinerDependsRhoSIE and SpinerDependsRhoT equations of state.
 - [[PR589]](https://github.com/lanl/singularity-eos/pull/589) InternalEnergyFromDensityPressure
 

From 3101b31ad68b8e74dd63bcffb524dc43c6c4954b Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 14:02:40 -0400
Subject: [PATCH 07/31] eospac implementation

---
 singularity-eos/eos/eos_eospac.hpp | 40 ++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 82f2f09e880..f39a0abc7ef 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -231,6 +231,11 @@ class EOSPAC : public EosBase<EOSPAC> {
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+  template <typename Lambda_t>
+  PORTABLE_INLINE_FUNCTION void
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real T,
+                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                             Real &drdT_P) const;
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -1710,6 +1715,41 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::DensityEnergyFromPressureTemperature(
 #endif // ON DEVICE
 }
 
+template <typename Lambda_t>
+PORTABLE_INLINE_FUNCTION void
+EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real press, const Real temp,
+                                   Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
+                                   Real &drdT_P) const {
+#if SINGULARITY_ON_DEVICE
+  PORTABLE_ALWAYS_ABORT("EOSPAC calls not supported on device\n");
+#else
+  PORTABLE_REQUIRE(split_ == TableSplit::Total,
+                   "Density of pressure and temperature only supported for total "
+                   "tables at this time");
+  using namespace EospacWrapper;
+  EOS_REAL P[1] = {pressureToSesame(press)};
+  EOS_REAL T[1] = {temperatureToSesame(temp)};
+  EOS_REAL R[1] = {rho};
+  EOS_REAL E[1] = {sieToSesame(sie)};
+  EOS_REAL z[1], dx[1], dy[1];
+  EOS_INTEGER nxypairs = 1;
+  EOS_INTEGER table;
+
+  table = PofRE_table_;
+  eosSafeInterpolate(&table, nxypairs, z, R, E, dx, dy, "PofRE", Verbosity::Quiet);
+  dedP_T = robust::ratio(1., dy[0]);
+
+  table = EofRT_table_;
+  eosSafeInterpolate(&table, nxypairs, z, R, T, dx, dy, "EofRT", Verbosity::Quiet);
+  dedT_P = dy[0];
+
+  table = RofPT_table_;
+  eosSafeInterpolate(&table, nxypairs, z, P, T, dx, dy, "RofPT", Verbosity::Quiet);
+  drdP_T = dx[0];
+  drdT_P = dx[1];
+#endif // ON DEVICE
+}
+
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 EOSPAC::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,

From fc5367df6f9c8dc671e14f8eb00a53b9ffebd463 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 16:07:32 -0400
Subject: [PATCH 08/31] add tests for modifiers

---
 test/test_eos_modifiers.cpp | 109 ++++++++++++++++++++++++++----------
 1 file changed, 80 insertions(+), 29 deletions(-)

diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp
index 859749f28ea..a585e04f815 100644
--- a/test/test_eos_modifiers.cpp
+++ b/test/test_eos_modifiers.cpp
@@ -122,6 +122,64 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
           eos_ramped.PrintParams();
         }
       }
+
+      WHEN("We compute PT derivatives from preferred for the modified EOS") {
+        // Choose a representative state
+        constexpr Real test_rho = 1.5;
+        constexpr Real test_T   = 2.0;
+
+        // Compute primitive quantities
+        const Real P   = eos.PressureFromDensityTemperature(test_rho, test_T);
+        const Real sie_val = eos.InternalEnergyFromDensityTemperature(test_rho, test_T);
+
+        // Derivatives from the EOS preferred interface
+        Real dedP_T, drdP_T, dedT_P, drdT_P;
+        eos.PTDerivativesFromPreferred(test_rho, sie_val, P, test_T,
+                                      static_cast<Real *>(nullptr),
+                                      dedP_T, drdP_T, dedT_P, drdT_P);
+
+        // Finite‑difference reference
+        constexpr Real derivative_eps = 3e-6;
+        Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+        singularity::eos_base::PTDerivativesByFiniteDifferences(
+            eos, test_rho, sie_val, P, test_T, derivative_eps,
+            static_cast<Real *>(nullptr),
+            dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
+
+        // Compare each component
+        THEN("dedP_T matches finite‑difference reference") {
+          bool close = isClose(dedP_T, dedP_T_fd);
+          if (!close) {
+            printf("dedP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
+                   dedP_T, dedP_T_fd, dedP_T - dedP_T_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("drdP_T matches finite‑difference reference") {
+          bool close = isClose(drdP_T, drdP_T_fd);
+          if (!close) {
+            printf("drdP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
+                   drdP_T, drdP_T_fd, drdP_T - drdP_T_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("dedT_P matches finite‑difference reference") {
+          bool close = isClose(dedT_P, dedT_P_fd);
+          if (!close) {
+            printf("dedT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
+                   dedT_P, dedT_P_fd, dedT_P - dedT_P_fd);
+          }
+          REQUIRE(close);
+        }
+        THEN("drdT_P matches finite‑difference reference") {
+          bool close = isClose(drdT_P, drdT_P_fd);
+          if (!close) {
+            printf("drdT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
+                   drdT_P, drdT_P_fd, drdT_P - drdT_P_fd);
+          }
+          REQUIRE(close);
+        }
+      }
     }
 #ifdef SINGULARITY_BUILD_CLOSURE
     WHEN("We construct a non-modifying modifier") {
@@ -170,56 +228,49 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
       // bmod (rho_t2)
       const Real bmodrt2 = rho_t2 * b / r0;
       THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") {
-        INFO("P_eos(alpha_0*rmid, T0): "
-             << P_armid
-             << " P_ramp(rmid, T0): " << igra.PressureFromDensityTemperature(rmid, T0));
+        printf("P_eos(alpha_0*rmid, T0): %.14e  P_ramp(rmid, T0): %.14e\n",
+               P_armid, igra.PressureFromDensityTemperature(rmid, T0));
         REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in P(rho) for rho <r0, [r0,rmid], [rmid,r1] "
            "and >r1") {
         // also check pressures on ramp
-        INFO("reference P((r0+rmid)/2, T0): "
-             << Prhot1 << " test P((r0+rmid)/2, T0): "
-             << igra.PressureFromDensityTemperature(rho_t1, T0));
+        printf("reference P((r0+rmid)/2, T0): %.14e  test P((r0+rmid)/2, T0): %.14e\n",
+               Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0));
         REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12));
-        INFO("reference P((rmid+r1)/2, T0): "
-             << Prhot2 << " test P((rmid+r1)/2, T0): "
-             << igra.PressureFromDensityTemperature(rho_t2, T0));
+        printf("reference P((rmid+r1)/2, T0): %.14e  test P((rmid+r1)/2, T0): %.14e\n",
+               Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0));
         REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check pressure below and beyond ramp matches unmodified ideal gas
-        INFO("reference P(0.8*r0, T0): "
-             << 0.8 * r0 * gm1 * Cv * T0 << " test P(0.8*r0, T0): "
-             << igra.PressureFromDensityTemperature(0.8 * r0, T0));
+        printf("reference P(0.8*r0, T0): %.14e  test P(0.8*r0, T0): %.14e\n",
+               0.8 * r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(0.8 * r0, T0));
         REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0,
                         igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12));
-        INFO("reference P(1.2*r1, T0): "
-             << 1.2 * r1 * gm1 * Cv * T0 << " test P(1.2*r1, T0): "
-             << igra.PressureFromDensityTemperature(1.2 * r1, T0));
-        REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0,
-                        igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12));
+        printf("reference P(1.2*r1, T0): %.14e  test P(1.2*r1, T0): %.14e\n",
+               1.2 * r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(1.2 * r0, T0));
+        REQUIRE(isClose(1.2 * r0 * gm1 * Cv * T0,
+                        igra.PressureFromDensityTemperature(1.2 * r0, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in bmod(rho) for rho <r0, [r0,rmid], "
            "[rmid,r1] and >r1") {
         // check bulk moduli on both pieces of ramp
-        INFO("reference bmod((r0+rmid)/2, T0): "
-             << bmodrt1 << " test bmod((r0+rmid)/2, T0): "
-             << igra.BulkModulusFromDensityTemperature(rho_t1, T0));
+        printf("reference bmod((r0+rmid)/2, T0): %.14e  test bmod((r0+rmid)/2, T0): %.14e\n",
+               bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0));
         REQUIRE(
             isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12));
-        INFO("reference bmod((rmid+r1)/2, T0): "
-             << bmodrt2 << " test bmod((rmid+r1)/2, T0): "
-             << igra.BulkModulusFromDensityTemperature(rho_t2, T0));
+        printf("reference bmod((rmid+r1)/2, T0): %.14e  test bmod((rmid+r1)/2, T0): %.14e\n",
+               bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0));
         REQUIRE(
             isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check bulk modulus below and beyond ramp matches unmodified ideal gas
-        INFO("reference bmod(0.8*r0, T0): "
-             << 0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(0.8*r0, T0): "
-             << igra.BulkModulusFromDensityTemperature(0.8 * r0, T0));
+        printf("reference bmod(0.8*r0, T0): %.14e  test bmod(0.8*r0, T0): %.14e\n",
+               0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0,
+               igra.BulkModulusFromDensityTemperature(0.8 * r0, T0));
         REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0,
                         igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12));
-        INFO("reference bmod(1.2*r1, T0): "
-             << 1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(1.2*r1, T0): "
-             << igra.BulkModulusFromDensityTemperature(1.2 * r1, T0));
+        printf("reference bmod(1.2*r1, T0): %.14e  test bmod(1.2*r1, T0): %.14e\n",
+               1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0,
+               igra.BulkModulusFromDensityTemperature(1.2 * r1, T0));
         REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0,
                         igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12));
       }

From ea7f14fda508f086e6555c619db8ccc64899dcd9 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 16:08:46 -0400
Subject: [PATCH 09/31] formatting

---
 singularity-eos/eos/eos_eospac.hpp            |  6 +--
 singularity-eos/eos/eos_ideal.hpp             |  2 +-
 .../eos/modifiers/floored_energy.hpp          |  6 +--
 test/test_eos_modifiers.cpp                   | 47 ++++++++++---------
 4 files changed, 32 insertions(+), 29 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index f39a0abc7ef..f72268a1d14 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1717,9 +1717,9 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::DensityEnergyFromPressureTemperature(
 
 template <typename Lambda_t>
 PORTABLE_INLINE_FUNCTION void
-EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real press, const Real temp,
-                                   Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
-                                   Real &drdT_P) const {
+EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real press,
+                                   const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                                   Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
 #if SINGULARITY_ON_DEVICE
   PORTABLE_ALWAYS_ABORT("EOSPAC calls not supported on device\n");
 #else
diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
index d057f2f564e..0885ae4658e 100644
--- a/singularity-eos/eos/eos_ideal.hpp
+++ b/singularity-eos/eos/eos_ideal.hpp
@@ -164,7 +164,7 @@ class IdealGas : public EosBase<IdealGas> {
 
     // rho(P, T) = (P / (gm1 cv)) T^{-1}
     // => drdT_P = -1 (P / (gm1 cv)) T^{-2}
-    drdT_P = - robust::ratio(P, _gm1 * _Cv * T *T);
+    drdT_P = -robust::ratio(P, _gm1 * _Cv * T * T);
   }
 
   template <typename Indexer_t = Real *>
diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
index b2700d4b7ca..e1d77764d07 100644
--- a/singularity-eos/eos/modifiers/floored_energy.hpp
+++ b/singularity-eos/eos/modifiers/floored_energy.hpp
@@ -129,9 +129,9 @@ class FlooredEnergy : public EosBase<FlooredEnergy<T>> {
 
   template <typename Lambda_t = Real *>
   PORTABLE_INLINE_FUNCTION void
-  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P, const Real temp,
-                             Lambda_t &&lambda, Real &dedP_T, Real &drdP_T, Real &dedT_P,
-                             Real &drdT_P) const {
+  PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
+                             const Real temp, Lambda_t &&lambda, Real &dedP_T,
+                             Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
     const Real min_sie = t_.MinInternalEnergyFromDensity(rho);
     t_.PTDerivativesFromPreferred(rho, std::max(sie, min_sie), P, temp, lambda, dedP_T,
                                   drdP_T, dedT_P, drdT_P);
diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp
index a585e04f815..4211de826db 100644
--- a/test/test_eos_modifiers.cpp
+++ b/test/test_eos_modifiers.cpp
@@ -126,56 +126,55 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
       WHEN("We compute PT derivatives from preferred for the modified EOS") {
         // Choose a representative state
         constexpr Real test_rho = 1.5;
-        constexpr Real test_T   = 2.0;
+        constexpr Real test_T = 2.0;
 
         // Compute primitive quantities
-        const Real P   = eos.PressureFromDensityTemperature(test_rho, test_T);
+        const Real P = eos.PressureFromDensityTemperature(test_rho, test_T);
         const Real sie_val = eos.InternalEnergyFromDensityTemperature(test_rho, test_T);
 
         // Derivatives from the EOS preferred interface
         Real dedP_T, drdP_T, dedT_P, drdT_P;
         eos.PTDerivativesFromPreferred(test_rho, sie_val, P, test_T,
-                                      static_cast<Real *>(nullptr),
-                                      dedP_T, drdP_T, dedT_P, drdT_P);
+                                       static_cast<Real *>(nullptr), dedP_T, drdP_T,
+                                       dedT_P, drdT_P);
 
         // Finite‑difference reference
         constexpr Real derivative_eps = 3e-6;
         Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
         singularity::eos_base::PTDerivativesByFiniteDifferences(
             eos, test_rho, sie_val, P, test_T, derivative_eps,
-            static_cast<Real *>(nullptr),
-            dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
+            static_cast<Real *>(nullptr), dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd);
 
         // Compare each component
         THEN("dedP_T matches finite‑difference reference") {
           bool close = isClose(dedP_T, dedP_T_fd);
           if (!close) {
-            printf("dedP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
-                   dedP_T, dedP_T_fd, dedP_T - dedP_T_fd);
+            printf("dedP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", dedP_T,
+                   dedP_T_fd, dedP_T - dedP_T_fd);
           }
           REQUIRE(close);
         }
         THEN("drdP_T matches finite‑difference reference") {
           bool close = isClose(drdP_T, drdP_T_fd);
           if (!close) {
-            printf("drdP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
-                   drdP_T, drdP_T_fd, drdP_T - drdP_T_fd);
+            printf("drdP_T mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", drdP_T,
+                   drdP_T_fd, drdP_T - drdP_T_fd);
           }
           REQUIRE(close);
         }
         THEN("dedT_P matches finite‑difference reference") {
           bool close = isClose(dedT_P, dedT_P_fd);
           if (!close) {
-            printf("dedT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
-                   dedT_P, dedT_P_fd, dedT_P - dedT_P_fd);
+            printf("dedT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", dedT_P,
+                   dedT_P_fd, dedT_P - dedT_P_fd);
           }
           REQUIRE(close);
         }
         THEN("drdT_P matches finite‑difference reference") {
           bool close = isClose(drdT_P, drdT_P_fd);
           if (!close) {
-            printf("drdT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n",
-                   drdT_P, drdT_P_fd, drdT_P - drdT_P_fd);
+            printf("drdT_P mismatch: eos = %.14e, fd = %.14e, diff = %.14e\n", drdT_P,
+                   drdT_P_fd, drdT_P - drdT_P_fd);
           }
           REQUIRE(close);
         }
@@ -228,8 +227,8 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
       // bmod (rho_t2)
       const Real bmodrt2 = rho_t2 * b / r0;
       THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") {
-        printf("P_eos(alpha_0*rmid, T0): %.14e  P_ramp(rmid, T0): %.14e\n",
-               P_armid, igra.PressureFromDensityTemperature(rmid, T0));
+        printf("P_eos(alpha_0*rmid, T0): %.14e  P_ramp(rmid, T0): %.14e\n", P_armid,
+               igra.PressureFromDensityTemperature(rmid, T0));
         REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in P(rho) for rho <r0, [r0,rmid], [rmid,r1] "
@@ -243,23 +242,27 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") {
         REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check pressure below and beyond ramp matches unmodified ideal gas
         printf("reference P(0.8*r0, T0): %.14e  test P(0.8*r0, T0): %.14e\n",
-               0.8 * r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(0.8 * r0, T0));
+               0.8 * r0 * gm1 * Cv * T0,
+               igra.PressureFromDensityTemperature(0.8 * r0, T0));
         REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0,
                         igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12));
         printf("reference P(1.2*r1, T0): %.14e  test P(1.2*r1, T0): %.14e\n",
-               1.2 * r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(1.2 * r0, T0));
+               1.2 * r0 * gm1 * Cv * T0,
+               igra.PressureFromDensityTemperature(1.2 * r0, T0));
         REQUIRE(isClose(1.2 * r0 * gm1 * Cv * T0,
                         igra.PressureFromDensityTemperature(1.2 * r0, T0), 1.e-12));
       }
       THEN("We obtain correct ramp behavior in bmod(rho) for rho <r0, [r0,rmid], "
            "[rmid,r1] and >r1") {
         // check bulk moduli on both pieces of ramp
-        printf("reference bmod((r0+rmid)/2, T0): %.14e  test bmod((r0+rmid)/2, T0): %.14e\n",
-               bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0));
+        printf(
+            "reference bmod((r0+rmid)/2, T0): %.14e  test bmod((r0+rmid)/2, T0): %.14e\n",
+            bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0));
         REQUIRE(
             isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12));
-        printf("reference bmod((rmid+r1)/2, T0): %.14e  test bmod((rmid+r1)/2, T0): %.14e\n",
-               bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0));
+        printf(
+            "reference bmod((rmid+r1)/2, T0): %.14e  test bmod((rmid+r1)/2, T0): %.14e\n",
+            bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0));
         REQUIRE(
             isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12));
         // check bulk modulus below and beyond ramp matches unmodified ideal gas

From 08ff08c439e9c70e0524933fc9ed6824f8f58172 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 16:11:23 -0400
Subject: [PATCH 10/31] formatting

---
 singularity-eos/eos/modifiers/shifted_eos.hpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index 40a7614203e..8b5db295131 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -179,7 +179,8 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
         singularity::mfuncname::member_func_name(typeid(ShiftedEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
     const auto shift_val = shift_;
-    portableFor(cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_val; });
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_val; });
   }
 
   template <typename LambdaIndexer>

From b5ff9e31e72b243499cb2c6796e9bffa0eb3efc0 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 18 Mar 2026 16:16:07 -0400
Subject: [PATCH 11/31] CC

---
 singularity-eos/eos/eos_base.hpp                   | 2 +-
 singularity-eos/eos/eos_eospac.hpp                 | 2 +-
 singularity-eos/eos/eos_variant.hpp                | 2 +-
 singularity-eos/eos/modifiers/eos_unitsystem.hpp   | 2 +-
 singularity-eos/eos/modifiers/floored_energy.hpp   | 2 +-
 singularity-eos/eos/modifiers/ramps_eos.hpp        | 2 +-
 singularity-eos/eos/modifiers/relativistic_eos.hpp | 2 +-
 singularity-eos/eos/modifiers/scaled_eos.hpp       | 2 +-
 singularity-eos/eos/modifiers/shifted_eos.hpp      | 2 +-
 singularity-eos/eos/modifiers/zsplit_eos.hpp       | 2 +-
 test/test_eos_ideal.cpp                            | 2 +-
 test/test_eos_modifiers.cpp                        | 2 +-
 12 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 00aa7f9f1d7..b0b3f9fb434 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.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
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index f72268a1d14..6a596606858 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.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
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index cb36ed8f5f5..463835ecc37 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.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
diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 9baf657c242..a0dfd53f113 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.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
diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
index e1d77764d07..53a2ec7ac99 100644
--- a/singularity-eos/eos/modifiers/floored_energy.hpp
+++ b/singularity-eos/eos/modifiers/floored_energy.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
diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp
index 4bec46dd247..26d94987181 100644
--- a/singularity-eos/eos/modifiers/ramps_eos.hpp
+++ b/singularity-eos/eos/modifiers/ramps_eos.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
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index a63b15ba0b2..6ca7606ae65 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.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
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index 11cce2cb99d..1d21825b671 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.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
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index 8b5db295131..4c856280c35 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.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
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index 9abc2cf1a2f..4e5c3bde589 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2024-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2024-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
diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp
index c86abf0c2ad..2faf1f6b33a 100644
--- a/test/test_eos_ideal.cpp
+++ b/test/test_eos_ideal.cpp
@@ -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
diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp
index 4211de826db..c9252b28c4a 100644
--- a/test/test_eos_modifiers.cpp
+++ b/test/test_eos_modifiers.cpp
@@ -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

From 4d3a781230333b8a1c59fb0c29ffadddff6d62a4 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 10:46:46 -0400
Subject: [PATCH 12/31] account for zero or negative pressures/temperatures

---
 singularity-eos/eos/eos_base.hpp | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index b0b3f9fb434..37d1738996d 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -196,8 +196,11 @@ PTDerivativesByFiniteDifferences(const EOS_t &eos, const Real rho, const Real si
                                  Real &dedT_P, Real &drdT_P) {
   Real r_pert, e_pert;
 
-  const Real dP = -P * derivative_eps;
-  const Real dT = T * derivative_eps;
+  const Real absP = std::max(std::abs(P), robust::EPS() / derivative_eps);
+  const Real dP = -robust::sgn(P) * absP * derivative_eps; // might be negative
+
+  const Real absT = std::max(std::abs(T), robust::EPS() / derivative_eps);
+  const Real dT = absT * derivative_eps; // can't be negative
 
   eos.DensityEnergyFromPressureTemperature(P + dP, T, lambda, r_pert, e_pert);
   drdP_T = robust::ratio(r_pert - rho, dP);

From 1ec6dfe61ad05dbd4585671c32a796da6a0df848 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 11:01:40 -0400
Subject: [PATCH 13/31] fix formula for eospac

---
 singularity-eos/eos/eos_eospac.hpp | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 6a596606858..3f4d16f3a16 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1735,18 +1735,17 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   EOS_INTEGER nxypairs = 1;
   EOS_INTEGER table;
 
-  table = PofRE_table_;
-  eosSafeInterpolate(&table, nxypairs, z, R, E, dx, dy, "PofRE", Verbosity::Quiet);
-  dedP_T = robust::ratio(1., dy[0]);
-
   table = EofRT_table_;
   eosSafeInterpolate(&table, nxypairs, z, R, T, dx, dy, "EofRT", Verbosity::Quiet);
-  dedT_P = dy[0];
+  const Real dedr_T = dx[0];
+  const Real dedT_r = dy[0];
 
   table = RofPT_table_;
   eosSafeInterpolate(&table, nxypairs, z, P, T, dx, dy, "RofPT", Verbosity::Quiet);
   drdP_T = dx[0];
   drdT_P = dx[1];
+  dedP_T = dedr_T * drdP_T;
+  dedT_P = dedT_r + dedr_T * drdT_P;
 #endif // ON DEVICE
 }
 

From bd5949083cc285d4cd57f4631e01407090dc91bc Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 12:07:07 -0400
Subject: [PATCH 14/31] oops forgot the inverse scale for sie

---
 singularity-eos/eos/modifiers/scaled_eos.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index 1d21825b671..aa60f0b3241 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -147,8 +147,8 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
   PTDerivativesFromPreferred(const Real rho, const Real sie, const Real P,
                              const Real temp, Lambda_t &&lambda, Real &dedP_T,
                              Real &drdP_T, Real &dedT_P, Real &drdT_P) const {
-    t_.PTDerivativesFromPreferred(scale_ * rho, sie, P, temp, lambda, dedP_T, drdP_T,
-                                  dedT_P, drdT_P);
+    t_.PTDerivativesFromPreferred(scale_ * rho, inv_scale_ * sie, P, temp, lambda, dedP_T,
+                                  drdP_T, dedT_P, drdT_P);
     dedP_T *= scale_;
     drdP_T *= inv_scale_;
     dedT_P *= scale_;

From a56cf66afc0d12e9f9296c2ef25cde8e2ee6f035 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 12:16:32 -0400
Subject: [PATCH 15/31] Test EOSPAC derivatives from preferred

---
 test/test_eos_tabulated.cpp | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index 0ad2c0d32ab..12b52bf22d6 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -97,6 +97,24 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
     auto steelEOS_host = steelEOS_host_polymorphic.get<SpinerEOSDependsRhoT>();
 
     EOS eospac = EOSPAC(steelID);
+    THEN("The EOSPAC model can provide derivatives from preferred") {
+      Real rho = 5;
+      Real T = 500;
+      Real P = eospac.PressureFromDensityTemperature(rho, T);
+      Real sie = eospac.InternalEnergyFromDensityTemperature(rho, T);
+
+      Real dedP_T, drdP_T, dedT_P, drdT_P;
+      Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
+      eospac.PTDerivativesFromPreferred(rho, sie, P, T, static_cast<Real *>(nullptr),
+                                        dedP_T, drdP_T, dedT_P, drdT_P);
+      singularity::eos_base::PTDerivativesByFiniteDifferences(
+          eospac, rho, sie, P, T, derivative_eps, static_cast<Real *>(nullptr), dedP_T_fd,
+          drdP_T_fd, dedT_P_fd, drdT_P_fd);
+      REQUIRE(isClose(dedP_T, dedP_T_fd));
+      REQUIRE(isClose(drdP_T, drdP_T_fd));
+      REQUIRE(isClose(dedT_P, dedT_P_fd));
+      REQUIRE(isClose(drdT_P, drdT_P_fd));
+    }
 
     THEN("The correct metadata is read in") {
       REQUIRE(steelEOS_host.matid() == steelID);

From 3dabe6f4dd1f6efe67168d886005131f44c0d482 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 12:40:08 -0400
Subject: [PATCH 16/31] need derivative eps

---
 test/test_eos_tabulated.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index 12b52bf22d6..d4bfb68fbd3 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -103,6 +103,7 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
       Real P = eospac.PressureFromDensityTemperature(rho, T);
       Real sie = eospac.InternalEnergyFromDensityTemperature(rho, T);
 
+      constexpr Real derivative_eps = 3e-6;
       Real dedP_T, drdP_T, dedT_P, drdT_P;
       Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
       eospac.PTDerivativesFromPreferred(rho, sie, P, T, static_cast<Real *>(nullptr),

From 84cd5240e411b3290949a7f8c9370c33e60257f1 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 14:39:23 -0400
Subject: [PATCH 17/31] OK Lets try this

---
 test/test_eos_tabulated.cpp | 28 ++++++++++++++++++++++++----
 1 file changed, 24 insertions(+), 4 deletions(-)

diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index d4bfb68fbd3..cce723b38c6 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -111,10 +111,30 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
       singularity::eos_base::PTDerivativesByFiniteDifferences(
           eospac, rho, sie, P, T, derivative_eps, static_cast<Real *>(nullptr), dedP_T_fd,
           drdP_T_fd, dedT_P_fd, drdT_P_fd);
-      REQUIRE(isClose(dedP_T, dedP_T_fd));
-      REQUIRE(isClose(drdP_T, drdP_T_fd));
-      REQUIRE(isClose(dedT_P, dedT_P_fd));
-      REQUIRE(isClose(drdT_P, drdT_P_fd));
+      bool is_close_dedp = isClose(dedP_T, dedP_T_fd);
+      if (!is_close_dedp) {
+        printf("dedP_T not close! %.14e %.14e %.14e\n", dedP_T, dedP_T_fd,
+               dedP_T - dedP_T_fd);
+      }
+      REQUIRE(is_close_dedp);
+      bool is_close_drdP = isClose(drdP_T, drdP_T_fd);
+      if (!is_close_drdP) {
+        printf("drdP_T not close! %.14e %.14e %.14e\n", drdP_T, drdP_T_fd,
+               drdP_T - drdP_T_fd);
+      }
+      REQUIRE(is_close_drdP);
+      bool is_close_dedT = isClose(dedT_P, dedT_P_fd);
+      if (!is_close_dedT) {
+        printf("dedT_P not close! %.14e %.14e %.14e\n", dedT_P, dedT_P_fd,
+               dedT_P - dedT_P_fd);
+      }
+      REQUIRE(is_close_dedT);
+      bool is_close_drdT = isClose(drdT_P, drdT_P_fd);
+      if (!is_close_drdT) {
+        printf("drdT_P not close! %.14e %.14e %.14e\n", drdT_P, drdT_P_fd,
+               drdT_P - drdT_P_fd);
+      }
+      REQUIRE(is_close_drdT);
     }
 
     THEN("The correct metadata is read in") {

From cae1dc71928b5d695219d75adadcc70d5ea51236 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 18:30:47 -0400
Subject: [PATCH 18/31] dx[1] should be dy[0]

---
 singularity-eos/eos/eos_eospac.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 3f4d16f3a16..407b9e961d6 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1743,7 +1743,7 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   table = RofPT_table_;
   eosSafeInterpolate(&table, nxypairs, z, P, T, dx, dy, "RofPT", Verbosity::Quiet);
   drdP_T = dx[0];
-  drdT_P = dx[1];
+  drdT_P = dy[0];
   dedP_T = dedr_T * drdP_T;
   dedT_P = dedT_r + dedr_T * drdT_P;
 #endif // ON DEVICE

From 02bbf9b1aa435947ab7042fe446edb2e152ea7c2 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 18:47:31 -0400
Subject: [PATCH 19/31] debug

---
 singularity-eos/eos/eos_eospac.hpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 407b9e961d6..10e504cb64b 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1744,6 +1744,8 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   eosSafeInterpolate(&table, nxypairs, z, P, T, dx, dy, "RofPT", Verbosity::Quiet);
   drdP_T = dx[0];
   drdT_P = dy[0];
+  printf("dedP_T, dedr_T, drdP_T = %.14e %.14e %.14e\n", // DEBUG
+         dedr_T, dedr_T, drdP_T);
   dedP_T = dedr_T * drdP_T;
   dedT_P = dedT_r + dedr_T * drdT_P;
 #endif // ON DEVICE

From 24090dc2bbda1f7da31736e09f6113939621a39b Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 19:06:13 -0400
Subject: [PATCH 20/31] oops ordering...

---
 singularity-eos/eos/eos_eospac.hpp | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 10e504cb64b..21a8dd817d4 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1736,12 +1736,17 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   EOS_INTEGER table;
 
   table = EofRT_table_;
-  eosSafeInterpolate(&table, nxypairs, z, R, T, dx, dy, "EofRT", Verbosity::Quiet);
+  eosSafeInterpolate(&table, nxypairs, R, T, z dx, dy, "EofRT", Verbosity::Quiet);
   const Real dedr_T = dx[0];
   const Real dedT_r = dy[0];
 
+  table = TofRP_table_;
+  eosSafeInterpolate(&table, nxypairs, R, P, z, dx, dy, "TofRP", Verbosity::Quiet);
+  const Real dTdr_P = dx[0];
+  const Real dTdP_r = dy[0];
+
   table = RofPT_table_;
-  eosSafeInterpolate(&table, nxypairs, z, P, T, dx, dy, "RofPT", Verbosity::Quiet);
+  eosSafeInterpolate(&table, nxypairs, P, T, z, dx, dy, "RofPT", Verbosity::Quiet);
   drdP_T = dx[0];
   drdT_P = dy[0];
   printf("dedP_T, dedr_T, drdP_T = %.14e %.14e %.14e\n", // DEBUG

From 9f88d518b41ada1dbf6b21c03175cdc50af07c41 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 19:27:55 -0400
Subject: [PATCH 21/31] gah missing comma

---
 singularity-eos/eos/eos_eospac.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 21a8dd817d4..745ff3eddfe 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1736,7 +1736,7 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   EOS_INTEGER table;
 
   table = EofRT_table_;
-  eosSafeInterpolate(&table, nxypairs, R, T, z dx, dy, "EofRT", Verbosity::Quiet);
+  eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "EofRT", Verbosity::Quiet);
   const Real dedr_T = dx[0];
   const Real dedT_r = dy[0];
 

From 755fa845daf02255d242aa866b94e6f6c0184331 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 19:48:01 -0400
Subject: [PATCH 22/31] bad printf

---
 singularity-eos/eos/eos_eospac.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 745ff3eddfe..58dbf64a0a9 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1750,7 +1750,7 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   drdP_T = dx[0];
   drdT_P = dy[0];
   printf("dedP_T, dedr_T, drdP_T = %.14e %.14e %.14e\n", // DEBUG
-         dedr_T, dedr_T, drdP_T);
+         dedT_T, dedr_T, drdP_T);
   dedP_T = dedr_T * drdP_T;
   dedT_P = dedT_r + dedr_T * drdT_P;
 #endif // ON DEVICE

From 00298ed18882620a0e4ffb28058bc798320c5f40 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 20:00:50 -0400
Subject: [PATCH 23/31] fix missing unit conversions

---
 singularity-eos/eos/eos_eospac.hpp | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 58dbf64a0a9..1532112d962 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1747,12 +1747,10 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
 
   table = RofPT_table_;
   eosSafeInterpolate(&table, nxypairs, P, T, z, dx, dy, "RofPT", Verbosity::Quiet);
-  drdP_T = dx[0];
-  drdT_P = dy[0];
-  printf("dedP_T, dedr_T, drdP_T = %.14e %.14e %.14e\n", // DEBUG
-         dedT_T, dedr_T, drdP_T);
-  dedP_T = dedr_T * drdP_T;
-  dedT_P = dedT_r + dedr_T * drdT_P;
+  drdP_T = pressureToSesame(dx[0]); // dividing by pressure means inverse conversion
+  drdT_P = temperatueToSesame(dy[0]); // ditto
+  dedP_T = sieFromSesame(pressureToSesame(dedr_T * drdP_T));
+  dedT_P = sieFromSesame(temperatureToSesame(dedT_r + dedr_T * drdT_P));
 #endif // ON DEVICE
 }
 

From c10a4686e806bcb47313da19f28ec89349c6b03c Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 20:08:01 -0400
Subject: [PATCH 24/31] Gah fix units and remove unused interp call

---
 singularity-eos/eos/eos_eospac.hpp | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 1532112d962..d0b98b55349 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1737,20 +1737,15 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
 
   table = EofRT_table_;
   eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "EofRT", Verbosity::Quiet);
-  const Real dedr_T = dx[0];
-  const Real dedT_r = dy[0];
-
-  table = TofRP_table_;
-  eosSafeInterpolate(&table, nxypairs, R, P, z, dx, dy, "TofRP", Verbosity::Quiet);
-  const Real dTdr_P = dx[0];
-  const Real dTdP_r = dy[0];
+  const Real dedr_T = sieFromSesame(dx[0]);
+  const Real dedT_r = sieFromSesame(temperatureToSesame(dy[0]));
 
   table = RofPT_table_;
   eosSafeInterpolate(&table, nxypairs, P, T, z, dx, dy, "RofPT", Verbosity::Quiet);
   drdP_T = pressureToSesame(dx[0]); // dividing by pressure means inverse conversion
-  drdT_P = temperatueToSesame(dy[0]); // ditto
-  dedP_T = sieFromSesame(pressureToSesame(dedr_T * drdP_T));
-  dedT_P = sieFromSesame(temperatureToSesame(dedT_r + dedr_T * drdT_P));
+  drdT_P = temperatureToSesame(dy[0]); // ditto
+  dedP_T = dedr_T * drdP_T;
+  dedT_P = dedT_r + dedr_T * drdT_P;
 #endif // ON DEVICE
 }
 

From 2264aedcbda477b14f33ca3a045202178ad3a443 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 20:31:02 -0400
Subject: [PATCH 25/31] make smaller finite differences in case issue is local
 slope

---
 test/test_eos_tabulated.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index cce723b38c6..c22f9033faa 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -103,7 +103,7 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
       Real P = eospac.PressureFromDensityTemperature(rho, T);
       Real sie = eospac.InternalEnergyFromDensityTemperature(rho, T);
 
-      constexpr Real derivative_eps = 3e-6;
+      constexpr Real derivative_eps = 1e-8;
       Real dedP_T, drdP_T, dedT_P, drdT_P;
       Real dedP_T_fd, drdP_T_fd, dedT_P_fd, drdT_P_fd;
       eospac.PTDerivativesFromPreferred(rho, sie, P, T, static_cast<Real *>(nullptr),

From 0f123ddabd13629d2afaa6cd0ec8b33ccefdf81e Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Mon, 23 Mar 2026 20:33:33 -0400
Subject: [PATCH 26/31] lets also move to a more gentle part of the phase
 space, where rational functions are more likely to agree with linear

---
 test/test_eos_tabulated.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp
index c22f9033faa..e61fb6b8185 100644
--- a/test/test_eos_tabulated.cpp
+++ b/test/test_eos_tabulated.cpp
@@ -98,7 +98,7 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") {
 
     EOS eospac = EOSPAC(steelID);
     THEN("The EOSPAC model can provide derivatives from preferred") {
-      Real rho = 5;
+      Real rho = 50;
       Real T = 500;
       Real P = eospac.PressureFromDensityTemperature(rho, T);
       Real sie = eospac.InternalEnergyFromDensityTemperature(rho, T);

From 655a01b9df7e6ad516193d3129affbf50272dae4 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonah.maxwell.miller@gmail.com>
Date: Wed, 25 Mar 2026 11:57:46 -0400
Subject: [PATCH 27/31] Update singularity-eos/eos/modifiers/eos_unitsystem.hpp

Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com>
---
 singularity-eos/eos/modifiers/eos_unitsystem.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 2a0a53974c5..42c351e957a 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -223,7 +223,7 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
                                   temp * temp_unit_, lambda, dedP_T, drdP_T, dedT_P,
                                   drdT_P);
     // scale outputs
-    dedP_T *= robust::ratio(press_unit_, temp_unit_);
+    dedP_T *= robust::ratio(press_unit_, energy_unit_);
     drdP_T *= robust::ratio(press_unit_, rho_unit_);
     dedT_P *= robust::ratio(temp_unit_, sie_unit_);
     drdT_P *= robust::ratio(temp_unit_, rho_unit_);

From 129f6264e7dc9c85b1f9d640f406b2c5a92036a1 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonah.maxwell.miller@gmail.com>
Date: Wed, 25 Mar 2026 11:59:16 -0400
Subject: [PATCH 28/31] Update doc/sphinx/src/using-eos.rst

Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com>
---
 doc/sphinx/src/using-eos.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index bd33a235236..b196e526f9a 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -1344,9 +1344,9 @@ The function
 
 computes the partial derivatives of density and specific internal
 energy with respect to pressure and temperature, with either pressure
-or temperature fixed. Density, energy, pressure, and temperature, must
-all be passed into this function, to maximize the ability for the
-underlying EOS to perform this calculation performantly. The intended
+or temperature fixed. Each EOS model expects consistent density, 
+energy, pressure and temperature values to be provided so that it
+can perform this calculation performantly. The intended
 use of this method is to compute the cell-averaged thermodynamic
 derivatives in a mixed cell in pressure-temperature equilibrium.
 

From a6ccf422b1754862255e67b3a5b255795879f724 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 25 Mar 2026 14:47:18 -0400
Subject: [PATCH 29/31] jhp suggestions

---
 singularity-eos/eos/eos_eospac.hpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index d0b98b55349..f050a48ac57 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1740,12 +1740,24 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
   const Real dedr_T = sieFromSesame(dx[0]);
   const Real dedT_r = sieFromSesame(temperatureToSesame(dy[0]));
 
+  table = PofRT_table_;
+  eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "PofRT", Verbosity::Quiet);
+  const Real dpdr_T = pressureFromSesame(dx[0]);
+  const Real dpdT_r = pressureFromSesame(temperatureToSesame(dy[0]));
+
+  dedP_T = robust::ratio(dedr_T, dpdr_T);
+  dedT_P = dedT_r - robust::ratio(dedr_T * dpdT_r, dPdr_T);
+  drdP_T = robust::ratio(1., dpdr_T);
+  drdT_P = -robust::ratio(dpdT_r, dpdr_T);
+  /*
+  // Alternative approach
   table = RofPT_table_;
   eosSafeInterpolate(&table, nxypairs, P, T, z, dx, dy, "RofPT", Verbosity::Quiet);
-  drdP_T = pressureToSesame(dx[0]); // dividing by pressure means inverse conversion
-  drdT_P = temperatureToSesame(dy[0]); // ditto
+  drdP_T = pressureToSesame(dx[0]);
+  drdT_P = temperatureToSesame(dy[0]);
   dedP_T = dedr_T * drdP_T;
   dedT_P = dedT_r + dedr_T * drdT_P;
+  */
 #endif // ON DEVICE
 }
 

From 178d51d873dd06c1ae1b7357c0d0fde56306abf4 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 25 Mar 2026 14:57:23 -0400
Subject: [PATCH 30/31] fix that...

---
 singularity-eos/eos/modifiers/eos_unitsystem.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 42c351e957a..a1b1bf4ffc6 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -223,7 +223,7 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
                                   temp * temp_unit_, lambda, dedP_T, drdP_T, dedT_P,
                                   drdT_P);
     // scale outputs
-    dedP_T *= robust::ratio(press_unit_, energy_unit_);
+    dedP_T *= robust::ratio(press_unit_, sie_unit_);
     drdP_T *= robust::ratio(press_unit_, rho_unit_);
     dedT_P *= robust::ratio(temp_unit_, sie_unit_);
     drdT_P *= robust::ratio(temp_unit_, rho_unit_);

From 194dbf7fcd55a0591516267a2636470a84d43728 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 26 Mar 2026 10:14:51 -0400
Subject: [PATCH 31/31] make P consistent capitalization

---
 singularity-eos/eos/eos_eospac.hpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index f050a48ac57..a2292743c70 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -1742,13 +1742,13 @@ EOSPAC::PTDerivativesFromPreferred(const Real rho, const Real sie, const Real pr
 
   table = PofRT_table_;
   eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "PofRT", Verbosity::Quiet);
-  const Real dpdr_T = pressureFromSesame(dx[0]);
-  const Real dpdT_r = pressureFromSesame(temperatureToSesame(dy[0]));
+  const Real dPdr_T = pressureFromSesame(dx[0]);
+  const Real dPdT_r = pressureFromSesame(temperatureToSesame(dy[0]));
 
-  dedP_T = robust::ratio(dedr_T, dpdr_T);
-  dedT_P = dedT_r - robust::ratio(dedr_T * dpdT_r, dPdr_T);
-  drdP_T = robust::ratio(1., dpdr_T);
-  drdT_P = -robust::ratio(dpdT_r, dpdr_T);
+  dedP_T = robust::ratio(dedr_T, dPdr_T);
+  dedT_P = dedT_r - robust::ratio(dedr_T * dPdT_r, dPdr_T);
+  drdP_T = robust::ratio(1., dPdr_T);
+  drdT_P = -robust::ratio(dPdT_r, dPdr_T);
   /*
   // Alternative approach
   table = RofPT_table_;