Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #ifdef ENABLE_NEK_COUPLING
20 :
21 : #include "NekSpatialBinUserObject.h"
22 : #include "NekInterface.h"
23 : #include "CardinalUtils.h"
24 :
25 : InputParameters
26 433 : NekSpatialBinUserObject::validParams()
27 : {
28 433 : InputParameters params = GeneralUserObject::validParams();
29 433 : params += NekBase::validParams();
30 433 : params += NekFieldInterface::validParams();
31 866 : params.addRequiredParam<std::vector<UserObjectName>>(
32 : "bins", "Userobjects providing a spatial bin given a point");
33 866 : params.addParam<unsigned int>(
34 : "interval",
35 866 : 1,
36 : "Frequency (in number of time steps) with which to execute this user object; user objects "
37 : "can be expensive and not necessary to evaluate on every single time step. NOTE: you "
38 : "probably want to match with 'time_step_interval' in the Output");
39 866 : params.addParam<bool>(
40 : "map_space_by_qp",
41 866 : false,
42 : "Whether to map the NekRS spatial domain to a bin according to the element centroids (true) "
43 : "or quadrature point locations (false).");
44 866 : params.addParam<bool>(
45 : "check_zero_contributions",
46 866 : true,
47 : "Whether to throw an error if no GLL points/element centroids in the NekRS mesh map to a "
48 : "spatial bin; this "
49 : "can be used to ensure that the bins are sufficiently big to get at least one contributing "
50 : "point from the NekRS mesh.");
51 433 : return params;
52 0 : }
53 :
54 220 : NekSpatialBinUserObject::NekSpatialBinUserObject(const InputParameters & parameters)
55 : : GeneralUserObject(parameters),
56 : NekBase(this, parameters),
57 : NekFieldInterface(this, parameters, true /* allow normal */),
58 219 : _interval(getParam<unsigned int>("interval")),
59 438 : _bin_names(getParam<std::vector<UserObjectName>>("bins")),
60 438 : _map_space_by_qp(getParam<bool>("map_space_by_qp")),
61 438 : _check_zero_contributions(getParam<bool>("check_zero_contributions")),
62 219 : _bin_values(nullptr),
63 219 : _bin_values_x(nullptr),
64 219 : _bin_values_y(nullptr),
65 219 : _bin_values_z(nullptr),
66 219 : _bin_volumes(nullptr),
67 219 : _bin_counts(nullptr),
68 219 : _bin_partial_values(nullptr),
69 220 : _bin_partial_counts(nullptr)
70 : {
71 219 : _fixed_mesh = !nekrs::hasMovingMesh();
72 :
73 219 : if (_bin_names.size() == 0)
74 0 : paramError("bins", "Length of vector must be greater than zero!");
75 :
76 570 : for (auto & b : _bin_names)
77 : {
78 : // first check that the user object exists
79 353 : if (!hasUserObjectByName<UserObject>(b))
80 2 : mooseError("Bin user object with name '" + b +
81 : "' not found in problem. The user objects "
82 1 : "in 'bins' must be listed before the '" +
83 0 : name() + "' user object.");
84 :
85 : // then check that it's the right type
86 352 : if (!hasUserObjectByName<SpatialBinUserObject>(b))
87 2 : mooseError("Bin user object with name '" + b +
88 : "' must inherit from SpatialBinUserObject.\n\n"
89 : "Volume options: HexagonalSubchannelBin, LayeredBin, RadialBin\n"
90 : "Side options: HexagonalSubchannelGapBin, LayeredGapBin");
91 :
92 351 : _bins.push_back(&getUserObjectByName<SpatialBinUserObject>(b));
93 : }
94 :
95 217 : _n_bins = num_bins();
96 :
97 217 : _bin_values = (double *)calloc(_n_bins, sizeof(double));
98 217 : _bin_volumes = (double *)calloc(_n_bins, sizeof(double));
99 217 : _bin_partial_values = (double *)calloc(_n_bins, sizeof(double));
100 217 : _bin_counts = (int *)calloc(_n_bins, sizeof(int));
101 217 : _bin_partial_counts = (int *)calloc(_n_bins, sizeof(int));
102 :
103 217 : if (_field == field::velocity_component)
104 : {
105 17 : _bin_values_x = (double *)calloc(_n_bins, sizeof(double));
106 17 : _bin_values_y = (double *)calloc(_n_bins, sizeof(double));
107 17 : _bin_values_z = (double *)calloc(_n_bins, sizeof(double));
108 : }
109 :
110 217 : _has_direction = {false, false, false};
111 217 : _bin_providing_direction.resize(3);
112 566 : for (unsigned int b = 0; b < _bins.size(); ++b)
113 : {
114 : const auto & bin = _bins[b];
115 :
116 : // directions provided by this bin
117 351 : auto bin_directions = bin->directions();
118 :
119 820 : for (const auto & d : bin_directions)
120 : {
121 471 : if (_has_direction[d])
122 : {
123 2 : const auto & bin_providing_d = _bins[_bin_providing_direction[d]];
124 2 : mooseError("Cannot combine multiple distributions in the same coordinate direction!\n"
125 2 : "Bin '" +
126 4 : bin->name() + "' conflicts with bin '" + bin_providing_d->name() + "'.");
127 : }
128 : else
129 : {
130 : _has_direction[d] = true;
131 469 : _bin_providing_direction[d] = b;
132 : }
133 : }
134 349 : }
135 :
136 : // initialize all points to (0, 0, 0)
137 44857 : for (unsigned int i = 0; i < _n_bins; ++i)
138 44642 : _points.push_back(Point(0.0, 0.0, 0.0));
139 :
140 : // we will at most have 3 separate distributions
141 215 : if (_bins.size() == 1)
142 102 : computePoints1D();
143 113 : else if (_bins.size() == 2)
144 96 : computePoints2D();
145 : else
146 17 : computePoints3D();
147 :
148 : // with a user-specified direction, the direction for each bin is the same
149 215 : if (_field == field::velocity_component && _velocity_component == component::user)
150 1016 : for (unsigned int i = 0; i < _n_bins; ++i)
151 1008 : _velocity_bin_directions.push_back(_velocity_direction);
152 215 : }
153 :
154 208 : NekSpatialBinUserObject::~NekSpatialBinUserObject()
155 : {
156 208 : freePointer(_bin_values);
157 208 : freePointer(_bin_volumes);
158 208 : freePointer(_bin_counts);
159 208 : freePointer(_bin_partial_values);
160 208 : freePointer(_bin_partial_counts);
161 :
162 208 : freePointer(_bin_values_x);
163 208 : freePointer(_bin_values_y);
164 208 : freePointer(_bin_values_z);
165 416 : }
166 :
167 : void
168 1432 : NekSpatialBinUserObject::execute()
169 : {
170 1432 : if (_fe_problem.timeStep() % _interval == 0)
171 216 : executeUserObject();
172 1432 : }
173 :
174 : Point
175 126587040 : NekSpatialBinUserObject::nekPoint(const int & local_elem_id, const int & local_node_id) const
176 : {
177 126587040 : if (_map_space_by_qp)
178 99589536 : return nekrs::gllPoint(local_elem_id, local_node_id);
179 : else
180 26997504 : return nekrs::centroid(local_elem_id);
181 : }
182 :
183 : void
184 459 : NekSpatialBinUserObject::resetPartialStorage()
185 : {
186 59638 : for (unsigned int i = 0; i < _n_bins; ++i)
187 : {
188 59179 : _bin_partial_values[i] = 0.0;
189 59179 : _bin_partial_counts[i] = 0;
190 : }
191 459 : }
192 :
193 : void
194 211 : NekSpatialBinUserObject::computeBinVolumes()
195 : {
196 211 : getBinVolumes();
197 :
198 211 : if (_check_zero_contributions)
199 : {
200 11327 : for (unsigned int i = 0; i < _n_bins; ++i)
201 : {
202 11142 : if (_bin_counts[i] == 0)
203 : {
204 3 : std::string map = _map_space_by_qp ? "GLL points" : "element centroids";
205 4 : mooseError(
206 4 : "Failed to map any " + map + " to bin " + Moose::stringify(i) +
207 : "!\n\n"
208 : "This can happen if the bins are much finer than the NekRS mesh or if the bins are "
209 : "defined in a way that results in bins entirely outside the NekRS domain. You can turn "
210 : "this error off by setting 'check_zero_contributions = false' (at your own risk!).");
211 : }
212 : }
213 : }
214 209 : }
215 :
216 : Real
217 1671078 : NekSpatialBinUserObject::spatialValue(const Point & p) const
218 : {
219 1671078 : return _bin_values[bin(p)];
220 : }
221 :
222 : const std::vector<unsigned int>
223 97088 : NekSpatialBinUserObject::unrolledBin(const unsigned int & total_bin_index) const
224 : {
225 : std::vector<unsigned int> local_bins;
226 97088 : local_bins.resize(_bins.size());
227 :
228 97088 : int running_index = total_bin_index;
229 269800 : for (int i = _bins.size() - 1; i >= 0; --i)
230 : {
231 172712 : local_bins[i] = running_index % _bins[i]->num_bins();
232 172712 : running_index -= local_bins[i];
233 172712 : running_index /= _bins[i]->num_bins();
234 : }
235 :
236 97088 : return local_bins;
237 0 : }
238 :
239 : const unsigned int
240 59791362 : NekSpatialBinUserObject::bin(const Point & p) const
241 : {
242 : // get the indices into each of the individual bin objects
243 : std::vector<unsigned int> indices;
244 169139690 : for (const auto & b : _bins)
245 109348328 : indices.push_back(b->bin(p));
246 :
247 : // convert to a total index into the multidimensional bin union
248 59791362 : unsigned int index = indices[0];
249 109348328 : for (unsigned int i = 1; i < _bins.size(); ++i)
250 49556966 : index = index * _bins[i]->num_bins() + indices[i];
251 :
252 59791362 : return index;
253 59791362 : }
254 :
255 : const unsigned int
256 8361 : NekSpatialBinUserObject::num_bins() const
257 : {
258 : unsigned int num_bins = 1;
259 25020 : for (const auto & b : _bins)
260 16659 : num_bins *= b->num_bins();
261 :
262 8361 : return num_bins;
263 : }
264 :
265 : void
266 102 : NekSpatialBinUserObject::computePoints1D()
267 : {
268 574 : for (unsigned int i = 0; i < _bins[0]->num_bins(); ++i)
269 : {
270 472 : std::vector<unsigned int> indices = {i};
271 472 : fillCoordinates(indices, _points[i]);
272 472 : }
273 102 : }
274 :
275 : void
276 96 : NekSpatialBinUserObject::computePoints2D()
277 : {
278 : int p = 0;
279 1778 : for (unsigned int i = 0; i < _bins[0]->num_bins(); ++i)
280 : {
281 35796 : for (unsigned int j = 0; j < _bins[1]->num_bins(); ++j, ++p)
282 : {
283 34114 : std::vector<unsigned int> indices = {i, j};
284 34114 : fillCoordinates(indices, _points[p]);
285 34114 : }
286 : }
287 96 : }
288 :
289 : void
290 17 : NekSpatialBinUserObject::computePoints3D()
291 : {
292 : int p = 0;
293 64 : for (unsigned int i = 0; i < _bins[0]->num_bins(); ++i)
294 : {
295 180 : for (unsigned int j = 0; j < _bins[1]->num_bins(); ++j)
296 : {
297 10189 : for (unsigned int k = 0; k < _bins[2]->num_bins(); ++k, ++p)
298 : {
299 10056 : std::vector<unsigned int> indices = {i, j, k};
300 10056 : fillCoordinates(indices, _points[p]);
301 10056 : }
302 : }
303 : }
304 17 : }
305 :
306 : void
307 44642 : NekSpatialBinUserObject::fillCoordinates(const std::vector<unsigned int> & indices, Point & p) const
308 : {
309 143510 : for (unsigned int b = 0; b < _bins.size(); ++b)
310 : {
311 : const auto & bin = _bins[b];
312 98868 : const auto & centers = bin->getBinCenters();
313 :
314 232118 : for (const auto & d : bin->directions())
315 232118 : p(d) = centers[indices[b]](d);
316 : }
317 44642 : }
318 :
319 : #endif
|
