Merge pull request #321 from LLNL/bugfix/PorosityDamage

Fix for porosity interaction with damage when starting with zero porosity
LLNL · Jan 3, 2025 · 69b1459 · 69b1459
2 parents 9a039f7 + bfb5c1d
commit 69b1459
Show file tree

Hide file tree

Showing 7 changed files with 106 additions and 77 deletions.
diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md
@@ -58,6 +58,7 @@ Notable changes include:
     * Bugfix for RZ solid CRKSPH with compatible energy.
     * Parsing of None string now always becomes None python type. Tests have been updated accordingly.
     * IO for checkpoints and visuzalization can now be properly turned off through SpheralController input options.
+    * Fixed porosity model interaction with damage for zero porosity case.
 
 Version v2024.06.1 -- Release date 2024-07-09
 ==============================================

diff --git a/src/Damage/ProbabilisticDamagePolicy.cc b/src/Damage/ProbabilisticDamagePolicy.cc
@@ -288,8 +288,8 @@ update(const KeyType& key,
           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;
-          CHECK(DD13Dt_p >= 0.0);
+          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);
         }
 

diff --git a/src/Porosity/ShadowPalphaPorosity.py b/src/Porosity/ShadowPalphaPorosity.py
@@ -79,7 +79,7 @@ def __init__(self,
                     last_alphae = 0.0
                     iter = 0
                     def DalphaDP_elastic(P, alpha):
-                        h = 1.0 + (alpha - 1.0)*(c0min - cS0)/(cS0*(alphae - 1.0))
+                        h = 1.0 + (alpha - 1.0)*(c0min - cS0)/max(1.0e-20, cS0*(alphae - 1.0))
                         return alpha*alpha/K0*(1.0 - 1.0/(h*h))
                     while abs(alphae - last_alphae) > 1.0e-10 and iter < 1000:
                         iter += 1

diff --git a/src/Porosity/computeHugoniotWithPorosity.py b/src/Porosity/computeHugoniotWithPorosity.py
@@ -76,15 +76,16 @@ def __init__(self,
 
         # Find alphae = alpha(Pe)
         self.alphae = alpha0   # Starting point
-        last_alphae = 0.0
-        iter = 0
-        while abs(self.alphae - last_alphae) > 1.0e-15 and iter < 1000:
-            iter += 1
-            last_alphae = self.alphae
-            self.alphae = scipy.integrate.solve_ivp(self.Dalpha_elasticDP,
-                                                    t_span = [self.P0, self.Pe],
-                                                    y0 = [self.alpha0],
-                                                    t_eval = [self.Pe]).y[0][0]
+        if alpha0 > 1.0:
+            last_alphae = 0.0
+            iter = 0
+            while abs(self.alphae - last_alphae) > 1.0e-15 and iter < 1000:
+                iter += 1
+                last_alphae = self.alphae
+                self.alphae = scipy.integrate.solve_ivp(self.Dalpha_elasticDP,
+                                                        t_span = [self.P0, self.Pe],
+                                                        y0 = [self.alpha0],
+                                                        t_eval = [self.Pe]).y[0][0]
 
         if self.alphat is None:
             self.alphat = self.alphae    # Reduces to Eq 8 in Jutzi 2008
@@ -96,8 +97,8 @@ def __init__(self,
         return
 
     def h(self, alpha):
-        assert self.alphae > 1.0 and self.c0 < self.cS0, "alphae={}, c0={}, cS0={}".format(self.alphae, self.c0, self.cS0)
-        return 1.0 + (alpha - 1.0)*(self.c0 - self.cS0)/(self.cS0*(self.alphae - 1.0))
+        assert self.alphae >= 1.0 and self.c0 <= self.cS0, "alphae={}, c0={}, cS0={}".format(self.alphae, self.c0, self.cS0)
+        return 1.0 + (alpha - 1.0)*(self.c0 - self.cS0)/max(1.0e-20, self.cS0*(self.alphae - 1.0))
 
     def Dalpha_elasticDP(self, P, alpha):
         return alpha*alpha/self.K0*(1.0 - 1.0/self.h(alpha)**2)

diff --git a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py
@@ -37,12 +37,12 @@
 # MFM
 #
 #ATS:mfm0 = test(         SELF, "--mfm True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories True --steps 100", label="Noh cylindrical MFM, nPerh=2.5", np=8, gsph=True)
-#ATS:mfm1 = testif(gsph0, SELF, "--mfm True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical MFM, nPerh=2.5, restart test", np=8, gsph=True)
+#ATS:mfm1 = testif(mfm0,  SELF, "--mfm True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical MFM, nPerh=2.5, restart test", np=8, gsph=True)
 #
 # MFV
 #
 #ATS:mfv0 = test(         SELF, "--mfv True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories True --steps 100", label="Noh cylindrical MFV, nPerh=2.5", np=8, gsph=True)
-#ATS:mfv1 = testif(gsph0, SELF, "--mfv True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical MFV, nPerh=2.5, restart test", np=8, gsph=True)
+#ATS:mfv1 = testif(mfv0,  SELF, "--mfv True --nRadial 100 --cfl 0.25 --nPerh 2.51 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical MFV, nPerh=2.5, restart test", np=8, gsph=True)
 
 #-------------------------------------------------------------------------------
 # The Cylindrical Noh test case run in 2-D.

diff --git a/tests/functional/Porosity/PlanarCompaction/PlanarCompaction-1d.py b/tests/functional/Porosity/PlanarCompaction/PlanarCompaction-1d.py
@@ -13,6 +13,11 @@
 #ATS:t1 = test(      SELF, "--graphics False --clearDirectories True  --checkError False --dataDirBase dumps-PlanarCompaction-1d-sph-restart --restartStep 100 --steps 200", label="Planar porous aluminum compaction problem -- 1-D (serial, restart test step 1)")
 #ATS:t2 = testif(t1, SELF, "--graphics False --clearDirectories False --checkError False --dataDirBase dumps-PlanarCompaction-1d-sph-restart --restartStep 100 --steps 100 --checkRestart True --restoreCycle 100 --postCleanup True", label="Planar porous aluminum compaction problem -- 1-D (serial, restart test step 2)")
 #
+# Ordinary SPH with no initial porosity
+#
+#ATS:t3 = testif(t0, SELF, "--alpha0 1.0 --useDamage False --goalTime 1.0 --graphics False --clearDirectories True  --checkError True  --dataDirBase dumps-PlanarCompaction-1d-sph --restartStep 100000 --postCleanup True", np=4, label="Planar porous aluminum compaction problem with no initial porosity -- 1-D (4 proc, no damage model)")
+#ATS:t4 = testif(t3, SELF, "--alpha0 1.0 --useDamage True  --goalTime 1.0 --graphics False --clearDirectories True  --checkError True  --dataDirBase dumps-PlanarCompaction-1d-sph --restartStep 100000 --postCleanup True", np=4, label="Planar porous aluminum compaction problem with no initial porosity -- 1-D (4 proc, with damage model)")
+#
 # FSISPH
 #
 #ATS:t10 = test(       SELF, "--graphics False --clearDirectories True  --checkError True  --hydroType FSISPH --dataDirBase dumps-PlanarCompaction-1d-fsisph --restartStep 100000 --postCleanup True", np=4, label="Planar porous aluminum compaction problem -- 1-D (FSISPH, 4 proc)", fsisph=True)
@@ -131,62 +136,81 @@
 #-------------------------------------------------------------------------------
 # The reference values for error norms checking for pass/fail
 #-------------------------------------------------------------------------------
-LnormRef = {"SPH": {"Mass density" : {"L1"   : 0.06784186300927694,
-                                      "L2"   : 0.012774373437299643,
-                                      "Linf" : 0.6245679444354701},
-                    "Spec Therm E" : {"L1"   : 0.0001200742460791407,
-                                      "L2"   : 2.2616105613742583e-05,
-                                      "Linf" : 0.0010923440797387786},
-                    "velocity    " : {"L1"   : 0.004921931042558655,
-                                      "L2"   : 0.0009173594117158436,
-                                      "Linf" : 0.0448725433453345},
-                    "pressure    " : {"L1"   : 0.0022217375280911347,
-                                      "L2"   : 0.00039479153550769805,
-                                      "Linf" : 0.018793913196205617},
-                    "alpha       " : {"L1"   : 0.0590391542763204,
-                                      "L2"   : 0.007963583413760916,
-                                      "Linf" : 0.2738180402369801},
-                    "h           " : {"L1"   : 0.00043261838627472803,
-                                      "L2"   : 8.062946952637553e-05,
-                                      "Linf" : 0.014201309070925212}},
-
-            "FSISPH": {"Mass density" : {"L1"   : 0.06781314493410028,
-                                         "L2"   : 0.012767602580471844,
-                                         "Linf" : 0.6245698198195724},
-                       "Spec Therm E" : {"L1"   : 0.00011965457154529887,
-                                         "L2"   : 2.254867764655585e-05,
-                                         "Linf" : 0.0010923441579020216},
-                       "velocity    " : {"L1"   : 0.004913235403786584,
-                                         "L2"   : 0.0009163181868064007,
-                                         "Linf" : 0.044872645505280966},
-                       "pressure    " : {"L1"   : 0.002210222289968727,
-                                         "L2"   : 0.00039387994606202237,
-                                         "Linf" : 0.018793847854203908},
-                       "alpha       " : {"L1"   : 0.05903530352469833,
-                                         "L2"   : 0.007960855343380124,
-                                         "Linf" : 0.2738175996911776},
-                       "h           " : {"L1"   : 0.0004319674641541397,
-                                         "L2"   : 8.05536967933465e-05,
-                                         "Linf" : 0.014201309071287523}},
-
-            "CRKSPH": {"Mass density" : {"L1"   : 0.0679707220773017,
-                                         "L2"   : 0.012783621322270034,
-                                         "Linf" : 0.6245687018640083},
-                       "Spec Therm E" : {"L1"   : 0.0001201303520771047,
-                                         "L2"   : 2.2622550331357265e-05,
-                                         "Linf" : 0.00109234444748564},
-                       "velocity    " : {"L1"   : 0.004926121368235753,
-                                         "L2"   : 0.0009173629106343205,
-                                         "Linf" : 0.044872623201390446},
-                       "pressure    " : {"L1"   : 0.0022258225868044238,
-                                         "L2"   : 0.00039496818856079414,
-                                         "Linf" : 0.018793890216913627},
-                       "alpha       " : {"L1"   : 0.05909030710035451,
-                                         "L2"   : 0.007965976598237856,
-                                         "Linf" : 0.2738173870953471},
-                       "h           " : {"L1"   : 0.00043274827065262975,
-                                         "L2"   : 8.062817025125132e-05,
-                                         "Linf" : 0.014201309070451522}},
+LnormRef = {("SPH", 1.275) : {"Mass density" : {"L1"   : 0.06784186300927694,
+                                                "L2"   : 0.012774373437299643,
+                                                "Linf" : 0.6245679444354701},
+                              "Spec Therm E" : {"L1"   : 0.0001200742460791407,
+                                                "L2"   : 2.2616105613742583e-05,
+                                                "Linf" : 0.0010923440797387786},
+                              "velocity    " : {"L1"   : 0.004921931042558655,
+                                                "L2"   : 0.0009173594117158436,
+                                                "Linf" : 0.0448725433453345},
+                              "pressure    " : {"L1"   : 0.0022217375280911347,
+                                                "L2"   : 0.00039479153550769805,
+                                                "Linf" : 0.018793913196205617},
+                              "alpha       " : {"L1"   : 0.0590391542763204,
+                                                "L2"   : 0.007963583413760916,
+                                                "Linf" : 0.2738180402369801},
+                              "h           " : {"L1"   : 0.00043261838627472803,
+                                                "L2"   : 8.062946952637553e-05,
+                                                "Linf" : 0.014201309070925212}},
+
+            ("SPH", 1.0) : {"Mass density" : {"L1"   : 0.005403187072834507,  
+                                              "L2"   : 0.001992916969258407,  
+                                              "Linf" : 0.22292338017813007},  
+                            "Spec Therm E" : {"L1"   : 2.9770274733200942e-05,
+                                              "L2"   : 1.038711319331483e-05, 
+                                              "Linf" : 0.0010487495498077324},
+                            "velocity    " : {"L1"   : 0.001085888733217598,  
+                                              "L2"   : 0.00040460303212683284,
+                                              "Linf" : 0.04537793164974422},  
+                            "pressure    " : {"L1"   : 0.0018074730403377138, 
+                                              "L2"   : 0.0006629679859824688, 
+                                              "Linf" : 0.0730479929478202},   
+                            "alpha       " : {"L1"   : 0.0,                   
+                                              "L2"   : 0.0,                   
+                                              "Linf" : 0.0},                  
+                            "h           " : {"L1"   : 6.116234877070288e-05, 
+                                              "L2"   : 3.0869281685704505e-05,
+                                              "Linf" : 0.014204020975568329}},
+
+            ("FSISPH", 1.275) : {"Mass density" : {"L1"   : 0.06781314493410028,
+                                                   "L2"   : 0.012767602580471844,
+                                                   "Linf" : 0.6245698198195724},
+                                 "Spec Therm E" : {"L1"   : 0.00011965457154529887,
+                                                   "L2"   : 2.254867764655585e-05,
+                                                   "Linf" : 0.0010923441579020216},
+                                 "velocity    " : {"L1"   : 0.004913235403786584,
+                                                   "L2"   : 0.0009163181868064007,
+                                                   "Linf" : 0.044872645505280966},
+                                 "pressure    " : {"L1"   : 0.002210222289968727,
+                                                   "L2"   : 0.00039387994606202237,
+                                                   "Linf" : 0.018793847854203908},
+                                 "alpha       " : {"L1"   : 0.05903530352469833,
+                                                   "L2"   : 0.007960855343380124,
+                                                   "Linf" : 0.2738175996911776},
+                                 "h           " : {"L1"   : 0.0004319674641541397,
+                                                   "L2"   : 8.05536967933465e-05,
+                                                   "Linf" : 0.014201309071287523}},
+
+            ("CRKSPH", 1.275) : {"Mass density" : {"L1"   : 0.0679707220773017,
+                                                   "L2"   : 0.012783621322270034,
+                                                   "Linf" : 0.6245687018640083},
+                                 "Spec Therm E" : {"L1"   : 0.0001201303520771047,
+                                                   "L2"   : 2.2622550331357265e-05,
+                                                   "Linf" : 0.00109234444748564},
+                                 "velocity    " : {"L1"   : 0.004926121368235753,
+                                                   "L2"   : 0.0009173629106343205,
+                                                   "Linf" : 0.044872623201390446},
+                                 "pressure    " : {"L1"   : 0.0022258225868044238,
+                                                   "L2"   : 0.00039496818856079414,
+                                                   "Linf" : 0.018793890216913627},
+                                 "alpha       " : {"L1"   : 0.05909030710035451,
+                                                   "L2"   : 0.007965976598237856,
+                                                   "Linf" : 0.2738173870953471},
+                                 "h           " : {"L1"   : 0.00043274827065262975,
+                                                   "L2"   : 8.062817025125132e-05,
+                                                   "Linf" : 0.014201309070451522}},
 }
 
 #-------------------------------------------------------------------------------
@@ -579,10 +603,10 @@ def epsX_from_alphaX(alphaX,
         Linf = Pn.gridpnorm("inf", xmin, xmax)
         print(f"{name}\t\t{L1} \t\t{L2} \t\t{Linf}")
 
-        if checkError and not (np.allclose(L1, LnormRef[hydroType][name]["L1"], tol, tol) and
-                               np.allclose(L2, LnormRef[hydroType][name]["L2"], tol, tol) and
-                               np.allclose(Linf, LnormRef[hydroType][name]["Linf"], tol, tol)):
-            print("Failing Lnorm tolerance for ", name, (L1, L2, Linf), LnormRef[hydroType][name])
+        if checkError and not (np.allclose(L1, LnormRef[(hydroType, alpha0)][name]["L1"], tol, tol) and
+                               np.allclose(L2, LnormRef[(hydroType, alpha0)][name]["L2"], tol, tol) and
+                               np.allclose(Linf, LnormRef[(hydroType, alpha0)][name]["Linf"], tol, tol)):
+            print("Failing Lnorm tolerance for ", name, (L1, L2, Linf), LnormRef[(hydroType, alpha0)][name])
             failure = True
 sys.stdout.flush()
 

diff --git a/tests/functional/Porosity/PlanarCompaction/PlanarCompactionSolution.py b/tests/functional/Porosity/PlanarCompaction/PlanarCompactionSolution.py
@@ -145,7 +145,10 @@ def soundSpeed_solution(self, t,
         us, rhos, epss, Ps, alphas, ue, rhoe, epse, Pe, alphae, xs, xe, v1, h1, v2, h2 = self.waveProperties(t)
         def _soundSpeed(rhoi, epsi, alphai):
             cS0 = self.eos.soundSpeed(alphai*rhoi, epsi)
-            return cS0 + (alphai - 1.0)/(self.alpha0 - 1.0)*(self.crushCurve.c0 - cS0)
+            if self.alpha0 > 1.0:
+                return cS0 + (alphai - 1.0)/(self.alpha0 - 1.0)*(self.crushCurve.c0 - cS0)
+            else:
+                return cS0
         c_s = _soundSpeed(rhos, epss, alphas)
         c_e = _soundSpeed(rhoe, epse, alphae)
         c_0 = self.crushCurve.c0