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 "MultiAppShapeEvaluationTransfer.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 "MooseAppCoordTransform.h"
19 :
20 : #include "libmesh/meshfree_interpolation.h"
21 : #include "libmesh/system.h"
22 : #include "libmesh/mesh_function.h"
23 : #include "libmesh/mesh_tools.h"
24 : #include "libmesh/parallel_algebra.h" // for communicator send and receive stuff
25 :
26 : // TIMPI includes
27 : #include "timpi/communicator.h"
28 : #include "timpi/parallel_sync.h"
29 :
30 : using namespace libMesh;
31 :
32 : registerMooseObjectDeprecated("MooseApp", MultiAppShapeEvaluationTransfer, "12/31/2024 24:00");
33 : registerMooseObjectRenamed("MooseApp",
34 : MultiAppMeshFunctionTransfer,
35 : "12/31/2023 24:00",
36 : MultiAppShapeEvaluationTransfer);
37 :
38 : InputParameters
39 29052 : MultiAppShapeEvaluationTransfer::validParams()
40 : {
41 29052 : InputParameters params = MultiAppConservativeTransfer::validParams();
42 29052 : params.addClassDescription(
43 : "Transfers field data at the MultiApp position using solution the finite element function "
44 : "from the main/parent application, via a 'libMesh::MeshFunction' object.");
45 :
46 87156 : params.addParam<bool>(
47 : "error_on_miss",
48 58104 : false,
49 : "Whether or not to error in the case that a target point is not found in the source domain.");
50 29052 : MultiAppTransfer::addBBoxFactorParam(params);
51 29052 : return params;
52 0 : }
53 :
54 261 : MultiAppShapeEvaluationTransfer::MultiAppShapeEvaluationTransfer(const InputParameters & parameters)
55 261 : : MultiAppConservativeTransfer(parameters), _error_on_miss(getParam<bool>("error_on_miss"))
56 : {
57 261 : mooseDeprecated("MultiAppShapeEvaluationTransfer is deprecated. Use "
58 : "MultiAppGeneralFieldShapeEvaluationTransfer instead and adapt the parameters");
59 :
60 261 : if (_to_var_names.size() == _from_var_names.size())
61 261 : _var_size = _to_var_names.size();
62 : else
63 0 : paramError("variable", "The number of variables to transfer to and from should be equal");
64 261 : }
65 :
66 : void
67 675 : MultiAppShapeEvaluationTransfer::execute()
68 : {
69 675 : TIME_SECTION("MultiAppShapeEvaluationTransfer::execute()",
70 : 5,
71 : "Transferring variables via finite element interpolation");
72 :
73 : // loop over the vector of variables and make the transfer one by one
74 1346 : for (unsigned int i = 0; i < _var_size; ++i)
75 675 : transferVariable(i);
76 :
77 671 : postExecute();
78 671 : }
79 :
80 : void
81 675 : MultiAppShapeEvaluationTransfer::transferVariable(unsigned int i)
82 : {
83 : mooseAssert(i < _var_size, "The variable of index " << i << " does not exist");
84 :
85 : /**
86 : * For every combination of global "from" problem and local "to" problem, find
87 : * which "from" bounding boxes overlap with which "to" elements. Keep track
88 : * of which processors own bounding boxes that overlap with which elements.
89 : * Build vectors of node locations/element centroids to send to other
90 : * processors for mesh function evaluations.
91 : */
92 :
93 : // Get the bounding boxes for the "from" domains.
94 675 : std::vector<BoundingBox> bboxes = getFromBoundingBoxes();
95 :
96 : // Figure out how many "from" domains each processor owns.
97 675 : std::vector<unsigned int> froms_per_proc = getFromsPerProc();
98 :
99 : // Point locations needed to send to from-domain
100 : // processor to points
101 675 : std::map<processor_id_type, std::vector<Point>> outgoing_points;
102 : // <processor, <system_id, node_i>> --> point_id
103 : std::map<processor_id_type, std::map<std::pair<unsigned int, dof_id_type>, dof_id_type>>
104 675 : point_index_map;
105 :
106 1350 : for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
107 : {
108 675 : System * to_sys = find_sys(*_to_es[i_to], _to_var_names[i]);
109 675 : unsigned int sys_num = to_sys->number();
110 675 : unsigned int var_num = to_sys->variable_number(_to_var_names[i]);
111 675 : MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
112 675 : auto & fe_type = to_sys->variable_type(var_num);
113 675 : bool is_constant = fe_type.order == CONSTANT;
114 675 : bool is_nodal = fe_type.family == LAGRANGE;
115 675 : const auto to_global_num = _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
116 675 : const auto & to_transform = *_to_transforms[to_global_num];
117 :
118 675 : if (fe_type.order > FIRST && !is_nodal)
119 0 : mooseError("We don't currently support second order or higher elemental variable.");
120 :
121 675 : if (is_nodal)
122 : {
123 85633 : for (const auto & node : to_mesh->local_node_ptr_range())
124 : {
125 : // Skip this node if the variable has no dofs at it.
126 42516 : if (node->n_dofs(sys_num, var_num) < 1)
127 0 : continue;
128 :
129 : // Loop over the "froms" on processor i_proc. If the node is found in
130 : // any of the "froms", add that node to the vector that will be sent to
131 : // i_proc.
132 42516 : unsigned int from0 = 0;
133 101157 : for (processor_id_type i_proc = 0; i_proc < n_processors();
134 58641 : from0 += froms_per_proc[i_proc], ++i_proc)
135 : {
136 58641 : bool point_found = false;
137 117282 : for (unsigned int i_from = from0; i_from < from0 + froms_per_proc[i_proc] && !point_found;
138 : ++i_from)
139 : {
140 58641 : auto transformed_node = to_transform(*node);
141 58641 : if (bboxes[i_from].contains_point(transformed_node))
142 : {
143 : // <system id, node id>
144 44425 : std::pair<unsigned int, dof_id_type> key(i_to, node->id());
145 : // map a tuple of pid, problem id and node id to point id
146 : // point id is counted from zero
147 44425 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
148 : // map pid to points
149 44425 : outgoing_points[i_proc].push_back(std::move(transformed_node));
150 44425 : point_found = true;
151 : }
152 : }
153 : }
154 601 : }
155 : }
156 : else // Elemental
157 : {
158 74 : std::vector<Point> points;
159 74 : std::vector<dof_id_type> point_ids;
160 35274 : for (auto & elem : as_range(to_mesh->local_elements_begin(), to_mesh->local_elements_end()))
161 : {
162 : // Skip this element if the variable has no dofs at it.
163 17600 : if (elem->n_dofs(sys_num, var_num) < 1)
164 0 : continue;
165 :
166 17600 : points.clear();
167 17600 : point_ids.clear();
168 : // grab sample points
169 : // for constant shape function, we take the element centroid
170 17600 : if (is_constant)
171 : {
172 1600 : points.push_back(elem->vertex_average());
173 1600 : point_ids.push_back(elem->id());
174 : }
175 :
176 : // for higher order method, we take all nodes of element
177 : // this works for the first order L2 Lagrange.
178 : else
179 80000 : for (auto & node : elem->node_ref_range())
180 : {
181 64000 : points.push_back(node);
182 64000 : point_ids.push_back(node.id());
183 : }
184 :
185 17600 : unsigned int offset = 0;
186 83200 : for (auto & point : points)
187 : {
188 : // Loop over the "froms" on processor i_proc. If the elem is found in
189 : // any of the "froms", add that elem to the vector that will be sent to
190 : // i_proc.
191 65600 : unsigned int from0 = 0;
192 151000 : for (processor_id_type i_proc = 0; i_proc < n_processors();
193 85400 : from0 += froms_per_proc[i_proc], ++i_proc)
194 : {
195 85400 : bool point_found = false;
196 170800 : for (unsigned int i_from = from0;
197 170800 : i_from < from0 + froms_per_proc[i_proc] && !point_found;
198 : ++i_from)
199 : {
200 85400 : auto transformed_point = to_transform(point);
201 85400 : if (bboxes[i_from].contains_point(transformed_point))
202 : {
203 66378 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
204 66378 : if (point_index_map[i_proc].find(key) != point_index_map[i_proc].end())
205 46540 : continue;
206 :
207 19838 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
208 19838 : outgoing_points[i_proc].push_back(std::move(transformed_point));
209 19838 : point_found = true;
210 : } // if
211 : } // i_from
212 : } // i_proc
213 65600 : offset++;
214 : } // point
215 :
216 74 : } // else
217 74 : }
218 : }
219 :
220 : // Get the local bounding boxes for current processor.
221 : // There could be more than one box because of the number of local apps
222 : // can be larger than one
223 675 : std::vector<BoundingBox> local_bboxes(froms_per_proc[processor_id()]);
224 : {
225 : // Find the index to the first of this processor's local bounding boxes.
226 675 : unsigned int local_start = 0;
227 858 : for (processor_id_type i_proc = 0; i_proc < n_processors() && i_proc != processor_id();
228 : ++i_proc)
229 : {
230 183 : local_start += froms_per_proc[i_proc];
231 : }
232 :
233 : // Extract the local bounding boxes.
234 1350 : for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; ++i_from)
235 : {
236 675 : local_bboxes[i_from] = bboxes[local_start + i_from];
237 : }
238 : }
239 :
240 : // Setup the local mesh functions.
241 675 : std::vector<MeshFunction> local_meshfuns;
242 675 : local_meshfuns.reserve(_from_problems.size());
243 1350 : for (unsigned int i_from = 0; i_from < _from_problems.size(); ++i_from)
244 : {
245 675 : FEProblemBase & from_problem = *_from_problems[i_from];
246 : MooseVariableFEBase & from_var =
247 675 : from_problem.getVariable(0,
248 675 : _from_var_names[i],
249 : Moose::VarKindType::VAR_ANY,
250 : Moose::VarFieldType::VAR_FIELD_STANDARD);
251 675 : System & from_sys = from_var.sys().system();
252 675 : unsigned int from_var_num = from_sys.variable_number(from_var.name());
253 :
254 1350 : local_meshfuns.emplace_back(getEquationSystem(from_problem, _displaced_source_mesh),
255 675 : *from_sys.current_local_solution,
256 : from_sys.get_dof_map(),
257 : from_var_num);
258 675 : local_meshfuns.back().init();
259 675 : local_meshfuns.back().enable_out_of_mesh_mode(OutOfMeshValue);
260 : }
261 :
262 : /**
263 : * Gather all of the evaluations, pick out the best ones for each point, and
264 : * apply them to the solution vector. When we are transferring from
265 : * multiapps, there may be multiple overlapping apps for a particular point.
266 : * In that case, we'll try to use the value from the app with the lowest id.
267 : */
268 :
269 : // Fill values and app ids for incoming points
270 : // We are responsible to compute values for these incoming points
271 : auto gather_functor =
272 915 : [this, &local_meshfuns, &local_bboxes](
273 : processor_id_type /*pid*/,
274 : const std::vector<Point> & incoming_points,
275 534535 : std::vector<std::pair<Real, unsigned int>> & vals_ids_for_incoming_points)
276 : {
277 915 : vals_ids_for_incoming_points.resize(incoming_points.size(), std::make_pair(OutOfMeshValue, 0));
278 65178 : for (MooseIndex(incoming_points.size()) i_pt = 0; i_pt < incoming_points.size(); ++i_pt)
279 : {
280 64263 : Point pt = incoming_points[i_pt];
281 :
282 : // Loop until we've found the lowest-ranked app that actually contains
283 : // the quadrature point.
284 128526 : for (MooseIndex(_from_problems.size()) i_from = 0;
285 192789 : i_from < _from_problems.size() &&
286 64263 : vals_ids_for_incoming_points[i_pt].first == OutOfMeshValue;
287 : ++i_from)
288 : {
289 64263 : if (local_bboxes[i_from].contains_point(pt))
290 : {
291 : const auto from_global_num =
292 64263 : _current_direction == TO_MULTIAPP ? 0 : _from_local2global_map[i_from];
293 : // Use mesh function to compute interpolation values
294 128526 : vals_ids_for_incoming_points[i_pt].first =
295 64263 : (local_meshfuns[i_from])(_from_transforms[from_global_num]->mapBack(pt));
296 : // Record problem ID as well
297 64263 : switch (_current_direction)
298 : {
299 20431 : case FROM_MULTIAPP:
300 20431 : vals_ids_for_incoming_points[i_pt].second = _from_local2global_map[i_from];
301 20431 : break;
302 43832 : case TO_MULTIAPP:
303 43832 : vals_ids_for_incoming_points[i_pt].second = _to_local2global_map[i_from];
304 43832 : break;
305 0 : default:
306 0 : mooseError("Unsupported direction");
307 : }
308 : }
309 : }
310 : }
311 915 : };
312 :
313 : // Incoming values and APP ids for outgoing points
314 675 : std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> incoming_vals_ids;
315 : // Copy data out to incoming_vals_ids
316 : auto action_functor =
317 915 : [&incoming_vals_ids](
318 : processor_id_type pid,
319 : const std::vector<Point> & /*my_outgoing_points*/,
320 1830 : const std::vector<std::pair<Real, unsigned int>> & vals_ids_for_outgoing_points)
321 : {
322 : // This lambda function might be called multiple times
323 915 : incoming_vals_ids[pid].reserve(vals_ids_for_outgoing_points.size());
324 : // Copy data for processor 'pid'
325 915 : std::copy(vals_ids_for_outgoing_points.begin(),
326 : vals_ids_for_outgoing_points.end(),
327 915 : std::back_inserter(incoming_vals_ids[pid]));
328 1590 : };
329 :
330 : // We assume incoming_vals_ids is ordered in the same way as outgoing_points
331 : // Hopefully, pull_parallel_vector_data will not mess up this
332 675 : const std::pair<Real, unsigned int> * ex = nullptr;
333 675 : libMesh::Parallel::pull_parallel_vector_data(
334 675 : comm(), outgoing_points, gather_functor, action_functor, ex);
335 :
336 1346 : for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
337 : {
338 675 : const auto to_global_num = _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
339 675 : System * to_sys = find_sys(*_to_es[i_to], _to_var_names[i]);
340 :
341 675 : unsigned int sys_num = to_sys->number();
342 675 : unsigned int var_num = to_sys->variable_number(_to_var_names[i]);
343 :
344 675 : NumericVector<Real> * solution = nullptr;
345 675 : switch (_current_direction)
346 : {
347 371 : case TO_MULTIAPP:
348 371 : solution = &getTransferVector(i_to, _to_var_names[i]);
349 371 : break;
350 304 : case FROM_MULTIAPP:
351 304 : solution = to_sys->solution.get();
352 304 : break;
353 0 : default:
354 0 : mooseError("Unknown direction");
355 : }
356 :
357 675 : MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
358 675 : auto & fe_type = to_sys->variable_type(var_num);
359 675 : bool is_constant = fe_type.order == CONSTANT;
360 675 : bool is_nodal = fe_type.family == LAGRANGE;
361 :
362 675 : if (is_nodal)
363 : {
364 84761 : for (const auto & node : to_mesh->local_node_ptr_range())
365 : {
366 : // Skip this node if the variable has no dofs at it.
367 42084 : if (node->n_dofs(sys_num, var_num) < 1)
368 0 : continue;
369 :
370 42084 : unsigned int lowest_app_rank = libMesh::invalid_uint;
371 42084 : Real best_val = 0.;
372 42084 : bool point_found = false;
373 95627 : for (auto & group : incoming_vals_ids)
374 : {
375 : // Skip this proc if the node wasn't in it's bounding boxes.
376 53543 : std::pair<unsigned int, dof_id_type> key(i_to, node->id());
377 : // Make sure point_index_map has data for corresponding pid
378 : mooseAssert(point_index_map.find(group.first) != point_index_map.end(),
379 : "Point index map does not have data for processor group.first");
380 53543 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
381 9972 : continue;
382 :
383 44209 : auto i_pt = point_index_map[group.first][key];
384 :
385 : // Ignore this proc if it's app has a higher rank than the
386 : // previously found lowest app rank.
387 44209 : if (_current_direction == FROM_MULTIAPP)
388 : {
389 19591 : if (group.second[i_pt].second >= lowest_app_rank)
390 0 : continue;
391 : }
392 :
393 : // Ignore this proc if the point was actually outside its meshes.
394 44209 : if (group.second[i_pt].first == OutOfMeshValue)
395 638 : continue;
396 :
397 43571 : best_val = group.second[i_pt].first;
398 43571 : point_found = true;
399 : }
400 :
401 42084 : if (_error_on_miss && !point_found)
402 4 : mooseError("Point not found in the reference space! ",
403 4 : (*_to_transforms[to_global_num])(*node));
404 :
405 42080 : dof_id_type dof = node->dof_number(sys_num, var_num, 0);
406 42080 : solution->set(dof, best_val);
407 597 : }
408 : }
409 : else // Elemental
410 : {
411 74 : std::vector<Point> points;
412 74 : std::vector<dof_id_type> point_ids;
413 35274 : for (auto & elem : as_range(to_mesh->local_elements_begin(), to_mesh->local_elements_end()))
414 : {
415 : // Skip this element if the variable has no dofs at it.
416 17600 : if (elem->n_dofs(sys_num, var_num) < 1)
417 0 : continue;
418 :
419 17600 : points.clear();
420 17600 : point_ids.clear();
421 : // grab sample points
422 : // for constant shape function, we take the element centroid
423 17600 : if (is_constant)
424 : {
425 1600 : points.push_back(elem->vertex_average());
426 1600 : point_ids.push_back(elem->id());
427 : }
428 : // for higher order method, we take all nodes of element
429 : // this works for the first order L2 Lagrange. Might not work
430 : // with something higher than the first order
431 : else
432 : {
433 80000 : for (auto & node : elem->node_ref_range())
434 : {
435 64000 : points.push_back(node);
436 64000 : point_ids.push_back(node.id());
437 : }
438 : }
439 :
440 17600 : auto n_points = points.size();
441 17600 : unsigned int n_comp = elem->n_comp(sys_num, var_num);
442 : // We assume each point corresponds to one component of elemental variable
443 17600 : if (n_points != n_comp)
444 0 : mooseError(" Number of points ",
445 : n_points,
446 : " does not equal to number of variable components ",
447 : n_comp);
448 83200 : for (unsigned int offset = 0; offset < n_points; offset++)
449 : {
450 65600 : unsigned int lowest_app_rank = libMesh::invalid_uint;
451 65600 : Real best_val = 0;
452 65600 : bool point_found = false;
453 147600 : for (auto & group : incoming_vals_ids)
454 : {
455 : // Skip this proc if the elem wasn't in it's bounding boxes.
456 82000 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
457 82000 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
458 17322 : continue;
459 :
460 66378 : unsigned int i_pt = point_index_map[group.first][key];
461 :
462 : // Ignore this proc if it's app has a higher rank than the
463 : // previously found lowest app rank.
464 66378 : if (_current_direction == FROM_MULTIAPP)
465 : {
466 840 : if (group.second[i_pt].second >= lowest_app_rank)
467 0 : continue;
468 : }
469 :
470 : // Ignore this proc if the point was actually outside its meshes.
471 66378 : if (group.second[i_pt].first == OutOfMeshValue)
472 1700 : continue;
473 :
474 64678 : best_val = group.second[i_pt].first;
475 64678 : point_found = true;
476 : }
477 :
478 65600 : if (_error_on_miss && !point_found)
479 0 : mooseError("Point not found in the reference space! ",
480 0 : (*_to_transforms[to_global_num])(elem->vertex_average()));
481 :
482 : // Get the value for a dof
483 65600 : dof_id_type dof = elem->dof_number(sys_num, var_num, offset);
484 65600 : solution->set(dof, best_val);
485 : } // point
486 74 : } // element
487 74 : }
488 671 : solution->close();
489 671 : to_sys->update();
490 : }
491 671 : }
|
