Merge pull request #394 from LLNL/rc-v2025.06.1

v2025.06.1 Patch Release

Author: mdavis36

RELEASE_NOTES.md

@@ -1,3 +1,15 @@
+Version v2025.06.1 -- Release date 2025-07-21
+==============================================
+  * Important Notes:
+    * This is a patch release for v2025.06.0.
+
+Notable changes include:
+
+  * Bug Fixes / improvements:
+    * The porosity models now restart the full field of initial sound speeds, rather than assuming they are set during problem initialization every time.
+    * Set limits on porosity contributions to damage while crushing out porosity.
+    * Added better limits on ANEOS internal eps(rho, T) interpolation, which seems to have improved ANEOS stability.
+
 Version v2025.06.0 -- Release date 2025-06-18
 ==============================================
   * Important Notes:

cmake/SpheralVersion.cmake

@@ -1 +1 @@
-set(SPHERAL_VERSION 2025.06.0)
+set(SPHERAL_VERSION 2025.06.1)

docs/conf.py

@@ -27,9 +27,9 @@
 author = 'J. Michael Owen'
 
 # The short X.Y version
-version = '2025.06.0'
+version = '2025.06.1'
 # The full version, including alpha/beta/rc tags
-release = '2025.06.0'
+release = '2025.06.1'
 
 
 # -- General configuration ---------------------------------------------------

docs/developer/dev/continuous_deployment.rst

@@ -9,7 +9,7 @@ to their environment.
 
 On LC systems we maintain:
   * ``Spheral/risky``
-  * ``Spheral/2025.06.0``
+  * ``Spheral/2025.06.1``
   * ``Spheral/2025.01.0``
   * ``Spheral/2024.06.1``
   * ``Spheral/2024.01.1``

src/Damage/ProbabilisticDamagePolicy.cc

@@ -279,10 +279,10 @@ update(const KeyType& key,
           const auto alpha0 = (*alpha0Ptr)(i);
           const auto alpha = (*alphaPtr)(i);
           const auto DalphaDti = std::min(0.0, (*DalphaDtPtr)(i));   // Only allowed to grow damage, not reduce it.
-          const auto phi0 = 1.0 - 1.0/alpha0;
-          const auto DD13Dt_p = -FastMath::CubeRootHalley2(phi0)*safeInv(3.0 * pow(1.0 - (alpha - 1.0)*safeInv(alpha0 - 1.0) + Dtiny, 2.0/3.0) * (alpha0 - 1.0))*Dtiny1*DalphaDti;
-          CHECK2(DD13Dt_p >= 0.0, "bad DD13Dt_p: " << DD13Dt_p);
-          D113 = std::min(1.0, D113 + multiplier*DD13Dt_p);
+          const auto phi0_13 = FastMath::CubeRootHalley2(1.0 - 1.0/alpha0);
+          const auto DD13Dt_p = -phi0_13/(3.0 * pow(1.0 - (alpha - 1.0)*safeInv(alpha0 - 1.0) + Dtiny, 2.0/3.0) * (alpha0 - 1.0))*Dtiny1*DalphaDti;
+          CHECK(DD13Dt_p >= 0.0);
+          D113 = std::min(1.0, D113 + std::min(phi0_13, multiplier*DD13Dt_p));
         }
 
         // Increment the damage tensor.

src/Damage/TensorDamagePolicy.cc

@@ -306,10 +306,10 @@ update(const KeyType& key,
           const auto alpha0 = (*alpha0Ptr)(i);
           const auto alpha = (*alphaPtr)(i);
           const auto DalphaDti = std::min(0.0, (*DalphaDtPtr)(i));   // Only allowed to grow damage, not reduce it.
-          const auto phi0 = 1.0 - 1.0/alpha0;
-          const auto DD13Dt_p = -FastMath::CubeRootHalley2(phi0)/(3.0 * pow(1.0 - (alpha - 1.0)*safeInv(alpha0 - 1.0) + Dtiny, 2.0/3.0) * (alpha0 - 1.0))*Dtiny1*DalphaDti;
+          const auto phi0_13 = FastMath::CubeRootHalley2(1.0 - 1.0/alpha0);
+          const auto DD13Dt_p = -phi0_13/(3.0 * pow(1.0 - (alpha - 1.0)*safeInv(alpha0 - 1.0) + Dtiny, 2.0/3.0) * (alpha0 - 1.0))*Dtiny1*DalphaDti;
           CHECK(DD13Dt_p >= 0.0);
-          D113 = std::min(1.0, D113 + multiplier*DD13Dt_p);
+          D113 = std::min(1.0, D113 + std::min(phi0_13, multiplier*DD13Dt_p));
         }
 
         // Increment the damage tensor.

src/Damage/TensorStrainPolicy.cc

@@ -11,6 +11,7 @@
 #include "Hydro/HydroFieldNames.hh"
 #include "DataBase/UpdatePolicyBase.hh"
 #include "DataBase/IncrementState.hh"
+#include "DataBase/IncrementBoundedState.hh"
 #include "DataBase/State.hh"
 #include "DataBase/StateDerivatives.hh"
 #include "Field/Field.hh"
@@ -85,26 +86,35 @@ update(const KeyType& key,
 
   // Check if a porosity model has registered a modifier for the deviatoric stress.
   // They should have added it as a dependency of this policy if so.
-  const auto usePorosity = state.registered(buildKey(SolidFieldNames::fDSjutzi));
+  const auto porosityScaling = state.registered(buildKey(SolidFieldNames::fDSjutzi));
   const Field<Dimension, Scalar>* fDSptr = nullptr;
-  if (usePorosity) {
+  const Field<Dimension, Scalar>* alphaPtr = nullptr;
+  const Field<Dimension, Scalar>* DalphaDtPtr = nullptr;
+  if (porosityScaling) {
     fDSptr = &state.field(buildKey(SolidFieldNames::fDSjutzi), 0.0);
+    alphaPtr = &state.field(buildKey(SolidFieldNames::porosityAlpha), 0.0);
+    DalphaDtPtr = &derivs.field(buildKey(IncrementBoundedState<Dimension, Scalar>::prefix() + SolidFieldNames::porosityAlpha), 0.0);
   }
 
   // Iterate over the internal nodes.
   const auto ni = stateField.numInternalElements();
 #pragma omp parallel for
   for (auto i = 0u; i < ni; ++i) {
-    const auto fDSi = (usePorosity ?
-                       (*fDSptr)(i) :
-                       1.0);
+    double fDSi = 1.0;
 
     // Begin the big bonanza of options!
 
     // PseudoPlasticStrain.
     if (mStrainType == TensorStrainAlgorithm::PseudoPlasticStrain) {
 
-      const auto DSDti = fDSi * DSDt(i);
+      auto DSDti = DSDt(i);
+      if (porosityScaling) {
+        fDSi = (*fDSptr)(i);
+        const auto alphai = (*alphaPtr)(i);
+        const auto DalphaDti = (*DalphaDtPtr)(i);
+        CHECK(alphai >= 1.0);
+        DSDti = (fDSi*DSDti - S(i)*DalphaDti/alphai)/alphai;
+      }
       strain(i) += multiplier*safeInv(mu(i), 1.0e-10)*DSDti;
       stateField(i) = strain(i);
 
@@ -127,7 +137,7 @@ update(const KeyType& key,
       switch(mStrainType) {
 
       case(TensorStrainAlgorithm::BenzAsphaugStrain):
-        CHECK(E(i) >= 0.0);
+        CHECK2(E(i) >= 0.0, "Bad Youngs modulus for " << stateField.nodeList().name() << " " << i << " : " << E(i));
         stateField(i) = (S(i) - P(i)*SymTensor::one)/(E(i) + tiny);
         break;
 

src/Porosity/PalphaPorosity.cc

@@ -158,50 +158,59 @@ evaluateDerivatives(const Scalar time,
     const auto DrhoDti = DrhoDt(i);
     const auto DuDti = DuDt(i);
 
-    // Time evolution of the solid pressure
-    const auto dPsdti = dPsdri*DrhoDti + dPsdui*DuDti;
+    // If fully compressed we're done
+    if (alphai == 1.0) {
 
-    // If we're above the solid transition pressure we should compress out all porosity
-    if (Pi >= mPs) {
-
-      DalphaDt(i) = (1.0 - alphai)*safeInvVar(dt);
+      DalphaDt(i) = 0.0;
       fDSnew(i) = 1.0;
 
     } else {
 
-      // First compute the derivative with respect to pressure
-      auto DalphaDpi = 0.0;
-      if (alphai > 1.0) {
-        if (Pi < mPe or dPsdti < 0.0) {
+      // Time evolution of the solid pressure
+      const auto dPsdti = dPsdri*DrhoDti + dPsdui*DuDti;
 
-          // Elastic
-          if (c0i != mcS0) {  // If initial porous sound speed is the same as solid phase, no elastic evolution
-            const auto halpha = 1.0 + (alphai - 1.0)*(c0i - mcS0)*safeInvVar(mcS0*(mAlphae - 1.0));
-            DalphaDpi = alphai*alphai/mKS0*(1.0 - safeInvVar(halpha*halpha));
-          }
+      // If we're above the solid transition pressure we should compress out all porosity
+      if (Pi >= mPs) {
 
-        } else {
+        DalphaDt(i) = (1.0 - alphai)*safeInvVar(dt);
+        fDSnew(i) = 1.0;
+
+      } else {
+
+        // First compute the derivative with respect to pressure
+        auto DalphaDpi = 0.0;
+        if (alphai > 1.0) {
+          if (Pi < mPe or dPsdti < 0.0) {
+
+            // Elastic
+            if (c0i != mcS0) {  // If initial porous sound speed is the same as solid phase, no elastic evolution
+              const auto halpha = 1.0 + (alphai - 1.0)*(c0i - mcS0)*safeInvVar(mcS0*(mAlphae - 1.0));
+              DalphaDpi = alphai*alphai/mKS0*(1.0 - safeInvVar(halpha*halpha));
+            }
 
-          // Plastic
-          DalphaDpi = (Pi < mPt ?
-                       mn1*(mAlphat - mAlphae)*pow((mPt - Pi)/(mPt - mPe), mn1)*safeInv(mPt - Pi) + mn2*(1.0 - mAlphat)*pow((mPs - Pi)/(mPs - mPe), mn2)*safeInv(mPs - Pi) :
-                       mn2*(1.0 - mAlphat)*pow((mPs - Pi)/(mPs - mPe), mn2)*safeInv(mPs - Pi));
+          } else {
 
+            // Plastic
+            DalphaDpi = (Pi < mPt ?
+                         mn1*(mAlphat - mAlphae)*pow((mPt - Pi)/(mPt - mPe), mn1)*safeInv(mPt - Pi) + mn2*(1.0 - mAlphat)*pow((mPs - Pi)/(mPs - mPe), mn2)*safeInv(mPs - Pi) :
+                         mn2*(1.0 - mAlphat)*pow((mPs - Pi)/(mPs - mPe), mn2)*safeInv(mPs - Pi));
+
+          }
         }
-      }
 
-      // Now we can compute the final time derivative
-      DalphaDpi = min(0.0, DalphaDpi);  // Keep things physical
-      const auto Ainv = safeInvVar(alphai + DalphaDpi*(Pi - rhoi*dPsdri));
-      const auto dPdti = (alphai*dPsdri*DrhoDti + dPsdui*DuDti)*Ainv;
-      DalphaDt(i) = DalphaDpi*dPdti;
+        // Now we can compute the final time derivative
+        DalphaDpi = min(0.0, DalphaDpi);  // Keep things physical
+        const auto Ainv = safeInvVar(alphai + DalphaDpi*(Pi - rhoi*dPsdri));
+        const auto dPdti = (alphai*dPsdri*DrhoDti + dPsdui*DuDti)*Ainv;
+        DalphaDt(i) = DalphaDpi*dPdti;
 
-      // Optionally update the deviatoric stress scaling term
-      if (mJutziStateUpdate) {
-        auto DalphaDrhoi = (Pi/(rhoi*rhoi)*dPsdui + alphai*dPsdri)*Ainv * DalphaDpi;
-        fDSnew(i) = std::max(0.0, std::min(1.0, 1.0 + DalphaDrhoi*rhoi/alphai));
-      } else {
-        fDSnew(i) = 1.0;
+        // Optionally update the deviatoric stress scaling term
+        if (mJutziStateUpdate) {
+          auto DalphaDrhoi = (Pi/(rhoi*rhoi)*dPsdui + alphai*dPsdri)*Ainv * DalphaDpi;
+          fDSnew(i) = std::max(0.0, std::min(1.0, 1.0 + DalphaDrhoi*rhoi/alphai));
+        } else {
+          fDSnew(i) = 1.0;
+        }
       }
     }
 
@@ -237,7 +246,6 @@ void
 PalphaPorosity<Dimension>::
 dumpState(FileIO& file, const string& pathName) const {
   PorosityModel<Dimension>::dumpState(file, pathName);
-  file.write(mc0, pathName + "/c0");
   file.write(mdPdU, pathName + "/dPdU");
   file.write(mdPdR, pathName + "/dPdR");
 }
@@ -250,7 +258,6 @@ void
 PalphaPorosity<Dimension>::
 restoreState(const FileIO& file, const string& pathName) {
   PorosityModel<Dimension>::restoreState(file, pathName);
-  file.read(mc0, pathName + "/c0");
   file.read(mdPdU, pathName + "/dPdU");
   file.read(mdPdR, pathName + "/dPdR");
 }

src/Porosity/PorosityModel.cc

@@ -16,6 +16,7 @@
 #include "DataBase/IncrementState.hh"
 #include "DataBase/IncrementBoundedState.hh"
 #include "DataBase/PureReplaceBoundedState.hh"
+#include "Utilities/SpheralMessage.hh"
 
 #include <string>
 
@@ -282,6 +283,7 @@ dumpState(FileIO& file, const string& pathName) const {
   file.write(mAlpha, pathName + "/alpha");
   file.write(mDalphaDt, pathName + "/DalphaDt");
   file.write(mSolidMassDensity, pathName + "/solidMassDensity");
+  file.write(mc0, pathName + "/c0");
   file.write(mfDS, pathName + "/fDS");
   file.write(mfDSnew, pathName + "/fDSnew");
 }
@@ -298,6 +300,7 @@ restoreState(const FileIO& file, const string& pathName) {
   file.read(mAlpha, pathName + "/alpha");
   file.read(mDalphaDt, pathName + "/DalphaDt");
   file.read(mSolidMassDensity, pathName + "/solidMassDensity");
+  if (file.readIfAvailable(mc0, pathName + "/c0") != 0) SpheralMessage("PorosityModel WARNING: unable to load variable c0");
   file.read(mfDS, pathName + "/fDS");
   file.read(mfDSnew, pathName + "/fDSnew");
 }

src/Porosity/StrainPorosity.cc

@@ -124,45 +124,49 @@ evaluateDerivatives(const Scalar /*time*/,
   const auto n = mNodeList.numInternalNodes();
 #pragma omp parallel for
   for (auto i = 0u; i < n; ++i) {
-    DstrainDt(i) = DvDx(i).Trace();
     const auto alphai = alpha(i);
-    const auto epsi = strain(i);
-    const auto DuDti = DuDt(i);
-    const auto strainRate = min(0.0, DstrainDt(i) - mGammaS0*safeInv(mcS0*mcS0)*DuDti);
-    const auto alpha0i = mAlpha0(i);
-    const auto c0i = mc0(i);
-    const auto csi = mcS0 + (alphai - 1.0)*safeInv(alpha0i - 1.0)*(c0i - mcS0);
-
-    // The time rate of change for alpha depends on if we're in the elastic, exponential, or
-    // power-law regime of compaction.
-    auto DalphaDepsi = 0.0;
-    if (-std::abs(epsi) >= mEpsE) {
-      // Elastic regime
-      // We assume compaction is reversible in this regime.
-      DalphaDepsi = alphai*(1.0 - FastMath::square(csi*safeInv(mcS0)));
-      DalphaDt(i) = DalphaDepsi*strainRate;
-    } else if (epsi < 0.0) { // The plastic regimes only apply for negative strains (compaction)
-      if (epsi >= mEpsX) {
-        // Exponential regime
-        DalphaDepsi = mKappa*alphai;
-        DalphaDt(i) = min(0.0, DalphaDepsi*strainRate);
-      } else {
-        // Power-law -- irreversible compaction.
-        const auto epsCi = 2.0*(1.0 - alpha0i*exp(mKappa*(mEpsX - mEpsE)))/(mKappa*alpha0i*exp(mKappa*(mEpsX - mEpsE))) + mEpsX;
-        const auto PLcoefi = 2.0*(1.0 - alpha0i*exp(mKappa*(mEpsX - mEpsE)))/FastMath::square(epsCi - mEpsX);
-        DalphaDepsi = PLcoefi*(mEpsE - epsi);
-        DalphaDt(i) = min(0.0, DalphaDepsi*strainRate);
+    if (alphai == 1.0) {
+      DstrainDt(i) = 0.0;
+      DalphaDt(i) = 0.0;
+      fDSnew(i) = 1.0;
+    } else {
+      DstrainDt(i) = DvDx(i).Trace();
+      const auto epsi = strain(i);
+      const auto DuDti = DuDt(i);
+      const auto strainRate = min(0.0, DstrainDt(i) - mGammaS0*safeInv(mcS0*mcS0)*DuDti);
+      const auto alpha0i = mAlpha0(i);
+      const auto c0i = mc0(i);
+      const auto csi = mcS0 + (alphai - 1.0)*safeInv(alpha0i - 1.0)*(c0i - mcS0);
+
+      // The time rate of change for alpha depends on if we're in the elastic, exponential, or
+      // power-law regime of compaction.
+      auto DalphaDepsi = 0.0;
+      if (-std::abs(epsi) >= mEpsE) {
+        // Elastic regime
+        // We assume compaction is reversible in this regime.
+        DalphaDepsi = alphai*(1.0 - FastMath::square(csi*safeInv(mcS0)));
+        DalphaDt(i) = DalphaDepsi*strainRate;
+      } else if (epsi < 0.0) { // The plastic regimes only apply for negative strains (compaction)
+        if (epsi >= mEpsX) {
+          // Exponential regime
+          DalphaDepsi = mKappa*alphai;
+          DalphaDt(i) = min(0.0, DalphaDepsi*strainRate);
+        } else {
+          // Power-law -- irreversible compaction.
+          const auto epsCi = 2.0*(1.0 - alpha0i*exp(mKappa*(mEpsX - mEpsE)))/(mKappa*alpha0i*exp(mKappa*(mEpsX - mEpsE))) + mEpsX;
+          const auto PLcoefi = 2.0*(1.0 - alpha0i*exp(mKappa*(mEpsX - mEpsE)))/FastMath::square(epsCi - mEpsX);
+          DalphaDepsi = PLcoefi*(mEpsE - epsi);
+          DalphaDt(i) = min(0.0, DalphaDepsi*strainRate);
+        }
       }
-    }
 
-    // Optionally update the deviatoric stress scaling term
-    if (mJutziStateUpdate) {
-      fDSnew(i) = (alphai == 1.0 ?
-                   1.0 :
-                   1.0 - DalphaDepsi/rhoS(i));
-                   // 1.0 + DalphaDt(i)/alphai * safeInv(DrhoDt(i)/rho(i), 0.01*csi*Dimension::nDim/H(i).Trace()));
-    } else {
-      fDSnew(i) = 1.0;
+      // Optionally update the deviatoric stress scaling term
+      if (mJutziStateUpdate) {
+        fDSnew(i) = std::max(0.0, std::min(1.0, 1.0 - DalphaDepsi/rhoS(i)));
+        // 1.0 + DalphaDt(i)/alphai * safeInv(DrhoDt(i)/rho(i), 0.01*csi*Dimension::nDim/H(i).Trace()));
+      } else {
+        fDSnew(i) = 1.0;
+      }
     }
 
 #pragma omp critical