Skip to content

Commit 7b5c0b8

Browse files
jmpearljmpearl
authored
Merge pull request #259 from LLNL/bugfix/FSISPH-ALE3D
Bugfix/fsisph ale3 d
2 parents ad4d3cd + a22b124 commit 7b5c0b8

5 files changed

+30
-3
lines changed

‎src/FSISPH/SolidFSISPHEvaluateDerivatives.cc

+13-3
Original file line numberDiff line numberDiff line change
@@ -146,6 +146,7 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
146146
auto DHDt = derivatives.fields(IncrementState<Dimension, SymTensor>::prefix() + HydroFieldNames::H, SymTensor::zero);
147147
auto Hideal = derivatives.fields(ReplaceBoundedState<Dimension, SymTensor>::prefix() + HydroFieldNames::H, SymTensor::zero);
148148
auto maxViscousPressure = derivatives.fields(HydroFieldNames::maxViscousPressure, 0.0);
149+
auto effViscousPressure = derivatives.fields(HydroFieldNames::effectiveViscousPressure, 0.0);
149150
auto XSPHWeightSum = derivatives.fields(HydroFieldNames::XSPHWeightSum, 0.0);
150151
auto XSPHDeltaV = derivatives.fields(HydroFieldNames::XSPHDeltaV, Vector::zero);
151152
auto weightedNeighborSum = derivatives.fields(HydroFieldNames::weightedNeighborSum, 0.0);
@@ -177,6 +178,7 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
177178
CHECK(DHDt.size() == numNodeLists);
178179
CHECK(Hideal.size() == numNodeLists);
179180
CHECK(maxViscousPressure.size() == numNodeLists);
181+
CHECK(effViscousPressure.size() == numNodeLists);
180182
CHECK(XSPHWeightSum.size() == numNodeLists);
181183
CHECK(XSPHDeltaV.size() == numNodeLists);
182184
CHECK(weightedNeighborSum.size() == numNodeLists);
@@ -222,6 +224,7 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
222224
auto weightedNeighborSum_thread = weightedNeighborSum.threadCopy(threadStack);
223225
auto massSecondMoment_thread = massSecondMoment.threadCopy(threadStack);
224226
auto maxViscousPressure_thread = maxViscousPressure.threadCopy(threadStack, ThreadReduction::MAX);
227+
auto effViscousPressure_thread = effViscousPressure.threadCopy(threadStack);
225228

226229
#pragma omp for
227230
for (auto kk = 0u; kk < numPairs; ++kk) {
@@ -274,6 +277,7 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
274277
auto& weightedNeighborSumi = weightedNeighborSum_thread(nodeListi, i);
275278
auto& massSecondMomenti = massSecondMoment_thread(nodeListi, i);
276279
auto& maxViscousPressurei = maxViscousPressure_thread(nodeListi, i);
280+
auto& effViscousPressurei = effViscousPressure_thread(nodeListi, i);
277281
auto& newInterfaceFlagsi = newInterfaceFlags_thread(nodeListi,i);
278282
auto& newInterfaceAreaVectorsi = newInterfaceAreaVectors_thread(nodeListi,i);
279283
auto& newInterfaceNormalsi = newInterfaceNormals_thread(nodeListi,i);
@@ -325,6 +329,7 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
325329
auto& weightedNeighborSumj = weightedNeighborSum_thread(nodeListj, j);
326330
auto& massSecondMomentj = massSecondMoment_thread(nodeListj, j);
327331
auto& maxViscousPressurej = maxViscousPressure_thread(nodeListj, j);
332+
auto& effViscousPressurej = effViscousPressure_thread(nodeListj, j);
328333
auto& newInterfaceFlagsj = newInterfaceFlags_thread(nodeListj,j);
329334
auto& newInterfaceAreaVectorsj = newInterfaceAreaVectors_thread(nodeListj,j);
330335
auto& newInterfaceNormalsj = newInterfaceNormals_thread(nodeListj,j);
@@ -492,9 +497,14 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
492497
}
493498

494499
// save our max pressure from the AV for each node
495-
maxViscousPressurei = max(maxViscousPressurei, rhoi*rhoj * QPiij.diagonalElements().maxAbsElement());
496-
maxViscousPressurej = max(maxViscousPressurej, rhoi*rhoj * QPiji.diagonalElements().maxAbsElement());
497-
500+
{
501+
const auto Qi = rhoi*rhoj * QPiij.diagonalElements().maxAbsElement();
502+
const auto Qj = rhoi*rhoj * QPiji.diagonalElements().maxAbsElement();
503+
maxViscousPressurei = max(maxViscousPressurei, Qi);
504+
maxViscousPressurej = max(maxViscousPressurej, Qj);
505+
effViscousPressurei += volj * Qi * Wi;
506+
effViscousPressurej += voli * Qj * Wj;
507+
}
498508
// stress tensor
499509
//{
500510
// apply yield pairwise

‎src/FSISPH/SolidFSISPHHydroBase.cc

+6
Original file line numberDiff line numberDiff line change
@@ -200,6 +200,7 @@ SolidFSISPHHydroBase(const SmoothingScaleBase<Dimension>& smoothingScaleMethod,
200200
mM(FieldStorageType::CopyFields),
201201
mLocalM(FieldStorageType::CopyFields),
202202
mMaxViscousPressure(FieldStorageType::CopyFields),
203+
mEffViscousPressure(FieldStorageType::CopyFields),
203204
mNormalization(FieldStorageType::CopyFields),
204205
mWeightedNeighborSum(FieldStorageType::CopyFields),
205206
mMassSecondMoment(FieldStorageType::CopyFields),
@@ -253,6 +254,7 @@ SolidFSISPHHydroBase(const SmoothingScaleBase<Dimension>& smoothingScaleMethod,
253254
mM = dataBase.newFluidFieldList(Tensor::zero, HydroFieldNames::M_SPHCorrection);
254255
mLocalM = dataBase.newFluidFieldList(Tensor::zero, "local " + HydroFieldNames::M_SPHCorrection);
255256
mMaxViscousPressure = dataBase.newFluidFieldList(0.0, HydroFieldNames::maxViscousPressure);
257+
mEffViscousPressure = dataBase.newFluidFieldList(0.0, HydroFieldNames::effectiveViscousPressure);
256258
mNormalization = dataBase.newFluidFieldList(0.0, HydroFieldNames::normalization);
257259
mWeightedNeighborSum = dataBase.newFluidFieldList(0.0, HydroFieldNames::weightedNeighborSum);
258260
mMassSecondMoment = dataBase.newFluidFieldList(SymTensor::zero, HydroFieldNames::massSecondMoment);
@@ -455,6 +457,7 @@ registerDerivatives(DataBase<Dimension>& dataBase,
455457
dataBase.resizeFluidFieldList(mM, Tensor::zero, HydroFieldNames::M_SPHCorrection, false);
456458
dataBase.resizeFluidFieldList(mLocalM, Tensor::zero, "local " + HydroFieldNames::M_SPHCorrection, false);
457459
dataBase.resizeFluidFieldList(mMaxViscousPressure, 0.0, HydroFieldNames::maxViscousPressure, false);
460+
dataBase.resizeFluidFieldList(mEffViscousPressure, 0.0, HydroFieldNames::effectiveViscousPressure, false);
458461
dataBase.resizeFluidFieldList(mNormalization, 0.0, HydroFieldNames::normalization, false);
459462
dataBase.resizeFluidFieldList(mWeightedNeighborSum, 0.0, HydroFieldNames::weightedNeighborSum, false);
460463
dataBase.resizeFluidFieldList(mMassSecondMoment, SymTensor::zero, HydroFieldNames::massSecondMoment, false);
@@ -492,6 +495,7 @@ registerDerivatives(DataBase<Dimension>& dataBase,
492495
derivs.enroll(mM);
493496
derivs.enroll(mLocalM);
494497
derivs.enroll(mMaxViscousPressure);
498+
derivs.enroll(mEffViscousPressure);
495499
derivs.enroll(mNormalization);
496500
derivs.enroll(mWeightedNeighborSum);
497501
derivs.enroll(mMassSecondMoment);
@@ -768,6 +772,7 @@ dumpState(FileIO& file, const string& pathName) const {
768772
file.write(mM, pathName + "/M");
769773
file.write(mLocalM, pathName + "/localM");
770774
file.write(mMaxViscousPressure, pathName + "/maxViscousPressure");
775+
file.write(mEffViscousPressure, pathName + "/effectiveViscousPressure");
771776
file.write(mNormalization, pathName + "/normalization");
772777
file.write(mWeightedNeighborSum, pathName + "/weightedNeighborSum");
773778
file.write(mMassSecondMoment, pathName + "/massSecondMoment");
@@ -819,6 +824,7 @@ restoreState(const FileIO& file, const string& pathName) {
819824
file.read(mM, pathName + "/M");
820825
file.read(mLocalM, pathName + "/localM");
821826
file.read(mMaxViscousPressure, pathName + "/maxViscousPressure");
827+
file.read(mEffViscousPressure, pathName + "/effectiveViscousPressure");
822828
file.read(mNormalization, pathName + "/normalization");
823829
file.read(mWeightedNeighborSum, pathName + "/weightedNeighborSum");
824830
file.read(mMassSecondMoment, pathName + "/massSecondMoment");

‎src/FSISPH/SolidFSISPHHydroBase.hh

+2
Original file line numberDiff line numberDiff line change
@@ -244,6 +244,7 @@ public:
244244
const FieldList<Dimension, Tensor>& M() const;
245245
const FieldList<Dimension, Tensor>& localM() const;
246246
const FieldList<Dimension, Scalar>& maxViscousPressure() const;
247+
const FieldList<Dimension, Scalar>& effectiveViscousPressure() const;
247248
const FieldList<Dimension, Scalar>& normalization() const;
248249
const FieldList<Dimension, Scalar>& weightedNeighborSum() const;
249250
const FieldList<Dimension, SymTensor>& massSecondMoment() const;
@@ -327,6 +328,7 @@ private:
327328
FieldList<Dimension, Tensor> mM;
328329
FieldList<Dimension, Tensor> mLocalM;
329330
FieldList<Dimension, Scalar> mMaxViscousPressure;
331+
FieldList<Dimension, Scalar> mEffViscousPressure;
330332
FieldList<Dimension, Scalar> mNormalization;
331333
FieldList<Dimension, Scalar> mWeightedNeighborSum;
332334
FieldList<Dimension, SymTensor> mMassSecondMoment;

‎src/FSISPH/SolidFSISPHHydroBaseInline.hh

+8
Original file line numberDiff line numberDiff line change
@@ -601,6 +601,14 @@ maxViscousPressure() const {
601601
return mMaxViscousPressure;
602602
}
603603

604+
template<typename Dimension>
605+
inline
606+
const FieldList<Dimension, typename Dimension::Scalar>&
607+
SolidFSISPHHydroBase<Dimension>::
608+
effectiveViscousPressure() const {
609+
return mEffViscousPressure;
610+
}
611+
604612
template<typename Dimension>
605613
inline
606614
const FieldList<Dimension, typename Dimension::Scalar>&

‎src/PYB11/FSISPH/SolidFSISPHHydroBase.py

+1
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,7 @@ def registerDerivatives(dataBase = "DataBase<%(Dimension)s>&",
148148
M = PYB11property("const FieldList<%(Dimension)s, Tensor>&", "M", returnpolicy="reference_internal")
149149
localM = PYB11property("const FieldList<%(Dimension)s, Tensor>&", "localM", returnpolicy="reference_internal")
150150
maxViscousPressure = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "maxViscousPressure", returnpolicy="reference_internal")
151+
effectiveViscousPressure = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "effectiveViscousPressure", returnpolicy="reference_internal")
151152
normalization = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "normalization", returnpolicy="reference_internal")
152153
weightedNeighborSum = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "weightedNeighborSum", returnpolicy="reference_internal")
153154
massSecondMoment = PYB11property("const FieldList<%(Dimension)s, SymTensor>&","massSecondMoment", returnpolicy="reference_internal")

0 commit comments

Comments
 (0)