From 5884d73b93dafffe8c2679e64d5c6eeee31cab54 Mon Sep 17 00:00:00 2001 From: Mike Owen Date: Thu, 2 Jan 2025 11:04:18 -0800 Subject: [PATCH 1/4] Robustness changes for setting phi=0 (alpha=1) in the porosity model --- src/Porosity/ShadowPalphaPorosity.py | 2 +- src/Porosity/computeHugoniotWithPorosity.py | 23 ++++++++++--------- .../PlanarCompactionSolution.py | 5 +++- 3 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/Porosity/ShadowPalphaPorosity.py b/src/Porosity/ShadowPalphaPorosity.py index d298d6f1e..93f6e06d0 100644 --- 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 index 5531b286d..00af13210 100644 --- 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/Porosity/PlanarCompaction/PlanarCompactionSolution.py b/tests/functional/Porosity/PlanarCompaction/PlanarCompactionSolution.py index 31bd677f0..2f769149f 100644 --- 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 From c46800052261f0189031cd7bc825cd070a92331e Mon Sep 17 00:00:00 2001 From: Mike Owen Date: Thu, 2 Jan 2025 11:19:29 -0800 Subject: [PATCH 2/4] Correcting porous damage evolution to handle limit of no initial porosity --- src/Damage/ProbabilisticDamagePolicy.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Damage/ProbabilisticDamagePolicy.cc b/src/Damage/ProbabilisticDamagePolicy.cc index 3f7c7a685..d6d74e8b3 100644 --- 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); } From 7ed15d9486473b0864fbf198948d5bd443e039e8 Mon Sep 17 00:00:00 2001 From: Mike Owen Date: Thu, 2 Jan 2025 13:52:55 -0800 Subject: [PATCH 3/4] Adding unit tests for porosity model and damage with zero initial porosity --- RELEASE_NOTES.md | 1 + .../PlanarCompaction/PlanarCompaction-1d.py | 144 ++++++++++-------- 2 files changed, 85 insertions(+), 60 deletions(-) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 5bd2124cc..a4e41d72a 100644 --- 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/tests/functional/Porosity/PlanarCompaction/PlanarCompaction-1d.py b/tests/functional/Porosity/PlanarCompaction/PlanarCompaction-1d.py index 6f2b2e314..086f8595d 100644 --- 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() From bfb5c1d9186e7b9588df33650c304b9213507fa1 Mon Sep 17 00:00:00 2001 From: Mike Owen Date: Thu, 2 Jan 2025 15:15:53 -0800 Subject: [PATCH 4/4] Copy/paste name dependency errors in ATS test names for Noh-cylindrical --- tests/functional/Hydro/Noh/Noh-cylindrical-2d.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py index 9b6309fc4..d38699676 100644 --- 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.