Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix for porosity interaction with damage when starting with zero porosity #321

Merged
merged 4 commits into from
Jan 3, 2025
Merged

Fix for porosity interaction with damage when starting with zero porosity #321

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5884d73
Robustness changes for setting phi=0 (alpha=1) in the porosity model
jmikeowen Jan 2, 2025
c468000
Correcting porous damage evolution to handle limit of no initial poro…
jmikeowen Jan 2, 2025
7ed15d9
Adding unit tests for porosity model and damage with zero initial por…
jmikeowen Jan 2, 2025
bfb5c1d
Copy/paste name dependency errors in ATS test names for Noh-cylindrical
jmikeowen Jan 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions RELEASE_NOTES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
==============================================
Expand Down
4 changes: 2 additions & 2 deletions src/Damage/ProbabilisticDamagePolicy.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
2 changes: 1 addition & 1 deletion src/Porosity/ShadowPalphaPorosity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
23 changes: 12 additions & 11 deletions src/Porosity/computeHugoniotWithPorosity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions tests/functional/Hydro/Noh/Noh-cylindrical-2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
144 changes: 84 additions & 60 deletions tests/functional/Porosity/PlanarCompaction/PlanarCompaction-1d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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}},
}

#-------------------------------------------------------------------------------
Expand Down Expand Up @@ -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()

Expand Down
5 changes: 4 additions & 1 deletion tests/functional/Porosity/PlanarCompaction/PlanarCompactionSolution.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down