Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "LayeredBase.h"
11 :
12 : // MOOSE includes
13 : #include "MooseEnum.h"
14 : #include "MooseMesh.h"
15 : #include "SubProblem.h"
16 : #include "UserObject.h"
17 :
18 : #include "libmesh/mesh_tools.h"
19 : #include "libmesh/point.h"
20 :
21 : InputParameters
22 310077 : LayeredBase::validParams()
23 : {
24 310077 : InputParameters params = emptyInputParameters();
25 310077 : MooseEnum directions("x y z");
26 :
27 310077 : params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
28 310077 : params.addParam<unsigned int>("num_layers", "The number of layers.");
29 310077 : params.addParam<std::vector<Real>>("bounds",
30 : "The 'bounding' positions of the layers i.e.: '0, "
31 : "1.2, 3.7, 4.2' will mean 3 layers between those "
32 : "positions.");
33 310077 : params.addRangeCheckedParam<unsigned int>("bound_uniform_splits",
34 : "bound_uniform_splits > 0",
35 : "The number of times the bins specified in 'bounds' "
36 : "should be split uniformly.");
37 :
38 310077 : MooseEnum sample_options("direct interpolate average", "direct");
39 310077 : params.addParam<MooseEnum>("sample_type",
40 : sample_options,
41 : "How to sample the layers. 'direct' means get the value of the layer "
42 : "the point falls in directly (or average if that layer has no value). "
43 : " 'interpolate' does a linear interpolation between the two closest "
44 : "layers. 'average' averages the two closest layers.");
45 :
46 930231 : params.addParam<unsigned int>("average_radius",
47 620154 : 1,
48 : "When using 'average' sampling this is how "
49 : "the number of values both above and below "
50 : "the layer that will be averaged.");
51 :
52 930231 : params.addParam<bool>(
53 : "cumulative",
54 620154 : false,
55 : "When true the value in each layer is the sum of the values up to and including that layer");
56 930231 : params.addParam<bool>(
57 : "positive_cumulative_direction",
58 620154 : true,
59 : "When 'cumulative' is true, whether the direction for summing the cumulative value "
60 : "is the positive direction or negative direction");
61 :
62 310077 : params.addParam<std::vector<SubdomainName>>(
63 : "block", "The list of block ids (SubdomainID) that this object will be applied");
64 :
65 310077 : params.addParam<std::vector<SubdomainName>>("layer_bounding_block",
66 : "List of block ids (SubdomainID) that are used to "
67 : "determine the upper and lower geometric bounds for "
68 : "all layers. If this is not specified, the ids "
69 : "specified in 'block' are used for this purpose.");
70 :
71 310077 : params.addParam<Real>("direction_min",
72 : "Minimum coordinate along 'direction' that bounds the layers");
73 310077 : params.addParam<Real>("direction_max",
74 : "Maximum coordinate along 'direction' that bounds the layers");
75 310077 : params.addParamNamesToGroup(
76 : "direction num_layers bounds direction_min direction_max bound_uniform_splits",
77 : "Layers extent and definition");
78 310077 : params.addParamNamesToGroup("sample_type average_radius cumulative positive_cumulative_direction",
79 : "Value sampling / aggregating");
80 620154 : return params;
81 310077 : }
82 :
83 5126 : LayeredBase::LayeredBase(const InputParameters & parameters)
84 10252 : : Restartable(parameters.getCheckedPointerParam<SubProblem *>("_subproblem")->getMooseApp(),
85 10252 : parameters.get<std::string>("_object_name") + "_layered_base",
86 : "LayeredBase",
87 5126 : parameters.get<THREAD_ID>("_tid")),
88 5126 : _layered_base_name(parameters.get<std::string>("_object_name")),
89 5126 : _layered_base_params(parameters),
90 5126 : _direction_enum(parameters.get<MooseEnum>("direction")),
91 5126 : _direction(_direction_enum),
92 5126 : _sample_type(parameters.get<MooseEnum>("sample_type")),
93 5126 : _average_radius(parameters.get<unsigned int>("average_radius")),
94 5126 : _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
95 5126 : _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
96 5126 : _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
97 5126 : _cumulative(parameters.get<bool>("cumulative")),
98 5126 : _positive_cumulative_direction(parameters.get<bool>("positive_cumulative_direction")),
99 5126 : _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
100 5126 : _layer_bounding_blocks(),
101 30756 : _has_direction_max_min(false)
102 : {
103 15324 : if (_layered_base_params.isParamValid("num_layers") &&
104 10198 : _layered_base_params.isParamValid("bounds"))
105 4 : mooseError("'bounds' and 'num_layers' cannot both be set");
106 :
107 5122 : if (!_cumulative && parameters.isParamSetByUser("positive_cumulative_direction"))
108 0 : mooseWarning(
109 : "The 'positive_cumulative_direction' parameter is unused when 'cumulative' is false");
110 :
111 5122 : if (_layered_base_params.isParamValid("num_layers"))
112 : {
113 5068 : _num_layers = _layered_base_params.get<unsigned int>("num_layers");
114 5068 : _interval_based = true;
115 : }
116 54 : else if (_layered_base_params.isParamValid("bounds"))
117 : {
118 50 : _interval_based = false;
119 :
120 50 : _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
121 :
122 50 : if (_layer_bounds.size() < 2)
123 4 : mooseError("At least two boundaries must be provided in 'bounds' to form layers!");
124 :
125 : // Make sure the bounds are sorted - we're going to depend on this
126 46 : std::sort(_layer_bounds.begin(), _layer_bounds.end());
127 :
128 : // If requested, we uniformly split the layers.
129 46 : if (_layered_base_params.isParamValid("bound_uniform_splits"))
130 : {
131 14 : const auto splits = _layered_base_params.get<unsigned int>("bound_uniform_splits");
132 28 : for (unsigned int s = 0; s < splits; ++s)
133 : {
134 14 : std::vector<Real> new_bnds;
135 14 : new_bnds.reserve(2 * (_layer_bounds.size() - 1) + 1);
136 56 : for (unsigned int i = 0; i < _layer_bounds.size() - 1; ++i)
137 : {
138 42 : new_bnds.emplace_back(_layer_bounds[i]);
139 42 : new_bnds.emplace_back(0.5 * (_layer_bounds[i] + _layer_bounds[i + 1]));
140 : }
141 14 : _layer_bounds = new_bnds;
142 14 : }
143 : }
144 :
145 46 : _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
146 46 : _direction_min = _layer_bounds.front();
147 46 : _direction_max = _layer_bounds.back();
148 46 : _has_direction_max_min = true;
149 : }
150 : else
151 4 : mooseError("One of 'bounds' or 'num_layers' must be specified");
152 :
153 5114 : if (!_interval_based && _sample_type == 1)
154 4 : mooseError("'sample_type = interpolate' not supported with 'bounds'");
155 :
156 5110 : bool has_layer_bounding_block = _layered_base_params.isParamValid("layer_bounding_block");
157 5110 : bool has_block = _layered_base_params.isParamValid("block");
158 5110 : bool has_direction_min = _layered_base_params.isParamValid("direction_min");
159 5110 : bool has_direction_max = _layered_base_params.isParamValid("direction_max");
160 :
161 5110 : if (_has_direction_max_min && has_direction_min)
162 0 : mooseWarning("'direction_min' is unused when providing 'bounds'");
163 :
164 5110 : if (_has_direction_max_min && has_direction_max)
165 0 : mooseWarning("'direction_max' is unused when providing 'bounds'");
166 :
167 : // can only specify one of layer_bounding_block or the pair direction_max/min
168 5110 : if (has_layer_bounding_block && (has_direction_min || has_direction_max))
169 4 : mooseError("Only one of 'layer_bounding_block' and the pair 'direction_max' and "
170 : "'direction_min' can be provided");
171 :
172 : // if either one of direction_min or direction_max is specified, must provide the other one
173 5106 : if (has_direction_min != has_direction_max)
174 4 : mooseError("If providing the layer max/min directions, both 'direction_max' and "
175 : "'direction_min' must be specified.");
176 :
177 5102 : if (has_layer_bounding_block)
178 84 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
179 56 : _layered_base_params.get<std::vector<SubdomainName>>("layer_bounding_block"));
180 5074 : else if (has_block)
181 936 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
182 624 : _layered_base_params.get<std::vector<SubdomainName>>("block"));
183 :
184 : // specifying the direction max/min overrides anything set with the 'block'
185 5102 : if (has_direction_min && has_direction_max)
186 : {
187 60 : _direction_min = parameters.get<Real>("direction_min");
188 60 : _direction_max = parameters.get<Real>("direction_max");
189 60 : _has_direction_max_min = true;
190 :
191 60 : if (_direction_max <= _direction_min)
192 4 : mooseError("'direction_max' must be larger than 'direction_min'");
193 : }
194 :
195 5098 : _layer_values.resize(_num_layers);
196 5098 : _layer_has_value.resize(_num_layers);
197 :
198 : // if we haven't already figured out the max/min in specified direction
199 : // (either with the 'bounds' or explicit specification from the user), do so
200 5098 : if (!_has_direction_max_min)
201 5000 : getBounds();
202 :
203 5098 : computeLayerCenters();
204 5098 : }
205 :
206 : Real
207 3963281 : LayeredBase::integralValue(const Point & p) const
208 : {
209 3963281 : unsigned int layer = getLayer(p);
210 :
211 3963281 : int higher_layer = -1;
212 3963281 : int lower_layer = -1;
213 :
214 4100201 : for (unsigned int i = layer; i < _layer_values.size(); i++)
215 : {
216 4058791 : if (_layer_has_value[i])
217 : {
218 3921871 : higher_layer = i;
219 3921871 : break;
220 : }
221 : }
222 :
223 4150479 : for (int i = layer - 1; i >= 0; i--)
224 : {
225 3127396 : if (_layer_has_value[i])
226 : {
227 2940198 : lower_layer = i;
228 2940198 : break;
229 : }
230 : }
231 :
232 3963281 : if (higher_layer == -1 && lower_layer == -1)
233 41248 : return 0; // TODO: We could error here but there are startup dependency problems
234 :
235 3922033 : switch (_sample_type)
236 : {
237 3915553 : case 0: // direct
238 : {
239 3915553 : if (higher_layer == -1) // Didn't find a higher layer
240 162 : return _layer_values[lower_layer];
241 :
242 3915391 : if (unsigned(higher_layer) == layer) // constant in a layer
243 3890855 : return _layer_values[higher_layer];
244 :
245 24536 : if (lower_layer == -1) // Didn't find a lower layer
246 3243 : return _layer_values[higher_layer];
247 :
248 21293 : return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
249 : }
250 0 : case 1: // interpolate
251 : {
252 0 : if (higher_layer == -1) // Didn't find a higher layer
253 0 : return _layer_values[lower_layer];
254 :
255 0 : Real layer_length = (_direction_max - _direction_min) / _num_layers;
256 0 : Real lower_coor = _direction_min;
257 0 : Real lower_value = 0;
258 0 : if (lower_layer != -1)
259 : {
260 0 : lower_coor += (lower_layer + 1) * layer_length;
261 0 : lower_value = _layer_values[lower_layer];
262 : }
263 :
264 : // Interpolate between the two points
265 0 : Real higher_value = _layer_values[higher_layer];
266 :
267 : // Linear interpolation
268 : return lower_value +
269 0 : (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
270 : }
271 6480 : case 2: // average
272 : {
273 6480 : Real total = 0;
274 6480 : unsigned int num_values = 0;
275 :
276 6480 : if (higher_layer != -1)
277 : {
278 18144 : for (const auto i : make_range(_average_radius))
279 : {
280 12960 : int current_layer = higher_layer + i;
281 :
282 12960 : if ((size_t)current_layer >= _layer_values.size())
283 1296 : break;
284 :
285 11664 : if (_layer_has_value[current_layer])
286 : {
287 11664 : total += _layer_values[current_layer];
288 11664 : num_values += 1;
289 : }
290 : }
291 : }
292 :
293 6480 : if (lower_layer != -1)
294 : {
295 14256 : for (const auto i : make_range(_average_radius))
296 : {
297 10368 : int current_layer = lower_layer - i;
298 :
299 10368 : if (current_layer < 0)
300 1296 : break;
301 :
302 9072 : if (_layer_has_value[current_layer])
303 : {
304 9072 : total += _layer_values[current_layer];
305 9072 : num_values += 1;
306 : }
307 : }
308 : }
309 :
310 6480 : return total / num_values;
311 : }
312 0 : default:
313 0 : mooseError("Unknown sample type!");
314 : }
315 : }
316 :
317 : Real
318 822348 : LayeredBase::getLayerValue(unsigned int layer) const
319 : {
320 822348 : if (layer >= _layer_values.size())
321 0 : mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
322 822348 : return _layer_values[layer];
323 : }
324 :
325 : void
326 10693 : LayeredBase::initialize()
327 : {
328 10693 : if (_using_displaced_mesh)
329 92 : getBounds();
330 :
331 74457 : for (const auto i : index_range(_layer_values))
332 : {
333 63764 : _layer_values[i] = 0.0;
334 63764 : _layer_has_value[i] = false;
335 : }
336 10693 : }
337 :
338 : void
339 9854 : LayeredBase::finalize()
340 : {
341 9854 : _layered_base_subproblem.comm().sum(_layer_values);
342 9854 : _layered_base_subproblem.comm().max(_layer_has_value);
343 :
344 9854 : if (_cumulative)
345 : {
346 456 : Real value = 0;
347 :
348 456 : if (_positive_cumulative_direction)
349 912 : for (const auto i : make_range(_num_layers))
350 : {
351 684 : value += getLayerValue(i);
352 684 : setLayerValue(i, value);
353 : }
354 : else
355 912 : for (int i = _num_layers - 1; i >= 0; i--)
356 : {
357 684 : value += getLayerValue(i);
358 684 : setLayerValue(i, value);
359 : }
360 : }
361 9854 : }
362 :
363 : void
364 813 : LayeredBase::threadJoin(const UserObject & y)
365 : {
366 813 : const auto & lb = dynamic_cast<const LayeredBase &>(y);
367 5925 : for (const auto i : index_range(_layer_values))
368 5112 : if (lb.layerHasValue(i))
369 2494 : setLayerValue(i, getLayerValue(i) + lb._layer_values[i]);
370 813 : }
371 :
372 : unsigned int
373 5433569 : LayeredBase::getLayer(const Point & p) const
374 : {
375 5433569 : Real direction_x = p(_direction);
376 :
377 5433569 : if (direction_x < _direction_min)
378 19341 : return 0;
379 :
380 5414228 : if (_interval_based)
381 : {
382 5403188 : unsigned int layer =
383 5403188 : std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
384 5403188 : static_cast<Real>(_num_layers));
385 :
386 5403188 : if (layer >= _num_layers)
387 84566 : layer = _num_layers - 1;
388 :
389 5403188 : return layer;
390 : }
391 : else // Figure out what layer we are in from the bounds
392 : {
393 : // This finds the first entry in the vector that is larger than what we're looking for
394 : std::vector<Real>::const_iterator one_higher =
395 11040 : std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
396 :
397 11040 : if (one_higher == _layer_bounds.end())
398 : {
399 : return static_cast<unsigned int>(
400 1620 : _layer_bounds.size() -
401 1620 : 2); // Just return the last layer. -2 because layers are "in-between" bounds
402 : }
403 9420 : else if (one_higher == _layer_bounds.begin())
404 0 : return 0; // Return the first layer
405 : else
406 : // The -1 is because the interval that we fall in is just _before_ the number that is bigger
407 : // (which is what we found
408 9420 : return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
409 : }
410 : }
411 :
412 : void
413 5098 : LayeredBase::computeLayerCenters()
414 : {
415 5098 : _layer_centers.resize(_num_layers);
416 :
417 5098 : if (_interval_based)
418 : {
419 5056 : Real dx = (_direction_max - _direction_min) / _num_layers;
420 :
421 44929 : for (const auto i : make_range(_num_layers))
422 39873 : _layer_centers[i] = (i + 0.5) * dx;
423 : }
424 : else
425 : {
426 182 : for (const auto i : make_range(_num_layers))
427 140 : _layer_centers[i] = 0.5 * (_layer_bounds[i + 1] + _layer_bounds[i]);
428 : }
429 5098 : }
430 :
431 : void
432 822588 : LayeredBase::setLayerValue(unsigned int layer, Real value)
433 : {
434 822588 : _layer_values[layer] = value;
435 822588 : _layer_has_value[layer] = true;
436 822588 : }
437 :
438 : void
439 5092 : LayeredBase::getBounds()
440 : {
441 5092 : if (_layer_bounding_blocks.size() == 0)
442 : {
443 4766 : BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
444 4766 : _direction_min = bounding_box.min()(_direction);
445 4766 : _direction_max = bounding_box.max()(_direction);
446 : }
447 : else
448 : {
449 326 : _direction_min = std::numeric_limits<Real>::infinity();
450 326 : _direction_max = -std::numeric_limits<Real>::infinity();
451 :
452 326 : MooseMesh & mesh = _layered_base_subproblem.mesh();
453 :
454 4788 : for (auto & elem_ptr : *mesh.getActiveLocalElementRange())
455 : {
456 4462 : auto subdomain_id = elem_ptr->subdomain_id();
457 :
458 4462 : if (std::find(_layer_bounding_blocks.begin(), _layer_bounding_blocks.end(), subdomain_id) ==
459 8924 : _layer_bounding_blocks.end())
460 2753 : continue;
461 :
462 7175 : for (auto & node : elem_ptr->node_ref_range())
463 : {
464 5466 : _direction_min = std::min(_direction_min, node(_direction));
465 5466 : _direction_max = std::max(_direction_max, node(_direction));
466 : }
467 : }
468 :
469 326 : mesh.comm().min(_direction_min);
470 326 : mesh.comm().max(_direction_max);
471 : }
472 5092 : }
|
