Skip to content

Commit 43c4d4f

Browse files
mcfadden8mcfadden8
authored
Merge pull request #488 from llnl/pr-474_476_477_478_480-combo
2 parents e125c29 + 0e3db62 commit 43c4d4f

20 files changed

Lines changed: 373 additions & 105 deletions

‎RELEASE_NOTES.md‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -550,7 +550,19 @@ Version vYYYY.MM.p -- Release date YYYY-MM-DD
550550
Notable changes include:
551551
552552
* New features / API changes:
553+
* Added verbose time step frequency to `SpheralController`, allowing time step information to be printed every N steps instead of every step.
554+
* Periodic work frequencies in `SpheralController` can now be callables, enabling dynamic frequency changes during a simulation.
555+
* Added `gradientPairs` to modernize the gradient calculation using node pairs.
556+
* Silo output now supports subdirectories within silo files.
557+
* Added support for 1D silo output, readable as curves in VisIt.
558+
* Added a damping factor to the ideal H iteration in both SPH and ASPH smoothing scales to improve convergence for difficult initial distributions.
559+
* Added a check in `SpheralController` that limits the initial neighbor count to a reasonable value before `iterateIdealH` begins.
560+
* `iterateIdealH` now supports extra packages during the initial H iteration and a flag to force Voronoi tessellation.
553561
554562
* Build changes / improvements:
563+
* Cleaned up the gradient instantiation template.
555564
556565
* Bug Fixes / improvements:
566+
* Moved axis boundary logic out of individual SPH and CRKSPH hydro constructors into `SpheralController`, where it is applied once during initialization. This allows the axis BC to work without hydro.
567+
* Added `allReduceLoc` utility so the controlling time step is printed once with its owning rank, rather than once per rank when time steps do not vary by processor.
568+

‎src/CRKSPH/CRKSPHHydros.py‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -104,8 +104,7 @@ def CRKSPH(dataBase,
104104

105105
# If we're using area-weighted RZ, we need to reflect from the axis
106106
if GeometryRegistrar.coords() == CoordinateType.RZ:
107-
result.zaxisBC = AxisBoundaryRZ(etaMinAxis)
108-
result.appendBoundary(result.zaxisBC)
107+
result.etaMinAxis = etaMinAxis
109108

110109
return result
111110

‎src/Distributed/allReduce.hh‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,22 +27,45 @@ namespace Spheral {
2727
#define SPHERAL_OP_PROD MPI_PROD
2828
#define SPHERAL_OP_LAND MPI_LAND
2929
#define SPHERAL_OP_LOR MPI_LOR
30+
#define SPHERAL_OP_MINLOC MPI_MINLOC
31+
#define SPHERAL_OP_MAXLOC MPI_MAXLOC
3032

3133
template<typename Value>
3234
constexpr Value
3335
allReduce(const Value& value, const MPI_Op op,
3436
const MPI_Comm comm = Communicator::communicator()) {
37+
CHECK(!(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC));
3538
Value tmp = value;
3639
Value result;
3740
MPI_Allreduce(&tmp, &result, 1,
3841
DataTypeTraits<Value>::MpiDataType(), op, comm);
3942
return result;
4043
}
4144

45+
template<typename Value>
46+
constexpr std::pair<Value, int>
47+
allReduceLoc(const Value value, const MPI_Op op,
48+
const MPI_Comm comm = Communicator::communicator()) {
49+
CHECK(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC);
50+
struct {
51+
Value val;
52+
int rank;
53+
} in, out;
54+
55+
MPI_Comm_rank(comm, &in.rank);
56+
in.val = value;
57+
58+
MPI_Allreduce(&in, &out, 1, DataTypeTraits<Value>::MpiLocDataType(), op, comm);
59+
60+
return {out.val, out.rank};
61+
}
62+
63+
4264
template<typename Value>
4365
constexpr Value
4466
distScan(const Value& value, const MPI_Op op,
4567
const MPI_Comm comm = Communicator::communicator()) {
68+
CHECK(!(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC));
4669
Value tmp = value;
4770
Value result;
4871
MPI_Scan(&tmp, &result, 1, DataTypeTraits<Value>::MpiDataType(), op, comm);
@@ -65,13 +88,22 @@ Barrier(const MPI_Comm comm = Communicator::communicator()) {
6588
#define SPHERAL_OP_PROD 4
6689
#define SPHERAL_OP_LAND 5
6790
#define SPHERAL_OP_LOR 6
91+
#define SPHERAL_OP_MINLOC 7
92+
#define SPHERAL_OP_MAXLOC 8
6893

6994
template<typename Value>
7095
constexpr Value
7196
allReduce(const Value& value, const int /*op*/, const int = 0) {
7297
return value;
7398
}
7499

100+
template<typename Value>
101+
inline std::pair<Value, int>
102+
allReduceLoc(const Value value, const int /*op*/,
103+
const int = 0) {
104+
return {value, 0};
105+
}
106+
75107
template<typename Value>
76108
constexpr Value
77109
distScan(const Value& value, const int /*op*/, const int = 0) {

‎src/FieldOperations/FieldListFunctions.hh‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88

99
namespace Spheral {
1010

11+
template<typename Dimension> class ConnectivityMap;
1112
template<typename Dimension> class TableKernel;
1213

1314
template<typename Dimension, typename DataType> class FieldList;
@@ -44,6 +45,16 @@ gradient(const FieldList<Dimension, std::vector<DataType>>& fieldList,
4445
const FieldList<Dimension, typename Dimension::SymTensor>& Hfield,
4546
const TableKernel<Dimension>& kernel);
4647

48+
template<typename Dimension, typename DataType>
49+
void
50+
gradientPairs(FieldList<Dimension, typename MathTraits<Dimension, DataType>::GradientType>& result,
51+
const FieldList<Dimension, DataType>& field,
52+
const FieldList<Dimension, typename Dimension::Vector>& position,
53+
const FieldList<Dimension, typename Dimension::Scalar>& weight,
54+
const FieldList<Dimension, typename Dimension::SymTensor>& H,
55+
const ConnectivityMap<Dimension>& conn,
56+
const TableKernel<Dimension>& kernel);
57+
4758
// Calculate the divergence of the given FieldList.
4859
template<typename Dimension, typename DataType>
4960
FieldList<Dimension, typename MathTraits<Dimension, DataType>::DivergenceType>

‎src/FieldOperations/gradient.cc‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "Field/FieldList.hh"
99
#include "Field/Field.hh"
1010
#include "Field/NodeIterators.hh"
11+
#include "Neighbor/ConnectivityMap.hh"
1112
#include "NodeList/NodeList.hh"
1213
#include "Neighbor/Neighbor.hh"
1314
#include "Kernel/TableKernel.hh"
@@ -277,6 +278,60 @@ gradient(const FieldList<Dimension, std::vector<DataType>>& fieldList,
277278
return result;
278279
}
279280

281+
//------------------------------------------------------------------------------
282+
// Calculate the gradient of a FieldList.
283+
//------------------------------------------------------------------------------
284+
template<typename Dimension, typename DataType>
285+
void
286+
gradientPairs(FieldList<Dimension, typename MathTraits<Dimension, DataType>::GradientType>& result,
287+
const FieldList<Dimension, DataType>& field,
288+
const FieldList<Dimension, typename Dimension::Vector>& position,
289+
const FieldList<Dimension, typename Dimension::Scalar>& weight,
290+
const FieldList<Dimension, typename Dimension::SymTensor>& H,
291+
const ConnectivityMap<Dimension>& conn,
292+
const TableKernel<Dimension>& kernel) {
293+
typedef typename MathTraits<Dimension, DataType>::GradientType GradientType;
294+
result = GradientType::zero();
295+
296+
const auto& pairs = conn.nodePairList();
297+
const auto npairs = pairs.size();
298+
299+
for (auto k = 0u; k < npairs; ++k) {
300+
const auto ni = pairs[k].i_list;
301+
const auto nj = pairs[k].j_list;
302+
const auto i = pairs[k].i_node;
303+
const auto j = pairs[k].j_node;
304+
305+
const auto& ri = position(ni, i);
306+
const auto& rj = position(nj, j);
307+
const auto& vi = weight(ni, i);
308+
const auto& vj = weight(nj, j);
309+
const auto& hi = H(ni, i);
310+
const auto& hj = H(nj, j);
311+
const auto& fi = field(ni, i);
312+
const auto& fj = field(nj, j);
313+
314+
const auto rij = ri - rj;
315+
const auto hdeti = hi.Determinant();
316+
const auto hdetj = hj.Determinant();
317+
const auto etai = hi * rij;
318+
const auto etaj = hj * rij;
319+
const auto etaMagi = etai.magnitude();
320+
const auto etaMagj = etaj.magnitude();
321+
const auto etaUniti = etai.unitVector();
322+
const auto etaUnitj = etaj.unitVector();
323+
const auto hetaUniti = hi * etaUniti;
324+
const auto hetaUnitj = hj * etaUnitj;
325+
326+
const auto dwi = hetaUniti * kernel.gradValue(etaMagi, hdeti);
327+
const auto dwj = hetaUnitj * kernel.gradValue(etaMagj, hdetj);
328+
const auto dwij = 0.5 * (dwi + dwj);
329+
330+
result(ni, i) += vj * (fj - fi) * dwij;
331+
result(nj, j) -= vi * (fi - fj) * dwij;
332+
}
333+
}
334+
280335
//------------------------------------------------------------------------------
281336
// The limiter method.
282337
//------------------------------------------------------------------------------

‎src/FieldOperations/gradientInst.cc.py‎

Lines changed: 42 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -6,62 +6,57 @@
66
#include "Geometry/Dimension.hh"
77
88
namespace Spheral {
9+
"""
910

10-
//============================== gradient() ==============================
11-
template
12-
FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>
13-
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& fieldList,
14-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
15-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
16-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
17-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
18-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
19-
const TableKernel< Dim< %(ndim)s > >& kernel);
11+
for Value in ("Scalar", "Vector"):
12+
text += """
2013
template
21-
FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>
22-
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& fieldList,
23-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
24-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
25-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
26-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
27-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
28-
const TableKernel< Dim< %(ndim)s > >& kernel);
14+
FieldList<%%(Dim)s, MathTraits<%%(Dim)s, %(Value)s>::GradientType>
15+
gradient<%%(Dim)s, %(Value)s>(const FieldList<%%(Dim)s, %(Value)s>& fieldList,
16+
const FieldList<%%(Dim)s, %%(Vector)s>& position,
17+
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
18+
const FieldList<%%(Dim)s, %%(Scalar)s>& mass,
19+
const FieldList<%%(Dim)s, %%(Scalar)s>& density,
20+
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
21+
const TableKernel< %%(Dim)s >& kernel);
2922
3023
template
31-
FieldList<Dim< %(ndim)s >, std::vector<MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>>
32-
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, std::vector<Dim< %(ndim)s >::Scalar>>& fieldList,
33-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
34-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
35-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
36-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
37-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
38-
const TableKernel< Dim< %(ndim)s > >& kernel);
39-
template
40-
FieldList<Dim< %(ndim)s >, std::vector<MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>>
41-
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, std::vector<Dim< %(ndim)s >::Vector>>& fieldList,
42-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
43-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
44-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
45-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
46-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
47-
const TableKernel< Dim< %(ndim)s > >& kernel);
24+
FieldList<%%(Dim)s, std::vector<MathTraits<%%(Dim)s, %(Value)s>::GradientType>>
25+
gradient<%%(Dim)s, %(Value)s>(const FieldList<%%(Dim)s, std::vector<%(Value)s>>& fieldList,
26+
const FieldList<%%(Dim)s, %%(Vector)s>& position,
27+
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
28+
const FieldList<%%(Dim)s, %%(Scalar)s>& mass,
29+
const FieldList<%%(Dim)s, %%(Scalar)s>& density,
30+
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
31+
const TableKernel< %%(Dim)s >& kernel);
4832
33+
template
34+
void
35+
gradientPairs<%%(Dim)s, %(Value)s>(FieldList<%%(Dim)s, MathTraits<%%(Dim)s, %(Value)s>::GradientType>& result,
36+
const FieldList<%%(Dim)s, %(Value)s>& field,
37+
const FieldList<%%(Dim)s, %%(Vector)s>& position,
38+
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
39+
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
40+
const ConnectivityMap<%%(Dim)s>& conn,
41+
const TableKernel<%%(Dim)s>& kernel);
42+
""" % {"Value" : "%(" + Value + ")s"}
4943

44+
text += """
5045
//============================== limiter() ==============================
5146
template
52-
FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>
53-
limiter<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& fieldList,
54-
const FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>& gradient,
55-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
56-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
57-
const TableKernel< Dim< %(ndim)s > >& kernel);
47+
FieldList<%(Dim)s, %(SymTensor)s>
48+
limiter<%(Dim)s, %(Scalar)s>(const FieldList<%(Dim)s, %(Scalar)s>& fieldList,
49+
const FieldList<%(Dim)s, MathTraits<%(Dim)s, %(Scalar)s>::GradientType>& gradient,
50+
const FieldList<%(Dim)s, %(Vector)s>& position,
51+
const FieldList<%(Dim)s, %(SymTensor)s>& Hfield,
52+
const TableKernel< %(Dim)s >& kernel);
5853
template
59-
FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>
60-
limiter<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& fieldList,
61-
const FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>& gradient,
62-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
63-
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
64-
const TableKernel< Dim< %(ndim)s > >& kernel);
54+
FieldList<%(Dim)s, %(SymTensor)s>
55+
limiter<%(Dim)s, %(Vector)s>(const FieldList<%(Dim)s, %(Vector)s>& fieldList,
56+
const FieldList<%(Dim)s, MathTraits<%(Dim)s, %(Vector)s>::GradientType>& gradient,
57+
const FieldList<%(Dim)s, %(Vector)s>& position,
58+
const FieldList<%(Dim)s, %(SymTensor)s>& Hfield,
59+
const TableKernel< %(Dim)s >& kernel);
6560
6661
}
6762
"""

‎src/Integrator/Integrator.cc‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@ Integrator(DataBase<Dimension>& dataBase,
4848
mCurrentTime(0.0),
4949
mCurrentCycle(0),
5050
mUpdateBoundaryFrequency(1),
51+
mVerboseStep(1),
5152
mVerbose(false),
5253
mAllowDtCheck(false),
5354
mRequireConnectivity(true),
@@ -150,16 +151,16 @@ selectDt(const typename Dimension::Scalar dtMin,
150151

151152
// In the parallel case we need to find the minimum timestep across all processors.
152153
#ifdef SPHERAL_ENABLE_GLOBALDT_REDUCTION
153-
const auto globalDt = allReduce(dt.first, SPHERAL_OP_MIN);
154+
const auto [globalDt, rank] = allReduceLoc(dt.first, SPHERAL_OP_MINLOC);
154155
#else
155156
const auto globalDt = dt.first;
157+
const auto rank = Process::getRank();
156158
#endif
157159

158160
// Are we verbose?
159-
if (dt.first == globalDt and
160-
(verbose() or globalDt < mDtMin)) {
161+
if (rank == Process::getRank() and (verbose() and currentCycle() % verboseStep() == 0)) {
161162
std::cout << "Selected timestep of "
162-
<< dt.first << std::endl
163+
<< dt.first << " on rank " << rank << std::endl
163164
<< dt.second << std::endl;
164165
}
165166
std::cout.flush();

‎src/Integrator/Integrator.hh‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -192,6 +192,10 @@ public:
192192
bool verbose() const { return mVerbose; }
193193
void verbose(const bool x) { mVerbose = x; }
194194

195+
// Select whether the integrator is verbose or not during a cycle.
196+
int verboseStep() const { return mVerboseStep; }
197+
void verboseStep(const int x) { mVerboseStep = x; }
198+
195199
// Should the integrator check interim timestep votes and abort steps?
196200
bool allowDtCheck() const { return mAllowDtCheck; }
197201
void allowDtCheck(const bool x) { mAllowDtCheck = x; }
@@ -235,7 +239,7 @@ protected:
235239
private:
236240
//--------------------------- Private Interface ---------------------------//
237241
Scalar mDtMin, mDtMax, mDtGrowth, mLastDt, mDtCheckFrac, mCurrentTime;
238-
int mCurrentCycle, mUpdateBoundaryFrequency;
242+
int mCurrentCycle, mUpdateBoundaryFrequency, mVerboseStep;
239243
bool mVerbose, mAllowDtCheck, mRequireConnectivity, mRequireGhostConnectivity, mRequireOverlapConnectivity, mRequireIntersectionConnectivity;
240244
std::reference_wrapper<DataBase<Dimension>> mDataBase;
241245
std::vector<Physics<Dimension>*> mPhysicsPackages;

0 commit comments

Comments
 (0)