From ce886441381559d32c57c932ce68f2586057b9a7 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 16 Oct 2021 10:49:27 -0700 Subject: [PATCH 1/6] removing outdated debug checks c=<0 is handled by the riemann solvers now --- src/GSPH/GSPHHydroBase.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/GSPH/GSPHHydroBase.cc b/src/GSPH/GSPHHydroBase.cc index 44d0e58e7..f425c7557 100644 --- a/src/GSPH/GSPHHydroBase.cc +++ b/src/GSPH/GSPHHydroBase.cc @@ -831,10 +831,6 @@ evaluateDerivatives(const typename Dimension::Scalar time, const auto Peffj = Pj + Rj; // Acceleration. - CHECK(rhoi > 0.0); - CHECK(rhoj > 0.0); - CHECK(ci > 0.0); - CHECK(cj > 0.0); const auto voli = mi/rhoi; const auto volj = mj/rhoj; From cd1276c6108c2040f2183875a208b72ac791d8e8 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 16 Oct 2021 12:30:37 -0700 Subject: [PATCH 2/6] adding GSPH option to gas-gas sod tests --- tests/Hydro/Sod/Sod-planar-1d.py | 101 +++++++++++++++++----------- tests/Hydro/Sod/Sod-planar-2d.py | 109 ++++++++++++++++++------------- tests/Hydro/Sod/Sod-planar-3d.py | 95 ++++++++++++++++----------- 3 files changed, 187 insertions(+), 118 deletions(-) diff --git a/tests/Hydro/Sod/Sod-planar-1d.py b/tests/Hydro/Sod/Sod-planar-1d.py index 39a885e11..8a089974f 100644 --- a/tests/Hydro/Sod/Sod-planar-1d.py +++ b/tests/Hydro/Sod/Sod-planar-1d.py @@ -14,6 +14,11 @@ #ATS:fsisph1 = test( SELF, "--crksph False --fsisph True --solid True --cfl 0.25 --graphics None --clearDirectories True --restartStep 20 --steps 40", label="Planar Sod problem with FSISPH -- 1-D (serial)") #ATS:fsisph2 = testif(fsisph1, SELF, "--crksph False --fsisph True --solid True --cfl 0.25 --graphics None --clearDirectories False --restartStep 20 --steps 20 --restoreCycle 20 --checkRestart True", label="Planar Sod problem with FSISPH -- 1-D (serial) RESTART CHECK") # +# GSPH +# +#ATS:gsph1 = test( SELF, "--gsph True --cfl 0.25 --graphics None --clearDirectories True --restartStep 20 --steps 40", label="Planar Sod problem with GSPH -- 1-D (serial)") +#ATS:gsph2 = testif(gsph1, SELF, "--gsph True --cfl 0.25 --graphics None --clearDirectories False --restartStep 20 --steps 20 --restoreCycle 20 --checkRestart True", label="Planar Sod problem with GSPH -- 1-D (serial) RESTART CHECK") +# import os, sys import shutil from SolidSpheral1d import * @@ -51,7 +56,8 @@ crksph = False, psph = False, fsisph = False, - + gsph = False, + evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy solid = False, # If true, use the fluid limit of the solid hydro option boolReduceViscosity = False, @@ -127,6 +133,7 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not (GSPH and (boolReduceViscosity or boolCullenViscosity)) assert not svph assert not (fsisph and not solid) assert numNodeLists in (1, 2) @@ -141,6 +148,8 @@ hydroname = "PSPH" elif fsisph: hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if solid: @@ -348,6 +357,22 @@ def specificEnergy(xi, rhoi, gammai): evolveTotalEnergy = evolveTotalEnergy, correctVelocityGradient = correctVelocityGradient, HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True,RiemannGradient) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) else: hydro = SPH(dataBase = db, W = WT, @@ -368,40 +393,41 @@ def specificEnergy(xi, rhoi, gammai): #------------------------------------------------------------------------------- # Tweak the artificial viscosity. #------------------------------------------------------------------------------- -q = hydro.Q -if not Cl is None: - q.Cl = Cl -if not Cq is None: - q.Cq = Cq -if not linearInExpansion is None: - q.linearInExpansion = linearInExpansion -if not quadraticInExpansion is None: - q.quadraticInExpansion = quadraticInExpansion -if not Qlimiter is None: - q.limiter = Qlimiter -if not epsilon2 is None: - q.epsilon2 = epsilon2 -if not etaCritFrac is None: - q.etaCritFrac = etaCritFrac -if not etaFoldFrac is None: - q.etaFoldFrac = etaFoldFrac -output("q") -output("q.Cl") -output("q.Cq") -output("q.limiter") -output("q.epsilon2") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + if not Cl is None: + q.Cl = Cl + if not Cq is None: + q.Cq = Cq + if not linearInExpansion is None: + q.linearInExpansion = linearInExpansion + if not quadraticInExpansion is None: + q.quadraticInExpansion = quadraticInExpansion + if not Qlimiter is None: + q.limiter = Qlimiter + if not epsilon2 is None: + q.epsilon2 = epsilon2 + if not etaCritFrac is None: + q.etaCritFrac = etaCritFrac + if not etaFoldFrac is None: + q.etaFoldFrac = etaFoldFrac + output("q") + output("q.Cl") + output("q.Cq") + output("q.limiter") + output("q.epsilon2") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Construct the Artificial Conduction physics object. @@ -635,10 +661,11 @@ def createList(x): plots += [(volPlot, "Sod-planar-vol.png"), (splot, "Sod-planar-surfacePoint.png")] - viscPlot = plotFieldList(hydro.maxViscousPressure, + if not gsph: + viscPlot = plotFieldList(hydro.maxViscousPressure, winTitle = "max($\\rho^2 \pi_{ij}$)", colorNodeLists = False) - plots.append((viscPlot, "Sod-planar-viscosity.png")) + plots.append((viscPlot, "Sod-planar-viscosity.png")) if boolCullenViscosity: cullAlphaPlot = plotFieldList(q.ClMultiplier, diff --git a/tests/Hydro/Sod/Sod-planar-2d.py b/tests/Hydro/Sod/Sod-planar-2d.py index e42447e78..5067fdef2 100644 --- a/tests/Hydro/Sod/Sod-planar-2d.py +++ b/tests/Hydro/Sod/Sod-planar-2d.py @@ -45,6 +45,7 @@ crksph = False, psph = False, fsisph = False, + gsph= False, evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy solid = False, # If true, use the fluid limit of the solid hydro option boolReduceViscosity = False, @@ -116,6 +117,7 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not (GSPH and (boolReduceViscosity or boolCullenViscosity)) assert numNodeLists in (1, 2) assert not (fsisph and not solid) assert not svph @@ -128,6 +130,8 @@ hydroname = "PSPH" elif fsisph: hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if solid: @@ -381,6 +385,22 @@ def specificEnergy(xi, rhoi, gammai): evolveTotalEnergy = evolveTotalEnergy, correctVelocityGradient = correctVelocityGradient, HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True,RiemannGradient) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) else: hydro = SPH(dataBase = db, W = WT, @@ -401,47 +421,48 @@ def specificEnergy(xi, rhoi, gammai): #------------------------------------------------------------------------------- # Tweak the artificial viscosity. #------------------------------------------------------------------------------- -q = hydro.Q -if not Cl is None: - q.Cl = Cl -if not Cq is None: - q.Cq = Cq -if not linearInExpansion is None: - q.linearInExpansion = linearInExpansion -if not quadraticInExpansion is None: - q.quadraticInExpansion = quadraticInExpansion -if not Qlimiter is None: - q.limiter = Qlimiter -if not epsilon2 is None: - q.epsilon2 = epsilon2 -if not etaCritFrac is None: - try: - q.etaCritFrac = etaCritFrac - except: - pass -if not etaFoldFrac is None: - try: - q.etaFoldFrac = etaFoldFrac - except: - pass -output("q") -output("q.Cl") -output("q.Cq") -output("q.limiter") -output("q.epsilon2") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + if not Cl is None: + q.Cl = Cl + if not Cq is None: + q.Cq = Cq + if not linearInExpansion is None: + q.linearInExpansion = linearInExpansion + if not quadraticInExpansion is None: + q.quadraticInExpansion = quadraticInExpansion + if not Qlimiter is None: + q.limiter = Qlimiter + if not epsilon2 is None: + q.epsilon2 = epsilon2 + if not etaCritFrac is None: + try: + q.etaCritFrac = etaCritFrac + except: + pass + if not etaFoldFrac is None: + try: + q.etaFoldFrac = etaFoldFrac + except: + pass + output("q") + output("q.Cl") + output("q.Cq") + output("q.limiter") + output("q.epsilon2") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Construct the Artificial Conduction physics object. @@ -639,13 +660,13 @@ def createList(x): (HPlot, "Sod-planar-h.png"), (csPlot, "Sod-planar-cs.png"), (Aplot, "Sod-planar-entropy.png")] - - viscPlot = plotFieldList(hydro.maxViscousPressure, + if not gsph: + viscPlot = plotFieldList(hydro.maxViscousPressure, winTitle = "max(rho^2 Piij)", plotStyle = "points", colorNodeLists = False, filterFunc = plotFilter) - plots.append((viscPlot, "Sod-planar-viscosity.png")) + plots.append((viscPlot, "Sod-planar-viscosity.png")) if crksph: volPlot = plotFieldList(control.RKCorrections.volume, winTitle = "volume", diff --git a/tests/Hydro/Sod/Sod-planar-3d.py b/tests/Hydro/Sod/Sod-planar-3d.py index 8fa13f1be..bc8c6329b 100644 --- a/tests/Hydro/Sod/Sod-planar-3d.py +++ b/tests/Hydro/Sod/Sod-planar-3d.py @@ -62,6 +62,7 @@ crksph = False, fsisph = False, psph = False, + gsph = False, asph = False, # Choose the H evolution -- works with all hydro options evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy solid = False, # If true, use the fluid limit of the solid hydro option @@ -134,6 +135,7 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not(GSPH and (boolCullenViscosity or boolReduceViscosity)) assert numNodeLists in (1, 2) assert not (fsisph and not solid) assert not svph @@ -149,6 +151,8 @@ hydroname = "PSPH" elif fsisph: hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if asph: @@ -363,6 +367,22 @@ def specificEnergy(xi, rhoi, gammai): evolveTotalEnergy = evolveTotalEnergy, correctVelocityGradient = correctVelocityGradient, HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True,RiemannGradient) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) else: hydro = SPH(dataBase = db, W = WT, @@ -383,40 +403,41 @@ def specificEnergy(xi, rhoi, gammai): #------------------------------------------------------------------------------- # Tweak the artificial viscosity. #------------------------------------------------------------------------------- -q = hydro.Q -if not Cl is None: - q.Cl = Cl -if not Cq is None: - q.Cq = Cq -if not linearInExpansion is None: - q.linearInExpansion = linearInExpansion -if not quadraticInExpansion is None: - q.quadraticInExpansion = quadraticInExpansion -if not Qlimiter is None: - q.limiter = Qlimiter -if not epsilon2 is None: - q.epsilon2 = epsilon2 -if not etaCritFrac is None: - q.etaCritFrac = etaCritFrac -if not etaFoldFrac is None: - q.etaFoldFrac = etaFoldFrac -output("q") -output("q.Cl") -output("q.Cq") -output("q.limiter") -output("q.epsilon2") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + if not Cl is None: + q.Cl = Cl + if not Cq is None: + q.Cq = Cq + if not linearInExpansion is None: + q.linearInExpansion = linearInExpansion + if not quadraticInExpansion is None: + q.quadraticInExpansion = quadraticInExpansion + if not Qlimiter is None: + q.limiter = Qlimiter + if not epsilon2 is None: + q.epsilon2 = epsilon2 + if not etaCritFrac is None: + q.etaCritFrac = etaCritFrac + if not etaFoldFrac is None: + q.etaFoldFrac = etaFoldFrac + output("q") + output("q.Cl") + output("q.Cq") + output("q.limiter") + output("q.epsilon2") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Construct the Artificial Conduction physics object. @@ -624,11 +645,11 @@ def createList(x): colorNodeLists = False) plots += [(volPlot, "Sod-planar-vol.png"), (splot, "Sod-planar-surfacePoint.png")] - - viscPlot = plotFieldList(hydro.maxViscousPressure, + if not gsph: + viscPlot = plotFieldList(hydro.maxViscousPressure, winTitle = "max($\\rho^2 \pi_{ij}$)", colorNodeLists = False) - plots.append((viscPlot, "Sod-planar-viscosity.png")) + plots.append((viscPlot, "Sod-planar-viscosity.png")) if boolCullenViscosity: cullAlphaPlot = plotFieldList(q.ClMultiplier, From 4ef426e1602d56855663aa26cfbb66f96d7d0fda Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 16 Oct 2021 14:39:20 -0700 Subject: [PATCH 3/6] gsph added to Noh --- tests/Hydro/Noh/Noh-cylindrical-2d.py | 16 +++- tests/Hydro/Noh/Noh-planar-1d.py | 2 + tests/Hydro/Noh/Noh-spherical-3d.py | 114 ++++++++++++++++++-------- 3 files changed, 95 insertions(+), 37 deletions(-) diff --git a/tests/Hydro/Noh/Noh-cylindrical-2d.py b/tests/Hydro/Noh/Noh-cylindrical-2d.py index db563f196..e0fb2afd3 100644 --- a/tests/Hydro/Noh/Noh-cylindrical-2d.py +++ b/tests/Hydro/Noh/Noh-cylindrical-2d.py @@ -137,6 +137,8 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not(fsisph and not solid) assert thetaFactor in (0.5, 1.0, 2.0) theta = thetaFactor * pi @@ -327,7 +329,19 @@ HUpdate = HUpdate, XSPH = XSPH, ASPH = asph) - +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) elif gsph: limiter = VanLeerLimiter() waveSpeed = DavisWaveSpeed() diff --git a/tests/Hydro/Noh/Noh-planar-1d.py b/tests/Hydro/Noh/Noh-planar-1d.py index 18eaec386..f81a14995 100644 --- a/tests/Hydro/Noh/Noh-planar-1d.py +++ b/tests/Hydro/Noh/Noh-planar-1d.py @@ -195,6 +195,8 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not(fsisph and not solid) if smallPressure: P0 = 1.0e-6 eps1 = P0/((gamma - 1.0)*rho1) diff --git a/tests/Hydro/Noh/Noh-spherical-3d.py b/tests/Hydro/Noh/Noh-spherical-3d.py index a99c818de..580cba239 100644 --- a/tests/Hydro/Noh/Noh-spherical-3d.py +++ b/tests/Hydro/Noh/Noh-spherical-3d.py @@ -47,6 +47,9 @@ svph = False, crksph = False, psph = False, + fsisph = False, + gsph = False, + asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. boolReduceViscosity = False, HopkinsConductivity = False, # For PSPH @@ -120,6 +123,9 @@ ) assert not(boolReduceViscosity and boolCullenViscosity) +assert not(gsph and (boolReduceViscosity and boolCullenViscosity)) +assert not(fsisph and not solid) + if smallPressure: P0 = 1.0e-6 eps0 = P0/((gamma - 1.0)*rho0) @@ -130,6 +136,10 @@ hydroname = "CRKSPH" elif psph: hydroname = "PSPH" +elif fsisph: + hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if asph: @@ -216,8 +226,8 @@ SPH = True) if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes3d - #from PeanoHilbertDistributeNodes import distributeNodes3d + #from VoronoiDistributeNodes import distributeNodes3d + from PeanoHilbertDistributeNodes import distributeNodes3d else: from DistributeNodes import distributeNodes3d @@ -272,6 +282,35 @@ densityUpdate = densityUpdate, HUpdate = HUpdate, ASPH = asph) +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True,RiemannGradient) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) elif psph: hydro = PSPH(dataBase = db, W = WT, @@ -302,49 +341,52 @@ ASPH = asph) output("hydro") output("hydro.kernel") -output("hydro.PiKernel") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") -output("hydro.densityUpdate") output("hydro.HEvolution") +if not (gsph or fsisph): + output("hydro.PiKernel") +if not fsisph: + output("hydro.densityUpdate") packages = [hydro] #------------------------------------------------------------------------------- # Set the artificial viscosity parameters. #------------------------------------------------------------------------------- -q = hydro.Q -if Cl: - q.Cl = Cl -if Cq: - q.Cq = Cq -if epsilon2: - q.epsilon2 = epsilon2 -if Qlimiter: - q.limiter = Qlimiter -if balsaraCorrection: - q.balsaraShearCorrection = balsaraCorrection -output("q") -output("q.Cl") -output("q.Cq") -output("q.epsilon2") -output("q.limiter") -output("q.balsaraShearCorrection") -try: - output("q.linearInExpansion") - output("q.quadraticInExpansion") -except: - pass - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nhQ,nhL,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + if Cl: + q.Cl = Cl + if Cq: + q.Cq = Cq + if epsilon2: + q.epsilon2 = epsilon2 + if Qlimiter: + q.limiter = Qlimiter + if balsaraCorrection: + q.balsaraShearCorrection = balsaraCorrection + output("q") + output("q.Cl") + output("q.Cq") + output("q.epsilon2") + output("q.limiter") + output("q.balsaraShearCorrection") + try: + output("q.linearInExpansion") + output("q.quadraticInExpansion") + except: + pass + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nhQ,nhL,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Create boundary conditions. From c1e48936685df98e462399be0a1f5ee892e71417 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 18 Oct 2021 09:56:58 -0700 Subject: [PATCH 4/6] updated Sedov scripts added GSPH and FSISPH hydros to sedov. Updated CRKSPH and PSPH for some input and name changes --- tests/Hydro/Sedov/Sedov-2d-ratio.py | 169 ++++++++++++++-------- tests/Hydro/Sedov/Sedov-cylindrical-2d.py | 149 +++++++++++++------ tests/Hydro/Sedov/Sedov-planar-1d.py | 168 +++++++++++++-------- tests/Hydro/Sedov/Sedov-spherical-3d.py | 140 ++++++++++++------ 4 files changed, 415 insertions(+), 211 deletions(-) diff --git a/tests/Hydro/Sedov/Sedov-2d-ratio.py b/tests/Hydro/Sedov/Sedov-2d-ratio.py index 8eaa14037..73bc13374 100644 --- a/tests/Hydro/Sedov/Sedov-2d-ratio.py +++ b/tests/Hydro/Sedov/Sedov-2d-ratio.py @@ -22,8 +22,6 @@ nTheta = 200, rmin = 0.0, rmax = 1.0, - nPerh = 1.51, - order = 5, rho0 = 1.0, eps0 = 0.0, @@ -34,15 +32,40 @@ smoothSpikeScale = 0.5, gamma = 5.0/3.0, mu = 1.0, + rhomin = 1e-10, + # kernel + HUpdate = IdealH, + nPerh = 1.51, + order = 5, + hmin = 1e-15, + hmax = 1.0, + + # hydros crksph = False, psph = False, - asph = False, # Only for H evolution, not hydro algorithm - correctionOrder = LinearOrder, - volumeType = RKSumVolume, + fsisph = False, + gsph = False, + + # hydro options + solid = False, + asph = False, + XSPH = False, + evolveTotalEnergy = False, + compatibleEnergy = True, + gradhCorrection = True, + correctVelocityGradient = True, densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - HUpdate = IdealH, filter = 0.0, + + # crksph parameters + correctionOrder = LinearOrder, + + # gsph parameters + RiemannGradientType = RiemannGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) + linearReconstruction = True, + + # Artficial Viscosity boolReduceViscosity = False, nh = 5.0, aMin = 0.1, @@ -56,17 +79,11 @@ betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - HopkinsConductivity = False, # For PSPH - evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy - - hmin = 1e-15, - hmax = 1.0, - cfl = 0.5, - useVelocityMagnitudeForDt = True, - XSPH = False, - rhomin = 1e-10, - + + # Integration IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.5, + useVelocityMagnitudeForDt = False, steps = None, goalTime = None, goalRadius = 0.8, @@ -74,23 +91,19 @@ dtMin = 1.0e-8, dtMax = None, dtGrowth = 2.0, - vizCycle = None, - vizTime = 0.1, maxSteps = None, statsStep = 1, smoothIters = 0, - HEvolution = IdealH, - compatibleEnergy = True, - gradhCorrection = True, - correctVelocityGradient = True, - + + # IO + vizCycle = None, + vizTime = 0.1, restoreCycle = -1, restartStep = 1000, - graphics = False, useVoronoiOutput = False, clearDirectories = False, - dataRoot = "dumps-cylindrical-Sedov", + dataDirBase = "dumps-cylindrical-Sedov", outputFile = "None", serialDump=True, @@ -104,7 +117,10 @@ print "WARNING: smallPressure specified, so setting eps0=%g" % eps0 assert not(boolReduceViscosity and boolCullenViscosity) +assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not(fsisph and not solid) assert thetaFactor in (0.5, 1.0, 2.0) + theta = thetaFactor * pi # Figure out what our goal time should be. @@ -136,12 +152,16 @@ hydroname = "CRKSPH" elif psph: hydroname = "PSPH" +elif fsisph: + hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if asph: hydroname = "A" + hydroname -dataDir = os.path.join(dataRoot, +dataDir = os.path.join(dataDirBase, hydroname, "nr=%i_nt=%i" % (nRadial, nTheta)) restartDir = os.path.join(dataDir, "restarts") @@ -177,7 +197,11 @@ #------------------------------------------------------------------------------- # Create a NodeList and associated Neighbor object. #------------------------------------------------------------------------------- -nodes1 = makeFluidNodeList("nodes1", eos, +if solid: + nodeListConstructor = makeSolidNodeList +else: + nodeListConstructor = makeFluidNodeList +nodes1 = nodeListConstructor("nodes1", eos, hmin = hmin, hmax = hmax, kernelExtent = kernelExtent, @@ -215,8 +239,8 @@ xlmax = xlmax) if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes2d - #from PeanoHilbertDistributeNodes import distributeNodes2d + #from VoronoiDistributeNodes import distributeNodes2d + from PeanoHilbertDistributeNodes import distributeNodes2d else: from DistributeNodes import distributeNodes2d @@ -225,10 +249,6 @@ output("mpi.reduce(nodes1.numInternalNodes, mpi.MAX)") output("mpi.reduce(nodes1.numInternalNodes, mpi.SUM)") - - - - # Set the point source of energy. Esum = 0.0 if smoothSpike or topHatSpike: @@ -285,24 +305,51 @@ #------------------------------------------------------------------------------- if crksph: hydro = CRKSPH(dataBase = db, - W = WT, + order = correctionOrder, filter = filter, cfl = cfl, compatibleEnergyEvolution = compatibleEnergy, XSPH = XSPH, - correctionOrder = correctionOrder, densityUpdate = densityUpdate, HUpdate = HUpdate, ASPH = asph) +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction,RiemannGradientType) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) elif psph: - hydro = PASPH(dataBase = db, + hydro = PSPH(dataBase = db, W = WT, filter = filter, cfl = cfl, useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, compatibleEnergyEvolution = compatibleEnergy, evolveTotalEnergy = evolveTotalEnergy, - HopkinsConductivity = HopkinsConductivity, correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, HUpdate = HUpdate, @@ -318,38 +365,36 @@ correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, XSPH = XSPH, - HUpdate = HEvolution, + HUpdate = HUpdate, ASPH = asph) + +packages = [hydro] + output("hydro") -output("hydro.kernel()") -output("hydro.PiKernel()") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.XSPH") -output("hydro.densityUpdate") output("hydro.HEvolution") -q = hydro.Q -output("q") -output("q.Cl") -output("q.Cq") -output("q.epsilon2") -output("q.limiter") -output("q.balsaraShearCorrection") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - -packages = [hydro] - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + output("q") + output("q.Cl") + output("q.Cq") + output("q.epsilon2") + output("q.limiter") + output("q.balsaraShearCorrection") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Create boundary conditions. diff --git a/tests/Hydro/Sedov/Sedov-cylindrical-2d.py b/tests/Hydro/Sedov/Sedov-cylindrical-2d.py index 19bb0bd4e..b588a003f 100644 --- a/tests/Hydro/Sedov/Sedov-cylindrical-2d.py +++ b/tests/Hydro/Sedov/Sedov-cylindrical-2d.py @@ -34,15 +34,39 @@ smoothSpikeScale = 0.5, gamma = 5.0/3.0, mu = 1.0, + rhomin = 1e-10, + + # kernel + HUpdate = IdealH, + hmin = 1e-15, + hmax = 1.0, + # hydros crksph = False, psph = False, - asph = False, # Only for H evolution, not hydro algorithm + fsisph = False, + gsph = False, + + # hydro options + solid = False, + asph = False, + XSPH = False, + evolveTotalEnergy = False, + compatibleEnergy = True, + gradhCorrection = True, + correctVelocityGradient = True, + densityUpdate = RigorousSumDensity, + filter = 0.0, + + # crksph options correctionOrder = LinearOrder, volumeType = RKSumVolume, - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - HUpdate = IdealH, - filter = 0.0, + + # gsph options + RiemannGradientType = RiemannGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) + linearReconstruction = True, + + # Artifical Viscosity boolReduceViscosity = False, nh = 5.0, aMin = 0.1, @@ -56,17 +80,10 @@ betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - HopkinsConductivity = False, # For PSPH - evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy - - hmin = 1e-15, - hmax = 1.0, - cfl = 0.25, - useVelocityMagnitudeForDt = True, - XSPH = False, - rhomin = 1e-10, - + + # Integration IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.25, steps = None, goalTime = None, goalRadius = 0.8, @@ -74,23 +91,21 @@ dtMin = 1.0e-8, dtMax = None, dtGrowth = 2.0, - vizCycle = None, - vizTime = 0.1, maxSteps = None, statsStep = 1, smoothIters = 0, - HEvolution = IdealH, - compatibleEnergy = True, - gradhCorrection = True, - correctVelocityGradient = True, + useVelocityMagnitudeForDt = False, + # IO + vizCycle = None, + vizTime = 0.1, restoreCycle = -1, restartStep = 1000, graphics = False, useVoronoiOutput = False, clearDirectories = False, - dataRoot = "dumps-cylindrical-Sedov", + dataDirBase = "dumps-cylindrical-Sedov", outputFile = "None", serialDump=True, ) @@ -102,6 +117,8 @@ assert not(boolReduceViscosity and boolCullenViscosity) assert thetaFactor in (0.5, 1.0, 2.0) +assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not(fsisph and not solid) theta = thetaFactor * pi # Figure out what our goal time should be. @@ -133,12 +150,16 @@ hydroname = "CRKSPH" elif psph: hydroname = "PSPH" +elif fsisph: + hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if asph: hydroname = "A" + hydroname -dataDir = os.path.join(dataRoot, +dataDir = os.path.join(dataDirBase, hydroname, "nperh=%4.2f" % nPerh, "XSPH=%s" % XSPH, @@ -180,7 +201,11 @@ #------------------------------------------------------------------------------- # Create a NodeList and associated Neighbor object. #------------------------------------------------------------------------------- -nodes1 = makeFluidNodeList("nodes1", eos, +if solid: + nodeListConstructor = makeSolidNodeList +else: + nodeListConstructor = makeFluidNodeList +nodes1 = nodeListConstructor("nodes1", eos, hmin = hmin, hmax = hmax, kernelExtent = kernelExtent, @@ -299,15 +324,43 @@ densityUpdate = densityUpdate, HUpdate = HUpdate, ASPH = asph) +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction,RiemannGradientType) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) elif psph: - hydro = PASPH(dataBase = db, + hydro = PSPH(dataBase = db, W = WT, filter = filter, cfl = cfl, useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, compatibleEnergyEvolution = compatibleEnergy, evolveTotalEnergy = evolveTotalEnergy, - HopkinsConductivity = HopkinsConductivity, correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, HUpdate = HUpdate, @@ -323,8 +376,11 @@ correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, XSPH = XSPH, - HUpdate = HEvolution, + HUpdate = HUpdate, ASPH = asph) + +packages = [hydro] + output("hydro") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") @@ -332,27 +388,26 @@ output("hydro.densityUpdate") output("hydro.HEvolution") -q = hydro.Q -output("q") -output("q.Cl") -output("q.Cq") -output("q.epsilon2") -output("q.limiter") -output("q.balsaraShearCorrection") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - -packages = [hydro] - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + output("q") + output("q.Cl") + output("q.Cq") + output("q.epsilon2") + output("q.limiter") + output("q.balsaraShearCorrection") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Create boundary conditions. diff --git a/tests/Hydro/Sedov/Sedov-planar-1d.py b/tests/Hydro/Sedov/Sedov-planar-1d.py index 2605ca414..4e4313854 100644 --- a/tests/Hydro/Sedov/Sedov-planar-1d.py +++ b/tests/Hydro/Sedov/Sedov-planar-1d.py @@ -12,31 +12,55 @@ #------------------------------------------------------------------------------- # Generic problem parameters #------------------------------------------------------------------------------- -commandLine(nRadial = 50, +commandLine(# discretization & domain + nRadial = 50, rmin = 0.0, rmax = 1.0, - nPerh = 1.51, - order = 5, + # material properties rho0 = 1.0, eps0 = 0.0, Espike = 1.0, - smoothSpike = True, + smoothSpike = False, topHatSpike = False, smoothSpikeScale = 0.5, gamma = 5.0/3.0, mu = 1.0, smallPressure = False, + rhomin = 1e-10, + #Kernel + HUpdate = IdealH, + nPerh = 1.51, + order = 5, + hmin = 1e-15, + hmax = 1.0, + + # hydro type crksph = False, psph = False, - asph = False, # Selects the H update algorithm -- can be used with CRK, PSPH, SPH, etc. - correctionOrder = LinearOrder, - volumeType = RKSumVolume, + gsph = False, + fsisph = False, + + # hydro options + solid = False, + asph = False, # Selects the H update algorithm -- can be used with CRK, PSPH, SPH, etc. + XSPH = False, + evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy + compatibleEnergy = True, + gradhCorrection = True, + correctVelocityGradient = True, densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - HUpdate = IdealH, filter = 0.0, + + # crksph options + correctionOrder = LinearOrder, + # gsph options + RiemannGradientType = RiemannGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) + linearReconstruction = True, + + # artificial viscosity options boolReduceViscosity = False, nh = 5.0, aMin = 0.1, @@ -49,16 +73,10 @@ betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - HopkinsConductivity = False, # For PSPH - evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy - - hmin = 1e-15, - hmax = 1.0, + + # integration cfl = 0.5, - useVelocityMagnitudeForDt = True, - XSPH = False, - rhomin = 1e-10, - + useVelocityMagnitudeForDt = False, steps = None, goalTime = None, goalRadius = 0.8, @@ -70,17 +88,14 @@ maxSteps = None, statsStep = 1, smoothIters = 0, - HEvolution = IdealH, - compatibleEnergy = True, - gradhCorrection = True, - correctVelocityGradient = True, + # IO restoreCycle = -1, restartStep = 1000, graphics = True, clearDirectories = False, - dataRoot = "dumps-planar-Sedov", + dataDirBase = "dumps-planar-Sedov", outputFile = "None", ) @@ -90,6 +105,9 @@ print "WARNING: smallPressure specified, so setting eps0=%g" % eps0 assert not(boolReduceViscosity and boolCullenViscosity) +assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not(fsisph and not solid) + # Figure out what our goal time should be. import SedovAnalyticSolution h0 = 1.0/nRadial*nPerh @@ -116,12 +134,16 @@ hydroname = "CRKSPH" elif psph: hydroname = "PSPH" +elif gsph: + hydroname = "GSPH" +elif fsisph: + hydroname = "FSISPH" else: hydroname = "SPH" if asph: hydroname = "A" + hydroname -dataDir = os.path.join(dataRoot, +dataDir = os.path.join(dataDirBase, hydroname, "nperh=%4.2f" % nPerh, "XSPH=%s" % XSPH, @@ -159,12 +181,17 @@ #------------------------------------------------------------------------------- # Create a NodeList and associated Neighbor object. #------------------------------------------------------------------------------- -nodes1 = makeFluidNodeList("nodes1", eos, - hmin = hmin, - hmax = hmax, - nPerh = nPerh, - kernelExtent = kernelExtent, - rhoMin = rhomin) +if solid: + nodeListConstructor = makeSolidNodeList +else: + nodeListConstructor = makeFluidNodeList + +nodes1 = nodeListConstructor("nodes1", eos, + hmin = hmin, + hmax = hmax, + nPerh = nPerh, + kernelExtent = kernelExtent, + rhoMin = rhomin) #------------------------------------------------------------------------------- # Set the node properties. @@ -236,24 +263,51 @@ #------------------------------------------------------------------------------- if crksph: hydro = CRKSPH(dataBase = db, - W = WT, + order = correctionOrder, filter = filter, cfl = cfl, compatibleEnergyEvolution = compatibleEnergy, XSPH = XSPH, - correctionOrder = correctionOrder, densityUpdate = densityUpdate, HUpdate = HUpdate, ASPH = asph) +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction,RiemannGradientType) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) elif psph: - hydro = PASPH(dataBase = db, + hydro = PSPH(dataBase = db, W = WT, filter = filter, cfl = cfl, useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, compatibleEnergyEvolution = compatibleEnergy, evolveTotalEnergy = evolveTotalEnergy, - HopkinsConductivity = HopkinsConductivity, correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, HUpdate = HUpdate, @@ -269,38 +323,36 @@ correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, XSPH = XSPH, - HUpdate = HEvolution, + HUpdate = HUpdate, ASPH = asph) + +packages = [hydro] output("hydro") -output("hydro.kernel()") -output("hydro.PiKernel()") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.XSPH") -output("hydro.densityUpdate") output("hydro.HEvolution") -q = hydro.Q -output("q") -output("q.Cl") -output("q.Cq") -output("q.epsilon2") -output("q.limiter") -output("q.balsaraShearCorrection") -output("q.linearInExpansion") -output("q.quadraticInExpansion") - -packages = [hydro] - -#------------------------------------------------------------------------------- -# Construct the MMRV physics object. -#------------------------------------------------------------------------------- -if boolReduceViscosity: - evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: - evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) - packages.append(evolveCullenViscosityMultiplier) +if not gsph: + q = hydro.Q + output("q") + output("q.Cl") + output("q.Cq") + output("q.epsilon2") + output("q.limiter") + output("q.balsaraShearCorrection") + output("q.linearInExpansion") + output("q.quadraticInExpansion") + + #------------------------------------------------------------------------------- + # Construct the MMRV physics object. + #------------------------------------------------------------------------------- + if boolReduceViscosity: + evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) + packages.append(evolveReducingViscosityMultiplier) + elif boolCullenViscosity: + evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) + packages.append(evolveCullenViscosityMultiplier) #------------------------------------------------------------------------------- # Create boundary conditions. diff --git a/tests/Hydro/Sedov/Sedov-spherical-3d.py b/tests/Hydro/Sedov/Sedov-spherical-3d.py index 6aae37dbb..7b13228a9 100644 --- a/tests/Hydro/Sedov/Sedov-spherical-3d.py +++ b/tests/Hydro/Sedov/Sedov-spherical-3d.py @@ -32,22 +32,46 @@ smoothSpikeScale = 0.5, gamma = 5.0/3.0, mu = 1.0, + rhomin = 1e-10, + + # kernel + HUpdate = IdealH, + hmin = 1e-15, + hmax = 1.0, + + # hydros + crksph = False, + psph = False, + gsph = False, + fsisph = False, + + # hydro parameters + asph = False, + solid = False, + XSPH = False, + evolveTotalEnergy = False, + compatibleEnergy = True, + gradhCorrection = True, + correctVelocityGradient = True, + densityUpdate = RigorousSumDensity, + filter = 0.0, + + # crksph parameters + correctionOrder = LinearOrder, + volumeType = RKSumVolume, + + # gsph parameters + RiemannGradientType = RiemannGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) + linearReconstruction = True, + # Artificial Viscosity + Qconstructor = MonaghanGingoldViscosity, Cl = 1.0, Cq = 0.75, epsilon2 = 1e-2, Qlimiter = False, balsaraCorrection = False, linearInExpansion = False, - - crksph = False, - psph = False, - asph = False, # Only for H evolution, not hydro algorithm - Qconstructor = MonaghanGingoldViscosity, - correctionOrder = LinearOrder, - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - HUpdate = IdealH, - filter = 0.0, boolReduceViscosity = False, nh = 5.0, aMin = 0.1, @@ -61,18 +85,11 @@ betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - HopkinsConductivity = False, # For PSPH - evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy - volumeType = RKSumVolume, # For CRK - - hmin = 1e-15, - hmax = 1.0, - cfl = 0.5, - useVelocityMagnitudeForDt = True, - XSPH = False, - rhomin = 1e-10, + # Integration IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.5, + useVelocityMagnitudeForDt = False, steps = None, goalTime = None, goalRadius = 0.8, @@ -80,22 +97,19 @@ dtMin = 1.0e-8, dtMax = None, dtGrowth = 2.0, - vizCycle = None, - vizTime = 0.1, maxSteps = None, statsStep = 1, smoothIters = 0, - HEvolution = IdealH, - compatibleEnergy = True, - gradhCorrection = True, - correctVelocityGradient = True, + # IO + vizCycle = None, + vizTime = 0.1, restoreCycle = -1, restartStep = 1000, graphics = True, clearDirectories = False, - dataRoot = "dumps-spherical-Sedov", + dataDirBase = "dumps-spherical-Sedov", outputFile = "None", ) @@ -132,12 +146,16 @@ hydroname += "A" if crksph: hydroname += "CRKSPH" +elif fsisph: + hydroname = "FSISPH" +elif gsph: + hydroname = "GSPH" elif psph: hydroname += "PSPH" else: hydroname += "SPH" -dataDir = os.path.join(dataRoot, +dataDir = os.path.join(dataDirBase, hydroname, "nperh=%4.2f" % nPerh, "XSPH=%s" % XSPH, @@ -182,7 +200,11 @@ #------------------------------------------------------------------------------- # Create a NodeList and associated Neighbor object. #------------------------------------------------------------------------------- -nodes1 = makeFluidNodeList("nodes1", eos, +if solid: + nodeListConstructor = makeSolidNodeList +else: + nodeListConstructor = makeFluidNodeList +nodes1 = nodeListConstructor("nodes1", eos, hmin = hmin, hmax = hmax, kernelExtent = kernelExtent, @@ -209,7 +231,8 @@ SPH = (not asph)) if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes3d + #from VoronoiDistributeNodes import distributeNodes3d + from PeanoHilbertDistributeNodes import distributeNodes3d else: from DistributeNodes import distributeNodes3d @@ -279,18 +302,19 @@ #------------------------------------------------------------------------------- # Construct the artificial viscosity. #------------------------------------------------------------------------------- -q = Qconstructor(Cl, Cq, linearInExpansion) -q.epsilon2 = epsilon2 -q.limiter = Qlimiter -q.balsaraShearCorrection = balsaraCorrection -output("q") -output("q.Cl") -output("q.Cq") -output("q.epsilon2") -output("q.limiter") -output("q.balsaraShearCorrection") -output("q.linearInExpansion") -output("q.quadraticInExpansion") +if not gsph: + q = Qconstructor(Cl, Cq, linearInExpansion) + q.epsilon2 = epsilon2 + q.limiter = Qlimiter + q.balsaraShearCorrection = balsaraCorrection + output("q") + output("q.Cl") + output("q.Cq") + output("q.epsilon2") + output("q.limiter") + output("q.balsaraShearCorrection") + output("q.linearInExpansion") + output("q.quadraticInExpansion") #------------------------------------------------------------------------------- # Construct the hydro physics object. @@ -302,9 +326,38 @@ cfl = cfl, compatibleEnergyEvolution = compatibleEnergy, XSPH = XSPH, - correctionOrder = correctionOrder, + order = correctionOrder, densityUpdate = densityUpdate, HUpdate = HUpdate) +elif fsisph: + hydro = FSISPH(dataBase = db, + W = WT, + filter = filter, + cfl = cfl, + interfaceMethod = ModulusInterface, + sumDensityNodeLists=[nodes1], + densityStabilizationCoefficient = 0.00, + useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + correctVelocityGradient = correctVelocityGradient, + HUpdate = HUpdate) +elif gsph: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction,RiemannGradientType) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) elif psph: hydro = PSPH(dataBase = db, W = WT, @@ -314,7 +367,6 @@ useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, compatibleEnergyEvolution = compatibleEnergy, evolveTotalEnergy = evolveTotalEnergy, - HopkinsConductivity = HopkinsConductivity, correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, HUpdate = HUpdate, @@ -330,7 +382,7 @@ correctVelocityGradient = correctVelocityGradient, densityUpdate = densityUpdate, XSPH = XSPH, - HUpdate = HEvolution) + HUpdate = HUpdate) output("hydro") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") From c9b7f159f87bf617e16fa46444fcb1b4917443ef Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 18 Oct 2021 11:24:56 -0700 Subject: [PATCH 5/6] bugfix: fsisph bool flag got deleted somehow --- tests/Hydro/Noh/Noh-cylindrical-2d.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/Hydro/Noh/Noh-cylindrical-2d.py b/tests/Hydro/Noh/Noh-cylindrical-2d.py index e0fb2afd3..b76da0ab1 100644 --- a/tests/Hydro/Noh/Noh-cylindrical-2d.py +++ b/tests/Hydro/Noh/Noh-cylindrical-2d.py @@ -16,8 +16,8 @@ # # 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) -#ATS:gsph1 = testif(sph0, SELF, "--gsph True --nRadial 100 --cfl 0.25 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical GSPH, nPerh=2.0, restart test", np=8) +#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) +#ATS:gsph1 = testif(gsph0, SELF, "--gsph True --nRadial 100 --cfl 0.25 --nPerh 2.01 --graphics False --restartStep 20 --clearDirectories False --steps 60 --restoreCycle 40 --checkRestart True", label="Noh cylindrical GSPH, nPerh=2.0, restart test", np=8) #------------------------------------------------------------------------------- @@ -65,6 +65,7 @@ svph = False, crksph = False, psph = False, + fsisph = False, gsph = False, asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. boolReduceViscosity = False, From c4e65e2fa814627c96bc36c410f33e8bfaa58c5d Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 19 Oct 2021 11:33:06 -0700 Subject: [PATCH 6/6] fixing smoothing spike edit --- tests/Hydro/Sedov/Sedov-planar-1d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/Hydro/Sedov/Sedov-planar-1d.py b/tests/Hydro/Sedov/Sedov-planar-1d.py index 4e4313854..bc45c7fbd 100644 --- a/tests/Hydro/Sedov/Sedov-planar-1d.py +++ b/tests/Hydro/Sedov/Sedov-planar-1d.py @@ -21,7 +21,7 @@ rho0 = 1.0, eps0 = 0.0, Espike = 1.0, - smoothSpike = False, + smoothSpike = True, topHatSpike = False, smoothSpikeScale = 0.5, gamma = 5.0/3.0,