Skip to content

Commit

Permalink
Merge pull request #274 from LLNL/bugfix/PalphaConstructor
Browse files Browse the repository at this point in the history 
Bugfixes for P-alpha constructor and ASPH algorithm
  • Loading branch information
@mdavis36
mdavis36 authored Jun 25, 2024
2 parents 28fb2fc + 69c5f2d commit b03d92f
Show file tree
Hide file tree
Showing 10 changed files with 95 additions and 74 deletions.
4 changes: 4 additions & 0 deletions RELEASE_NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ Notable changes include:

* New features/ API changes:
* added MFV hydro from Hopkins 2015 with extension for ALE options
* Adding optional user specified smoothing scale method for SPH, FSISPH, and CRKSPH

* Build changes / improvements:
* PYBind11 libraries no longer depend on the structure of the PYB11 source directory.
Expand All @@ -24,6 +25,9 @@ Notable changes include:
* Versions for python dependencies in the Spheral spack recipe are fixed and updated (in some cases).

* Bug Fixes / improvements:
* Corrected an erroneous VERIFY in the P-alpha porosity constructor (with Fields of porosity and sound speed) that forced runs to stop even with correct input parameters
* Fixed a bug in the standard ASPH hydros (ASPH, SolidASPH, and RZ varieties) that gave incorrect results. FSI ad CRK models with ASPH smoothing scales were OK, but standard
SPH using ASPH smoothing scales were simply incorrect for non-unit aspect ratio H's. Also added ATS tests to help catch such errors going forward.

Version v2024.01.1 -- Release date 2024-02-17
==============================================
Expand Down
12 changes: 7 additions & 5 deletions src/CRKSPH/CRKSPHHydros.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@ def CRKSPH(dataBase,
damageRelieveRubble = False,
ASPH = False,
etaMinAxis = 0.1,
crktype = "default"):
crktype = "default",
smoothingScaleMethod = None):

# We use the provided DataBase to sniff out what sort of NodeLists are being
# used, and based on this determine which SPH object to build.
Expand Down Expand Up @@ -62,10 +63,11 @@ def CRKSPH(dataBase,
Q = eval("LimitedMonaghanGingoldViscosity%id(Clinear=%g, Cquadratic=%g)" % (ndim, Cl, Cq))

# Smoothing scale update
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)
if smoothingScaleMethod is None:
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)

# Build the constructor arguments
kwargs = {"smoothingScaleMethod" : smoothingScaleMethod,
Expand Down
63 changes: 33 additions & 30 deletions src/FSISPH/FSISPHHydros.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,32 +4,34 @@
dims = spheralDimensions()

def FSISPH(dataBase,
W,
Q = None,
slides=None,
cfl = 0.35,
surfaceForceCoefficient=0.0,
densityStabilizationCoefficient=0.1,
specificThermalEnergyDiffusionCoefficient=0.1,
xsphCoefficient=0.0,
interfaceMethod=HLLCInterface,
kernelAveragingMethod = NeverAverageKernels,
sumDensityNodeLists=[],
useVelocityMagnitudeForDt = False,
compatibleEnergyEvolution = True,
evolveTotalEnergy = False,
linearCorrectGradients = True,
planeStrain = False,
interfacePmin = 0.0,
interfaceNeighborAngleThreshold=0.707,
HUpdate = IdealH,
densityUpdate = FSISumMassDensity,
epsTensile = 0.0,
nTensile = 4.0,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
ASPH = False,
RZ = False):
W,
Q = None,
slides=None,
cfl = 0.35,
surfaceForceCoefficient=0.0,
densityStabilizationCoefficient=0.1,
specificThermalEnergyDiffusionCoefficient=0.1,
xsphCoefficient=0.0,
interfaceMethod=HLLCInterface,
kernelAveragingMethod = NeverAverageKernels,
sumDensityNodeLists=[],
useVelocityMagnitudeForDt = False,
compatibleEnergyEvolution = True,
evolveTotalEnergy = False,
linearCorrectGradients = True,
planeStrain = False,
interfacePmin = 0.0,
interfaceNeighborAngleThreshold=0.707,
HUpdate = IdealH,
densityUpdate = FSISumMassDensity,
epsTensile = 0.0,
nTensile = 4.0,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
ASPH = False,
RZ = False,
smoothingScaleMethod = None):

######################################################################
# some of these parameters are inactive and possible on there was out.
# strengthInDamage and damageRelieveRubble are old switches and are not
Expand Down Expand Up @@ -86,10 +88,11 @@ def FSISPH(dataBase,
slides = eval("SlideSurface%id(dataBase,contactTypes)" % ndim)

# Smoothing scale update
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)
if smoothingScaleMethod is None:
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)

# Build the constructor arguments
xmin = (ndim,) + xmin
Expand Down
3 changes: 0 additions & 3 deletions src/Porosity/PalphaPorosity.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,6 @@ PalphaPorosity(const SolidNodeList<Dimension>& nodeList,
for (auto i = 0u; i < n; ++i) {
mc0[i] = c0[i];
}
const auto alpha0_min = mAlpha0.min();
VERIFY2((1.0 <= mAlphae) and (mAlphae <= mAlphat) and (mAlphat <= alpha0_min),
"PalphaPorosity input ERROR : require 1.0 <= alphae <= alphat <= alpha0, (alphae, alphat, alpha0) = " << mAlphae << ", " << mAlphat << ", " << alpha0_min);
}
}

Expand Down
11 changes: 6 additions & 5 deletions src/SPH/SPHHydroBase.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -843,15 +843,16 @@ evaluateDerivatives(const typename Dimension::Scalar time,
const auto etaj = Hj*rij;
const auto etaMagi = etai.magnitude();
const auto etaMagj = etaj.magnitude();
const auto etaUnit = etai*safeInvVar(etaMagi);
const auto etaiUnit = etai*safeInvVar(etaMagi);
const auto etajUnit = etaj*safeInvVar(etaMagj);
CHECK(etaMagi >= 0.0);
CHECK(etaMagj >= 0.0);

// Symmetrized kernel weight and gradient.
W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi);
W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj);
gradWi = gWi*Hi*etaUnit;
gradWj = gWj*Hj*etaUnit;
gradWi = gWi*Hi*etaiUnit;
gradWj = gWj*Hj*etajUnit;
if (oneKernel) {
WQi = Wi;
WQj = Wj;
Expand All @@ -860,8 +861,8 @@ evaluateDerivatives(const typename Dimension::Scalar time,
} else {
WQ.kernelAndGradValue(etaMagi, Hdeti, WQi, gWQi);
WQ.kernelAndGradValue(etaMagj, Hdetj, WQj, gWQj);
gradWQi = gWQi*Hi*etaUnit;
gradWQj = gWQj*Hj*etaUnit;
gradWQi = gWQi*Hi*etaiUnit;
gradWQj = gWQj*Hj*etajUnit;
}

// Zero'th and second moment of the node distribution -- used for the
Expand Down
11 changes: 6 additions & 5 deletions src/SPH/SPHHydroBaseRZ.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,15 +392,16 @@ evaluateDerivatives(const Dim<2>::Scalar /*time*/,
const auto etaj = Hj*xij;
const auto etaMagi = etai.magnitude();
const auto etaMagj = etaj.magnitude();
const auto etaUnit = etai*safeInvVar(etaMagi);
const auto etaiUnit = etai*safeInvVar(etaMagi);
const auto etajUnit = etaj*safeInvVar(etaMagj);
CHECK(etaMagi >= 0.0);
CHECK(etaMagj >= 0.0);

// Symmetrized kernel weight and gradient.
W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi);
W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj);
gradWi = gWi*Hi*etaUnit;
gradWj = gWj*Hj*etaUnit;
gradWi = gWi*Hi*etaiUnit;
gradWj = gWj*Hj*etajUnit;
if (oneKernel) {
WQi = Wi;
WQj = Wj;
Expand All @@ -409,8 +410,8 @@ evaluateDerivatives(const Dim<2>::Scalar /*time*/,
} else {
WQ.kernelAndGradValue(etaMagi, Hdeti, WQi, gWQi);
WQ.kernelAndGradValue(etaMagj, Hdetj, WQj, gWQj);
gradWQi = gWQi*Hi*etaUnit;
gradWQj = gWQj*Hj*etaUnit;
gradWQi = gWQi*Hi*etaiUnit;
gradWQj = gWQj*Hj*etajUnit;
}

// Zero'th and second moment of the node distribution -- used for the
Expand Down
18 changes: 7 additions & 11 deletions src/SPH/SPHHydros.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@ def SPH(W,
xmin = (-1e100, -1e100, -1e100),
xmax = ( 1e100, 1e100, 1e100),
etaMinAxis = 0.1,
ASPH = False):
ASPH = False,
smoothingScaleMethod = None):

# Check if we're running solid or fluid hydro
nfluid = dataBase.numFluidNodeLists
Expand Down Expand Up @@ -86,16 +87,11 @@ def SPH(W,
Q = eval("LimitedMonaghanGingoldViscosity%id(Clinear=%g, Cquadratic=%g)" % (ndim, Cl, Cq))

# Smoothing scale update
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)

# Smoothing scale update
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)
if smoothingScaleMethod is None:
if ASPH:
smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim)
else:
smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim)

# Build the constructor arguments
xmin = (ndim,) + xmin
Expand Down
15 changes: 8 additions & 7 deletions src/SPH/SolidSPHHydroBase.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -548,15 +548,16 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/,
const auto etaj = Hj*rij;
const auto etaMagi = etai.magnitude();
const auto etaMagj = etaj.magnitude();
const auto etaUnit = etai*safeInvVar(etaMagi);
const auto etaiUnit = etai*safeInvVar(etaMagi);
const auto etajUnit = etaj*safeInvVar(etaMagj);
CHECK(etaMagi >= 0.0);
CHECK(etaMagj >= 0.0);

// Symmetrized kernel weight and gradient.
W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi);
W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj);
gradWi = gWi*Hi*etaUnit;
gradWj = gWj*Hj*etaUnit;
gradWi = gWi*Hi*etaiUnit;
gradWj = gWj*Hj*etajUnit;
if (oneKernelQ) {
WQi = Wi;
WQj = Wj;
Expand All @@ -565,15 +566,15 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/,
} else {
WQ.kernelAndGradValue(etaMagi, Hdeti, WQi, gWQi);
WQ.kernelAndGradValue(etaMagj, Hdetj, WQj, gWQj);
gradWQi = gWQi*Hi*etaUnit;
gradWQj = gWQj*Hj*etaUnit;
gradWQi = gWQi*Hi*etaiUnit;
gradWQj = gWQj*Hj*etajUnit;
}
if (oneKernelG) {
gradWGi = gradWi;
gradWGj = gradWj;
} else {
gradWGi = Hi*etaUnit * WG.gradValue(etaMagi, Hdeti);
gradWGj = Hj*etaUnit * WG.gradValue(etaMagj, Hdetj);
gradWGi = Hi*etaiUnit * WG.gradValue(etaMagi, Hdeti);
gradWGj = Hj*etajUnit * WG.gradValue(etaMagj, Hdetj);
}

// Zero'th and second moment of the node distribution -- used for the
Expand Down
15 changes: 8 additions & 7 deletions src/SPH/SolidSPHHydroBaseRZ.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,15 +469,16 @@ evaluateDerivatives(const Dim<2>::Scalar /*time*/,
const auto etaj = Hj*xij;
const auto etaMagi = etai.magnitude();
const auto etaMagj = etaj.magnitude();
const auto etaUnit = etai*safeInvVar(etaMagi);
const auto etaiUnit = etai*safeInvVar(etaMagi);
const auto etajUnit = etaj*safeInvVar(etaMagj);
CHECK(etaMagi >= 0.0);
CHECK(etaMagj >= 0.0);

// Symmetrized kernel weight and gradient.
W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi);
W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj);
gradWi = gWi*Hi*etaUnit;
gradWj = gWj*Hj*etaUnit;
gradWi = gWi*Hi*etaiUnit;
gradWj = gWj*Hj*etajUnit;
if (oneKernelQ) {
WQi = Wi;
WQj = Wj;
Expand All @@ -486,15 +487,15 @@ evaluateDerivatives(const Dim<2>::Scalar /*time*/,
} else {
WQ.kernelAndGradValue(etaMagi, Hdeti, WQi, gWQi);
WQ.kernelAndGradValue(etaMagj, Hdetj, WQj, gWQj);
gradWQi = gWQi*Hi*etaUnit;
gradWQj = gWQj*Hj*etaUnit;
gradWQi = gWQi*Hi*etaiUnit;
gradWQj = gWQj*Hj*etajUnit;
}
if (oneKernelG) {
gradWGi = gradWi;
gradWGj = gradWj;
} else {
gradWGi = Hi*etaUnit * WG.gradValue(etaMagi, Hdeti);
gradWGj = Hj*etaUnit * WG.gradValue(etaMagj, Hdetj);
gradWGi = Hi*etaiUnit * WG.gradValue(etaMagi, Hdeti);
gradWGj = Hj*etajUnit * WG.gradValue(etaMagj, Hdetj);
}

// Determine how we're applying damage.
Expand Down
17 changes: 16 additions & 1 deletion tests/functional/Hydro/Noh/Noh-cylindrical-2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,11 @@
#ATS:sph0 = test( SELF, "--crksph False --nRadial 100 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories True --steps 100", label="Noh cylindrical SPH, nPerh=2.0", np=8)
#ATS:sph1 = testif(sph0, SELF, "--crksph False --nRadial 100 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical SPH, nPerh=2.0, restart test", np=8)
#
# ASPH
#
#ATS:asph0 = test( SELF, "--crksph False --asph True --nRadial 100 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories True --steps 100", label="Noh cylindrical ASPH, nPerh=2.0", np=8)
#ATS:asph1 = testif(sph0, SELF, "--crksph False --asph True --nRadial 100 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical ASPH, nPerh=2.0, restart test", np=8)
#
# CRK (SumVolume)
#
#ATS:crk0 = test( SELF, "--crksph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKSumVolume --clearDirectories True --steps 50", label="Noh cylindrical CRK (sum vol), nPerh=2.0", np=2)
Expand All @@ -14,6 +19,16 @@
#ATS:crk2 = test( SELF, "--crksph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKVoronoiVolume --clearDirectories True --steps 50", label="Noh cylindrical CRK (Voronoi vol), nPerh=2.0", np=2)
#ATS:crk3 = testif(crk2, SELF, "--crksph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKVoronoiVolume --clearDirectories False --steps 10 --restoreCycle 40 --checkRestart True", label="Noh cylindrical CRK (Voronoi vol) , nPerh=2.0, restart test", np=2)
#
# ACRK (SumVolume)
#
#ATS:acrk0 = test( SELF, "--crksph True --asph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKSumVolume --clearDirectories True --steps 50", label="Noh cylindrical ACRK (sum vol), nPerh=2.0", np=2)
#ATS:acrk1 = testif(acrk0, SELF, "--crksph True --asph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKSumVolume --clearDirectories False --steps 10 --restoreCycle 40 --checkRestart True", label="Noh cylindrical ACRK (sum vol), nPerh=2.0, restart test", np=2)
#
# ACRK (VoroniVolume)
#
#ATS:acrk2 = test( SELF, "--crksph True --asph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKVoronoiVolume --clearDirectories True --steps 50", label="Noh cylindrical ACRK (Voronoi vol), nPerh=2.0", np=2)
#ATS:acrk3 = testif(acrk2, SELF, "--crksph True --asph True --nRadial 20 --cfl 0.25 --Cl 1.0 --Cq 1.0 --filter 0.0 --nPerh 2.01 --graphics False --restartStep 20 --volumeType RKVoronoiVolume --clearDirectories False --steps 10 --restoreCycle 40 --checkRestart True", label="Noh cylindrical ACRK (Voronoi vol) , nPerh=2.0, restart test", np=2)
#
# GSPH
#
#ATS:gsph0 = test( SELF, "--gsph True --nRadial 100 --cfl 0.25 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories True --steps 100", label="Noh cylindrical GSPH, nPerh=2.0", np=8)
Expand Down Expand Up @@ -551,7 +566,7 @@
vizTime = vizTime,
vizDerivs = vizDerivs,
#skipInitialPeriodicWork = SVPH,
SPH = True, # Only for iterating H
SPH = not asph, # Only for iterating H
)
output("control")

Expand Down

0 comments on commit b03d92f

Please sign in to comment.