Fix Bug in PTE Solver and Several Code Quality Improvements

CHANGELOG.md

@@ -26,6 +26,7 @@
 - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface
 - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions
 - [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust
+- [[PR383]](https://github.com/lanl/singularity-eos/pull/383) Fix bug in step scaling for PTE solver
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls

sesame2spiner/examples/unit_tests/copper.dat

@@ -0,0 +1,20 @@
+# ======================================================================
+#  Example input deck for sesame2spiner,
+#  a tool for converting eospac to spiner
+#  Author: Jonah Miller ([email protected])
+#  © 2021-2023. 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
+#  Nuclear Security Administration. All rights in the program are
+#  reserved by Triad National Security, LLC, and the U.S. Department of
+#  Energy/National Nuclear Security Administration. The Government is
+#  granted for itself and others acting on its behalf a nonexclusive,
+#  paid-up, irrevocable worldwide license in this material to reproduce,
+#  prepare derivative works, distribute copies to the public, perform
+#  publicly and display publicly, and to permit others to do so.
+# ======================================================================
+
+matid=3337
+name=copper
+rhomin=0.

singularity-eos/closure/mixed_cell_models.hpp

@@ -170,9 +170,9 @@ class PTESolverBase {
   PORTABLE_INLINE_FUNCTION
   virtual void Fixup() const {
     Real vsum = 0;
-    for (int m = 0; m < nmat; ++m)
+    for (int m = 0; m < nmat; ++m) {
       vsum += vfrac[m];
-    PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive");
+    }
     for (int m = 0; m < nmat; ++m)
       vfrac[m] *= robust::ratio(vfrac_total, vsum);
   }
@@ -182,7 +182,6 @@ class PTESolverBase {
   void Finalize() {
     for (int m = 0; m < nmat; m++) {
       temp[m] *= Tnorm;
-      PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned");
       u[m] *= uscale;
       press[m] *= uscale;
     }
@@ -219,8 +218,9 @@ class PTESolverBase {
     // material m averaged over the full PTE volume
     rho_total = 0.0;
     for (int m = 0; m < nmat; ++m) {
-      PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction");
-      PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density");
+      PORTABLE_REQUIRE(vfrac[m] > 0.,
+                       "Non-positive volume fraction provided to PTE solver");
+      PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density provided to PTE solver");
       rhobar[m] = rho[m] * vfrac[m];
       rho_total += rhobar[m];
     }
@@ -353,12 +353,12 @@ class PTESolverBase {
       Asum += vfrac[m] * robust::ratio(press[m], temp[m]);
       rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]);
     }
-    PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess");
+    PORTABLE_REQUIRE(Tnorm > 0., "Non-positive Ideal PTE temperature guess");
     Asum *= uscale / Tnorm;
     rhoBsum /= Tnorm;
-    PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density");
+    PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive Ideal PTE energy density guess");
     Tideal = uscale / rhoBsum / Tnorm;
-    PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum");
+    PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive Ideal PTE volume fraction sum");
     Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total);
   }
 
@@ -671,6 +671,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   PORTABLE_INLINE_FUNCTION
   Real ScaleDx() const {
     using namespace mix_params;
+    // Each check reduces the scale further if necessary
     Real scale = 1.0;
     // control how big of a step toward vfrac = 0 is allowed
     for (int m = 0; m < nmat; ++m) {
@@ -684,8 +685,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
       const Real rho_min =
           std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew));
       const Real alpha_max = robust::ratio(rhobar[m], rho_min);
+      if (alpha_max < vfrac[m]) {
+        // Despite our best efforts, we're already in the unstable regime (i.e.
+        // dPdV_T > 0) so we would actually want to *increase* the step instead
+        // of decreasing it. As a result, this code doesn't work as intended and
+        // should be skipped.
+        continue;
+      }
       if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) {
-        scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]);
+        scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]);
       }
     }
     // control how big of a step toward T = 0 is allowed
@@ -701,8 +709,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   // 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) {
-    if (scale == 1.0) {
+  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;
       for (int m = 0; m < nmat; ++m)
         vtemp[m] = vfrac[m];
@@ -879,6 +889,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   PORTABLE_INLINE_FUNCTION
   Real ScaleDx() const {
     using namespace mix_params;
+    // Each check reduces the scale further if necessary
     Real scale = 1.0;
     // control how big of a step toward vfrac = 0 is allowed
     for (int m = 0; m < nmat; ++m) {
@@ -890,6 +901,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
     for (int m = 0; m < nmat; m++) {
       const Real rho_min = eos[m].RhoPmin(Tequil);
       const Real alpha_max = robust::ratio(rhobar[m], rho_min);
+      if (alpha_max < vfrac[m]) {
+        // Despite our best efforts, we're already in the unstable regime (i.e.
+        // dPdV_T > 0) so we would actually want to *increase* the step instead
+        // of decreasing it. As a result, this code doesn't work as intended and
+        // should be skipped.
+        continue;
+      }
       if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) {
         scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]);
       }
@@ -903,8 +921,10 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   // 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) {
-    if (scale == 1.0) {
+  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
       for (int m = 0; m < nmat; ++m)
         vtemp[m] = vfrac[m];
     }
@@ -1093,6 +1113,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   PORTABLE_INLINE_FUNCTION
   Real ScaleDx() const {
     using namespace mix_params;
+    // Each check reduces the scale further if necessary
     Real scale = 1.0;
     // control how big of a step toward vfrac = 0 is allowed
     for (int m = 0; m < nmat; ++m) {
@@ -1106,6 +1127,13 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
       const Real rho_min =
           std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew));
       const Real alpha_max = robust::ratio(rhobar[m], rho_min);
+      if (alpha_max < vfrac[m]) {
+        // Despite our best efforts, we're already in the unstable regime (i.e.
+        // dPdV_T > 0) so we would actually want to *increase* the step instead
+        // of decreasing it. As a result, this code doesn't work as intended and
+        // should be skipped.
+        continue;
+      }
       if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) {
         scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]);
       }
@@ -1123,8 +1151,10 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
   // 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) {
-    if (scale == 1.0) {
+  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;
       for (int m = 0; m < nmat; ++m)
         vtemp[m] = vfrac[m];
@@ -1332,6 +1362,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   PORTABLE_INLINE_FUNCTION
   Real ScaleDx() const {
     using namespace mix_params;
+    // Each check reduces the scale further if necessary
     Real scale = 1.0;
     for (int m = 0; m < nmat; ++m) {
       // control how big of a step toward vfrac = 0 is allowed
@@ -1348,6 +1379,13 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
           std::max(eos[m].RhoPmin(Tnorm * temp[m]), eos[m].RhoPmin(Tnorm * tt));
       const Real alpha_max = robust::ratio(rhobar[m], rho_min);
       // control how big of a step toward rho = rho(Pmin) is allowed
+      if (alpha_max < vfrac[m]) {
+        // Despite our best efforts, we're already in the unstable regime (i.e.
+        // dPdV_T > 0) so we would actually want to *increase* the step instead
+        // of decreasing it. As a result, this code doesn't work as intended and
+        // should be skipped.
+        continue;
+      }
       if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) {
         scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]);
       }
@@ -1361,8 +1399,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
   // 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) const {
-    if (scale == 1.0) {
+  Real TestUpdate(const Real scale, bool const cache_state = false) const {
+    if (cache_state) {
+      // Store the current state in temp variables for first iteration of line
+      // search
       for (int m = 0; m < nmat; ++m) {
         vtemp[m] = vfrac[m];
         utemp[m] = u[m];
@@ -1417,25 +1457,26 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) {
 
     // possibly scale the update to stay within reasonable bounds
     Real scale = s.ScaleDx();
-    // const Real scale_save = scale;
+    PORTABLE_REQUIRE(scale <= 1.0, "PTE Solver is attempting to increase the step size");
 
     // Line search
     Real gradfdx = -2.0 * scale * err;
-    scale = 1.0;
+    scale = 1.0; // New scale for line search
     Real err_old = err;
-    err = s.TestUpdate(scale);
+    // Test the update and reset the cache the current state
+    err = s.TestUpdate(scale, true /* cache_state */);
     if (err > err_old + line_search_alpha * gradfdx) {
-      // backtrack
-      Real err_mid = s.TestUpdate(0.5);
+      // backtrack to middle of step
+      scale = 0.5;
+      Real err_mid = s.TestUpdate(scale);
       if (err_mid < err && err_mid < err_old) {
+        // try a larger step since a half step reduces the error
         scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old);
-      } else {
-        scale = line_search_fac;
       }
-
       for (int line_iter = 0; line_iter < line_search_max_iter; line_iter++) {
         err = s.TestUpdate(scale);
         if (err < err_old + line_search_alpha * scale * gradfdx) break;
+        // shrink the step if the error isn't reduced
         scale *= line_search_fac;
       }
     }

singularity-eos/eos/get_sg_eos_functors.hpp

@@ -63,32 +63,37 @@ struct init_functor {
 
   PORTABLE_INLINE_FUNCTION
   void operator()(const int i, const int tid, double &mass_sum, int &npte,
-                  const Real t_mult, const Real s_mult, const Real p_mult) const {
+                  double &vfrac_sum, const Real t_mult, const Real s_mult,
+                  const Real p_mult) const {
     /* normalize mass fractions */
     /* first find the mass sum */
     /* also set idxs as the decrement of the eos offsets */
     /* to take into account 1 based indexing in fortran */
     for (int m = 0; m < nmat; ++m) {
       mass_sum += frac_mass_v(i, m);
       pte_idxs(tid, m) = eos_offsets_v(m) - 1;
-      frac_vol_v(i, m) = 0.0;
+      // Assume non-participating materials are all at cell density. PTE volumes
+      // will be overwritten
+      frac_vol_v(i, m) = frac_mass_v(i, m);
+      // Initialize all materials to the cell sie (if provided). PTE sie values
+      // will be overwritten. This also means that the sie for the PTE
+      // materials is the same as the cell PTE after removing small materials
+      frac_ie_v(i, m) = sie_v(i) * frac_mass_v(i, m) * s_mult;
     }
     for (int m = 0; m < nmat; ++m) {
       frac_mass_v(i, m) /= mass_sum;
     }
     check_all_vals(i);
-    // count the number of participating materials and zero the inputs
+    // count the number of participating materials
     npte = 0;
     for (int m = 0; m < nmat; ++m) {
       if (frac_mass_v(i, m) > mass_frac_cutoff) {
         // participating materials are those with non-negligible mass fractions
         pte_idxs(tid, npte) = eos_offsets_v(m) - 1;
         pte_mats(tid, npte) = m;
         npte += 1;
-      } else {
-        frac_ie_v(i, m) = 0.0;
       }
-      // zero the inputs
+      // zero the PTE solver material inputs
       vfrac_pte(tid, m) = 0.0;
       sie_pte(tid, m) = 0.0;
       temp_pte(tid, m) = 0.0;
@@ -98,10 +103,12 @@ struct init_functor {
     // NOTE: the volume fractions and densities need to be consistent with the
     // total specific volume since they are used to calculate internal
     // quantities for the PTE solver
+    vfrac_sum = 0.;
     for (int mp = 0; mp < npte; ++mp) {
       const int m = pte_mats(tid, mp);
       // Need to guess volume fractions
       vfrac_pte(tid, mp) = frac_mass_v(i, m);
+      vfrac_sum += vfrac_pte(tid, mp);
       // Calculate densities to be consistent with these volume fractions
       rho_pte(tid, mp) = frac_mass_v(i, m) / spvol_v(i) / vfrac_pte(tid, mp);
       temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult;

singularity-eos/eos/get_sg_eos_rho_e.cpp

@@ -37,8 +37,9 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v,
         const int32_t tid{small_loop ? iloop : token};
         double mass_sum{0.0};
         int npte{0};
+        double vfrac_sum{0.0};
         // initialize values for solver / lookup
-        i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0);
+        i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 1.0, 0.0);
         // need to initialize the scratch before it's used to avoid undefined behavior
         for (int idx = 0; idx < solver_scratch.extent(1); ++idx) {
           solver_scratch(tid, idx) = 0.0;
@@ -54,7 +55,7 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v,
           singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
           // reset inputs
           PTESolverRhoT<singularity::EOSAccessor_, Real *, Real **> method(
-              npte, eos_inx, 1.0, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
+              npte, eos_inx, vfrac_sum, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
               &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache,
               &solver_scratch(tid, 0));
           pte_converged = PTESolver(method);

singularity-eos/eos/get_sg_eos_rho_p.cpp

@@ -38,8 +38,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
         const int32_t tid{small_loop ? iloop : token};
         double mass_sum{0.0};
         int npte{0};
+        double vfrac_sum{0.0};
         // initialize values for solver / lookup
-        i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0);
+        i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 0.0, 1.0);
         Real sie_tot_true{0.0};
         // need to initialize the scratch before it's used to avoid undefined behavior
         for (int idx = 0; idx < solver_scratch.extent(1); ++idx) {
@@ -54,9 +55,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
           // eos accessor
           singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
           PTESolverFixedP<singularity::EOSAccessor_, Real *, Real *> method(
-              npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
-              &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0],
-              &solver_scratch(tid, 0));
+              npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0),
+              &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0),
+              cache[0], &solver_scratch(tid, 0));
           pte_converged = PTESolver(method);
           // calculate total sie
           for (int mp = 0; mp < npte; ++mp) {

singularity-eos/eos/get_sg_eos_rho_t.cpp

@@ -38,8 +38,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v,
         const int32_t tid{small_loop ? iloop : token};
         double mass_sum{0.0};
         int npte{0};
+        double vfrac_sum{0.0};
         // initialize values for solver / lookup
-        i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0);
+        i_func(i, tid, mass_sum, npte, vfrac_sum, 1.0, 0.0, 0.0);
         // calculate pte condition (lookup for 1 mat cell)
         Real sie_tot_true{0.0};
         // need to initialize the scratch before it's used to avoid undefined behavior
@@ -55,9 +56,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v,
           // eos accessor
           singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
           PTESolverFixedT<singularity::EOSAccessor_, Real *, Real **> method(
-              npte, eos_inx, 1.0, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
-              &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache,
-              &solver_scratch(tid, 0));
+              npte, eos_inx, vfrac_sum, temp_pte(tid, 0), &rho_pte(tid, 0),
+              &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0),
+              cache, &solver_scratch(tid, 0));
           pte_converged = PTESolver(method);
           // calculate total internal energy
           for (int mp = 0; mp < npte; ++mp) {

spack-repo/packages/singularity-eos/package.py

@@ -112,6 +112,7 @@ class SingularityEos(CMakePackage, CudaPackage):
 
     depends_on("[email protected],1.5.2:", when="@:1.7.0")
     depends_on("[email protected]:", when="@1.7.1:")
+    depends_on("[email protected]:", when="@1.8.1:")
     # request HEAD of main branch
     depends_on("ports-of-call@main", when="@main")