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 6636 : MultiAppShapeEvaluationTransfer::validParams()
40 : {
41 6636 : InputParameters params = MultiAppConservativeTransfer::validParams();
42 13272 : 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 13272 : params.addParam<bool>(
47 : "error_on_miss",
48 13272 : false,
49 : "Whether or not to error in the case that a target point is not found in the source domain.");
50 6636 : MultiAppTransfer::addBBoxFactorParam(params);
51 6636 : return params;
52 0 : }
53 :
54 257 : MultiAppShapeEvaluationTransfer::MultiAppShapeEvaluationTransfer(const InputParameters & parameters)
55 514 : : MultiAppConservativeTransfer(parameters), _error_on_miss(getParam<bool>("error_on_miss"))
56 : {
57 257 : mooseDeprecated("MultiAppShapeEvaluationTransfer is deprecated. Use "
58 : "MultiAppGeneralFieldShapeEvaluationTransfer instead and adapt the parameters");
59 :
60 257 : if (_to_var_names.size() == _from_var_names.size())
61 257 : _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 257 : }
65 :
66 : void
67 675 : MultiAppShapeEvaluationTransfer::execute()
68 : {
69 3375 : 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 1347 : for (unsigned int i = 0; i < _var_size; ++i)
75 675 : transferVariable(i);
76 :
77 672 : postExecute();
78 672 : }
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 43361 : for (const auto & node : to_mesh->local_node_ptr_range())
124 : {
125 : // Skip this node if the variable has no dofs at it.
126 42758 : 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 42758 : unsigned int from0 = 0;
133 101641 : for (processor_id_type i_proc = 0; i_proc < n_processors();
134 58883 : from0 += froms_per_proc[i_proc], ++i_proc)
135 : {
136 58883 : bool point_found = false;
137 117766 : for (unsigned int i_from = from0; i_from < from0 + froms_per_proc[i_proc] && !point_found;
138 : ++i_from)
139 : {
140 58883 : auto transformed_node = to_transform(*node);
141 58883 : if (bboxes[i_from].contains_point(transformed_node))
142 : {
143 : // <system id, node id>
144 44722 : 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 44722 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
148 : // map pid to points
149 44722 : outgoing_points[i_proc].push_back(std::move(transformed_node));
150 44722 : point_found = true;
151 : }
152 : }
153 : }
154 603 : }
155 : }
156 : else // Elemental
157 : {
158 72 : std::vector<Point> points;
159 72 : std::vector<dof_id_type> point_ids;
160 16872 : 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 16800 : if (elem->n_dofs(sys_num, var_num) < 1)
164 0 : continue;
165 :
166 16800 : points.clear();
167 16800 : point_ids.clear();
168 : // grab sample points
169 : // for constant shape function, we take the element centroid
170 16800 : 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 76000 : for (auto & node : elem->node_ref_range())
180 : {
181 60800 : points.push_back(node);
182 60800 : point_ids.push_back(node.id());
183 : }
184 :
185 16800 : unsigned int offset = 0;
186 79200 : 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 62400 : unsigned int from0 = 0;
192 144600 : for (processor_id_type i_proc = 0; i_proc < n_processors();
193 82200 : from0 += froms_per_proc[i_proc], ++i_proc)
194 : {
195 82200 : bool point_found = false;
196 164400 : for (unsigned int i_from = from0;
197 164400 : i_from < from0 + froms_per_proc[i_proc] && !point_found;
198 : ++i_from)
199 : {
200 82200 : auto transformed_point = to_transform(point);
201 82200 : if (bboxes[i_from].contains_point(transformed_point))
202 : {
203 63178 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
204 63178 : if (point_index_map[i_proc].find(key) != point_index_map[i_proc].end())
205 44222 : continue;
206 :
207 18956 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
208 18956 : outgoing_points[i_proc].push_back(std::move(transformed_point));
209 18956 : point_found = true;
210 : } // if
211 : } // i_from
212 : } // i_proc
213 62400 : offset++;
214 : } // point
215 :
216 72 : } // else
217 72 : }
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 : 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 64593 : for (MooseIndex(incoming_points.size()) i_pt = 0; i_pt < incoming_points.size(); ++i_pt)
279 : {
280 63678 : 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 127356 : for (MooseIndex(_from_problems.size()) i_from = 0;
285 191034 : i_from < _from_problems.size() &&
286 63678 : vals_ids_for_incoming_points[i_pt].first == OutOfMeshValue;
287 : ++i_from)
288 : {
289 63678 : if (local_bboxes[i_from].contains_point(pt))
290 : {
291 : // Use mesh function to compute interpolation values
292 127356 : vals_ids_for_incoming_points[i_pt].first = (local_meshfuns[i_from])(
293 191034 : getPointInSourceAppFrame(pt, i_from, "Shape evaluation transfer"));
294 : // Record problem ID as well
295 63678 : switch (_current_direction)
296 : {
297 20552 : case FROM_MULTIAPP:
298 20552 : vals_ids_for_incoming_points[i_pt].second = _from_local2global_map[i_from];
299 20552 : break;
300 43126 : case TO_MULTIAPP:
301 43126 : vals_ids_for_incoming_points[i_pt].second = _to_local2global_map[i_from];
302 43126 : break;
303 0 : default:
304 0 : mooseError("Unsupported direction");
305 : }
306 : }
307 : }
308 : }
309 915 : };
310 :
311 : // Incoming values and APP ids for outgoing points
312 675 : std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> incoming_vals_ids;
313 : // Copy data out to incoming_vals_ids
314 : auto action_functor =
315 915 : [&incoming_vals_ids](
316 : processor_id_type pid,
317 : const std::vector<Point> & /*my_outgoing_points*/,
318 : const std::vector<std::pair<Real, unsigned int>> & vals_ids_for_outgoing_points)
319 : {
320 : // This lambda function might be called multiple times
321 915 : incoming_vals_ids[pid].reserve(vals_ids_for_outgoing_points.size());
322 : // Copy data for processor 'pid'
323 915 : std::copy(vals_ids_for_outgoing_points.begin(),
324 : vals_ids_for_outgoing_points.end(),
325 915 : std::back_inserter(incoming_vals_ids[pid]));
326 1590 : };
327 :
328 : // We assume incoming_vals_ids is ordered in the same way as outgoing_points
329 : // Hopefully, pull_parallel_vector_data will not mess up this
330 675 : const std::pair<Real, unsigned int> * ex = nullptr;
331 675 : libMesh::Parallel::pull_parallel_vector_data(
332 675 : comm(), outgoing_points, gather_functor, action_functor, ex);
333 :
334 1347 : for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
335 : {
336 675 : const auto to_global_num = _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
337 675 : System * to_sys = find_sys(*_to_es[i_to], _to_var_names[i]);
338 :
339 675 : unsigned int sys_num = to_sys->number();
340 675 : unsigned int var_num = to_sys->variable_number(_to_var_names[i]);
341 :
342 675 : NumericVector<Real> * solution = nullptr;
343 675 : switch (_current_direction)
344 : {
345 370 : case TO_MULTIAPP:
346 370 : solution = &getTransferVector(i_to, _to_var_names[i]);
347 370 : break;
348 305 : case FROM_MULTIAPP:
349 305 : solution = to_sys->solution.get();
350 305 : break;
351 0 : default:
352 0 : mooseError("Unknown direction");
353 : }
354 :
355 675 : MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
356 675 : auto & fe_type = to_sys->variable_type(var_num);
357 675 : bool is_constant = fe_type.order == CONSTANT;
358 675 : bool is_nodal = fe_type.family == LAGRANGE;
359 :
360 675 : if (is_nodal)
361 : {
362 43034 : for (const auto & node : to_mesh->local_node_ptr_range())
363 : {
364 : // Skip this node if the variable has no dofs at it.
365 42434 : if (node->n_dofs(sys_num, var_num) < 1)
366 0 : continue;
367 :
368 42434 : unsigned int lowest_app_rank = libMesh::invalid_uint;
369 42434 : Real best_val = 0.;
370 42434 : bool point_found = false;
371 96327 : for (auto & group : incoming_vals_ids)
372 : {
373 : // Skip this proc if the node wasn't in it's bounding boxes.
374 53893 : std::pair<unsigned int, dof_id_type> key(i_to, node->id());
375 : // Make sure point_index_map has data for corresponding pid
376 : mooseAssert(point_index_map.find(group.first) != point_index_map.end(),
377 : "Point index map does not have data for processor group.first");
378 53893 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
379 9975 : continue;
380 :
381 44560 : auto i_pt = point_index_map[group.first][key];
382 :
383 : // Ignore this proc if it's app has a higher rank than the
384 : // previously found lowest app rank.
385 44560 : if (_current_direction == FROM_MULTIAPP)
386 : {
387 19712 : if (group.second[i_pt].second >= lowest_app_rank)
388 0 : continue;
389 : }
390 :
391 : // Ignore this proc if the point was actually outside its meshes.
392 44560 : if (group.second[i_pt].first == OutOfMeshValue)
393 642 : continue;
394 :
395 43918 : best_val = group.second[i_pt].first;
396 43918 : point_found = true;
397 : }
398 :
399 42434 : if (_error_on_miss && !point_found)
400 3 : mooseError("Point not found in the reference space! ",
401 3 : (*_to_transforms[to_global_num])(*node));
402 :
403 42431 : dof_id_type dof = node->dof_number(sys_num, var_num, 0);
404 42431 : solution->set(dof, best_val);
405 600 : }
406 : }
407 : else // Elemental
408 : {
409 72 : std::vector<Point> points;
410 72 : std::vector<dof_id_type> point_ids;
411 16872 : for (auto & elem : as_range(to_mesh->local_elements_begin(), to_mesh->local_elements_end()))
412 : {
413 : // Skip this element if the variable has no dofs at it.
414 16800 : if (elem->n_dofs(sys_num, var_num) < 1)
415 0 : continue;
416 :
417 16800 : points.clear();
418 16800 : point_ids.clear();
419 : // grab sample points
420 : // for constant shape function, we take the element centroid
421 16800 : if (is_constant)
422 : {
423 1600 : points.push_back(elem->vertex_average());
424 1600 : point_ids.push_back(elem->id());
425 : }
426 : // for higher order method, we take all nodes of element
427 : // this works for the first order L2 Lagrange. Might not work
428 : // with something higher than the first order
429 : else
430 : {
431 76000 : for (auto & node : elem->node_ref_range())
432 : {
433 60800 : points.push_back(node);
434 60800 : point_ids.push_back(node.id());
435 : }
436 : }
437 :
438 16800 : auto n_points = points.size();
439 16800 : unsigned int n_comp = elem->n_comp(sys_num, var_num);
440 : // We assume each point corresponds to one component of elemental variable
441 16800 : if (n_points != n_comp)
442 0 : mooseError(" Number of points ",
443 : n_points,
444 : " does not equal to number of variable components ",
445 : n_comp);
446 79200 : for (unsigned int offset = 0; offset < n_points; offset++)
447 : {
448 62400 : unsigned int lowest_app_rank = libMesh::invalid_uint;
449 62400 : Real best_val = 0;
450 62400 : bool point_found = false;
451 141200 : for (auto & group : incoming_vals_ids)
452 : {
453 : // Skip this proc if the elem wasn't in it's bounding boxes.
454 78800 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
455 78800 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
456 17336 : continue;
457 :
458 63178 : unsigned int i_pt = point_index_map[group.first][key];
459 :
460 : // Ignore this proc if it's app has a higher rank than the
461 : // previously found lowest app rank.
462 63178 : if (_current_direction == FROM_MULTIAPP)
463 : {
464 840 : if (group.second[i_pt].second >= lowest_app_rank)
465 0 : continue;
466 : }
467 :
468 : // Ignore this proc if the point was actually outside its meshes.
469 63178 : if (group.second[i_pt].first == OutOfMeshValue)
470 1714 : continue;
471 :
472 61464 : best_val = group.second[i_pt].first;
473 61464 : point_found = true;
474 : }
475 :
476 62400 : if (_error_on_miss && !point_found)
477 0 : mooseError("Point not found in the reference space! ",
478 0 : (*_to_transforms[to_global_num])(elem->vertex_average()));
479 :
480 : // Get the value for a dof
481 62400 : dof_id_type dof = elem->dof_number(sys_num, var_num, offset);
482 62400 : solution->set(dof, best_val);
483 : } // point
484 72 : } // element
485 72 : }
486 672 : solution->close();
487 672 : to_sys->update();
488 : }
489 672 : }
|
