diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index fe35d2250..ea2435e59 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -6,6 +6,7 @@ Version vYYYY.MM.p -- Release date YYYY-MM-DD Notable changes include: * New features/ API changes: + * added MFV hydro from Hopkins 2015 with extension for ALE options * Build changes / improvements: * tpl-manager.py will no longer use generic x86_64 configs for non LC systems. Users will be required to supply their own configs for pointing spack at external packages. diff --git a/src/GSPH/CMakeLists.txt b/src/GSPH/CMakeLists.txt index 9b095b88e..09d04055b 100644 --- a/src/GSPH/CMakeLists.txt +++ b/src/GSPH/CMakeLists.txt @@ -4,23 +4,30 @@ set(GSPH_inst computeSPHVolume computeSumVolume computeMFMDensity + initializeGradients GenericRiemannHydro GSPHHydroBase MFMHydroBase + MFVHydroBase + Policies/MassFluxPolicy Policies/ReplaceWithRatioPolicy + Policies/MFVIncrementVelocityPolicy + Policies/MFVIncrementSpecificThermalEnergyPolicy + Policies/CompatibleMFVSpecificThermalEnergyPolicy Limiters/LimiterBase Limiters/VanLeerLimiter Limiters/SuperbeeLimiter Limiters/MinModLimiter Limiters/VanAlbaLimiter Limiters/OspreLimiter + Limiters/BarthJespersenLimiter WaveSpeeds/WaveSpeedBase WaveSpeeds/AcousticWaveSpeed WaveSpeeds/DavisWaveSpeed WaveSpeeds/EinfeldtWaveSpeed RiemannSolvers/RiemannSolverBase RiemannSolvers/HLLC - RiemannSolvers/GHLLC) + RiemannSolvers/SecondOrderArtificialViscosity) set(GSPH_sources GSPHFieldNames.cc) @@ -29,24 +36,31 @@ set(GSPH_headers computeSPHVolume.hh computeSumVolume.hh computeMFMDensity.hh + initializeGradients.hh GSPHFieldNames.hh GenericRiemannHydro.hh GSPHHydroBase.hh MFMHydroBase.hh + MFVHydroBase.hh + Policies/MassFluxPolicy.hh Policies/ReplaceWithRatioPolicy.hh + Policies/MFVIncrementVelocityPolicy.hh + Policies/MFVIncrementSpecificThermalEnergyPolicy.hh + Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh Limiters/LimiterBase.hh Limiters/VanLeerLimiter.hh Limiters/SuperbeeLimiter.hh Limiters/MinModLimiter.hh Limiters/VanAlbaLimiter.hh Limiters/OspreLimiter.hh + Limiters/BarthJespersenLimiter.hh WaveSpeeds/WaveSpeedBase.hh WaveSpeeds/AcousticWaveSpeed.hh WaveSpeeds/DavisWaveSpeed.hh WaveSpeeds/EinfeldtWaveSpeed.hh RiemannSolvers/RiemannSolverBase.hh RiemannSolvers/HLLC.hh - RiemannSolvers/GHLLC.hh) + RiemannSolvers/SecondOrderArtificialViscosity.hh) add_subdirectory(Limiters) add_subdirectory(WaveSpeeds) diff --git a/src/GSPH/GSPHEvaluateDerivatives.cc b/src/GSPH/GSPHEvaluateDerivatives.cc index 33e16a1a9..3f464732e 100644 --- a/src/GSPH/GSPHEvaluateDerivatives.cc +++ b/src/GSPH/GSPHEvaluateDerivatives.cc @@ -67,6 +67,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, // Derivative FieldLists. const auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + const auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); auto normalization = derivatives.fields(HydroFieldNames::normalization, 0.0); auto DxDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::position, Vector::zero); auto DrhoDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::massDensity, 0.0); @@ -83,6 +84,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); + CHECK(DrhoDx.size() == numNodeLists); CHECK(M.size() == numNodeLists); CHECK(normalization.size() == numNodeLists); CHECK(DxDt.size() == numNodeLists); @@ -139,7 +141,6 @@ evaluateDerivatives(const typename Dimension::Scalar time, const auto& vi = velocity(nodeListi, i); const auto& rhoi = massDensity(nodeListi, i); const auto& voli = volume(nodeListi, i); - //const auto& epsi = specificThermalEnergy(nodeListi, i); const auto& Pi = pressure(nodeListi, i); const auto& Hi = H(nodeListi, i); const auto& ci = soundSpeed(nodeListi, i); @@ -158,6 +159,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& weightedNeighborSumi = weightedNeighborSum_thread(nodeListi, i); auto& massSecondMomenti = massSecondMoment_thread(nodeListi, i); auto& XSPHDeltaVi = XSPHDeltaV_thread(nodeListi,i); + const auto& gradRhoi = DrhoDx(nodeListi,i); const auto& Mi = M(nodeListi,i); @@ -169,7 +171,6 @@ evaluateDerivatives(const typename Dimension::Scalar time, const auto& vj = velocity(nodeListj, j); const auto& rhoj = massDensity(nodeListj, j); const auto& volj = volume(nodeListj, j); - //const auto& epsj = specificThermalEnergy(nodeListj, j); const auto& Pj = pressure(nodeListj, j); const auto& Hj = H(nodeListj, j); const auto& cj = soundSpeed(nodeListj, j); @@ -188,8 +189,9 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& weightedNeighborSumj = weightedNeighborSum_thread(nodeListj, j); auto& massSecondMomentj = massSecondMoment_thread(nodeListj, j); auto& XSPHDeltaVj = XSPHDeltaV_thread(nodeListj,j); + const auto& gradRhoj = DrhoDx(nodeListj,j); const auto& Mj = M(nodeListj,j); - + // Node displacement. const auto rij = ri - rj; const auto rhatij =rij.unitVector(); @@ -233,26 +235,28 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto gradPj = riemannDpDxj; auto gradVi = riemannDvDxi; auto gradVj = riemannDvDxj; - if (gradType==GradientType::SPHSameTimeGradient){ - gradPi = newRiemannDpDxi; - gradPj = newRiemannDpDxj; - gradVi = newRiemannDvDxi; - gradVj = newRiemannDvDxj; + if (gradType==GradientType::SPHSameTimeGradient or + gradType==GradientType::SPHUncorrectedGradient){ + gradPi = newRiemannDpDx(nodeListi,i); + gradPj = newRiemannDpDx(nodeListj,j); + gradVi = newRiemannDvDx(nodeListi,i); + gradVj = newRiemannDvDx(nodeListj,j); } - riemannSolver.interfaceState(i, j, - nodeListi, nodeListj, - ri, rj, + + riemannSolver.interfaceState(ri, rj, + Hi, Hj, rhoi, rhoj, ci, cj, Peffi, Peffj, vi, vj, + gradRhoi, gradRhoj, gradPi, gradPj, gradVi, gradVj, - Pstar, - vstar, - rhostari, - rhostarj); - + Pstar, //output + vstar, //output + rhostari, //output + rhostarj); //output + // get our basis function and interface area vectors //-------------------------------------------------------- psii = volj*Wi; @@ -266,17 +270,16 @@ evaluateDerivatives(const typename Dimension::Scalar time, // acceleration //------------------------------------------------------ const auto deltaDvDt = Pstar*(Ai+Aj); - DvDti -= deltaDvDt; - DvDtj += deltaDvDt; + DvDti -= deltaDvDt/mi; + DvDtj += deltaDvDt/mj; // energy //------------------------------------------------------ const auto deltaDepsDti = 2.0*Pstar*Ai.dot(vi-vstar); const auto deltaDepsDtj = 2.0*Pstar*Aj.dot(vstar-vj); + DepsDti += deltaDepsDti/mi; + DepsDtj += deltaDepsDtj/mj; - DepsDti += deltaDepsDti; - DepsDtj += deltaDepsDtj; - if(compatibleEnergy){ const auto invmij = 1.0/(mi*mj); pairAccelerations[kk] = deltaDvDt*invmij; @@ -376,9 +379,6 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& weightedNeighborSumi = weightedNeighborSum(nodeListi, i); auto& massSecondMomenti = massSecondMoment(nodeListi, i); - DvDti /= mi; - DepsDti /= mi; - normi += voli*Hdeti*W0; DrhoDti = - rhoi * DvDxi.Trace(); @@ -434,23 +434,26 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, const DataBase& dataBase, const State& state, StateDerivatives& derivatives) const { - // The kernels and such. const auto& W = this->kernel(); - const auto gradType = this->gradientType(); - + const auto calcSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient + or this->gradientType() == GradientType::SPHUncorrectedGradient); + const auto correctSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient); + // The connectivity. const auto& connectivityMap = dataBase.connectivityMap(); const auto& nodeLists = connectivityMap.nodeLists(); const auto numNodeLists = nodeLists.size(); // Get the state and derivative FieldLists. + const auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); const auto volume = state.fields(HydroFieldNames::volume, 0.0); const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); const auto pressure = state.fields(HydroFieldNames::pressure, 0.0); const auto position = state.fields(HydroFieldNames::position, Vector::zero); const auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + CHECK(massDensity.size() == numNodeLists); CHECK(volume.size() == numNodeLists); CHECK(velocity.size() == numNodeLists); CHECK(pressure.size() == numNodeLists); @@ -458,10 +461,12 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(H.size() == numNodeLists); auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); CHECK(M.size() == numNodeLists); + CHECK(DrhoDx.size() == numNodeLists); CHECK(newRiemannDpDx.size() == numNodeLists); CHECK(newRiemannDvDx.size() == numNodeLists); @@ -476,6 +481,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, typename SpheralThreads::FieldListStack threadStack; auto M_thread = M.threadCopy(threadStack); + auto DrhoDx_thread = DrhoDx.threadCopy(threadStack); auto newRiemannDpDx_thread = newRiemannDpDx.threadCopy(threadStack); auto newRiemannDvDx_thread = newRiemannDvDx.threadCopy(threadStack); @@ -487,6 +493,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, nodeListj = pairs[kk].j_list; // Get the state for node i. + const auto& rhoi = massDensity(nodeListi, i); const auto& ri = position(nodeListi, i); const auto& voli = volume(nodeListi, i); const auto& Hi = H(nodeListi, i); @@ -495,8 +502,10 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(Hdeti > 0.0); auto& Mi = M_thread(nodeListi, i); + auto& DrhoDxi = DrhoDx_thread(nodeListi, i); // Get the state for node j + const auto& rhoj = massDensity(nodeListj, j); const auto& rj = position(nodeListj, j); const auto& volj = volume(nodeListj, j); const auto& Hj = H(nodeListj, j); @@ -505,6 +514,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(Hdetj > 0.0); auto& Mj = M_thread(nodeListj, j); + auto& DrhoDxj = DrhoDx_thread(nodeListj, j); const auto rij = ri - rj; @@ -530,12 +540,17 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, Mi -= rij.dyad(gradPsii); Mj -= rij.dyad(gradPsij); + DrhoDxi -= (rhoi - rhoj) * gradPsii; + DrhoDxj -= (rhoi - rhoj) * gradPsij; + // // based on nodal values - if (gradType == GradientType::SPHSameTimeGradient){ + if (calcSpatialGradients){ + const auto& vi = velocity(nodeListi, i); const auto& Pi = pressure(nodeListi, i); const auto& vj = velocity(nodeListj, j); const auto& Pj = pressure(nodeListj, j); + auto& newRiemannDpDxi = newRiemannDpDx_thread(nodeListi, i); auto& newRiemannDvDxi = newRiemannDvDx_thread(nodeListi, i); auto& newRiemannDpDxj = newRiemannDpDx_thread(nodeListj, j); @@ -546,6 +561,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, newRiemannDvDxi -= (vi-vj).dyad(gradPsii); newRiemannDvDxj -= (vi-vj).dyad(gradPsij); + } } // loop over pairs @@ -554,7 +570,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, } // OpenMP parallel region - // Finish up the spatial gradient calculation + // loop the nodes to finish up the spatial gradient calculation for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { const auto& nodeList = M[nodeListi]->nodeList(); const auto ni = nodeList.numInternalNodes(); @@ -562,6 +578,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, for (auto i = 0u; i < ni; ++i) { const auto numNeighborsi = connectivityMap.numNeighborsForNode(nodeListi, i); auto& Mi = M(nodeListi, i); + auto& DrhoDxi = DrhoDx(nodeListi, i); const auto Mdeti = std::abs(Mi.Determinant()); @@ -570,25 +587,28 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, Mi = ( goodM ? Mi.Inverse() : Tensor::one); - if (gradType == GradientType::SPHSameTimeGradient){ + DrhoDxi = Mi.Transpose()*DrhoDxi; + + if (correctSpatialGradients){ + auto& newRiemannDpDxi = newRiemannDpDx(nodeListi, i); auto& newRiemannDvDxi = newRiemannDvDx(nodeListi, i); - newRiemannDpDxi = Mi.Transpose()*newRiemannDpDxi; newRiemannDvDxi = newRiemannDvDxi*Mi; - } - } - - } + + } // if correctSpatialGradients + } // for each node + } // for each node list for (ConstBoundaryIterator boundItr = this->boundaryBegin(); boundItr != this->boundaryEnd(); ++boundItr)(*boundItr)->applyFieldListGhostBoundary(M); - if (gradType == GradientType::SPHSameTimeGradient){ + if (calcSpatialGradients){ for (ConstBoundaryIterator boundItr = this->boundaryBegin(); boundItr != this->boundaryEnd(); ++boundItr){ + (*boundItr)->applyFieldListGhostBoundary(DrhoDx); (*boundItr)->applyFieldListGhostBoundary(newRiemannDpDx); (*boundItr)->applyFieldListGhostBoundary(newRiemannDvDx); } @@ -597,6 +617,5 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, boundaryItr != this->boundaryEnd(); ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary(); -} - +} // MC correction method } // spheral namespace diff --git a/src/GSPH/GSPHFieldNames.cc b/src/GSPH/GSPHFieldNames.cc index 46d5b6b83..5b7d9dbdb 100644 --- a/src/GSPH/GSPHFieldNames.cc +++ b/src/GSPH/GSPHFieldNames.cc @@ -6,9 +6,13 @@ #include "GSPHFieldNames.hh" +const std::string Spheral::GSPHFieldNames::nodalVelocity = "velocity of node"; +const std::string Spheral::GSPHFieldNames::momentum = "momentum"; +const std::string Spheral::GSPHFieldNames::thermalEnergy = "thermal energy"; const std::string Spheral::GSPHFieldNames::densityGradient = "density gradient"; const std::string Spheral::GSPHFieldNames::pressureGradient = "pressure gradient"; const std::string Spheral::GSPHFieldNames::deviatoricStressTensorGradient = "deviatoric stress tensor gradient"; const std::string Spheral::GSPHFieldNames::RiemannPressureGradient = "Riemann solvers pressure gradient"; const std::string Spheral::GSPHFieldNames::RiemannVelocityGradient = "Riemann solvers velocity gradient"; -const std::string Spheral::GSPHFieldNames::RiemannDeviatoricStressTensorGradient = "Riemann solvers deviatoric stress tensor gradient"; \ No newline at end of file +const std::string Spheral::GSPHFieldNames::RiemannDeviatoricStressTensorGradient = "Riemann solvers deviatoric stress tensor gradient"; +const std::string Spheral::GSPHFieldNames::pairMassFlux = "pairwise mass flux"; \ No newline at end of file diff --git a/src/GSPH/GSPHFieldNames.hh b/src/GSPH/GSPHFieldNames.hh index 07421567f..d87846995 100644 --- a/src/GSPH/GSPHFieldNames.hh +++ b/src/GSPH/GSPHFieldNames.hh @@ -12,12 +12,16 @@ namespace Spheral { struct GSPHFieldNames { + static const std::string nodalVelocity; + static const std::string momentum; + static const std::string thermalEnergy; static const std::string densityGradient; static const std::string pressureGradient; static const std::string deviatoricStressTensorGradient; static const std::string RiemannPressureGradient; static const std::string RiemannVelocityGradient; static const std::string RiemannDeviatoricStressTensorGradient; + static const std::string pairMassFlux; }; } diff --git a/src/GSPH/GSPHHydroBase.cc b/src/GSPH/GSPHHydroBase.cc index e4e1679eb..ee809189c 100644 --- a/src/GSPH/GSPHHydroBase.cc +++ b/src/GSPH/GSPHHydroBase.cc @@ -1,5 +1,8 @@ //---------------------------------Spheral++----------------------------------// -// GSPHHydroBase -- The Godunov SPH hydrodynamic package for Spheral++. +// GSPHHydroBase -- A Riemann-solver-based implementation of SPH. Compared to +// MFM/MFV this approach requires a larger neighbor set. 2.5 +// nodes per kernel extent instead of 2-2.25 for MFM/MFV but +// does perform better on certain tests (Noh implosion) // // J.M. Pearl 2021 //----------------------------------------------------------------------------// @@ -207,11 +210,8 @@ initialize(const typename Dimension::Scalar time, State& state, StateDerivatives& derivs) { TIME_BEGIN("GSPHinitialize"); - GenericRiemannHydro::initialize(time,dt,dataBase,state,derivs); - TIME_END("GSPHinitialize"); - } //------------------------------------------------------------------------------ @@ -239,9 +239,7 @@ GSPHHydroBase:: applyGhostBoundaries(State& state, StateDerivatives& derivs) { TIME_BEGIN("GSPHghostBounds"); - GenericRiemannHydro::applyGhostBoundaries(state,derivs); - TIME_END("GSPHghostBounds"); } @@ -254,9 +252,7 @@ GSPHHydroBase:: enforceBoundaries(State& state, StateDerivatives& derivs) { TIME_BEGIN("GSPHenforceBounds"); - GenericRiemannHydro::enforceBoundaries(state,derivs); - TIME_END("GSPHenforceBounds"); } diff --git a/src/GSPH/GSPHHydroBase.hh b/src/GSPH/GSPHHydroBase.hh index 964c70472..0137c768b 100644 --- a/src/GSPH/GSPHHydroBase.hh +++ b/src/GSPH/GSPHHydroBase.hh @@ -1,5 +1,8 @@ //---------------------------------Spheral++----------------------------------// -// GSPHHydroBase -- The Godunov SPH hydrodynamic package for Spheral++. +// GSPHHydroBase -- A Riemann-solver-based implementation of SPH. Compared to +// MFM/MFV this approach requires a larger neighbor set. 2.5 +// nodes per kernel extent instead of 2-2.25 for MFM/MFV but +// does perform better on certain tests (Noh implosion) // // J.M. Pearl 2021 //----------------------------------------------------------------------------// diff --git a/src/GSPH/GSPHHydros.py b/src/GSPH/GSPHHydros.py index 618871580..31166132b 100644 --- a/src/GSPH/GSPHHydros.py +++ b/src/GSPH/GSPHHydros.py @@ -3,118 +3,6 @@ from spheralDimensions import spheralDimensions dims = spheralDimensions() -#------------------------------------------------------------------------------- -# density-based GSPH factory string -#------------------------------------------------------------------------------- -GSPHHydroFactoryString = """ -class %(classname)s%(dim)s(GSPHHydroBase%(dim)s): - - def __init__(self, - dataBase, - riemannSolver, - W, - epsDiffusionCoeff = 0.0, - cfl = 0.25, - useVelocityMagnitudeForDt = False, - compatibleEnergyEvolution = True, - evolveTotalEnergy = False, - XSPH = True, - correctVelocityGradient = True, - gradientType = HydroAccelerationGradient, - densityUpdate = IntegrateDensity, - HUpdate = IdealH, - epsTensile = 0.0, - nTensile = 4.0, - xmin = Vector%(dim)s(-1e100, -1e100, -1e100), - xmax = Vector%(dim)s( 1e100, 1e100, 1e100)): - self._smoothingScaleMethod = %(smoothingScaleMethod)s%(dim)s() - GSPHHydroBase%(dim)s.__init__(self, - self._smoothingScaleMethod, - dataBase, - riemannSolver, - W, - epsDiffusionCoeff, - cfl, - useVelocityMagnitudeForDt, - compatibleEnergyEvolution, - evolveTotalEnergy, - XSPH, - correctVelocityGradient, - gradientType, - densityUpdate, - HUpdate, - epsTensile, - nTensile, - xmin, - xmax) - return -""" - -#------------------------------------------------------------------------------- -# volume-based GSPH factory string (MFM) -#------------------------------------------------------------------------------- -MFMHydroFactoryString = """ -class %(classname)s%(dim)s(MFMHydroBase%(dim)s): - - def __init__(self, - dataBase, - riemannSolver, - W, - epsDiffusionCoeff = 0.0, - cfl = 0.25, - useVelocityMagnitudeForDt = False, - compatibleEnergyEvolution = True, - evolveTotalEnergy = False, - XSPH = True, - correctVelocityGradient = True, - gradientType = HydroAccelerationGradient, - densityUpdate = IntegrateDensity, - HUpdate = IdealH, - epsTensile = 0.0, - nTensile = 4.0, - xmin = Vector%(dim)s(-1e100, -1e100, -1e100), - xmax = Vector%(dim)s( 1e100, 1e100, 1e100)): - self._smoothingScaleMethod = %(smoothingScaleMethod)s%(dim)s() - MFMHydroBase%(dim)s.__init__(self, - self._smoothingScaleMethod, - dataBase, - riemannSolver, - W, - epsDiffusionCoeff, - cfl, - useVelocityMagnitudeForDt, - compatibleEnergyEvolution, - evolveTotalEnergy, - XSPH, - correctVelocityGradient, - gradientType, - densityUpdate, - HUpdate, - epsTensile, - nTensile, - xmin, - xmax) - return -""" - -#------------------------------------------------------------------------------- -# Make 'em. -#------------------------------------------------------------------------------- -for dim in dims: - exec(GSPHHydroFactoryString % {"dim" : "%id" % dim, - "classname" : "GSPHHydro", - "smoothingScaleMethod" : "SPHSmoothingScale"}) - exec(GSPHHydroFactoryString % {"dim" : "%id" % dim, - "classname" : "AGSPHHydro", - "smoothingScaleMethod" : "ASPHSmoothingScale"}) - - exec(MFMHydroFactoryString % {"dim" : "%id" % dim, - "classname" : "MFMHydro", - "smoothingScaleMethod" : "SPHSmoothingScale"}) - exec(MFMHydroFactoryString % {"dim" : "%id" % dim, - "classname" : "AMFMHydro", - "smoothingScaleMethod" : "ASPHSmoothingScale"}) - #------------------------------------------------------------------------------- # GSPH convience wrapper function #------------------------------------------------------------------------------- @@ -155,10 +43,13 @@ def GSPH(dataBase, print(" which will result in fluid behaviour for those nodes.") raise RuntimeError("Cannot mix solid and fluid NodeLists.") + Constructor = eval("GSPHHydroBase%id" % ndim) + + # Smoothing scale update if ASPH: - Constructor = eval("AGSPHHydro%id" % ndim) + smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim) else: - Constructor = eval("GSPHHydro%id" % ndim) + smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim) if riemannSolver is None: waveSpeedMethod = eval("DavisWaveSpeed%id()" % (ndim)) @@ -170,7 +61,9 @@ def GSPH(dataBase, xmin = (ndim,) + xmin xmax = (ndim,) + xmax - kwargs = {"riemannSolver" : riemannSolver, + kwargs = {"smoothingScaleMethod" : smoothingScaleMethod, + "dataBase" : dataBase, + "riemannSolver" : riemannSolver, "W" : W, "epsDiffusionCoeff" : specificThermalEnergyDiffusionCoefficient, "dataBase" : dataBase, @@ -180,7 +73,7 @@ def GSPH(dataBase, "evolveTotalEnergy" : evolveTotalEnergy, "XSPH" : XSPH, "correctVelocityGradient" : correctVelocityGradient, - "gradientType" : gradientType, + "gradType" : gradientType, "densityUpdate" : densityUpdate, "HUpdate" : HUpdate, "epsTensile" : epsTensile, @@ -191,8 +84,10 @@ def GSPH(dataBase, # Build and return the thing. result = Constructor(**kwargs) - return result + result._smoothingScaleMethod = smoothingScaleMethod + return result + #------------------------------------------------------------------------------- # MFM convience wrapper function #------------------------------------------------------------------------------- @@ -233,10 +128,14 @@ def MFM(dataBase, print(" which will result in fluid behaviour for those nodes.") raise RuntimeError("Cannot mix solid and fluid NodeLists.") + Constructor = eval("MFMHydroBase%id" % ndim) + + # Smoothing scale update if ASPH: - Constructor = eval("AMFMHydro%id" % ndim) + smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim) else: - Constructor = eval("MFMHydro%id" % ndim) + smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim) + if riemannSolver is None: waveSpeedMethod = eval("DavisWaveSpeed%id()" % (ndim)) @@ -248,7 +147,9 @@ def MFM(dataBase, xmin = (ndim,) + xmin xmax = (ndim,) + xmax - kwargs = {"riemannSolver" : riemannSolver, + kwargs = {"smoothingScaleMethod" : smoothingScaleMethod, + "dataBase" : dataBase, + "riemannSolver" : riemannSolver, "W" : W, "epsDiffusionCoeff" : specificThermalEnergyDiffusionCoefficient, "dataBase" : dataBase, @@ -258,7 +159,7 @@ def MFM(dataBase, "evolveTotalEnergy" : evolveTotalEnergy, "XSPH" : XSPH, "correctVelocityGradient" : correctVelocityGradient, - "gradientType" : gradientType, + "gradType" : gradientType, "densityUpdate" : densityUpdate, "HUpdate" : HUpdate, "epsTensile" : epsTensile, @@ -269,5 +170,99 @@ def MFM(dataBase, # Build and return the thing. result = Constructor(**kwargs) + result._smoothingScaleMethod = smoothingScaleMethod + + return result + + + +#------------------------------------------------------------------------------- +# MFM convience wrapper function +#------------------------------------------------------------------------------- +def MFV(dataBase, + W, + riemannSolver=None, + specificThermalEnergyDiffusionCoefficient = 0.0, + cfl = 0.25, + nodeMotionCoefficient = 0.2, + nodeMotionType = NodeMotionType.Lagrangian, + gradientType = HydroAccelerationGradient, + densityUpdate = IntegrateDensity, + useVelocityMagnitudeForDt = False, + compatibleEnergyEvolution = True, + evolveTotalEnergy = False, + XSPH = False, + correctVelocityGradient = False, + HUpdate = IdealH, + epsTensile = 0.0, + nTensile = 4.0, + damageRelieveRubble = False, + negativePressureInDamage = False, + strengthInDamage = False, + xmin = (-1e100, -1e100, -1e100), + xmax = ( 1e100, 1e100, 1e100), + ASPH = False, + RZ = False): + + # for now we'll just piggy back off this enum + assert densityUpdate in (RigorousSumDensity,IntegrateDensity) + + # 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. + ndim = dataBase.nDim + nfluid = dataBase.numFluidNodeLists + nsolid = dataBase.numSolidNodeLists + if nsolid > 0 and nsolid != nfluid: + print("MFM Error: you have provided both solid and fluid NodeLists, which is currently not supported.") + print(" If you want some fluids active, provide SolidNodeList without a strength option specfied,") + print(" which will result in fluid behaviour for those nodes.") + raise RuntimeError("Cannot mix solid and fluid NodeLists.") + + Constructor = eval("MFVHydroBase%id" % ndim) + + # Smoothing scale update + if ASPH: + smoothingScaleMethod = eval("ASPHSmoothingScale%id()" % ndim) + else: + smoothingScaleMethod = eval("SPHSmoothingScale%id()" % ndim) + + + if riemannSolver is None: + waveSpeedMethod = eval("DavisWaveSpeed%id()" % (ndim)) + slopeLimiter = eval("VanLeerLimiter%id()" % (ndim)) + linearReconstruction=True + riemannSolver = eval("HLLC%id(slopeLimiter,waveSpeedMethod,linearReconstruction)" % (ndim)) + + # Build the constructor arguments + xmin = (ndim,) + xmin + xmax = (ndim,) + xmax + + kwargs = {"smoothingScaleMethod" : smoothingScaleMethod, + "dataBase" : dataBase, + "riemannSolver" : riemannSolver, + "W" : W, + "epsDiffusionCoeff" : specificThermalEnergyDiffusionCoefficient, + "dataBase" : dataBase, + "cfl" : cfl, + "useVelocityMagnitudeForDt" : useVelocityMagnitudeForDt, + "compatibleEnergyEvolution" : compatibleEnergyEvolution, + "evolveTotalEnergy" : evolveTotalEnergy, + "XSPH" : XSPH, + "correctVelocityGradient" : correctVelocityGradient, + "nodeMotionCoefficient" : nodeMotionCoefficient, + "nodeMotionType" : nodeMotionType, + "gradType" : gradientType, + "densityUpdate" : densityUpdate, + "HUpdate" : HUpdate, + "epsTensile" : epsTensile, + "nTensile" : nTensile, + "xmin" : eval("Vector%id(%g, %g, %g)" % xmin), + "xmax" : eval("Vector%id(%g, %g, %g)" % xmax)} + + #print(nodeMotionType) + # Build and return the thing. + result = Constructor(**kwargs) + result._smoothingScaleMethod = smoothingScaleMethod + #print(result.nodeMotionType) return result diff --git a/src/GSPH/GenericRiemannHydro.cc b/src/GSPH/GenericRiemannHydro.cc index 6c9ea80c9..b072e6f04 100644 --- a/src/GSPH/GenericRiemannHydro.cc +++ b/src/GSPH/GenericRiemannHydro.cc @@ -14,6 +14,7 @@ #include "DataBase/StateDerivatives.hh" #include "DataBase/IncrementState.hh" #include "DataBase/ReplaceState.hh" +#include "DataBase/PureReplaceState.hh" #include "DataBase/IncrementBoundedState.hh" #include "DataBase/ReplaceBoundedState.hh" #include "DataBase/updateStateFields.hh" @@ -35,6 +36,7 @@ #include "GSPH/GSPHFieldNames.hh" #include "GSPH/GenericRiemannHydro.hh" #include "GSPH/computeSPHVolume.hh" +#include "GSPH/initializeGradients.hh" #include "GSPH/RiemannSolvers/RiemannSolverBase.hh" #ifdef _OPENMP @@ -50,6 +52,7 @@ using std::string; using std::pair; using std::to_string; using std::make_pair; +using std::make_shared; namespace { @@ -119,9 +122,10 @@ GenericRiemannHydro(const SmoothingScaleBase& smoothingScaleMethod, mXSPHDeltaV(FieldStorageType::CopyFields), mM(FieldStorageType::CopyFields), mDxDt(FieldStorageType::CopyFields), - mDvDt(FieldStorageType::CopyFields), - mDspecificThermalEnergyDt(FieldStorageType::CopyFields), - mDHDt(FieldStorageType::CopyFields), + mDvDt(FieldStorageType::CopyFields), // move up one layer + mDspecificThermalEnergyDt(FieldStorageType::CopyFields), // move up one layer + mDHDt(FieldStorageType::CopyFields), + mDrhoDx(FieldStorageType::CopyFields), mDvDx(FieldStorageType::CopyFields), mRiemannDpDx(FieldStorageType::CopyFields), mRiemannDvDx(FieldStorageType::CopyFields), @@ -146,6 +150,7 @@ GenericRiemannHydro(const SmoothingScaleBase& smoothingScaleMethod, mDvDt = dataBase.newFluidFieldList(Vector::zero, HydroFieldNames::hydroAcceleration); mDspecificThermalEnergyDt = dataBase.newFluidFieldList(0.0, IncrementState::prefix() + HydroFieldNames::specificThermalEnergy); mDHDt = dataBase.newFluidFieldList(SymTensor::zero, IncrementState::prefix() + HydroFieldNames::H); + mDrhoDx = dataBase.newFluidFieldList(Vector::zero,GSPHFieldNames::densityGradient); mDvDx = dataBase.newFluidFieldList(Tensor::zero, HydroFieldNames::velocityGradient); mRiemannDpDx = dataBase.newFluidFieldList(Vector::zero,GSPHFieldNames::RiemannPressureGradient); mRiemannDvDx = dataBase.newFluidFieldList(Tensor::zero,GSPHFieldNames::RiemannVelocityGradient); @@ -174,16 +179,43 @@ initializeProblemStartupDependencies(DataBase& dataBase, State& state, StateDerivatives& derivs) { - // Set the moduli. + const auto& connectivityMap = dataBase.connectivityMap(); + const auto mass = dataBase.fluidMass(); + const auto massDensity = dataBase.fluidMassDensity(); + const auto position = dataBase.fluidPosition(); + const auto H = dataBase.fluidHfield(); + auto velocity = dataBase.fluidVelocity(); + updateStateFields(HydroFieldNames::pressure, state, derivs); updateStateFields(HydroFieldNames::soundSpeed, state, derivs); - // for now initialize with SPH volume to make sure things are defined - const auto mass = dataBase.fluidMass(); - const auto massDensity = dataBase.fluidMassDensity(); computeSPHVolume(mass,massDensity,mVolume); + + for (ConstBoundaryIterator boundItr = this->boundaryBegin(); + boundItr != this->boundaryEnd(); + ++boundItr){ + (*boundItr)->applyFieldListGhostBoundary(mVolume); + (*boundItr)->applyFieldListGhostBoundary(velocity); + (*boundItr)->applyFieldListGhostBoundary(mPressure); + } + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary(); + + initializeGradients(connectivityMap, + this->kernel(), + position, + H, + mVolume, + mPressure, + velocity, + mM, + mRiemannDpDx, + mRiemannDvDx); + } + //------------------------------------------------------------------------------ // Register the state we need/are going to evolve. //------------------------------------------------------------------------------ @@ -227,38 +259,42 @@ registerState(DataBase& dataBase, VERIFY2(false, "SPH ERROR: Unknown Hevolution option "); } } + + auto positionPolicy = make_shared>(); + auto pressurePolicy = make_shared>(); + auto csPolicy = make_shared>(); + auto pressureGradientPolicy = make_shared>(); + auto velocityGradientPolicy = make_shared>(); + auto velocityPolicy = make_policy>({HydroFieldNames::position,HydroFieldNames::specificThermalEnergy},true); // normal state variables state.enroll(mTimeStepMask); state.enroll(mVolume); state.enroll(mass); state.enroll(massDensity); - state.enroll(position, std::make_shared>()); - state.enroll(mPressure, std::make_shared>()); - state.enroll(mSoundSpeed, std::make_shared>()); - state.enroll(mRiemannDpDx); - state.enroll(mRiemannDvDx); + state.enroll(position, positionPolicy); + state.enroll(mPressure, pressurePolicy); + state.enroll(mSoundSpeed, csPolicy); + state.enroll(velocity, velocityPolicy); + + if (mRiemannSolver.linearReconstruction()){ + state.enroll(mRiemannDpDx, pressureGradientPolicy); + state.enroll(mRiemannDvDx, velocityGradientPolicy); + }else{ + state.enroll(mRiemannDpDx); + state.enroll(mRiemannDvDx); + } // conditional for energy method if (mCompatibleEnergyEvolution) { - - state.enroll(specificThermalEnergy, std::make_shared>(dataBase)); - state.enroll(velocity, make_policy>({HydroFieldNames::position, - HydroFieldNames::specificThermalEnergy}, - true)); - } else if (mEvolveTotalEnergy) { - - state.enroll(specificThermalEnergy, std::make_shared>()); - state.enroll(velocity, make_policy>({HydroFieldNames::position, - HydroFieldNames::specificThermalEnergy}, - true)); - + auto thermalEnergyPolicy = make_shared>(dataBase); + state.enroll(specificThermalEnergy, thermalEnergyPolicy); + }else if (mEvolveTotalEnergy) { + auto thermalEnergyPolicy = make_shared>(); + state.enroll(specificThermalEnergy, thermalEnergyPolicy); } else { - - state.enroll(specificThermalEnergy, std::make_shared>()); - state.enroll(velocity, make_policy>({HydroFieldNames::position}, - true)); - + auto thermalEnergyPolicy = make_shared>(); + state.enroll(specificThermalEnergy, thermalEnergyPolicy); } } @@ -285,6 +321,7 @@ registerDerivatives(DataBase& dataBase, dataBase.resizeFluidFieldList(mDspecificThermalEnergyDt, 0.0, IncrementState::prefix() + HydroFieldNames::specificThermalEnergy, false); dataBase.resizeFluidFieldList(mDHDt, SymTensor::zero, IncrementState::prefix() + HydroFieldNames::H, false); dataBase.resizeFluidFieldList(mDvDx, Tensor::zero, HydroFieldNames::velocityGradient, false); + dataBase.resizeFluidFieldList(mDrhoDx, Vector::zero, GSPHFieldNames::densityGradient, false); dataBase.resizeFluidFieldList(mM, Tensor::zero, HydroFieldNames::M_SPHCorrection, false); // Check if someone already registered DxDt. @@ -294,6 +331,7 @@ registerDerivatives(DataBase& dataBase, } // Check that no-one else is trying to control the hydro vote for DvDt. CHECK(not derivs.registered(mDvDt)); + derivs.enroll(mDrhoDx); derivs.enroll(mNewRiemannDpDx); derivs.enroll(mNewRiemannDvDx); derivs.enroll(mDvDt); @@ -488,51 +526,6 @@ initialize(const typename Dimension::Scalar time, const DataBase& dataBase, State& state, StateDerivatives& derivs) { - - auto& riemannSolver = this->riemannSolver(); - // const auto& W = this->kernel(); - - // riemannSolver.initialize(dataBase, - // state, - // derivs, - // this->boundaryBegin(), - // this->boundaryEnd(), - // time, - // dt, - // W); - - if(riemannSolver.linearReconstruction()){ - const auto& connectivityMap = dataBase.connectivityMap(); - const auto& nodeLists = connectivityMap.nodeLists(); - const auto numNodeLists = nodeLists.size(); - - // copy from previous time step - for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { - const auto& nodeList = nodeLists[nodeListi]; - const auto ni = nodeList->numInternalNodes(); - #pragma omp parallel for - for (auto i = 0u; i < ni; ++i) { - const auto DvDxi = mNewRiemannDvDx(nodeListi,i); - const auto DpDxi = mNewRiemannDpDx(nodeListi,i); - - mRiemannDvDx(nodeListi,i) = DvDxi; - mRiemannDpDx(nodeListi,i) = DpDxi; - - } - } - - for (auto boundItr =this->boundaryBegin(); - boundItr != this->boundaryEnd(); - ++boundItr) { - (*boundItr)->applyFieldListGhostBoundary(mRiemannDvDx); - (*boundItr)->applyFieldListGhostBoundary(mRiemannDpDx); - } - - for (auto boundItr = this->boundaryBegin(); - boundItr != this->boundaryEnd(); - ++boundItr) (*boundItr)->finalizeGhostBoundary(); - - } // if LinearReconstruction } @@ -568,8 +561,8 @@ template void GenericRiemannHydro:: applyGhostBoundaries(State& state, - StateDerivatives& /*derivs*/) { - // Apply boundary conditions to the basic fluid state Fields. + StateDerivatives& derivs) { + auto volume = state.fields(HydroFieldNames::volume, 0.0); auto mass = state.fields(HydroFieldNames::mass, 0.0); auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); @@ -577,8 +570,6 @@ applyGhostBoundaries(State& state, auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); auto pressure = state.fields(HydroFieldNames::pressure, 0.0); auto soundSpeed = state.fields(HydroFieldNames::soundSpeed, 0.0); - - // our store vars in the riemann solver auto DpDx = state.fields(GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto DvDx = state.fields(GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); @@ -605,9 +596,8 @@ template void GenericRiemannHydro:: enforceBoundaries(State& state, - StateDerivatives& /*derivs*/) { + StateDerivatives& derivs) { - // Enforce boundary conditions on the fluid state Fields. auto volume = state.fields(HydroFieldNames::volume, 0.0); auto mass = state.fields(HydroFieldNames::mass, 0.0); auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); @@ -615,12 +605,9 @@ enforceBoundaries(State& state, auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); auto pressure = state.fields(HydroFieldNames::pressure, 0.0); auto soundSpeed = state.fields(HydroFieldNames::soundSpeed, 0.0); - - // our store vars in the riemann solver auto DpDx = state.fields(GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto DvDx = state.fields(GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); - for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); boundaryItr != this->boundaryEnd(); ++boundaryItr) { @@ -646,6 +633,7 @@ GenericRiemannHydro:: dumpState(FileIO& file, const string& pathName) const { file.write(mTimeStepMask, pathName + "/timeStepMask"); + file.write(mVolume, pathName + "/volume"); file.write(mPressure, pathName + "/pressure"); file.write(mSoundSpeed, pathName + "/soundSpeed"); @@ -682,6 +670,7 @@ GenericRiemannHydro:: restoreState(const FileIO& file, const string& pathName) { file.read(mTimeStepMask, pathName + "/timeStepMask"); + file.read(mVolume, pathName + "/volume"); file.read(mPressure, pathName + "/pressure"); file.read(mSoundSpeed, pathName + "/soundSpeed"); diff --git a/src/GSPH/GenericRiemannHydro.hh b/src/GSPH/GenericRiemannHydro.hh index 971370eba..e4bd66862 100644 --- a/src/GSPH/GenericRiemannHydro.hh +++ b/src/GSPH/GenericRiemannHydro.hh @@ -1,6 +1,6 @@ //---------------------------------Spheral++----------------------------------// // GenericRiemannHydro -- pure virtual class for hydros using a Riemann -// solver +// solver // // J.M. Pearl 2022 //----------------------------------------------------------------------------// @@ -20,7 +20,15 @@ enum class GradientType { HydroAccelerationGradient = 1, SPHGradient = 2, MixedMethodGradient = 3, - SPHSameTimeGradient = 4 + SPHSameTimeGradient = 4, + SPHUncorrectedGradient = 5, + NoGradient = 6 +}; + +enum class GSPHEvolutionType { + IdealH = 0, + IntegrateH = 1, + constantNeighborCount = 2 }; template class State; @@ -206,6 +214,7 @@ public: const std::vector& pairAccelerations() const; const std::vector& pairDepsDt() const; + const FieldList& DrhoDx() const; const FieldList& riemannDpDx() const; const FieldList& riemannDvDx() const; const FieldList& newRiemannDpDx() const; @@ -232,7 +241,7 @@ private: MassDensityType mDensityUpdate; HEvolutionType mHEvolution; - // A bunch of switches. + // A bunch of switches. bool mCompatibleEnergyEvolution; bool mEvolveTotalEnergy; bool mXSPH; @@ -266,6 +275,7 @@ private: FieldList mDspecificThermalEnergyDt; FieldList mDHDt; + FieldList mDrhoDx; FieldList mDvDx; FieldList mRiemannDpDx; FieldList mRiemannDvDx; diff --git a/src/GSPH/GenericRiemannHydroInline.hh b/src/GSPH/GenericRiemannHydroInline.hh index ab296c78f..e3264732a 100644 --- a/src/GSPH/GenericRiemannHydroInline.hh +++ b/src/GSPH/GenericRiemannHydroInline.hh @@ -643,6 +643,14 @@ specificThermalEnergyDiffusionCoefficient() const { } +template +inline +const FieldList& +GenericRiemannHydro:: +DrhoDx() const { + return mDrhoDx; +} + template inline const FieldList& diff --git a/src/GSPH/Limiters/BarthJespersenLimiter.cc b/src/GSPH/Limiters/BarthJespersenLimiter.cc new file mode 100644 index 000000000..4c6259af5 --- /dev/null +++ b/src/GSPH/Limiters/BarthJespersenLimiter.cc @@ -0,0 +1,40 @@ +//---------------------------------Spheral++----------------------------------// +// BarthJespersenLimiter +// J. Barth, D. C. Jespersen, The design and application of upwind schemes +// on unstructured meshes, in: 27th Aerospace Sciences Meetings, AIAA Paper +// 89-0366, Reno, NV, 1989. doi:10.2514/6.1989-366 +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "BarthJespersenLimiter.hh" + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +template +BarthJespersenLimiter:: +BarthJespersenLimiter(): + LimiterBase(true,true){ +} + +//------------------------------------------------------------------------------ +// Destructor +//------------------------------------------------------------------------------ +template +BarthJespersenLimiter:: +~BarthJespersenLimiter(){} + +//------------------------------------------------------------------------------ +// slope limiter +//------------------------------------------------------------------------------ +template +typename Dimension::Scalar +BarthJespersenLimiter:: +fluxLimiter(const typename Dimension::Scalar x) const { + return std::min(std::min(0.5*(x+1),2.0),2.0*x) ; +} + +} \ No newline at end of file diff --git a/src/GSPH/Limiters/BarthJespersenLimiter.hh b/src/GSPH/Limiters/BarthJespersenLimiter.hh new file mode 100644 index 000000000..ed4ea0f39 --- /dev/null +++ b/src/GSPH/Limiters/BarthJespersenLimiter.hh @@ -0,0 +1,44 @@ +//---------------------------------Spheral++----------------------------------// +// BarthJespersenLimiter +// J. Barth, D. C. Jespersen, The design and application of upwind schemes +// on unstructured meshes, in: 27th Aerospace Sciences Meetings, AIAA Paper +// 89-0366, Reno, NV, 1989. doi:10.2514/6.1989-366 +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#ifndef __Spheral_BarthJespersenLimiter_hh__ +#define __Spheral_BarthJespersenLimiter_hh__ + +#include "LimiterBase.hh" + +namespace Spheral { + +template +class BarthJespersenLimiter : public LimiterBase { + +public: + + typedef typename Dimension::Scalar Scalar; + + BarthJespersenLimiter(); + + ~BarthJespersenLimiter(); + + virtual + Scalar fluxLimiter(const Scalar) const override; + +}; + + +} + +#else + +// Forward declaration. +namespace Spheral { + template class BarthJespersenLimiter; +} + +#endif + diff --git a/src/GSPH/Limiters/BarthJespersenLimiterInst.cc.py b/src/GSPH/Limiters/BarthJespersenLimiterInst.cc.py new file mode 100644 index 000000000..b7ec5c8d8 --- /dev/null +++ b/src/GSPH/Limiters/BarthJespersenLimiterInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "Geometry/Dimension.hh" +#include "GSPH/Limiters/BarthJespersenLimiter.cc" + +namespace Spheral { + template class BarthJespersenLimiter >; +} +""" diff --git a/src/GSPH/MFMEvaluateDerivatives.cc b/src/GSPH/MFMEvaluateDerivatives.cc index 75a5846a2..4e5098708 100644 --- a/src/GSPH/MFMEvaluateDerivatives.cc +++ b/src/GSPH/MFMEvaluateDerivatives.cc @@ -66,6 +66,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, // Derivative FieldLists. const auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + const auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); auto normalization = derivatives.fields(HydroFieldNames::normalization, 0.0); auto DxDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::position, Vector::zero); auto DvolDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::volume, 0.0); @@ -82,6 +83,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); + CHECK(DrhoDx.size() == numNodeLists); CHECK(M.size() == numNodeLists); CHECK(normalization.size() == numNodeLists); CHECK(DxDt.size() == numNodeLists); @@ -158,6 +160,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& massSecondMomenti = massSecondMoment_thread(nodeListi, i); auto& XSPHDeltaVi = XSPHDeltaV_thread(nodeListi,i); const auto& Mi = M(nodeListi,i); + const auto& gradRhoi = DrhoDx(nodeListi,i); // Get the state for node j @@ -188,10 +191,12 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& massSecondMomentj = massSecondMoment_thread(nodeListj, j); auto& XSPHDeltaVj = XSPHDeltaV_thread(nodeListj,j); const auto& Mj = M(nodeListj,j); + const auto& gradRhoj = DrhoDx(nodeListj,j); // Node displacement. const auto rij = ri - rj; const auto rhatij =rij.unitVector(); + //const auto rMagij = rij.magnitude2(); const auto vij = vi - vj; const auto etai = Hi*rij; const auto etaj = Hj*rij; @@ -218,7 +223,8 @@ evaluateDerivatives(const typename Dimension::Scalar time, weightedNeighborSumj += std::abs(gWj); massSecondMomenti += gradWi.magnitude2()*thpt; massSecondMomentj += gradWj.magnitude2()*thpt; - + //massSecondMomenti -= voli*rij.selfdyad()*gWi/rMagij;//.magnitude2()*thpt; + //massSecondMomentj -= volj*rij.selfdyad()*gWj/rMagij;//.magnitude2()*thpt; // Determine an effective pressure including a term to fight the tensile instability. //const auto fij = epsTensile*pow(Wi/(Hdeti*WnPerh), nTensile); const auto fij = epsTensile*FastMath::pow4(Wi/(Hdeti*WnPerh)); @@ -233,25 +239,26 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto gradPj = riemannDpDxj; auto gradVi = riemannDvDxi; auto gradVj = riemannDvDxj; - if (gradType==GradientType::SPHSameTimeGradient){ - gradPi = newRiemannDpDxi; - gradPj = newRiemannDpDxj; - gradVi = newRiemannDvDxi; - gradVj = newRiemannDvDxj; + if (gradType==GradientType::SPHSameTimeGradient or + gradType==GradientType::SPHUncorrectedGradient){ + gradPi = newRiemannDpDx(nodeListi,i); + gradPj = newRiemannDpDx(nodeListj,j); + gradVi = newRiemannDvDx(nodeListi,i); + gradVj = newRiemannDvDx(nodeListj,j); } - riemannSolver.interfaceState(i, j, - nodeListi, nodeListj, - ri, rj, + riemannSolver.interfaceState(ri, rj, + Hi, Hj, rhoi, rhoj, ci, cj, Peffi, Peffj, vi, vj, + gradRhoi, gradRhoj, gradPi, gradPj, gradVi, gradVj, - Pstar, - vstar, - rhostari, - rhostarj); + Pstar, //output + vstar, //output + rhostari, //output + rhostarj); //output // get our basis function and interface area vectors //-------------------------------------------------------- @@ -265,16 +272,16 @@ evaluateDerivatives(const typename Dimension::Scalar time, // acceleration //------------------------------------------------------ const auto deltaDvDt = Pstar*Astar; - DvDti -= deltaDvDt; - DvDtj += deltaDvDt; + DvDti -= deltaDvDt/mi; + DvDtj += deltaDvDt/mj; // energy //------------------------------------------------------ const auto deltaDepsDti = Pstar*Astar.dot(vi-vstar); const auto deltaDepsDtj = Pstar*Astar.dot(vstar-vj); - DepsDti += deltaDepsDti; - DepsDtj += deltaDepsDtj; + DepsDti += deltaDepsDti/mi; + DepsDtj += deltaDepsDtj/mj; if(compatibleEnergy){ const auto invmij = 1.0/(mi*mj); @@ -345,7 +352,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, const auto hmax = nodeList.hmax(); const auto hminratio = nodeList.hminratio(); const auto nPerh = nodeList.nodesPerSmoothingScale(); - + //const auto kernelExtent = nodeList.neighbor().kernelExtent(); const auto ni = nodeList.numInternalNodes(); #pragma omp parallel for for (auto i = 0u; i < ni; ++i) { @@ -373,9 +380,6 @@ evaluateDerivatives(const typename Dimension::Scalar time, auto& weightedNeighborSumi = weightedNeighborSum(nodeListi, i); auto& massSecondMomenti = massSecondMoment(nodeListi, i); - DvDti /= mi; - DepsDti /= mi; - normi += voli*Hdeti*W0; DvolDti = voli * DvDxi.Trace() ; @@ -401,6 +405,7 @@ evaluateDerivatives(const typename Dimension::Scalar time, hmax, hminratio, nPerh); + Hideali = smoothingScale.newSmoothingScale(Hi, ri, weightedNeighborSumi, @@ -431,9 +436,11 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, const State& state, StateDerivatives& derivatives) const { + const auto calcSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient + or this->gradientType() == GradientType::SPHUncorrectedGradient); + const auto correctSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient); // The kernels and such. const auto& W = this->kernel(); - const auto gradType = this->gradientType(); // The connectivity. const auto& connectivityMap = dataBase.connectivityMap(); @@ -441,11 +448,14 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, const auto numNodeLists = nodeLists.size(); // Get the state and derivative FieldLists. + const auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); const auto volume = state.fields(HydroFieldNames::volume, 0.0); const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); const auto pressure = state.fields(HydroFieldNames::pressure, 0.0); const auto position = state.fields(HydroFieldNames::position, Vector::zero); const auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + + CHECK(massDensity.size() == numNodeLists); CHECK(volume.size() == numNodeLists); CHECK(velocity.size() == numNodeLists); CHECK(pressure.size() == numNodeLists); @@ -453,10 +463,12 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(H.size() == numNodeLists); auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); CHECK(M.size() == numNodeLists); + CHECK(DrhoDx.size() == numNodeLists); CHECK(newRiemannDpDx.size() == numNodeLists); CHECK(newRiemannDvDx.size() == numNodeLists); @@ -471,6 +483,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, typename SpheralThreads::FieldListStack threadStack; auto M_thread = M.threadCopy(threadStack); + auto DrhoDx_thread = DrhoDx.threadCopy(threadStack); auto newRiemannDpDx_thread = newRiemannDpDx.threadCopy(threadStack); auto newRiemannDvDx_thread = newRiemannDvDx.threadCopy(threadStack); @@ -482,6 +495,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, nodeListj = pairs[kk].j_list; // Get the state for node i. + const auto& rhoi = massDensity(nodeListi, i); const auto& ri = position(nodeListi, i); const auto& voli = volume(nodeListi, i); const auto& Hi = H(nodeListi, i); @@ -489,9 +503,11 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(voli > 0.0); CHECK(Hdeti > 0.0); + auto& DrhoDxi = DrhoDx_thread(nodeListi, i); auto& Mi = M_thread(nodeListi, i); // Get the state for node j + const auto& rhoj = massDensity(nodeListj, j); const auto& rj = position(nodeListj, j); const auto& volj = volume(nodeListj, j); const auto& Hj = H(nodeListj, j); @@ -499,6 +515,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, CHECK(volj > 0.0); CHECK(Hdetj > 0.0); + auto& DrhoDxj = DrhoDx_thread(nodeListj, j); auto& Mj = M_thread(nodeListj, j); const auto rij = ri - rj; @@ -524,12 +541,17 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, Mi -= rij.dyad(gradPsii); Mj -= rij.dyad(gradPsij); + DrhoDxi -= (rhoi - rhoj) * gradPsii; + DrhoDxj -= (rhoi - rhoj) * gradPsij; + // // based on nodal values - if (gradType == GradientType::SPHSameTimeGradient){ + if (calcSpatialGradients){ + const auto& vi = velocity(nodeListi, i); const auto& Pi = pressure(nodeListi, i); const auto& vj = velocity(nodeListj, j); const auto& Pj = pressure(nodeListj, j); + auto& newRiemannDpDxi = newRiemannDpDx_thread(nodeListi, i); auto& newRiemannDvDxi = newRiemannDvDx_thread(nodeListi, i); auto& newRiemannDpDxj = newRiemannDpDx_thread(nodeListj, j); @@ -540,6 +562,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, newRiemannDvDxi -= (vi-vj).dyad(gradPsii); newRiemannDvDxj -= (vi-vj).dyad(gradPsij); + } } // loop over pairs @@ -564,12 +587,17 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, Mi = ( goodM ? Mi.Inverse() : Tensor::one); - if (gradType == GradientType::SPHSameTimeGradient){ + auto& DrhoDxi = DrhoDx(nodeListi, i); + DrhoDxi = Mi.Transpose()*DrhoDxi; + + if (correctSpatialGradients){ + auto& newRiemannDpDxi = newRiemannDpDx(nodeListi, i); auto& newRiemannDvDxi = newRiemannDvDx(nodeListi, i); newRiemannDpDxi = Mi.Transpose()*newRiemannDpDxi; newRiemannDvDxi = newRiemannDvDxi*Mi; + } } @@ -579,10 +607,11 @@ computeMCorrection(const typename Dimension::Scalar /*time*/, boundItr != this->boundaryEnd(); ++boundItr)(*boundItr)->applyFieldListGhostBoundary(M); - if (gradType == GradientType::SPHSameTimeGradient){ + if (calcSpatialGradients){ for (ConstBoundaryIterator boundItr = this->boundaryBegin(); boundItr != this->boundaryEnd(); ++boundItr){ + (*boundItr)->applyFieldListGhostBoundary(DrhoDx); (*boundItr)->applyFieldListGhostBoundary(newRiemannDpDx); (*boundItr)->applyFieldListGhostBoundary(newRiemannDvDx); } diff --git a/src/GSPH/MFMHydroBase.cc b/src/GSPH/MFMHydroBase.cc index ea00b8748..6b540aec3 100644 --- a/src/GSPH/MFMHydroBase.cc +++ b/src/GSPH/MFMHydroBase.cc @@ -87,7 +87,6 @@ MFMHydroBase(const SmoothingScaleBase& smoothingScaleMethod, mDvolumeDt = dataBase.newFluidFieldList(0.0, IncrementState::prefix() + HydroFieldNames::volume); - } //------------------------------------------------------------------------------ @@ -198,7 +197,7 @@ initialize(const typename Dimension::Scalar time, const DataBase& dataBase, State& state, StateDerivatives& derivs) { - GenericRiemannHydro::initialize(time,dt,dataBase,state,derivs); + GenericRiemannHydro::initialize(time,dt,dataBase,state,derivs); } //------------------------------------------------------------------------------ diff --git a/src/GSPH/MFVEvaluateDerivatives.cc b/src/GSPH/MFVEvaluateDerivatives.cc new file mode 100644 index 000000000..f2a5efc6a --- /dev/null +++ b/src/GSPH/MFVEvaluateDerivatives.cc @@ -0,0 +1,770 @@ +namespace Spheral { + +template +void +MFVHydroBase:: +evaluateDerivatives(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const { + this->firstDerivativesLoop(time,dt,dataBase,state,derivatives); + this->secondDerivativesLoop(time,dt,dataBase,state,derivatives); + //this->setH(time,dt,dataBase,state,derivatves) +} +//------------------------------------------------------------------------------ +// Determine the principle derivatives. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +secondDerivativesLoop(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const { + + const auto& riemannSolver = this->riemannSolver(); + const auto& smoothingScale = this->smoothingScaleMethod(); + + // A few useful constants we'll use in the following loop. + const auto tiny = std::numeric_limits::epsilon(); + const auto epsTensile = this->epsilonTensile(); + const auto compatibleEnergy = this->compatibleEnergyEvolution(); + //const auto totalEnergy = this->evolveTotalEnergy(); + const auto gradType = this->gradientType(); + //const auto correctVelocityGradient = this->correctVelocityGradient(); + + // The connectivity. + const auto& connectivityMap = dataBase.connectivityMap(); + const auto& nodeLists = connectivityMap.nodeLists(); + const auto& pairs = connectivityMap.nodePairList(); + const auto npairs = pairs.size(); + const auto numNodeLists = nodeLists.size(); + const auto nPerh = nodeLists[0]->nodesPerSmoothingScale(); + + // kernel + const auto& W = this->kernel(); + const auto WnPerh = W(1.0/nPerh, 1.0); + //const auto W0 = W(0.0, 1.0); + + // Get the state and derivative FieldLists. + // State FieldLists. + const auto mass = state.fields(HydroFieldNames::mass, 0.0); + const auto position = state.fields(HydroFieldNames::position, Vector::zero); + const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); + //const auto nodalVelocity = state.fields(GSPHFieldNames::nodalVelocity, Vector::zero); + const auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); + const auto volume = state.fields(HydroFieldNames::volume, 0.0); + const auto specificThermalEnergy = state.fields(HydroFieldNames::specificThermalEnergy, 0.0); + const auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + const auto pressure = state.fields(HydroFieldNames::pressure, 0.0); + const auto soundSpeed = state.fields(HydroFieldNames::soundSpeed, 0.0); + const auto riemannDpDx = state.fields(GSPHFieldNames::RiemannPressureGradient,Vector::zero); + const auto riemannDvDx = state.fields(GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); + + //CHECK(nodalVelocity.size() == numNodeLists); + CHECK(mass.size() == numNodeLists); + CHECK(position.size() == numNodeLists); + CHECK(velocity.size() == numNodeLists); + CHECK(massDensity.size() == numNodeLists); + CHECK(volume.size() == numNodeLists); + CHECK(specificThermalEnergy.size() == numNodeLists); + CHECK(H.size() == numNodeLists); + CHECK(pressure.size() == numNodeLists); + CHECK(soundSpeed.size() == numNodeLists); + CHECK(riemannDpDx.size() == numNodeLists); + CHECK(riemannDvDx.size() == numNodeLists); + + // Derivative FieldLists. + const auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + const auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); + auto normalization = derivatives.fields(HydroFieldNames::normalization, 0.0); + auto DxDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::position, Vector::zero); + auto DvolDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::volume, 0.0); + auto DmDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::mass, 0.0); + auto DEDt = derivatives.fields(IncrementState::prefix() + GSPHFieldNames::thermalEnergy, 0.0); + auto DpDt = derivatives.fields(IncrementState::prefix() + GSPHFieldNames::momentum, Vector::zero); + auto DvDx = derivatives.fields(HydroFieldNames::velocityGradient, Tensor::zero); + auto DHDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::H, SymTensor::zero); + auto Hideal = derivatives.fields(ReplaceBoundedState::prefix() + HydroFieldNames::H, SymTensor::zero); + auto XSPHDeltaV = derivatives.fields(HydroFieldNames::XSPHDeltaV, Vector::zero); + auto weightedNeighborSum = derivatives.fields(HydroFieldNames::weightedNeighborSum, 0.0); + auto massSecondMoment = derivatives.fields(HydroFieldNames::massSecondMoment, SymTensor::zero); + //auto HStretchTensor = derivatives.fields("HStretchTensor", SymTensor::zero); + auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); + auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); + auto& pairAccelerations = derivatives.getAny(HydroFieldNames::pairAccelerations, vector()); + auto& pairDepsDt = derivatives.getAny(HydroFieldNames::pairWork, vector()); + auto& pairMassFlux = derivatives.getAny(GSPHFieldNames::pairMassFlux, vector()); + + CHECK(DrhoDx.size() == numNodeLists); + CHECK(M.size() == numNodeLists); + CHECK(normalization.size() == numNodeLists); + CHECK(DxDt.size() == numNodeLists); + CHECK(DvolDt.size() == numNodeLists); + CHECK(DEDt.size() == numNodeLists); + CHECK(DpDt.size() == numNodeLists); + CHECK(DvDx.size() == numNodeLists); + CHECK(DHDt.size() == numNodeLists); + CHECK(Hideal.size() == numNodeLists); + //CHECK(XSPHDeltaV.size() == numNodeLists); + CHECK(weightedNeighborSum.size() == numNodeLists); + CHECK(massSecondMoment.size() == numNodeLists); + //CHECK(HStretchTensor.size() == numNodeLists); + CHECK(newRiemannDpDx.size() == numNodeLists); + CHECK(newRiemannDvDx.size() == numNodeLists); + + if (compatibleEnergy){ + pairAccelerations.resize(npairs); + pairDepsDt.resize(2*npairs); + pairMassFlux.resize(npairs); + } + + // Walk all the interacting pairs. +#pragma omp parallel + { + // Thread private scratch variables + int i, j, nodeListi, nodeListj; + Scalar Wi, Wj, gWi, gWj, Pstar, rhostari, rhostarj; + Vector vstar; + + typename SpheralThreads::FieldListStack threadStack; + //auto weightedNeighborSum_thread = weightedNeighborSum.threadCopy(threadStack); + //auto massSecondMoment_thread = massSecondMoment.threadCopy(threadStack); + auto DvolDt_thread = DvolDt.threadCopy(threadStack); + auto DmDt_thread = DmDt.threadCopy(threadStack); + auto DEDt_thread = DEDt.threadCopy(threadStack); + auto DpDt_thread = DpDt.threadCopy(threadStack); + auto DvDx_thread = DvDx.threadCopy(threadStack); + auto newRiemannDpDx_thread = newRiemannDpDx.threadCopy(threadStack); + auto newRiemannDvDx_thread = newRiemannDvDx.threadCopy(threadStack); + auto XSPHDeltaV_thread = XSPHDeltaV.threadCopy(threadStack); + //auto normalization_thread = normalization.threadCopy(threadStack); + //auto HStretchTensor_thread = HStretchTensor.threadCopy(threadStack); + + // this is kind of criminal and should be fixed, but for testing purposes + // I'm going to say its allowable. We're going to zero out the thread + // copy of the Hstretch Tensor so that we can zero it out then replace + // the original with the smoothed version. + // HStretchTensor_thread.Zero(); + +#pragma omp for + for (auto kk = 0u; kk < npairs; ++kk) { + i = pairs[kk].i_node; + j = pairs[kk].j_node; + nodeListi = pairs[kk].i_list; + nodeListj = pairs[kk].j_list; + + const auto& mi = mass(nodeListi, i); + const auto& mj = mass(nodeListj, j); + + if( mi >tiny or mj > tiny){ + + // Get the state for node i. + //const auto& ui = nodalVelocity(nodeListi,i); + const auto& riemannDpDxi = riemannDpDx(nodeListi, i); + const auto& riemannDvDxi = riemannDvDx(nodeListi, i); + const auto& ri = position(nodeListi, i); + const auto& vi = velocity(nodeListi, i); + const auto& rhoi = massDensity(nodeListi, i); + const auto& voli = volume(nodeListi, i); + const auto& epsi = specificThermalEnergy(nodeListi, i); + const auto& Pi = pressure(nodeListi, i); + const auto& Hi = H(nodeListi, i); + const auto& ci = soundSpeed(nodeListi, i); + const auto Hdeti = Hi.Determinant(); + CHECK(voli > 0.0); + CHECK(rhoi > 0.0); + CHECK(Hdeti > 0.0); + + //auto& normi = normalization_thread(nodeListi,i); + auto& DvolDti = DvolDt_thread(nodeListi,i); + auto& DmDti = DmDt_thread(nodeListi, i); + auto& DEDti = DEDt_thread(nodeListi, i); + auto& DpDti = DpDt_thread(nodeListi, i); + const auto& DxDti = DxDt(nodeListi,i); + auto& newRiemannDpDxi = newRiemannDpDx_thread(nodeListi, i); + auto& newRiemannDvDxi = newRiemannDvDx_thread(nodeListi, i); + auto& DvDxi = DvDx_thread(nodeListi, i); + //auto& weightedNeighborSumi = weightedNeighborSum_thread(nodeListi, i); + //auto& massSecondMomenti = massSecondMoment(nodeListi, i); + auto& XSPHDeltaVi = XSPHDeltaV_thread(nodeListi,i); + const auto& gradRhoi = DrhoDx(nodeListi, i); + const auto& Mi = M(nodeListi,i); + + + // Get the state for node j + //const auto& uj = nodalVelocity(nodeListj,j); + const auto& riemannDpDxj = riemannDpDx(nodeListj, j); + const auto& riemannDvDxj = riemannDvDx(nodeListj, j); + const auto& rj = position(nodeListj, j); + const auto& vj = velocity(nodeListj, j); + const auto& rhoj = massDensity(nodeListj, j); + const auto& volj = volume(nodeListj, j); + const auto& epsj = specificThermalEnergy(nodeListj, j); + const auto& Pj = pressure(nodeListj, j); + const auto& Hj = H(nodeListj, j); + const auto& cj = soundSpeed(nodeListj, j); + const auto Hdetj = Hj.Determinant(); + CHECK(rhoj > 0.0); + CHECK(volj > 0.0); + CHECK(Hdetj > 0.0); + + //auto& normj = normalization_thread(nodeListj,j); + auto& DvolDtj = DvolDt_thread(nodeListj,j); + auto& DmDtj = DmDt_thread(nodeListj, j); + auto& DEDtj = DEDt_thread(nodeListj, j); + auto& DpDtj = DpDt_thread(nodeListj, j); + const auto& DxDtj = DxDt(nodeListj,j); + auto& newRiemannDpDxj = newRiemannDpDx_thread(nodeListj,j); + auto& newRiemannDvDxj = newRiemannDvDx_thread(nodeListj,j); + auto& DvDxj = DvDx_thread(nodeListj, j); + //auto& weightedNeighborSumj = weightedNeighborSum_thread(nodeListj, j); + //auto& massSecondMomentj = massSecondMoment(nodeListj, j); + auto& XSPHDeltaVj = XSPHDeltaV_thread(nodeListj,j); + const auto& gradRhoj = DrhoDx(nodeListj, j); + const auto& Mj = M(nodeListj,j); + + // Node displacement. + const auto rij = ri - rj; + //const auto rMagij = rij.magnitude(); + //const auto rMagij2 = rij.magnitude2(); + const auto rhatij =rij.unitVector(); + const auto vij = vi - vj; + const auto etai = Hi*rij; + const auto etaj = Hj*rij; + const auto etaMagi = etai.magnitude(); + const auto etaMagj = etaj.magnitude(); + CHECK(etaMagi >= 0.0); + CHECK(etaMagj >= 0.0); + + + // Symmetrized kernel weight and gradient. + //------------------------------------------------------ + W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi); + const auto Hetai = Hi*etai.unitVector(); + const auto gradWi = gWi*Hetai; + + W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj); + const auto Hetaj = Hj*etaj.unitVector(); + const auto gradWj = gWj*Hetaj; + + const auto gradPsii = voli * Mi.Transpose()*gradWi; + const auto gradPsij = volj * Mj.Transpose()*gradWj; + + const auto Astar = voli*gradPsii + volj*gradPsij; + + // Determine an effective pressure including a term to fight the tensile instability. + //const auto fij = epsTensile*pow(Wi/(Hdeti*WnPerh), nTensile); + const auto fij = epsTensile*FastMath::pow4(Wi/(Hdeti*WnPerh)); + const auto Ri = fij*(Pi < 0.0 ? -Pi : 0.0); + const auto Rj = fij*(Pj < 0.0 ? -Pj : 0.0); + const auto Peffi = Pi + Ri; + const auto Peffj = Pj + Rj; + + // Reimann Solver and Fluxes + //------------------------------------------------------ + // we'll clean this up when we have a gradient + // implementation we're in love with + auto gradPi = riemannDpDxi; + auto gradPj = riemannDpDxj; + auto gradVi = riemannDvDxi; + auto gradVj = riemannDvDxj; + if (gradType==GradientType::SPHSameTimeGradient or + gradType==GradientType::SPHUncorrectedGradient){ + gradPi = newRiemannDpDx(nodeListi,i); + gradPj = newRiemannDpDx(nodeListj,j); + gradVi = newRiemannDvDx(nodeListi,i); + gradVj = newRiemannDvDx(nodeListj,j); + } + // need grad rho and grad eps + riemannSolver.interfaceState(ri, rj, + Hi, Hj, + rhoi, rhoj, + ci, cj, + Peffi, Peffj, + vi, vj, + gradRhoi, gradRhoj, + gradPi, gradPj, + gradVi, gradVj, + Pstar, //output + vstar, //output + rhostari, //output + rhostarj); //output + + const auto fluxSwitch = 1.0;//(nodeListi==nodeListj ? 1.0 : 0.0); + const auto vframe = (DxDti+DxDtj)*0.5; + const auto vflux = vstar-vframe; + const auto fluxTowardsNodei = vflux.dot(rhatij) > 0; + const auto rhostar = (fluxTowardsNodei ? rhostarj : rhostari); // we'll need to fix these later + const auto epsstar = (fluxTowardsNodei ? epsj : epsi); // we'll need to fix these later + + const auto massFlux = fluxSwitch * rhostar * vflux.dot(Astar); + const auto momentumFlux = massFlux * vstar; + const auto energyFlux = massFlux * epsstar; + + // mass + //------------------------------------------------------ + DmDti -= massFlux; + DmDtj += massFlux; + + // momentum + //------------------------------------------------------ + const auto deltaDvDt = Pstar*Astar + momentumFlux; + DpDti -= deltaDvDt; + DpDtj += deltaDvDt; + + // energy + //------------------------------------------------------ + const auto deltaDepsDti = Pstar*Astar.dot(vi-vstar) - energyFlux; + const auto deltaDepsDtj = Pstar*Astar.dot(vstar-vj) + energyFlux; + + DEDti += deltaDepsDti; + DEDtj += deltaDepsDtj; + + if(compatibleEnergy){ + pairMassFlux[kk] = massFlux; + pairAccelerations[kk] = deltaDvDt; + pairDepsDt[2*kk] = deltaDepsDti; + pairDepsDt[2*kk+1] = deltaDepsDtj; + } + + // volume change based on nodal velocity + //----------------------------------------------------- + DvolDti -= (DxDti-DxDtj).dot(gradPsii); + DvolDtj -= (DxDti-DxDtj).dot(gradPsij); + + // gradients + //------------------------------------------------------ + const auto deltaDvDxi = 2.0*(vi-vstar).dyad(gradPsii); + const auto deltaDvDxj = 2.0*(vstar-vj).dyad(gradPsij); + + XSPHDeltaVi -= voli*gradWi; + XSPHDeltaVj += volj*gradWj; + + // based on riemann soln + DvDxi -= deltaDvDxi; + DvDxj -= deltaDvDxj; + + // while we figure out what we want ... + switch(gradType){ + case GradientType::RiemannGradient: // default grad based on riemann soln + newRiemannDvDxi -= deltaDvDxi; + newRiemannDvDxj -= deltaDvDxj; + newRiemannDpDxi -= 2.0*(Pi-Pstar)*gradPsii; + newRiemannDpDxj -= 2.0*(Pstar-Pj)*gradPsij; + break; + case GradientType::HydroAccelerationGradient: // based on hydro accel for DpDx + newRiemannDvDxi -= deltaDvDxi; + newRiemannDvDxj -= deltaDvDxj; + newRiemannDpDxi += Pstar*Astar/voli; + newRiemannDpDxj -= Pstar*Astar/volj; + break; + case GradientType::SPHGradient: // raw gradients + newRiemannDvDxi -= (vi-vj).dyad(gradPsii); + newRiemannDvDxj -= (vi-vj).dyad(gradPsij); + newRiemannDpDxi -= (Pi-Pj)*gradPsii; + newRiemannDpDxj -= (Pi-Pj)*gradPsij; + break; + case GradientType::MixedMethodGradient: // raw gradient for P riemann gradient for v + newRiemannDvDxi -= deltaDvDxi; + newRiemannDvDxj -= deltaDvDxj; + newRiemannDpDxi -= (Pi-Pj)*gradPsii; + newRiemannDpDxj -= (Pi-Pj)*gradPsij; + break; + default: + break; + } + + } //if statement + } // loop over pairs + threadReduceFieldLists(threadStack); + } // OpenMP parallel region + + + // Finish up the derivatives for each point. + for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { + const auto& nodeList = mass[nodeListi]->nodeList(); + const auto hmin = nodeList.hmin(); + const auto hmax = nodeList.hmax(); + const auto hminratio = nodeList.hminratio(); + const auto nPerh = nodeList.nodesPerSmoothingScale(); + //const auto kernelExtent = nodeList.neighbor().kernelExtent(); + const auto ni = nodeList.numInternalNodes(); + +#pragma omp parallel for + for (auto i = 0u; i < ni; ++i) { + + // Get the state for node i. + const auto& ri = position(nodeListi, i); + const auto& voli = volume(nodeListi,i); + //const auto& ui = nodalVelocity(nodeListi,i); + //const auto& vi = velocity(nodeListi,i); + //const auto& ci = soundSpeed(nodeListi,i); + const auto& Hi = H(nodeListi, i); + const auto Hdeti = Hi.Determinant(); + CHECK(voli > 0.0); + CHECK(Hdeti > 0.0); + + //auto& normi = normalization(nodeListi, i); + //auto& DxDti = DxDt(nodeListi, i); + auto& DvolDti = DvolDt(nodeListi, i); + auto& DvDxi = DvDx(nodeListi, i); + auto& DHDti = DHDt(nodeListi, i); + auto& Hideali = Hideal(nodeListi, i); + auto& XSPHDeltaVi = XSPHDeltaV(nodeListi, i); + const auto& weightedNeighborSumi = weightedNeighborSum(nodeListi, i); + const auto& massSecondMomenti = massSecondMoment(nodeListi, i); + //const auto& HStretchTensori = HStretchTensor(nodeListi, i); + + XSPHDeltaVi /= Dimension::rootnu(Hdeti); + DvolDti *= voli; + + // If needed finish the total energy derivative. + //if (totalEnergy) DepsDti = mi*(vi.dot(DvDti) + DepsDti); + + // ----------------------------------------------- + // TODO: + // this makes ui be vi from the previous timestep. We might need a special update method for hthis + // We culd also just take care of these in the primary loop and make the node velocity a deriv + // ----------------------------------------------- + //if(true){ + DHDti = smoothingScale.smoothingScaleDerivative(Hi, + ri, + DvDxi, + hmin, + hmax, + hminratio, + nPerh); + Hideali = smoothingScale.newSmoothingScale(Hi, // Hi + ri, // ri + weightedNeighborSumi, // ? + massSecondMomenti, // Hstretch tensor + W, // W + hmin, // hmin + hmax, // hmax + hminratio, // hminratio + nPerh, // Ngb + connectivityMap, // connectivityMap + nodeListi, // nodeListi + i); // i + // }else{ + // // smoothing scale construction + // const auto Ngb_target = (Dimension::nDim == 3 ? 32 : + // (Dimension::nDim == 2 ? 16 : + // 4)); + // const auto stretchFactor = 0.00; + + // // set on construction + // const auto C = (Dimension::nDim == 3 ? 1.33333*3.1415 : + // (Dimension::nDim == 2 ? 3.1415 : + // 1.0)); + + // // pass + // const auto Ngb = C /(Hdeti*voli) * pow(kernelExtent,Dimension::nDim); + + // const auto Hstretch = ((1.00-stretchFactor)* SymTensor::one + + // stretchFactor * HStretchTensori)*Hi; + + // const auto scaleFactor = (1.0+0.5*(Ngb - Ngb_target)/Ngb_target); + // Hideali = std::min(std::max(scaleFactor,0.8),1.2) * Hstretch; + + // DHDti = 0.25*(Hideali-Hi)/dt; + // } + } // nodes loop + } // nodeLists loop +} // eval derivs method + +//------------------------------------------------------------------------------ +// EvalDerivs subroutine for spatial derivs +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +firstDerivativesLoop(const typename Dimension::Scalar /*time*/, + const typename Dimension::Scalar /*dt*/, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const { + + const auto tiny = std::numeric_limits::epsilon(); + //const auto epsTensile = this->epsilonTensile(); + const auto nodeMotionCoeff = this->nodeMotionCoefficient(); + + const auto calcSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient + or this->gradientType() == GradientType::SPHUncorrectedGradient); + const auto correctSpatialGradients = (this->gradientType() == GradientType::SPHSameTimeGradient); + + const auto nodeMotion = this->nodeMotionType(); + const auto xsphMotion = (nodeMotion == NodeMotionType::XSPH); + const auto ficianMotion = (nodeMotion == NodeMotionType::Fician); + const auto noMotion = (nodeMotion == NodeMotionType::Eulerian); + + // The connectivity. + const auto& connectivityMap = dataBase.connectivityMap(); + const auto& nodeLists = connectivityMap.nodeLists(); + const auto numNodeLists = nodeLists.size(); + const auto& pairs = connectivityMap.nodePairList(); + const auto npairs = pairs.size(); + //const auto nPerh = nodeLists[0]->nodesPerSmoothingScale(); + + // kernel + const auto& W = this->kernel(); + //const auto WnPerh = W(1.0/nPerh, 1.0); + const auto W0 = W(0.0, 1.0); + + // Get the state and derivative FieldLists. + const auto soundSpeed = state.fields(HydroFieldNames::soundSpeed, 0.0); + const auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); + const auto volume = state.fields(HydroFieldNames::volume, 0.0); + const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); + const auto pressure = state.fields(HydroFieldNames::pressure, 0.0); + const auto position = state.fields(HydroFieldNames::position, Vector::zero); + const auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + + CHECK(soundSpeed.size() == numNodeLists); + CHECK(massDensity.size() == numNodeLists); + CHECK(volume.size() == numNodeLists); + CHECK(velocity.size() == numNodeLists); + CHECK(pressure.size() == numNodeLists); + CHECK(position.size() == numNodeLists); + CHECK(H.size() == numNodeLists); + + auto M = derivatives.fields(HydroFieldNames::M_SPHCorrection, Tensor::zero); + auto DxDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::position, Vector::zero); + auto DrhoDx = derivatives.fields(GSPHFieldNames::densityGradient, Vector::zero); + auto newRiemannDpDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero); + auto newRiemannDvDx = derivatives.fields(ReplaceState::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero); + auto massSecondMoment = derivatives.fields(HydroFieldNames::massSecondMoment, SymTensor::zero); + auto weightedNeighborSum = derivatives.fields(HydroFieldNames::weightedNeighborSum, 0.0); + //auto HStretchTensor = derivatives.fields("HStretchTensor", SymTensor::zero); + auto normalization = derivatives.fields(HydroFieldNames::normalization, 0.0); + + CHECK(M.size() == numNodeLists); + CHECK(DrhoDx.size() == numNodeLists); + CHECK(DxDt.size() == numNodeLists); + CHECK(newRiemannDpDx.size() == numNodeLists); + CHECK(newRiemannDvDx.size() == numNodeLists); + CHECK(massSecondMoment.size() == numNodeLists) + CHECK(weightedNeighborSum.size() == numNodeLists) + CHECK(normalization.size() == numNodeLists) + //CHECK(HStretchTensor.size() == numNodeLists) + +#pragma omp parallel + { + // Thread private scratch variables + int i, j, nodeListi, nodeListj; + Scalar Wi, Wj, gWi, gWj; + + typename SpheralThreads::FieldListStack threadStack; + auto M_thread = M.threadCopy(threadStack); + auto DrhoDx_thread = DrhoDx.threadCopy(threadStack); + auto newRiemannDpDx_thread = newRiemannDpDx.threadCopy(threadStack); + auto newRiemannDvDx_thread = newRiemannDvDx.threadCopy(threadStack); + auto DxDt_thread = DxDt.threadCopy(threadStack); + auto massSecondMoment_thread = massSecondMoment.threadCopy(threadStack); + auto weightedNeighborSum_thread = weightedNeighborSum.threadCopy(threadStack); + //auto HStretchTensor_thread = HStretchTensor.threadCopy(threadStack); + auto normalization_thread = normalization.threadCopy(threadStack); + +#pragma omp for + for (auto kk = 0u; kk < npairs; ++kk) { + i = pairs[kk].i_node; + j = pairs[kk].j_node; + nodeListi = pairs[kk].i_list; + nodeListj = pairs[kk].j_list; + + // Get the state for node i. + const auto& vi = velocity(nodeListi, i); + const auto& Pi = pressure(nodeListi, i); + const auto& ci = soundSpeed(nodeListi, i); + const auto& rhoi = massDensity(nodeListi, i); + const auto& ri = position(nodeListi, i); + const auto& voli = volume(nodeListi, i); + const auto& Hi = H(nodeListi, i); + const auto Hdeti = Hi.Determinant(); + CHECK(voli > 0.0); + CHECK(Hdeti > 0.0); + + auto& DxDti = DxDt_thread(nodeListi,i); + //auto& HStretchTensori = HStretchTensor_thread(nodeListi,i); + auto& weightedNeighborSumi = weightedNeighborSum_thread(nodeListi,i); + auto& massSecondMomenti = massSecondMoment_thread(nodeListi, i); + auto& normi = normalization(nodeListi,i); + auto& DrhoDxi = DrhoDx_thread(nodeListi, i); + auto& newRiemannDpDxi = newRiemannDpDx_thread(nodeListi, i); + auto& newRiemannDvDxi = newRiemannDvDx_thread(nodeListi, i); + auto& Mi = M_thread(nodeListi, i); + + // Get the state for node j + const auto& vj = velocity(nodeListj, j); + const auto& Pj = pressure(nodeListj, j); + const auto& cj = soundSpeed(nodeListj, j); + const auto& rhoj = massDensity(nodeListj, j); + const auto& rj = position(nodeListj, j); + const auto& volj = volume(nodeListj, j); + const auto& Hj = H(nodeListj, j); + const auto Hdetj = Hj.Determinant(); + CHECK(volj > 0.0); + CHECK(Hdetj > 0.0); + + auto& DxDtj = DxDt_thread(nodeListj,j); + //auto& HStretchTensorj = HStretchTensor_thread(nodeListj,j); + auto& weightedNeighborSumj = weightedNeighborSum_thread(nodeListj,j); + auto& massSecondMomentj = massSecondMoment_thread(nodeListj, j); + auto& normj = normalization(nodeListj,j); + auto& DrhoDxj = DrhoDx_thread(nodeListj, j); + auto& newRiemannDpDxj = newRiemannDpDx_thread(nodeListj, j); + auto& newRiemannDvDxj = newRiemannDvDx_thread(nodeListj, j); + auto& Mj = M_thread(nodeListj, j); + + const auto rij = ri - rj; + //const auto rMagij = safeInv(rij.magnitude()); + + const auto etai = Hi*rij; + const auto etaj = Hj*rij; + const auto etaMagi = etai.magnitude(); + const auto etaMagj = etaj.magnitude(); + CHECK(etaMagi >= 0.0); + CHECK(etaMagj >= 0.0); + + W.kernelAndGradValue(etaMagi, Hdeti, Wi, gWi); + const auto Hetai = Hi*etai.unitVector(); + const auto gradWi = gWi*Hetai; + + W.kernelAndGradValue(etaMagj, Hdetj, Wj, gWj); + const auto Hetaj = Hj*etaj.unitVector(); + const auto gradWj = gWj*Hetaj; + + const auto psii = voli*Wi; + const auto psij = volj*Wj; + const auto gradPsii = voli*gradWi; + const auto gradPsij = volj*gradWj; + + weightedNeighborSumi += std::abs(gWi); + weightedNeighborSumj += std::abs(gWj); + + //HStretchTensori -= voli*rij.selfdyad()*gWi*rMagij; + //HStretchTensorj -= volj*rij.selfdyad()*gWj*rMagij; + + const auto rij2 = rij.magnitude2(); + const auto thpt = rij.selfdyad()*safeInvVar(rij2*rij2*rij2); + massSecondMomenti += gradWi.magnitude2()*thpt; + massSecondMomentj += gradWj.magnitude2()*thpt; + + // gradients + Mi -= rij.dyad(gradPsii); + Mj -= rij.dyad(gradPsij); + + DrhoDxi -= (rhoi - rhoj) * gradPsii; + DrhoDxj -= (rhoi - rhoj) * gradPsij; + + if (calcSpatialGradients){ + newRiemannDpDxi -= (Pi-Pj)*gradPsii; + newRiemannDpDxj -= (Pi-Pj)*gradPsij; + + newRiemannDvDxi -= (vi-vj).dyad(gradPsii); + newRiemannDvDxj -= (vi-vj).dyad(gradPsij); + } + + // node motion relative to fluid + //----------------------------------------------------------- + if (xsphMotion) { + const auto cij = 0.5*(ci+cj); + const auto wij = cij * safeInv(max(10*(vi-vj).magnitude(),cij)); + DxDti -= wij*psii*(vi-vj); + DxDtj -= wij*psij*(vj-vi); + } + if(ficianMotion){ + //const auto fi = FastMath::pow4(Wi/(Hdeti*WnPerh)); + //const auto fj = FastMath::pow4(Wj/(Hdetj*WnPerh)); + DxDti -= -rij*psii; + DxDtj += -rij*psij; + } + + normi += psii;//voli*gradWi.magnitude(); + normj += psij;//volj*gradWj.magnitude(); + + } // loop over pairs + + // Reduce the thread values to the master. + threadReduceFieldLists(threadStack); + + } // OpenMP parallel region + + // Finish up the spatial gradient calculation + for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { + const auto& nodeList = M[nodeListi]->nodeList(); + const auto ni = nodeList.numInternalNodes(); + +#pragma omp parallel for + for (auto i = 0u; i < ni; ++i) { + + const auto numNeighborsi = connectivityMap.numNeighborsForNode(nodeListi, i); + + const auto& ci = soundSpeed(nodeListi,i); + const auto& vi = velocity(nodeListi,i); + const auto& voli = volume(nodeListi,i); + const auto& Hi = H(nodeListi, i); + const auto Hdeti = Hi.Determinant(); + + auto& DxDti = DxDt(nodeListi,i); + auto& Mi = M(nodeListi, i); + auto& massSecondMomenti = massSecondMoment(nodeListi,i); + auto& weightedNeighborSumi = weightedNeighborSum(nodeListi,i); + //auto& HStretchTensori = HStretchTensor(nodeListi,i); + auto& normi = normalization(nodeListi, i); + const auto Mdeti = std::abs(Mi.Determinant()); + + normi += voli*Hdeti*W0; + weightedNeighborSumi = Dimension::rootnu(max(0.0, weightedNeighborSumi/Hdeti)); + //HStretchTensori /= Dimension::rootnu(max(HStretchTensori.Determinant(),tiny)); + massSecondMomenti /= Hdeti*Hdeti; + + const auto enoughNeighbors = numNeighborsi > Dimension::pownu(2); + const auto goodM = (Mdeti > 1e-2 and enoughNeighbors); + + Mi = ( goodM ? Mi.Inverse() : Tensor::one); + + if (correctSpatialGradients){ + auto& newRiemannDpDxi = newRiemannDpDx(nodeListi, i); + auto& newRiemannDvDxi = newRiemannDvDx(nodeListi, i); + auto& DrhoDxi = DrhoDx(nodeListi,i); + + DrhoDxi = Mi.Transpose()*DrhoDxi; + newRiemannDpDxi = Mi.Transpose()*newRiemannDpDxi; + newRiemannDvDxi = newRiemannDvDxi*Mi; + } + + if (xsphMotion) DxDti *= nodeMotionCoeff/max(tiny, normi); + if(ficianMotion) DxDti *= nodeMotionCoeff * ci * ci * Dimension::rootnu(Hdeti) * + safeInv( max(10.0*DxDti.magnitude(),ci)); + if(!noMotion) DxDti += vi; + } + + } + + for (ConstBoundaryIterator boundItr = this->boundaryBegin(); + boundItr != this->boundaryEnd(); + ++boundItr){ + (*boundItr)->applyFieldListGhostBoundary(M); + (*boundItr)->applyFieldListGhostBoundary(DxDt); + } + + if (calcSpatialGradients){ + for (ConstBoundaryIterator boundItr = this->boundaryBegin(); + boundItr != this->boundaryEnd(); + ++boundItr){ + (*boundItr)->applyFieldListGhostBoundary(DrhoDx); + (*boundItr)->applyFieldListGhostBoundary(newRiemannDpDx); + (*boundItr)->applyFieldListGhostBoundary(newRiemannDvDx); + } + } + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary(); + +} + +} // spheral namespace \ No newline at end of file diff --git a/src/GSPH/MFVHydroBase.cc b/src/GSPH/MFVHydroBase.cc new file mode 100644 index 000000000..f2da58ce8 --- /dev/null +++ b/src/GSPH/MFVHydroBase.cc @@ -0,0 +1,385 @@ +//---------------------------------Spheral++----------------------------------// +// MFVHydroBase -- This is an Arbitrary Eulerian-Lagrangian extension of the +// MFV approach of Hopkins 2015. Its got several node-motion +// approaches which promote more regular particle distributions. +// +// Each of the ALE options defines the velocity of the nodes +// differently. The flux that results from the difference +// between the nodes velocities and the fluid velocity. +// The velocities are defined as follows for the +// NodeMotionTypes: +// +// 1) Eulerian ---- static Nodes +// 2) Lagrangian -- nodal velocity = fluid velocity. (This is +// a spheralized version of MFV so there +// is some flux between nodes) +// 3) Fician ------ nodal velocity = fluid velocity + Fician +// PST correction +// 4) XSPH -------- nodal velocity = xsph velocity +// +// Hopkins P.F. (2015) "A New Class of Accurate, Mesh-Free Hydrodynamic +// Simulation Methods," MNRAS, 450(1):53-110 +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// +// TODO: +// 1 backpressure and fician particle shifting +// 2 Eulerian model will still crash on the Noh implosion due to void particles +// 3 Good implementation of Ngb update +// 4 treatment for material interfaces +//---------------------------------------------------------------------------// + +#include "FileIO/FileIO.hh" +#include "NodeList/SmoothingScaleBase.hh" +#include "Hydro/HydroFieldNames.hh" + +#include "DataBase/DataBase.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" +#include "DataBase/IncrementState.hh" +#include "DataBase/ReplaceState.hh" +#include "DataBase/PureReplaceState.hh" +#include "DataBase/ReplaceBoundedState.hh" +#include "DataBase/IncrementBoundedState.hh" + +#include "Field/FieldList.hh" +#include "Field/NodeIterators.hh" +#include "Boundary/Boundary.hh" +#include "Neighbor/ConnectivityMap.hh" + +#include "GSPH/MFVHydroBase.hh" +#include "GSPH/GSPHFieldNames.hh" +#include "GSPH/computeSumVolume.hh" +#include "GSPH/computeMFMDensity.hh" +#include "GSPH/Policies/MassFluxPolicy.hh" +#include "GSPH/Policies/ReplaceWithRatioPolicy.hh" +#include "GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.hh" +#include "GSPH/Policies/MFVIncrementVelocityPolicy.hh" +#include "GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh" +#include "GSPH/RiemannSolvers/RiemannSolverBase.hh" + +#ifdef _OPENMP +#include "omp.h" +#endif + +#include + +using std::string; +using std::min; +using std::max; + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor. +//------------------------------------------------------------------------------ +template +MFVHydroBase:: +MFVHydroBase(const SmoothingScaleBase& smoothingScaleMethod, + DataBase& dataBase, + RiemannSolverBase& riemannSolver, + const TableKernel& W, + const Scalar epsDiffusionCoeff, + const double cfl, + const bool useVelocityMagnitudeForDt, + const bool compatibleEnergyEvolution, + const bool evolveTotalEnergy, + const bool XSPH, + const bool correctVelocityGradient, + const double nodeMotionCoefficient, + const NodeMotionType nodeMotionType, + const GradientType gradType, + const MassDensityType densityUpdate, + const HEvolutionType HUpdate, + const double epsTensile, + const double nTensile, + const Vector& xmin, + const Vector& xmax): + GenericRiemannHydro(smoothingScaleMethod, + dataBase, + riemannSolver, + W, + epsDiffusionCoeff, + cfl, + useVelocityMagnitudeForDt, + compatibleEnergyEvolution, + evolveTotalEnergy, + XSPH, + correctVelocityGradient, + gradType, + densityUpdate, + HUpdate, + epsTensile, + nTensile, + xmin, + xmax), + mNodeMotionCoefficient(nodeMotionCoefficient), + mNodeMotionType(nodeMotionType), + mNodalVelocity(FieldStorageType::CopyFields), + mDmassDt(FieldStorageType::CopyFields), + mDthermalEnergyDt(FieldStorageType::CopyFields), + mDmomentumDt(FieldStorageType::CopyFields), + mDvolumeDt(FieldStorageType::CopyFields), + //mHStretchTensor(FieldStorageType::CopyFields), + mPairMassFlux(){ + mNodalVelocity = dataBase.newFluidFieldList(Vector::zero, GSPHFieldNames::nodalVelocity); + mDmassDt = dataBase.newFluidFieldList(0.0, IncrementState::prefix() + HydroFieldNames::mass); + mDthermalEnergyDt = dataBase.newFluidFieldList(0.0, IncrementState::prefix() + GSPHFieldNames::thermalEnergy); + mDmomentumDt = dataBase.newFluidFieldList(Vector::zero, IncrementState::prefix() + GSPHFieldNames::momentum); + mDvolumeDt = dataBase.newFluidFieldList(0.0, IncrementState::prefix() + HydroFieldNames::volume); + //mHStretchTensor = dataBase.newFluidFieldList(SymTensor::zero, "HStretchTensor"); + mPairMassFlux.clear(); +} + +//------------------------------------------------------------------------------ +// Destructor +//------------------------------------------------------------------------------ +template +MFVHydroBase:: +~MFVHydroBase() { +} + +//------------------------------------------------------------------------------ +// On problem start up, we need to initialize our internal data. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +initializeProblemStartup(DataBase& dataBase) { + GenericRiemannHydro::initializeProblemStartup(dataBase); +} + +//------------------------------------------------------------------------------ +// Register the state we need/are going to evolve. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +registerState(DataBase& dataBase, + State& state) { + + GenericRiemannHydro::registerState(dataBase,state); + + dataBase.resizeFluidFieldList(mNodalVelocity, Vector::zero, GSPHFieldNames::nodalVelocity,false); + + auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); + auto position = state.fields(HydroFieldNames::position,Vector::zero); + auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); + auto volume = state.fields(HydroFieldNames::volume, 0.0); + auto mass = state.fields(HydroFieldNames::mass, 0.0); + auto specificThermalEnergy = state.fields(HydroFieldNames::specificThermalEnergy, 0.0); + + // We use the thermal energy to update the specific thermal energy + state.removePolicy(specificThermalEnergy,false); + + CHECK(position.numFields() == dataBase.numFluidNodeLists()); + CHECK(velocity.numFields() == dataBase.numFluidNodeLists()); + CHECK(volume.numFields() == dataBase.numFluidNodeLists()); + CHECK(mass.numFields() == dataBase.numFluidNodeLists()); + CHECK(specificThermalEnergy.numFields() == dataBase.numFluidNodeLists()); + + auto nodeListi = 0u; + for (auto itr = dataBase.fluidNodeListBegin(); + itr < dataBase.fluidNodeListEnd(); + ++itr, ++nodeListi) { + auto& massi = (*itr)->mass(); + auto minVolume = massi.min()/(*itr)->rhoMax(); + auto maxVolume = massi.max()/(*itr)->rhoMin(); + state.enroll(*volume[nodeListi], make_policy>(minVolume, + maxVolume)); + } + + + state.enroll(massDensity, make_policy>({HydroFieldNames::mass, + HydroFieldNames::volume}, + HydroFieldNames::mass, + HydroFieldNames::volume)); + + state.enroll(mass, make_policy>({HydroFieldNames::velocity, + HydroFieldNames::specificThermalEnergy})); + + state.enroll(velocity, + make_policy>({HydroFieldNames::specificThermalEnergy})); + + + if (this->compatibleEnergyEvolution()) { + auto thermalEnergyPolicy = make_policy>(dataBase); + state.enroll(specificThermalEnergy, thermalEnergyPolicy); + }else if (this->evolveTotalEnergy()) { + std::cout <<"evolve total energy not implemented for MFV" << std::endl; + } else { + auto thermalEnergyPolicy = make_policy>(); + state.enroll(specificThermalEnergy,thermalEnergyPolicy); + } + +} + +//------------------------------------------------------------------------------ +// Register the state derivative fields. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +registerDerivatives(DataBase& dataBase, + StateDerivatives& derivs) { + GenericRiemannHydro::registerDerivatives(dataBase,derivs); + dataBase.resizeFluidFieldList(mDmassDt, 0.0, IncrementState::prefix() + HydroFieldNames::mass, false); + dataBase.resizeFluidFieldList(mDthermalEnergyDt, 0.0, IncrementState::prefix() + GSPHFieldNames::thermalEnergy, false); + dataBase.resizeFluidFieldList(mDmomentumDt, Vector::zero, IncrementState::prefix() + GSPHFieldNames::momentum, false); + dataBase.resizeFluidFieldList(mDvolumeDt, 0.0, IncrementState::prefix() + HydroFieldNames::volume, false); + //dataBase.resizeFluidFieldList(mHStretchTensor,SymTensor::zero, "HStretchTensor", false); + derivs.enroll(mDmassDt); + derivs.enroll(mDthermalEnergyDt); + derivs.enroll(mDmomentumDt); + derivs.enroll(mDvolumeDt); + //derivs.enroll(mHStretchTensor); + derivs.enrollAny(GSPHFieldNames::pairMassFlux, mPairMassFlux); +} + +//------------------------------------------------------------------------------ +// This method is called once at the beginning of a timestep, after all state registration. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +preStepInitialize(const DataBase& dataBase, + State& state, + StateDerivatives& derivs) { + GenericRiemannHydro::preStepInitialize(dataBase,state,derivs); + + if(this->densityUpdate() == MassDensityType::RigorousSumDensity){ + + const auto position = state.fields(HydroFieldNames::position, Vector::zero); + const auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + const auto mass = state.fields(HydroFieldNames::mass, 0.0); + auto massDensity = state.fields(HydroFieldNames::massDensity, 0.0); + auto volume = state.fields(HydroFieldNames::volume, 0.0); + + computeSumVolume(dataBase.connectivityMap(),this->kernel(),position,H,volume); + computeMFMDensity(mass,volume,massDensity); + + for (auto boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr){ + (*boundaryItr)->applyFieldListGhostBoundary(volume); + (*boundaryItr)->applyFieldListGhostBoundary(massDensity); + } + for (auto boundaryItr = this->boundaryBegin(); + boundaryItr < this->boundaryEnd(); + ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary(); + } +} + +//------------------------------------------------------------------------------ +// Initialize the hydro before calling evaluateDerivatives +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +initialize(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + State& state, + StateDerivatives& derivs) { + GenericRiemannHydro::initialize(time,dt,dataBase,state,derivs); +} + +//------------------------------------------------------------------------------ +// Finalize the derivatives. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +finalizeDerivatives(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivs) const { + // hackish solution and I should be ashamed. + if (this->compatibleEnergyEvolution()) { + auto DpDt = derivs.fields(IncrementState::prefix() + GSPHFieldNames::momentum, Vector::zero); + auto DmDt = derivs.fields(IncrementState::prefix() + HydroFieldNames::mass, 0.0); + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr){ + (*boundaryItr)->applyFieldListGhostBoundary(DpDt); + (*boundaryItr)->applyFieldListGhostBoundary(DmDt); + } + + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary(); + } +} + +//------------------------------------------------------------------------------ +// Apply the ghost boundary conditions for hydro state fields. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +applyGhostBoundaries(State& state, + StateDerivatives& derivs) { + GenericRiemannHydro::applyGhostBoundaries(state,derivs); + + auto nodalVelocity = state.fields(GSPHFieldNames::nodalVelocity, Vector::zero); + + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr) { + (*boundaryItr)->applyFieldListGhostBoundary(nodalVelocity); + } +} + +//------------------------------------------------------------------------------ +// Enforce the boundary conditions for hydro state fields. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +enforceBoundaries(State& state, + StateDerivatives& derivs) { + GenericRiemannHydro::enforceBoundaries(state,derivs); + + auto nodalVelocity = state.fields(GSPHFieldNames::nodalVelocity, Vector::zero); + + for (ConstBoundaryIterator boundaryItr = this->boundaryBegin(); + boundaryItr != this->boundaryEnd(); + ++boundaryItr) { + (*boundaryItr)->enforceFieldListBoundary(nodalVelocity); + } +} + + +//------------------------------------------------------------------------------ +// Dump the current state to the given file. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +dumpState(FileIO& file, const string& pathName) const { + GenericRiemannHydro::dumpState(file,pathName); + file.write(mNodalVelocity, pathName + "/nodalVelocity"); + file.write(mDmassDt, pathName + "/DmassDt"); + file.write(mDthermalEnergyDt, pathName + "/DthermalEnergyDt"); + file.write(mDmomentumDt, pathName + "/DmomentumDt"); + file.write(mDvolumeDt, pathName + "/DvolumeDt"); +} + +//------------------------------------------------------------------------------ +// Restore the state from the given file. +//------------------------------------------------------------------------------ +template +void +MFVHydroBase:: +restoreState(const FileIO& file, const string& pathName) { + GenericRiemannHydro::restoreState(file,pathName); + file.read(mNodalVelocity, pathName + "/nodalVelocity"); + file.read(mDmassDt, pathName + "/DmassDt"); + file.read(mDthermalEnergyDt, pathName + "/DthermalEnergyDt"); + file.read(mDmomentumDt, pathName + "/DmomentumDt"); + file.read(mDvolumeDt, pathName + "/DvolumeDt"); +} + +} + diff --git a/src/GSPH/MFVHydroBase.hh b/src/GSPH/MFVHydroBase.hh new file mode 100644 index 000000000..593fde4bf --- /dev/null +++ b/src/GSPH/MFVHydroBase.hh @@ -0,0 +1,219 @@ +//---------------------------------Spheral++----------------------------------// +// MFVHydroBase -- This is an Arbitrary Eulerian-Lagrangian extension of the +// MFV approach of Hopkins 2015. Its got several node-motion +// approaches which promote more regular particle distributions. +// +// Each of the ALE options defines the velocity of the nodes +// differently. The flux that results from the difference +// between the nodes velocities and the fluid velocity. +// The velocities are defined as follows for the +// NodeMotionTypes: +// +// 1) Eulerian ---- static Nodes +// 2) Lagrangian -- nodal velocity = fluid velocity. (This is +// a spheralized version of MFV so there +// is some flux between nodes) +// 3) Fician ------ nodal velocity = fluid velocity + Fician +// PST correction +// 4) XSPH -------- nodal velocity = xsph velocity +// +// Hopkins P.F. (2015) "A New Class of Accurate, Mesh-Free Hydrodynamic +// Simulation Methods," MNRAS, 450(1):53-110 +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// +// TODO: +// 1 backpressure and fician particle shifting +// 2 Eulerian model will still crash on the Noh implosion due to void particles +// 3 Good implementation of Ngb update +// 4 treatment for material interfaces +//---------------------------------------------------------------------------// + +#ifndef __Spheral_MFVHydroBase_hh__ +#define __Spheral_MFVHydroBase_hh__ + +#include + +#include "GSPH/GenericRiemannHydro.hh" + +namespace Spheral { + +enum class NodeMotionType { + Lagrangian = 0, + Eulerian = 1, + Fician = 2, + XSPH = 3, + BackgroundPressure = 4, +}; + +template class State; +template class StateDerivatives; +template class SmoothingScaleBase; +template class TableKernel; +template class RiemannSolverBase; +template class DataBase; +template class Field; +template class FieldList; +class FileIO; + +template +class MFVHydroBase: public GenericRiemannHydro { + +public: + //--------------------------- Public Interface ---------------------------// + typedef typename Dimension::Scalar Scalar; + typedef typename Dimension::Vector Vector; + typedef typename Dimension::Tensor Tensor; + typedef typename Dimension::SymTensor SymTensor; + typedef typename Dimension::ThirdRankTensor ThirdRankTensor; + + typedef typename GenericRiemannHydro::TimeStepType TimeStepType; + typedef typename GenericRiemannHydro::ConstBoundaryIterator ConstBoundaryIterator; + + // Constructors. + MFVHydroBase(const SmoothingScaleBase& smoothingScaleMethod, + DataBase& dataBase, + RiemannSolverBase& riemannSolver, + const TableKernel& W, + const Scalar epsDiffusionCoeff, + const double cfl, + const bool useVelocityMagnitudeForDt, + const bool compatibleEnergyEvolution, + const bool evolveTotalEnergy, + const bool XSPH, + const bool correctVelocityGradient, + const double nodeMotionCoefficient, + const NodeMotionType nodeMotionType, + const GradientType gradType, + const MassDensityType densityUpdate, + const HEvolutionType HUpdate, + const double epsTensile, + const double nTensile, + const Vector& xmin, + const Vector& xmax); + + // Destructor. + virtual ~MFVHydroBase(); + + // Tasks we do once on problem startup. + virtual + void initializeProblemStartup(DataBase& dataBase) override; + + // Register the state Hydro expects to use and evolve. + virtual + void registerState(DataBase& dataBase, + State& state) override; + + // Register the derivatives/change fields for updating state. + virtual + void registerDerivatives(DataBase& dataBase, + StateDerivatives& derivs) override; + + + // This method is called once at the beginning of a timestep, after all state registration. + virtual void preStepInitialize(const DataBase& dataBase, + State& state, + StateDerivatives& derivs) override; + + // Initialize the Hydro before we start a derivative evaluation. + virtual + void initialize(const Scalar time, + const Scalar dt, + const DataBase& dataBase, + State& state, + StateDerivatives& derivs) override; + + // Evaluate the derivatives for the principle hydro variables: + // mass density, velocity, and specific thermal energy. + virtual + void evaluateDerivatives(const Scalar time, + const Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const override; + void + firstDerivativesLoop(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const; + + void + secondDerivativesLoop(const typename Dimension::Scalar time, + const typename Dimension::Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivatives) const; + + // Finalize the derivatives. + virtual + void finalizeDerivatives(const Scalar time, + const Scalar dt, + const DataBase& dataBase, + const State& state, + StateDerivatives& derivs) const override; + + // Apply boundary conditions to the physics specific fields. + virtual + void applyGhostBoundaries(State& state, + StateDerivatives& derivs) override; + + // Enforce boundary conditions for the physics specific fields. + virtual + void enforceBoundaries(State& state, + StateDerivatives& derivs) override; + + Scalar nodeMotionCoefficient() const; + void nodeMotionCoefficient(const Scalar x); + + NodeMotionType nodeMotionType() const; + void nodeMotionType(NodeMotionType x); + + const FieldList& nodalVelocity() const; + const FieldList& DmassDt() const; + const FieldList& DthermalEnergyDt() const; + const FieldList& DmomentumDt() const; + const FieldList& DvolumeDt() const; + //const FieldList& HStretchTensor() const; + + const std::vector& pairMassFlux() const; + + //**************************************************************************** + // Methods required for restarting. + virtual std::string label() const override { return "MFVHydroBase" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + //**************************************************************************** +private: + + Scalar mNodeMotionCoefficient; + + NodeMotionType mNodeMotionType; + + FieldList mNodalVelocity; + FieldList mDmassDt; + FieldList mDthermalEnergyDt; + FieldList mDmomentumDt; + FieldList mDvolumeDt; + //FieldList mHStretchTensor; + + std::vector mPairMassFlux; + + // No default constructor, copying, or assignment. + MFVHydroBase(); + MFVHydroBase(const MFVHydroBase&); + MFVHydroBase& operator=(const MFVHydroBase&); +}; + +} + +#include "MFVHydroBaseInline.hh" + +#else + +// Forward declaration. +namespace Spheral { + template class MFVHydroBase; +} + +#endif diff --git a/src/GSPH/MFVHydroBaseInline.hh b/src/GSPH/MFVHydroBaseInline.hh new file mode 100644 index 000000000..7e945fa09 --- /dev/null +++ b/src/GSPH/MFVHydroBaseInline.hh @@ -0,0 +1,94 @@ +namespace Spheral { + + +//------------------------------------------------------------------------------ +// set/get for mesh motion coefficient +//------------------------------------------------------------------------------ +template +inline +typename Dimension::Scalar +MFVHydroBase::nodeMotionCoefficient() const { + return mNodeMotionCoefficient; +} + +template +inline +void +MFVHydroBase:: +nodeMotionCoefficient(typename Dimension::Scalar x) { + mNodeMotionCoefficient = x; +} + +//------------------------------------------------------------------------------ +// set/get mesh motion type +//------------------------------------------------------------------------------ +template +inline +NodeMotionType +MFVHydroBase:: +nodeMotionType() const { + return mNodeMotionType; +} + +template +inline +void +MFVHydroBase:: +nodeMotionType(NodeMotionType x) { + mNodeMotionType=x; +} + +//------------------------------------------------------------------------------ +// The internal state field lists. +//------------------------------------------------------------------------------ +template +inline +const FieldList& +MFVHydroBase:: +nodalVelocity() const { + return mNodalVelocity; +} + +template +inline +const FieldList& +MFVHydroBase:: +DmassDt() const { + return mDmassDt; +} +template +inline +const FieldList& +MFVHydroBase:: +DthermalEnergyDt() const { + return mDthermalEnergyDt; +} +template +inline +const FieldList& +MFVHydroBase:: +DmomentumDt() const { + return mDmomentumDt; +} +template +inline +const FieldList& +MFVHydroBase:: +DvolumeDt() const { + return mDvolumeDt; +} +template +inline +const std::vector& +MFVHydroBase:: +pairMassFlux() const { + return mPairMassFlux; +} +// template +// inline +// const FieldList& +// MFVHydroBase:: +// HStretchTensor() const { +// return mHStretchTensor; +// } +} \ No newline at end of file diff --git a/src/GSPH/MFVHydroBaseInst.cc.py b/src/GSPH/MFVHydroBaseInst.cc.py new file mode 100644 index 000000000..4f1510ced --- /dev/null +++ b/src/GSPH/MFVHydroBaseInst.cc.py @@ -0,0 +1,12 @@ +text = """ +//------------------------------------------------------------------------------ +// Explict instantiation. +//------------------------------------------------------------------------------ +#include "GSPH/MFVHydroBase.cc" +#include "GSPH/MFVEvaluateDerivatives.cc" +#include "Geometry/Dimension.hh" + +namespace Spheral { + template class MFVHydroBase< Dim< %(ndim)s > >; +} +""" diff --git a/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.cc b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.cc new file mode 100644 index 000000000..eb7a33a27 --- /dev/null +++ b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.cc @@ -0,0 +1,226 @@ +//---------------------------------Spheral++----------------------------------// +// CompatibleMFVSpecificThermalEnergyPolicy -- This is a generalization of the +// Lagrangian compatible energy scheme to ALE-based scheme with mass flux +// between nodes. +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh" +#include "GSPH/GSPHFieldNames.hh" + +#include "Hydro/HydroFieldNames.hh" + +#include "NodeList/NodeList.hh" +#include "NodeList/FluidNodeList.hh" + +#include "DataBase/IncrementState.hh" +#include "DataBase/DataBase.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" + +#include "Neighbor/ConnectivityMap.hh" + +#include "Field/Field.hh" +#include "Field/FieldList.hh" + +#include "Utilities/DBC.hh" +#include "Utilities/safeInv.hh" +#include "Utilities/SpheralFunctions.hh" + +#include +#include +using std::vector; +using std::numeric_limits; +using std::abs; +using std::min; +using std::max; + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor. +//------------------------------------------------------------------------------ +template +CompatibleMFVSpecificThermalEnergyPolicy:: +CompatibleMFVSpecificThermalEnergyPolicy(const DataBase& dataBase): + UpdatePolicyBase(), + mDataBasePtr(&dataBase){ +} + +//------------------------------------------------------------------------------ +// Destructor. +//------------------------------------------------------------------------------ +template +CompatibleMFVSpecificThermalEnergyPolicy:: +~CompatibleMFVSpecificThermalEnergyPolicy() { +} + +//------------------------------------------------------------------------------ +// Update the field. +//------------------------------------------------------------------------------ +template +void +CompatibleMFVSpecificThermalEnergyPolicy:: +update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double /*t*/, + const double /*dt*/) { + + KeyType fieldKey, nodeListKey; + StateBase::splitFieldKey(key, fieldKey, nodeListKey); + REQUIRE(fieldKey == HydroFieldNames::specificThermalEnergy and + nodeListKey == UpdatePolicyBase::wildcard()); + auto eps = state.fields(fieldKey, Scalar()); + const auto numFields = eps.numFields(); + + // constant we'll need for the weighting scheme + const auto tiny = numeric_limits::epsilon(); + + // Get the state fields. + const auto mass = state.fields(HydroFieldNames::mass, Scalar()); + const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); + const auto DmassDt = derivs.fields(IncrementState::prefix() + HydroFieldNames::mass, 0.0); + const auto DmomentumDt = derivs.fields(IncrementState::prefix() + GSPHFieldNames::momentum, Vector::zero); + const auto& pairAccelerations = derivs.getAny(HydroFieldNames::pairAccelerations, vector()); + const auto& pairDepsDt = derivs.getAny(HydroFieldNames::pairWork, vector()); + const auto& pairMassFlux = derivs.getAny(GSPHFieldNames::pairMassFlux, vector()); + + const auto& connectivityMap = mDataBasePtr->connectivityMap(); + const auto& pairs = connectivityMap.nodePairList(); + const auto npairs = pairs.size(); + + CHECK(pairAccelerations.size() == npairs); + CHECK(pairMassFlux.size() == npairs); + CHECK(pairDepsDt.size() == 2*npairs); + + auto DepsDt = derivs.fields(IncrementState::prefix() + GSPHFieldNames::thermalEnergy, 0.0); + DepsDt.Zero(); + + const auto hdt = 0.5*multiplier; + + // Walk all pairs and figure out the discrete work for each point +#pragma omp parallel + { + // Thread private variables + auto DepsDt_thread = DepsDt.threadCopy(); + +#pragma omp for + for (auto kk = 0u; kk < npairs; ++kk) { + const auto i = pairs[kk].i_node; + const auto j = pairs[kk].j_node; + const auto nodeListi = pairs[kk].i_list; + const auto nodeListj = pairs[kk].j_list; + + const auto& paccij = pairAccelerations[kk]; + const auto& DepsDt0i = pairDepsDt[2*kk]; + const auto& DepsDt0j = pairDepsDt[2*kk+1]; + const auto& massFlux = pairMassFlux[kk]; + + const auto mi = mass(nodeListi, i); + const auto pi = mi*velocity(nodeListi, i); + const auto& DPDti = DmomentumDt(nodeListi, i); + const auto& DmDti = DmassDt(nodeListi,i); + + const auto mj = mass(nodeListj, j); + const auto pj = mj*velocity(nodeListj, j); + const auto& DPDtj = DmomentumDt(nodeListj, j); + const auto& DmDtj = DmassDt(nodeListj,j); + + // half-step momenta + const auto pi12 = pi + DPDti*hdt; + const auto pj12 = pj + DPDtj*hdt; + //const auto pij = pi12 - pj12; + + // weighting scheme + const auto weighti = abs(DepsDt0i) + tiny; + const auto weightj = abs(DepsDt0j) + tiny; + const auto wi = weighti/(weighti+weightj); + + // safeInv + const auto invmi0 = safeInv(mi); + const auto invmj0 = safeInv(mj); + const auto invmi1 = safeInv(mi+DmDti*multiplier); + const auto invmj1 = safeInv(mj+DmDtj*multiplier); + + const Scalar delta_duij = (pi12*invmi1 - pj12*invmj1).dot(paccij) + + (pj.dot(pj)*invmj0*invmj1 - pi.dot(pi)*invmi0*invmi1) * massFlux*0.5 + - DepsDt0i-DepsDt0j; + + CHECK(wi >= 0.0 and wi <= 1.0); + CHECK(invmi0 >= 0.0); + CHECK(invmj0 >= 1.0); + + const auto depsi = (wi *delta_duij+DepsDt0i); + const auto depsj = ((1.0-wi)*delta_duij+DepsDt0j); + + // make conservative + DepsDt_thread(nodeListi, i) += depsi; + DepsDt_thread(nodeListj, j) += depsj; + + } + +#pragma omp critical + { + DepsDt_thread.threadReduce(); + } + } + +// // Now we can update the energy. + for (auto nodeListi = 0u; nodeListi < numFields; ++nodeListi) { + const auto n = eps[nodeListi]->numInternalElements(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + const auto m1 = mass(nodeListi,i)+DmassDt(nodeListi,i)*multiplier; + if (m1 > tiny) eps(nodeListi, i) += (DepsDt(nodeListi, i) - DmassDt(nodeListi, i)*eps(nodeListi, i)) * multiplier * safeInv(m1); + } + } + +} + +//------------------------------------------------------------------------------ +// Update the field using increments +//------------------------------------------------------------------------------ +template +void +CompatibleMFVSpecificThermalEnergyPolicy:: +updateAsIncrement(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) { + + KeyType fieldKey, nodeListKey; + StateBase::splitFieldKey(key, fieldKey, nodeListKey); + REQUIRE(fieldKey == HydroFieldNames::specificThermalEnergy and + nodeListKey == UpdatePolicyBase::wildcard()); + auto eps = state.fields(fieldKey, Scalar()); + + // Build an increment policy to use. + IncrementState fpolicy; + + // Do the deed for each of our Fields. + for (auto fptr: eps) { + fpolicy.updateAsIncrement(State::key(*fptr), + state, derivs, multiplier, t, dt); + } +} + +//------------------------------------------------------------------------------ +// Equivalence operator. +//------------------------------------------------------------------------------ +template +bool +CompatibleMFVSpecificThermalEnergyPolicy:: +operator==(const UpdatePolicyBase& rhs) const { + + // We're only equal if the other guy is also an increment operator. + const auto* rhsPtr = dynamic_cast*>(&rhs); + return rhsPtr != nullptr; +} + +} + diff --git a/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh new file mode 100644 index 000000000..2165468cd --- /dev/null +++ b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.hh @@ -0,0 +1,76 @@ +//---------------------------------Spheral++----------------------------------// +// CompatibleMFVSpecificThermalEnergyPolicy -- This is a generalization of the +// Lagrangian compatible energy scheme to ALE-based scheme with mass flux +// between nodes. +// +// J.M. Pearl 2024 +//----------------------------------------------------------------------------// + +#ifndef __Spheral_CompatibleMFVSpecificThermalEnergyPolicy_hh__ +#define __Spheral_CompatibleMFVSpecificThermalEnergyPolicy_hh__ + +#include "DataBase/UpdatePolicyBase.hh" + +#include + +namespace Spheral { + +// Forward declarations. +template class State; +template class StateDerivatives; +template class FluidNodeList; +template class FieldList; +template class DataBase; + +template +class CompatibleMFVSpecificThermalEnergyPolicy: + public UpdatePolicyBase { +public: + //--------------------------- Public Interface ---------------------------// + // Useful typedefs + using Scalar = typename Dimension::Scalar; + using Vector = typename Dimension::Vector; + using KeyType = typename UpdatePolicyBase::KeyType; + + // Constructors, destructor. + CompatibleMFVSpecificThermalEnergyPolicy(const DataBase& db); + virtual ~CompatibleMFVSpecificThermalEnergyPolicy(); + + // Overload the methods describing how to update Fields. + virtual void update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // If the derivative stored values for the pair-accelerations has not been updated, + // we need to just time advance normally. + virtual void updateAsIncrement(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // Equivalence. + virtual bool operator==(const UpdatePolicyBase& rhs) const override; + +private: + //--------------------------- Private Interface ---------------------------// + const DataBase* mDataBasePtr; + + CompatibleMFVSpecificThermalEnergyPolicy(const CompatibleMFVSpecificThermalEnergyPolicy& rhs); + CompatibleMFVSpecificThermalEnergyPolicy& operator=(const CompatibleMFVSpecificThermalEnergyPolicy& rhs); +}; + +} + +#else + +// Forward declaration. +namespace Spheral { + template class CompatibleMFVSpecificThermalEnergyPolicy; +} + +#endif diff --git a/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicyInst.cc.py b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicyInst.cc.py new file mode 100644 index 000000000..2559684fc --- /dev/null +++ b/src/GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicyInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "Geometry/Dimension.hh" +#include "GSPH/Policies/CompatibleMFVSpecificThermalEnergyPolicy.cc" + +namespace Spheral { + template class CompatibleMFVSpecificThermalEnergyPolicy >; +} +""" diff --git a/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.cc b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.cc new file mode 100644 index 000000000..1d240cede --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.cc @@ -0,0 +1,94 @@ +//---------------------------------Spheral++----------------------------------// +// MFVIncrementSpecificThermalEnergyPolicy -- This is a specialized increment +// policy for the specific thermal energy for schemes that allow +// for flux between nodes. The specific thermal energy is updated +// based on the time derivative of thermal energy. The mass and +// time derivative are needed to got from thermal to specific +// thermal. +// +// J.M. Pearl 2022 +//----------------------------------------------------------------------------// +// TODO: the edge case handing for m->0 needs to be improved to robustly +// handle void when full Eulerian. +//----------------------------------------------------------------------------// + +#include "GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.hh" +#include "GSPH/GSPHFieldNames.hh" +#include "DataBase/IncrementState.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" +#include "Field/Field.hh" +#include "Utilities/DBC.hh" +#include "Hydro/HydroFieldNames.hh" + +#include + +namespace Spheral { +//------------------------------------------------------------------------------ +// Constructors. +//------------------------------------------------------------------------------ +template +MFVIncrementSpecificThermalEnergyPolicy:: +MFVIncrementSpecificThermalEnergyPolicy(std::initializer_list depends): + FieldUpdatePolicy(depends){ +} + +//------------------------------------------------------------------------------ +// Destructor. +//------------------------------------------------------------------------------ +template +MFVIncrementSpecificThermalEnergyPolicy:: +~MFVIncrementSpecificThermalEnergyPolicy() { +} + +//------------------------------------------------------------------------------ +// Update the field. +//------------------------------------------------------------------------------ +template +void +MFVIncrementSpecificThermalEnergyPolicy:: +update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) { + + const auto tiny = std::numeric_limits::epsilon(); + + KeyType fieldKey, nodeListKey; + StateBase::splitFieldKey(key, fieldKey, nodeListKey); + + const auto massKey = StateBase::buildFieldKey(HydroFieldNames::mass, nodeListKey); + const auto derivFieldKey = StateBase::buildFieldKey(prefix() + GSPHFieldNames::thermalEnergy, nodeListKey); + + const auto& m = state.field(massKey, Scalar()); + auto& eps = state.field(key, Scalar()); + + const auto& DmDt = derivs.field(prefix() + massKey, Scalar()); + const auto& DmepsDt = derivs.field(derivFieldKey, Scalar()); + + const auto n = m.numInternalElements(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + const auto m1 = m(i)+DmDt(i)*multiplier; + if (m1 > tiny) eps(i) += (DmepsDt(i) - DmDt(i)*eps(i)) * multiplier * safeInv(m1); + } + +} + +//------------------------------------------------------------------------------ +// Equivalence operator. +//------------------------------------------------------------------------------ +template +bool +MFVIncrementSpecificThermalEnergyPolicy:: +operator==(const UpdatePolicyBase& rhs) const { + + // We're only equal if the other guy is also an replace operator. + const auto* rhsPtr = dynamic_cast*>(&rhs); + return rhsPtr != nullptr; +} + +} + diff --git a/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.hh b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.hh new file mode 100644 index 000000000..703237b7b --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.hh @@ -0,0 +1,70 @@ +//---------------------------------Spheral++----------------------------------// +// MFVIncrementSpecificThermalEnergyPolicy -- This is a specialized increment +// policy for the specific thermal energy for schemes that allow +// for flux between nodes. The specific thermal energy is updated +// based on the time derivative of thermal energy. The mass and +// time derivative are needed to got from thermal to specific +// thermal. +// +// J.M. Pearl 2022 +//----------------------------------------------------------------------------// +// TODO: the edge case handing for m->0 needs to be improved to robustly +// handle void when full Eulerian. +//----------------------------------------------------------------------------// + +#ifndef __Spheral_MFVIncrementSpecificThermalEnergyPolicy_hh__ +#define __Spheral_MFVIncrementSpecificThermalEnergyPolicy_hh__ + +#include "DataBase/FieldUpdatePolicy.hh" + +#include + +namespace Spheral { + +template +class MFVIncrementSpecificThermalEnergyPolicy: public FieldUpdatePolicy { +public: + + //--------------------------- Public Interface ---------------------------// + // Useful typedefs + using Scalar = typename Dimension::Scalar; + using Vector = typename Dimension::Vector; + using KeyType = typename FieldUpdatePolicy::KeyType; + + // Constructors, destructor. + MFVIncrementSpecificThermalEnergyPolicy(std::initializer_list depends={}); + ~MFVIncrementSpecificThermalEnergyPolicy(); + + // Overload the methods describing how to update FieldLists. + virtual void update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // Equivalence. + virtual bool operator==(const UpdatePolicyBase& rhs) const override; + + static const std::string prefix() { return "delta "; } + +private: + + const std::string mStateKey; + const std::string mDerivativeKey; + + //--------------------------- Private Interface ---------------------------// + MFVIncrementSpecificThermalEnergyPolicy(const MFVIncrementSpecificThermalEnergyPolicy& rhs); + MFVIncrementSpecificThermalEnergyPolicy& operator=(const MFVIncrementSpecificThermalEnergyPolicy& rhs); +}; + +} + +#else + +// Forward declaration. +namespace Spheral { + template class MFVIncrementSpecificThermalEnergyPolicy; +} + +#endif diff --git a/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicyInst.cc.py b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicyInst.cc.py new file mode 100644 index 000000000..975b8a38b --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicyInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "Geometry/Dimension.hh" +#include "GSPH/Policies/MFVIncrementSpecificThermalEnergyPolicy.cc" + +namespace Spheral { + template class MFVIncrementSpecificThermalEnergyPolicy>; +} +""" diff --git a/src/GSPH/Policies/MFVIncrementVelocityPolicy.cc b/src/GSPH/Policies/MFVIncrementVelocityPolicy.cc new file mode 100644 index 000000000..1083c332e --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementVelocityPolicy.cc @@ -0,0 +1,94 @@ +//---------------------------------Spheral++----------------------------------// +// MFVIncrementVelocityPolicy -- specialized policy for hydros that allow for mass +// flux between nodes. The momentum time derivative +// is used to update the velocity. The "hydro acceleration" +// is also added in to be compatible w/ phys packages +// that apply a pure acceleration. +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// +// TODO : HydroAcceleration needs to be added in +//----------------------------------------------------------------------------// + +#include "GSPH/Policies/MFVIncrementVelocityPolicy.hh" +#include "GSPH/GSPHFieldNames.hh" +#include "DataBase/IncrementState.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" +#include "Field/Field.hh" +#include "Utilities/DBC.hh" +#include "Hydro/HydroFieldNames.hh" + +#include + +namespace Spheral { +//------------------------------------------------------------------------------ +// Constructors. +//------------------------------------------------------------------------------ +template +MFVIncrementVelocityPolicy:: +MFVIncrementVelocityPolicy(std::initializer_list depends): + FieldUpdatePolicy(depends){ +} + +//------------------------------------------------------------------------------ +// Destructor. +//------------------------------------------------------------------------------ +template +MFVIncrementVelocityPolicy:: +~MFVIncrementVelocityPolicy() { +} + +//------------------------------------------------------------------------------ +// Update the field. +//------------------------------------------------------------------------------ +template +void +MFVIncrementVelocityPolicy:: +update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) { + + const auto tiny = std::numeric_limits::epsilon(); + + KeyType fieldKey, nodeListKey; + StateBase::splitFieldKey(key, fieldKey, nodeListKey); + + const auto massKey = StateBase::buildFieldKey(HydroFieldNames::mass, nodeListKey); + const auto momDerivFieldKey = StateBase::buildFieldKey(prefix() + GSPHFieldNames::momentum, nodeListKey); + //const auto accDerivFieldKey = StateBase::buildFieldKey(HydroFieldNames::hydroAcceleration, nodeListKey); + + const auto& m = state.field(massKey, Scalar()); + auto& v = state.field(key, Vector()); + + const auto& DmDt = derivs.field(prefix() + massKey, Scalar()); + const auto& DpDt = derivs.field(momDerivFieldKey, Vector()); + //const auto& DvDt = derivs.field(accDerivFieldKey, Vector()); + + const auto n = m.numInternalElements(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + const auto m1 = m(i)+DmDt(i)*multiplier; + const auto DpDti = DpDt(i); + if (m1 > tiny) v(i) += (DpDti - DmDt(i)*v(i)) * multiplier * safeInv(m1); + } +} + +//------------------------------------------------------------------------------ +// Equivalence operator. +//------------------------------------------------------------------------------ +template +bool +MFVIncrementVelocityPolicy:: +operator==(const UpdatePolicyBase& rhs) const { + + // We're only equal if the other guy is also an replace operator. + const auto* rhsPtr = dynamic_cast*>(&rhs); + return rhsPtr != nullptr; +} + +} + diff --git a/src/GSPH/Policies/MFVIncrementVelocityPolicy.hh b/src/GSPH/Policies/MFVIncrementVelocityPolicy.hh new file mode 100644 index 000000000..979499cbf --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementVelocityPolicy.hh @@ -0,0 +1,65 @@ +//---------------------------------Spheral++----------------------------------// +// MFVIncrementVelocityPolicy -- specialized policy for hydros that allow for mass +// flux between nodes. The momentum time derivative +// is used to update the velocity. The "hydro acceleration" +// is also added in to be compatible w/ phys packages +// that apply a pure acceleration. +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// +// TODO : HydroAcceleration needs to be added in +//----------------------------------------------------------------------------// + +#ifndef __Spheral_MFVIncrementVelocityPolicy_hh__ +#define __Spheral_MFVIncrementVelocityPolicy_hh__ + +#include "DataBase/FieldUpdatePolicy.hh" + +#include + +namespace Spheral { + +template +class MFVIncrementVelocityPolicy: public FieldUpdatePolicy { +public: + + //--------------------------- Public Interface ---------------------------// + // Useful typedefs + using Scalar = typename Dimension::Scalar; + using Vector = typename Dimension::Vector; + using KeyType = typename FieldUpdatePolicy::KeyType; + + // Constructors, destructor. + MFVIncrementVelocityPolicy(std::initializer_list depends={}); + ~MFVIncrementVelocityPolicy(); + + // Overload the methods describing how to update FieldLists. + virtual void update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // Equivalence. + virtual bool operator==(const UpdatePolicyBase& rhs) const override; + + static const std::string prefix() { return "delta "; } + +private: + + //--------------------------- Private Interface ---------------------------// + MFVIncrementVelocityPolicy(const MFVIncrementVelocityPolicy& rhs); + MFVIncrementVelocityPolicy& operator=(const MFVIncrementVelocityPolicy& rhs); +}; + +} + +#else + +// Forward declaration. +namespace Spheral { + template class MFVIncrementVelocityPolicy; +} + +#endif diff --git a/src/GSPH/Policies/MFVIncrementVelocityPolicyInst.cc.py b/src/GSPH/Policies/MFVIncrementVelocityPolicyInst.cc.py new file mode 100644 index 000000000..76a744576 --- /dev/null +++ b/src/GSPH/Policies/MFVIncrementVelocityPolicyInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "Geometry/Dimension.hh" +#include "GSPH/Policies/MFVIncrementVelocityPolicy.cc" + +namespace Spheral { + template class MFVIncrementVelocityPolicy>; +} +""" diff --git a/src/GSPH/Policies/MassFluxPolicy.cc b/src/GSPH/Policies/MassFluxPolicy.cc new file mode 100644 index 000000000..448be74e1 --- /dev/null +++ b/src/GSPH/Policies/MassFluxPolicy.cc @@ -0,0 +1,73 @@ +//---------------------------------Spheral++----------------------------------// +// MassFluxPolicy -- update method for ALE - based hydro schemes that allow +// for mass flux between nodes. +// +// J. M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "MassFluxPolicy.hh" +#include "Hydro/HydroFieldNames.hh" +#include "DataBase/IncrementState.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" +#include "Field/Field.hh" +#include "Utilities/DBC.hh" + +namespace Spheral { + + +//------------------------------------------------------------------------------ +// Constructor. +//------------------------------------------------------------------------------ +template +MassFluxPolicy:: +MassFluxPolicy(std::initializer_list depends): + IncrementState(depends) { +} + +//------------------------------------------------------------------------------ +// Destructor. +//------------------------------------------------------------------------------ +template +MassFluxPolicy:: +~MassFluxPolicy() { +} + +//------------------------------------------------------------------------------ +// Update the field. +//------------------------------------------------------------------------------ +template +void +MassFluxPolicy:: +update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double /*t*/, + const double /*dt*/) { + + auto& m = state.field(key, 0.0); + const auto& dmdt = derivs.field(this->prefix() + key, 0.0); + + const auto n = m.numInternalElements(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + m(i) += multiplier*(dmdt(i)); + } +} + +//------------------------------------------------------------------------------ +// Equivalence operator. +//------------------------------------------------------------------------------ +template +bool +MassFluxPolicy:: +operator==(const UpdatePolicyBase& rhs) const { + + // We're only equal if the other guy is also an increment operator. + const auto* rhsPtr = dynamic_cast*>(&rhs); + return rhsPtr != nullptr; +} + +} + diff --git a/src/GSPH/Policies/MassFluxPolicy.hh b/src/GSPH/Policies/MassFluxPolicy.hh new file mode 100644 index 000000000..5cf4a3ce3 --- /dev/null +++ b/src/GSPH/Policies/MassFluxPolicy.hh @@ -0,0 +1,63 @@ +//---------------------------------Spheral++----------------------------------// +// MassFluxPolicy -- update method for ALE - based hydro schemes that allow +// for mass flux between nodes. +// +// J. M. Pearl 2023 +//----------------------------------------------------------------------------// + +#ifndef __Spheral_MassFluxPolicy_hh__ +#define __Spheral_MassFluxPolicy_hh__ + +#include "DataBase/IncrementState.hh" + +#include + +namespace Spheral { + +// Forward declarations. +template class State; +template class StateDerivatives; +template class FluidNodeList; +template class FieldList; + +template +class MassFluxPolicy: + public IncrementState { +public: + //--------------------------- Public Interface ---------------------------// + // Useful typedefs + using Scalar = typename Dimension::Scalar; + using Vector = typename Dimension::Vector; + using KeyType = typename IncrementState::KeyType; + + // Constructors, destructor. + MassFluxPolicy(std::initializer_list depends = {}); + virtual ~MassFluxPolicy(); + + // Overload the methods describing how to update Fields. + virtual void update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // Equivalence. + virtual bool operator==(const UpdatePolicyBase& rhs) const override; + +private: + //--------------------------- Private Interface ---------------------------// + MassFluxPolicy(const MassFluxPolicy& rhs); + MassFluxPolicy& operator=(const MassFluxPolicy& rhs); +}; + +} + +#else + +// Forward declaration. +namespace Spheral { + template class MassFluxPolicy; +} + +#endif diff --git a/src/GSPH/RiemannSolvers/GHLLCInst.cc.py b/src/GSPH/Policies/MassFluxPolicyInst.cc.py similarity index 73% rename from src/GSPH/RiemannSolvers/GHLLCInst.cc.py rename to src/GSPH/Policies/MassFluxPolicyInst.cc.py index db174ede1..79a4346cb 100644 --- a/src/GSPH/RiemannSolvers/GHLLCInst.cc.py +++ b/src/GSPH/Policies/MassFluxPolicyInst.cc.py @@ -2,10 +2,10 @@ //------------------------------------------------------------------------------ // Explicit instantiation. //------------------------------------------------------------------------------ +#include "GSPH/Policies/MassFluxPolicy.cc" #include "Geometry/Dimension.hh" -#include "GSPH/RiemannSolvers/GHLLC.cc" namespace Spheral { - template class GHLLC >; + template class MassFluxPolicy >; } """ diff --git a/src/GSPH/RiemannSolvers/GHLLC.cc b/src/GSPH/RiemannSolvers/GHLLC.cc deleted file mode 100644 index e26aebcfd..000000000 --- a/src/GSPH/RiemannSolvers/GHLLC.cc +++ /dev/null @@ -1,182 +0,0 @@ -//---------------------------------Spheral++----------------------------------// -// GHLLC -- HLLC with gravitational source term -// -// J.M. Pearl 2021 -//----------------------------------------------------------------------------// -#include "FileIO/FileIO.hh" -#include "DataBase/DataBase.hh" -#include "DataBase/State.hh" -#include "DataBase/StateDerivatives.hh" - -#include "Field/FieldList.hh" -#include "Neighbor/ConnectivityMap.hh" - -#include "Hydro/HydroFieldNames.hh" -#include "GSPH/GSPHFieldNames.hh" - -#include "GSPH/WaveSpeeds/WaveSpeedBase.hh" -#include "GSPH/Limiters/LimiterBase.hh" -#include "GSPH/RiemannSolvers/HLLC.hh" -#include "GSPH/RiemannSolvers/GHLLC.hh" - -#include - -namespace Spheral { - -//------------------------------------------------------------------------------ -// Constructor -//------------------------------------------------------------------------------ -template -GHLLC:: -GHLLC(LimiterBase& slopeLimiter, - WaveSpeedBase& waveSpeed, - const bool linearReconstruction, - const typename Dimension::Vector gravitationalAcceleration): - HLLC(slopeLimiter, - waveSpeed, - linearReconstruction), - mGravitationalAcceleration(gravitationalAcceleration){ - -} - -//------------------------------------------------------------------------------ -// Destructor -//------------------------------------------------------------------------------ -template -GHLLC:: -~GHLLC(){} - - -//------------------------------------------------------------------------------ -// Interface State fluid hydro -//------------------------------------------------------------------------------ -// template -// void -// GHLLC:: -// interfaceState(const int i, -// const int j, -// const int nodelisti, -// const int nodelistj, -// const typename Dimension::Vector& ri, -// const typename Dimension::Vector& rj, -// const typename Dimension::Scalar& rhoi, -// const typename Dimension::Scalar& rhoj, -// const typename Dimension::Scalar& ci, -// const typename Dimension::Scalar& cj, -// const typename Dimension::Scalar& Pi, -// const typename Dimension::Scalar& Pj, -// const typename Dimension::Vector& vi, -// const typename Dimension::Vector& vj, -// typename Dimension::Scalar& Pstar, -// typename Dimension::Vector& vstar, -// typename Dimension::Scalar& rhostari, -// typename Dimension::Scalar& rhostarj) const{ - -// // pressure + linear grav contribution -// const auto rhogh = 0.5*(rhoi+rhoj)*mGravitationalAcceleration.dot(ri-rj); -// const auto p1i = Pi - rhogh; -// const auto p1j = Pj + rhogh; -// HLLC::interfaceState(i, -// j, -// nodelisti, -// nodelistj, -// ri, -// rj, -// rhoi, -// rhoj, -// ci, -// cj, -// p1i, -// p1j, -// vi, -// vj, -// Pstar, -// vstar, -// rhostari, -// rhostarj); -// }// Scalar interface class - - -template -void -GHLLC:: -interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const typename Dimension::Vector& ri, - const typename Dimension::Vector& rj, - const typename Dimension::Scalar& rhoi, - const typename Dimension::Scalar& rhoj, - const typename Dimension::Scalar& ci, - const typename Dimension::Scalar& cj, - const typename Dimension::Scalar& Pi, - const typename Dimension::Scalar& Pj, - const typename Dimension::Vector& vi, - const typename Dimension::Vector& vj, - const typename Dimension::Vector& DpDxi, - const typename Dimension::Vector& DpDxj, - const typename Dimension::Tensor& DvDxi, - const typename Dimension::Tensor& DvDxj, - typename Dimension::Scalar& Pstar, - typename Dimension::Vector& vstar, - typename Dimension::Scalar& rhostari, - typename Dimension::Scalar& rhostarj) const{ - - // pressure + linear grav contribution - const auto rhogh = 0.5*(rhoi+rhoj)*mGravitationalAcceleration.dot(ri-rj); - const auto p1i = Pi - rhogh; - const auto p1j = Pj + rhogh; - HLLC::interfaceState(i, - j, - nodelisti, - nodelistj, - ri, - rj, - rhoi, - rhoj, - ci, - cj, - p1i, - p1j, - vi, - vj, - DpDxi, - DpDxj, - DvDxi, - DvDxj, - Pstar, - vstar, - rhostari, - rhostarj); -}// Scalar interface class - -template -void -GHLLC:: -interfaceState(const int /*i*/, - const int /*j*/, - const int /*nodelisti*/, - const int /*nodelistj*/, - const Vector& /*ri*/, - const Vector& /*rj*/, - const Scalar& /*rhoi*/, - const Scalar& /*rhoj*/, - const Scalar& /*ci*/, - const Scalar& /*cj*/, - const Scalar& /*Pi*/, - const Scalar& /*Pj*/, - const Vector& /*vi*/, - const Vector& /*vj*/, - const SymTensor& /*Si*/, - const SymTensor& /*Sj*/, - const Tensor& /*Di*/, - const Tensor& /*Dj*/, - Vector& /*Tstar*/, - Vector& /*vstar*/) const{ - - - -} - -} // spheral namespace \ No newline at end of file diff --git a/src/GSPH/RiemannSolvers/GHLLCInline.hh b/src/GSPH/RiemannSolvers/GHLLCInline.hh deleted file mode 100644 index ad222696e..000000000 --- a/src/GSPH/RiemannSolvers/GHLLCInline.hh +++ /dev/null @@ -1,21 +0,0 @@ -namespace Spheral { - - -//-------------------------------------------------------- -// setter/getter for the gravitational acceleration -//-------------------------------------------------------- -template -typename Dimension::Vector -GHLLC:: -gravitationalAcceleration() const{ - return mGravitationalAcceleration; -} - -template -void -GHLLC:: -gravitationalAcceleration(const typename Dimension::Vector g) { - mGravitationalAcceleration = g; -} - -} \ No newline at end of file diff --git a/src/GSPH/RiemannSolvers/HLLC.cc b/src/GSPH/RiemannSolvers/HLLC.cc index f108e8a83..17d824199 100644 --- a/src/GSPH/RiemannSolvers/HLLC.cc +++ b/src/GSPH/RiemannSolvers/HLLC.cc @@ -1,19 +1,11 @@ //---------------------------------Spheral++----------------------------------// // HLLC -- approximate riemann solver -// Toro E.F., Spruce M., Speares W., (1994) "Restoration of the Contact Surface in -// the HLL-Riemann Solver," Shock Waves, 4:25-34 +// Toro E.F., Spruce M., Speares W., (1994) "Restoration of the Contact +// Surface in the HLL-Riemann Solver," Shock Waves, 4:25-34 // // J.M. Pearl 2021 //----------------------------------------------------------------------------// -#include "FileIO/FileIO.hh" -#include "DataBase/DataBase.hh" -#include "DataBase/State.hh" -#include "DataBase/StateDerivatives.hh" - -#include "Field/FieldList.hh" -#include "Neighbor/ConnectivityMap.hh" - #include "Hydro/HydroFieldNames.hh" #include "GSPH/GSPHFieldNames.hh" @@ -47,134 +39,16 @@ HLLC:: ~HLLC(){} -//------------------------------------------------------------------------------ -// Interface State scalar -//------------------------------------------------------------------------------ -// template -// void -// HLLC:: -// interfaceState(const int i, -// const int j, -// const int nodelisti, -// const int nodelistj, -// const typename Dimension::Vector& ri, -// const typename Dimension::Vector& rj, -// const typename Dimension::Scalar& rhoi, -// const typename Dimension::Scalar& rhoj, -// const typename Dimension::Scalar& ci, -// const typename Dimension::Scalar& cj, -// const typename Dimension::Scalar& Pi, -// const typename Dimension::Scalar& Pj, -// const typename Dimension::Vector& vi, -// const typename Dimension::Vector& vj, -// typename Dimension::Scalar& Pstar, -// typename Dimension::Vector& vstar, -// typename Dimension::Scalar& /*rhostari*/, -// typename Dimension::Scalar& /*rhostarj*/) const{ - -// Scalar Si, Sj; - -// const auto tiny = std::numeric_limits::epsilon(); - -// const auto& limiter = this->limiter(); -// const auto& waveSpeedObject = this->waveSpeed(); - -// const auto& DpDx0 = this->DpDx(); -// const auto& DvDx0 = this->DvDx(); - -// const auto rij = ri - rj; -// const auto rhatij = rij.unitVector(); - -// vstar = 0.5*(vi+vj); -// Pstar = 0.5*(Pi+Pj); - -// if (ci > tiny or cj > tiny){ - - -// // default to nodal values -// auto v1i = vi; -// auto v1j = vj; - -// auto p1i = Pi; -// auto p1j = Pj; - -// // linear reconstruction -// if(this->linearReconstruction()){ - -// // gradients -// const auto DvDxi = DvDx0(nodelisti,i); -// const auto DvDxj = DvDx0(nodelistj,j); -// const auto DpDxi = DpDx0(nodelisti,i); -// const auto DpDxj = DpDx0(nodelistj,j); - -// // gradients along line of action -// if (true){ -// this->linearReconstruction(ri,rj, Pi,Pj, DpDxi,DpDxj, -// p1i,p1j); -// this->linearReconstruction(ri,rj, vi,vj, DvDxi,DvDxj, -// v1i,v1j); -// }else{ - -// const auto xij = 0.5*(rij); -// const auto Dpi = DpDxi.dot(xij); -// const auto Dpj = DpDxj.dot(xij); -// const auto Dvi = DvDxi.dot(xij); -// const auto Dvj = DvDxj.dot(xij); -// const auto Dui = Dvi.dot(rhatij); -// const auto Duj = Dvj.dot(rhatij); -// //const auto Dp0 = 0.5*(Pi-Pj); -// //const auto Du0 = 0.5*(vi-vj).dot(rhatij); -// const auto rui = Dui/(sgn(Duj)*std::max(tiny, abs(Duj))); -// const auto ruj = Duj/(sgn(Dui)*std::max(tiny, abs(Dui))); -// const auto xu = std::min(rui,ruj); -// const auto phiu = limiter.slopeLimiter(xu); - -// const auto rpi = Dpi/(sgn(Dpj)*std::max(tiny, abs(Dpj))); -// const auto rpj = Dpj/(sgn(Dpi)*std::max(tiny, abs(Dpi))); -// const auto xp = std::min(rpi,rpj); -// const auto phip = limiter.slopeLimiter(xp); - -// v1i = vi - phiu * Dvi; -// v1j = vj + phiu * Dvj; -// p1i = Pi - phip * Dpi; -// p1j = Pj + phip * Dpj; -// } - -// } - -// const auto ui = v1i.dot(rhatij); -// const auto uj = v1j.dot(rhatij); -// const auto wi = v1i - ui*rhatij; -// const auto wj = v1j - uj*rhatij; - -// waveSpeedObject.waveSpeed(rhoi,rhoj,ci,cj,ui,uj,Si,Sj); - -// const auto denom = safeInv(Si - Sj); - -// const auto ustar = (Si*ui - Sj*uj - p1i + p1j )*denom; -// const auto wstar = (Si*wi - Sj*wj)*denom; -// vstar = ustar*rhatij + wstar; -// Pstar = Sj * (ustar-uj) + p1j; - -// }else{ // if ci & cj too small punt to normal av -// const auto uij = std::min((vi-vj).dot(rhatij),0.0); -// Pstar += 0.25 * (rhoi+rhoj) * (uij*uij); -// } -// }// Scalar interface class - - //------------------------------------------------------------------------------ // Interface State fluid hydro //------------------------------------------------------------------------------ template void HLLC:: -interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const typename Dimension::Vector& ri, +interfaceState(const typename Dimension::Vector& ri, const typename Dimension::Vector& rj, + const typename Dimension::SymTensor& Hi, + const typename Dimension::SymTensor& Hj, const typename Dimension::Scalar& rhoi, const typename Dimension::Scalar& rhoj, const typename Dimension::Scalar& ci, @@ -183,14 +57,16 @@ interfaceState(const int i, const typename Dimension::Scalar& Pj, const typename Dimension::Vector& vi, const typename Dimension::Vector& vj, + const typename Dimension::Vector& DrhoDxi, + const typename Dimension::Vector& DrhoDxj, const typename Dimension::Vector& DpDxi, const typename Dimension::Vector& DpDxj, const typename Dimension::Tensor& DvDxi, const typename Dimension::Tensor& DvDxj, typename Dimension::Scalar& Pstar, typename Dimension::Vector& vstar, - typename Dimension::Scalar& /*rhostari*/, - typename Dimension::Scalar& /*rhostarj*/) const{ + typename Dimension::Scalar& rhostari, + typename Dimension::Scalar& rhostarj) const{ Scalar Si, Sj; @@ -203,6 +79,8 @@ interfaceState(const int i, vstar = 0.5*(vi+vj); Pstar = 0.5*(Pi+Pj); + rhostari = rhoi; + rhostarj = rhoj; if (ci > tiny or cj > tiny){ @@ -214,14 +92,19 @@ interfaceState(const int i, auto p1i = Pi; auto p1j = Pj; + //auto rho1i = rhoi; + //auto rho1j = rhoj; + // linear reconstruction if(this->linearReconstruction()){ // gradients along line of action - this->linearReconstruction(ri,rj, Pi,Pj, DpDxi,DpDxj, - p1i,p1j); - this->linearReconstruction(ri,rj, vi,vj, DvDxi,DvDxj, - v1i,v1j); + //this->linearReconstruction(ri,rj, rhoi,rhoj,DrhoDxi,DrhoDxj, //inputs + // rho1i,rho1j); //outputs + this->linearReconstruction(ri,rj, Pi,Pj,DpDxi,DpDxj, //inputs + p1i,p1j); //outputs + this->linearReconstruction(ri,rj, vi,vj,DvDxi,DvDxj, //inputs + v1i,v1j); //outputs } @@ -230,7 +113,8 @@ interfaceState(const int i, const auto wi = v1i - ui*rhatij; const auto wj = v1j - uj*rhatij; - waveSpeedObject.waveSpeed(rhoi,rhoj,ci,cj,ui,uj,Si,Sj); + waveSpeedObject.waveSpeed(rhoi,rhoj,ci,cj,ui,uj, //inputs + Si,Sj); //outputs const auto denom = safeInv(Si - Sj); @@ -238,23 +122,24 @@ interfaceState(const int i, const auto wstar = (Si*wi - Sj*wj)*denom; vstar = ustar*rhatij + wstar; Pstar = Sj * (ustar-uj) + p1j; + //rhostari = rho1i;// * (Si - ui)*safeInv(Si-ustar); + //rhostarj = rho1j;// * (Sj - uj)*safeInv(Sj-ustar); }else{ // if ci & cj too small punt to normal av const auto uij = std::min((vi-vj).dot(rhatij),0.0); Pstar += 0.25 * (rhoi+rhoj) * (uij*uij); } + }// Scalar interface class template void HLLC:: -interfaceState(const int /*i*/, - const int /*j*/, - const int /*nodelisti*/, - const int /*nodelistj*/, - const Vector& /*ri*/, +interfaceState(const Vector& /*ri*/, const Vector& /*rj*/, + const SymTensor& /*Hi*/, + const SymTensor& /*Hj*/, const Scalar& /*rhoi*/, const Scalar& /*rhoj*/, const Scalar& /*ci*/, diff --git a/src/GSPH/RiemannSolvers/HLLC.hh b/src/GSPH/RiemannSolvers/HLLC.hh index 8432a2fef..a401cc8a9 100644 --- a/src/GSPH/RiemannSolvers/HLLC.hh +++ b/src/GSPH/RiemannSolvers/HLLC.hh @@ -37,34 +37,11 @@ public: ~HLLC(); - // virtual - // void interfaceState(const int i, - // const int j, - // const int nodelisti, - // const int nodelistj, - // const Vector& ri, - // const Vector& rj, - // const Scalar& rhoi, - // const Scalar& rhoj, - // const Scalar& ci, - // const Scalar& cj, - // const Scalar& sigmai, - // const Scalar& sigmaj, - // const Vector& vi, - // const Vector& vj, - // Scalar& Pstar, - // Vector& vstar, - // Scalar& rhostari, - // Scalar& rhostarj) const override; - - // ^ temporary class to wrap the above ^ virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, @@ -73,6 +50,8 @@ public: const Scalar& sigmaj, const Vector& vi, const Vector& vj, + const Vector& DrhoDxi, + const Vector& DrhoDxj, const Vector& DpDxi, const Vector& DpDxj, const Tensor& DvDxi, @@ -84,12 +63,10 @@ public: virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, diff --git a/src/GSPH/RiemannSolvers/RiemannSolverBase.cc b/src/GSPH/RiemannSolvers/RiemannSolverBase.cc index 4e68e2f91..6fe5490c5 100644 --- a/src/GSPH/RiemannSolvers/RiemannSolverBase.cc +++ b/src/GSPH/RiemannSolvers/RiemannSolverBase.cc @@ -35,8 +35,6 @@ RiemannSolverBase(LimiterBase& slopeLimiter, mSlopeLimiter(slopeLimiter), mWaveSpeed(waveSpeed), mLinearReconstruction(linearReconstruction){ - //mDpDx(FieldStorageType::CopyFields), - //mDvDx(FieldStorageType::CopyFields){ } //------------------------------------------------------------------------------ @@ -52,92 +50,15 @@ RiemannSolverBase:: template void RiemannSolverBase:: -initialize(const DataBase& dataBase, - const State& state, - const StateDerivatives& derivs, - typename RiemannSolverBase::ConstBoundaryIterator boundaryBegin, - typename RiemannSolverBase::ConstBoundaryIterator boundaryEnd, +initialize(const DataBase& /*dataBase*/, + const State& /*state*/, + const StateDerivatives& /*derivs*/, + typename RiemannSolverBase::ConstBoundaryIterator /*boundaryBegin*/, + typename RiemannSolverBase::ConstBoundaryIterator /*boundaryEnd*/, const typename Dimension::Scalar /*time*/, const typename Dimension::Scalar /*dt*/, const TableKernel& /*W*/){ - // if(mLinearReconstruction){ - // dataBase.resizeFluidFieldList(mDpDx,Vector::zero,GSPHFieldNames::RiemannPressureGradient0,true); - // dataBase.resizeFluidFieldList(mDvDx,Tensor::zero,GSPHFieldNames::RiemannVelocityGradient0,true); - - // //const auto& DpDx0 = derivs.fields( GSPHFieldNames::pressureGradient, Vector::zero); - // //const auto& DpDxRaw0 = derivs.fields( GSPHFieldNames::pressureGradient+"RAW", Vector::zero); - // const auto& DvDx0 = derivs.fields( HydroFieldNames::velocityGradient,Tensor::zero); - // //const auto& localDvDx0 = derivs.fields( HydroFieldNames::internalVelocityGradient,Tensor::zero); - // //const auto& DvDxRaw0 = derivs.fields( HydroFieldNames::velocityGradient+"RAW",Tensor::zero); - // const auto& DvDt0 = derivs.fields(HydroFieldNames::hydroAcceleration, Vector::zero); - // const auto& rho0 = state.fields(HydroFieldNames::massDensity, 0.0); - - // const auto& connectivityMap = dataBase.connectivityMap(); - // const auto& nodeLists = connectivityMap.nodeLists(); - // const auto numNodeLists = nodeLists.size(); - - // // copy from previous time step - // for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { - // const auto& nodeList = nodeLists[nodeListi]; - // const auto ni = nodeList->numInternalNodes(); - // #pragma omp parallel for - // for (auto i = 0u; i < ni; ++i) { - // //const auto localDvDxi = localDvDx0(nodeListi,i); - // const auto DvDxi = DvDx0(nodeListi,i); - // //const auto DvDxRawi = DvDxRaw0(nodeListi,i); - // //const auto DpDxi = DpDx0(nodeListi,i); - // //const auto DpDxRawi = DpDxRaw0(nodeListi,i); - // const auto DvDti = DvDt0(nodeListi,i); - // const auto rhoi = rho0(nodeListi,i); - - // // this'll need some cleaning - // // switch(mGradientType){ - // // case GradientType::RiemannGradient: // default grad based on riemann soln - // mDvDx(nodeListi,i) = DvDxi; - // mDpDx(nodeListi,i) = -rhoi*DvDti; - // // break; - // // case GradientType::HydroAccelerationGradient: // based on hydro accel for DpDx - // // mDvDx(nodeListi,i) = DvDxi; - // // mDpDx(nodeListi,i) = -rhoi*DvDti; - // // break; - // // case GradientType::SPHGradient: // raw gradients - // // mDvDx(nodeListi,i) = DvDxRawi; - // // mDpDx(nodeListi,i) = DpDxRawi; - // // break; - // // case GradientType::MixedMethodGradient: // raw gradient for P riemann gradient for v - // // mDvDx(nodeListi,i) = DvDxi; - // // mDpDx(nodeListi,i) = DpDxRawi; - // // break; - // // case GradientType::OnlyDvDxGradient: // raw gradients - // // mDvDx(nodeListi,i) = DvDxi; - // // mDpDx(nodeListi,i) = Vector::zero; - // // break; - // // case GradientType::LocalDvDxGradient: // local velocity gradient - // // mDvDx(nodeListi,i) = localDvDxi; - // // mDpDx(nodeListi,i) = DpDxi; - // // break; - // // default : - // // mDvDx(nodeListi,i) = Tensor::zero; - // // mDpDx(nodeListi,i) = Vector::zero; - - // // } - // } - // } - - // for (auto boundItr = boundaryBegin; - // boundItr != boundaryEnd; - // ++boundItr) { - // (*boundItr)->applyFieldListGhostBoundary(mDpDx); - // (*boundItr)->applyFieldListGhostBoundary(mDvDx); - // } - - // for (auto boundItr = boundaryBegin; - // boundItr != boundaryEnd; - // ++boundItr) (*boundItr)->finalizeGhostBoundary(); - - // } // if LinearReconstruction - } // initialize method @@ -227,39 +148,15 @@ linearReconstruction(const typename Dimension::Vector& ri, //------------------------------------------------------------------------------ // default to non-op //------------------------------------------------------------------------------ -// template -// void -// RiemannSolverBase:: -// interfaceState(const int /*i*/, -// const int /*j*/, -// const int /*nodelisti*/, -// const int /*nodelistj*/, -// const Vector& /*ri*/, -// const Vector& /*rj*/, -// const Scalar& /*rhoi*/, -// const Scalar& /*rhoj*/, -// const Scalar& /*ci*/, -// const Scalar& /*cj*/, -// const Scalar& /*sigmai*/, -// const Scalar& /*sigmaj*/, -// const Vector& /*vi*/, -// const Vector& /*vj*/, -// Scalar& /*Pstar*/, -// Vector& /*vstar*/, -// Scalar& /*rhostari*/, -// Scalar& /*rhostarj*/) const{ -// } template void RiemannSolverBase:: -interfaceState(const int /*i*/, - const int /*j*/, - const int /*nodelisti*/, - const int /*nodelistj*/, - const Vector& /*ri*/, +interfaceState(const Vector& /*ri*/, const Vector& /*rj*/, + const SymTensor& /*Hi*/, + const SymTensor& /*Hj*/, const Scalar& /*rhoi*/, const Scalar& /*rhoj*/, const Scalar& /*ci*/, @@ -268,6 +165,8 @@ interfaceState(const int /*i*/, const Scalar& /*Pj*/, const Vector& /*vi*/, const Vector& /*vj*/, + const Vector& /*DrhoDxi*/, + const Vector& /*DrhoDxj*/, const Vector& /*DpDxi*/, const Vector& /*DpDxj*/, const Tensor& /*DvDxi*/, @@ -285,12 +184,10 @@ interfaceState(const int /*i*/, template void RiemannSolverBase:: -interfaceState(const int /*i*/, - const int /*j*/, - const int /*nodelisti*/, - const int /*nodelistj*/, - const Vector& /*ri*/, +interfaceState(const Vector& /*ri*/, const Vector& /*rj*/, + const SymTensor& /*Hi*/, + const SymTensor& /*Hj*/, const Scalar& /*rhoi*/, const Scalar& /*rhoj*/, const Scalar& /*ci*/, diff --git a/src/GSPH/RiemannSolvers/RiemannSolverBase.hh b/src/GSPH/RiemannSolvers/RiemannSolverBase.hh index fbdd996d2..e5aacffdf 100644 --- a/src/GSPH/RiemannSolvers/RiemannSolverBase.hh +++ b/src/GSPH/RiemannSolvers/RiemannSolverBase.hh @@ -46,33 +46,12 @@ public: const Scalar dt, const TableKernel& W); - // virtual - // void interfaceState(const int i, - // const int j, - // const int nodelisti, - // const int nodelistj, - // const Vector& ri, - // const Vector& rj, - // const Scalar& rhoi, - // const Scalar& rhoj, - // const Scalar& ci, - // const Scalar& cj, - // const Scalar& Pi, - // const Scalar& Pj, - // const Vector& vi, - // const Vector& vj, - // Scalar& Pstar, - // Vector& vstar, - // Scalar& rhostari, - // Scalar& rhostarj) const; virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, @@ -81,6 +60,8 @@ public: const Scalar& Pj, const Vector& vi, const Vector& vj, + const Vector& DrhoDxi, + const Vector& DrhoDxj, const Vector& DpDxi, const Vector& DpDxj, const Tensor& DvDxi, @@ -92,12 +73,10 @@ public: virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, @@ -118,9 +97,6 @@ public: bool linearReconstruction() const; void linearReconstruction(bool x); - - // GradientType gradientType() const; - // void gradientType(GradientType x); virtual void linearReconstruction(const Vector& ri, @@ -142,22 +118,13 @@ public: Vector& ytildei, Vector& ytildej) const; - // we'll want the ability to modify these (make better) - // FieldList& DpDx(); - // FieldList& DvDx(); - // const FieldList& DpDx() const; - // const FieldList& DvDx() const; private: LimiterBase& mSlopeLimiter; WaveSpeedBase& mWaveSpeed; bool mLinearReconstruction; - //GradientType mGradientType; - - // FieldList mDpDx; - // FieldList mDvDx; }; diff --git a/src/GSPH/RiemannSolvers/RiemannSolverBaseInline.hh b/src/GSPH/RiemannSolvers/RiemannSolverBaseInline.hh index 6df4c6f3a..d04dec1d2 100644 --- a/src/GSPH/RiemannSolvers/RiemannSolverBaseInline.hh +++ b/src/GSPH/RiemannSolvers/RiemannSolverBaseInline.hh @@ -41,44 +41,4 @@ linearReconstruction(bool x) { mLinearReconstruction=x; } - -//------------------------------------------------------------------------------ -// field getters -//------------------------------------------------------------------------------ - - -// template -// inline -// FieldList& -// RiemannSolverBase:: -// DpDx() { -// return mDpDx; -// } - -// template -// inline -// FieldList& -// RiemannSolverBase:: -// DvDx() { -// return mDvDx; -// } - - -// template -// inline -// const FieldList& -// RiemannSolverBase:: -// DpDx() const { -// return mDpDx; -// } - -// template -// inline -// const FieldList& -// RiemannSolverBase:: -// DvDx() const { -// return mDvDx; -// } - - } \ No newline at end of file diff --git a/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.cc b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.cc new file mode 100644 index 000000000..a67dbfd26 --- /dev/null +++ b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.cc @@ -0,0 +1,135 @@ +//---------------------------------Spheral++----------------------------------// +// SecondOrderArtificialViscosity +// Frontiere, Raskin, Owen (2017) "CRKSPH:- A Conservative Reproducing Kernel +// Smoothed Particle Hydrodynamics Scheme," J. Comp. Phys. +// +// This is a reimplementation of the LimitedArtificialViscosity class as a +// derivative of RiemannSolverBase so it can be used with GSPH derived +// classes +// +// J.M. Pearl 2021 +//----------------------------------------------------------------------------// + +#include "Hydro/HydroFieldNames.hh" +#include "GSPH/GSPHFieldNames.hh" + +#include "GSPH/WaveSpeeds/WaveSpeedBase.hh" +#include "GSPH/Limiters/LimiterBase.hh" +#include "GSPH/RiemannSolvers/SecondOrderArtificialViscosity.hh" + +#include + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +template +SecondOrderArtificialViscosity:: +SecondOrderArtificialViscosity(const Scalar Cl, + const Scalar Cq, + LimiterBase& slopeLimiter, + WaveSpeedBase& waveSpeed, + const bool linearReconstruction): + RiemannSolverBase(slopeLimiter, + waveSpeed, + linearReconstruction), + mCl(Cl), + mCq(Cq){ + +} + +//------------------------------------------------------------------------------ +// Destructor +//------------------------------------------------------------------------------ +template +SecondOrderArtificialViscosity:: +~SecondOrderArtificialViscosity(){} + + +//------------------------------------------------------------------------------ +// Interface State fluid hydro +//------------------------------------------------------------------------------ +template +void +SecondOrderArtificialViscosity:: +interfaceState(const typename Dimension::Vector& ri, + const typename Dimension::Vector& rj, + const typename Dimension::SymTensor& Hi, + const typename Dimension::SymTensor& Hj, + const typename Dimension::Scalar& rhoi, + const typename Dimension::Scalar& rhoj, + const typename Dimension::Scalar& ci, + const typename Dimension::Scalar& cj, + const typename Dimension::Scalar& Pi, + const typename Dimension::Scalar& Pj, + const typename Dimension::Vector& vi, + const typename Dimension::Vector& vj, + const typename Dimension::Vector& /*DrhoDxi*/, + const typename Dimension::Vector& /*DrhoDxj*/, + const typename Dimension::Vector& /*DpDxi*/, + const typename Dimension::Vector& /*DpDxj*/, + const typename Dimension::Tensor& DvDxi, + const typename Dimension::Tensor& DvDxj, + typename Dimension::Scalar& Pstar, + typename Dimension::Vector& vstar, + typename Dimension::Scalar& /*rhostari*/, + typename Dimension::Scalar& /*rhostarj*/) const{ + + + const auto tiny = std::numeric_limits::epsilon(); + + const Vector rij = ri - rj; + const Vector rhatij = rij.unitVector(); + const Vector etaij = 0.5*(Hi+Hj)*rij; + + // default to nodal values + Vector v1i = vi; + Vector v1j = vj; + + // linear reconstruction + if(this->linearReconstruction()){ + + this->linearReconstruction(ri,rj, vi,vj,DvDxi,DvDxj, //inputs + v1i,v1j); //outputs + + } + const Vector vij = v1i-v1j; + const Scalar muij = std::max(0.0,-vij.dot(etaij)/(etaij.magnitude2() + tiny)); + const Scalar cij = 0.5*(ci+cj); + const Scalar rhoij = 2*rhoi*rhoj/(rhoi+rhoj); + Pstar = 0.5*(Pi+Pj) + + rhoij*muij*(this->Cl()*cij + +this->Cq()*muij); + vstar = 0.5*(vi+vj); + +}// Scalar interface class + + +template +void +SecondOrderArtificialViscosity:: +interfaceState(const Vector& /*ri*/, + const Vector& /*rj*/, + const SymTensor& /*Hi*/, + const SymTensor& /*Hj*/, + const Scalar& /*rhoi*/, + const Scalar& /*rhoj*/, + const Scalar& /*ci*/, + const Scalar& /*cj*/, + const Scalar& /*Pi*/, + const Scalar& /*Pj*/, + const Vector& /*vi*/, + const Vector& /*vj*/, + const SymTensor& /*Si*/, + const SymTensor& /*Sj*/, + const Tensor& /*Di*/, + const Tensor& /*Dj*/, + Vector& /*Tstar*/, + Vector& /*vstar*/) const{ + + + +} + +} // spheral namespace \ No newline at end of file diff --git a/src/GSPH/RiemannSolvers/GHLLC.hh b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.hh similarity index 56% rename from src/GSPH/RiemannSolvers/GHLLC.hh rename to src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.hh index 951233334..3b7d91a0b 100644 --- a/src/GSPH/RiemannSolvers/GHLLC.hh +++ b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosity.hh @@ -1,10 +1,19 @@ //---------------------------------Spheral++----------------------------------// -// GHLLC -- HLLC solver w/ gravitational source term +// SecondOrderArtificialViscosity +// Frontiere, Raskin, Owen (2017) "CRKSPH:- A Conservative Reproducing Kernel +// Smoothed Particle Hydrodynamics Scheme," J. Comp. Phys. +// +// This is a reimplementation of the LimitedArtificialViscosity class as a +// derivative of RiemannSolverBase so it can be used with GSPH derived +// classes +// +// J.M. Pearl 2021 //----------------------------------------------------------------------------// -#ifndef __Spheral_GHLLC_hh__ -#define __Spheral_GHLLC_hh__ -#include "HLLC.hh" +#ifndef __Spheral_SecondOrderArtificialViscosity_hh__ +#define __Spheral_SecondOrderArtificialViscosity_hh__ + +#include "RiemannSolverBase.hh" namespace Spheral { @@ -17,7 +26,7 @@ template class Field; template class FieldList; template -class GHLLC : public HLLC { +class SecondOrderArtificialViscosity : public RiemannSolverBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; @@ -26,40 +35,20 @@ class GHLLC : public HLLC { public: - GHLLC(LimiterBase& slopeLimiter, - WaveSpeedBase& waveSpeedBase, - const bool linearReconstruction, - const Vector gravitationalAcceleration); - - ~GHLLC(); - - // virtual - // void interfaceState(const int i, - // const int j, - // const int nodelisti, - // const int nodelistj, - // const Vector& ri, - // const Vector& rj, - // const Scalar& rhoi, - // const Scalar& rhoj, - // const Scalar& ci, - // const Scalar& cj, - // const Scalar& sigmai, - // const Scalar& sigmaj, - // const Vector& vi, - // const Vector& vj, - // Scalar& Pstar, - // Vector& vstar, - // Scalar& rhostari, - // Scalar& rhostarj) const override; + SecondOrderArtificialViscosity(const Scalar Cl, + const Scalar Cq, + LimiterBase& slopeLimiter, + WaveSpeedBase& waveSpeedBase, + const bool linearReconstruction); + + ~SecondOrderArtificialViscosity(); + virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, @@ -68,6 +57,8 @@ public: const Scalar& sigmaj, const Vector& vi, const Vector& vj, + const Vector& DrhoDxi, + const Vector& DrhoDxj, const Vector& DpDxi, const Vector& DpDxj, const Tensor& DvDxi, @@ -77,13 +68,12 @@ public: Scalar& rhostari, Scalar& rhostarj) const override; + virtual - void interfaceState(const int i, - const int j, - const int nodelisti, - const int nodelistj, - const Vector& ri, + void interfaceState(const Vector& ri, const Vector& rj, + const SymTensor& Hi, + const SymTensor& Hj, const Scalar& rhoi, const Scalar& rhoj, const Scalar& ci, @@ -99,17 +89,21 @@ public: Vector& Tstar, Vector& vstar) const override; - - void gravitationalAcceleration(const Vector g); - Vector gravitationalAcceleration() const; + Scalar Cl() const; + void Cl(const Scalar x); + + Scalar Cq() const; + void Cq(const Scalar x); private: - Vector mGravitationalAcceleration; + Scalar mCl; + Scalar mCq; + }; } -#include "GHLLCInline.hh" +#include "SecondOrderArtificialViscosityInline.hh" #endif diff --git a/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInline.hh b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInline.hh new file mode 100644 index 000000000..baa102bd7 --- /dev/null +++ b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInline.hh @@ -0,0 +1,40 @@ +namespace Spheral { + +//------------------------------------------------------------------------------ +// set/get linear reconstruction switch +//------------------------------------------------------------------------------ +template +inline +typename Dimension::Scalar +SecondOrderArtificialViscosity:: +Cl() const { + return mCl; +} + +template +inline +void +SecondOrderArtificialViscosity:: +Cl(typename Dimension::Scalar x) { + mCl=x; +} + + +template +inline +typename Dimension::Scalar +SecondOrderArtificialViscosity:: +Cq() const { + return mCq; +} + +template +inline +void +SecondOrderArtificialViscosity:: +Cq(typename Dimension::Scalar x) { + mCq=x; +} + + +} \ No newline at end of file diff --git a/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInst.cc.py b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInst.cc.py new file mode 100644 index 000000000..4aebeff1f --- /dev/null +++ b/src/GSPH/RiemannSolvers/SecondOrderArtificialViscosityInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "Geometry/Dimension.hh" +#include "GSPH/RiemannSolvers/SecondOrderArtificialViscosity.cc" + +namespace Spheral { + template class SecondOrderArtificialViscosity >; +} +""" diff --git a/src/GSPH/computeMFMDensity.cc b/src/GSPH/computeMFMDensity.cc index 41e1c0a6c..c617172aa 100644 --- a/src/GSPH/computeMFMDensity.cc +++ b/src/GSPH/computeMFMDensity.cc @@ -1,5 +1,7 @@ //---------------------------------Spheral++----------------------------------// // Compute the density from m/V +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #include "GSPH/computeMFMDensity.hh" diff --git a/src/GSPH/computeMFMDensity.hh b/src/GSPH/computeMFMDensity.hh index ebfc8fae2..99f25822d 100644 --- a/src/GSPH/computeMFMDensity.hh +++ b/src/GSPH/computeMFMDensity.hh @@ -1,5 +1,7 @@ //---------------------------------Spheral++----------------------------------// // Compute the density from m/V +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #ifndef __Spheral__computeMFMDensity__ diff --git a/src/GSPH/computeSPHVolume.cc b/src/GSPH/computeSPHVolume.cc index 0fba6577a..ce803c322 100644 --- a/src/GSPH/computeSPHVolume.cc +++ b/src/GSPH/computeSPHVolume.cc @@ -1,5 +1,7 @@ //---------------------------------Spheral++----------------------------------// // Compute the volume from m/rho +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #include "GSPH/computeSPHVolume.hh" diff --git a/src/GSPH/computeSPHVolume.hh b/src/GSPH/computeSPHVolume.hh index c0c856faa..536101c43 100644 --- a/src/GSPH/computeSPHVolume.hh +++ b/src/GSPH/computeSPHVolume.hh @@ -1,5 +1,7 @@ //---------------------------------Spheral++----------------------------------// // Compute the volume from m/rho +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #ifndef __Spheral__computeSPHVolume__ diff --git a/src/GSPH/computeSumVolume.cc b/src/GSPH/computeSumVolume.cc index 60e527051..5dd4126dc 100644 --- a/src/GSPH/computeSumVolume.cc +++ b/src/GSPH/computeSumVolume.cc @@ -1,5 +1,10 @@ //---------------------------------Spheral++----------------------------------// -// Compute volume from inverse of the kernel summation +// Compute volume from inverse of the kernel summation. +// +// Hopkins P.F. (2015) "A New Class of Accurate, Mesh-Free Hydrodynamic +// Simulation Methods," MNRAS, 450(1):53-110 +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #include "GSPH/computeSumVolume.hh" diff --git a/src/GSPH/computeSumVolume.hh b/src/GSPH/computeSumVolume.hh index 1eec78308..283a3e19f 100644 --- a/src/GSPH/computeSumVolume.hh +++ b/src/GSPH/computeSumVolume.hh @@ -1,5 +1,10 @@ //---------------------------------Spheral++----------------------------------// -// Compute volume from inverse of the kernel summation +// Compute volume from inverse of the kernel summation. +// +// Hopkins P.F. (2015) "A New Class of Accurate, Mesh-Free Hydrodynamic +// Simulation Methods," MNRAS, 450(1):53-110 +// +// J.M. Pearl 2022 //----------------------------------------------------------------------------// #ifndef __Spheral__computeSumVolume__ diff --git a/src/GSPH/initializeGradients.cc b/src/GSPH/initializeGradients.cc new file mode 100644 index 000000000..baced1087 --- /dev/null +++ b/src/GSPH/initializeGradients.cc @@ -0,0 +1,155 @@ +//---------------------------------Spheral++----------------------------------// +// initializes the pressure and velocity gradients for Riemann solver - based +// SPH varients +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "GSPH/initializeGradients.hh" +#include "Field/FieldList.hh" +#include "Neighbor/ConnectivityMap.hh" +#include "Kernel/TableKernel.hh" +#include "NodeList/NodeList.hh" + +namespace Spheral{ + +template +void +initializeGradients(const ConnectivityMap& connectivityMap, + const TableKernel& W, + const FieldList& position, + const FieldList& H, + const FieldList& volume, + const FieldList& pressure, + const FieldList& velocity, + FieldList& M, + FieldList& DpDx, + FieldList& DvDx) { + + typedef typename Dimension::Tensor Tensor; + + const auto& nodeLists = connectivityMap.nodeLists(); + const auto& pairs = connectivityMap.nodePairList(); + const auto npairs = pairs.size(); + const auto numNodeLists = nodeLists.size(); + + REQUIRE(volume.size() == numNodeLists); + REQUIRE(velocity.size() == numNodeLists); + REQUIRE(pressure.size() == numNodeLists); + REQUIRE(position.size() == numNodeLists); + REQUIRE(H.size() == numNodeLists); + + REQUIRE(DpDx.size() == numNodeLists); + REQUIRE(DvDx.size() == numNodeLists); + REQUIRE(M.size() == numNodeLists); + +#pragma omp parallel + { + // Thread private scratch variables + int i, j, nodeListi, nodeListj; + + typename SpheralThreads::FieldListStack threadStack; + auto M_thread = M.threadCopy(threadStack); + auto DpDx_thread = DpDx.threadCopy(threadStack); + auto DvDx_thread = DvDx.threadCopy(threadStack); + +#pragma omp for + for (auto kk = 0u; kk < npairs; ++kk) { + i = pairs[kk].i_node; + j = pairs[kk].j_node; + nodeListi = pairs[kk].i_list; + nodeListj = pairs[kk].j_list; + + // Get the state for node i. + const auto& vi = velocity(nodeListi, i); + const auto& Pi = pressure(nodeListi, i); + const auto& ri = position(nodeListi, i); + const auto& voli = volume(nodeListi, i); + const auto& Hi = H(nodeListi, i); + const auto Hdeti = Hi.Determinant(); + + CHECK(voli > 0.0); + CHECK(Hdeti > 0.0); + + auto& DpDxi = DpDx(nodeListi, i); + auto& DvDxi = DvDx(nodeListi, i); + auto& Mi = M(nodeListi, i); + + // Get the state for node j + const auto& vj = velocity(nodeListj, j); + const auto& Pj = pressure(nodeListj, j); + const auto& rj = position(nodeListj, j); + const auto& volj = volume(nodeListj, j); + const auto& Hj = H(nodeListj, j); + const auto Hdetj = Hj.Determinant(); + + CHECK(volj > 0.0); + CHECK(Hdetj > 0.0); + + auto& DpDxj = DpDx(nodeListj, j); + auto& DvDxj = DvDx(nodeListj, j); + auto& Mj = M(nodeListj, j); + + const auto rij = ri - rj; + const auto vij = vi - vj; + const auto Pij = Pi - Pj; + + const auto etai = Hi*rij; + const auto etaj = Hj*rij; + const auto etaMagi = etai.magnitude(); + const auto etaMagj = etaj.magnitude(); + + CHECK(etaMagi >= 0.0); + CHECK(etaMagj >= 0.0); + + const auto gWi = W.gradValue(etaMagi, Hdeti); + const auto Hetai = Hi*etai.unitVector(); + const auto gradWi = gWi*Hetai; + + const auto gWj = W.gradValue(etaMagj, Hdetj); + const auto Hetaj = Hj*etaj.unitVector(); + const auto gradWj = gWj*Hetaj; + + const auto gradPsii = volj*gradWi; + const auto gradPsij = voli*gradWj; + + Mi -= rij.dyad(gradPsii); + Mj -= rij.dyad(gradPsij); + + DpDxi -= Pij*gradPsii; + DpDxj -= Pij*gradPsij; + + DvDxi -= vij.dyad(gradPsii); + DvDxj -= vij.dyad(gradPsij); + + } // loop over pairs + // Reduce the thread values to the master. + threadReduceFieldLists(threadStack); + } // OpenMP parallel region + + // Finish up the spatial gradient calculation + for (auto nodeListi = 0u; nodeListi < numNodeLists; ++nodeListi) { + const auto& nodeList = volume[nodeListi]->nodeList(); + const auto ni = nodeList.numInternalNodes(); +#pragma omp parallel for + for (auto i = 0u; i < ni; ++i) { + const auto numNeighborsi = connectivityMap.numNeighborsForNode(nodeListi, i); + + auto& Mi = M(nodeListi, i); + auto& DpDxi = DpDx(nodeListi, i); + auto& DvDxi = DvDx(nodeListi, i); + + const auto Mdeti = std::abs(Mi.Determinant()); + + const auto enoughNeighbors = numNeighborsi > Dimension::pownu(2); + const auto goodM = (Mdeti > 1e-2 and enoughNeighbors); + + Mi = ( goodM ? Mi.Inverse() : Tensor::one); + + DpDxi = Mi.Transpose()*DpDxi; + DvDxi = DvDxi*Mi; + + } // loop nodes + } // loop nodelists +} // function +} // spheral namespace diff --git a/src/GSPH/initializeGradients.hh b/src/GSPH/initializeGradients.hh new file mode 100644 index 000000000..8930dd563 --- /dev/null +++ b/src/GSPH/initializeGradients.hh @@ -0,0 +1,37 @@ +//---------------------------------Spheral++----------------------------------// +// initializes the pressure and velocity gradients for Riemann solver - based +// SPH varients +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#ifndef __Spheral__initializeGradients__ +#define __Spheral__initializeGradients__ + +#include + +namespace Spheral { + + // Forward declarations. + template class ConnectivityMap; + template class TableKernel; + template class FieldList; + + +template +void +initializeGradients(const ConnectivityMap& connectivityMap, + const TableKernel& W, + const FieldList& position, + const FieldList& H, + const FieldList& volume, + const FieldList& pressure, + const FieldList& velocity, + FieldList& M, + FieldList& DpDx, + FieldList& DvDx); + +} + + + #endif \ No newline at end of file diff --git a/src/GSPH/initializeGradientsInst.cc.py b/src/GSPH/initializeGradientsInst.cc.py new file mode 100644 index 000000000..d27852d9d --- /dev/null +++ b/src/GSPH/initializeGradientsInst.cc.py @@ -0,0 +1,22 @@ +text = """ +//------------------------------------------------------------------------------ +// Explicit instantiation. +//------------------------------------------------------------------------------ +#include "GSPH/initializeGradients.cc" +#include "Geometry/Dimension.hh" + +namespace Spheral { + + + template void initializeGradients(const ConnectivityMap >&, + const TableKernel >&, + const FieldList, Dim< %(ndim)s >::Vector>&, + const FieldList, Dim< %(ndim)s >::SymTensor>&, + const FieldList, Dim< %(ndim)s >::Scalar>&, + const FieldList, Dim< %(ndim)s >::Scalar>&, + const FieldList, Dim< %(ndim)s >::Vector>&, + FieldList, Dim< %(ndim)s >::Tensor>&, + FieldList, Dim< %(ndim)s >::Vector>&, + FieldList, Dim< %(ndim)s >::Tensor>&); +} +""" diff --git a/src/Hydro/CompatibleDifferenceSpecificThermalEnergyPolicy.cc b/src/Hydro/CompatibleDifferenceSpecificThermalEnergyPolicy.cc index 511390859..154900215 100644 --- a/src/Hydro/CompatibleDifferenceSpecificThermalEnergyPolicy.cc +++ b/src/Hydro/CompatibleDifferenceSpecificThermalEnergyPolicy.cc @@ -1,12 +1,12 @@ //---------------------------------Spheral++----------------------------------// -// CompatibleDifferenceSpecificThermalEnergyPolicy -- An implementation of UpdatePolicyBase -// specialized for the updating the specific thermal energy as a dependent -// quantity. +// CompatibleDifferenceSpecificThermalEnergyPolicy -- An implementation of +// UpdatePolicyBase specialized for the updating the specific thermal energy +// as a dependent quantity. // // This version is specialized for materials with different properties. A // compatible energy discretization in which pairwise work allows for opposite // sign pair-wise work. DepsDti and DepsDtj are used as weights and the -// difference between the conservative and consistent formulations is added +// difference between the conservative and consistent formulations is added // back in. //----------------------------------------------------------------------------// #include "Hydro/CompatibleDifferenceSpecificThermalEnergyPolicy.hh" diff --git a/src/NodeGenerators/CubicNodeGenerator.py b/src/NodeGenerators/CubicNodeGenerator.py index b7cca64ae..56edacbfb 100644 --- a/src/NodeGenerators/CubicNodeGenerator.py +++ b/src/NodeGenerators/CubicNodeGenerator.py @@ -59,9 +59,9 @@ def __init__(self, assert nxdomains * nydomains == mpi.procs # The number of nodes per domain. - nxperdomain = nx / nxdomains + nxperdomain = nx // nxdomains nxremainder = nx % nxdomains - nyperdomain = ny / nydomains + nyperdomain = ny // nydomains nyremainder = ny % nydomains assert nxremainder < nxdomains assert nyremainder < nydomains @@ -84,7 +84,7 @@ def __init__(self, # Compute our domain indicies. ixdomain = mpi.rank % nxdomains - iydomain = mpi.rank / nxdomains + iydomain = mpi.rank // nxdomains ixmin = nodeindex(ixdomain, nxperdomain, nxremainder) ixmax = nodeindex(ixdomain + 1, nxperdomain, nxremainder) iymin = nodeindex(iydomain, nyperdomain, nyremainder) diff --git a/src/PYB11/GSPH/GSPHHydroBase.py b/src/PYB11/GSPH/GSPHHydroBase.py index 49c33231e..a8e8c2806 100644 --- a/src/PYB11/GSPH/GSPHHydroBase.py +++ b/src/PYB11/GSPH/GSPHHydroBase.py @@ -7,6 +7,7 @@ @PYB11template("Dimension") @PYB11module("SpheralGSPH") +@PYB11dynamic_attr class GSPHHydroBase(GenericRiemannHydro): PYB11typedefs = """ diff --git a/src/PYB11/GSPH/GSPH_PYB11.py b/src/PYB11/GSPH/GSPH_PYB11.py index c465528a4..f753cae9b 100644 --- a/src/PYB11/GSPH/GSPH_PYB11.py +++ b/src/PYB11/GSPH/GSPH_PYB11.py @@ -12,6 +12,7 @@ from GenericRiemannHydro import * from GSPHHydroBase import * from MFMHydroBase import * +from MFVHydroBase import * from WaveSpeeds import * from Limiters import * from RiemannSolvers import * @@ -22,6 +23,7 @@ PYB11includes += ['"GSPH/GenericRiemannHydro.hh"', '"GSPH/GSPHHydroBase.hh"', '"GSPH/MFMHydroBase.hh"', + '"GSPH/MFVHydroBase.hh"', '"GSPH/WaveSpeeds/WaveSpeedBase.hh"', '"GSPH/WaveSpeeds/AcousticWaveSpeed.hh"', '"GSPH/WaveSpeeds/DavisWaveSpeed.hh"', @@ -32,9 +34,10 @@ '"GSPH/Limiters/VanAlbaLimiter.hh"', '"GSPH/Limiters/SuperbeeLimiter.hh"', '"GSPH/Limiters/OspreLimiter.hh"', + '"GSPH/Limiters/BarthJespersenLimiter.hh"', '"GSPH/RiemannSolvers/RiemannSolverBase.hh"', '"GSPH/RiemannSolvers/HLLC.hh"', - '"GSPH/RiemannSolvers/GHLLC.hh"', + '"GSPH/RiemannSolvers/SecondOrderArtificialViscosity.hh"', '"FileIO/FileIO.hh"'] #------------------------------------------------------------------------------- @@ -49,7 +52,15 @@ "HydroAccelerationGradient", "SPHGradient", "MixedMethodGradient", - "SPHSameTimeGradient"), export_values = True) + "SPHSameTimeGradient", + "SPHUncorrectedGradient", + "NoGradient"), export_values = True) + +NodeMotionType = PYB11enum(("Lagrangian", + "Eulerian", + "Fician", + "XSPH", + "BackgroundPressure"), export_values = False) #------------------------------------------------------------------------------- # Instantiate our types @@ -59,6 +70,7 @@ GenericRiemannHydro%(ndim)id = PYB11TemplateClass(GenericRiemannHydro, template_parameters="%(Dimension)s") GSPHHydroBase%(ndim)id = PYB11TemplateClass(GSPHHydroBase, template_parameters="%(Dimension)s") MFMHydroBase%(ndim)id = PYB11TemplateClass(MFMHydroBase, template_parameters="%(Dimension)s") +MFVHydroBase%(ndim)id = PYB11TemplateClass(MFVHydroBase, template_parameters="%(Dimension)s") WaveSpeedBase%(ndim)id = PYB11TemplateClass(WaveSpeedBase, template_parameters="%(Dimension)s") AcousticWaveSpeed%(ndim)id = PYB11TemplateClass(AcousticWaveSpeed, template_parameters="%(Dimension)s") DavisWaveSpeed%(ndim)id = PYB11TemplateClass(DavisWaveSpeed, template_parameters="%(Dimension)s") @@ -69,9 +81,10 @@ VanAlbaLimiter%(ndim)id = PYB11TemplateClass(VanAlbaLimiter, template_parameters="%(Dimension)s") SuperbeeLimiter%(ndim)id = PYB11TemplateClass(SuperbeeLimiter, template_parameters="%(Dimension)s") OspreLimiter%(ndim)id = PYB11TemplateClass(OspreLimiter, template_parameters="%(Dimension)s") +BarthJespersenLimiter%(ndim)id = PYB11TemplateClass(BarthJespersenLimiter, template_parameters="%(Dimension)s") RiemannSolverBase%(ndim)id = PYB11TemplateClass(RiemannSolverBase, template_parameters="%(Dimension)s") HLLC%(ndim)id = PYB11TemplateClass(HLLC, template_parameters="%(Dimension)s") -GHLLC%(ndim)id = PYB11TemplateClass(GHLLC, template_parameters="%(Dimension)s") +SecondOrderArtificialViscosity%(ndim)id = PYB11TemplateClass(SecondOrderArtificialViscosity, template_parameters="%(Dimension)s") ''' % {"ndim" : ndim, "Dimension" : "Dim<" + str(ndim) + ">"}) diff --git a/src/PYB11/GSPH/GenericRiemannHydro.py b/src/PYB11/GSPH/GenericRiemannHydro.py index 6db20d7e5..57dcfd80d 100644 --- a/src/PYB11/GSPH/GenericRiemannHydro.py +++ b/src/PYB11/GSPH/GenericRiemannHydro.py @@ -7,6 +7,7 @@ @PYB11template("Dimension") @PYB11module("SpheralGSPH") +@PYB11dynamic_attr class GenericRiemannHydro(Physics): PYB11typedefs = """ @@ -146,7 +147,7 @@ def enforceBoundaries(state = "State<%(Dimension)s>&", cfl = PYB11property("Scalar", "cfl", "cfl", doc="The Courant-Friedrichs-Lewy timestep limit multiplier") specificThermalEnergyDiffusionCoefficient = PYB11property("Scalar", "specificThermalEnergyDiffusionCoefficient", "specificThermalEnergyDiffusionCoefficient", doc="coefficient used to diffuse specificThermalEnergy amongst like nodes.") - riemannSolver = PYB11property("RiemannSolverBase<%(Dimension)s>&", "riemannSolver",returnpolicy="reference_internal",doc="The object defining the interface state construction.") + riemannSolver = PYB11property("RiemannSolverBase<%(Dimension)s>&", "riemannSolver", doc="The object defining the interface state construction.") kernel = PYB11property("const TableKernel<%(Dimension)s>&", "kernel", doc="The interpolation kernel") gradientType = PYB11property("GradientType", "gradientType", "gradientType", doc="Enum to selecting different gradients we can use") @@ -178,6 +179,7 @@ def enforceBoundaries(state = "State<%(Dimension)s>&", timeStepMask = PYB11property("const FieldList<%(Dimension)s, int>&", "timeStepMask", returnpolicy="reference_internal") pressure = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "pressure", returnpolicy="reference_internal") + volume = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "volume", returnpolicy="reference_internal") soundSpeed = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "soundSpeed", returnpolicy="reference_internal") Hideal = PYB11property("const FieldList<%(Dimension)s, SymTensor>&","Hideal", returnpolicy="reference_internal") normalization = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "normalization", returnpolicy="reference_internal") diff --git a/src/PYB11/GSPH/Limiters.py b/src/PYB11/GSPH/Limiters.py index d2cb8d647..b3e85d8f8 100644 --- a/src/PYB11/GSPH/Limiters.py +++ b/src/PYB11/GSPH/Limiters.py @@ -177,4 +177,32 @@ def fluxLimiter(self, "slope limiter from flux limiter." return "Scalar" +#------------------------------------------------------------------------------- +# Barth Jespersen limiter +#------------------------------------------------------------------------------- +@PYB11template("Dimension") +@PYB11module("SpheralGSPH") +class BarthJespersenLimiter(LimiterBase): + + PYB11typedefs = """ + typedef typename %(Dimension)s::Scalar Scalar; + """ + + def pyinit(): + "Barth Jespersen slope limiter constructor" + + @PYB11virtual + @PYB11const + def slopeLimiter(self, + x = "const Scalar"): + "slope limiter from flux limiter." + return "Scalar" + + @PYB11virtual + @PYB11const + def fluxLimiter(self, + x = "const Scalar"): + "slope limiter from flux limiter." + return "Scalar" + diff --git a/src/PYB11/GSPH/MFMHydroBase.py b/src/PYB11/GSPH/MFMHydroBase.py index b713433ee..50e413746 100644 --- a/src/PYB11/GSPH/MFMHydroBase.py +++ b/src/PYB11/GSPH/MFMHydroBase.py @@ -7,6 +7,7 @@ @PYB11template("Dimension") @PYB11module("SpheralGSPH") +@PYB11dynamic_attr class MFMHydroBase(GenericRiemannHydro): PYB11typedefs = """ diff --git a/src/PYB11/GSPH/MFVHydroBase.py b/src/PYB11/GSPH/MFVHydroBase.py new file mode 100644 index 000000000..61b857f06 --- /dev/null +++ b/src/PYB11/GSPH/MFVHydroBase.py @@ -0,0 +1,120 @@ +#------------------------------------------------------------------------------- +# GSPHHydroBase +#------------------------------------------------------------------------------- +from PYB11Generator import * +from GenericRiemannHydro import * +from RestartMethods import * + +@PYB11template("Dimension") +@PYB11module("SpheralGSPH") +@PYB11dynamic_attr +class MFVHydroBase(GenericRiemannHydro): + + PYB11typedefs = """ + typedef typename %(Dimension)s::Scalar Scalar; + typedef typename %(Dimension)s::Vector Vector; + typedef typename %(Dimension)s::Tensor Tensor; + typedef typename %(Dimension)s::SymTensor SymTensor; + typedef typename Physics<%(Dimension)s>::TimeStepType TimeStepType; +""" + + def pyinit(smoothingScaleMethod = "const SmoothingScaleBase<%(Dimension)s>&", + dataBase = "DataBase<%(Dimension)s>&", + riemannSolver = "RiemannSolverBase<%(Dimension)s>&", + W = "const TableKernel<%(Dimension)s>&", + epsDiffusionCoeff = "const Scalar", + cfl = "const double", + useVelocityMagnitudeForDt = "const bool", + compatibleEnergyEvolution = "const bool", + evolveTotalEnergy = "const bool", + XSPH = "const bool", + correctVelocityGradient = "const bool", + nodeMotionCoefficient = "const double", + nodeMotionType = "const NodeMotionType", + gradType = "const GradientType", + densityUpdate = "const MassDensityType", + HUpdate = "const HEvolutionType", + epsTensile = "const double", + nTensile = "const double", + xmin = "const Vector&", + xmax = "const Vector&"): + "GSPHHydroBase constructor" + + #........................................................................... + # Virtual methods + + @PYB11virtual + def initializeProblemStartup(dataBase = "DataBase<%(Dimension)s>&"): + "Tasks we do once on problem startup." + return "void" + + @PYB11virtual + def registerState(dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&"): + "Register the state Hydro expects to use and evolve." + return "void" + + @PYB11virtual + def registerDerivatives(dataBase = "DataBase<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Register the derivatives/change fields for updating state." + return "void" + + @PYB11virtual + def preStepInitialize(self, + dataBase = "const DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Optional hook to be called at the beginning of a time step." + return "void" + + @PYB11virtual + def initialize(time = "const Scalar", + dt = "const Scalar", + dataBase = "const DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Initialize the Hydro before we start a derivative evaluation." + return "void" + + @PYB11virtual + @PYB11const + def evaluateDerivatives(time = "const Scalar", + dt = "const Scalar", + dataBase = "const DataBase<%(Dimension)s>&", + state = "const State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + """Evaluate the derivatives for the principle hydro +mass density, velocity, and specific thermal energy.""" + return "void" + + @PYB11virtual + @PYB11const + def finalizeDerivatives(time = "const Scalar", + dt = "const Scalar", + dataBase = "const DataBase<%(Dimension)s>&", + state = "const State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Finalize the derivatives." + return "void" + + @PYB11virtual + def applyGhostBoundaries(state = "State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Apply boundary conditions to the physics specific fields." + return "void" + + @PYB11virtual + def enforceBoundaries(state = "State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Enforce boundary conditions for the physics specific fields." + return "void" + + DvolumeDt = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "DvolumeDt", returnpolicy="reference_internal") + nodeMotionCoefficient = PYB11property("double", "nodeMotionCoefficient", "nodeMotionCoefficient",doc="multiplier for XSPH and Fician node motion schemes.") + nodeMotionType = PYB11property("NodeMotionType","nodeMotionType","nodeMotionType") + +#------------------------------------------------------------------------------- +# Inject methods +#------------------------------------------------------------------------------- +PYB11inject(RestartMethods, MFVHydroBase) diff --git a/src/PYB11/GSPH/RiemannSolvers.py b/src/PYB11/GSPH/RiemannSolvers.py index 3c1c4a1db..1f3206327 100644 --- a/src/PYB11/GSPH/RiemannSolvers.py +++ b/src/PYB11/GSPH/RiemannSolvers.py @@ -25,11 +25,6 @@ def pyinit(slopeLimiter = "LimiterBase<%(Dimension)s>&", waveSpeed = PYB11property("WaveSpeedBase<%(Dimension)s>&", "waveSpeed",returnpolicy="reference_internal", doc="wave speed object") limiter = PYB11property("LimiterBase<%(Dimension)s>&", "limiter",returnpolicy="reference_internal", doc="slope limiter object") - #DpDx = PYB11property("const FieldList<%(Dimension)s, Vector>&", "DpDx",returnpolicy="reference_internal") - #DvDx = PYB11property("const FieldList<%(Dimension)s, Tensor>&", "DvDx",returnpolicy="reference_internal") - -#PYB11inject(RiemannSolverBaseAbstractMethods, RiemannSolverBase, pure_virtual=True) - #------------------------------------------------------------------------------- # HLLC Approximate Riemann Solver #------------------------------------------------------------------------------- @@ -47,10 +42,11 @@ def pyinit(slopeLimiter = "LimiterBase<%(Dimension)s>&", linearReconstruction = "const bool"): "slope limiter constructor" + #------------------------------------------------------------------------------- -# HLLC Approximate Riemann Solver with constant grav acceleration +# HLLC Approximate Riemann Solver #------------------------------------------------------------------------------- -class GHLLC(HLLC): +class SecondOrderArtificialViscosity(RiemannSolverBase): PYB11typedefs = """ typedef typename %(Dimension)s::Scalar Scalar; @@ -59,11 +55,12 @@ class GHLLC(HLLC): typedef typename %(Dimension)s::SymTensor SymTensor; """ - def pyinit(slopeLimiter = "LimiterBase<%(Dimension)s>&", + def pyinit(Cl = "const Scalar", + Cq = "const Scalar", + slopeLimiter = "LimiterBase<%(Dimension)s>&", waveSpeed = "WaveSpeedBase<%(Dimension)s>&", - linearReconstruction = "const bool", - gravitationalAcceleration = "const Vector"): + linearReconstruction = "const bool"): "slope limiter constructor" - - gravitationalAcceleration = PYB11property("Vector", "gravitationalAcceleration", "gravitationalAcceleration", doc="constant gravitational acceleration vector") + Cl = PYB11property("Scalar", "Cl", "Cl", doc="linear artificial viscosity coefficient") + Cq = PYB11property("Scalar", "Cq", "Cq", doc="quadratic artificial viscosity coefficient") diff --git a/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d.py b/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d.py index f1e064d36..51d7ae2ed 100644 --- a/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d.py +++ b/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d.py @@ -8,7 +8,7 @@ from math import * from Spheral2d import * from SpheralTestUtilities import * -from SpheralGnuPlotUtilities import * +#from SpheralGnuPlotUtilities import * from findLastRestart import * from GenerateNodeDistribution2d import * from CompositeNodeDistribution import * @@ -28,9 +28,9 @@ # Generic problem parameters #------------------------------------------------------------------------------- commandLine(nx1 = 100, - ny1 = 50, + ny1 = 50, nx2 = 100, - ny2 = 50, + ny2 = 50, rho1 = 2.0, rho2 = 1.0, @@ -51,9 +51,9 @@ # kernel HUpdate = IdealH, - nPerh = 1.51, - KernelConstructor = BSplineKernel, - order = 5, + nPerh = 3.0, + KernelConstructor = WendlandC2Kernel, + order = 3, hmin = 0.0001, hmax = 0.5, hminratio = 0.1, @@ -64,6 +64,8 @@ crksph = False, fsisph = False, gsph = False, + mfm = False, + mfv = False, # hydro options solid = False, # fluid limit of the solid hydro @@ -78,6 +80,7 @@ evolveTotalEnergy = False, # integrate total instead of specific energy gradhCorrection = True, # correct for temporal variation in h correctVelocityGradient = True, # linear exact velocity gradient (M correction) (corrected kernesl for GSPH and FSISPH) + # SVPH parameters fcentroidal = 0.0, fcellPressure = 0.0, @@ -90,14 +93,17 @@ fsiInterfaceMethod = HLLCInterface, # (HLLCInterface, ModulusInterface) fsiKernelMethod = NeverAverageKernels, # (NeverAverageKernels, AlwaysAverageKernels, AverageInterfaceKernels) - # GSPH parameters + # GSPH/MFM/MFV parameters gsphEpsDiffuseCoeff = 0.0, gsphLinearCorrect = True, + LimiterConstructor = VanLeerLimiter, + WaveSpeedConstructor = DavisWaveSpeed, + nodeMotionCoefficient = 0.2, # artificial viscosity + Qconstructor = LimitedMonaghanGingoldViscosity, Cl = 1.0, Cq = 1.0, - Qconstructor = MonaghanGingoldViscosity, linearConsistent = False, boolReduceViscosity = False, nh = 5.0, @@ -114,7 +120,7 @@ arCondAlpha = 0.5, # integrator - cfl = 0.5, + cfl = 0.25, IntegratorConstructor = CheapSynchronousRK2Integrator, goalTime = 2.0, steps = None, @@ -150,21 +156,32 @@ assert not svph assert not (compatibleEnergy and evolveTotalEnergy) -assert sum([fsisph,psph,gsph,crksph,svph])<=1 +assert sum([fsisph,psph,gsph,crksph,svph,mfm])<=1 assert not (fsisph and not solid) +assert not ((mfm or gsph or mfv) and ( boolReduceViscosity)) # Decide on our hydro algorithm. hydroname = 'SPH' +useArtificialViscosity=True + if svph: hydroname = "SVPH" elif crksph: hydroname = "CRK"+hydroname + Qconstructor = LimitedMonaghanGingoldViscosity elif psph: hydroname = "P"+hydroname elif fsisph: hydroname = "FSI"+hydroname elif gsph: hydroname = "G"+hydroname + useArtificialViscosity=False +elif mfm: + hydroname = "MFM" + useArtificialViscosity=False +elif mfv: + hydroname = "MFV" + useArtificialViscosity=False if asph: hydorname = "A"+hydroname if solid: @@ -186,12 +203,6 @@ restartBaseName = os.path.join(restartDir, "KelvinHelmholtz-2d") vizBaseName = "KelvinHelmholtz-2d" -#------------------------------------------------------------------------------- -# CRKSPH Switches to ensure consistency -#------------------------------------------------------------------------------- -if crksph or fsisph: - Qconstructor = LimitedMonaghanGingoldViscosity - #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- @@ -218,14 +229,12 @@ #------------------------------------------------------------------------------- # Interpolation kernels. #------------------------------------------------------------------------------- -if KernelConstructor=="NBSplineKernel": +if KernelConstructor == NBSplineKernel: WT = TableKernel(NBSplineKernel(order), 1000) - WTPi = TableKernel(NBSplineKernel(order), 1000, Qhmult) else: WT = TableKernel(KernelConstructor(), 1000) - WTPi = TableKernel(KernelConstructor(), 1000, Qhmult) output("WT") -output("WTPi") + kernelExtent = WT.kernelExtent #------------------------------------------------------------------------------- @@ -332,7 +341,7 @@ def vy(ri): #------------------------------------------------------------------------------- # Construct the artificial viscosity. #------------------------------------------------------------------------------- -if not gsph: +if useArtificialViscosity: q = Qconstructor(Cl, Cq, linearInExpansion) q.epsilon2 = epsilon2 q.limiter = Qlimiter @@ -351,6 +360,7 @@ def vy(ri): #------------------------------------------------------------------------------- if crksph: hydro = CRKSPH(dataBase = db, + Q=q, cfl = cfl, filter = filter, epsTensile = epsilonTensile, @@ -391,18 +401,54 @@ def vy(ri): ASPH = asph, epsTensile = epsilonTensile) elif gsph: - limiter = VanLeerLimiter() - waveSpeed = DavisWaveSpeed() + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) hydro = GSPH(dataBase = db, riemannSolver = solver, W = WT, cfl=cfl, - specificThermalEnergyDiffusionCoefficient = gsphEpsDiffuseCoeff, compatibleEnergyEvolution = compatibleEnergy, correctVelocityGradient= correctVelocityGradient, evolveTotalEnergy = evolveTotalEnergy, densityUpdate=densityUpdate, + gradientType = SPHSameTimeGradient, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfm: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = MFM(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + gradientType = SPHSameTimeGradient, + densityUpdate=densityUpdate, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfv: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = MFV(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + nodeMotionCoefficient = nodeMotionCoefficient, + nodeMotionType = NodeMotionType.Lagrangian, + gradientType = SPHSameTimeGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, XSPH = xsph, ASPH = asph, epsTensile = epsilonTensile, @@ -433,7 +479,7 @@ def vy(ri): #------------------------------------------------------------------------------- # Construct the MMRV physics object. #------------------------------------------------------------------------------- -if boolReduceViscosity: +if boolReduceViscosity and useArtificialViscosity: evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) packages.append(evolveReducingViscosityMultiplier) @@ -494,6 +540,7 @@ def vy(ri): # import SpheralPointmeshSiloDump # vizMethod = SpheralPointmeshSiloDump.dumpPhysicsState control = SpheralController(integrator, WT, + initializeDerivatives=True, statsStep = statsStep, restartStep = restartStep, restartBaseName = restartBaseName, diff --git a/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d_McNally.py b/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d_McNally.py index 05feab1b9..b0f2ce737 100644 --- a/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d_McNally.py +++ b/tests/functional/Hydro/KelvinHelmholtz/KelvinHelmholtz-2d_McNally.py @@ -64,12 +64,13 @@ psph = False, fsisph = False, gsph = False, + mfm = False, # hydro options solid = False, asph = False, # Just for choosing the H algorithm - useVelocityMagnitudeForDt = False, XSPH = False, + useVelocityMagnitudeForDt = False, epsilonTensile = 0.0, nTensile = 8, filter = 0.0, @@ -80,6 +81,8 @@ evolveTotalEnergy = False, # artificial viscosity + Cl = None, + Cq = None, linearConsistent = False, boolReduceViscosity = False, nh = 5.0, @@ -94,8 +97,6 @@ betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - Cl = None, - Cq = None, linearInExpansion = False, Qlimiter = None, balsaraCorrection = False, @@ -120,6 +121,8 @@ # GSPH parameters gsphEpsDiffuseCoeff = 0.0, gsphLinearCorrect = True, + LimiterConstructor = VanLeerLimiter, + WaveSpeedConstructor = DavisWaveSpeed, ## integrator cfl = 0.5, @@ -164,10 +167,12 @@ assert numNodeLists in (1, 2) assert not svph assert not (compatibleEnergy and evolveTotalEnergy) -assert sum([fsisph,psph,gsph,crksph,svph])<=1 +assert sum([fsisph,psph,gsph,crksph,svph,mfm])<=1 assert not (fsisph and not solid) +assert not ((mfm or gsph) and (boolCullenViscosity or boolReduceViscosity)) # hydro algorithm label +useArtificialViscosity = True if svph: hydroname = "SVPH" elif crksph: @@ -180,11 +185,19 @@ hydroname = "FSISPH" elif gsph: hydroname = "GSPH" + useArtificialViscosity=False +elif mfm: + hydroname = "MFM" + useArtificialViscosity=False else: hydroname = "SPH" + if asph: hydroname = "A" + hydroname +if solid: + hydroname = "solid" + hydroname + dataDir = os.path.join(dataDir, "rho1=%g-rho2=%g" % (rho1, rho2), "vx1=%g-vx2=%g" % (abs(vx1), abs(vx2)), @@ -429,14 +442,29 @@ ASPH = asph, epsTensile = epsilonTensile) elif gsph: - limiter = VanLeerLimiter() - waveSpeed = DavisWaveSpeed() + limiter = LimiterConstructor + waveSpeed = WaveSpeedConstructor + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, + XSPH = XSPH, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfm: + limiter = LimiterConstructor + waveSpeed = WaveSpeedConstructor solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) hydro = GSPH(dataBase = db, riemannSolver = solver, W = WT, cfl=cfl, - specificThermalEnergyDiffusionCoefficient = gsphEpsDiffuseCoeff, compatibleEnergyEvolution = compatibleEnergy, correctVelocityGradient= correctVelocityGradient, evolveTotalEnergy = evolveTotalEnergy, @@ -471,7 +499,7 @@ #------------------------------------------------------------------------------- # Set the artificial viscosity parameters. #------------------------------------------------------------------------------- -if not gsph: +if useArtificialViscosity: q = hydro.Q if not Cl is None: q.Cl = Cl @@ -498,10 +526,10 @@ #------------------------------------------------------------------------------- # Construct the MMRV physics object. #------------------------------------------------------------------------------- -if boolReduceViscosity: +if boolReduceViscosity and useArtificialViscosity: evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: +elif boolCullenViscosity and useArtificialViscosity: evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WT,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) packages.append(evolveCullenViscosityMultiplier) diff --git a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py index 636f4332d..b11576f72 100644 --- a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py +++ b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py @@ -25,24 +25,27 @@ # # W.F. Noh 1987, JCP, 72, 78-120. #------------------------------------------------------------------------------- -import os, shutil +import os, shutil, mpi, sys from math import * + from SolidSpheral2d import * from SpheralTestUtilities import * from GenerateNodeDistribution2d import * from CubicNodeGenerator import GenerateSquareNodeDistribution from CentroidalVoronoiRelaxation import * -import mpi -import DistributeNodes +if mpi.procs > 1: + from VoronoiDistributeNodes import distributeNodes2d + #from PeanoHilbertDistributeNodes import distributeNodes2d +else: + from DistributeNodes import distributeNodes2d title("2-D integrated hydro test -- cylindrical Noh problem") #------------------------------------------------------------------------------- # Generic problem parameters #------------------------------------------------------------------------------- -commandLine(order = 5, - seed = "constantDTheta", +commandLine(seed = "constantDTheta", thetaFactor = 0.5, azimuthalOffsetFraction = 0.0, @@ -50,7 +53,6 @@ nTheta = 50, rmin = 0.0, rmax = 1.0, - nPerh = 2.01, rho0 = 1.0, eps0 = 0.0, smallPressure = False, @@ -60,22 +62,53 @@ gamma = 5.0/3.0, mu = 1.0, + # hydro type (only one!) + svph = False, + crksph = False, # high order conservative formulation of SPH + psph = False, # pressure-based formulation of SPH + fsisph = False, # formulation for multimaterial problems + gsph = False, # godunov SPH + mfm = False, # moving finite mass of Hopkins 2015 + mfv=False, # moving finite volume of Hopkins 2015 + asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. solid = False, # If true, use the fluid limit of the solid hydro option - svph = False, - crksph = False, - psph = False, - fsisph = False, - gsph = False, - mfm = False, - asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. - boolReduceViscosity = False, + # general hydro options + densityUpdate = RigorousSumDensity, # (IntegrateDensity) + evolveTotalEnergy = False, # evolve total rather than specific energy + compatibleEnergy = True, # evolve specific in a energy conserving manner + gradhCorrection = True, # only for SPH, PSPH (correction for time evolution of h) + correctVelocityGradient = True, # linear gradient correction + XSPH = False, # monaghan's xsph -- move w/ averaged velocity + epsilonTensile = 0.0, # coefficient for the tensile correction + nTensile = 8, # exponent for tensile correction + filter = 0.0, + + # PSPH options HopkinsConductivity = False, # For PSPH + + #CRKSPH options + correctionOrder = LinearOrder, + volumeType = RKSumVolume, + + # MFV + nodeMotion = NodeMotionType.Lagrangian, + + # artificial viscosity + Cl = None, + Cq = None, + linearInExpansion = None, + Qlimiter = None, + balsaraCorrection = None, + epsilon2 = None, + + boolReduceViscosity = False, # morris-monaghan reducing AV nhQ = 5.0, nhL = 10.0, aMin = 0.1, aMax = 2.0, - boolCullenViscosity = False, + + boolCullenViscosity = False, # cullen dehnen AV limiter alphMax = 2.0, alphMin = 0.02, betaC = 0.7, @@ -85,22 +118,18 @@ boolHopkinsCorrection = True, linearConsistent = False, - Cl = None, - Cq = None, - linearInExpansion = None, - Qlimiter = None, - balsaraCorrection = None, - epsilon2 = None, + # kernel options + KernelConstructor = NBSplineKernel, #(NBSplineKernel,WendlandC2Kernel,WendlandC4Kernel,WendlandC6Kernel) + nPerh = 2.01, + HUpdate = IdealH, + order = 5, hmin = 0.0001, hmax = 0.5, hminratio = 0.1, - cfl = 0.25, - XSPH = False, - epsilonTensile = 0.0, - nTensile = 8, - filter = 0.0, + # integrator options IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.25, goalTime = 0.6, steps = None, vizCycle = None, @@ -112,19 +141,11 @@ maxSteps = None, statsStep = 10, smoothIters = 0, - HUpdate = IdealH, - correctionOrder = LinearOrder, - volumeType = RKSumVolume, domainIndependent = False, rigorousBoundaries = False, dtverbose = False, - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - evolveTotalEnergy = False, # Only for SPH variants -- evolve total rather than specific energy - compatibleEnergy = True, - gradhCorrection = True, - correctVelocityGradient = True, - + # output options useVoronoiOutput = True, clearDirectories = False, vizDerivs = False, @@ -141,6 +162,7 @@ assert not(boolReduceViscosity and boolCullenViscosity) assert not((gsph or mfm) and (boolReduceViscosity or boolCullenViscosity)) assert not(fsisph and not solid) +assert sum([crksph,psph,fsisph,svph,gsph,mfm,mfv])<=1 assert thetaFactor in (0.5, 1.0, 2.0) theta = thetaFactor * pi @@ -163,10 +185,14 @@ hydroname = os.path.join("CRKSPH", str(correctionOrder), str(volumeType)) +elif fsisph: + hydroname = "FSISPH" elif gsph: hydroname = "GSPH" elif mfm: hydroname = "MFM" +elif mfv: + hydroname = "MFV" elif psph: hydroname = "PSPH" else: @@ -183,6 +209,7 @@ "compatibleEnergy=%s" % compatibleEnergy, "Cullen=%s" % boolCullenViscosity, "filter=%f" % filter, + "%s" % nodeMotion, "nrad=%i_ntheta=%i" % (nRadial, nTheta)) restartDir = os.path.join(dataDir, "restarts") restartBaseName = os.path.join(restartDir, "Noh-cylindrical-2d-%ix%i" % (nRadial, nTheta)) @@ -196,7 +223,6 @@ #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- -import os, sys if mpi.rank == 0: if clearDirectories and os.path.exists(dataDir): shutil.rmtree(dataDir) @@ -214,7 +240,10 @@ #------------------------------------------------------------------------------- # Interpolation kernels. #------------------------------------------------------------------------------- -WT = TableKernel(NBSplineKernel(order), 1000) +if KernelConstructor==NBSplineKernel: + WT = TableKernel(KernelConstructor(order), 1000) +else: + WT = TableKernel(KernelConstructor(), 1000) output("WT") kernelExtent = WT.kernelExtent @@ -250,8 +279,8 @@ generator = GenerateSquareNodeDistribution(nRadial, nTheta, rho0, - xmin, - xmax, + xmin=xmin, + xmax=xmax, nNodePerh = nPerh, SPH = not asph) else: @@ -265,12 +294,6 @@ nNodePerh = nPerh, SPH = not asph) -if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes2d - #from PeanoHilbertDistributeNodes import distributeNodes2d -else: - from DistributeNodes import distributeNodes2d - distributeNodes2d((nodes1, generator)) output("mpi.reduce(nodes1.numInternalNodes, mpi.MIN)") output("mpi.reduce(nodes1.numInternalNodes, mpi.MAX)") @@ -339,7 +362,7 @@ cfl = cfl, interfaceMethod = HLLCModulus, sumDensityNodeLists=[nodes1], - densityStabilizationCoefficient = 0.00, + densityStabilizationCoefficient = 0.1, compatibleEnergyEvolution = compatibleEnergy, evolveTotalEnergy = evolveTotalEnergy, linearCorrectGradients = correctVelocityGradient, @@ -358,7 +381,7 @@ evolveTotalEnergy = evolveTotalEnergy, XSPH = XSPH, ASPH = asph, - gradientType = RiemannGradient, + gradientType = SPHSameTimeGradient, densityUpdate=densityUpdate, HUpdate = HUpdate, epsTensile = epsilonTensile, @@ -377,7 +400,28 @@ evolveTotalEnergy = evolveTotalEnergy, XSPH = XSPH, ASPH = asph, - gradientType = RiemannGradient, + gradientType = HydroAccelerationGradient, + densityUpdate=densityUpdate, + HUpdate = HUpdate, + epsTensile = epsilonTensile, + nTensile = nTensile) + +elif mfv: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True) + hydro = MFV(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + nodeMotionType=nodeMotion, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + ASPH = asph, + gradientType = SPHSameTimeGradient, densityUpdate=densityUpdate, HUpdate = HUpdate, epsTensile = epsilonTensile, @@ -410,7 +454,7 @@ #------------------------------------------------------------------------------- # Set the artificial viscosity parameters. #------------------------------------------------------------------------------- -if not (gsph or mfm): +if not (gsph or mfm or mfv): q = hydro.Q if Cl: q.Cl = Cl @@ -716,7 +760,11 @@ comparisonFile = os.path.join(dataDir, comparisonFile) import filecmp assert filecmp.cmp(outputFile, comparisonFile) + + +Masserror = (control.conserve.massHistory[-1] - control.conserve.massHistory[0])/max(1.0e-30, control.conserve.massHistory[0]) Eerror = (control.conserve.EHistory[-1] - control.conserve.EHistory[0])/max(1.0e-30, control.conserve.EHistory[0]) +print("Total mass error: %g" % Masserror) print("Total energy error: %g" % Eerror) if compatibleEnergy and abs(Eerror) > 1e-13: raise ValueError("Energy error outside allowed bounds.") diff --git a/tests/functional/Hydro/Noh/Noh-planar-1d.py b/tests/functional/Hydro/Noh/Noh-planar-1d.py index c261b5b09..0d2f7d1bd 100644 --- a/tests/functional/Hydro/Noh/Noh-planar-1d.py +++ b/tests/functional/Hydro/Noh/Noh-planar-1d.py @@ -50,17 +50,17 @@ # # GSPH # -#ATS:t500 = test( SELF, "--gsph True --gsphReconstructionGradient=RiemannGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and RiemannGradient -- 1-D (serial)") -#ATS:t501 = testif(t500, SELF, "--gsph True --gsphReconstructionGradient=RiemannGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and RiemannGradient -- 1-D (serial) RESTART CHECK") -#ATS:t502 = test( SELF, "--gsph True --gsphReconstructionGradient=HydroAccelerationGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and and HydroAccelerationGradient -- 1-D (serial)") -#ATS:t503 = testif(t502, SELF, "--gsph True --gsphReconstructionGradient=HydroAccelerationGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and HydroAccelerationGradient -- 1-D (serial) RESTART CHECK") -#ATS:t504 = test( SELF, "--gsph True --gsphReconstructionGradient=SPHGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and SPHGradient -- 1-D (serial)") -#ATS:t505 = testif(t504, SELF, "--gsph True --gsphReconstructionGradient=SPHGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and SPHGradient -- 1-D (serial) RESTART CHECK") +#ATS:t500 = test( SELF, "--gsph True --gsphReconstructionGradient RiemannGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and RiemannGradient -- 1-D (serial)") +#ATS:t501 = testif(t500, SELF, "--gsph True --gsphReconstructionGradient RiemannGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and RiemannGradient -- 1-D (serial) RESTART CHECK") +#ATS:t502 = test( SELF, "--gsph True --gsphReconstructionGradient HydroAccelerationGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and and HydroAccelerationGradient -- 1-D (serial)") +#ATS:t503 = testif(t502, SELF, "--gsph True --gsphReconstructionGradient HydroAccelerationGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and HydroAccelerationGradient -- 1-D (serial) RESTART CHECK") +#ATS:t504 = test( SELF, "--gsph True --gsphReconstructionGradient SPHGradient --graphics None --clearDirectories True --checkError True --restartStep 20", label="Planar Noh problem with GSPH and SPHGradient -- 1-D (serial)") +#ATS:t505 = testif(t504, SELF, "--gsph True --gsphReconstructionGradient SPHGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with GSPH and SPHGradient -- 1-D (serial) RESTART CHECK") # # MFM # -#ATS:t600 = test( SELF, "--mfm True --gsphReconstructionGradient=RiemannGradient --graphics None --clearDirectories True --checkError False --restartStep 20", label="Planar Noh problem with MFM -- 1-D (serial)") -#ATS:t601 = testif(t600, SELF, "--mfm True --gsphReconstructionGradient=RiemannGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with MFM -- 1-D (serial) RESTART CHECK") +#ATS:t600 = test( SELF, "--mfm True --gsphReconstructionGradient RiemannGradient --graphics None --clearDirectories True --checkError False --restartStep 20", label="Planar Noh problem with MFM -- 1-D (serial)") +#ATS:t601 = testif(t600, SELF, "--mfm True --gsphReconstructionGradient RiemannGradient --graphics None --clearDirectories False --checkError False --restartStep 20 --restoreCycle 20 --steps 20 --checkRestart True", label="Planar Noh problem with MFM -- 1-D (serial) RESTART CHECK") import os, shutil, sys from SolidSpheral1d import * diff --git a/tests/functional/Hydro/Noh/Noh-spherical-3d.py b/tests/functional/Hydro/Noh/Noh-spherical-3d.py index d6d58a6a8..43226af1a 100644 --- a/tests/functional/Hydro/Noh/Noh-spherical-3d.py +++ b/tests/functional/Hydro/Noh/Noh-spherical-3d.py @@ -55,6 +55,7 @@ fsisph = False, gsph = False, mfm = False, + mfv = False, asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. boolReduceViscosity = False, @@ -144,8 +145,10 @@ hydroname = "FSISPH" elif gsph: hydroname = "GSPH" -elif gsph: +elif mfm: hydroname = "MFM" +elif mfv: + hydroname = "MFV" else: hydroname = "SPH" if asph: @@ -327,6 +330,23 @@ HUpdate = IdealH, epsTensile = epsilonTensile, nTensile = nTensile) +elif mfv: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,True) + hydro = MFV(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + gradientType = RiemannGradient, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) elif psph: hydro = PSPH(dataBase = db, W = WT, @@ -360,7 +380,7 @@ output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.HEvolution") -if not (gsph or fsisph): +if not (gsph or mfm or mfv or fsisph): output("hydro.PiKernel") if not fsisph: output("hydro.densityUpdate") @@ -370,7 +390,7 @@ #------------------------------------------------------------------------------- # Set the artificial viscosity parameters. #------------------------------------------------------------------------------- -if not (gsph or mfm): +if not (gsph or mfm or mfv): q = hydro.Q if Cl: q.Cl = Cl diff --git a/tests/functional/Hydro/RayleighTaylor/RT-2d.py b/tests/functional/Hydro/RayleighTaylor/RT-2d.py index 150846102..c0203ce2e 100644 --- a/tests/functional/Hydro/RayleighTaylor/RT-2d.py +++ b/tests/functional/Hydro/RayleighTaylor/RT-2d.py @@ -1,7 +1,7 @@ #------------------------------------------------------------------------------- # This is the basic Rayleigh-Taylor Problem #------------------------------------------------------------------------------- -import shutil +import shutil, os, sys, mpi from math import * from Spheral2d import * from SpheralTestUtilities import * @@ -10,13 +10,16 @@ from GenerateNodeDistribution2d import * from CompositeNodeDistribution import * from CentroidalVoronoiRelaxation import * +from HydrostaticReflectingBoundary import HydrostaticReflectingBoundary2d as HydrostaticReflectingBoundary -import mpi -import DistributeNodes +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes2d +else: + from DistributeNodes import distributeNodes2d title("Rayleigh-Taylor test problem in 2D") -class ExponentialDensity: +class ExponentialProfile: def __init__(self, y1, rho0, @@ -26,6 +29,8 @@ def __init__(self, self.alpha = alpha return def __call__(self, r): + #if r.y > 1.0: + # print self.rho0*exp(self.alpha*(r.y - self.y1)) return self.rho0*exp(self.alpha*(r.y - self.y1)) #------------------------------------------------------------------------------- @@ -35,6 +40,9 @@ def __call__(self, r): ny1 = 100, nx2 = 100, ny2 = 100, + refineFactor = 1, + nybound = 8, # number of layers in const node bc + rho0 = 1.0, eps0 = 1.0, x0 = 0.0, @@ -47,53 +55,94 @@ def __call__(self, r): vx1 = 0.0, vx2 = 0.0, freq = 1.0, - alpha = 0.01, # amplitude of displacement + alpha = 0.0025, # amplitude of displacement beta = 5.0, # speed at which displacement decays away from midline - S = 3.0, # density jump at surface - g0 = -2.0, - w0 = 0.1, + S = 10.0, # density jump at surface + g0 = -0.5, + w0 = 0.005, sigma = 0.05/sqrt(2.0), gamma = 5.0/3.0, mu = 1.0, - nPerh = 1.51, + useHydrostaticBoundary = True, + + # kernel options + HUpdate = IdealH, + KernelConstructor = WendlandC2Kernel, + order = 3, + nPerh = 3.0, + hmin = 0.0001, + hmax = 0.5, + hminratio = 0.1, + + # hydro type + svph = False, + crksph = False, + psph = False, + fsisph = False, + gsph = False, + mfm = False, + + # hydro options + asph = False, + xsph = False, + solid = False, + filter = 0.0, + densityUpdate = IntegrateDensity, + compatibleEnergy = True, + evolveTotalEnergy = False, + useVelocityMagnitudeForDt = False, + correctVelocityGradient = False, + epsilonTensile = 0.0, + nTensile = 8, + + # SPH/PSPH options + gradhCorrection = False, - SVPH = False, - CRKSPH = False, - ASPH = False, - SPH = True, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. - filter = 0.0, # CRKSPH filtering - Qconstructor = MonaghanGingoldViscosity, - #Qconstructor = TensorMonaghanGingoldViscosity, - linearConsistent = False, + # svph options fcentroidal = 0.0, fcellPressure = 0.0, - boolReduceViscosity = False, - nh = 5.0, - aMin = 0.1, - aMax = 2.0, - Qhmult = 1.0, + linearConsistent = False, + + # FSISPH parameters + fsiSurfaceCoefficient = 0.00, # adds additional repulsive force to material interfaces) + fsiRhoStabilizeCoeff = 0.1, # coefficient that smooths the density field + fsiEpsDiffuseCoeff = 0.1, # explicit diiffusion of the thermal energy + fsiXSPHCoeff = 0.00, # fsi uses multiplier for XSPH instead of binary switch + fsiInterfaceMethod = ModulusInterface, # (HLLCInterface, ModulusInterface) + fsiKernelMethod = NeverAverageKernels, # (NeverAverageKernels, AlwaysAverageKernels, AverageInterfaceKernels) + + # GSPH/MFM parameters + gsphEpsDiffuseCoeff = 0.0, + gsphLinearCorrect = True, + LimiterConstructor = VanLeerLimiter, + WaveSpeedConstructor = DavisWaveSpeed, + + # artificial viscosity options + Qconstructor = LimitedMonaghanGingoldViscosity, Cl = 1.0, Cq = 1.0, linearInExpansion = False, Qlimiter = False, balsaraCorrection = False, epsilon2 = 1e-2, - hmin = 0.0001, - hmax = 0.5, - hminratio = 0.1, - cfl = 0.5, - useVelocityMagnitudeForDt = False, - XSPH = False, - epsilonTensile = 0.0, - nTensile = 8, + + boolReduceViscosity = False, + nh = 5.0, + aMin = 0.1, + aMax = 2.0, + Qhmult = 1.0, + # artificial conduction options + bArtificialConduction = False, + arCondAlpha = 0.5, + + # integrator & options IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.25, goalTime = 2.0, steps = None, - vizCycle = None, - vizTime = 0.01, dt = 0.0001, dtMin = 1.0e-8, dtMax = 0.1, @@ -101,54 +150,62 @@ def __call__(self, r): maxSteps = None, statsStep = 10, smoothIters = 0, - HUpdate = IdealH, domainIndependent = False, rigorousBoundaries = False, dtverbose = False, - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - compatibleEnergy = True, # <--- Important! rigorousBoundaries does not work with the compatibleEnergy algorithm currently. - gradhCorrection = False, - + # output options + vizCycle = None, + vizTime = 0.1, useVoronoiOutput = False, clearDirectories = False, restoreCycle = None, restartStep = 100, - redistributeStep = 500, + redistributeStep = 50000, checkRestart = False, dataDir = "dumps-Rayleigh-Taylor-2d", outputFile = "None", comparisonFile = "None", serialDump = False, #whether to dump a serial ascii file at the end for viz - - bArtificialConduction = False, - arCondAlpha = 0.5, ) + +assert not (compatibleEnergy and evolveTotalEnergy) +assert sum([fsisph,psph,gsph,crksph,svph,mfm])<=1 +assert not (fsisph and not solid) +assert not ((mfm or gsph) and (boolReduceViscosity)) + # Decide on our hydro algorithm. -if SVPH: - if ASPH: - HydroConstructor = ASVPHFacetedHydro - else: - HydroConstructor = SVPHFacetedHydro -elif CRKSPH: - if ASPH: - HydroConstructor = ACRKSPHHydro - else: - HydroConstructor = CRKSPHHydro -else: - if ASPH: - HydroConstructor = ASPHHydro - else: - HydroConstructor = SPHHydro +hydroname = 'SPH' +useArtificialViscosity=True + +if svph: + hydroname = "SVPH" +elif crksph: + Qconstructor = LimitedMonaghanGingoldViscosity + hydroname = "CRK"+hydroname +elif psph: + hydroname = "P"+hydroname +elif fsisph: + hydroname = "FSI"+hydroname +elif gsph: + hydroname = "G"+hydroname + useArtificialViscosity=False +elif mfm: + hydroname = "MFM" + useArtificialViscosity=False +if asph: + hydorname = "A"+hydroname +if solid: + hydroname = "solid"+hydroname dataDir = os.path.join(dataDir, "S=%g" % (S), "vx1=%g-vx2=%g" % (abs(vx1), abs(vx2)), - str(HydroConstructor).split("'")[1].split(".")[-1], + hydroname, "densityUpdate=%s" % (densityUpdate), - "XSPH=%s" % XSPH, + "XSPH=%s" % xsph, "filter=%s" % filter, "%s-Cl=%g-Cq=%g" % (str(Qconstructor).split("'")[1].split(".")[-1], Cl, Cq), "%ix%i" % (nx1, ny1 + ny2), @@ -158,16 +215,10 @@ def __call__(self, r): restartBaseName = os.path.join(restartDir, "Rayleigh-Taylor-2d") vizBaseName = "Rayleigh-Taylor-2d" -#------------------------------------------------------------------------------- -# CRKSPH Switches to ensure consistency -#------------------------------------------------------------------------------- -if CRKSPH: - Qconstructor = LimitedMonaghanGingoldViscosity #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- -import os, sys if mpi.rank == 0: if clearDirectories and os.path.exists(dataDir): shutil.rmtree(dataDir) @@ -191,19 +242,27 @@ def __call__(self, r): #------------------------------------------------------------------------------- # Interpolation kernels. #------------------------------------------------------------------------------- -WT = TableKernel(BSplineKernel(), 1000) +if KernelConstructor==NBSplineKernel: + WT = TableKernel(NBSplineKernel(order), 1000) +else: + WT = TableKernel(KernelConstructor(), 1000) output("WT") kernelExtent = WT.kernelExtent #------------------------------------------------------------------------------- # Make the NodeList. #------------------------------------------------------------------------------- -nodes1 = makeFluidNodeList("High density gas", eos, +if solid: + nodeListConstructor=makeSolidNodeList +else: + nodeListConstructor=makeFluidNodeList + +nodes1 = nodeListConstructor("Low density gas", eos, hmin = hmin, hmax = hmax, hminratio = hminratio, nPerh = nPerh) -nodes2 = makeFluidNodeList("Low density gas", eos, +nodes2 = nodeListConstructor("High density gas", eos, hmin = hmin, hmax = hmax, hminratio = hminratio, @@ -216,35 +275,47 @@ def __call__(self, r): output("nodes.hminratio") output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# functions for ICs. +#------------------------------------------------------------------------------- +eps1=eps0 +eps2=eps1/S + +lowerDensity = ExponentialProfile(y1, + rho0/S, + g0/((gamma - 1.0)*eps1)) +upperDensity = ExponentialProfile(y1, + rho0, + g0/((gamma - 1.0)*eps2)) + + #------------------------------------------------------------------------------- # Set the node properties. #------------------------------------------------------------------------------- + if restoreCycle is None: - generator1 = GenerateNodeDistribution2d(nx1, ny1, - rho = ExponentialDensity(y1, - rho0/S, - g0/((gamma - 1.0)*eps0)), - distributionType = "lattice", - xmin = (x0,y0), + nx1 *= refineFactor + ny1 *= refineFactor + nx2 *= refineFactor + ny2 *= refineFactor + dy = (y1 - y0)/ny1 + + generator1 = GenerateNodeDistribution2d(nx1, ny1+nybound, + rho = lowerDensity, + distributionType = "xstaggeredLattice", + xmin = (x0,y0-nybound*dy), xmax = (x1,y1), - nNodePerh = nPerh, SPH = SPH) - generator2 = GenerateNodeDistribution2d(nx2, ny2, - rho = ExponentialDensity(y1, - rho0, - g0*S/((gamma - 1.0)*eps0)), - distributionType = "lattice", + generator2 = GenerateNodeDistribution2d(nx2, ny2+nybound, + rho = upperDensity, + distributionType = "xstaggeredLattice", xmin = (x0,y1), - xmax = (x1,y2), + xmax = (x1,y2+nybound*dy), nNodePerh = nPerh, SPH = SPH) - if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes2d - else: - from DistributeNodes import distributeNodes2d - distributeNodes2d((nodes1, generator1), (nodes2, generator2)) @@ -275,64 +346,114 @@ def dy(ri): #------------------------------------------------------------------------------- # 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 useArtificialViscosity: + 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. #------------------------------------------------------------------------------- -if SVPH: - hydro = HydroConstructor(W = WT, - Q = q, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - densityUpdate = densityUpdate, - XSVPH = XSPH, - linearConsistent = linearConsistent, - generateVoid = False, - HUpdate = HUpdate, - fcentroidal = fcentroidal, - fcellPressure = fcellPressure, - xmin = Vector(-2.0, -2.0), - xmax = Vector(3.0, 3.0)) -# xmin = Vector(x0 - 0.5*(x2 - x0), y0 - 0.5*(y2 - y0)), -# xmax = Vector(x2 + 0.5*(x2 - x0), y2 + 0.5*(y2 - y0))) -elif CRKSPH: - hydro = HydroConstructor(W = WT, - Q = q, - filter = filter, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - XSPH = XSPH, - densityUpdate = densityUpdate, - HUpdate = HUpdate) +if crksph: + hydro = CRKSPH(dataBase = db, + cfl = cfl, + filter = filter, + epsTensile = epsilonTensile, + nTensile = nTensile, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = xsph, + densityUpdate = densityUpdate, + HUpdate = HUpdate) +elif psph: + hydro = PSPH(dataBase = db, + cfl = cfl, + W = WT, + Q = q, + filter = filter, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate = densityUpdate, + HUpdate = HUpdate, + XSPH = xsph, + ASPH = asph) +if fsisph: + sumDensityNodeListSwitch =[nodes1,nodes2] + hydro = FSISPH(dataBase = db, + Q=q, + W = WT, + cfl = cfl, + surfaceForceCoefficient = fsiSurfaceCoefficient, + densityStabilizationCoefficient = fsiRhoStabilizeCoeff, + specificThermalEnergyDiffusionCoefficient = fsiEpsDiffuseCoeff, + xsphCoefficient = fsiXSPHCoeff, + interfaceMethod = fsiInterfaceMethod, + kernelAveragingMethod = fsiKernelMethod, + sumDensityNodeLists = sumDensityNodeListSwitch, + correctVelocityGradient = correctVelocityGradient, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + ASPH = asph, + epsTensile = epsilonTensile) +elif gsph: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfm: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = MFM(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) else: - hydro = HydroConstructor(W = WT, - Q = q, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - gradhCorrection = gradhCorrection, - XSPH = XSPH, - densityUpdate = densityUpdate, - HUpdate = HUpdate, - epsTensile = epsilonTensile, - nTensile = nTensile) + hydro = SPH(dataBase = db, + cfl = cfl, + W = WT, + Q = q, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + gradhCorrection = gradhCorrection, + correctVelocityGradient = correctVelocityGradient, + XSPH = xsph, + ASPH = asph, + densityUpdate = densityUpdate, + HUpdate = HUpdate, + epsTensile = epsilonTensile, + nTensile = nTensile) + output("hydro") output("hydro.kernel()") -output("hydro.PiKernel()") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.densityUpdate") @@ -344,9 +465,8 @@ def dy(ri): # Construct the MMRV physics object. #------------------------------------------------------------------------------- -if boolReduceViscosity: +if boolReduceViscosity and useArtificialViscosity: evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) - packages.append(evolveReducingViscosityMultiplier) #------------------------------------------------------------------------------- @@ -354,9 +474,7 @@ def dy(ri): #------------------------------------------------------------------------------- if bArtificialConduction: - #q.reducingViscosityCorrection = True - ArtyCond = ArtificialConduction(WT,arCondAlpha) - + ArtyCond = ArtificialConduction(WT,arCondAlpha) packages.append(ArtyCond) #------------------------------------------------------------------------------- @@ -370,9 +488,6 @@ def dy(ri): for i in range(nodes2.numInternalNodes): nodeIndicies2.append(i) -#nodeIndicies1.extend(range(nodes1.numInternalNodes)) -#nodeIndicies2.extend(range(nodes2.numInternalNodes)) - gravity1 = ConstantAcceleration2d(Vector2d(0.0, g0), nodes1, nodeIndicies1) @@ -391,11 +506,23 @@ def dy(ri): yp1 = Plane(Vector(x0, y0), Vector(0.0, 1.0)) yp2 = Plane(Vector(x0, y2), Vector(0.0, -1.0)) xbc = PeriodicBoundary(xp1, xp2) -#ybc = PeriodicBoundary(yp1, yp2) -ybc1 = ReflectingBoundary(yp1) -ybc2 = ReflectingBoundary(yp2) -bcSet = [xbc, ybc1, ybc2] -#bcSet = [xbc,ybc1] + +pos = nodes1.positions() +ylow, yhigh = vector_of_int(), vector_of_int() +for i in xrange(nodes1.numInternalNodes): + if pos[i].y < y0: + ylow.append(i) + +pos = nodes2.positions() +for i in xrange(nodes2.numInternalNodes): + if pos[i].y > y2: + yhigh.append(i) + +print(yhigh) +ybc1 = ConstantBoundary(db, nodes1, ylow, yp1) +ybc2 = ConstantBoundary(db, nodes2, yhigh, yp2) + +bcSet = [ybc1, ybc2, xbc] for bc in bcSet: for p in packages: @@ -443,6 +570,7 @@ def dy(ri): redistributeStep = redistributeStep, vizMethod = vizMethod, vizBaseName = vizBaseName, + vizGhosts=True, vizDir = vizDir, vizStep = vizCycle, vizTime = vizTime, diff --git a/tests/functional/Hydro/RayleighTaylor/RT-2d_Hopkins.py b/tests/functional/Hydro/RayleighTaylor/RT-2d_Hopkins.py index 5946445ec..3175acd48 100644 --- a/tests/functional/Hydro/RayleighTaylor/RT-2d_Hopkins.py +++ b/tests/functional/Hydro/RayleighTaylor/RT-2d_Hopkins.py @@ -5,7 +5,7 @@ #------------------------------------------------------------------------------- # This is the basic Rayleigh-Taylor Problem #------------------------------------------------------------------------------- -import shutil +import shutil, os, sys from math import * from Spheral2d import * from SpheralTestUtilities import * @@ -18,7 +18,10 @@ from RTMixLength import RTMixLength import mpi -import DistributeNodes +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes2d +else: + from DistributeNodes import distributeNodes2d class ExponentialDensity: def __init__(self, @@ -46,25 +49,72 @@ def __call__(self, r): y0 = 0.0, y1 = 1.0, gval = -0.5, - w0 = 0.025, + w0 = 0.025, delta = 0.025, gamma = 1.4, mu = 1.0, - nPerh = 1.01, - - SVPH = False, - CRKSPH = False, - PSPH = False, - SPH = True, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. - filter = 0.0, # CRKSPH filtering - Qconstructor = MonaghanGingoldViscosity, - #Qconstructor = TensorMonaghanGingoldViscosity, - KernelConstructor = BSplineKernel, + # kernel options + KernelConstructor = WendlandC2Kernel, + HUpdate = IdealH, + nPerh = 3.01, order = 5, + hmin = 0.0001, + hmax = 0.5, + hminratio = 0.1, + + # hydros + svph = False, + crksph = False, + psph = False, + fsisph = False, + gsph = False, + mfm = False, + + # general hydro options + asph = False, + xsph = False, + solid = False, + filter = 0.0, + epsilonTensile = 0.0, + nTensile = 8, + useVelocityMagnitudeForDt = False, + densityUpdate = IntegrateDensity, + compatibleEnergy = True, + correctVelocityGradient = True, + evolveTotalEnergy= False, + + # SPH/PSPH options + gradhCorrection = True, + HopkinsConductivity=False, + + # svph options linearConsistent = False, fcentroidal = 0.0, fcellPressure = 0.0, + + # crksph options + correctionOrder = LinearOrder, + + # FSISPH parameters + fsiSurfaceCoefficient = 0.00, # adds additional repulsive force to material interfaces) + fsiRhoStabilizeCoeff = 0.1, # coefficient that smooths the density field + fsiEpsDiffuseCoeff = 0.1, # explicit diiffusion of the thermal energy + fsiXSPHCoeff = 0.00, # fsi uses multiplier for XSPH instead of binary switch + fsiInterfaceMethod = ModulusInterface, # (HLLCInterface, ModulusInterface) + fsiKernelMethod = NeverAverageKernels, # (NeverAverageKernels, AlwaysAverageKernels, AverageInterfaceKernels) + + # GSPH/MFM parameters + gsphEpsDiffuseCoeff = 0.0, + gsphLinearCorrect = True, + LimiterConstructor = VanLeerLimiter, + WaveSpeedConstructor = DavisWaveSpeed, + riemannGradientType = HydroAccelerationGradient, + + # artificial viscosity + Qconstructor = LimitedMonaghanGingoldViscosity, + Cl = 1.0, + Cq = 1.0, boolReduceViscosity = False, nh = 5.0, aMin = 0.1, @@ -78,22 +128,13 @@ def __call__(self, r): betaE = 1.0, fKern = 1.0/3.0, boolHopkinsCorrection = True, - Cl = 1.0, - Cq = 1.0, linearInExpansion = False, Qlimiter = False, balsaraCorrection = False, epsilon2 = 1e-2, - hmin = 0.0001, - hmax = 0.5, - hminratio = 0.1, - cfl = 0.5, - useVelocityMagnitudeForDt = False, - XSPH = False, - epsilonTensile = 0.0, - nTensile = 8, IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.5, goalTime = 4.0, steps = None, vizCycle = None, @@ -105,18 +146,9 @@ def __call__(self, r): maxSteps = None, statsStep = 10, smoothIters = 0, - HUpdate = IdealH, domainIndependent = False, rigorousBoundaries = False, dtverbose = False, - - correctionOrder = LinearOrder, - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - compatibleEnergy = True, # <--- Important! rigorousBoundaries does not work with the compatibleEnergy algorithm currently. - gradhCorrection = True, - correctVelocityGradient = True, - evolveTotalEnergy= False, - HopkinsConductivity=False, clearDirectories = False, restoreCycle = -1, @@ -134,39 +166,44 @@ def __call__(self, r): arCondAlpha = 0.5, ) +assert not (compatibleEnergy and evolveTotalEnergy) +assert sum([fsisph,psph,gsph,crksph,svph,mfm])<=1 +assert not (fsisph and not solid) +assert not ((mfm or gsph) and (boolReduceViscosity or boolReduceViscosity)) assert not(boolReduceViscosity and boolCullenViscosity) -# Decide on our hydro algorithm. -if SVPH: - if SPH: - HydroConstructor = SVPHFacetedHydro - else: - HydroConstructor = ASVPHFacetedHydro -elif CRKSPH: + +hydroname = 'SPH' +useArtificialViscosity=True + +if svph: + hydroname = "SVPH" +elif crksph: Qconstructor = LimitedMonaghanGingoldViscosity - if SPH: - HydroConstructor = CRKSPHHydro - else: - HydroConstructor = ACRKSPHHydro -elif PSPH: - if SPH: - HydroConstructor = PSPHHydro - else: - HydroConstructor = APSPHHydro -else: - if SPH: - HydroConstructor = SPHHydro - else: - HydroConstructor = ASPHHydro + hydroname = "CRK"+hydroname +elif psph: + hydroname = "P"+hydroname +elif fsisph: + hydroname = "FSI"+hydroname +elif gsph: + hydroname = "G"+hydroname + useArtificialViscosity=False +elif mfm: + hydroname = "MFM" + useArtificialViscosity=False +if asph: + hydorname = "A"+hydroname +if solid: + hydroname = "solid"+hydroname dataDir = os.path.join(dataDir, "gval=%g" % (gval), "w0=%g" % w0, - HydroConstructor.__name__, + hydroname, Qconstructor.__name__, KernelConstructor.__name__, "densityUpdate=%s" % (densityUpdate), "correctionOrder=%s" % (correctionOrder), - "XSPH=%s" % XSPH, + "XSPH=%s" % xsph, "filter=%s" % filter, "compatible=%s" % compatibleEnergy, "Cullen=%s" % boolCullenViscosity, @@ -181,7 +218,6 @@ def __call__(self, r): #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- -import os, sys if mpi.rank == 0: if clearDirectories and os.path.exists(dataDir): shutil.rmtree(dataDir) @@ -201,18 +237,20 @@ def __call__(self, r): #------------------------------------------------------------------------------- if KernelConstructor==NBSplineKernel: WT = TableKernel(NBSplineKernel(order), 10000) - WTPi = TableKernel(NBSplineKernel(order), 10000, Qhmult) else: WT = TableKernel(KernelConstructor(), 10000) - WTPi = TableKernel(KernelConstructor(), 10000, Qhmult) output("WT") -output("WTPi") kernelExtent = WT.kernelExtent #------------------------------------------------------------------------------- # Make the NodeList. #------------------------------------------------------------------------------- -nodes = makeFluidNodeList("High density gas", eos, +if solid: + nodeListConstructor=makeSolidNodeList +else: + nodeListConstructor=makeFluidNodeList + +nodes = nodeListConstructor("High density gas", eos, hmin = hmin, hmax = hmax, hminratio = hminratio, @@ -240,11 +278,6 @@ def __call__(self, r): nNodePerh = nPerh, SPH = SPH) -if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes2d -else: - from DistributeNodes import distributeNodes2d - distributeNodes2d((nodes, generator)) #Set IC @@ -276,78 +309,117 @@ def __call__(self, r): #------------------------------------------------------------------------------- # 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 useArtificialViscosity: + 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. #------------------------------------------------------------------------------- -if SVPH: - hydro = HydroConstructor(W = WT, - Q = q, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - densityUpdate = densityUpdate, - XSVPH = XSPH, - linearConsistent = linearConsistent, - generateVoid = False, - HUpdate = HUpdate, - fcentroidal = fcentroidal, - fcellPressure = fcellPressure, - xmin = Vector(-2.0, -2.0), - xmax = Vector(3.0, 3.0)) -elif CRKSPH: - hydro = HydroConstructor(W = WT, - WPi = WTPi, - Q = q, - filter = filter, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - XSPH = XSPH, - correctionOrder = correctionOrder, - densityUpdate = densityUpdate, - HUpdate = HUpdate) -elif PSPH: - hydro = HydroConstructor(W = WT, - Q = q, - filter = filter, - cfl = cfl, - compatibleEnergyEvolution = compatibleEnergy, - evolveTotalEnergy = evolveTotalEnergy, - HopkinsConductivity = HopkinsConductivity, - densityUpdate = densityUpdate, - correctVelocityGradient = correctVelocityGradient, - HUpdate = HUpdate, - XSPH = XSPH) +if crksph: + hydro = CRKSPH(dataBase = db, + cfl = cfl, + filter = filter, + epsTensile = epsilonTensile, + nTensile = nTensile, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = xsph, + densityUpdate = densityUpdate, + HUpdate = HUpdate) +elif psph: + hydro = PSPH(dataBase = db, + cfl = cfl, + W = WT, + Q = q, + filter = filter, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate = densityUpdate, + HUpdate = HUpdate, + XSPH = xsph, + ASPH = asph) +if fsisph: + sumDensityNodeListSwitch =[nodes] + hydro = FSISPH(dataBase = db, + Q=q, + W = WT, + cfl = cfl, + surfaceForceCoefficient = fsiSurfaceCoefficient, + densityStabilizationCoefficient = fsiRhoStabilizeCoeff, + specificThermalEnergyDiffusionCoefficient = fsiEpsDiffuseCoeff, + xsphCoefficient = fsiXSPHCoeff, + interfaceMethod = fsiInterfaceMethod, + kernelAveragingMethod = fsiKernelMethod, + sumDensityNodeLists = sumDensityNodeListSwitch, + correctVelocityGradient = correctVelocityGradient, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + ASPH = asph, + epsTensile = epsilonTensile) +elif gsph: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + gradientType = riemannGradientType, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfm: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,gsphLinearCorrect) + hydro = MFM(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + gradientType = riemannGradientType, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + densityUpdate=densityUpdate, + XSPH = xsph, + ASPH = asph, + epsTensile = epsilonTensile, + nTensile = nTensile) else: - hydro = HydroConstructor(W = WT, - WPi = WTPi, - Q = q, - cfl = cfl, - useVelocityMagnitudeForDt = useVelocityMagnitudeForDt, - compatibleEnergyEvolution = compatibleEnergy, - gradhCorrection = gradhCorrection, - correctVelocityGradient = correctVelocityGradient, - evolveTotalEnergy = evolveTotalEnergy, - XSPH = XSPH, - densityUpdate = densityUpdate, - HUpdate = HUpdate, - epsTensile = epsilonTensile) + hydro = SPH(dataBase = db, + cfl = cfl, + W = WT, + Q = q, + compatibleEnergyEvolution = compatibleEnergy, + evolveTotalEnergy = evolveTotalEnergy, + gradhCorrection = gradhCorrection, + correctVelocityGradient = correctVelocityGradient, + XSPH = xsph, + ASPH = asph, + densityUpdate = densityUpdate, + HUpdate = HUpdate, + epsTensile = epsilonTensile, + nTensile = nTensile) + output("hydro") output("hydro.kernel()") -output("hydro.PiKernel()") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.densityUpdate") @@ -358,10 +430,10 @@ def __call__(self, r): #------------------------------------------------------------------------------- # Construct the MMRV physics object. #------------------------------------------------------------------------------- -if boolReduceViscosity: +if boolReduceViscosity and useArtificialViscosity: evolveReducingViscosityMultiplier = MorrisMonaghanReducingViscosity(q,nh,aMin,aMax) packages.append(evolveReducingViscosityMultiplier) -elif boolCullenViscosity: +elif boolCullenViscosity and useArtificialViscosity: evolveCullenViscosityMultiplier = CullenDehnenViscosity(q,WTPi,alphMax,alphMin,betaC,betaD,betaE,fKern,boolHopkinsCorrection) packages.append(evolveCullenViscosityMultiplier) @@ -407,8 +479,8 @@ def __call__(self, r): ylow.append(i) elif pos[i].y > y1: yhigh.append(i) -ybc1 = ConstantBoundary(nodes, ylow, yp1) -ybc2 = ConstantBoundary(nodes, yhigh, yp2) +ybc1 = ConstantBoundary(db,nodes, ylow, yp1) +ybc2 = ConstantBoundary(db,nodes, yhigh, yp2) bcSet = [ybc1, ybc2, xbc] # <-- ybc should be first! @@ -465,7 +537,7 @@ def __call__(self, r): restartBaseName = restartBaseName, restoreCycle = restoreCycle, redistributeStep = None, - vizMethod = vizMethod, + #vizMethod = vizMethod, vizBaseName = vizBaseName, vizDir = vizDir, vizStep = vizCycle, diff --git a/tests/functional/Hydro/Sedov/Sedov-cylindrical-2d.py b/tests/functional/Hydro/Sedov/Sedov-cylindrical-2d.py index e63fb8d82..29a6472f4 100644 --- a/tests/functional/Hydro/Sedov/Sedov-cylindrical-2d.py +++ b/tests/functional/Hydro/Sedov/Sedov-cylindrical-2d.py @@ -1,14 +1,19 @@ #------------------------------------------------------------------------------- # The cylindrical Sedov test case (2-D). #------------------------------------------------------------------------------- -import os, sys, shutil +import os, sys, shutil, mpi from Spheral2d import * from SpheralTestUtilities import * #from SpheralGnuPlotUtilities import * from GenerateNodeDistribution2d import * from CubicNodeGenerator import GenerateSquareNodeDistribution -import mpi +if mpi.procs > 1: + from VoronoiDistributeNodes import distributeNodes2d + #from PeanoHilbertDistributeNodes import distributeNodes2d +else: + from DistributeNodes import distributeNodes2d + title("2-D integrated hydro test -- planar Sedov problem") #------------------------------------------------------------------------------- @@ -22,8 +27,8 @@ nTheta = 50, rmin = 0.0, rmax = 1.0, - nPerh = 1.51, - order = 5, + nPerh = 1.00, + order = 3, rho0 = 1.0, eps0 = 0.0, @@ -46,6 +51,8 @@ psph = False, fsisph = False, gsph = False, + mfm = False, + mfv = False, # hydro options solid = False, @@ -63,7 +70,7 @@ volumeType = RKSumVolume, # gsph options - RiemannGradientType = RiemannGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) + RiemannGradientType = SPHGradient, # (RiemannGradient,SPHGradient,HydroAccelerationGradient,OnlyDvDxGradient,MixedMethodGradient) linearReconstruction = True, # Artifical Viscosity @@ -95,9 +102,11 @@ statsStep = 1, smoothIters = 0, useVelocityMagnitudeForDt = False, + dtverbose = False, # IO vizCycle = None, + vizDerivs = False, vizTime = 0.1, restoreCycle = -1, restartStep = 1000, @@ -117,7 +126,7 @@ assert not(boolReduceViscosity and boolCullenViscosity) assert thetaFactor in (0.5, 1.0, 2.0) -assert not(gsph and (boolReduceViscosity or boolCullenViscosity)) +assert not((gsph or mfm or mfv) and (boolReduceViscosity or boolCullenViscosity)) assert not(fsisph and not solid) theta = thetaFactor * pi @@ -154,6 +163,10 @@ hydroname = "FSISPH" elif gsph: hydroname = "GSPH" +elif mfm: + hydroname = "MFM" +elif mfv: + hydroname = "MFV" else: hydroname = "SPH" if asph: @@ -175,7 +188,6 @@ #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- -import os, sys if mpi.rank == 0: if clearDirectories and os.path.exists(dataDir): shutil.rmtree(dataDir) @@ -249,12 +261,6 @@ nNodePerh = nPerh, SPH = (not ASPH)) -if mpi.procs > 1: - from VoronoiDistributeNodes import distributeNodes2d - #from PeanoHilbertDistributeNodes import distributeNodes2d -else: - from DistributeNodes import distributeNodes2d - distributeNodes2d((nodes1, generator)) output("mpi.reduce(nodes1.numInternalNodes, mpi.MIN)") output("mpi.reduce(nodes1.numInternalNodes, mpi.MAX)") @@ -348,6 +354,42 @@ compatibleEnergyEvolution = compatibleEnergy, correctVelocityGradient= correctVelocityGradient, evolveTotalEnergy = evolveTotalEnergy, + gradientType = RiemannGradientType, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) +elif mfm: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction) + hydro = MFM(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + gradientType = RiemannGradientType, + XSPH = XSPH, + ASPH = asph, + densityUpdate=densityUpdate, + HUpdate = HUpdate) +elif mfv: + limiter = VanLeerLimiter() + waveSpeed = DavisWaveSpeed() + solver = HLLC(limiter,waveSpeed,linearReconstruction) + hydro = MFV(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + specificThermalEnergyDiffusionCoefficient = 0.00, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient= correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + gradientType = RiemannGradientType, + nodeMotionType=NodeMotionType.Fician, XSPH = XSPH, ASPH = asph, densityUpdate=densityUpdate, @@ -386,8 +428,10 @@ output("hydro.densityUpdate") output("hydro.HEvolution") -if not gsph: +if not (gsph or mfm or mfv): q = hydro.Q + q.Cq = 2 + q.Cl = 2 output("q") output("q.Cl") output("q.Cq") @@ -433,6 +477,7 @@ if dtMax: integrator.dtMax = dtMax integrator.dtGrowth = dtGrowth +integrator.verbose = dtverbose integrator.allowDtCheck = True output("integrator") output("integrator.havePhysicsPackage(hydro)") @@ -457,6 +502,7 @@ restoreCycle = restoreCycle, vizMethod = vizMethod, vizBaseName = "Sedov-cylindrical-2d-%ix%i" % (nRadial, nTheta), + vizDerivs=vizDerivs, vizDir = vizDir, vizStep = vizCycle, vizTime = vizTime, @@ -579,6 +625,7 @@ # Plot the final state. #------------------------------------------------------------------------------- if graphics: + from SpheralMatplotlib import * rPlot = plotNodePositions2d(db, colorNodeLists=0, colorDomains=1) rhoPlot, velPlot, epsPlot, PPlot, HPlot = plotRadialState(db) plotAnswer(answer, control.time(), @@ -591,30 +638,30 @@ (HPlot, "Sedov-cylindrical-h.png")] # Plot the specific entropy. - AsimData = Gnuplot.Data(xprof, A, - with_ = "points", - title = "Simulation", - inline = True) - AansData = Gnuplot.Data(xprof, Aans, - with_ = "lines", - title = "Solution", - inline = True) + # AsimData = Gnuplot.Data(xprof, A, + # with_ = "points", + # title = "Simulation", + # inline = True) + # AansData = Gnuplot.Data(xprof, Aans, + # with_ = "lines", + # title = "Solution", + # inline = True) - Aplot = generateNewGnuPlot() - Aplot.plot(AsimData) - Aplot.replot(AansData) - Aplot.title("Specific entropy") - Aplot.refresh() - plots.append((Aplot, "Sedov-cylindrical-entropy.png")) - - if boolCullenViscosity: - cullAlphaPlot = plotFieldList(q.ClMultiplier(), - xFunction = "%s.magnitude()", - plotStyle = "points", - winTitle = "Cullen alpha") - plots += [(cullAlphaPlot, "Sedov-planar-Cullen-alpha.png")] + # Aplot = generateNewGnuPlot() + # Aplot.plot(AsimData) + # Aplot.replot(AansData) + # Aplot.title("Specific entropy") + # Aplot.refresh() + # plots.append((Aplot, "Sedov-cylindrical-entropy.png")) + + # if boolCullenViscosity: + # cullAlphaPlot = plotFieldList(q.ClMultiplier(), + # xFunction = "%s.magnitude()", + # plotStyle = "points", + # winTitle = "Cullen alpha") + # plots += [(cullAlphaPlot, "Sedov-planar-Cullen-alpha.png")] # Make hardcopies of the plots. for p, filename in plots: - p.hardcopy(os.path.join(dataDir, filename), terminal="png") + p.figure.savefig(os.path.join(dataDir, filename)) diff --git a/tests/functional/Hydro/YeeVortex/YeeVortex.py b/tests/functional/Hydro/YeeVortex/YeeVortex.py index 404443b6c..0f7b51829 100644 --- a/tests/functional/Hydro/YeeVortex/YeeVortex.py +++ b/tests/functional/Hydro/YeeVortex/YeeVortex.py @@ -5,7 +5,7 @@ #------------------------------------------------------------------------------- # The Yee-Vortex Test #------------------------------------------------------------------------------- -import shutil +import shutil, os, sys, mpi from math import * from Spheral2d import * from SpheralTestUtilities import * @@ -13,9 +13,7 @@ from findLastRestart import * from GenerateNodeDistribution2d import * from CubicNodeGenerator import GenerateSquareNodeDistribution -from CentroidalVoronoiRelaxation import * - -import mpi +from CentroidalVoronoiRelaxation import * import DistributeNodes class YeeDensity: @@ -53,16 +51,17 @@ def __call__(self, r): #Center and radius of Vortex xc=0.0, yc=0.0, - rmax = 5.0, + rmax = 6.0, # How far should we measure the error norms? rmaxnorm = 5.0, # The number of radial points on the outside to force with constant BC - nbcrind = 10, + nbcrind = 6, #Vortex strength beta = 5.0, + #Tempurature at inf temp_inf = 1.0, @@ -70,19 +69,53 @@ def __call__(self, r): nRadial = 64, seed = "constantDTheta", - nPerh = 1.51, + # kernel options + KernelConstructor = WendlandC2Kernel, + nPerh = 3.01, + order = 7, + hmin = 1e-5, + hmax = 0.5, + hminratio = 0.1, + # hydros svph = False, crksph = False, fsisph = False, psph = False, + gsph = False, + mfm = False, + + # general hydro options asph = False, solid = False, + XSPH = False, + epsilonTensile = 0.0, + nTensile = 8, + densityUpdate = RigorousSumDensity, # VolumeScaledDensity, + compatibleEnergy = True, + evolveTotalEnergy = False, + correctVelocityGradient = True, + + # default SPH options + gradhCorrection = True, + + # CRKSPH options filter = 0.0, # For CRKSPH - KernelConstructor = NBSplineKernel, - order = 5, - Qconstructor = MonaghanGingoldViscosity, - #Qconstructor = TensorMonaghanGingoldViscosity, + + # PSPH options + HopkinsConductivity = False, + XPSH=False, + + # MFM/GSPH options + WaveSpeedConstructor = DavisWaveSpeed, # Einfeldt, Acoustic + LimiterConstructor = VanLeerLimiter, # VanLeer, Opsre, MinMod, VanAlba, Superbee + riemannLinearReconstruction = True, + riemannGradientType = SPHSameTimeGradient, # HydroAccelerationGradient, SPHGradient, RiemannGradient, MixedMethodGradient, SPHSameTimeGradient + + # artificial viscosity + Qconstructor = LimitedMonaghanGingoldViscosity, # TensorMonaghanGingoldViscosity, + Cl = 1.0, + Cq = 1.0, boolReduceViscosity = False, nhQ = 5.0, nhL = 10.0, @@ -99,25 +132,18 @@ def __call__(self, r): linearConsistent = False, fcentroidal = 0.0, fcellPressure = 0.0, - Cl = 1.0, - Cq = 0.75, linearInExpansion = False, Qlimiter = False, balsaraCorrection = False, epsilon2 = 1e-2, - hmin = 1e-5, - hmax = 0.5, - hminratio = 0.1, - cfl = 0.5, - XSPH = False, - epsilonTensile = 0.0, - nTensile = 8, - + + # integrator IntegratorConstructor = CheapSynchronousRK2Integrator, + cfl = 0.25, goalTime = 8.0, steps = None, - vizCycle = 20, - vizTime = 0.1, + vizCycle = None, + vizTime = 2.0, dt = 0.0001, dtMin = 1.0e-5, dtMax = 1.0, @@ -128,23 +154,17 @@ def __call__(self, r): domainIndependent = False, rigorousBoundaries = False, dtverbose = False, - - densityUpdate = RigorousSumDensity, # VolumeScaledDensity, - compatibleEnergy = True, - gradhCorrection = True, - HopkinsConductivity = False, # For PSPH - correctVelocityGradient = True, - evolveTotalEnergy = False, - XPSH=False, - + + # output useVoronoiOutput = False, clearDirectories = False, restoreCycle = -1, restartStep = 200, dataDir = "dumps-yeevortex-xy", graphics = True, - smooth = None, - outputFile = "None", + smooth = False, + outputFileBase = ".out", + convergenceFileBase = "xstaglattice_converge.txt", ) assert not(boolReduceViscosity and boolCullenViscosity) @@ -159,7 +179,12 @@ def __call__(self, r): elif psph: hydroname = "PSPH" elif fsisph: + Qconstructor = LimitedMonaghanGingoldViscosity hydroname = "FSISPH" +elif mfm: + hydroname = "MFM" +elif gsph: + hydroname = "GSPH" else: hydroname = "SPH" if asph: @@ -167,6 +192,13 @@ def __call__(self, r): if solid: hydroname = "solid"+hydroname + +if mfm or gsph: + convergenceFile = hydroname + "_" + str(densityUpdate) + "_" + str(riemannGradientType) + "_" + convergenceFileBase + outputFile = hydroname + "_" + str(densityUpdate) + "_" + str(riemannGradientType) + "_" + str(nRadial) + ".out" +else: + convergenceFile = hydroname+"_"+str(densityUpdate) + "_" + convergenceFileBase + outputFile = hydroname + "_" + str(densityUpdate) + "_" + str(nRadial) + ".out" #------------------------------------------------------------------------------- # Build our directory paths. #------------------------------------------------------------------------------- @@ -176,16 +208,18 @@ def __call__(self, r): SumVoronoiCellDensity : "SumVoronoiCellDensity"} baseDir = os.path.join(dataDir, hydroname, - Qconstructor.__name__, - KernelConstructor.__name__, - "Cl=%g_Cq=%g" % (Cl, Cq), - densityUpdateLabel[densityUpdate], - "compatibleEnergy=%s" % compatibleEnergy, - "Cullen=%s" % boolCullenViscosity, - "nPerh=%3.1f" % nPerh, - "fcentroidal=%f" % max(fcentroidal, filter), - "fcellPressure=%f" % fcellPressure, - "seed=%s" % seed, + #str(riemannGradientType), + #LimiterConstructor.__name__, + #Qconstructor.__name__, + #KernelConstructor.__name__, + #"Cl=%g_Cq=%g" % (Cl, Cq), + #densityUpdateLabel[densityUpdate], + #"compatibleEnergy=%s" % compatibleEnergy, + #"Cullen=%s" % boolCullenViscosity, + #"nPerh=%3.1f" % nPerh, + #"fcentroidal=%f" % max(fcentroidal, filter), + #"fcellPressure=%f" % fcellPressure, + #"seed=%s" % seed, str(nRadial)) restartDir = os.path.join(baseDir, "restarts") restartBaseName = os.path.join(restartDir, "yeevortex-xy-%i" % nRadial) @@ -199,7 +233,6 @@ def __call__(self, r): #------------------------------------------------------------------------------- # Check if the necessary output directories exist. If not, create them. #------------------------------------------------------------------------------- -import os, sys if mpi.rank == 0: if clearDirectories and os.path.exists(baseDir): shutil.rmtree(baseDir) @@ -216,7 +249,7 @@ def __call__(self, r): K = 1.0 eos = GammaLawGasMKS(gamma, mu) #eos = PolytropicEquationOfStateMKS(K,gamma,mu) - +YeeDensityFunc = YeeDensity(xc,yc,gamma,beta,temp_inf) #------------------------------------------------------------------------------- # Interpolation kernels. #------------------------------------------------------------------------------- @@ -250,10 +283,11 @@ def __call__(self, r): # Set the node properties. #------------------------------------------------------------------------------- rmaxbound = rmax + rmax/nRadial*nbcrind +print rmaxbound nr1 = nRadial + nbcrind -if seed == "lattice": +if seed == "lattice" or "xstaggeredLattice": generator = GenerateNodeDistribution2d(2*nr1, 2*nr1, - rho = YeeDensity(xc,yc,gamma,beta,temp_inf), + rho = YeeDensityFunc, distributionType = seed, xmin = (-rmaxbound, -rmaxbound), xmax = (rmaxbound, rmaxbound), @@ -263,8 +297,9 @@ def __call__(self, r): nNodePerh = nPerh, SPH = SPH) else: + generator = GenerateNodeDistribution2d(nr1, nr1, - rho = YeeDensity(xc,yc,gamma,beta,temp_inf), + rho = YeeDensityFunc, distributionType = seed, xmin = (-rmaxbound, -rmaxbound), xmax = (rmaxbound, rmaxbound), @@ -358,10 +393,15 @@ def __call__(self, r): densityUpdate = densityUpdate, HUpdate = HUpdate) elif fsisph: + if densityUpdate==RigorousSumDensity: + sumDensityNodeLists = [nodes] + else: + sumDensityNodeLists = [] hydro = FSISPH(dataBase = db, Q=q, W = WT, - cfl = cfl, + cfl = cfl, + sumDensityNodeLists = sumDensityNodeLists, densityStabilizationCoefficient = 0.1, specificThermalEnergyDiffusionCoefficient = 0.1, linearCorrectGradients = correctVelocityGradient, @@ -370,6 +410,40 @@ def __call__(self, r): ASPH = asph, epsTensile = epsilonTensile, nTensile = nTensile) +elif gsph: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,riemannLinearReconstruction) + hydro = GSPH(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + gradientType = riemannGradientType, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) +elif mfm: + limiter = LimiterConstructor() + waveSpeed = WaveSpeedConstructor() + solver = HLLC(limiter,waveSpeed,riemannLinearReconstruction) + hydro = MFM(dataBase = db, + riemannSolver = solver, + W = WT, + cfl=cfl, + gradientType = riemannGradientType, + compatibleEnergyEvolution = compatibleEnergy, + correctVelocityGradient=correctVelocityGradient, + evolveTotalEnergy = evolveTotalEnergy, + XSPH = XSPH, + densityUpdate=densityUpdate, + HUpdate = IdealH, + epsTensile = epsilonTensile, + nTensile = nTensile) elif psph: hydro = PSPH(dataBase=db, W=WT, @@ -487,19 +561,23 @@ def __call__(self, r): #------------------------------------------------------------------------------- # Make the problem controller. #------------------------------------------------------------------------------- -# if useVoronoiOutput: -# import SpheralVoronoiSiloDump -# vizMethod = SpheralVoronoiSiloDump.dumpPhysicsState -# else: -# vizMethod = None # default +if useVoronoiOutput: + import SpheralVoronoiSiloDump + vizMethod = SpheralVoronoiSiloDump.dumpPhysicsState +else: + vizMethod = None # default +from SpheralPointmeshSiloDump import dumpPhysicsState + control = SpheralController(integrator, WT, initializeDerivatives = True, statsStep = statsStep, restartStep = restartStep, restartBaseName = restartBaseName, restoreCycle = restoreCycle, - #vizMethod = vizMethod, + #vizMethod = dumpPhysicsState, + #vizGhosts=True, vizBaseName = vizBaseName, + vizDerivs=True, vizDir = vizDir, vizStep = vizCycle, vizTime = vizTime, @@ -540,6 +618,7 @@ def __call__(self, r): if mpi.rank == 0: import numpy as np from Pnorm import Pnorm + rprof = np.array([sqrt(xi*xi + yi*yi) for xi, yi in zip(xprof, yprof)]) multiSort(rprof, mo, xprof, yprof, rhoprof, Pprof, vprof, epsprof, hprof,velx,vely) epsans = [] @@ -560,9 +639,13 @@ def __call__(self, r): rhoans.append(rhoi) velans.append(Vector(velxans,velyans).magnitude()) Pans.append(temp*rhoi) + L1rho2 =sum([abs(rhoprof[i]-rhoans[i]) for i in range(len(rhoans))])/len(rhoans) + Linfrho2 =max([abs(rhoprof[i]-rhoans[i]) for i in range(len(rhoans))])/len(rhoans) L1rho = Pnorm(rhoprof, rprof, rhoans).pnorm(1, rmin=0.0, rmax=rmaxnorm) + print (L1rho2,L1rho) L2rho = Pnorm(rhoprof, rprof, rhoans).pnorm(2, rmin=0.0, rmax=rmaxnorm) Linfrho = Pnorm(rhoprof, rprof, rhoans).pnorm("inf", rmin=0.0, rmax=rmaxnorm) + print (Linfrho2,Linfrho) L1eps = Pnorm(epsprof, rprof, epsans).pnorm(1, rmin=0.0, rmax=rmaxnorm) L2eps = Pnorm(epsprof, rprof, epsans).pnorm(2, rmin=0.0, rmax=rmaxnorm) Linfeps = Pnorm(epsprof, rprof, epsans).pnorm("inf", rmin=0.0, rmax=rmaxnorm) @@ -572,8 +655,11 @@ def __call__(self, r): L1P = Pnorm(Pprof, rprof, Pans).pnorm(1, rmin=0.0, rmax=rmaxnorm) L2P = Pnorm(Pprof, rprof, Pans).pnorm(2, rmin=0.0, rmax=rmaxnorm) LinfP = Pnorm(Pprof, rprof, velans).pnorm("inf", rmin=0.0, rmax=rmaxnorm) - with open("converge-CRK-%s-cullen-%s-PSPH-%s.txt" % (CRKSPH,boolCullenViscosity,PSPH), "a") as myfile: - myfile.write(("#" + 14*"%16s\t " + "%16s\n") % ("nRadial", "L1rho", "L1eps", "L1vel", "L2rho", "L2eps", "L2vel", "Linfrho", "Linfeps", "Linfvel", "L1P", "L2P", "LinfP", "cycles", "runtime")) + + isNewFile = not os.path.exists(convergenceFile) + with open(convergenceFile, "a") as myfile: + if isNewFile: + myfile.write(("#" + 14*"%16s\t " + "%16s\n") % ("nRadial", "L1rho", "L1eps", "L1vel", "L2rho", "L2eps", "L2vel", "Linfrho", "Linfeps", "Linfvel", "L1P", "L2P", "LinfP", "cycles", "runtime")) myfile.write((14*"%16s\t " + "%16s\n") % (nRadial, L1rho, L1eps, L1vel, L2rho, L2eps, L2vel, Linfrho, Linfeps, Linfvel, L1P, L2P, LinfP, control.totalSteps, control.stepTimer.elapsedTime)) f = open(outputFile, "w") f.write(("# " + 19*"%15s " + "\n") % ("r", "x", "y", "rho", "P", "v", "eps", "h", "mortonOrder", "rhoans", "epsans", "velans",