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