diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index deb1e5079..55b8d4ed6 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -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. @@ -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 ============================================== diff --git a/src/CRKSPH/CRKSPHHydros.py b/src/CRKSPH/CRKSPHHydros.py index 8e398606f..b5ddd0569 100644 --- a/src/CRKSPH/CRKSPHHydros.py +++ b/src/CRKSPH/CRKSPHHydros.py @@ -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. @@ -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, diff --git a/src/FSISPH/FSISPHHydros.py b/src/FSISPH/FSISPHHydros.py index 4a6ca1db5..044dacef6 100644 --- a/src/FSISPH/FSISPHHydros.py +++ b/src/FSISPH/FSISPHHydros.py @@ -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 @@ -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 diff --git a/src/Porosity/PalphaPorosity.cc b/src/Porosity/PalphaPorosity.cc index c4a3e6d55..99f3736d4 100644 --- a/src/Porosity/PalphaPorosity.cc +++ b/src/Porosity/PalphaPorosity.cc @@ -109,9 +109,6 @@ PalphaPorosity(const SolidNodeList& 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); } } diff --git a/src/SPH/SPHHydroBase.cc b/src/SPH/SPHHydroBase.cc index 13e0c3e6f..0a9c5f6a2 100644 --- a/src/SPH/SPHHydroBase.cc +++ b/src/SPH/SPHHydroBase.cc @@ -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; @@ -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 diff --git a/src/SPH/SPHHydroBaseRZ.cc b/src/SPH/SPHHydroBaseRZ.cc index c9789551b..59dc54262 100644 --- a/src/SPH/SPHHydroBaseRZ.cc +++ b/src/SPH/SPHHydroBaseRZ.cc @@ -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; @@ -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 diff --git a/src/SPH/SPHHydros.py b/src/SPH/SPHHydros.py index 0a7a527cd..cfe8daa0c 100644 --- a/src/SPH/SPHHydros.py +++ b/src/SPH/SPHHydros.py @@ -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 @@ -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 diff --git a/src/SPH/SolidSPHHydroBase.cc b/src/SPH/SolidSPHHydroBase.cc index 72a1164cc..2295f311a 100644 --- a/src/SPH/SolidSPHHydroBase.cc +++ b/src/SPH/SolidSPHHydroBase.cc @@ -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; @@ -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 diff --git a/src/SPH/SolidSPHHydroBaseRZ.cc b/src/SPH/SolidSPHHydroBaseRZ.cc index f03222aad..3757aa5a2 100644 --- a/src/SPH/SolidSPHHydroBaseRZ.cc +++ b/src/SPH/SolidSPHHydroBaseRZ.cc @@ -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; @@ -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. diff --git a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py index b11576f72..8be099298 100644 --- a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py +++ b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py @@ -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) @@ -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) @@ -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")