Skip to content

Commit dcd1080

Browse files
jmikeowenjmikeowen
authored
Merge pull request #264 from LLNL/bugfix/domainParallelism
Fixes some initialization problems with the damage models volume definitions.
2 parents d3c5737 + ccc2d61 commit dcd1080

16 files changed

+137
-90
lines changed

‎RELEASE_NOTES.md

+1
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ Notable changes include:
1515

1616
* Bug Fixes / improvements:
1717
* Fixed bug with ConstantBoundary in the presence of porosity with the new porosity models introduced in v2024.01.00.
18+
* Initial volumes for damage models were incorrectly not taking into account pore space when computing failure statistics for seeding flaws. Fixed.
1819

1920
Version v2024.01.00 -- Release date 2024-01-19
2021
==============================================

‎src/Damage/GradyKippTensorDamage.py

+38-40
Original file line numberDiff line numberDiff line change
@@ -153,34 +153,33 @@ def __init__(self, *args, **kwargs):
153153
if weibull_kwargs["mask"] is None:
154154
weibull_kwargs["mask"] = IntField("mask", damage_kwargs["nodeList"], 1)
155155

156-
# Build the flaw distribution.
157-
damage_kwargs["flaws"] = weibullFlawDistributionBenzAsphaug(**weibull_kwargs)
156+
# Preserve the input for constructing the flaws
157+
self.weibull_kwargs = weibull_kwargs
158158

159159
# Invoke the parent constructor.
160160
TensorDamageModel.__init__(self, **damage_kwargs)
161161
return
162162

163+
########################################################################
164+
def initializeProblemStartupDependencies(self,
165+
dataBase,
166+
state,
167+
derivs):
168+
# Set the flaws
169+
self.weibull_kwargs["state"] = state
170+
self.flaws = weibullFlawDistributionBenzAsphaug(**self.weibull_kwargs)
171+
172+
TensorDamageModel.initializeProblemStartupDependencies(self,
173+
dataBase,
174+
state,
175+
derivs)
176+
177+
return
178+
179+
########################################################################
163180
def label(self):
164181
return "GradyKippTensorDamageBenzAsphaug"
165182

166-
# def dumpState(self,
167-
# file,
168-
# pathName):
169-
# TensorDamageModel.dumpState(self, file, pathName)
170-
# #file.writeObject(self.kWeibull, pathName + "/kWeibull")
171-
# #file.writeObject(self.mWeibull, pathName + "/mWeibull")
172-
# #file.writeObject(self.seed, pathName + "/seed")
173-
# return
174-
175-
# def restoreState(self,
176-
# file,
177-
# pathName):
178-
# TensorDamageModel.restoreState(self, file, pathName)
179-
# #self.kWeibull = file.readObject(pathName + "/kWeibull")
180-
# #self.mWeibull = file.readObject(pathName + "/mWeibull")
181-
# #self.seed = file.readObject(pathName + "/seed")
182-
# return
183-
184183
return GradyKippTensorDamageBenzAsphaug
185184

186185
#-------------------------------------------------------------------------------
@@ -324,35 +323,34 @@ def __init__(self, *args, **kwargs):
324323
if weibull_kwargs["mask"] is None:
325324
weibull_kwargs["mask"] = IntField("mask", damage_kwargs["nodeList"], 1)
326325

327-
# Build the flaw distribution.
328-
damage_kwargs["flaws"] = weibullFlawDistributionOwen(**weibull_kwargs)
326+
# Preserve the input for constructing the flaws, and build dummy flaws for now
327+
self.weibull_kwargs = weibull_kwargs
329328

330329
# Invoke the parent constructor.
331330
TensorDamageModel.__init__(self, **damage_kwargs)
332331

333332
return
334333

334+
########################################################################
335+
def initializeProblemStartupDependencies(self,
336+
dataBase,
337+
state,
338+
derivs):
339+
# Set the flaws
340+
self.weibull_kwargs["state"] = state
341+
self.flaws = weibullFlawDistributionOwen(**self.weibull_kwargs)
342+
343+
TensorDamageModel.initializeProblemStartupDependencies(self,
344+
dataBase,
345+
state,
346+
derivs)
347+
348+
return
349+
350+
########################################################################
335351
def label(self):
336352
return "GradyKippTensorDamageOwen"
337353

338-
# def dumpState(self,
339-
# file,
340-
# pathName):
341-
# TensorDamageModel.dumpState(self, file, pathName)
342-
# # file.writeObject(self.kWeibull, pathName + "/kWeibull")
343-
# # file.writeObject(self.mWeibull, pathName + "/mWeibull")
344-
# # file.writeObject(self.seed, pathName + "/seed")
345-
# return
346-
347-
# def restoreState(self,
348-
# file,
349-
# pathName):
350-
# TensorDamageModel.restoreState(self, file, pathName)
351-
# # self.kWeibull = file.readObject(pathName + "/kWeibull")
352-
# # self.mWeibull = file.readObject(pathName + "/mWeibull")
353-
# # self.seed = file.readObject(pathName + "/seed")
354-
# return
355-
356354
return GradyKippTensorDamageOwen
357355

358356
#-------------------------------------------------------------------------------

‎src/Damage/ProbabilisticDamageModel.cc

+9-15
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ ProbabilisticDamageModel<Dimension>::
101101
}
102102

103103
//------------------------------------------------------------------------------
104-
// initializeProblemStartup
104+
// initializeProblemStartupDependencies
105105
//
106106
// After all initial state has been initialize (node positions, masses, etc),
107107
// but before we try to run any physics cycles. This is when we initialize a
@@ -110,7 +110,9 @@ ProbabilisticDamageModel<Dimension>::
110110
template<typename Dimension>
111111
void
112112
ProbabilisticDamageModel<Dimension>::
113-
initializeProblemStartup(DataBase<Dimension>& dataBase) {
113+
initializeProblemStartupDependencies(DataBase<Dimension>& dataBase,
114+
State<Dimension>& state,
115+
StateDerivatives<Dimension>& derivs) {
114116

115117
// How many points are actually being damaged?
116118
// We have to be careful to use an unsigned size_t here due to overflow
@@ -132,8 +134,11 @@ initializeProblemStartup(DataBase<Dimension>& dataBase) {
132134
// Compute the initial volumes and random seeds for each node. We hash the
133135
// morton index of each point with the global seed value to create unique
134136
// but reproducible seeds for each points random number generator.
135-
const auto& mass = nodes.mass();
136-
const auto& rho = nodes.massDensity();
137+
auto buildKey = [&](const std::string& fkey) -> std::string { return StateBase<Dimension>::buildFieldKey(fkey, nodes.name()); };
138+
const auto& mass = state.field(buildKey(HydroFieldNames::mass), 0.0);
139+
const auto& rho = (state.registered(buildKey(SolidFieldNames::porositySolidDensity)) ?
140+
state.field(buildKey(SolidFieldNames::porositySolidDensity), 0.0) :
141+
state.field(buildKey(HydroFieldNames::massDensity), 0.0));
137142
const auto nlocal = nodes.numInternalNodes();
138143
vector<uniform_random> randomGenerators(nlocal);
139144
#pragma omp parallel for
@@ -213,17 +218,6 @@ initializeProblemStartup(DataBase<Dimension>& dataBase) {
213218
<< " Avg Neff/Nflaws : " << numFlawsRatio << endl;
214219
}
215220
}
216-
}
217-
218-
//------------------------------------------------------------------------------
219-
// On problem start up, we need to initialize our internal data.
220-
//------------------------------------------------------------------------------
221-
template<typename Dimension>
222-
void
223-
ProbabilisticDamageModel<Dimension>::
224-
initializeProblemStartupDependencies(DataBase<Dimension>& dataBase,
225-
State<Dimension>& state,
226-
StateDerivatives<Dimension>& derivs) {
227221

228222
// Set the moduli.
229223
updateStateFields(HydroFieldNames::pressure, state, derivs);

‎src/Damage/ProbabilisticDamageModel.hh

-6
Original file line numberDiff line numberDiff line change
@@ -50,12 +50,6 @@ public:
5050
//............................................................................
5151
// Override the Physics package interface.
5252

53-
// An optional hook to initialize once when the problem is starting up.
54-
// Typically this is used to size arrays once all the materials and NodeLists have
55-
// been created. It is assumed after this method has been called it is safe to
56-
// call Physics::registerState for instance to create full populated State objects.
57-
virtual void initializeProblemStartup(DataBase<Dimension>& dataBase) override;
58-
5953
// A second optional method to be called on startup, after Physics::initializeProblemStartup has
6054
// been called.
6155
// One use for this hook is to fill in dependendent state using the State object, such as

‎src/Damage/TensorDamageModel.cc

+24
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,30 @@ TensorDamageModel(SolidNodeList<Dimension>& nodeList,
7373
mDamageInCompression(damageInCompression) {
7474
}
7575

76+
//------------------------------------------------------------------------------
77+
// Constructor.
78+
//------------------------------------------------------------------------------
79+
template<typename Dimension>
80+
TensorDamageModel<Dimension>::
81+
TensorDamageModel(SolidNodeList<Dimension>& nodeList,
82+
const TensorStrainAlgorithm strainAlgorithm,
83+
const DamageCouplingAlgorithm damageCouplingAlgorithm,
84+
const TableKernel<Dimension>& W,
85+
const double crackGrowthMultiplier,
86+
const double criticalDamageThreshold,
87+
const bool damageInCompression):
88+
DamageModel<Dimension>(nodeList, W, crackGrowthMultiplier, damageCouplingAlgorithm),
89+
mFlaws(SolidFieldNames::flaws, nodeList),
90+
mYoungsModulus(SolidFieldNames::YoungsModulus, nodeList),
91+
mLongitudinalSoundSpeed(SolidFieldNames::longitudinalSoundSpeed, nodeList),
92+
mStrain(SolidFieldNames::strainTensor, nodeList),
93+
mEffectiveStrain(SolidFieldNames::effectiveStrainTensor, nodeList),
94+
mDdamageDt(TensorDamagePolicy<Dimension>::prefix() + SolidFieldNames::scalarDamage, nodeList),
95+
mStrainAlgorithm(strainAlgorithm),
96+
mCriticalDamageThreshold(criticalDamageThreshold),
97+
mDamageInCompression(damageInCompression) {
98+
}
99+
76100
//------------------------------------------------------------------------------
77101
// Destructor.
78102
//------------------------------------------------------------------------------

‎src/Damage/TensorDamageModel.hh

+9
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,13 @@ public:
7373
const double criticalDamageThreshold,
7474
const bool damageInCompression,
7575
const FlawStorageType& flaws);
76+
TensorDamageModel(SolidNodeList<Dimension>& nodeList,
77+
const TensorStrainAlgorithm strainAlgorithm,
78+
const DamageCouplingAlgorithm damageCouplingAlgorithm,
79+
const TableKernel<Dimension>& W,
80+
const double crackGrowthMultiplier,
81+
const double criticalDamageThreshold,
82+
const bool damageInCompression);
7683
virtual ~TensorDamageModel();
7784

7885
//...........................................................................
@@ -137,6 +144,8 @@ public:
137144
const FlawStorageType& flaws() const;
138145
FlawStorageType& flaws();
139146

147+
void flaws(const FlawStorageType& x);
148+
140149
// The algorithms to update the strain.
141150
TensorStrainAlgorithm strainAlgorithm() const;
142151

‎src/Damage/TensorDamageModelInline.hh

+8
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,14 @@ flaws() {
6767
return mFlaws;
6868
}
6969

70+
template<typename Dimension>
71+
inline
72+
void
73+
TensorDamageModel<Dimension>::
74+
flaws(const typename TensorDamageModel<Dimension>::FlawStorageType& x) {
75+
mFlaws = x;
76+
}
77+
7078
//------------------------------------------------------------------------------
7179
// The strain update algorithm.
7280
//------------------------------------------------------------------------------

‎src/Damage/weibullFlawDistributionBenzAsphaug.cc

+16-9
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,9 @@
1515
#include "NodeList/FluidNodeList.hh"
1616
#include "Field/Field.hh"
1717
#include "Field/FieldList.hh"
18-
#include "DataBase/DataBase.hh"
18+
#include "Hydro/HydroFieldNames.hh"
19+
#include "Strength/SolidFieldNames.hh"
20+
#include "DataBase/State.hh"
1921
#include "Distributed/Communicator.hh"
2022
#include "Utilities/allReduce.hh"
2123

@@ -45,6 +47,7 @@ weibullFlawDistributionBenzAsphaug(double volume,
4547
const double kWeibull,
4648
const double mWeibull,
4749
const FluidNodeList<Dimension>& nodeList,
50+
const State<Dimension>& state,
4851
const int minFlawsPerNode,
4952
const int minTotalFlaws,
5053
const Field<Dimension, int>& mask) {
@@ -58,8 +61,6 @@ weibullFlawDistributionBenzAsphaug(double volume,
5861
REQUIRE(minTotalFlaws > 0);
5962
REQUIRE(mask.nodeListPtr() == &nodeList);
6063

61-
typedef typename Dimension::Scalar Scalar;
62-
6364
// Prepare the result.
6465
Field<Dimension, vector<double> > flaws("Weibull flaw distribution",
6566
nodeList);
@@ -70,7 +71,9 @@ weibullFlawDistributionBenzAsphaug(double volume,
7071

7172
// Prepare a table to faciliate looking local IDs from global.
7273
unordered_map<unsigned, unsigned> global2local;
73-
for (unsigned i = 0; i != nodeList.numInternalNodes(); ++i) global2local[globalIDs(i)] = i;
74+
const auto nlocal = nodeList.numInternalNodes();
75+
#pragma omp parallel for
76+
for (auto i = 0u; i < nlocal; i++) global2local[globalIDs(i)] = i;
7477
CHECK(global2local.size() == nodeList.numInternalNodes());
7578

7679
// Prepare an int per *each* node, so that each process can keep track of how many
@@ -84,11 +87,15 @@ weibullFlawDistributionBenzAsphaug(double volume,
8487
if (n > 0) {
8588

8689
// If the user did not speicify a volume, we compute it from the information
87-
// in the NodeList.
90+
// in the State.
8891
if (volume == 0.0) {
89-
const Field<Dimension, Scalar>& mass = nodeList.mass();
90-
const Field<Dimension, Scalar>& rho = nodeList.massDensity();
91-
for (auto i = 0u; i != nodeList.numInternalNodes(); ++i) {
92+
auto buildKey = [&](const std::string& fkey) -> std::string { return StateBase<Dimension>::buildFieldKey(fkey, nodeList.name()); };
93+
const auto& mass = state.field(buildKey(HydroFieldNames::mass), 0.0);
94+
const auto& rho = (state.registered(buildKey(SolidFieldNames::porositySolidDensity)) ?
95+
state.field(buildKey(SolidFieldNames::porositySolidDensity), 0.0) :
96+
state.field(buildKey(HydroFieldNames::massDensity), 0.0));
97+
#pragma omp parallel for
98+
for (auto i = 0u; i < nlocal; i++) {
9299
CHECK(rho(i) > 0.0);
93100
volume += mass(i)/rho(i);
94101
}
@@ -148,7 +155,7 @@ weibullFlawDistributionBenzAsphaug(double volume,
148155
unsigned totalNumFlaws = 0;
149156
double epsMax = 0.0;
150157
double sumFlaws = 0.0;
151-
for (auto i = 0u; i != nodeList.numInternalNodes(); ++i) {
158+
for (auto i = 0u; i < nodeList.numInternalNodes(); ++i) {
152159
minNumFlaws = min(minNumFlaws, unsigned(flaws(i).size()));
153160
maxNumFlaws = max(maxNumFlaws, unsigned(flaws(i).size()));
154161
totalNumFlaws += flaws(i).size();

‎src/Damage/weibullFlawDistributionBenzAsphaug.hh

+2
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111
namespace Spheral {
1212
template<typename Dimension> class FluidNodeList;
1313
template<typename Dimension, typename DataType> class Field;
14+
template<typename Dimension> class State;
1415
}
1516

1617
namespace Spheral {
@@ -27,6 +28,7 @@ weibullFlawDistributionBenzAsphaug(double volume,
2728
const double kWeibull,
2829
const double mWeibull,
2930
const FluidNodeList<Dimension>& nodeList,
31+
const State<Dimension>& state,
3032
const int minFlawsPerNode,
3133
const int minTotalFlaws,
3234
const Field<Dimension, int>& mask);

‎src/Damage/weibullFlawDistributionBenzAsphaugInst.cc.py

+1
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
const double,
1313
const double,
1414
const Spheral::FluidNodeList<Spheral::Dim< %(ndim)s > >&,
15+
const Spheral::State<Spheral::Dim< %(ndim)s > >&,
1516
const int,
1617
const int,
1718
const Spheral::Field<Spheral::Dim<%(ndim)s>, int>&);

0 commit comments

Comments
 (0)