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 309213 : LayeredBase::validParams()
23 : {
24 309213 : InputParameters params = emptyInputParameters();
25 309213 : MooseEnum directions("x y z");
26 :
27 309213 : params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
28 309213 : params.addParam<unsigned int>("num_layers", "The number of layers.");
29 309213 : 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 :
34 309213 : MooseEnum sample_options("direct interpolate average", "direct");
35 309213 : params.addParam<MooseEnum>("sample_type",
36 : sample_options,
37 : "How to sample the layers. 'direct' means get the value of the layer "
38 : "the point falls in directly (or average if that layer has no value). "
39 : " 'interpolate' does a linear interpolation between the two closest "
40 : "layers. 'average' averages the two closest layers.");
41 :
42 927639 : params.addParam<unsigned int>("average_radius",
43 618426 : 1,
44 : "When using 'average' sampling this is how "
45 : "the number of values both above and below "
46 : "the layer that will be averaged.");
47 :
48 927639 : params.addParam<bool>(
49 : "cumulative",
50 618426 : false,
51 : "When true the value in each layer is the sum of the values up to and including that layer");
52 927639 : params.addParam<bool>(
53 : "positive_cumulative_direction",
54 618426 : true,
55 : "When 'cumulative' is true, whether the direction for summing the cumulative value "
56 : "is the positive direction or negative direction");
57 :
58 309213 : params.addParam<std::vector<SubdomainName>>(
59 : "block", "The list of block ids (SubdomainID) that this object will be applied");
60 :
61 309213 : params.addParam<std::vector<SubdomainName>>("layer_bounding_block",
62 : "List of block ids (SubdomainID) that are used to "
63 : "determine the upper and lower geometric bounds for "
64 : "all layers. If this is not specified, the ids "
65 : "specified in 'block' are used for this purpose.");
66 :
67 309213 : params.addParam<Real>("direction_min",
68 : "Minimum coordinate along 'direction' that bounds the layers");
69 309213 : params.addParam<Real>("direction_max",
70 : "Maximum coordinate along 'direction' that bounds the layers");
71 309213 : params.addParamNamesToGroup("direction num_layers bounds direction_min direction_max",
72 : "Layers extent and definition");
73 309213 : params.addParamNamesToGroup("sample_type average_radius cumulative positive_cumulative_direction",
74 : "Value sampling / aggregating");
75 618426 : return params;
76 309213 : }
77 :
78 4722 : LayeredBase::LayeredBase(const InputParameters & parameters)
79 9444 : : Restartable(parameters.getCheckedPointerParam<SubProblem *>("_subproblem")->getMooseApp(),
80 9444 : parameters.get<std::string>("_object_name") + "_layered_base",
81 : "LayeredBase",
82 4722 : parameters.get<THREAD_ID>("_tid")),
83 4722 : _layered_base_name(parameters.get<std::string>("_object_name")),
84 4722 : _layered_base_params(parameters),
85 4722 : _direction_enum(parameters.get<MooseEnum>("direction")),
86 4722 : _direction(_direction_enum),
87 4722 : _sample_type(parameters.get<MooseEnum>("sample_type")),
88 4722 : _average_radius(parameters.get<unsigned int>("average_radius")),
89 4722 : _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
90 4722 : _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
91 4722 : _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
92 4722 : _cumulative(parameters.get<bool>("cumulative")),
93 4722 : _positive_cumulative_direction(parameters.get<bool>("positive_cumulative_direction")),
94 4722 : _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
95 4722 : _layer_bounding_blocks(),
96 28332 : _has_direction_max_min(false)
97 : {
98 14132 : if (_layered_base_params.isParamValid("num_layers") &&
99 9410 : _layered_base_params.isParamValid("bounds"))
100 4 : mooseError("'bounds' and 'num_layers' cannot both be set");
101 :
102 4718 : if (!_cumulative && parameters.isParamSetByUser("positive_cumulative_direction"))
103 0 : mooseWarning(
104 : "The 'positive_cumulative_direction' parameter is unused when 'cumulative' is false");
105 :
106 4718 : if (_layered_base_params.isParamValid("num_layers"))
107 : {
108 4684 : _num_layers = _layered_base_params.get<unsigned int>("num_layers");
109 4684 : _interval_based = true;
110 : }
111 34 : else if (_layered_base_params.isParamValid("bounds"))
112 : {
113 30 : _interval_based = false;
114 :
115 30 : _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
116 :
117 : // Make sure the bounds are sorted - we're going to depend on this
118 30 : std::sort(_layer_bounds.begin(), _layer_bounds.end());
119 :
120 30 : _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
121 30 : _direction_min = _layer_bounds.front();
122 30 : _direction_max = _layer_bounds.back();
123 30 : _has_direction_max_min = true;
124 : }
125 : else
126 4 : mooseError("One of 'bounds' or 'num_layers' must be specified");
127 :
128 4714 : if (!_interval_based && _sample_type == 1)
129 4 : mooseError("'sample_type = interpolate' not supported with 'bounds'");
130 :
131 4710 : bool has_layer_bounding_block = _layered_base_params.isParamValid("layer_bounding_block");
132 4710 : bool has_block = _layered_base_params.isParamValid("block");
133 4710 : bool has_direction_min = _layered_base_params.isParamValid("direction_min");
134 4710 : bool has_direction_max = _layered_base_params.isParamValid("direction_max");
135 :
136 4710 : if (_has_direction_max_min && has_direction_min)
137 0 : mooseWarning("'direction_min' is unused when providing 'bounds'");
138 :
139 4710 : if (_has_direction_max_min && has_direction_max)
140 0 : mooseWarning("'direction_max' is unused when providing 'bounds'");
141 :
142 : // can only specify one of layer_bounding_block or the pair direction_max/min
143 4710 : if (has_layer_bounding_block && (has_direction_min || has_direction_max))
144 4 : mooseError("Only one of 'layer_bounding_block' and the pair 'direction_max' and "
145 : "'direction_min' can be provided");
146 :
147 : // if either one of direction_min or direction_max is specified, must provide the other one
148 4706 : if (has_direction_min != has_direction_max)
149 4 : mooseError("If providing the layer max/min directions, both 'direction_max' and "
150 : "'direction_min' must be specified.");
151 :
152 4702 : if (has_layer_bounding_block)
153 78 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
154 52 : _layered_base_params.get<std::vector<SubdomainName>>("layer_bounding_block"));
155 4676 : else if (has_block)
156 873 : _layer_bounding_blocks = _layered_base_subproblem.mesh().getSubdomainIDs(
157 582 : _layered_base_params.get<std::vector<SubdomainName>>("block"));
158 :
159 : // specifying the direction max/min overrides anything set with the 'block'
160 4702 : if (has_direction_min && has_direction_max)
161 : {
162 56 : _direction_min = parameters.get<Real>("direction_min");
163 56 : _direction_max = parameters.get<Real>("direction_max");
164 56 : _has_direction_max_min = true;
165 :
166 56 : if (_direction_max <= _direction_min)
167 4 : mooseError("'direction_max' must be larger than 'direction_min'");
168 : }
169 :
170 4698 : _layer_values.resize(_num_layers);
171 4698 : _layer_has_value.resize(_num_layers);
172 :
173 : // if we haven't already figured out the max/min in specified direction
174 : // (either with the 'bounds' or explicit specification from the user), do so
175 4698 : if (!_has_direction_max_min)
176 4620 : getBounds();
177 :
178 4698 : computeLayerCenters();
179 4698 : }
180 :
181 : Real
182 3209115 : LayeredBase::integralValue(const Point & p) const
183 : {
184 3209115 : unsigned int layer = getLayer(p);
185 :
186 3209115 : int higher_layer = -1;
187 3209115 : int lower_layer = -1;
188 :
189 3372086 : for (unsigned int i = layer; i < _layer_values.size(); i++)
190 : {
191 3326498 : if (_layer_has_value[i])
192 : {
193 3163527 : higher_layer = i;
194 3163527 : break;
195 : }
196 : }
197 :
198 3406989 : for (int i = layer - 1; i >= 0; i--)
199 : {
200 2558846 : if (_layer_has_value[i])
201 : {
202 2360972 : lower_layer = i;
203 2360972 : break;
204 : }
205 : }
206 :
207 3209115 : if (higher_layer == -1 && lower_layer == -1)
208 45444 : return 0; // TODO: We could error here but there are startup dependency problems
209 :
210 3163671 : switch (_sample_type)
211 : {
212 3157911 : case 0: // direct
213 : {
214 3157911 : if (higher_layer == -1) // Didn't find a higher layer
215 144 : return _layer_values[lower_layer];
216 :
217 3157767 : if (unsigned(higher_layer) == layer) // constant in a layer
218 3135785 : return _layer_values[higher_layer];
219 :
220 21982 : if (lower_layer == -1) // Didn't find a lower layer
221 2884 : return _layer_values[higher_layer];
222 :
223 19098 : return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
224 : }
225 0 : case 1: // interpolate
226 : {
227 0 : if (higher_layer == -1) // Didn't find a higher layer
228 0 : return _layer_values[lower_layer];
229 :
230 0 : Real layer_length = (_direction_max - _direction_min) / _num_layers;
231 0 : Real lower_coor = _direction_min;
232 0 : Real lower_value = 0;
233 0 : if (lower_layer != -1)
234 : {
235 0 : lower_coor += (lower_layer + 1) * layer_length;
236 0 : lower_value = _layer_values[lower_layer];
237 : }
238 :
239 : // Interpolate between the two points
240 0 : Real higher_value = _layer_values[higher_layer];
241 :
242 : // Linear interpolation
243 : return lower_value +
244 0 : (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
245 : }
246 5760 : case 2: // average
247 : {
248 5760 : Real total = 0;
249 5760 : unsigned int num_values = 0;
250 :
251 5760 : if (higher_layer != -1)
252 : {
253 16128 : for (const auto i : make_range(_average_radius))
254 : {
255 11520 : int current_layer = higher_layer + i;
256 :
257 11520 : if ((size_t)current_layer >= _layer_values.size())
258 1152 : break;
259 :
260 10368 : if (_layer_has_value[current_layer])
261 : {
262 10368 : total += _layer_values[current_layer];
263 10368 : num_values += 1;
264 : }
265 : }
266 : }
267 :
268 5760 : if (lower_layer != -1)
269 : {
270 12672 : for (const auto i : make_range(_average_radius))
271 : {
272 9216 : int current_layer = lower_layer - i;
273 :
274 9216 : if (current_layer < 0)
275 1152 : break;
276 :
277 8064 : if (_layer_has_value[current_layer])
278 : {
279 8064 : total += _layer_values[current_layer];
280 8064 : num_values += 1;
281 : }
282 : }
283 : }
284 :
285 5760 : return total / num_values;
286 : }
287 0 : default:
288 0 : mooseError("Unknown sample type!");
289 : }
290 : }
291 :
292 : Real
293 694627 : LayeredBase::getLayerValue(unsigned int layer) const
294 : {
295 694627 : if (layer >= _layer_values.size())
296 0 : mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
297 694627 : return _layer_values[layer];
298 : }
299 :
300 : void
301 9749 : LayeredBase::initialize()
302 : {
303 9749 : if (_using_displaced_mesh)
304 83 : getBounds();
305 :
306 67650 : for (const auto i : index_range(_layer_values))
307 : {
308 57901 : _layer_values[i] = 0.0;
309 57901 : _layer_has_value[i] = false;
310 : }
311 9749 : }
312 :
313 : void
314 8920 : LayeredBase::finalize()
315 : {
316 8920 : _layered_base_subproblem.comm().sum(_layer_values);
317 8920 : _layered_base_subproblem.comm().max(_layer_has_value);
318 :
319 8920 : if (_cumulative)
320 : {
321 422 : Real value = 0;
322 :
323 422 : if (_positive_cumulative_direction)
324 844 : for (const auto i : make_range(_num_layers))
325 : {
326 633 : value += getLayerValue(i);
327 633 : setLayerValue(i, value);
328 : }
329 : else
330 844 : for (int i = _num_layers - 1; i >= 0; i--)
331 : {
332 633 : value += getLayerValue(i);
333 633 : setLayerValue(i, value);
334 : }
335 : }
336 8920 : }
337 :
338 : void
339 805 : LayeredBase::threadJoin(const UserObject & y)
340 : {
341 805 : const auto & lb = dynamic_cast<const LayeredBase &>(y);
342 5866 : for (const auto i : index_range(_layer_values))
343 5061 : if (lb.layerHasValue(i))
344 2471 : setLayerValue(i, getLayerValue(i) + lb._layer_values[i]);
345 805 : }
346 :
347 : unsigned int
348 4441451 : LayeredBase::getLayer(const Point & p) const
349 : {
350 4441451 : Real direction_x = p(_direction);
351 :
352 4441451 : if (direction_x < _direction_min)
353 15695 : return 0;
354 :
355 4425756 : if (_interval_based)
356 : {
357 4420764 : unsigned int layer =
358 4420764 : std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
359 4420764 : static_cast<Real>(_num_layers));
360 :
361 4420764 : if (layer >= _num_layers)
362 68955 : layer = _num_layers - 1;
363 :
364 4420764 : return layer;
365 : }
366 : else // Figure out what layer we are in from the bounds
367 : {
368 : // This finds the first entry in the vector that is larger than what we're looking for
369 : std::vector<Real>::const_iterator one_higher =
370 4992 : std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
371 :
372 4992 : if (one_higher == _layer_bounds.end())
373 : {
374 : return static_cast<unsigned int>(
375 0 : _layer_bounds.size() -
376 0 : 2); // Just return the last layer. -2 because layers are "in-between" bounds
377 : }
378 4992 : else if (one_higher == _layer_bounds.begin())
379 0 : return 0; // Return the first layer
380 : else
381 : // The -1 is because the interval that we fall in is just _before_ the number that is bigger
382 : // (which is what we found
383 4992 : return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
384 : }
385 : }
386 :
387 : void
388 4698 : LayeredBase::computeLayerCenters()
389 : {
390 4698 : _layer_centers.resize(_num_layers);
391 :
392 4698 : if (_interval_based)
393 : {
394 4672 : Real dx = (_direction_max - _direction_min) / _num_layers;
395 :
396 41354 : for (const auto i : make_range(_num_layers))
397 36682 : _layer_centers[i] = (i + 0.5) * dx;
398 : }
399 : else
400 : {
401 91 : for (const auto i : make_range(_num_layers))
402 65 : _layer_centers[i] = 0.5 * (_layer_bounds[i + 1] + _layer_bounds[i]);
403 : }
404 4698 : }
405 :
406 : void
407 694847 : LayeredBase::setLayerValue(unsigned int layer, Real value)
408 : {
409 694847 : _layer_values[layer] = value;
410 694847 : _layer_has_value[layer] = true;
411 694847 : }
412 :
413 : void
414 4703 : LayeredBase::getBounds()
415 : {
416 4703 : if (_layer_bounding_blocks.size() == 0)
417 : {
418 4399 : BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
419 4399 : _direction_min = bounding_box.min()(_direction);
420 4399 : _direction_max = bounding_box.max()(_direction);
421 : }
422 : else
423 : {
424 304 : _direction_min = std::numeric_limits<Real>::infinity();
425 304 : _direction_max = -std::numeric_limits<Real>::infinity();
426 :
427 304 : MooseMesh & mesh = _layered_base_subproblem.mesh();
428 :
429 4358 : for (auto & elem_ptr : *mesh.getActiveLocalElementRange())
430 : {
431 4054 : auto subdomain_id = elem_ptr->subdomain_id();
432 :
433 4054 : if (std::find(_layer_bounding_blocks.begin(), _layer_bounding_blocks.end(), subdomain_id) ==
434 8108 : _layer_bounding_blocks.end())
435 2495 : continue;
436 :
437 6565 : for (auto & node : elem_ptr->node_ref_range())
438 : {
439 5006 : _direction_min = std::min(_direction_min, node(_direction));
440 5006 : _direction_max = std::max(_direction_max, node(_direction));
441 : }
442 : }
443 :
444 304 : mesh.comm().min(_direction_min);
445 304 : mesh.comm().max(_direction_max);
446 : }
447 4703 : }
|
