pt space pte solver

.gitlab/build_and_test.sh

@@ -272,6 +272,7 @@ cmake_test() {
   (
   source ${BUILD_ENV}
   export CTEST_OUTPUT_ON_FAILURE=1
+  export OMP_PROC_BIND=false
   if ${BUILD_WITH_CTEST}; then
     ctest -V -S .gitlab/build_and_test.cmake,Test,$REPORT_ERRORS
   else

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR453]](https://github.com/lanl/singularity-eos/pull/453) A PT space PTE solver
 - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS
 
 ### Fixed (Repair bugs, etc)

doc/sphinx/src/using-closures.rst

@@ -278,6 +278,42 @@ choice of the second independent variable is discussed below and has
 implications for both the number of additional unknowns and the stability of the
 method.
 
+.. _pressure-temperature-formulation:
+The Pressure-Temperature Formulation
+`````````````````````````````````````
+
+An obvious choice is to treat the independent variables as pressure
+and temperature. Then one has only two equations and two unknowns. The
+residual contains only the volume fraction and energy summmation rules
+described above. Taylor expanding these residuals about fixed
+temeprature and pressure points leads to two residual equations of the
+form
+
+.. math::
+
+  0 = (T^* - T) \sum_{i = 0}^{N-1} \left(\frac{\partial f_i}{\partial T}\right)_P + (P^* - P) \sum_{i = 0}^{N-1} \left(\frac{\partial f_i}{\partial P}\right)_T\\
+  0 = (T^* - T) \sum_{i = 0}^{N-1} \left(\frac{\partial \epsilon_i}{\partial T}\right)_P + (P^* - P) \sum_{i = 0}^{N-1} \left(\frac{\partial \epsilon_i}{\partial P}\right)_T
+
+However, derivatives in the volume fraction are not easily
+accessible. They can, however, be cast in terms of the density,
+resulting in the residual:
+
+.. math::
+
+  0 = (T^* - T) \sum_{i = 0}^{N-1} \frac{\rho_{\mathrm{tot}}}{\rho_i^2} \left(\frac{\partial \rho_i}{\partial T}\right)_P + (P^* - P) \sum_{i = 0}^{N-1} \frac{\rho_{\mathrm{tot}}}{\rho_i^2} \left(\frac{\partial \rho_i}{\partial P}\right)_T\\
+  0 = (T^* - T) \sum_{i = 0}^{N-1} \left(\frac{\partial \epsilon_i}{\partial T}\right)_P + (P^* - P) \sum_{i = 0}^{N-1} \left(\frac{\partial \epsilon_i}{\partial P}\right)_T
+
+where :math:`\rho_{\mathrm{tot}}` is the sum of densities over all materials.
+
+The primary advantage of the pressure-temperature space solver is that
+it has only two independent variables and two unknowns, meaning the
+cost scales only linearly with the number of materials, not
+quadratically (or worse). The primary disadvantage, is that most
+equations of state are not formulated in terms of pressure and
+temperature, meaning additional inversions are required.
+
+In the code, this method is referred to as ``PTESolverPT``.
+
 .. _density-energy-formalism:
 The Density-Energy Formulation
 ''''''''''''''''''''''''''''''

singularity-eos/closure/mixed_cell_models.hpp

@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2025. 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
@@ -163,6 +163,9 @@ class CacheAccessor {
   Real *cache_;
 };
 
+// ======================================================================
+// Base class
+// ======================================================================
 template <typename EOSIndexer, typename RealIndexer>
 class PTESolverBase {
  public:
@@ -530,6 +533,9 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU(
   return FastMath::pow2((1.0 - alpha) * lTlo + alpha * lThi);
 }
 
+// ======================================================================
+// PTE Solver RhoT
+// ======================================================================
 inline int PTESolverRhoTRequiredScratch(const std::size_t nmat) {
   std::size_t neq = nmat + 1;
   return neq * neq                 // jacobian
@@ -647,10 +653,9 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
       const Real vf_pert = vfrac[m] + dv;
       const Real rho_pert = robust::ratio(rhobar[m], vf_pert);
 
-      Real p_pert{};
       Real e_pert =
           eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]);
-      p_pert =
+      Real p_pert =
           robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil,
                                                        e_pert, Cache[m], false),
                         uscale);
@@ -757,7 +762,232 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   Real Tequil, Ttemp;
 };
 
+// ======================================================================
+// PT space solver
+// ======================================================================
+inline int PTESolverPTRequiredScratch(const int nmat) {
+  int neq = 2;
+  return neq * neq                 // jacobian
+         + 4 * neq                 // dx, residual, and sol_scratch
+         + 2 * nmat                // all the nmat sized arrays
+         + MAX_NUM_LAMBDAS * nmat; // the cache
+}
+inline size_t PTESolverPTRequiredScratchInBytes(const int nmat) {
+  return PTESolverPTRequiredScratch(nmat) * sizeof(Real);
+}
+// TODO(JMM): make sure normalizations are correct everywhere
+template <typename EOSIndexer, typename RealIndexer, typename LambdaIndexer>
+class PTESolverPT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::InitBase;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::AssignIncrement;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::nmat;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::neq;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::ResidualNorm;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::vfrac_total;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::sie_total;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::eos;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rho;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::vfrac;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::sie;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::temp;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::press;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rho_total;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::uscale;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::utotal_scale;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::MatIndex;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::TryIdealPTE;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::jacobian;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::residual;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::dx;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::u;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rhobar;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::Cache;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::Tnorm;
+  using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::params_;
+
+  enum RES { RV = 0, RSIE = 1 };
+
+ public:
+  // template the ctor to get type deduction/universal references prior to c++17
+  template <typename EOS_t, typename Real_t, typename Lambda_t>
+  PORTABLE_INLINE_FUNCTION
+  PTESolverPT(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real sie_tot,
+              Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, Real_t &&press,
+              Lambda_t &&lambda, Real *scratch, const Real Tnorm = 0.0,
+              const MixParams &params = MixParams())
+      : mix_impl::PTESolverBase<EOSIndexer, RealIndexer>(nmat, 2, eos, vfrac_tot, sie_tot,
+                                                         rho, vfrac, sie, temp, press,
+                                                         scratch, Tnorm, params) {
+    // TODO(JCD): use whatever lambdas are passed in
+    /*for (int m = 0; m < nmat; m++) {
+      if (!variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
+    }*/
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Init() {
+    InitBase();
+    Residual();
+    // calculate an initial equilibrium pressure. Volume-fraction
+    // weighting pressures seems unwise. Use minimum positive pressure
+    // accross initial pressures.
+    Pequil = 1e100; // not actual infinity in case we run with fast math
+    for (int m = 0; m < nmat; ++m) {
+      if ((press[m] < Pequil) && (press[m] > 0)) {
+        Pequil = press[m];
+      }
+    }
+    if (Pequil >= 1e100) Pequil = 1; // note includes uscale
+    // all normalized temps set to 1 so no averaging necessary here.
+    Tequil = 1; // Because it's = Tnorm = initial guess
+    // Leave this in for now, but comment out because I'm not sure it's a good idea
+    // TryIdealPTE(this);
+    // Set the current guess for the equilibrium temperature.  Note that this is already
+    // scaled.
+    return ResidualNorm();
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  void Residual() const {
+    Real vsum = 0.0;
+    Real esum = 0.0;
+    for (int m = 0; m < nmat; ++m) {
+      vsum += vfrac[m];
+      esum += u[m];
+    }
+    residual[RV] = vfrac_total - vsum;
+    residual[RSIE] = utotal_scale - esum;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  bool CheckPTE() const {
+    Real error_v = std::abs(residual[RV]);
+    Real error_u = std::abs(residual[RSIE]);
+    // this may not be quite right as we may need an energy scaling factor
+    // Check for convergence
+    bool converged_u =
+        (error_u < params_.pte_rel_tolerance_e || error_u < params_.pte_abs_tolerance_e);
+    bool converged_v =
+        (error_v < params_.pte_rel_tolerance_e || error_v < params_.pte_abs_tolerance_e);
+    return converged_v && converged_u;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  void Jacobian() const {
+    // sum_m d e_m / dT )_P
+    Real dedT_P_sum = 0.0;
+    // sum_m d e_m / dP )_T
+    Real dedP_T_sum = 0.0;
+    // - sum_m rho_tot / rho_m^2 * d rho_m / dT )_P
+    Real rtor2_dr_dT_P_sum = 0.0;
+    // - sum_m rho_tot / rho_m^2 d rho_m / dP )_T
+    Real rtor2_dr_dP_T_sum = 0.0;
+
+    // JMM: Note rescaling for u, P, T
+
+    // TODO(JMM): Should we use the thermodynamic derivatives rather
+    // than finite differences?
+    for (int m = 0; m < nmat; ++m) {
+      Real r_pert = rho[m]; // provide initial guesses
+      Real e_pert = sie[m];
+
+      //////////////////////////////
+      // perturb pressures
+      //////////////////////////////
+      Real dp = -Pequil * params_.derivative_eps; // always move towards phase transition
+      eos[m].DensityEnergyFromPressureTemperature(uscale * (Pequil + dp), Tnorm * Tequil,
+                                                  Cache[m], r_pert, e_pert);
+      Real drdp = robust::ratio(r_pert - rho[m], dp);
+      // JMM: Note uscale here. Both sides divided by rhobar
+      Real dedp = robust::ratio(rhobar[m] * e_pert - u[m], rhobar[m] * uscale * dp);
+
+      rtor2_dr_dP_T_sum += robust::ratio(rho_total, rho[m] * rho[m]) * drdp;
+      dedP_T_sum += dedp;
+
+      //////////////////////////////
+      // perturb temperatures
+      //////////////////////////////
+      Real dT = Tequil * params_.derivative_eps;
+      eos[m].DensityEnergyFromPressureTemperature(uscale * Pequil, Tnorm * (Tequil + dT),
+                                                  Cache[m], r_pert, e_pert);
+      Real drdT = robust::ratio(r_pert - rho[m], dT);
+      Real dedT = robust::ratio(
+          robust::ratio(e_pert, uscale) - robust::ratio(sie[m], uscale), dT);
+
+      rtor2_dr_dT_P_sum += robust::ratio(rho_total, rho[m] * rho[m]) * drdT;
+      dedT_P_sum += dedT;
+    }
+
+    // Fill in the Jacobian
+    jacobian[0] = -rtor2_dr_dT_P_sum; // TODO(JMM): Check positions
+    jacobian[1] = -rtor2_dr_dP_T_sum;
+    jacobian[2] = dedT_P_sum;
+    jacobian[3] = dedP_T_sum;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real ScaleDx() const {
+    Real scale = 1.0;
+    if (scale * dx[0] < -0.95 * Tequil) {
+      scale = robust::ratio(-0.95 * Tequil, dx[0]);
+    }
+    auto bounded = [=](Real Pbound, Real delta) {
+      return robust::ratio(robust::ratio(Pbound, uscale) - Pequil, delta);
+    };
+
+    for (int m = 0; m < nmat; ++m) {
+      Real Pmin = eos[m].MinimumPressure();
+      scale = std::min(std::abs(scale), std::abs(0.95 * bounded(Pmin, dx[1])));
+    }
+    for (int m = 0; m < nmat; ++m) {
+      Real Ttest = (Tequil + scale * dx[0]) * Tnorm;
+      Real Pmax = eos[m].MaximumPressureAtTemperature(Ttest);
+      scale = std::min(std::abs(scale), std::abs(0.95 * bounded(Pmax, dx[0])));
+    }
+
+    for (int i = 0; i < neq; ++i) {
+      dx[i] *= scale;
+    }
+    return scale;
+  }
+
+  // Update the solution and return new residual.  Possibly called repeatedly with
+  // different scale factors as part of a line search
+  PORTABLE_INLINE_FUNCTION
+  Real TestUpdate(const Real scale, bool const cache_state = false) {
+    if (cache_state) {
+      // Store the current state in temp variables for first iteration of line
+      // search
+      Ttemp = Tequil;
+      Ptemp = Pequil;
+    }
+    Tequil = Ttemp + scale * dx[0];
+    Pequil = Ptemp + scale * dx[1];
+    for (int m = 0; m < nmat; ++m) {
+      eos[m].DensityEnergyFromPressureTemperature(Pequil * uscale, Tequil * Tnorm,
+                                                  Cache[m], rho[m], sie[m]);
+      vfrac[m] = robust::ratio(rhobar[m], rho[m]);
+      u[m] = robust::ratio(sie[m] * rhobar[m], uscale);
+      temp[m] = Tequil;
+    }
+    Residual();
+    return ResidualNorm();
+  }
+
+ private:
+  // TODO(JMM): Should these have trailing underscores?
+  // Current P, T state
+  Real Pequil;
+  Real Tequil;
+  // Scratch states for test update
+  Real Ptemp;
+  Real Ttemp;
+  // TODO(JMM): Should there be a P norm as well as a Tnorm?
+};
+
+// ======================================================================
 // fixed temperature solver
+// ======================================================================
 inline std::size_t PTESolverFixedTRequiredScratch(const std::size_t nmat) {
   std::size_t neq = nmat;
   return neq * neq                 // jacobian
@@ -963,7 +1193,9 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   Real Tequil, Ttemp;
 };
 
+// ======================================================================
 // fixed P solver
+// ======================================================================
 inline std::size_t PTESolverFixedPRequiredScratch(const std::size_t nmat) {
   std::size_t neq = nmat + 1;
   return neq * neq                 // jacobian
@@ -1198,6 +1430,9 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   Real Tequil, Ttemp, Pequil;
 };
 
+// ======================================================================
+// RhoU Solver
+// ======================================================================
 inline std::size_t PTESolverRhoURequiredScratch(const std::size_t nmat) {
   std::size_t neq = 2 * nmat;
   return neq * neq                 // jacobian
@@ -1443,6 +1678,9 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   Real *dpdv, *dtdv, *dpde, *dtde, *vtemp, *utemp;
 };
 
+// ======================================================================
+// Solver loop
+// ======================================================================
 template <class System>
 PORTABLE_INLINE_FUNCTION SolverStatus PTESolver(System &s) {
   SolverStatus status;

test/CMakeLists.txt

@@ -24,6 +24,7 @@ add_executable(
   test_eos_vinet.cpp
   test_eos_mgusup.cpp
   test_eos_powermg.cpp
+  test_pte_ideal.cpp
   test_eos_zsplit.cpp
   )