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 489 : NekSpatialBinUserObject::validParams()
27 : {
28 489 : InputParameters params = GeneralUserObject::validParams();
29 489 : params += NekBase::validParams();
30 489 : params += NekFieldInterface::validParams();
31 978 : params.addRequiredParam<std::vector<UserObjectName>>(
32 : "bins", "Userobjects providing a spatial bin given a point");
33 978 : params.addParam<unsigned int>(
34 : "interval",
35 978 : 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 978 : params.addParam<bool>(
40 : "map_space_by_qp",
41 978 : 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 978 : params.addParam<bool>(
45 : "check_zero_contributions",
46 978 : 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 489 : return params;
52 0 : }
53 :
54 248 : NekSpatialBinUserObject::NekSpatialBinUserObject(const InputParameters & parameters)
55 : : GeneralUserObject(parameters),
56 : NekBase(this, parameters),
57 : NekFieldInterface(this, parameters, true /* allow normal */),
58 247 : _interval(getParam<unsigned int>("interval")),
59 494 : _bin_names(getParam<std::vector<UserObjectName>>("bins")),
60 494 : _map_space_by_qp(getParam<bool>("map_space_by_qp")),
61 494 : _check_zero_contributions(getParam<bool>("check_zero_contributions")),
62 247 : _bin_values(nullptr),
63 247 : _bin_values_x(nullptr),
64 247 : _bin_values_y(nullptr),
65 247 : _bin_values_z(nullptr),
66 247 : _bin_volumes(nullptr),
67 247 : _bin_counts(nullptr),
68 247 : _bin_partial_values(nullptr),
69 495 : _bin_partial_counts(nullptr)
70 : {
71 247 : _fixed_mesh = !nekrs::hasMovingMesh();
72 :
73 247 : if (_bin_names.size() == 0)
74 0 : paramError("bins", "Length of vector must be greater than zero!");
75 :
76 652 : for (auto & b : _bin_names)
77 : {
78 : // first check that the user object exists
79 407 : 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 406 : 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 405 : _bins.push_back(&getUserObjectByName<SpatialBinUserObject>(b));
93 : }
94 :
95 245 : _n_bins = num_bins();
96 :
97 245 : _bin_values = (double *)calloc(_n_bins, sizeof(double));
98 245 : _bin_volumes = (double *)calloc(_n_bins, sizeof(double));
99 245 : _bin_partial_values = (double *)calloc(_n_bins, sizeof(double));
100 245 : _bin_counts = (int *)calloc(_n_bins, sizeof(int));
101 245 : _bin_partial_counts = (int *)calloc(_n_bins, sizeof(int));
102 :
103 245 : if (_field == field::velocity_component)
104 : {
105 25 : _bin_values_x = (double *)calloc(_n_bins, sizeof(double));
106 25 : _bin_values_y = (double *)calloc(_n_bins, sizeof(double));
107 25 : _bin_values_z = (double *)calloc(_n_bins, sizeof(double));
108 : }
109 :
110 245 : _has_direction = {false, false, false};
111 245 : _bin_providing_direction.resize(3);
112 648 : for (unsigned int b = 0; b < _bins.size(); ++b)
113 : {
114 : const auto & bin = _bins[b];
115 :
116 : // directions provided by this bin
117 405 : auto bin_directions = bin->directions();
118 :
119 948 : for (const auto & d : bin_directions)
120 : {
121 545 : 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 543 : _bin_providing_direction[d] = b;
132 : }
133 : }
134 : }
135 :
136 : // initialize all points to (0, 0, 0)
137 50361 : for (unsigned int i = 0; i < _n_bins; ++i)
138 50118 : _points.push_back(Point(0.0, 0.0, 0.0));
139 :
140 : // we will at most have 3 separate distributions
141 243 : if (_bins.size() == 1)
142 104 : computePoints1D();
143 139 : else if (_bins.size() == 2)
144 122 : computePoints2D();
145 : else
146 17 : computePoints3D();
147 :
148 : // with a user-specified direction, the direction for each bin is the same
149 243 : 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 243 : }
153 :
154 236 : NekSpatialBinUserObject::~NekSpatialBinUserObject()
155 : {
156 236 : freePointer(_bin_values);
157 236 : freePointer(_bin_volumes);
158 236 : freePointer(_bin_counts);
159 236 : freePointer(_bin_partial_values);
160 236 : freePointer(_bin_partial_counts);
161 :
162 236 : freePointer(_bin_values_x);
163 236 : freePointer(_bin_values_y);
164 236 : freePointer(_bin_values_z);
165 236 : }
166 :
167 : void
168 1964 : NekSpatialBinUserObject::execute()
169 : {
170 1964 : if (_fe_problem.timeStep() % _interval == 0)
171 348 : executeUserObject();
172 1964 : }
173 :
174 : Point
175 245866656 : NekSpatialBinUserObject::nekPoint(const int & local_elem_id, const int & local_node_id) const
176 : {
177 245866656 : if (_map_space_by_qp)
178 218869152 : return nekrs::gllPoint(local_elem_id, local_node_id);
179 : else
180 26997504 : return nekrs::centroid(local_elem_id);
181 : }
182 :
183 : void
184 635 : NekSpatialBinUserObject::resetPartialStorage()
185 : {
186 90894 : for (unsigned int i = 0; i < _n_bins; ++i)
187 : {
188 90259 : _bin_partial_values[i] = 0.0;
189 90259 : _bin_partial_counts[i] = 0;
190 : }
191 635 : }
192 :
193 : void
194 239 : NekSpatialBinUserObject::computeBinVolumes()
195 : {
196 239 : getBinVolumes();
197 :
198 239 : if (_check_zero_contributions)
199 : {
200 17023 : for (unsigned int i = 0; i < _n_bins; ++i)
201 : {
202 16798 : 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 237 : }
215 :
216 : Real
217 2136608 : NekSpatialBinUserObject::spatialValue(const Point & p) const
218 : {
219 2136608 : return _bin_values[bin(p)];
220 : }
221 :
222 : const std::vector<unsigned int>
223 127448 : NekSpatialBinUserObject::unrolledBin(const unsigned int & total_bin_index) const
224 : {
225 : std::vector<unsigned int> local_bins;
226 127448 : local_bins.resize(_bins.size());
227 :
228 127448 : int running_index = total_bin_index;
229 360880 : for (int i = _bins.size() - 1; i >= 0; --i)
230 : {
231 233432 : local_bins[i] = running_index % _bins[i]->num_bins();
232 233432 : running_index -= local_bins[i];
233 233432 : running_index /= _bins[i]->num_bins();
234 : }
235 :
236 127448 : return local_bins;
237 : }
238 :
239 : const unsigned int
240 159219476 : 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 465941188 : for (const auto & b : _bins)
245 306721712 : indices.push_back(b->bin(p));
246 :
247 : // convert to a total index into the multidimensional bin union
248 159219476 : unsigned int index = indices[0];
249 306721712 : for (unsigned int i = 1; i < _bins.size(); ++i)
250 147502236 : index = index * _bins[i]->num_bins() + indices[i];
251 :
252 159219476 : return index;
253 : }
254 :
255 : const unsigned int
256 32981 : NekSpatialBinUserObject::num_bins() const
257 : {
258 : unsigned int num_bins = 1;
259 98638 : for (const auto & b : _bins)
260 65657 : num_bins *= b->num_bins();
261 :
262 32981 : return num_bins;
263 : }
264 :
265 : void
266 104 : NekSpatialBinUserObject::computePoints1D()
267 : {
268 586 : for (unsigned int i = 0; i < _bins[0]->num_bins(); ++i)
269 : {
270 482 : std::vector<unsigned int> indices = {i};
271 482 : fillCoordinates(indices, _points[i]);
272 : }
273 104 : }
274 :
275 : void
276 122 : NekSpatialBinUserObject::computePoints2D()
277 : {
278 : int p = 0;
279 2398 : for (unsigned int i = 0; i < _bins[0]->num_bins(); ++i)
280 : {
281 41856 : for (unsigned int j = 0; j < _bins[1]->num_bins(); ++j, ++p)
282 : {
283 39580 : std::vector<unsigned int> indices = {i, j};
284 39580 : fillCoordinates(indices, _points[p]);
285 : }
286 : }
287 122 : }
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 : }
302 : }
303 : }
304 17 : }
305 :
306 : void
307 50118 : NekSpatialBinUserObject::fillCoordinates(const std::vector<unsigned int> & indices, Point & p) const
308 : {
309 159928 : for (unsigned int b = 0; b < _bins.size(); ++b)
310 : {
311 : const auto & bin = _bins[b];
312 109810 : const auto & centers = bin->getBinCenters();
313 :
314 259438 : for (const auto & d : bin->directions())
315 149628 : p(d) = centers[indices[b]](d);
316 : }
317 50118 : }
318 :
319 : #endif
|
