Peano 4
Loading...
Searching...
No Matches
MassAccumulator.cpp
Go to the documentation of this file.
2
3#include <cmath>
4
5#include "tarch/Assertions.h"
6#include "tarch/la/Vector.h"
7#include "tarch/mpi/Rank.h"
11
12
14 "applications::exahype2::euler::sphericalaccretion::MassAccumulator"
15);
16
17
19 double shellWidth, double threshold
20):
21 _shellWidth(shellWidth),
22 _threshold(threshold) {}
23
24
26 return (int)(std::floor(radius / _shellWidth));
27}
28
29
31 tarch::multicore::Lock lock(_semaphore);
32 int bucket = getBucket(radius);
33 while (_currentMass.size() <= bucket) {
34 _currentMass.push_back(0.0);
35 }
36 _currentMass[bucket] += mass;
37}
38
39
41 double radius
42) const {
43 int bucket = getBucket(radius);
44 if (_previousMassAccumulated.empty()) {
45 return 0.0;
46 } else if (_previousMassAccumulated.size() <= bucket) {
47 return _previousMassAccumulated.at(_previousMassAccumulated.size() - 1);
48 } else {
49 return _previousMassAccumulated.at(bucket);
50 }
51}
52
53
55) const {
56 int bucket = getBucket(radius);
57
58 if (_previousMassAccumulated.empty()) {
59 return 0.0;
60 } else if (_previousMassAccumulated.size() <= bucket) {
61 return _previousMassAccumulated.back();
62 } else {
63 int leftBucket = bucket;
64 int rightBucket = bucket;
65 double rightWeight = 1.0;
66
67 if (radius - bucket * _shellWidth < _shellWidth / 2.0 and bucket > 0) {
68 leftBucket--;
69 rightWeight = (radius - (leftBucket + 0.5) * _shellWidth) / _shellWidth;
70 } else if (radius - bucket * _shellWidth < _shellWidth / 2.0 and bucket == 0) {
71 rightWeight = 1.0;
72 } else if (radius - bucket * _shellWidth >= _shellWidth / 2.0 and bucket < static_cast<int>(_previousMassAccumulated.size()) - 1) {
73 rightBucket++;
74 rightWeight = (radius - (leftBucket + 0.5) * _shellWidth) / _shellWidth;
75 }
76 else if (radius-bucket*_shellWidth >= _shellWidth/2.0 and bucket==static_cast<int>(_previousMassAccumulated.size())-1) {
77 rightWeight = 0.0;
78 } else {
79 assertion(false);
80 }
81
83 tarch::la::greaterEquals(rightWeight, 0.0), rightWeight, leftBucket, rightBucket, bucket, radius, _shellWidth
84 );
86 tarch::la::smallerEquals(rightWeight, 1.0), rightWeight, leftBucket, rightBucket, bucket, radius, _shellWidth
87 );
88
89 return rightWeight * _previousMassAccumulated.at(rightBucket)
90 + (1.0 - rightWeight) * _previousMassAccumulated.at(leftBucket);
91 }
92}
93
94
96#ifdef Parallel
97 int globalMaxNumberOfBuckets = _currentMass.size();
98 int localMaxNumberOfBuckets = _currentMass.size();
99 logDebug("finishAccumulation()", "have " << localMaxNumberOfBuckets << " shells/buckets on rank");
101 .allReduce(&localMaxNumberOfBuckets, &globalMaxNumberOfBuckets, 1, MPI_INT, MPI_MAX, [&]() -> void {
103 });
104 logDebug("finishAccumulation()", "have " << globalMaxNumberOfBuckets << " shells/buckets globally");
105
106 assertion2(globalMaxNumberOfBuckets >= _currentMass.size(), globalMaxNumberOfBuckets, _currentMass.size());
107
108 double* globalBucketValues = new double[globalMaxNumberOfBuckets];
109 double* localBucketValues = new double[globalMaxNumberOfBuckets];
110 for (int bucket = 0; bucket < globalMaxNumberOfBuckets; bucket++) {
111 localBucketValues[bucket] = bucket < static_cast<int>(_currentMass.size()) ? _currentMass[bucket] : 0.0;
112 globalBucketValues[bucket] = bucket < static_cast<int>(_currentMass.size()) ? _currentMass[bucket] : 0.0;
113 logDebug("finishAccumulation()", "bucket " << bucket << " hosts value " << globalBucketValues[bucket]);
114 }
115
117 .allReduce(localBucketValues, globalBucketValues, globalMaxNumberOfBuckets, MPI_DOUBLE, MPI_SUM, [&]() -> void {
119 });
120
121 for (int bucket = 0; bucket < globalMaxNumberOfBuckets; bucket++) {
122 if (bucket < _currentMass.size()) {
123 _currentMass[bucket] = globalBucketValues[bucket];
124 } else {
125 _currentMass.push_back(globalBucketValues[bucket]);
126 }
127 }
128 delete[] globalBucketValues;
129 delete[] localBucketValues;
130#endif
131
132 logDebug("finishAccumulation()", "old relevant radius equals r=" << getMaxRelevantRadius());
133
134 _previousMassAccumulated.clear();
135 _previousMassAccumulated.insert(_previousMassAccumulated.end(), _currentMass.begin(), _currentMass.end());
136 _currentMass.clear();
137
138 // accumulate
139 for (int i = 1; i < static_cast<int>(_previousMassAccumulated.size()); i++) {
140 _previousMassAccumulated[i] += _previousMassAccumulated[i - 1];
141 }
142
143 int largestRelevantShell = 0;
144 for (int i = 0; i < static_cast<int>(_previousMassAccumulated.size()) - 1; i++) {
145 if (not tarch::la::equals(_previousMassAccumulated[i], 0.0) and std::abs(_previousMassAccumulated[i + 1] / _previousMassAccumulated[i]) >= 1.0 + _threshold) {
146 largestRelevantShell = std::max(largestRelevantShell, i + 1);
147 }
148 }
149 largestRelevantShell = std::max(
150 0, std::min(largestRelevantShell + 1, static_cast<int>(_previousMassAccumulated.size()))
151 );
152 logDebug(
153 "finishAccumulation()",
154 "roll over "
155 << largestRelevantShell << " out of " << _previousMassAccumulated.size()
156 << " buckets (remaining buckets make no significant contribution)"
157 );
158 _previousMassAccumulated.resize(largestRelevantShell);
159
160
161 double currentShellWidth = 0.0;
162 for (auto& p : _previousMassAccumulated) {
163 currentShellWidth += _shellWidth;
165 logInfo("finishAccumulation()", "bucket (r<=" << currentShellWidth << ") = " << p);
166 }
167 }
168 if (tarch::mpi::Rank::getInstance().isGlobalMaster()) {
169 logInfo("finishAccumulation()", "new relevant radius equals r=" << getMaxRelevantRadius());
170 }
171}
172
173
175 return _previousMassAccumulated.empty()
176 ? std::numeric_limits<double>::max()
177 : (_previousMassAccumulated.size() + 1.0) * _shellWidth;
178}
179
180
182 return _previousMassAccumulated.empty() ? 0.0 : _previousMassAccumulated.back();
183}
#define assertion2(expr, param0, param1)
#define assertion(expr)
#define assertion6(expr, param0, param1, param2, param3, param4, param5)
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#define logInfo(methodName, logMacroMessageStream)
Wrapper macro around tarch::tarch::logging::Log to improve logging.
Definition Log.h:411
int getBucket(double radius) const
Identifies the right bucket.
double getMaxRelevantRadius() const
Return the biggest radius which contributes towards mass.
void finishAccumulation()
Finalise data accumulation for next time step.
double getMass_piecewiseConstantInterpolation(double radius) const
If you try to read the mass for a distance where we had no accumulation yet, then this routine return...
double getTotalMass() const
This routine returns the total mass.
MassAccumulator(double shellWidth=0.01, double threshold=0.01)
Log Device.
Definition Log.h:516
bool isGlobalMaster() const
Is this node the global master process, i.e.
Definition Rank.cpp:419
static Rank & getInstance()
This operation returns the singleton instance.
Definition Rank.cpp:538
void allReduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, std::function< void()> waitor=[]() -> void {})
Definition Rank.cpp:279
Create a lock around a boolean semaphore region.
Definition Lock.h:19
virtual void receiveDanglingMessages() override
Answer to MPI Messages.
static ServiceRepository & getInstance()
bool greaterEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.