Skip to content

Commit

Permalink
Merge pull request #273 from LLNL/solidGSPH
Browse files Browse the repository at this point in the history 
GSPH update
  • Loading branch information
@jmpearl
jmpearl authored Jun 6, 2024
2 parents 426d5af + 22b202d commit 97ddcbe
Show file tree
Hide file tree
Showing 71 changed files with 4,377 additions and 1,404 deletions.
1 change: 1 addition & 0 deletions RELEASE_NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
18 changes: 16 additions & 2 deletions src/GSPH/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
95 changes: 57 additions & 38 deletions src/GSPH/GSPHEvaluateDerivatives.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Dimension, Vector>::prefix() + HydroFieldNames::position, Vector::zero);
auto DrhoDt = derivatives.fields(IncrementState<Dimension, Scalar>::prefix() + HydroFieldNames::massDensity, 0.0);
Expand All @@ -83,6 +84,7 @@ evaluateDerivatives(const typename Dimension::Scalar time,
auto newRiemannDpDx = derivatives.fields(ReplaceState<Dimension, Scalar>::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero);
auto newRiemannDvDx = derivatives.fields(ReplaceState<Dimension, Scalar>::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero);

CHECK(DrhoDx.size() == numNodeLists);
CHECK(M.size() == numNodeLists);
CHECK(normalization.size() == numNodeLists);
CHECK(DxDt.size() == numNodeLists);
Expand Down Expand Up @@ -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);
Expand All @@ -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);


Expand All @@ -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);
Expand All @@ -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();
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -434,34 +434,39 @@ computeMCorrection(const typename Dimension::Scalar /*time*/,
const DataBase<Dimension>& dataBase,
const State<Dimension>& state,
StateDerivatives<Dimension>& 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);
CHECK(position.size() == numNodeLists);
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<Dimension, Scalar>::prefix() + GSPHFieldNames::RiemannPressureGradient,Vector::zero);
auto newRiemannDvDx = derivatives.fields(ReplaceState<Dimension, Scalar>::prefix() + GSPHFieldNames::RiemannVelocityGradient,Tensor::zero);

CHECK(M.size() == numNodeLists);
CHECK(DrhoDx.size() == numNodeLists);
CHECK(newRiemannDpDx.size() == numNodeLists);
CHECK(newRiemannDvDx.size() == numNodeLists);

Expand All @@ -476,6 +481,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/,

typename SpheralThreads<Dimension>::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);

Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -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;

Expand All @@ -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);
Expand All @@ -546,6 +561,7 @@ computeMCorrection(const typename Dimension::Scalar /*time*/,

newRiemannDvDxi -= (vi-vj).dyad(gradPsii);
newRiemannDvDxj -= (vi-vj).dyad(gradPsij);

}
} // loop over pairs

Expand All @@ -554,14 +570,15 @@ 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();
#pragma omp parallel for
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());

Expand All @@ -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);
}
Expand All @@ -597,6 +617,5 @@ computeMCorrection(const typename Dimension::Scalar /*time*/,
boundaryItr != this->boundaryEnd();
++boundaryItr) (*boundaryItr)->finalizeGhostBoundary();

}

} // MC correction method
} // spheral namespace
6 changes: 5 additions & 1 deletion src/GSPH/GSPHFieldNames.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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";
const std::string Spheral::GSPHFieldNames::RiemannDeviatoricStressTensorGradient = "Riemann solvers deviatoric stress tensor gradient";
const std::string Spheral::GSPHFieldNames::pairMassFlux = "pairwise mass flux";
4 changes: 4 additions & 0 deletions src/GSPH/GSPHFieldNames.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};

}
Expand Down
Loading

0 comments on commit 97ddcbe

Please sign in to comment.