pt derivatives from preferred

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
 
 ### Fixed (Repair bugs, etc)
 

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. 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.
+
+.. 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*>

singularity-eos/eos/eos_base.hpp

@@ -206,6 +206,33 @@ 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 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);
+  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.
@@ -288,6 +315,23 @@ 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 {
+    constexpr Real derivative_eps = 3.0e-6;
+    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
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void

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
@@ -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,52 @@ 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 = EofRT_table_;
+  eosSafeInterpolate(&table, nxypairs, R, T, z, dx, dy, "EofRT", Verbosity::Quiet);
+  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]);
+  drdT_P = temperatureToSesame(dy[0]);
+  dedP_T = dedr_T * drdP_T;
+  dedT_P = dedT_r + dedr_T * drdT_P;
+  */
+#endif // ON DEVICE
+}
+
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 EOSPAC::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,

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,

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,

singularity-eos/eos/eos_spiner_rho_temp.hpp

@@ -196,6 +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,
@@ -1005,6 +1013,42 @@ 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 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?
+  } 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,

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
@@ -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,

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_, 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_);
+  }
+
   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,

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,

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,

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,

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, inv_scale_ * 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,

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,