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 "MultiAppGeneralFieldTransfer.h"
11 :
12 : // MOOSE includes
13 : #include "DisplacedProblem.h"
14 : #include "FEProblem.h"
15 : #include "MooseMesh.h"
16 : #include "MooseTypes.h"
17 : #include "MooseVariableFE.h"
18 : #include "MeshDivision.h"
19 : #include "Positions.h"
20 : #include "Factory.h"
21 : #include "MooseAppCoordTransform.h"
22 : #include "MultiAppGeneralFieldFunctorTransfer.h"
23 :
24 : // libmesh includes
25 : #include "libmesh/point_locator_base.h"
26 : #include "libmesh/enum_point_locator_type.h"
27 :
28 : // TIMPI includes
29 : #include "timpi/communicator.h"
30 : #include "timpi/parallel_sync.h"
31 :
32 : using namespace libMesh;
33 :
34 : namespace GeneralFieldTransfer
35 : {
36 : Number OutOfMeshValue = std::numeric_limits<Real>::infinity();
37 : }
38 :
39 : InputParameters
40 30274 : MultiAppGeneralFieldTransfer::validParams()
41 : {
42 30274 : InputParameters params = MultiAppConservativeTransfer::validParams();
43 : // Expansion default to make a node in the target mesh overlapping with a node in the origin
44 : // mesh always register as being inside the origin application bounding box. The contains_point
45 : // bounding box checks uses exact comparisons
46 151370 : params.addRangeCheckedParam<Real>("bbox_factor",
47 60548 : 1.00000001,
48 : "bbox_factor>0",
49 : "Factor to inflate or deflate the source app bounding boxes");
50 181644 : params.addRangeCheckedParam<std::vector<Real>>(
51 : "fixed_bounding_box_size",
52 : "fixed_bounding_box_size >= 0",
53 : "Override source app bounding box size(s) for searches. App bounding boxes will still be "
54 : "centered on the same coordinates. Only non-zero components passed will override.");
55 90822 : params.addParam<Real>(
56 : "extrapolation_constant",
57 60548 : 0,
58 : "Constant to use when no source app can provide a valid value for a target location.");
59 121096 : MooseEnum extrap_options("none nearest-valid-target", "none");
60 121096 : params.addParam<MooseEnum>("post_transfer_extrapolation",
61 : extrap_options,
62 : "Post treatment to apply to the field after the transfer");
63 :
64 : // Block restrictions
65 121096 : params.addParam<std::vector<SubdomainName>>(
66 : "from_blocks",
67 : "Subdomain restriction to transfer from (defaults to all the origin app domain)");
68 121096 : params.addParam<std::vector<SubdomainName>>(
69 : "to_blocks", "Subdomain restriction to transfer to, (defaults to all the target app domain)");
70 :
71 : // Boundary restrictions
72 121096 : params.addParam<std::vector<BoundaryName>>(
73 : "from_boundaries",
74 : "The boundary we are transferring from (if not specified, whole domain is used).");
75 121096 : params.addParam<std::vector<BoundaryName>>(
76 : "to_boundaries",
77 : "The boundary we are transferring to (if not specified, whole domain is used).");
78 121096 : MooseEnum nodes_or_sides("nodes sides", "sides");
79 121096 : params.addParam<MooseEnum>("elemental_boundary_restriction",
80 : nodes_or_sides,
81 : "Whether elemental variable boundary restriction is considered by "
82 : "element side or element nodes");
83 :
84 : // Mesh division restriction
85 121096 : params.addParam<MeshDivisionName>("from_mesh_division",
86 : "Mesh division object on the origin application");
87 121096 : params.addParam<MeshDivisionName>("to_mesh_division",
88 : "Mesh division object on the target application");
89 : MooseEnum mesh_division_uses("spatial_restriction matching_division matching_subapp_index none",
90 121096 : "none");
91 121096 : params.addParam<MooseEnum>("from_mesh_division_usage",
92 : mesh_division_uses,
93 : "How to use the source mesh division in the transfer. See object "
94 : "documentation for description of each option");
95 121096 : params.addParam<MooseEnum>("to_mesh_division_usage",
96 : mesh_division_uses,
97 : "How to use the target mesh division in the transfer. See object "
98 : "documentation for description of each option");
99 :
100 : // Array and vector variables
101 90822 : params.addParam<std::vector<unsigned int>>("source_variable_components",
102 60548 : std::vector<unsigned int>(),
103 : "The source array or vector variable component(s).");
104 90822 : params.addParam<std::vector<unsigned int>>("target_variable_components",
105 60548 : std::vector<unsigned int>(),
106 : "The target array or vector variable component(s).");
107 :
108 : // Search options
109 90822 : params.addParam<bool>(
110 : "greedy_search",
111 60548 : false,
112 : "Whether or not to send a point to all the domains. If true, all the processors will be "
113 : "checked for a given point."
114 : "The code will be slow if this flag is on but it will give a better solution.");
115 90822 : params.addParam<bool>(
116 : "error_on_miss",
117 60548 : true,
118 : "Whether or not to error in the case that a target point is not found in the source domain.");
119 90822 : params.addParam<bool>("use_bounding_boxes",
120 60548 : true,
121 : "When set to false, bounding boxes will not be used to restrict the source "
122 : "of the transfer. Either source applications must be set using the "
123 : "from_mesh_division parameter, or a greedy search must be used.");
124 90822 : params.addParam<bool>(
125 : "use_nearest_app",
126 60548 : false,
127 : "When True, transfers from a child application will work by finding the nearest (using "
128 : "the `position` + mesh centroid) sub-app and query that app for the value to transfer.");
129 121096 : params.addParam<PositionsName>(
130 : "use_nearest_position",
131 : "Name of the the Positions object (in main app) such that transfers to/from a child "
132 : "application will work by finding the nearest position to a target and query only the "
133 : "app / points closer to this position than to any other position for the value to transfer.");
134 90822 : params.addParam<bool>(
135 : "from_app_must_contain_point",
136 60548 : false,
137 : "Wether on not the origin mesh must contain the point to evaluate data at. If false, this "
138 : "allows for interpolation between origin app meshes. Origin app bounding boxes are still "
139 : "considered so you may want to increase them with 'fixed_bounding_box_size'");
140 90822 : params.addParam<bool>("search_value_conflicts",
141 60548 : false,
142 : "Whether to look for potential conflicts between two valid and different "
143 : "source values for any target point");
144 90822 : params.addParam<unsigned int>(
145 : "value_conflicts_output",
146 60548 : 10,
147 : "Maximum number of conflicts to output if value-conflicts, from equidistant sources to a "
148 : "given transfer target location, search is turned on");
149 :
150 121096 : params.addParamNamesToGroup(
151 : "to_blocks from_blocks to_boundaries from_boundaries elemental_boundary_restriction "
152 : "from_mesh_division from_mesh_division_usage to_mesh_division to_mesh_division_usage",
153 : "Transfer spatial restriction");
154 121096 : params.addParamNamesToGroup(
155 : "greedy_search use_bounding_boxes use_nearest_app use_nearest_position "
156 : "search_value_conflicts",
157 : "Search algorithm");
158 121096 : params.addParamNamesToGroup("error_on_miss from_app_must_contain_point extrapolation_constant",
159 : "Extrapolation behavior");
160 90822 : params.addParamNamesToGroup("bbox_factor fixed_bounding_box_size", "Source app bounding box");
161 :
162 : // We need a level of ghosting to move elemental value one layer out when using post transfer
163 : // extrapolation. We need geometric for the element vertex average, and algebraic for the dof
164 : // value NOTE: we need this on the target mesh!
165 90822 : params.addRelationshipManager(
166 : "ElementSideNeighborLayers",
167 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC,
168 0 : [](const InputParameters & obj_params, InputParameters & rm_params)
169 : {
170 98223 : if (!obj_params.isParamValid("from_multi_app") &&
171 31173 : obj_params.get<MooseEnum>("post_transfer_extrapolation") == "nearest-valid-target")
172 0 : rm_params.set<unsigned short>("layers") = 2;
173 44700 : rm_params.set<bool>("use_displaced_mesh") = obj_params.get<bool>("displaced_target_mesh");
174 22350 : });
175 :
176 60548 : return params;
177 30274 : }
178 :
179 7486 : MultiAppGeneralFieldTransfer::MultiAppGeneralFieldTransfer(const InputParameters & parameters)
180 : : MultiAppConservativeTransfer(parameters),
181 7486 : _from_var_components(getParam<std::vector<unsigned int>>("source_variable_components")),
182 14972 : _to_var_components(getParam<std::vector<unsigned int>>("target_variable_components")),
183 14972 : _use_bounding_boxes(getParam<bool>("use_bounding_boxes")),
184 14972 : _use_nearest_app(getParam<bool>("use_nearest_app")),
185 7486 : _nearest_positions_obj(
186 14972 : isParamValid("use_nearest_position")
187 8032 : ? &_fe_problem.getPositionsObject(getParam<PositionsName>("use_nearest_position"))
188 : : nullptr),
189 14972 : _source_app_must_contain_point(getParam<bool>("from_app_must_contain_point")),
190 14972 : _from_mesh_division_behavior(getParam<MooseEnum>("from_mesh_division_usage")),
191 14972 : _to_mesh_division_behavior(getParam<MooseEnum>("to_mesh_division_usage")),
192 7486 : _elemental_boundary_restriction_on_sides(
193 14972 : getParam<MooseEnum>("elemental_boundary_restriction") == "sides"),
194 14972 : _greedy_search(getParam<bool>("greedy_search")),
195 14972 : _search_value_conflicts(getParam<bool>("search_value_conflicts")),
196 7486 : _already_output_search_value_conflicts(false),
197 14972 : _search_value_conflicts_max_log(getParam<unsigned int>("value_conflicts_output")),
198 14972 : _post_transfer_extrapolation(getParam<MooseEnum>("post_transfer_extrapolation")),
199 14972 : _error_on_miss(getParam<bool>("error_on_miss")),
200 14972 : _default_extrapolation_value(getParam<Real>("extrapolation_constant")),
201 14972 : _bbox_factor(getParam<Real>("bbox_factor")),
202 37430 : _fixed_bbox_size(isParamValid("fixed_bounding_box_size")
203 7486 : ? getParam<std::vector<Real>>("fixed_bounding_box_size")
204 37430 : : std::vector<Real>(3, 0))
205 : {
206 7486 : _var_size = _to_var_names.size();
207 11680 : if (_to_var_names.size() != _from_var_names.size() && !parameters.isPrivate("source_variable"))
208 6 : paramError("variable", "The number of variables to transfer to and from should be equal");
209 :
210 : // Check the parameters of the components of the array / vector variable
211 7483 : if (_from_var_names.size() != _from_var_components.size() && _from_var_components.size() > 0)
212 0 : paramError("source_variable_components",
213 : "This parameter must be equal to the number of source variables");
214 7483 : if (_to_var_names.size() != _to_var_components.size() && _to_var_components.size() > 0)
215 0 : paramError("target_variable_components",
216 : "This parameter must be equal to the number of target variables");
217 :
218 : // Make simple 'use_nearest_app' parameter rely on Positions
219 7483 : if (_use_nearest_app)
220 : {
221 141 : if (_nearest_positions_obj)
222 6 : paramError("use_nearest_app", "Cannot use nearest-app and nearest-position together");
223 138 : if (!hasFromMultiApp())
224 6 : paramError("use_nearest_app",
225 : "Should have a 'from_multiapp' when using the nearest-app informed search");
226 270 : auto pos_params = _app.getFactory().getValidParams("MultiAppPositions");
227 540 : pos_params.set<std::vector<MultiAppName>>("multiapps") = {getFromMultiApp()->name()};
228 405 : _fe_problem.addReporter("MultiAppPositions", "_created_for_" + name(), pos_params);
229 135 : _nearest_positions_obj = &_fe_problem.getPositionsObject("_created_for_" + name());
230 135 : }
231 :
232 : // Dont let users get wrecked by bounding boxes if it looks like they are trying to extrapolate
233 14873 : if (!_source_app_must_contain_point && _use_bounding_boxes &&
234 28861 : (_nearest_positions_obj || isParamSetByUser("from_app_must_contain_point")))
235 1938 : if (!isParamSetByUser("bbox_factor") && !isParamSetByUser("fixed_bounding_box_size"))
236 3 : mooseWarning(
237 : "Extrapolation (nearest-source options, outside-app source) parameters have "
238 : "been passed, but no subapp bounding box expansion parameters have been passed.");
239 :
240 7549 : if (!_use_bounding_boxes &&
241 7849 : (isParamValid("fixed_bounding_box_size") || isParamSetByUser("bbox_factor")))
242 0 : paramError("use_bounding_boxes",
243 : "Cannot pass additional bounding box parameters (sizes, expansion, etc) if we are "
244 : "not using bounding boxes");
245 7744 : }
246 :
247 : void
248 7309 : MultiAppGeneralFieldTransfer::initialSetup()
249 : {
250 7309 : MultiAppConservativeTransfer::initialSetup();
251 :
252 : // Use IDs for block and boundary restriction
253 : // Loop over all source problems
254 17278 : for (const auto i_from : index_range(_from_problems))
255 : {
256 9996 : const auto & from_moose_mesh = _from_problems[i_from]->mesh(_displaced_source_mesh);
257 29988 : if (isParamValid("from_blocks"))
258 : {
259 3438 : const auto & block_names = getParam<std::vector<SubdomainName>>("from_blocks");
260 :
261 3747 : for (const auto & b : block_names)
262 2031 : if (!MooseMeshUtils::hasSubdomainName(from_moose_mesh.getMesh(), b))
263 6 : paramError("from_blocks", "The block '", b, "' was not found in the mesh");
264 :
265 1716 : if (!block_names.empty())
266 : {
267 1716 : const auto ids = from_moose_mesh.getSubdomainIDs(block_names);
268 1716 : _from_blocks.insert(ids.begin(), ids.end());
269 1716 : }
270 : }
271 :
272 29979 : if (isParamValid("from_boundaries"))
273 : {
274 354 : const auto & boundary_names = getParam<std::vector<BoundaryName>>("from_boundaries");
275 507 : for (const auto & bn : boundary_names)
276 333 : if (!MooseMeshUtils::hasBoundaryName(from_moose_mesh.getMesh(), bn))
277 6 : paramError("from_boundaries", "The boundary '", bn, "' was not found in the mesh");
278 :
279 174 : if (!boundary_names.empty())
280 : {
281 174 : const auto boundary_ids = from_moose_mesh.getBoundaryIDs(boundary_names);
282 174 : _from_boundaries.insert(boundary_ids.begin(), boundary_ids.end());
283 174 : }
284 : }
285 :
286 29970 : if (isParamValid("from_mesh_division"))
287 : {
288 708 : const auto & mesh_div_name = getParam<MeshDivisionName>("from_mesh_division");
289 354 : _from_mesh_divisions.push_back(&_from_problems[i_from]->getMeshDivision(mesh_div_name));
290 : // Check that the behavior set makes sense
291 354 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::RESTRICTION)
292 : {
293 144 : if (_from_mesh_divisions[i_from]->coversEntireMesh())
294 0 : mooseInfo("'from_mesh_division_usage' is set to use a spatial restriction but the "
295 0 : "'from_mesh_division' for source app of global index " +
296 0 : std::to_string(getGlobalSourceAppIndex(i_from)) +
297 : " covers the entire mesh. Do not expect any restriction from a mesh "
298 : "division that covers the entire mesh");
299 : }
300 309 : else if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX &&
301 507 : !isParamValid("to_mesh_division"))
302 0 : paramError("to_mesh_division_usage",
303 : "Source mesh division cannot match target mesh division if no target mesh "
304 : "division is specified");
305 210 : else if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX)
306 : {
307 111 : if (!hasToMultiApp())
308 6 : paramError("from_mesh_division_usage",
309 : "Cannot match source mesh division index to target subapp index if there is "
310 : "only one target: the parent app (not a subapp)");
311 216 : else if (getToMultiApp()->numGlobalApps() !=
312 108 : _from_mesh_divisions[i_from]->getNumDivisions())
313 3 : mooseWarning("Attempting to match target subapp index with the number of source mesh "
314 3 : "divisions, which is " +
315 6 : std::to_string(_from_mesh_divisions[i_from]->getNumDivisions()) +
316 6 : " while there are " + std::to_string(getToMultiApp()->numGlobalApps()) +
317 : " target subapps");
318 105 : if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX)
319 : // We do not support it because it would require sending the point + target app index +
320 : // target app division index, and we only send the Point + one number
321 6 : paramError("from_mesh_division_usage",
322 : "We do not support using target subapp index for source division behavior and "
323 : "matching the division index for the target mesh division behavior.");
324 : }
325 99 : else if (_from_mesh_division_behavior == "none")
326 0 : paramError("from_mesh_division_usage", "User must specify a 'from_mesh_division_usage'");
327 : }
328 9636 : else if (_from_mesh_division_behavior != "none")
329 6 : paramError("from_mesh_division",
330 : "'from_mesh_division' must be specified if the usage method is specified");
331 : }
332 :
333 : // Loop over all target problems
334 17237 : for (const auto i_to : index_range(_to_problems))
335 : {
336 9973 : const auto & to_moose_mesh = _to_problems[i_to]->mesh(_displaced_target_mesh);
337 29919 : if (isParamValid("to_blocks"))
338 : {
339 2910 : const auto & block_names = getParam<std::vector<SubdomainName>>("to_blocks");
340 3219 : for (const auto & b : block_names)
341 1767 : if (!MooseMeshUtils::hasSubdomainName(to_moose_mesh.getMesh(), b))
342 6 : paramError("to_blocks", "The block '", b, "' was not found in the mesh");
343 :
344 1452 : if (!block_names.empty())
345 : {
346 1452 : const auto ids = to_moose_mesh.getSubdomainIDs(block_names);
347 1452 : _to_blocks.insert(ids.begin(), ids.end());
348 1452 : }
349 : }
350 :
351 29910 : if (isParamValid("to_boundaries"))
352 : {
353 3282 : const auto & boundary_names = getParam<std::vector<BoundaryName>>("to_boundaries");
354 5493 : for (const auto & bn : boundary_names)
355 3855 : if (!MooseMeshUtils::hasBoundaryName(to_moose_mesh.getMesh(), bn))
356 6 : paramError("to_boundaries", "The boundary '", bn, "' was not found in the mesh");
357 :
358 1638 : if (!boundary_names.empty())
359 : {
360 1638 : const auto boundary_ids = to_moose_mesh.getBoundaryIDs(boundary_names);
361 1638 : _to_boundaries.insert(boundary_ids.begin(), boundary_ids.end());
362 1638 : }
363 : }
364 :
365 29901 : if (isParamValid("to_mesh_division"))
366 : {
367 702 : const auto & mesh_div_name = getParam<MeshDivisionName>("to_mesh_division");
368 351 : _to_mesh_divisions.push_back(&_to_problems[i_to]->getMeshDivision(mesh_div_name));
369 : // Check that the behavior set makes sense
370 351 : if (_to_mesh_division_behavior == MeshDivisionTransferUse::RESTRICTION)
371 : {
372 144 : if (_to_mesh_divisions[i_to]->coversEntireMesh())
373 0 : mooseInfo("'to_mesh_division_usage' is set to use a spatial restriction but the "
374 0 : "'to_mesh_division' for target application of global index " +
375 0 : std::to_string(getGlobalSourceAppIndex(i_to)) +
376 : " covers the entire mesh. Do not expect any restriction from a mesh "
377 : "division that covers the entire mesh");
378 : }
379 207 : else if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX)
380 : {
381 297 : if (!isParamValid("from_mesh_division"))
382 0 : paramError("to_mesh_division_usage",
383 : "Target mesh division cannot match source mesh division if no source mesh "
384 : "division is specified");
385 198 : else if ((*_from_mesh_divisions.begin())->getNumDivisions() !=
386 99 : _to_mesh_divisions[i_to]->getNumDivisions())
387 3 : mooseWarning("Source and target mesh divisions do not have the same number of bins. If "
388 : "this is what you expect, please reach out to a MOOSE or app developer to "
389 : "ensure appropriate use");
390 : }
391 108 : else if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX)
392 : {
393 108 : if (!hasFromMultiApp())
394 6 : paramError(
395 : "to_mesh_division_usage",
396 : "Cannot match target mesh division index to source subapp index if there is only one "
397 : "source: the parent app (not a subapp)");
398 105 : else if (getFromMultiApp()->numGlobalApps() != _to_mesh_divisions[i_to]->getNumDivisions())
399 3 : mooseWarning("Attempting to match source subapp index with the number of target mesh "
400 3 : "divisions, which is " +
401 6 : std::to_string(_to_mesh_divisions[i_to]->getNumDivisions()) +
402 6 : " while there are " + std::to_string(getFromMultiApp()->numGlobalApps()) +
403 : " source subapps");
404 102 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX)
405 0 : paramError(
406 : "from_mesh_division_usage",
407 : "We do not support using source subapp index for the target division behavior and "
408 : "matching the division index for the source mesh division behavior.");
409 : }
410 0 : else if (_to_mesh_division_behavior == "none")
411 0 : paramError("to_mesh_division_usage", "User must specify a 'to_mesh_division_usage'");
412 : }
413 9616 : else if (_to_mesh_division_behavior != "none")
414 3 : paramError("to_mesh_division",
415 3 : "'to_mesh_division' must be specified if usage method '" +
416 6 : Moose::stringify(_to_mesh_division_behavior) + "' is specified");
417 : }
418 :
419 : // Check if components are set correctly if using an array variable
420 17206 : for (const auto i_from : index_range(_from_problems))
421 : {
422 17268 : for (const auto var_index : make_range(_from_var_names.size()))
423 : {
424 : MooseVariableFieldBase & from_var =
425 14652 : _from_problems[i_from]->getVariable(0,
426 7326 : _from_var_names[var_index],
427 : Moose::VarKindType::VAR_ANY,
428 : Moose::VarFieldType::VAR_FIELD_ANY);
429 7326 : if (from_var.count() > 1 && _from_var_components.empty())
430 6 : paramError("source_variable_components", "Component must be passed for an array variable");
431 7323 : if (_from_var_components.size() && from_var.count() < _from_var_components[var_index])
432 6 : paramError("source_variable_components",
433 : "Component passed is larger than size of variable");
434 : }
435 : }
436 17201 : for (const auto i_to : index_range(_to_problems))
437 : {
438 20183 : for (const auto var_index : make_range(_to_var_names.size()))
439 : {
440 : MooseVariableFieldBase & to_var =
441 20480 : _to_problems[i_to]->getVariable(0,
442 10240 : _to_var_names[var_index],
443 : Moose::VarKindType::VAR_ANY,
444 : Moose::VarFieldType::VAR_FIELD_ANY);
445 10240 : if (to_var.count() > 1 && _to_var_components.empty())
446 6 : paramError("target_variable_components", "Component must be passed for an array variable");
447 10237 : if (_to_var_components.size() && to_var.count() < _to_var_components[var_index])
448 6 : paramError("target_variable_components",
449 : "Component passed is larger than size of variable");
450 : }
451 : }
452 :
453 : // Cache some quantities to avoid having to get them on every transferred point
454 7252 : if (_to_problems.size())
455 : {
456 7234 : _to_variables.resize(_to_var_names.size());
457 14636 : for (const auto i_var : index_range(_to_var_names))
458 14804 : _to_variables[i_var] = &_to_problems[0]->getVariable(
459 7402 : 0, _to_var_names[i_var], Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_ANY);
460 : }
461 7252 : }
462 :
463 : void
464 59651 : MultiAppGeneralFieldTransfer::getAppInfo()
465 : {
466 59651 : MultiAppFieldTransfer::getAppInfo();
467 :
468 : // Create the point locators to locate evaluation points in the origin mesh(es)
469 59651 : _from_point_locators.resize(_from_problems.size());
470 124361 : for (const auto i_from : index_range(_from_problems))
471 : {
472 64710 : const auto & from_moose_mesh = _from_problems[i_from]->mesh(_displaced_source_mesh);
473 64710 : _from_point_locators[i_from] =
474 129420 : PointLocatorBase::build(TREE_LOCAL_ELEMENTS, from_moose_mesh.getMesh());
475 64710 : _from_point_locators[i_from]->enable_out_of_mesh_mode();
476 : }
477 59651 : }
478 :
479 : void
480 52342 : MultiAppGeneralFieldTransfer::execute()
481 : {
482 52342 : TIME_SECTION(
483 : "MultiAppGeneralFieldTransfer::execute()_" + name(), 5, "Transfer execution " + name());
484 52342 : getAppInfo();
485 :
486 : // Set up bounding boxes, etc
487 52342 : prepareToTransfer();
488 :
489 : // loop over the vector of variables and make the transfer one by one
490 104778 : for (const auto i : make_range(_var_size))
491 52476 : transferVariable(i);
492 :
493 52302 : postExecute();
494 52302 : }
495 :
496 : void
497 52342 : MultiAppGeneralFieldTransfer::prepareToTransfer()
498 : {
499 : // Get the bounding boxes for the "from" domains.
500 : // Clean up _from_bboxes from the previous transfer execution
501 52342 : _from_bboxes.clear();
502 :
503 : // NOTE: This ignores the app's bounding box inflation and padding
504 52342 : _from_bboxes = getRestrictedFromBoundingBoxes();
505 :
506 : // Expand bounding boxes. Some desired points might be excluded
507 : // without an expansion
508 52342 : extendBoundingBoxes(_bbox_factor, _from_bboxes);
509 :
510 : // Figure out how many "from" domains each processor owns.
511 52342 : _froms_per_proc.clear();
512 52342 : _froms_per_proc = getFromsPerProc();
513 :
514 : // Get the index for the first source app every processor owns
515 52342 : _global_app_start_per_proc = getGlobalStartAppPerProc();
516 :
517 : // No need to keep searching for conflicts if the mesh has not changed
518 52342 : if (_already_output_search_value_conflicts && !_displaced_source_mesh && !_displaced_target_mesh)
519 0 : _search_value_conflicts = false;
520 52342 : }
521 :
522 : void
523 52302 : MultiAppGeneralFieldTransfer::postExecute()
524 : {
525 52302 : MultiAppConservativeTransfer::postExecute();
526 52302 : if (_search_value_conflicts)
527 259 : _already_output_search_value_conflicts = true;
528 52302 : }
529 :
530 : void
531 52476 : MultiAppGeneralFieldTransfer::transferVariable(unsigned int i)
532 : {
533 : mooseAssert(i < _var_size, "The variable of index " << i << " does not exist");
534 :
535 : // Find outgoing target points
536 : // We need to know what points we need to send which processors
537 : // One processor will receive many points from many processors
538 : // One point may go to different processors
539 52476 : ProcessorToPointVec outgoing_points;
540 52476 : extractOutgoingPoints(i, outgoing_points);
541 :
542 52470 : if (_from_var_names.size() || dynamic_cast<MultiAppGeneralFieldFunctorTransfer *>(this))
543 51258 : prepareEvaluationOfInterpValues(i);
544 : else
545 1212 : prepareEvaluationOfInterpValues(-1);
546 :
547 : // Fill values and app ids for incoming points
548 : // We are responsible to compute values for these incoming points
549 : auto gather_functor =
550 81100 : [this, &i](processor_id_type /*pid*/,
551 : const std::vector<std::pair<Point, unsigned int>> & incoming_locations,
552 : std::vector<std::pair<Real, Real>> & outgoing_vals)
553 : {
554 81100 : outgoing_vals.resize(
555 : incoming_locations.size(),
556 : {GeneralFieldTransfer::OutOfMeshValue, GeneralFieldTransfer::OutOfMeshValue});
557 : // Evaluate interpolation values for these incoming points
558 81100 : evaluateInterpValues(i, incoming_locations, outgoing_vals);
559 133570 : };
560 :
561 104940 : DofobjectToInterpValVec dofobject_to_valsvec(_to_problems.size());
562 52470 : InterpCaches interp_caches(_to_problems.size(), getMaxToProblemsBBoxDimensions());
563 52470 : InterpCaches distance_caches(_to_problems.size(), getMaxToProblemsBBoxDimensions());
564 :
565 : // Copy data out to incoming_vals_ids
566 81100 : auto action_functor = [this, &i, &dofobject_to_valsvec, &interp_caches, &distance_caches](
567 : processor_id_type pid,
568 : const std::vector<std::pair<Point, unsigned int>> & my_outgoing_points,
569 : const std::vector<std::pair<Real, Real>> & incoming_vals)
570 : {
571 81100 : auto & pointInfoVec = _processor_to_pointInfoVec[pid];
572 :
573 : // Cache interpolation values for each dof object / points
574 81100 : cacheIncomingInterpVals(pid,
575 : i,
576 : pointInfoVec,
577 : my_outgoing_points,
578 : incoming_vals,
579 : dofobject_to_valsvec,
580 : interp_caches,
581 : distance_caches);
582 133570 : };
583 :
584 : // We assume incoming_vals_ids is ordered in the same way as outgoing_points
585 : // Hopefully, pull_parallel_vector_data will not mess up this
586 52470 : const std::pair<Real, Real> * ex = nullptr;
587 52470 : libMesh::Parallel::pull_parallel_vector_data(
588 52470 : comm(), outgoing_points, gather_functor, action_functor, ex);
589 :
590 : // Check for conflicts and overlaps from the maps that were built during the transfer
591 52470 : if (_search_value_conflicts)
592 287 : outputValueConflicts(i, dofobject_to_valsvec, distance_caches);
593 :
594 : // Set cached values into solution vector
595 52442 : setSolutionVectorValues(i, dofobject_to_valsvec, interp_caches);
596 :
597 : // Modify solution vector values (notably extrapolation options in functor transfer)
598 52436 : correctSolutionVectorValues(i, dofobject_to_valsvec, interp_caches);
599 52436 : }
600 :
601 : void
602 5360367 : MultiAppGeneralFieldTransfer::locatePointReceivers(const Point point,
603 : std::set<processor_id_type> & processors)
604 : {
605 : // Check which processors have apps that may include or be near this point
606 : // A point may be close enough to several problems, hosted on several processes
607 5360367 : bool found = false;
608 :
609 : // Additional process-restriction techniques we could use (TODOs):
610 : // - create a heuristic for using nearest-positions
611 : // - from_mesh_divisions could be polled for which divisions they possess on each
612 : // process, depending on the behavior chosen. This could limit potential senders.
613 : // This should be done ahead of this function call, for all points at once
614 :
615 : // Determine the apps which will be receiving points (then sending values) using various
616 : // heuristics
617 5360367 : if (_use_nearest_app)
618 : {
619 : // Find the nearest position for the point
620 7919 : const bool initial = _fe_problem.getCurrentExecuteOnFlag() == EXEC_INITIAL;
621 : // The apps form the nearest positions here, this is the index of the nearest app
622 7919 : const auto nearest_index = _nearest_positions_obj->getNearestPositionIndex(point, initial);
623 :
624 : // Find the apps that are nearest to the same position
625 : // Global search over all applications
626 18559 : for (processor_id_type i_proc = 0; i_proc < n_processors(); ++i_proc)
627 : {
628 : // We need i_from to correspond to the global app index
629 10640 : unsigned int from0 = _global_app_start_per_proc[i_proc];
630 30014 : for (unsigned int i_from = from0; i_from < from0 + _froms_per_proc[i_proc]; ++i_from)
631 : {
632 19374 : if (_greedy_search || _search_value_conflicts || i_from == nearest_index)
633 : {
634 10534 : processors.insert(i_proc);
635 10534 : found = true;
636 : }
637 : mooseAssert(i_from < getFromMultiApp()->numGlobalApps(), "We should not reach this");
638 : }
639 : }
640 : mooseAssert((getFromMultiApp()->numGlobalApps() < n_processors() || processors.size() == 1) ||
641 : _greedy_search || _search_value_conflicts,
642 : "Should only be one source processor when using more processors than source apps");
643 : }
644 5352448 : else if (_use_bounding_boxes)
645 : {
646 : // We examine all (global) bounding boxes and find the minimum of the maximum distances within a
647 : // bounding box from the point. This creates a sphere around the point of interest. Any app
648 : // with a bounding box that intersects this sphere (with a bboxMinDistance <
649 : // nearest_max_distance) will be considered a potential source
650 : // NOTE: This is a heuristic. We could try others
651 : // NOTE: from_bboxes are in the reference space, as is the point.
652 5350165 : Real nearest_max_distance = std::numeric_limits<Real>::max();
653 13335380 : for (const auto & bbox : _from_bboxes)
654 : {
655 7985215 : Real distance = bboxMaxDistance(point, bbox);
656 7985215 : if (distance < nearest_max_distance)
657 6589925 : nearest_max_distance = distance;
658 : }
659 :
660 5350165 : unsigned int from0 = 0;
661 12785553 : for (processor_id_type i_proc = 0; i_proc < n_processors();
662 7435388 : from0 += _froms_per_proc[i_proc], ++i_proc)
663 : // i_from here is a hybrid index based on the cumulative sum of the apps per processor
664 15420603 : for (unsigned int i_from = from0; i_from < from0 + _froms_per_proc[i_proc]; ++i_from)
665 : {
666 7985215 : Real distance = bboxMinDistance(point, _from_bboxes[i_from]);
667 : // We will not break here because we want to send a point to all possible source domains
668 8164185 : if (_greedy_search || distance <= nearest_max_distance ||
669 178970 : _from_bboxes[i_from].contains_point(point))
670 : {
671 7806963 : processors.insert(i_proc);
672 7806963 : found = true;
673 : }
674 : }
675 : }
676 : // Greedy search will contact every single processor. It's not scalable, but if there's valid data
677 : // on any subapp on any process, it will find it
678 2283 : else if (_greedy_search)
679 : {
680 2080 : found = true;
681 4940 : for (const auto i_proc : make_range(n_processors()))
682 2860 : processors.insert(i_proc);
683 : }
684 : // Since we indicated that we only wanted values from a subapp with the same global index as the
685 : // target mesh division, we might as well only communicate with the process that owns this app
686 403 : else if (!_to_mesh_divisions.empty() &&
687 200 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX)
688 : {
689 : // The target point could have a different index in each target mesh division. So on paper, we
690 : // would need to check all of them.
691 200 : auto saved_target_div = MooseMeshDivision::INVALID_DIVISION_INDEX;
692 400 : for (const auto i_to : index_range(_to_meshes))
693 : {
694 400 : const auto target_div = _to_mesh_divisions[i_to]->divisionIndex(
695 200 : _to_transforms[getGlobalTargetAppIndex(i_to)]->mapBack(point));
696 : // If it's the same division index, do not redo the search
697 200 : if (target_div == saved_target_div)
698 0 : continue;
699 : else
700 200 : saved_target_div = target_div;
701 :
702 : // Look for the processors owning a source-app with an index equal to the target mesh division
703 475 : for (const auto i_proc : make_range(n_processors()))
704 1075 : for (const auto i_from : make_range(_froms_per_proc[i_proc]))
705 800 : if (target_div == _global_app_start_per_proc[i_proc] + i_from)
706 : {
707 200 : processors.insert(i_proc);
708 200 : found = true;
709 : }
710 : }
711 : }
712 : else
713 3 : mooseError("No algorithm were selected to find which processes may send value data "
714 : "for a each target point. Please either specify using bounding boxes, "
715 : "greedy search, or to_mesh_division-based parameters");
716 :
717 : // Error out if we could not find this point when ask us to do so
718 5360364 : if (!found && _error_on_miss)
719 0 : mooseError(
720 : "Cannot find a source application to provide a value at point: ",
721 : point,
722 : " \n ",
723 : "It must be that mismatched meshes, between the source and target application, are being "
724 : "used.\nIf you are using the bounding boxes or nearest-app heuristics, or mesh-divisions, "
725 : "please consider using the greedy_search to confirm. Then consider choosing a different "
726 : "transfer type.\nThis check can be turned off by setting 'error_on_miss' to false. The "
727 : "'extrapolation_constant' parameter will be used to set the local value at missed points.");
728 5360364 : }
729 :
730 : void
731 5360367 : MultiAppGeneralFieldTransfer::cacheOutgoingPointInfo(const Point point,
732 : const dof_id_type dof_object_id,
733 : const unsigned int problem_id,
734 : ProcessorToPointVec & outgoing_points)
735 : {
736 5360367 : std::set<processor_id_type> processors;
737 : // Find which processors will receive point data so they can send back value data
738 : // The list can be larger than needed, depending on the heuristic / algorithm used to make
739 : // the call on whether a processor (and the apps it runs) should be involved
740 5360367 : processors.clear();
741 5360367 : locatePointReceivers(point, processors);
742 :
743 : // We need to send this location data to these processors so they can send back values
744 12766729 : for (const auto pid : processors)
745 : {
746 : // Select which from_mesh_division the source data must come from for this point
747 7406365 : unsigned int required_source_division = 0;
748 7406365 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX)
749 9452 : required_source_division = getGlobalTargetAppIndex(problem_id);
750 7396913 : else if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
751 14790966 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
752 7394053 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX)
753 7866 : required_source_division = _to_mesh_divisions[problem_id]->divisionIndex(
754 7866 : _to_transforms[getGlobalTargetAppIndex(problem_id)]->mapBack(point));
755 :
756 : // Skip if we already know we don't want the point
757 7406365 : if (required_source_division == MooseMeshDivision::INVALID_DIVISION_INDEX)
758 0 : continue;
759 :
760 : // Store outgoing information for every source process
761 7406365 : outgoing_points[pid].push_back(std::pair<Point, unsigned int>(point, required_source_division));
762 :
763 : // Store point information locally for processing received data
764 : // We can use these information when inserting values into the solution vector
765 : PointInfo pointinfo;
766 7406365 : pointinfo.problem_id = problem_id;
767 7406365 : pointinfo.dof_object_id = dof_object_id;
768 7406365 : _processor_to_pointInfoVec[pid].push_back(pointinfo);
769 : }
770 5360364 : }
771 :
772 : void
773 52476 : MultiAppGeneralFieldTransfer::extractOutgoingPoints(const unsigned int var_index,
774 : ProcessorToPointVec & outgoing_points)
775 : {
776 : // Get the variable name, with the accommodation for array/vector names
777 52476 : const auto & var_name = getToVarName(var_index);
778 :
779 : // Clean up the map from processor to pointInfo vector
780 : // This map should be consistent with outgoing_points
781 52476 : _processor_to_pointInfoVec.clear();
782 :
783 : // Loop over all problems
784 107298 : for (const auto i_to : index_range(_to_problems))
785 : {
786 54828 : const auto global_i_to = getGlobalTargetAppIndex(i_to);
787 :
788 : // libMesh EquationSystems
789 54828 : auto & es = getEquationSystem(*_to_problems[i_to], _displaced_target_mesh);
790 : // libMesh system that has this variable
791 54828 : System * to_sys = find_sys(es, var_name);
792 54828 : auto sys_num = to_sys->number();
793 54828 : auto var_num = _to_variables[var_index]->number();
794 54828 : auto & fe_type = _to_variables[var_index]->feType();
795 54828 : bool is_nodal = _to_variables[var_index]->isNodal();
796 :
797 : // Moose mesh
798 54828 : const auto & to_moose_mesh = _to_problems[i_to]->mesh(_displaced_target_mesh);
799 54828 : const auto & to_mesh = to_moose_mesh.getMesh();
800 :
801 : // We support more general variables via libMesh GenericProjector
802 54828 : if (fe_type.order > CONSTANT && !is_nodal)
803 : {
804 404 : GeneralFieldTransfer::RecordRequests<Number> f;
805 404 : GeneralFieldTransfer::RecordRequests<Gradient> g;
806 404 : GeneralFieldTransfer::NullAction<Number> nullsetter;
807 404 : const std::vector<unsigned int> varvec(1, var_num);
808 :
809 : libMesh::GenericProjector<GeneralFieldTransfer::RecordRequests<Number>,
810 : GeneralFieldTransfer::RecordRequests<Gradient>,
811 : Number,
812 : GeneralFieldTransfer::NullAction<Number>>
813 404 : request_gather(*to_sys, f, &g, nullsetter, varvec);
814 :
815 : // Defining only boundary values will not be enough to describe the variable, disallow it
816 404 : if (_to_boundaries.size() && (_to_variables[var_index]->getContinuity() == DISCONTINUOUS))
817 3 : mooseError("Higher order discontinuous elemental variables are not supported for "
818 : "target-boundary "
819 : "restricted transfers");
820 :
821 : // Not implemented as the target mesh division could similarly be cutting elements in an
822 : // arbitrary way with not enough requested points to describe the target variable
823 401 : if (!_to_mesh_divisions.empty() && !_to_mesh_divisions[i_to]->coversEntireMesh())
824 0 : mooseError("Higher order variable support not implemented for target mesh division "
825 : "unless the mesh is fully covered / indexed in the mesh division. This must be "
826 : "set programmatically in the MeshDivision object used.");
827 :
828 : // We dont look at boundary restriction, not supported for higher order target variables
829 : // Same for mesh divisions
830 401 : const auto & to_begin = _to_blocks.empty()
831 662 : ? to_mesh.active_local_elements_begin()
832 401 : : to_mesh.active_local_subdomain_set_elements_begin(_to_blocks);
833 :
834 401 : const auto & to_end = _to_blocks.empty()
835 662 : ? to_mesh.active_local_elements_end()
836 401 : : to_mesh.active_local_subdomain_set_elements_end(_to_blocks);
837 :
838 401 : ConstElemRange to_elem_range(to_begin, to_end);
839 :
840 401 : request_gather.project(to_elem_range);
841 :
842 401 : dof_id_type point_id = 0;
843 283601 : for (const Point & p : f.points_requested())
844 : // using the point number as a "dof_object_id" will serve to identify the point if we ever
845 : // rework interp/distance_cache into the dof_id_to_value maps
846 283200 : this->cacheOutgoingPointInfo(
847 283200 : (*_to_transforms[global_i_to])(p), point_id++, i_to, outgoing_points);
848 :
849 : // This is going to require more complicated transfer work
850 401 : if (!g.points_requested().empty())
851 0 : mooseError("We don't currently support variables with gradient degrees of freedom");
852 401 : }
853 54424 : else if (is_nodal)
854 : {
855 5271772 : for (const auto & node : to_mesh.local_node_ptr_range())
856 : {
857 : // Skip this node if the variable has no dofs at it.
858 5221047 : if (node->n_dofs(sys_num, var_num) < 1)
859 301176 : continue;
860 :
861 : // Skip if it is a block restricted transfer and current node does not have
862 : // specified blocks
863 4919871 : if (!_to_blocks.empty() && !inBlocks(_to_blocks, to_moose_mesh, node))
864 17076 : continue;
865 :
866 4902795 : if (!_to_boundaries.empty() && !onBoundaries(_to_boundaries, to_moose_mesh, node))
867 215568 : continue;
868 :
869 : // Skip if the node does not meet the target mesh division behavior
870 : // We cannot know from which app the data will come from so we cannot know
871 : // the source mesh division index and the source app global index
872 4687227 : if (!_to_mesh_divisions.empty() && _to_mesh_divisions[i_to]->divisionIndex(*node) ==
873 : MooseMeshDivision::INVALID_DIVISION_INDEX)
874 4400 : continue;
875 :
876 : // Cache point information
877 : // We will use this information later for setting values back to solution vectors
878 4682827 : cacheOutgoingPointInfo(
879 4682827 : (*_to_transforms[global_i_to])(*node), node->id(), i_to, outgoing_points);
880 50725 : }
881 : }
882 : else // Elemental, constant monomial
883 : {
884 3696 : for (const auto & elem :
885 462384 : as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
886 : {
887 : // Skip this element if the variable has no dofs at it.
888 454992 : if (elem->n_dofs(sys_num, var_num) < 1)
889 2512 : continue;
890 :
891 : // Skip if the element is not inside the block restriction
892 452480 : if (!_to_blocks.empty() && !inBlocks(_to_blocks, elem))
893 5668 : continue;
894 :
895 : // Skip if the element does not have a side on the boundary
896 446812 : if (!_to_boundaries.empty() && !onBoundaries(_to_boundaries, to_moose_mesh, elem))
897 49480 : continue;
898 :
899 : // Skip if the element is not indexed within the mesh division
900 397332 : if (!_to_mesh_divisions.empty() && _to_mesh_divisions[i_to]->divisionIndex(*elem) ==
901 : MooseMeshDivision::INVALID_DIVISION_INDEX)
902 2992 : continue;
903 :
904 : // Cache point information
905 : // We will use this information later for setting values back to solution vectors
906 394340 : cacheOutgoingPointInfo((*_to_transforms[global_i_to])(elem->vertex_average()),
907 394340 : elem->id(),
908 : i_to,
909 : outgoing_points);
910 3696 : } // for
911 : } // else
912 : } // for
913 52470 : }
914 :
915 : void
916 3934 : MultiAppGeneralFieldTransfer::extractLocalFromBoundingBoxes(std::vector<BoundingBox> & local_bboxes)
917 : {
918 3934 : local_bboxes.resize(_froms_per_proc[processor_id()]);
919 : // Find the index to the first of this processor's local bounding boxes.
920 3934 : unsigned int local_start = 0;
921 4984 : for (processor_id_type i_proc = 0; i_proc < n_processors() && i_proc != processor_id(); ++i_proc)
922 1050 : local_start += _froms_per_proc[i_proc];
923 :
924 : // Extract the local bounding boxes.
925 9078 : for (const auto i_from : make_range(_froms_per_proc[processor_id()]))
926 5144 : local_bboxes[i_from] = _from_bboxes[local_start + i_from];
927 3934 : }
928 :
929 : void
930 81100 : MultiAppGeneralFieldTransfer::cacheIncomingInterpVals(
931 : processor_id_type pid,
932 : const unsigned int var_index,
933 : std::vector<PointInfo> & pointInfoVec,
934 : const std::vector<std::pair<Point, unsigned int>> & point_requests,
935 : const std::vector<std::pair<Real, Real>> & incoming_vals,
936 : DofobjectToInterpValVec & dofobject_to_valsvec,
937 : InterpCaches & interp_caches,
938 : InterpCaches & distance_caches)
939 : {
940 : mooseAssert(pointInfoVec.size() == incoming_vals.size(),
941 : "Number of dof objects does not equal to the number of incoming values");
942 :
943 81100 : dof_id_type val_offset = 0;
944 7487465 : for (const auto & pointinfo : pointInfoVec)
945 : {
946 : // Retrieve target information from cached point infos
947 7406365 : const auto problem_id = pointinfo.problem_id;
948 7406365 : const auto dof_object_id = pointinfo.dof_object_id;
949 :
950 7406365 : auto & fe_type = _to_variables[var_index]->feType();
951 7406365 : bool is_nodal = _to_variables[var_index]->isNodal();
952 :
953 : // In the higher order elemental variable case, we receive point values, not nodal or
954 : // elemental. We use an InterpCache to store the values. The distance_cache is necessary to
955 : // choose between multiple origin problems sending values. This code could be unified with the
956 : // lower order order case by using the dofobject_to_valsvec
957 7406365 : if (fe_type.order > CONSTANT && !is_nodal)
958 : {
959 : // Cache solution on target mesh in its local frame of reference
960 383613 : InterpCache & value_cache = interp_caches[problem_id];
961 383613 : InterpCache & distance_cache = distance_caches[problem_id];
962 383613 : Point p = _to_transforms[getGlobalTargetAppIndex(problem_id)]->mapBack(
963 383613 : point_requests[val_offset].first);
964 383613 : const Number val = incoming_vals[val_offset].first;
965 :
966 : // Initialize distance to be able to compare
967 383613 : if (!distance_cache.hasKey(p))
968 283200 : distance_cache[p] = std::numeric_limits<Real>::max();
969 :
970 : // We should only have one closest value for each variable at any given point.
971 : // While there are shared Qps, on vertices for higher order variables usually,
972 : // the generic projector only queries each point once
973 383613 : if (_search_value_conflicts && !GeneralFieldTransfer::isOutOfMeshValue(val) &&
974 383613 : value_cache.hasKey(p) != 0 && !MooseUtils::absoluteFuzzyEqual(value_cache[p], val) &&
975 0 : MooseUtils::absoluteFuzzyEqual(distance_cache[p], incoming_vals[val_offset].second))
976 0 : registerConflict(problem_id, dof_object_id, p, incoming_vals[val_offset].second, false);
977 :
978 : // if we use the nearest app, even if the value is bad we want to save the distance because
979 : // it's the distance to the app, if that's the closest app then so be it with the bad value
980 467038 : if ((!GeneralFieldTransfer::isOutOfMeshValue(val) || _use_nearest_app) &&
981 83425 : MooseUtils::absoluteFuzzyGreaterThan(distance_cache[p], incoming_vals[val_offset].second))
982 : {
983 : // NOTE: We store the distance as well as the value. We really only need the
984 : // value to construct the variable, but the distance is used to make decisions in nearest
985 : // node schemes on which value to use
986 75700 : value_cache[p] = val;
987 75700 : distance_cache[p] = incoming_vals[val_offset].second;
988 : }
989 : }
990 : else
991 : {
992 : // Using the dof object pointer, so we can handle
993 : // both element and node using the same code
994 : #ifndef NDEBUG
995 : auto var_num = _to_variables[var_index]->number();
996 : auto & to_sys = _to_variables[var_index]->sys();
997 :
998 : const MeshBase & to_mesh = _to_problems[problem_id]->mesh(_displaced_target_mesh).getMesh();
999 : const DofObject * dof_object_ptr = nullptr;
1000 : const auto sys_num = to_sys.number();
1001 : // It is a node
1002 : if (is_nodal)
1003 : dof_object_ptr = to_mesh.node_ptr(dof_object_id);
1004 : // It is an element
1005 : else
1006 : dof_object_ptr = to_mesh.elem_ptr(dof_object_id);
1007 :
1008 : // We should only be supporting nodal and constant elemental
1009 : // variables in this code path; if we see multiple DoFs on one
1010 : // object we should have been using GenericProjector
1011 : mooseAssert(dof_object_ptr->n_dofs(sys_num, var_num) == 1,
1012 : "Unexpectedly found " << dof_object_ptr->n_dofs(sys_num, var_num)
1013 : << "dofs instead of 1");
1014 : #endif
1015 :
1016 7022752 : auto & dofobject_to_val = dofobject_to_valsvec[problem_id];
1017 :
1018 : // Check if we visited this dof object earlier
1019 7022752 : auto values_ptr = dofobject_to_val.find(dof_object_id);
1020 : // We did not visit this
1021 7022752 : if (values_ptr == dofobject_to_val.end())
1022 : {
1023 : // Values for this dof object
1024 5077164 : auto & val = dofobject_to_val[dof_object_id];
1025 : // Interpolation value
1026 5077164 : val.interp = incoming_vals[val_offset].first;
1027 : // Where this value came from
1028 5077164 : val.pid = pid;
1029 : // Distance
1030 5077164 : val.distance = incoming_vals[val_offset].second;
1031 : }
1032 : else
1033 : {
1034 1945588 : auto & val = values_ptr->second;
1035 :
1036 : // Look for value conflicts
1037 1945588 : if (detectConflict(val.interp,
1038 1945588 : incoming_vals[val_offset].first,
1039 : val.distance,
1040 1945588 : incoming_vals[val_offset].second))
1041 : {
1042 : // Keep track of distance and value
1043 : const Point p =
1044 0 : getPointInTargetAppFrame(point_requests[val_offset].first,
1045 : problem_id,
1046 : "Registration of received equi-distant value conflict");
1047 0 : registerConflict(problem_id, dof_object_id, p, incoming_vals[val_offset].second, false);
1048 : }
1049 :
1050 : // We adopt values that are, in order of priority
1051 : // - valid (or from nearest app)
1052 : // - closest distance
1053 : // - the smallest rank with the same distance
1054 : // It is debatable whether we want invalid values from the nearest app. It could just be
1055 : // that the app position was closer but the extent of another child app was large enough
1056 3891176 : if ((!GeneralFieldTransfer::isOutOfMeshValue(incoming_vals[val_offset].first) ||
1057 3667601 : _use_nearest_app) &&
1058 1722013 : (MooseUtils::absoluteFuzzyGreaterThan(val.distance, incoming_vals[val_offset].second) ||
1059 1382629 : ((val.pid > pid) &&
1060 594040 : MooseUtils::absoluteFuzzyEqual(val.distance, incoming_vals[val_offset].second))))
1061 : {
1062 342941 : val.interp = incoming_vals[val_offset].first;
1063 342941 : val.pid = pid;
1064 342941 : val.distance = incoming_vals[val_offset].second;
1065 : }
1066 : }
1067 : }
1068 :
1069 : // Move it to next position
1070 7406365 : val_offset++;
1071 : }
1072 81100 : }
1073 :
1074 : void
1075 363 : MultiAppGeneralFieldTransfer::registerConflict(
1076 : unsigned int problem, dof_id_type dof_id, Point p, Real dist, bool local)
1077 : {
1078 : // NOTE We could be registering the same conflict several times, we could count them instead
1079 363 : if (local)
1080 363 : _local_conflicts.push_back(std::make_tuple(problem, dof_id, p, dist));
1081 : else
1082 0 : _received_conflicts.push_back(std::make_tuple(problem, dof_id, p, dist));
1083 363 : }
1084 :
1085 : void
1086 287 : MultiAppGeneralFieldTransfer::examineReceivedValueConflicts(
1087 : const unsigned int var_index,
1088 : const DofobjectToInterpValVec & dofobject_to_valsvec,
1089 : const InterpCaches & distance_caches)
1090 : {
1091 287 : const auto var_name = getToVarName(var_index);
1092 : // We must check a posteriori because we could have two
1093 : // equidistant points with different values from two different problems, but a third point from
1094 : // another problem is actually closer, so there is no conflict because only that last one
1095 : // matters We check here whether the potential conflicts actually were the nearest points Loop
1096 : // over potential conflicts
1097 287 : for (auto conflict_it = _received_conflicts.begin(); conflict_it != _received_conflicts.end();)
1098 : {
1099 0 : const auto potential_conflict = *conflict_it;
1100 0 : bool overlap_found = false;
1101 :
1102 : // Extract info for the potential conflict
1103 0 : const unsigned int problem_id = std::get<0>(potential_conflict);
1104 0 : const dof_id_type dof_object_id = std::get<1>(potential_conflict);
1105 0 : const Point p = std::get<2>(potential_conflict);
1106 0 : const Real distance = std::get<3>(potential_conflict);
1107 :
1108 : // Extract target variable info
1109 0 : auto & es = getEquationSystem(*_to_problems[problem_id], _displaced_target_mesh);
1110 0 : System * to_sys = find_sys(es, var_name);
1111 0 : auto var_num = to_sys->variable_number(var_name);
1112 0 : auto & fe_type = to_sys->variable_type(var_num);
1113 0 : bool is_nodal = _to_variables[var_index]->isNodal();
1114 :
1115 : // Higher order elemental
1116 0 : if (fe_type.order > CONSTANT && !is_nodal)
1117 : {
1118 0 : auto cached_distance = distance_caches[problem_id].find(p);
1119 0 : if (cached_distance == distance_caches[problem_id].end())
1120 0 : mooseError("Conflict point was not found in the map of all origin-target distances");
1121 : // Distance is still the distance when we detected a potential overlap
1122 0 : if (MooseUtils::absoluteFuzzyEqual(cached_distance->second, distance))
1123 0 : overlap_found = true;
1124 : }
1125 : // Nodal and const monomial variable
1126 0 : else if (MooseUtils::absoluteFuzzyEqual(
1127 0 : dofobject_to_valsvec[problem_id].find(dof_object_id)->second.distance, distance))
1128 0 : overlap_found = true;
1129 :
1130 : // Map will only keep the actual overlaps
1131 0 : if (!overlap_found)
1132 0 : _received_conflicts.erase(conflict_it);
1133 : else
1134 0 : ++conflict_it;
1135 : }
1136 287 : }
1137 :
1138 : void
1139 287 : MultiAppGeneralFieldTransfer::examineLocalValueConflicts(
1140 : const unsigned int var_index,
1141 : const DofobjectToInterpValVec & dofobject_to_valsvec,
1142 : const InterpCaches & distance_caches)
1143 : {
1144 287 : const auto var_name = getToVarName(var_index);
1145 : // We must check a posteriori because we could have:
1146 : // - two equidistant points with different values from two different problems
1147 : // - two (or more) equidistant points with different values from the same problem
1148 : // but a third point/value couple from another problem is actually closer, so there is no
1149 : // conflict because only that last one matters. We check here whether the potential conflicts
1150 : // actually were the nearest points. We use several global reductions. If there are not too many
1151 : // potential conflicts (and there should not be in a well-posed problem) it should be manageably
1152 : // expensive
1153 :
1154 : // Move relevant conflict info (location, distance) to a smaller data structure
1155 287 : std::vector<std::tuple<Point, Real>> potential_conflicts;
1156 287 : potential_conflicts.reserve(_local_conflicts.size());
1157 :
1158 : // Loop over potential conflicts to broadcast all the conflicts
1159 650 : for (auto conflict_it = _local_conflicts.begin(); conflict_it != _local_conflicts.end();
1160 363 : ++conflict_it)
1161 : {
1162 : // Extract info for the potential conflict
1163 363 : const auto potential_conflict = *conflict_it;
1164 363 : const unsigned int i_from = std::get<0>(potential_conflict);
1165 363 : Point p = std::get<2>(potential_conflict);
1166 363 : const Real distance = std::get<3>(potential_conflict);
1167 : // If not using nearest-positions: potential conflict was saved in the source frame
1168 : // If using nearest-positions: potential conflict was saved in the reference frame
1169 363 : if (!_nearest_positions_obj)
1170 : {
1171 327 : const auto from_global_num = getGlobalSourceAppIndex(i_from);
1172 327 : p = (*_from_transforms[from_global_num])(p);
1173 : }
1174 :
1175 : // Send data in the global frame of reference
1176 363 : potential_conflicts.push_back(std::make_tuple(p, distance));
1177 : }
1178 287 : _communicator.allgather(potential_conflicts, false);
1179 : // conflicts could have been reported multiple times within a tolerance
1180 287 : std::sort(potential_conflicts.begin(), potential_conflicts.end());
1181 287 : potential_conflicts.erase(unique(potential_conflicts.begin(),
1182 : potential_conflicts.end(),
1183 335 : [](auto l, auto r)
1184 : {
1185 368 : return std::get<0>(l).absolute_fuzzy_equals(std::get<0>(r)) &&
1186 368 : std::abs(std::get<1>(l) - std::get<1>(r)) < TOLERANCE;
1187 : }),
1188 287 : potential_conflicts.end());
1189 :
1190 287 : std::vector<std::tuple<Point, Real>> real_conflicts;
1191 287 : real_conflicts.reserve(potential_conflicts.size());
1192 :
1193 : // For each potential conflict, we need to identify what problem asked for that value
1194 617 : for (auto conflict_it = potential_conflicts.begin(); conflict_it != potential_conflicts.end();
1195 330 : ++conflict_it)
1196 : {
1197 : // Extract info for the potential conflict
1198 330 : auto potential_conflict = *conflict_it;
1199 330 : const Point p = std::get<0>(potential_conflict);
1200 330 : const Real distance = std::get<1>(potential_conflict);
1201 :
1202 : // Check all the problems to try to find this requested point in the data structures filled
1203 : // with the received information
1204 330 : bool target_found = false;
1205 330 : bool conflict_real = false;
1206 660 : for (const auto i_to : index_range(_to_problems))
1207 : {
1208 : // Extract variable info
1209 330 : auto & es = getEquationSystem(*_to_problems[i_to], _displaced_target_mesh);
1210 330 : System * to_sys = find_sys(es, var_name);
1211 330 : auto var_num = to_sys->variable_number(var_name);
1212 330 : auto & fe_type = to_sys->variable_type(var_num);
1213 330 : bool is_nodal = _to_variables[var_index]->isNodal();
1214 :
1215 : // Move to the local frame of reference for the target problem
1216 : Point local_p =
1217 660 : getPointInTargetAppFrame(p, i_to, "Resolution of local value conflicts detected");
1218 :
1219 : // Higher order elemental
1220 330 : if (fe_type.order > CONSTANT && !is_nodal)
1221 : {
1222 : // distance_caches finds use a binned floating point search
1223 0 : auto cached_distance = distance_caches[i_to].find(local_p);
1224 0 : if (cached_distance != distance_caches[i_to].end())
1225 : {
1226 0 : target_found = true;
1227 : // Distance between source & target is still the distance we found in the sending
1228 : // process when we detected a potential overlap while gathering values to send
1229 0 : if (MooseUtils::absoluteFuzzyEqual(cached_distance->second, distance))
1230 0 : conflict_real = true;
1231 : }
1232 : }
1233 : // Nodal-value-dof-only and const monomial variable
1234 : else
1235 : {
1236 : // Find the dof id for the variable to be set
1237 330 : dof_id_type dof_object_id = std::numeric_limits<dof_id_type>::max();
1238 330 : auto pl = _to_problems[i_to]->mesh().getPointLocator();
1239 330 : pl->enable_out_of_mesh_mode();
1240 330 : if (is_nodal)
1241 : {
1242 309 : auto node = pl->locate_node(local_p);
1243 309 : if (node)
1244 : // this is not the dof_id for the variable, but the dof_object_id
1245 309 : dof_object_id = node->id();
1246 : }
1247 : else
1248 : {
1249 21 : auto elem = (*pl)(local_p);
1250 21 : if (elem)
1251 21 : dof_object_id = elem->id();
1252 : }
1253 330 : pl->disable_out_of_mesh_mode();
1254 :
1255 : // point isn't even in mesh
1256 330 : if (dof_object_id == std::numeric_limits<dof_id_type>::max())
1257 0 : continue;
1258 :
1259 : // this dof was not requested by this problem on this process
1260 330 : if (dofobject_to_valsvec[i_to].find(dof_object_id) == dofobject_to_valsvec[i_to].end())
1261 0 : continue;
1262 :
1263 330 : target_found = true;
1264 : // Check the saved distance in the vector of saved results. If the same, then the local
1265 : // conflict we detected with that distance is still an issue after receiving all values
1266 660 : if (MooseUtils::absoluteFuzzyEqual(
1267 330 : dofobject_to_valsvec[i_to].find(dof_object_id)->second.distance, distance))
1268 330 : conflict_real = true;
1269 330 : }
1270 : }
1271 : // Only keep the actual conflicts / overlaps
1272 330 : if (target_found && conflict_real)
1273 330 : real_conflicts.push_back(potential_conflict);
1274 : }
1275 :
1276 : // Communicate real conflicts to all so they can be checked by every process
1277 287 : _communicator.allgather(real_conflicts, false);
1278 :
1279 : // Delete potential conflicts that were resolved
1280 : // Each local list of conflicts will now be updated. It's important to keep conflict lists local
1281 : // so we can give more context like the sending processor id (the domain of which can be
1282 : // inspected by the user)
1283 650 : for (auto conflict_it = _local_conflicts.begin(); conflict_it != _local_conflicts.end();)
1284 : {
1285 : // Extract info for the potential conflict
1286 363 : const auto potential_conflict = *conflict_it;
1287 363 : const unsigned int i_from = std::get<0>(potential_conflict);
1288 363 : Point p = std::get<2>(potential_conflict);
1289 363 : const Real distance = std::get<3>(potential_conflict);
1290 363 : if (!_nearest_positions_obj)
1291 : {
1292 327 : const auto from_global_num = getGlobalSourceAppIndex(i_from);
1293 327 : p = (*_from_transforms[from_global_num])(p);
1294 : }
1295 :
1296 : // If not in the vector of real conflicts, was not real so delete it
1297 363 : if (std::find_if(real_conflicts.begin(),
1298 : real_conflicts.end(),
1299 5094 : [p, distance](const auto & item)
1300 : {
1301 5457 : return std::get<0>(item).absolute_fuzzy_equals(p) &&
1302 5457 : std::abs(std::get<1>(item) - distance) < TOLERANCE;
1303 726 : }) == real_conflicts.end())
1304 0 : _local_conflicts.erase(conflict_it);
1305 : else
1306 363 : ++conflict_it;
1307 : }
1308 287 : }
1309 :
1310 : void
1311 287 : MultiAppGeneralFieldTransfer::outputValueConflicts(
1312 : const unsigned int var_index,
1313 : const DofobjectToInterpValVec & dofobject_to_valsvec,
1314 : const InterpCaches & distance_caches)
1315 : {
1316 : // Remove potential conflicts that did not materialize, the value did not end up being used
1317 287 : examineReceivedValueConflicts(var_index, dofobject_to_valsvec, distance_caches);
1318 287 : examineLocalValueConflicts(var_index, dofobject_to_valsvec, distance_caches);
1319 :
1320 : // Output the conflicts from the selection of local values (evaluateInterpValues-type routines)
1321 : // to send in response to value requests at target points
1322 287 : const std::string rank_str = std::to_string(_communicator.rank());
1323 287 : if (_local_conflicts.size())
1324 : {
1325 28 : unsigned int num_outputs = 0;
1326 56 : std::string local_conflicts_string = "";
1327 : std::string potential_reasons =
1328 : "Are some points in target mesh equidistant from the sources "
1329 28 : "(nodes/centroids/apps/positions, depending on transfer) in origin mesh(es)?\n";
1330 28 : if (hasFromMultiApp() && _from_problems.size() > 1)
1331 22 : potential_reasons += "Are multiple subapps overlapping?\n";
1332 391 : for (const auto & conflict : _local_conflicts)
1333 : {
1334 363 : const unsigned int problem_id = std::get<0>(conflict);
1335 363 : Point p = std::get<2>(conflict);
1336 363 : num_outputs++;
1337 :
1338 363 : std::string origin_domain_message;
1339 363 : if (hasFromMultiApp() && !_nearest_positions_obj)
1340 : {
1341 : // NOTES:
1342 : // - The origin app for a conflict may not be unique.
1343 : // - The conflicts vectors only store the conflictual points, not the original one
1344 : // The original value found with a given distance could be retrieved from the main
1345 : // caches
1346 327 : const auto app_id = _from_local2global_map[problem_id];
1347 327 : origin_domain_message = "In source child app " + std::to_string(app_id) + " mesh,";
1348 : }
1349 : // We can't locate the source app when considering nearest positions, so we saved the data
1350 : // in the reference space. So we return the conflict location in the target app (parent or
1351 : // sibling) instead
1352 36 : else if (hasFromMultiApp() && _nearest_positions_obj)
1353 : {
1354 36 : if (_to_problems.size() == 1 || _skip_coordinate_collapsing)
1355 : {
1356 36 : p = (*_to_transforms[0])(p);
1357 36 : origin_domain_message = "In target app mesh,";
1358 : }
1359 : else
1360 0 : origin_domain_message = "In reference (post-coordinate collapse) mesh,";
1361 : }
1362 : else
1363 0 : origin_domain_message = "In source parent app mesh,";
1364 :
1365 363 : if (num_outputs < _search_value_conflicts_max_log)
1366 330 : local_conflicts_string += origin_domain_message + " point: (" + std::to_string(p(0)) +
1367 660 : ", " + std::to_string(p(1)) + ", " + std::to_string(p(2)) +
1368 660 : "), equi-distance: " + std::to_string(std::get<3>(conflict)) +
1369 165 : "\n";
1370 198 : else if (num_outputs == _search_value_conflicts_max_log)
1371 : local_conflicts_string +=
1372 : "Maximum output of the search for value conflicts has been reached. Further conflicts "
1373 6 : "will not be output.\nIncrease 'search_value_conflicts_max_log' to output more.";
1374 363 : }
1375 : // Explicitly name source to give more context
1376 28 : const std::string source_str = getDataSourceName(var_index);
1377 :
1378 56 : mooseWarning("On rank " + rank_str +
1379 : ", multiple valid values from equidistant points were "
1380 28 : "found in the origin mesh for source " +
1381 56 : source_str + " for " + std::to_string(_local_conflicts.size()) +
1382 28 : " target points.\n" + potential_reasons + "Conflicts detected at :\n" +
1383 : local_conflicts_string);
1384 0 : }
1385 :
1386 : // Output the conflicts discovered when receiving values from multiple origin problems
1387 259 : if (_received_conflicts.size())
1388 : {
1389 0 : unsigned int num_outputs = 0;
1390 0 : std::string received_conflicts_string = "";
1391 : std::string potential_reasons =
1392 : "Are some points in target mesh equidistant from the sources "
1393 0 : "(nodes/centroids/apps/positions, depending on transfer) in origin mesh(es)?\n";
1394 0 : if (hasToMultiApp() && _to_problems.size() > 1)
1395 0 : potential_reasons += "Are multiple subapps overlapping?\n";
1396 0 : for (const auto & conflict : _received_conflicts)
1397 : {
1398 : // Extract info for the potential overlap
1399 0 : const unsigned int problem_id = std::get<0>(conflict);
1400 0 : const Point p = std::get<2>(conflict);
1401 0 : num_outputs++;
1402 :
1403 0 : std::string target_domain_message;
1404 0 : if (hasToMultiApp())
1405 : {
1406 0 : const auto app_id = _to_local2global_map[problem_id];
1407 0 : target_domain_message = "In target child app " + std::to_string(app_id) + " mesh,";
1408 : }
1409 : else
1410 0 : target_domain_message = "In target parent app mesh,";
1411 :
1412 0 : if (num_outputs < _search_value_conflicts_max_log)
1413 0 : received_conflicts_string += target_domain_message + " point: (" + std::to_string(p(0)) +
1414 0 : ", " + std::to_string(p(1)) + ", " + std::to_string(p(2)) +
1415 0 : "), equi-distance: " + std::to_string(std::get<3>(conflict)) +
1416 0 : "\n";
1417 0 : else if (num_outputs == _search_value_conflicts_max_log)
1418 : received_conflicts_string +=
1419 : "Maximum output of the search for value conflicts has been reached. Further conflicts "
1420 0 : "will not be output.\nIncrease 'search_value_conflicts_max_log' to output more.";
1421 0 : }
1422 0 : mooseWarning("On rank " + rank_str +
1423 : ", multiple valid values from equidistant points were "
1424 0 : "received for target variable '" +
1425 0 : getToVarName(var_index) + "' for " + std::to_string(_received_conflicts.size()) +
1426 0 : " target points.\n" + potential_reasons + "Conflicts detected at :\n" +
1427 : received_conflicts_string);
1428 0 : }
1429 :
1430 259 : if (_local_conflicts.empty() && _received_conflicts.empty())
1431 : {
1432 777 : if (isParamSetByUser("search_value_conflict"))
1433 0 : mooseInfo("Automated diagnosis did not detect floating point indetermination in transfer");
1434 259 : else if (_to_problems.size() > 10 || _from_problems.size() > 10 || _communicator.size() > 10)
1435 0 : mooseInfo(
1436 : "Automated diagnosis did not detect any floating point indetermination in "
1437 : "the transfer. You may consider turning it off using `search_value_conflicts=false` "
1438 : "to improve performance/scalability.");
1439 : }
1440 :
1441 : // Reset the conflicts vectors, to be used for checking conflicts when transferring the next
1442 : // variable
1443 259 : _local_conflicts.clear();
1444 259 : _received_conflicts.clear();
1445 259 : }
1446 :
1447 : void
1448 52442 : MultiAppGeneralFieldTransfer::setSolutionVectorValues(
1449 : const unsigned int var_index,
1450 : const DofobjectToInterpValVec & dofobject_to_valsvec,
1451 : const InterpCaches & interp_caches)
1452 : {
1453 : // Get the variable name, with the accommodation for array/vector names
1454 52442 : const auto & var_name = getToVarName(var_index);
1455 :
1456 107230 : for (const auto problem_id : index_range(_to_problems))
1457 : {
1458 54794 : auto & dofobject_to_val = dofobject_to_valsvec[problem_id];
1459 :
1460 : // libMesh EquationSystems
1461 : // NOTE: we would expect to set variables from the displaced equation system here
1462 54794 : auto & es = getEquationSystem(*_to_problems[problem_id], false);
1463 :
1464 : // libMesh system
1465 54794 : System * to_sys = find_sys(es, var_name);
1466 :
1467 : // libMesh mesh
1468 54794 : const MeshBase & to_mesh = _to_problems[problem_id]->mesh(_displaced_target_mesh).getMesh();
1469 54794 : auto var_num = to_sys->variable_number(var_name);
1470 54794 : auto sys_num = to_sys->number();
1471 :
1472 54794 : auto & fe_type = _to_variables[var_index]->feType();
1473 54794 : bool is_nodal = _to_variables[var_index]->isNodal();
1474 :
1475 54794 : if (fe_type.order > CONSTANT && !is_nodal)
1476 : {
1477 : // We may need to use existing data values in places where the
1478 : // from app domain doesn't overlap
1479 401 : MeshFunction to_func(es, *to_sys->current_local_solution, to_sys->get_dof_map(), var_num);
1480 401 : to_func.init();
1481 :
1482 : GeneralFieldTransfer::CachedData<Number> f(
1483 401 : interp_caches[problem_id], to_func, _default_extrapolation_value);
1484 401 : libMesh::VectorSetAction<Number> setter(*to_sys->solution);
1485 401 : const std::vector<unsigned int> varvec(1, var_num);
1486 :
1487 : libMesh::GenericProjector<GeneralFieldTransfer::CachedData<Number>,
1488 : GeneralFieldTransfer::CachedData<Gradient>,
1489 : Number,
1490 : libMesh::VectorSetAction<Number>>
1491 401 : set_solution(*to_sys, f, nullptr, setter, varvec);
1492 :
1493 : // We dont look at boundary restriction, not supported for higher order target variables
1494 401 : const auto & to_begin = _to_blocks.empty()
1495 662 : ? to_mesh.active_local_elements_begin()
1496 401 : : to_mesh.active_local_subdomain_set_elements_begin(_to_blocks);
1497 :
1498 401 : const auto & to_end = _to_blocks.empty()
1499 662 : ? to_mesh.active_local_elements_end()
1500 401 : : to_mesh.active_local_subdomain_set_elements_end(_to_blocks);
1501 :
1502 401 : ConstElemRange active_local_elem_range(to_begin, to_end);
1503 :
1504 401 : set_solution.project(active_local_elem_range);
1505 401 : }
1506 : else
1507 : {
1508 5127704 : for (const auto & val_pair : dofobject_to_val)
1509 : {
1510 5073317 : const auto dof_object_id = val_pair.first;
1511 :
1512 5073317 : const DofObject * dof_object = nullptr;
1513 5073317 : if (is_nodal)
1514 4679754 : dof_object = to_mesh.node_ptr(dof_object_id);
1515 : else
1516 393563 : dof_object = to_mesh.elem_ptr(dof_object_id);
1517 :
1518 5073317 : const auto dof = dof_object->dof_number(sys_num, var_num, 0);
1519 5073317 : const auto val = val_pair.second.interp;
1520 :
1521 : // This will happen if meshes are mismatched
1522 5073317 : if (_error_on_miss && GeneralFieldTransfer::isOutOfMeshValue(val))
1523 : {
1524 : const auto target_location =
1525 6 : hasToMultiApp()
1526 6 : ? " on target app " + std::to_string(getGlobalTargetAppIndex(problem_id))
1527 15 : : " on parent app";
1528 6 : const auto info_msg = "\nThis check can be turned off by setting 'error_on_miss' to "
1529 : "false. The 'extrapolation_constant' parameter will be used to set "
1530 : "the local value at missed points.";
1531 6 : if (is_nodal)
1532 3 : mooseError("No source value for node ",
1533 : dof_object_id,
1534 : target_location,
1535 : " could be located. Node details:\n",
1536 3 : _to_meshes[problem_id]->nodePtr(dof_object_id)->get_info(),
1537 : "\n",
1538 : info_msg);
1539 : else
1540 3 : mooseError("No source value for element ",
1541 : dof_object_id,
1542 : target_location,
1543 : " could be located. Element details:\n",
1544 3 : _to_meshes[problem_id]->elemPtr(dof_object_id)->get_info(),
1545 : "\n",
1546 : info_msg);
1547 0 : }
1548 :
1549 : // We should not put garbage into our solution vector
1550 : // but it can be that we want to set it to a different value than what was already there
1551 : // for example: the source app has been displaced and was sending an indicator of its
1552 : // position
1553 5073311 : if (GeneralFieldTransfer::isOutOfMeshValue(val))
1554 : {
1555 566504 : if (!GeneralFieldTransfer::isOutOfMeshValue(_default_extrapolation_value))
1556 : {
1557 : // For nearest-valid-target, keep the out-of-mesh sentinel in the solution so
1558 : // that correctSolutionVectorValues can reliably identify which DOFs still need
1559 : // extrapolation. Writing _default_extrapolation_value here instead would make it
1560 : // impossible to distinguish a legitimately-transferred value that happens to equal
1561 : // the extrapolation constant from a DOF that never received data.
1562 566056 : const auto missing_value = _post_transfer_extrapolation == "nearest-valid-target"
1563 566056 : ? GeneralFieldTransfer::OutOfMeshValue
1564 566056 : : _default_extrapolation_value;
1565 566056 : to_sys->solution->set(dof, missing_value);
1566 : }
1567 566504 : continue;
1568 566504 : }
1569 4506807 : to_sys->solution->set(dof, val);
1570 : }
1571 : }
1572 :
1573 54788 : to_sys->solution->close();
1574 : // Sync local solutions
1575 54788 : to_sys->update();
1576 : }
1577 52436 : }
1578 :
1579 : bool
1580 1825856 : MultiAppGeneralFieldTransfer::acceptPointInOriginMesh(unsigned int i_from,
1581 : const std::vector<BoundingBox> & local_bboxes,
1582 : const Point & pt,
1583 : const unsigned int only_from_mesh_div,
1584 : Real & distance) const
1585 : {
1586 1825856 : if (_use_bounding_boxes && !local_bboxes[i_from].contains_point(pt))
1587 1484626 : return false;
1588 : else
1589 : {
1590 341230 : auto * pl = _from_point_locators[i_from].get();
1591 341230 : const auto from_global_num = getGlobalSourceAppIndex(i_from);
1592 : const auto transformed_pt =
1593 682460 : getPointInSourceAppFrame(pt, i_from, "Source point acceptance check");
1594 :
1595 : // Check point against source block restriction
1596 341230 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, pl, transformed_pt))
1597 12313 : return false;
1598 :
1599 : // Check point against source boundary restriction. Block restriction will speed up the search
1600 341034 : if (!_from_boundaries.empty() &&
1601 0 : !onBoundaries(_from_boundaries, _from_blocks, *_from_meshes[i_from], pl, transformed_pt))
1602 0 : return false;
1603 :
1604 : // Check point against the source mesh division
1605 347798 : if ((!_from_mesh_divisions.empty() || !_to_mesh_divisions.empty()) &&
1606 6764 : !acceptPointMeshDivision(transformed_pt, i_from, only_from_mesh_div))
1607 3796 : return false;
1608 :
1609 : // Get nearest position (often a subapp position) for the target point
1610 : // We want values from the child app that is closest to the same position as the target
1611 337238 : Point nearest_position_source;
1612 337238 : if (_nearest_positions_obj)
1613 : {
1614 14399 : const bool initial = _fe_problem.getCurrentExecuteOnFlag() == EXEC_INITIAL;
1615 : // The search for the nearest position is done in the reference frame
1616 14399 : const Point nearest_position = _nearest_positions_obj->getNearestPosition(pt, initial);
1617 14399 : nearest_position_source = _nearest_positions_obj->getNearestPosition(
1618 14399 : (*_from_transforms[from_global_num])(Point(0, 0, 0)), initial);
1619 :
1620 14399 : if (!_skip_coordinate_collapsing &&
1621 0 : _from_transforms[from_global_num]->hasNonTranslationTransformation())
1622 0 : mooseError("Rotation and scaling currently unsupported with nearest positions transfer.");
1623 :
1624 : // Compute distance to nearest position and nearest position source
1625 14399 : const Real distance_to_position_nearest_source = (pt - nearest_position_source).norm();
1626 14399 : const Real distance_to_nearest_position = (pt - nearest_position).norm();
1627 :
1628 : // Source (usually app position) is not closest to the same positions as the target, dont
1629 : // send values. We check the distance instead of the positions because if they are the same
1630 : // that means there's two equidistant positions and we would want to capture that as a "value
1631 : // conflict"
1632 14399 : if (!MooseUtils::absoluteFuzzyEqual(distance_to_position_nearest_source,
1633 : distance_to_nearest_position))
1634 8201 : return false;
1635 :
1636 : // Set the distance as the distance from the nearest position to the target point
1637 6198 : distance = distance_to_position_nearest_source;
1638 : }
1639 :
1640 : // Check that the app actually contains the origin point
1641 : // We dont need to check if we already found it in a block or a boundary
1642 329157 : if (_from_blocks.empty() && _from_boundaries.empty() && _source_app_must_contain_point &&
1643 120 : !inMesh(pl, transformed_pt))
1644 120 : return false;
1645 : }
1646 328917 : return true;
1647 : }
1648 :
1649 : void
1650 52436 : MultiAppGeneralFieldTransfer::correctSolutionVectorValues(
1651 : const unsigned int var_index,
1652 : const DofobjectToInterpValVec & dofobject_to_valsvec,
1653 : const InterpCaches & /*interp_caches*/)
1654 : {
1655 : // TODO: variable component support
1656 :
1657 : // Get the variable name, with the accommodation for array/vector names
1658 52436 : const auto & var_name = getToVarName(var_index);
1659 :
1660 107224 : for (const auto problem_id : index_range(_to_problems))
1661 : {
1662 54788 : auto & dofobject_to_val = dofobject_to_valsvec[problem_id];
1663 :
1664 : // libMesh EquationSystems
1665 : // NOTE: we would expect to set variables from the displaced equation system here
1666 54788 : auto & es = getEquationSystem(*_to_problems[problem_id], false);
1667 :
1668 : // libMesh system
1669 54788 : System * to_sys = find_sys(es, var_name);
1670 :
1671 : // libMesh mesh
1672 54788 : const MeshBase & to_mesh = _to_problems[problem_id]->mesh(_displaced_target_mesh).getMesh();
1673 54788 : auto var_num = to_sys->variable_number(var_name);
1674 54788 : auto sys_num = to_sys->number();
1675 :
1676 54788 : auto & fe_type = getToVariable(var_index)->feType();
1677 54788 : bool is_nodal = getToVariable(var_index)->isNodal();
1678 :
1679 : // We might need the synchronization of values that update provides
1680 : // to find the nearest target value
1681 : // NOTE: we are checking the buffers still for the values transfered, so we actually don't gain
1682 : // anything from ghosting We have to still work with buffers, how else do we know the source
1683 : // (from transferred buffers) or target (from all the points listed in buffers) are met
1684 54788 : if (_post_transfer_extrapolation == "nearest-valid-target")
1685 : {
1686 44 : if (fe_type.order > CONSTANT && !is_nodal)
1687 0 : paramError("post_transfer_extrapolation",
1688 : "Nearest-valid-target is not implemented for higher order elemental variables");
1689 : const auto & node_to_elem_map =
1690 44 : _to_problems[problem_id]->mesh(_displaced_target_mesh).nodeToElemMap();
1691 :
1692 732 : for (const auto & val_pair : dofobject_to_val)
1693 : {
1694 688 : const auto dof_object_id = val_pair.first;
1695 :
1696 : // Check that the value was out of bounds
1697 688 : const DofObject * dof_object = nullptr;
1698 688 : if (is_nodal)
1699 416 : dof_object = to_mesh.node_ptr(dof_object_id);
1700 : else
1701 272 : dof_object = to_mesh.elem_ptr(dof_object_id);
1702 688 : const auto dof = dof_object->dof_number(sys_num, var_num, 0);
1703 688 : const auto val = val_pair.second.interp;
1704 688 : if (GeneralFieldTransfer::isOutOfMeshValue(val))
1705 : {
1706 328 : Real nearest_value = 0.;
1707 328 : dof_id_type min_dist_id = std::numeric_limits<dof_id_type>::max();
1708 :
1709 : // Find the nearest valid value
1710 328 : if (is_nodal)
1711 : {
1712 184 : const auto node = to_mesh.node_ptr(dof_object_id);
1713 : // Find nearest node
1714 : // NOTE: we have access to a bunch of values now here, we could interpolate!
1715 184 : Real min_distance_sq = std::numeric_limits<Real>::max();
1716 568 : for (const auto & elem_id : libmesh_map_find(node_to_elem_map, node->id()))
1717 : {
1718 384 : const auto elem = to_mesh.elem_ptr(elem_id);
1719 1920 : for (const auto & elem_node : elem->node_ref_range())
1720 : {
1721 1536 : Real distance_sq = (Point(elem_node) - Point(*node)).norm_sq();
1722 : // Avoid using another bad value from a node which did not receive data
1723 : // Note: if the node is on another process ID, we can't obtain the value from a
1724 : // buffer here Note: we could seek from the solution vector instead BUT if we do
1725 : // that we may be ignoring source restrictions set to the transfer.
1726 : // Note: Target mesh restrictions are fine since we are picking from
1727 : // dofobject_to_val
1728 1536 : if (distance_sq < min_distance_sq && elem_node.id() != node->id())
1729 : {
1730 851 : if (auto it = dofobject_to_val.find(elem_node.id());
1731 1652 : it != dofobject_to_val.end() &&
1732 801 : !GeneralFieldTransfer::isOutOfMeshValue(it->second.interp))
1733 : {
1734 172 : min_distance_sq = distance_sq;
1735 172 : min_dist_id = elem_node.id();
1736 172 : nearest_value = it->second.interp;
1737 : }
1738 679 : else if (elem_node.n_dofs(sys_num, var_num) > 0)
1739 : {
1740 679 : const auto other_dof = elem_node.dof_number(sys_num, var_num, 0);
1741 : try
1742 : {
1743 : // setSolutionVectorValues leaves DOFs that did not receive a transfer
1744 : // value marked with OutOfMeshValue, so isOutOfMeshValue is sufficient
1745 : // to reject them here. DOFs that did receive data (even if the value
1746 : // equals _default_extrapolation_value) are accepted correctly.
1747 679 : if (const auto sol_val = (*to_sys->current_local_solution)(other_dof);
1748 679 : !GeneralFieldTransfer::isOutOfMeshValue(sol_val))
1749 : {
1750 20 : min_distance_sq = distance_sq;
1751 20 : min_dist_id = elem_node.id();
1752 20 : nearest_value = sol_val;
1753 : }
1754 : }
1755 0 : catch (...)
1756 : {
1757 : // Access in ghosted vector failed, just keep going
1758 0 : }
1759 : }
1760 : }
1761 : }
1762 : }
1763 : }
1764 : else
1765 : {
1766 144 : const auto elem = to_mesh.elem_ptr(dof_object_id);
1767 144 : Real min_distance_sq = std::numeric_limits<Real>::max();
1768 720 : for (const auto neigh : elem->neighbor_ptr_range())
1769 : {
1770 576 : if (!neigh || neigh == libMesh::remote_elem)
1771 168 : continue;
1772 408 : Real distance_sq = (neigh->vertex_average() - elem->vertex_average()).norm_sq();
1773 408 : if (distance_sq < min_distance_sq)
1774 : {
1775 352 : if (auto it = dofobject_to_val.find(neigh->id());
1776 681 : it != dofobject_to_val.end() &&
1777 329 : !GeneralFieldTransfer::isOutOfMeshValue(it->second.interp))
1778 : {
1779 79 : min_distance_sq = distance_sq;
1780 79 : min_dist_id = neigh->id();
1781 79 : nearest_value = it->second.interp;
1782 : }
1783 : // Access into ghosted solution vector. See comments for node
1784 273 : else if (neigh->n_dofs(sys_num, var_num) > 0)
1785 : {
1786 273 : const auto other_dof = neigh->dof_number(sys_num, var_num, 0);
1787 : try
1788 : {
1789 : // Same reasoning as the nodal branch: DOFs without transfer data carry
1790 : // OutOfMeshValue, so isOutOfMeshValue is the correct rejection criterion.
1791 273 : if (const auto sol_val = (*to_sys->current_local_solution)(other_dof);
1792 273 : !GeneralFieldTransfer::isOutOfMeshValue(sol_val))
1793 : {
1794 9 : nearest_value = sol_val;
1795 9 : min_distance_sq = distance_sq;
1796 9 : min_dist_id = neigh->id();
1797 : }
1798 : }
1799 0 : catch (...)
1800 : {
1801 : // Access in ghosted vector failed, just keep going
1802 0 : }
1803 : }
1804 : }
1805 : }
1806 : }
1807 328 : nearest_value = (min_dist_id != std::numeric_limits<dof_id_type>::max())
1808 328 : ? nearest_value
1809 : : _default_extrapolation_value;
1810 :
1811 328 : if (min_dist_id != std::numeric_limits<dof_id_type>::max())
1812 224 : to_sys->solution->set(dof, nearest_value);
1813 : else
1814 : {
1815 : // No valid neighbor was found; replace the out-of-mesh sentinel with the
1816 : // fallback value so the solution vector does not retain an invalid sentinel.
1817 104 : to_sys->solution->set(dof, _default_extrapolation_value);
1818 126 : flagSolutionWarning(
1819 : "Search for the valid target nearest from a target point for which no "
1820 : "values were found (and thus extrapolation is required) failed. This warning will "
1821 : "not be repeated on the console for further failures.");
1822 : }
1823 : }
1824 : }
1825 44 : to_sys->solution->close();
1826 : // Sync local solutions
1827 44 : to_sys->update();
1828 : }
1829 : }
1830 52436 : }
1831 :
1832 : bool
1833 120 : MultiAppGeneralFieldTransfer::inMesh(const PointLocatorBase * const pl, const Point & point) const
1834 : {
1835 : // Note: we do not take advantage of a potential block restriction of the mesh here. This is
1836 : // because we can avoid this routine by calling inBlocks() instead
1837 120 : const Elem * elem = (*pl)(point);
1838 120 : return (elem != nullptr);
1839 : }
1840 :
1841 : bool
1842 49736 : MultiAppGeneralFieldTransfer::inBlocks(const std::set<SubdomainID> & blocks,
1843 : const Elem * elem) const
1844 : {
1845 49736 : return blocks.find(elem->subdomain_id()) != blocks.end();
1846 : }
1847 :
1848 : bool
1849 30400 : MultiAppGeneralFieldTransfer::inBlocks(const std::set<SubdomainID> & blocks,
1850 : const MooseMesh & /* mesh */,
1851 : const Elem * elem) const
1852 : {
1853 30400 : return inBlocks(blocks, elem);
1854 : }
1855 :
1856 : bool
1857 45172 : MultiAppGeneralFieldTransfer::inBlocks(const std::set<SubdomainID> & blocks,
1858 : const MooseMesh & mesh,
1859 : const Node * node) const
1860 : {
1861 45172 : const auto & node_blocks = mesh.getNodeBlockIds(*node);
1862 45172 : std::set<SubdomainID> u;
1863 45172 : std::set_intersection(blocks.begin(),
1864 : blocks.end(),
1865 : node_blocks.begin(),
1866 : node_blocks.end(),
1867 : std::inserter(u, u.begin()));
1868 90344 : return !u.empty();
1869 45172 : }
1870 :
1871 : bool
1872 9324 : MultiAppGeneralFieldTransfer::inBlocks(const std::set<SubdomainID> & blocks,
1873 : const PointLocatorBase * const pl,
1874 : const Point & point) const
1875 : {
1876 9324 : const Elem * elem = (*pl)(point, &blocks);
1877 9324 : return (elem != nullptr);
1878 : }
1879 :
1880 : bool
1881 512656 : MultiAppGeneralFieldTransfer::onBoundaries(const std::set<BoundaryID> & boundaries,
1882 : const MooseMesh & mesh,
1883 : const Node * node) const
1884 : {
1885 512656 : const BoundaryInfo & bnd_info = mesh.getMesh().get_boundary_info();
1886 512656 : std::vector<BoundaryID> vec_to_fill;
1887 512656 : bnd_info.boundary_ids(node, vec_to_fill);
1888 512656 : std::set<BoundaryID> vec_to_fill_set(vec_to_fill.begin(), vec_to_fill.end());
1889 512656 : std::set<BoundaryID> u;
1890 512656 : std::set_intersection(boundaries.begin(),
1891 : boundaries.end(),
1892 : vec_to_fill_set.begin(),
1893 : vec_to_fill_set.end(),
1894 : std::inserter(u, u.begin()));
1895 1025312 : return !u.empty();
1896 512656 : }
1897 :
1898 : bool
1899 88612 : MultiAppGeneralFieldTransfer::onBoundaries(const std::set<BoundaryID> & boundaries,
1900 : const MooseMesh & mesh,
1901 : const Elem * elem) const
1902 : {
1903 : // Get all boundaries each side of the element is part of
1904 88612 : const BoundaryInfo & bnd_info = mesh.getMesh().get_boundary_info();
1905 88612 : std::vector<BoundaryID> vec_to_fill;
1906 88612 : std::vector<BoundaryID> vec_to_fill_temp;
1907 88612 : if (_elemental_boundary_restriction_on_sides)
1908 536284 : for (const auto side : make_range(elem->n_sides()))
1909 : {
1910 459672 : bnd_info.boundary_ids(elem, side, vec_to_fill_temp);
1911 459672 : vec_to_fill.insert(vec_to_fill.end(), vec_to_fill_temp.begin(), vec_to_fill_temp.end());
1912 : }
1913 : else
1914 108000 : for (const auto node_index : make_range(elem->n_nodes()))
1915 : {
1916 96000 : bnd_info.boundary_ids(elem->node_ptr(node_index), vec_to_fill_temp);
1917 96000 : vec_to_fill.insert(vec_to_fill.end(), vec_to_fill_temp.begin(), vec_to_fill_temp.end());
1918 : }
1919 88612 : std::set<BoundaryID> vec_to_fill_set(vec_to_fill.begin(), vec_to_fill.end());
1920 :
1921 : // Look for a match between the boundaries from the restriction and those near the element
1922 88612 : std::set<BoundaryID> u;
1923 88612 : std::set_intersection(boundaries.begin(),
1924 : boundaries.end(),
1925 : vec_to_fill_set.begin(),
1926 : vec_to_fill_set.end(),
1927 : std::inserter(u, u.begin()));
1928 177224 : return !u.empty();
1929 88612 : }
1930 :
1931 : bool
1932 0 : MultiAppGeneralFieldTransfer::onBoundaries(const std::set<BoundaryID> & boundaries,
1933 : const std::set<SubdomainID> & block_restriction,
1934 : const MooseMesh & mesh,
1935 : const PointLocatorBase * const pl,
1936 : const Point & point) const
1937 : {
1938 : // Find the element containing the point and use the block restriction if known for speed
1939 : const Elem * elem;
1940 0 : if (block_restriction.empty())
1941 0 : elem = (*pl)(point);
1942 : else
1943 0 : elem = (*pl)(point, &block_restriction);
1944 :
1945 0 : if (!elem)
1946 0 : return false;
1947 0 : return onBoundaries(boundaries, mesh, elem);
1948 : }
1949 :
1950 : bool
1951 6764 : MultiAppGeneralFieldTransfer::acceptPointMeshDivision(
1952 : const Point & pt, const unsigned int i_local, const unsigned int only_from_this_mesh_div) const
1953 : {
1954 : // This routine can also be called to examine if the to_mesh_division index matches the current
1955 : // source subapp index
1956 6764 : unsigned int source_mesh_div = MooseMeshDivision::INVALID_DIVISION_INDEX - 1;
1957 6764 : if (!_from_mesh_divisions.empty())
1958 4824 : source_mesh_div = _from_mesh_divisions[i_local]->divisionIndex(pt);
1959 :
1960 : // If the point is not indexed in the source division
1961 6764 : if (!_from_mesh_divisions.empty() && source_mesh_div == MooseMeshDivision::INVALID_DIVISION_INDEX)
1962 2896 : return false;
1963 : // If the point is not the at the same index in the target and the origin meshes, reject
1964 3868 : else if ((_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
1965 3868 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX) &&
1966 : source_mesh_div != only_from_this_mesh_div)
1967 0 : return false;
1968 : // If the point is at a certain division index that is not the same as the index of the subapp
1969 : // we wanted the information to be from for that point, reject
1970 3868 : else if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
1971 : source_mesh_div != only_from_this_mesh_div)
1972 462 : return false;
1973 4044 : else if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
1974 638 : only_from_this_mesh_div != getGlobalSourceAppIndex(i_local))
1975 438 : return false;
1976 : else
1977 2968 : return true;
1978 : }
1979 :
1980 : bool
1981 146289 : MultiAppGeneralFieldTransfer::closestToPosition(unsigned int pos_index, const Point & pt) const
1982 : {
1983 : mooseAssert(_nearest_positions_obj, "Should not be here without a positions object");
1984 146289 : if (!_skip_coordinate_collapsing)
1985 0 : paramError("skip_coordinate_collapsing", "Coordinate collapsing not implemented");
1986 146289 : bool initial = _fe_problem.getCurrentExecuteOnFlag() == EXEC_INITIAL;
1987 146289 : if (!_search_value_conflicts)
1988 : // Faster to just compare the index
1989 17017 : return pos_index == _nearest_positions_obj->getNearestPositionIndex(pt, initial);
1990 : else
1991 : {
1992 : // Get the distance to the position and see if we are missing a value just because the position
1993 : // is not officially the closest, but it is actually at the same distance
1994 129272 : const auto nearest_position = _nearest_positions_obj->getNearestPosition(pt, initial);
1995 129272 : const auto nearest_position_at_index = _nearest_positions_obj->getPosition(pos_index, initial);
1996 129272 : Real distance_to_position_at_index = (pt - nearest_position_at_index).norm();
1997 129272 : const Real distance_to_nearest_position = (pt - nearest_position).norm();
1998 :
1999 129272 : if (!MooseUtils::absoluteFuzzyEqual(distance_to_position_at_index,
2000 : distance_to_nearest_position))
2001 84053 : return false;
2002 : // Actually the same position (point)
2003 45219 : else if (nearest_position == nearest_position_at_index)
2004 45219 : return true;
2005 : else
2006 : {
2007 0 : mooseWarning("Two equidistant positions ",
2008 : nearest_position,
2009 : " and ",
2010 : nearest_position_at_index,
2011 : " detected near point ",
2012 : pt);
2013 0 : return true;
2014 : }
2015 : }
2016 : }
2017 :
2018 : Real
2019 7985215 : MultiAppGeneralFieldTransfer::bboxMaxDistance(const Point & p, const BoundingBox & bbox) const
2020 : {
2021 7985215 : std::array<Point, 2> source_points = {{bbox.first, bbox.second}};
2022 :
2023 7985215 : std::array<Point, 8> all_points;
2024 23955645 : for (unsigned int x = 0; x < 2; x++)
2025 47911290 : for (unsigned int y = 0; y < 2; y++)
2026 95822580 : for (unsigned int z = 0; z < 2; z++)
2027 63881720 : all_points[x + 2 * y + 4 * z] =
2028 127763440 : Point(source_points[x](0), source_points[y](1), source_points[z](2));
2029 :
2030 7985215 : Real max_distance = 0.;
2031 :
2032 71866935 : for (unsigned int i = 0; i < 8; i++)
2033 : {
2034 63881720 : Real distance = (p - all_points[i]).norm();
2035 63881720 : if (distance > max_distance)
2036 17095062 : max_distance = distance;
2037 : }
2038 :
2039 7985215 : return max_distance;
2040 : }
2041 :
2042 : Real
2043 7985215 : MultiAppGeneralFieldTransfer::bboxMinDistance(const Point & p, const BoundingBox & bbox) const
2044 : {
2045 7985215 : std::array<Point, 2> source_points = {{bbox.first, bbox.second}};
2046 :
2047 7985215 : std::array<Point, 8> all_points;
2048 23955645 : for (unsigned int x = 0; x < 2; x++)
2049 47911290 : for (unsigned int y = 0; y < 2; y++)
2050 95822580 : for (unsigned int z = 0; z < 2; z++)
2051 63881720 : all_points[x + 2 * y + 4 * z] =
2052 127763440 : Point(source_points[x](0), source_points[y](1), source_points[z](2));
2053 :
2054 7985215 : Real min_distance = std::numeric_limits<Real>::max();
2055 :
2056 71866935 : for (unsigned int i = 0; i < 8; i++)
2057 : {
2058 63881720 : Real distance = (p - all_points[i]).norm();
2059 63881720 : if (distance < min_distance)
2060 17632467 : min_distance = distance;
2061 : }
2062 :
2063 7985215 : return min_distance;
2064 : }
2065 :
2066 : std::vector<BoundingBox>
2067 52342 : MultiAppGeneralFieldTransfer::getRestrictedFromBoundingBoxes() const
2068 : {
2069 52342 : std::vector<std::pair<Point, Point>> bb_points(_from_meshes.size());
2070 52342 : const Real min_r = std::numeric_limits<Real>::lowest();
2071 52342 : const Real max_r = std::numeric_limits<Real>::max();
2072 :
2073 107032 : for (const auto j : make_range(_from_meshes.size()))
2074 : {
2075 54690 : Point min(max_r, max_r, max_r);
2076 54690 : Point max(min_r, min_r, min_r);
2077 54690 : bool at_least_one = false;
2078 54690 : const auto & from_mesh = _from_problems[j]->mesh(_displaced_source_mesh);
2079 :
2080 109380 : for (const auto & elem : as_range(from_mesh.getMesh().local_elements_begin(),
2081 4760836 : from_mesh.getMesh().local_elements_end()))
2082 : {
2083 4596766 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, from_mesh, elem))
2084 20496 : continue;
2085 :
2086 26492626 : for (const auto & node : elem->node_ref_range())
2087 : {
2088 21916356 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, &node))
2089 161592 : continue;
2090 :
2091 21754764 : at_least_one = true;
2092 87019056 : for (const auto i : make_range(LIBMESH_DIM))
2093 : {
2094 65264292 : min(i) = std::min(min(i), node(i));
2095 65264292 : max(i) = std::max(max(i), node(i));
2096 : }
2097 : }
2098 54690 : }
2099 :
2100 : // For 2D RZ problems, we need to amend the bounding box to cover the whole XYZ projection
2101 : // - The XYZ-Y axis is assumed aligned with the RZ-Z axis
2102 : // - RZ systems also cover negative coordinates hence the use of the maximum R
2103 : // NOTE: We will only support the case where there is only one coordinate system
2104 54690 : if ((from_mesh.getUniqueCoordSystem() == Moose::COORD_RZ) && (LIBMESH_DIM == 3))
2105 : {
2106 119 : min(0) = -max(0);
2107 119 : min(2) = -max(0);
2108 119 : max(2) = max(0);
2109 : }
2110 :
2111 54690 : BoundingBox bbox(min, max);
2112 54690 : if (!at_least_one)
2113 31 : bbox.min() = max; // If we didn't hit any nodes, this will be _the_ minimum bbox
2114 : else
2115 : {
2116 : // Translate the bounding box to the from domain's position. We may have rotations so we
2117 : // must be careful in constructing the new min and max (first and second)
2118 54659 : const auto from_global_num = getGlobalSourceAppIndex(j);
2119 54659 : transformBoundingBox(bbox, *_from_transforms[from_global_num]);
2120 : }
2121 :
2122 : // Cast the bounding box into a pair of points (so it can be put through
2123 : // MPI communication).
2124 54690 : bb_points[j] = static_cast<std::pair<Point, Point>>(bbox);
2125 : }
2126 :
2127 : // Serialize the bounding box points.
2128 52342 : _communicator.allgather(bb_points);
2129 :
2130 : // Recast the points back into bounding boxes and return.
2131 52342 : std::vector<BoundingBox> bboxes(bb_points.size());
2132 136401 : for (const auto i : make_range(bb_points.size()))
2133 84059 : bboxes[i] = static_cast<BoundingBox>(bb_points[i]);
2134 :
2135 : // TODO move up
2136 : // Check for a user-set fixed bounding box size and modify the sizes as appropriate
2137 52342 : if (_fixed_bbox_size != std::vector<Real>(3, 0))
2138 1260 : for (const auto i : make_range(LIBMESH_DIM))
2139 945 : if (!MooseUtils::absoluteFuzzyEqual(_fixed_bbox_size[i], 0))
2140 1986 : for (const auto j : make_range(bboxes.size()))
2141 : {
2142 1245 : const auto current_width = (bboxes[j].second - bboxes[j].first)(i);
2143 1245 : bboxes[j].first(i) -= (_fixed_bbox_size[i] - current_width) / 2;
2144 1245 : bboxes[j].second(i) += (_fixed_bbox_size[i] - current_width) / 2;
2145 : }
2146 :
2147 104684 : return bboxes;
2148 52342 : }
2149 :
2150 : std::vector<unsigned int>
2151 52342 : MultiAppGeneralFieldTransfer::getGlobalStartAppPerProc() const
2152 : {
2153 52342 : std::vector<unsigned int> global_app_start_per_proc(1, -1);
2154 52342 : if (_from_local2global_map.size())
2155 27537 : global_app_start_per_proc[0] = _from_local2global_map[0];
2156 52342 : _communicator.allgather(global_app_start_per_proc, true);
2157 52342 : return global_app_start_per_proc;
2158 0 : }
2159 :
2160 : std::string
2161 19 : MultiAppGeneralFieldTransfer::getDataSourceName(unsigned int var_index) const
2162 : {
2163 : mooseAssert(var_index < _from_var_names.size(), "No source variable at this index");
2164 19 : return "variable '" + getFromVarName(var_index) + "'";
2165 : }
2166 :
2167 : VariableName
2168 54851 : MultiAppGeneralFieldTransfer::getFromVarName(unsigned int var_index) const
2169 : {
2170 : mooseAssert(var_index < _from_var_names.size(), "No source variable at this index");
2171 54851 : VariableName var_name = _from_var_names[var_index];
2172 54851 : if (_from_var_components.size())
2173 280 : var_name += "_" + std::to_string(_from_var_components[var_index]);
2174 54851 : return var_name;
2175 0 : }
2176 :
2177 : VariableName
2178 157928 : MultiAppGeneralFieldTransfer::getToVarName(unsigned int var_index)
2179 : {
2180 : mooseAssert(var_index < _to_var_names.size(), "No target variable at this index");
2181 157928 : VariableName var_name = _to_var_names[var_index];
2182 157928 : if (_to_var_components.size())
2183 528 : var_name += "_" + std::to_string(_to_var_components[var_index]);
2184 157928 : return var_name;
2185 0 : }
2186 :
2187 : Point
2188 104940 : MultiAppGeneralFieldTransfer::getMaxToProblemsBBoxDimensions() const
2189 : {
2190 : Point max_dimension = {std::numeric_limits<Real>::min(),
2191 : std::numeric_limits<Real>::min(),
2192 104940 : std::numeric_limits<Real>::min()};
2193 :
2194 214584 : for (const auto & to_mesh : _to_meshes)
2195 : {
2196 109644 : const auto bbox = to_mesh->getInflatedProcessorBoundingBox();
2197 438576 : for (const auto dim : make_range(LIBMESH_DIM))
2198 328932 : max_dimension(dim) = std::max(
2199 657864 : max_dimension(dim), std::max(std::abs(bbox.first(dim)), std::abs(bbox.second(dim))));
2200 : }
2201 :
2202 104940 : return max_dimension;
2203 : }
2204 :
2205 : bool
2206 2274505 : MultiAppGeneralFieldTransfer::detectConflict(Real current_value,
2207 : Real new_value,
2208 : Real current_distance,
2209 : Real new_distance) const
2210 : {
2211 : // No conflict if we're not looking for them
2212 2274505 : if (_search_value_conflicts)
2213 : // Only consider conflicts if the values are valid and different
2214 4545 : if (current_value != GeneralFieldTransfer::OutOfMeshValue &&
2215 7066 : new_value != GeneralFieldTransfer::OutOfMeshValue &&
2216 2521 : !MooseUtils::absoluteFuzzyEqual(current_value, new_value))
2217 : // Conflict only occurs if the origin points are equidistant
2218 2269 : if (MooseUtils::absoluteFuzzyEqual(current_distance, new_distance))
2219 57 : return true;
2220 2274448 : return false;
2221 : }
|
