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 75620 : LayeredBase::validParams()
23 : {
24 75620 : InputParameters params = emptyInputParameters();
25 302480 : MooseEnum directions("x y z");
26 :
27 302480 : params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
28 302480 : params.addParam<unsigned int>("num_layers", "The number of layers.");
29 302480 : 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 453720 : 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 302480 : params.addParam<std::vector<unsigned int>>(
38 : "bound_splits", "Number of uniform splits of all layers in 'bounds' (default to all ones)");
39 :
40 302480 : MooseEnum sample_options("direct interpolate average", "direct");
41 302480 : params.addParam<MooseEnum>("sample_type",
42 : sample_options,
43 : "How to sample the layers. 'direct' means get the value of the layer "
44 : "the point falls in directly (or average if that layer has no value). "
45 : " 'interpolate' does a linear interpolation between the two closest "
46 : "layers. 'average' averages the two closest layers.");
47 :
48 226860 : params.addParam<unsigned int>("average_radius",
49 151240 : 1,
50 : "When using 'average' sampling this is how "
51 : "the number of values both above and below "
52 : "the layer that will be averaged.");
53 :
54 226860 : params.addParam<bool>(
55 : "cumulative",
56 151240 : false,
57 : "When true the value in each layer is the sum of the values up to and including that layer");
58 226860 : params.addParam<bool>(
59 : "positive_cumulative_direction",
60 151240 : true,
61 : "When 'cumulative' is true, whether the direction for summing the cumulative value "
62 : "is the positive direction or negative direction");
63 :
64 302480 : params.addParam<std::vector<SubdomainName>>(
65 : "block", "The list of block ids (SubdomainID) that this object will be applied");
66 :
67 302480 : params.addParam<std::vector<SubdomainName>>("layer_bounding_block",
68 : "List of block ids (SubdomainID) that are used to "
69 : "determine the upper and lower geometric bounds for "
70 : "all layers. If this is not specified, the ids "
71 : "specified in 'block' are used for this purpose.");
72 :
73 302480 : params.addParam<Real>("direction_min",
74 : "Minimum coordinate along 'direction' that bounds the layers");
75 302480 : params.addParam<Real>("direction_max",
76 : "Maximum coordinate along 'direction' that bounds the layers");
77 302480 : params.addParamNamesToGroup(
78 : "direction num_layers bounds direction_min direction_max bound_uniform_splits bound_splits",
79 : "Layers extent and definition");
80 226860 : params.addParamNamesToGroup("sample_type average_radius cumulative positive_cumulative_direction",
81 : "Value sampling / aggregating");
82 151240 : return params;
83 75620 : }
84 :
85 5614 : LayeredBase::LayeredBase(const InputParameters & parameters)
86 22456 : : Restartable(parameters.getCheckedPointerParam<SubProblem *>("_subproblem")->getMooseApp(),
87 11228 : parameters.getObjectName() + "_layered_base",
88 : "LayeredBase",
89 5614 : parameters.get<THREAD_ID>("_tid")),
90 5614 : _layered_base_name(parameters.getObjectName()),
91 5614 : _layered_base_params(parameters),
92 5614 : _direction_enum(parameters.get<MooseEnum>("direction")),
93 5614 : _direction(_direction_enum),
94 5614 : _sample_type(parameters.get<MooseEnum>("sample_type")),
95 5614 : _average_radius(parameters.get<unsigned int>("average_radius")),
96 5614 : _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
97 11228 : _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
98 11228 : _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
99 5614 : _cumulative(parameters.get<bool>("cumulative")),
100 5614 : _positive_cumulative_direction(parameters.get<bool>("positive_cumulative_direction")),
101 16842 : _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
102 5614 : _layer_bounding_blocks(),
103 56140 : _has_direction_max_min(false)
104 : {
105 28009 : if (_layered_base_params.isParamValid("num_layers") &&
106 22273 : _layered_base_params.isParamValid("bounds"))
107 3 : mooseError("'bounds' and 'num_layers' cannot both be set");
108 :
109 16781 : if (!_cumulative && parameters.isParamSetByUser("positive_cumulative_direction"))
110 0 : mooseWarning(
111 : "The 'positive_cumulative_direction' parameter is unused when 'cumulative' is false");
112 :
113 16833 : if (_layered_base_params.isParamValid("num_layers"))
114 : {
115 5550 : _num_layers = _layered_base_params.get<unsigned int>("num_layers");
116 5550 : _interval_based = true;
117 : }
118 183 : else if (_layered_base_params.isParamValid("bounds"))
119 : {
120 58 : _interval_based = false;
121 :
122 58 : _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
123 :
124 58 : if (_layer_bounds.size() < 2)
125 3 : mooseError("At least two boundaries must be provided in 'bounds' to form layers!");
126 :
127 : // Make sure the bounds are sorted - we're going to depend on this
128 55 : std::sort(_layer_bounds.begin(), _layer_bounds.end());
129 :
130 233 : if (_layered_base_params.isParamValid("bound_splits") &&
131 94 : _layered_base_params.isParamValid("bound_uniform_splits"))
132 0 : mooseError("Parameters 'bound_splits' and 'bound_uniform_splits' cannot be both supplied");
133 :
134 : // If requested, we uniformly split the layers.
135 165 : if (_layered_base_params.isParamValid("bound_uniform_splits"))
136 : {
137 13 : const auto splits = _layered_base_params.get<unsigned int>("bound_uniform_splits");
138 26 : for (unsigned int s = 0; s < splits; ++s)
139 : {
140 13 : std::vector<Real> new_bnds;
141 13 : new_bnds.reserve(2 * (_layer_bounds.size() - 1) + 1);
142 52 : for (unsigned int i = 0; i < _layer_bounds.size() - 1; ++i)
143 : {
144 39 : new_bnds.emplace_back(_layer_bounds[i]);
145 39 : new_bnds.emplace_back(0.5 * (_layer_bounds[i] + _layer_bounds[i + 1]));
146 : }
147 13 : new_bnds.emplace_back(_layer_bounds.back());
148 13 : _layer_bounds = new_bnds;
149 13 : }
150 : }
151 :
152 165 : if (_layered_base_params.isParamValid("bound_splits"))
153 : {
154 : // expand layers with subdivisions
155 13 : const auto nsubs = _layered_base_params.get<std::vector<unsigned int>>("bound_splits");
156 13 : if (nsubs.size() != _layer_bounds.size() - 1)
157 0 : mooseError("'bound_splits' size must be equal to the size of 'bounds' minus one");
158 :
159 13 : std::vector<Real> new_layers;
160 13 : new_layers.reserve(std::accumulate(nsubs.begin(), nsubs.end(), (unsigned int)0));
161 13 : new_layers.push_back(_layer_bounds[0]);
162 26 : for (const auto i : index_range(nsubs))
163 : {
164 13 : Real dz = (_layer_bounds[i + 1] - _layer_bounds[i]) / nsubs[i];
165 52 : for (const auto j : make_range(nsubs[i]))
166 : {
167 39 : libmesh_ignore(j);
168 39 : new_layers.push_back(new_layers.back() + dz);
169 : }
170 : }
171 13 : _layer_bounds = new_layers;
172 13 : }
173 :
174 55 : _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
175 55 : _direction_min = _layer_bounds.front();
176 55 : _direction_max = _layer_bounds.back();
177 55 : _has_direction_max_min = true;
178 : }
179 : else
180 3 : mooseError("One of 'bounds' or 'num_layers' must be specified");
181 :
182 5605 : if (!_interval_based && _sample_type == 1)
183 3 : mooseError("'sample_type = interpolate' not supported with 'bounds'");
184 :
185 11204 : bool has_layer_bounding_block = _layered_base_params.isParamValid("layer_bounding_block");
186 11204 : bool has_block = _layered_base_params.isParamValid("block");
187 11204 : bool has_direction_min = _layered_base_params.isParamValid("direction_min");
188 11204 : bool has_direction_max = _layered_base_params.isParamValid("direction_max");
189 :
190 5602 : if (_has_direction_max_min && has_direction_min)
191 0 : mooseWarning("'direction_min' is unused when providing 'bounds'");
192 :
193 5602 : if (_has_direction_max_min && has_direction_max)
194 0 : mooseWarning("'direction_max' is unused when providing 'bounds'");
195 :
196 : // can only specify one of layer_bounding_block or the pair direction_max/min
197 5602 : if (has_layer_bounding_block && (has_direction_min || has_direction_max))
198 3 : mooseError("Only one of 'layer_bounding_block' and the pair 'direction_max' and "
199 : "'direction_min' can be provided");
200 :
201 : // if either one of direction_min or direction_max is specified, must provide the other one
202 5599 : if (has_direction_min != has_direction_max)
203 3 : mooseError("If providing the layer max/min directions, both 'direction_max' and "
204 : "'direction_min' must be specified.");
205 :
206 5596 : if (has_layer_bounding_block)
207 78 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
208 52 : _layered_base_params.get<std::vector<SubdomainName>>("layer_bounding_block"));
209 5570 : else if (has_block)
210 858 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
211 572 : _layered_base_params.get<std::vector<SubdomainName>>("block"));
212 :
213 : // specifying the direction max/min overrides anything set with the 'block'
214 5596 : if (has_direction_min && has_direction_max)
215 : {
216 55 : _direction_min = parameters.get<Real>("direction_min");
217 55 : _direction_max = parameters.get<Real>("direction_max");
218 55 : _has_direction_max_min = true;
219 :
220 55 : if (_direction_max <= _direction_min)
221 3 : mooseError("'direction_max' must be larger than 'direction_min'");
222 : }
223 :
224 5593 : _layer_values.resize(_num_layers);
225 5593 : _layer_has_value.resize(_num_layers);
226 :
227 : // if we haven't already figured out the max/min in specified direction
228 : // (either with the 'bounds' or explicit specification from the user), do so
229 5593 : if (!_has_direction_max_min)
230 5489 : getBounds();
231 :
232 5593 : computeLayerCenters();
233 5593 : }
234 :
235 : Real
236 3652048 : LayeredBase::integralValue(const Point & p) const
237 : {
238 3652048 : unsigned int layer = getLayer(p);
239 :
240 3652048 : int higher_layer = -1;
241 3652048 : int lower_layer = -1;
242 :
243 3771350 : for (unsigned int i = layer; i < _layer_values.size(); i++)
244 : {
245 3736940 : if (_layer_has_value[i])
246 : {
247 3617638 : higher_layer = i;
248 3617638 : break;
249 : }
250 : }
251 :
252 3847213 : for (int i = layer - 1; i >= 0; i--)
253 : {
254 2916399 : if (_layer_has_value[i])
255 : {
256 2721234 : lower_layer = i;
257 2721234 : break;
258 : }
259 : }
260 :
261 3652048 : if (higher_layer == -1 && lower_layer == -1)
262 34410 : return 0; // TODO: We could error here but there are startup dependency problems
263 :
264 3617638 : switch (_sample_type)
265 : {
266 3611878 : case 0: // direct
267 : {
268 3611878 : if (higher_layer == -1) // Didn't find a higher layer
269 0 : return _layer_values[lower_layer];
270 :
271 3611878 : if (unsigned(higher_layer) == layer) // constant in a layer
272 3578107 : return _layer_values[higher_layer];
273 :
274 33771 : if (lower_layer == -1) // Didn't find a lower layer
275 2486 : return _layer_values[higher_layer];
276 :
277 31285 : return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
278 : }
279 0 : case 1: // interpolate
280 : {
281 0 : if (higher_layer == -1) // Didn't find a higher layer
282 0 : return _layer_values[lower_layer];
283 :
284 0 : Real layer_length = (_direction_max - _direction_min) / _num_layers;
285 0 : Real lower_coor = _direction_min;
286 0 : Real lower_value = 0;
287 0 : if (lower_layer != -1)
288 : {
289 0 : lower_coor += (lower_layer + 1) * layer_length;
290 0 : lower_value = _layer_values[lower_layer];
291 : }
292 :
293 : // Interpolate between the two points
294 0 : Real higher_value = _layer_values[higher_layer];
295 :
296 : // Linear interpolation
297 : return lower_value +
298 0 : (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
299 : }
300 5760 : case 2: // average
301 : {
302 5760 : Real total = 0;
303 5760 : unsigned int num_values = 0;
304 :
305 5760 : if (higher_layer != -1)
306 : {
307 16128 : for (const auto i : make_range(_average_radius))
308 : {
309 11520 : int current_layer = higher_layer + i;
310 :
311 11520 : if ((size_t)current_layer >= _layer_values.size())
312 1152 : break;
313 :
314 10368 : if (_layer_has_value[current_layer])
315 : {
316 10368 : total += _layer_values[current_layer];
317 10368 : num_values += 1;
318 : }
319 : }
320 : }
321 :
322 5760 : if (lower_layer != -1)
323 : {
324 12672 : for (const auto i : make_range(_average_radius))
325 : {
326 9216 : int current_layer = lower_layer - i;
327 :
328 9216 : if (current_layer < 0)
329 1152 : break;
330 :
331 8064 : if (_layer_has_value[current_layer])
332 : {
333 8064 : total += _layer_values[current_layer];
334 8064 : num_values += 1;
335 : }
336 : }
337 : }
338 :
339 5760 : return total / num_values;
340 : }
341 0 : default:
342 0 : mooseError("Unknown sample type!");
343 : }
344 : }
345 :
346 : Real
347 807343 : LayeredBase::getLayerValue(unsigned int layer) const
348 : {
349 807343 : if (layer >= _layer_values.size())
350 0 : mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
351 807343 : return _layer_values[layer];
352 : }
353 :
354 : void
355 10976 : LayeredBase::initialize()
356 : {
357 10976 : if (_using_displaced_mesh)
358 83 : getBounds();
359 :
360 79030 : for (const auto i : index_range(_layer_values))
361 : {
362 68054 : _layer_values[i] = 0.0;
363 68054 : _layer_has_value[i] = false;
364 : }
365 10976 : }
366 :
367 : void
368 10030 : LayeredBase::finalize()
369 : {
370 10030 : _layered_base_subproblem.comm().sum(_layer_values);
371 10030 : _layered_base_subproblem.comm().max(_layer_has_value);
372 :
373 10030 : if (_cumulative)
374 : {
375 422 : Real value = 0;
376 :
377 422 : if (_positive_cumulative_direction)
378 844 : for (const auto i : make_range(_num_layers))
379 : {
380 633 : value += getLayerValue(i);
381 633 : setLayerValue(i, value);
382 : }
383 : else
384 844 : for (int i = _num_layers - 1; i >= 0; i--)
385 : {
386 633 : value += getLayerValue(i);
387 633 : setLayerValue(i, value);
388 : }
389 : }
390 10030 : }
391 :
392 : void
393 922 : LayeredBase::threadJoin(const UserObject & y)
394 : {
395 922 : const auto & lb = dynamic_cast<const LayeredBase &>(y);
396 7006 : for (const auto i : index_range(_layer_values))
397 6084 : if (lb.layerHasValue(i))
398 3089 : setLayerValue(i, getLayerValue(i) + lb._layer_values[i]);
399 922 : }
400 :
401 : unsigned int
402 5092612 : LayeredBase::getLayer(const Point & p) const
403 : {
404 5092612 : Real direction_x = p(_direction);
405 :
406 5092612 : if (direction_x < _direction_min)
407 17021 : return 0;
408 :
409 5075591 : if (_interval_based)
410 : {
411 5059262 : unsigned int layer =
412 5059262 : std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
413 5059262 : static_cast<Real>(_num_layers));
414 :
415 5059262 : if (layer >= _num_layers)
416 75763 : layer = _num_layers - 1;
417 :
418 5059262 : return layer;
419 : }
420 : else // Figure out what layer we are in from the bounds
421 : {
422 : // This finds the first entry in the vector that is larger than what we're looking for
423 : std::vector<Real>::const_iterator one_higher =
424 16329 : std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
425 :
426 16329 : if (one_higher == _layer_bounds.end())
427 : {
428 : return static_cast<unsigned int>(
429 0 : _layer_bounds.size() -
430 0 : 2); // Just return the last layer. -2 because layers are "in-between" bounds
431 : }
432 16329 : else if (one_higher == _layer_bounds.begin())
433 0 : return 0; // Return the first layer
434 : else
435 : // The -1 is because the interval that we fall in is just _before_ the number that is bigger
436 : // (which is what we found
437 32658 : return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
438 : }
439 : }
440 :
441 : void
442 5593 : LayeredBase::computeLayerCenters()
443 : {
444 5593 : _layer_centers.resize(_num_layers);
445 :
446 5593 : if (_interval_based)
447 : {
448 5541 : Real dx = (_direction_max - _direction_min) / _num_layers;
449 :
450 51143 : for (const auto i : make_range(_num_layers))
451 45602 : _layer_centers[i] = (i + 0.5) * dx;
452 : }
453 : else
454 : {
455 234 : for (const auto i : make_range(_num_layers))
456 182 : _layer_centers[i] = 0.5 * (_layer_bounds[i + 1] + _layer_bounds[i]);
457 : }
458 5593 : }
459 :
460 : void
461 807563 : LayeredBase::setLayerValue(unsigned int layer, Real value)
462 : {
463 807563 : _layer_values[layer] = value;
464 807563 : _layer_has_value[layer] = true;
465 807563 : }
466 :
467 : void
468 5572 : LayeredBase::getBounds()
469 : {
470 5572 : if (_layer_bounding_blocks.size() == 0)
471 : {
472 5273 : BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
473 5273 : _direction_min = bounding_box.min()(_direction);
474 5273 : _direction_max = bounding_box.max()(_direction);
475 : }
476 : else
477 : {
478 299 : _direction_min = std::numeric_limits<Real>::infinity();
479 299 : _direction_max = -std::numeric_limits<Real>::infinity();
480 :
481 299 : MooseMesh & mesh = _layered_base_subproblem.mesh();
482 :
483 4301 : for (auto & elem_ptr : *mesh.getActiveLocalElementRange())
484 : {
485 4002 : auto subdomain_id = elem_ptr->subdomain_id();
486 :
487 4002 : if (std::find(_layer_bounding_blocks.begin(), _layer_bounding_blocks.end(), subdomain_id) ==
488 8004 : _layer_bounding_blocks.end())
489 2487 : continue;
490 :
491 6361 : for (auto & node : elem_ptr->node_ref_range())
492 : {
493 4846 : _direction_min = std::min(_direction_min, node(_direction));
494 4846 : _direction_max = std::max(_direction_max, node(_direction));
495 : }
496 : }
497 :
498 299 : mesh.comm().min(_direction_min);
499 299 : mesh.comm().max(_direction_max);
500 : }
501 5572 : }
|
