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 29096 : MultiAppShapeEvaluationTransfer::validParams()
40 : {
41 29096 : InputParameters params = MultiAppConservativeTransfer::validParams();
42 29096 : 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 87288 : params.addParam<bool>(
47 : "error_on_miss",
48 58192 : false,
49 : "Whether or not to error in the case that a target point is not found in the source domain.");
50 29096 : MultiAppTransfer::addBBoxFactorParam(params);
51 29096 : return params;
52 0 : }
53 :
54 283 : MultiAppShapeEvaluationTransfer::MultiAppShapeEvaluationTransfer(const InputParameters & parameters)
55 283 : : MultiAppConservativeTransfer(parameters), _error_on_miss(getParam<bool>("error_on_miss"))
56 : {
57 283 : mooseDeprecated("MultiAppShapeEvaluationTransfer is deprecated. Use "
58 : "MultiAppGeneralFieldShapeEvaluationTransfer instead and adapt the parameters");
59 :
60 283 : if (_to_var_names.size() == _from_var_names.size())
61 283 : _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 283 : }
65 :
66 : void
67 740 : MultiAppShapeEvaluationTransfer::execute()
68 : {
69 740 : 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 1476 : for (unsigned int i = 0; i < _var_size; ++i)
75 740 : transferVariable(i);
76 :
77 736 : postExecute();
78 736 : }
79 :
80 : void
81 740 : 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 740 : std::vector<BoundingBox> bboxes = getFromBoundingBoxes();
95 :
96 : // Figure out how many "from" domains each processor owns.
97 740 : std::vector<unsigned int> froms_per_proc = getFromsPerProc();
98 :
99 : // Point locations needed to send to from-domain
100 : // processor to points
101 740 : 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 740 : point_index_map;
105 :
106 1480 : for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
107 : {
108 740 : System * to_sys = find_sys(*_to_es[i_to], _to_var_names[i]);
109 740 : unsigned int sys_num = to_sys->number();
110 740 : unsigned int var_num = to_sys->variable_number(_to_var_names[i]);
111 740 : MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
112 740 : auto & fe_type = to_sys->variable_type(var_num);
113 740 : bool is_constant = fe_type.order == CONSTANT;
114 740 : bool is_nodal = fe_type.family == LAGRANGE;
115 740 : const auto to_global_num = _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
116 740 : const auto & to_transform = *_to_transforms[to_global_num];
117 :
118 740 : if (fe_type.order > FIRST && !is_nodal)
119 0 : mooseError("We don't currently support second order or higher elemental variable.");
120 :
121 740 : if (is_nodal)
122 : {
123 97167 : for (const auto & node : to_mesh->local_node_ptr_range())
124 : {
125 : // Skip this node if the variable has no dofs at it.
126 48254 : 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 48254 : unsigned int from0 = 0;
133 112633 : for (processor_id_type i_proc = 0; i_proc < n_processors();
134 64379 : from0 += froms_per_proc[i_proc], ++i_proc)
135 : {
136 64379 : bool point_found = false;
137 128758 : for (unsigned int i_from = from0; i_from < from0 + froms_per_proc[i_proc] && !point_found;
138 : ++i_from)
139 : {
140 64379 : auto transformed_node = to_transform(*node);
141 64379 : if (bboxes[i_from].contains_point(transformed_node))
142 : {
143 : // <system id, node id>
144 50163 : 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 50163 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
148 : // map pid to points
149 50163 : outgoing_points[i_proc].push_back(std::move(transformed_node));
150 50163 : point_found = true;
151 : }
152 : }
153 : }
154 659 : }
155 : }
156 : else // Elemental
157 : {
158 81 : std::vector<Point> points;
159 81 : std::vector<dof_id_type> point_ids;
160 39681 : 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 19800 : if (elem->n_dofs(sys_num, var_num) < 1)
164 0 : continue;
165 :
166 19800 : points.clear();
167 19800 : point_ids.clear();
168 : // grab sample points
169 : // for constant shape function, we take the element centroid
170 19800 : if (is_constant)
171 : {
172 1800 : points.push_back(elem->vertex_average());
173 1800 : 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 90000 : for (auto & node : elem->node_ref_range())
180 : {
181 72000 : points.push_back(node);
182 72000 : point_ids.push_back(node.id());
183 : }
184 :
185 19800 : unsigned int offset = 0;
186 93600 : 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 73800 : unsigned int from0 = 0;
192 167400 : for (processor_id_type i_proc = 0; i_proc < n_processors();
193 93600 : from0 += froms_per_proc[i_proc], ++i_proc)
194 : {
195 93600 : bool point_found = false;
196 187200 : for (unsigned int i_from = from0;
197 187200 : i_from < from0 + froms_per_proc[i_proc] && !point_found;
198 : ++i_from)
199 : {
200 93600 : auto transformed_point = to_transform(point);
201 93600 : if (bboxes[i_from].contains_point(transformed_point))
202 : {
203 74347 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
204 74347 : if (point_index_map[i_proc].find(key) != point_index_map[i_proc].end())
205 52184 : continue;
206 :
207 22163 : point_index_map[i_proc][key] = outgoing_points[i_proc].size();
208 22163 : outgoing_points[i_proc].push_back(std::move(transformed_point));
209 22163 : point_found = true;
210 : } // if
211 : } // i_from
212 : } // i_proc
213 73800 : offset++;
214 : } // point
215 :
216 81 : } // else
217 81 : }
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 740 : 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 740 : unsigned int local_start = 0;
227 923 : 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 1480 : for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; ++i_from)
235 : {
236 740 : local_bboxes[i_from] = bboxes[local_start + i_from];
237 : }
238 : }
239 :
240 : // Setup the local mesh functions.
241 740 : std::vector<MeshFunction> local_meshfuns;
242 740 : local_meshfuns.reserve(_from_problems.size());
243 1480 : for (unsigned int i_from = 0; i_from < _from_problems.size(); ++i_from)
244 : {
245 740 : FEProblemBase & from_problem = *_from_problems[i_from];
246 : MooseVariableFEBase & from_var =
247 740 : from_problem.getVariable(0,
248 740 : _from_var_names[i],
249 : Moose::VarKindType::VAR_ANY,
250 : Moose::VarFieldType::VAR_FIELD_STANDARD);
251 740 : System & from_sys = from_var.sys().system();
252 740 : unsigned int from_var_num = from_sys.variable_number(from_var.name());
253 :
254 1480 : local_meshfuns.emplace_back(getEquationSystem(from_problem, _displaced_source_mesh),
255 740 : *from_sys.current_local_solution,
256 : from_sys.get_dof_map(),
257 : from_var_num);
258 740 : local_meshfuns.back().init();
259 740 : 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 980 : [this, &local_meshfuns, &local_bboxes](
273 : processor_id_type /*pid*/,
274 : const std::vector<Point> & incoming_points,
275 601647 : std::vector<std::pair<Real, unsigned int>> & vals_ids_for_incoming_points)
276 : {
277 980 : vals_ids_for_incoming_points.resize(incoming_points.size(), std::make_pair(OutOfMeshValue, 0));
278 73306 : for (MooseIndex(incoming_points.size()) i_pt = 0; i_pt < incoming_points.size(); ++i_pt)
279 : {
280 72326 : 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 144652 : for (MooseIndex(_from_problems.size()) i_from = 0;
285 216978 : i_from < _from_problems.size() &&
286 72326 : vals_ids_for_incoming_points[i_pt].first == OutOfMeshValue;
287 : ++i_from)
288 : {
289 72326 : if (local_bboxes[i_from].contains_point(pt))
290 : {
291 : const auto from_global_num =
292 72326 : _current_direction == TO_MULTIAPP ? 0 : _from_local2global_map[i_from];
293 : // Use mesh function to compute interpolation values
294 144652 : vals_ids_for_incoming_points[i_pt].first =
295 72326 : (local_meshfuns[i_from])(_from_transforms[from_global_num]->mapBack(pt));
296 : // Record problem ID as well
297 72326 : switch (_current_direction)
298 : {
299 23039 : case FROM_MULTIAPP:
300 23039 : vals_ids_for_incoming_points[i_pt].second = _from_local2global_map[i_from];
301 23039 : break;
302 49287 : case TO_MULTIAPP:
303 49287 : vals_ids_for_incoming_points[i_pt].second = _to_local2global_map[i_from];
304 49287 : break;
305 0 : default:
306 0 : mooseError("Unsupported direction");
307 : }
308 : }
309 : }
310 : }
311 980 : };
312 :
313 : // Incoming values and APP ids for outgoing points
314 740 : 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 980 : [&incoming_vals_ids](
318 : processor_id_type pid,
319 : const std::vector<Point> & /*my_outgoing_points*/,
320 1960 : const std::vector<std::pair<Real, unsigned int>> & vals_ids_for_outgoing_points)
321 : {
322 : // This lambda function might be called multiple times
323 980 : incoming_vals_ids[pid].reserve(vals_ids_for_outgoing_points.size());
324 : // Copy data for processor 'pid'
325 980 : std::copy(vals_ids_for_outgoing_points.begin(),
326 : vals_ids_for_outgoing_points.end(),
327 980 : std::back_inserter(incoming_vals_ids[pid]));
328 1720 : };
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 740 : const std::pair<Real, unsigned int> * ex = nullptr;
333 740 : libMesh::Parallel::pull_parallel_vector_data(
334 740 : comm(), outgoing_points, gather_functor, action_functor, ex);
335 :
336 1476 : for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
337 : {
338 740 : const auto to_global_num = _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
339 740 : System * to_sys = find_sys(*_to_es[i_to], _to_var_names[i]);
340 :
341 740 : unsigned int sys_num = to_sys->number();
342 740 : unsigned int var_num = to_sys->variable_number(_to_var_names[i]);
343 :
344 740 : NumericVector<Real> * solution = nullptr;
345 740 : switch (_current_direction)
346 : {
347 407 : case TO_MULTIAPP:
348 407 : solution = &getTransferVector(i_to, _to_var_names[i]);
349 407 : break;
350 333 : case FROM_MULTIAPP:
351 333 : solution = to_sys->solution.get();
352 333 : break;
353 0 : default:
354 0 : mooseError("Unknown direction");
355 : }
356 :
357 740 : MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
358 740 : auto & fe_type = to_sys->variable_type(var_num);
359 740 : bool is_constant = fe_type.order == CONSTANT;
360 740 : bool is_nodal = fe_type.family == LAGRANGE;
361 :
362 740 : if (is_nodal)
363 : {
364 96295 : for (const auto & node : to_mesh->local_node_ptr_range())
365 : {
366 : // Skip this node if the variable has no dofs at it.
367 47822 : if (node->n_dofs(sys_num, var_num) < 1)
368 0 : continue;
369 :
370 47822 : unsigned int lowest_app_rank = libMesh::invalid_uint;
371 47822 : Real best_val = 0.;
372 47822 : bool point_found = false;
373 107103 : for (auto & group : incoming_vals_ids)
374 : {
375 : // Skip this proc if the node wasn't in it's bounding boxes.
376 59281 : 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 59281 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
381 9972 : continue;
382 :
383 49947 : 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 49947 : if (_current_direction == FROM_MULTIAPP)
388 : {
389 22099 : 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 49947 : if (group.second[i_pt].first == OutOfMeshValue)
395 638 : continue;
396 :
397 49309 : best_val = group.second[i_pt].first;
398 49309 : point_found = true;
399 : }
400 :
401 47822 : 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 47818 : dof_id_type dof = node->dof_number(sys_num, var_num, 0);
406 47818 : solution->set(dof, best_val);
407 655 : }
408 : }
409 : else // Elemental
410 : {
411 81 : std::vector<Point> points;
412 81 : std::vector<dof_id_type> point_ids;
413 39681 : 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 19800 : if (elem->n_dofs(sys_num, var_num) < 1)
417 0 : continue;
418 :
419 19800 : points.clear();
420 19800 : point_ids.clear();
421 : // grab sample points
422 : // for constant shape function, we take the element centroid
423 19800 : if (is_constant)
424 : {
425 1800 : points.push_back(elem->vertex_average());
426 1800 : 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 90000 : for (auto & node : elem->node_ref_range())
434 : {
435 72000 : points.push_back(node);
436 72000 : point_ids.push_back(node.id());
437 : }
438 : }
439 :
440 19800 : auto n_points = points.size();
441 19800 : unsigned int n_comp = elem->n_comp(sys_num, var_num);
442 : // We assume each point corresponds to one component of elemental variable
443 19800 : 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 93600 : for (unsigned int offset = 0; offset < n_points; offset++)
449 : {
450 73800 : unsigned int lowest_app_rank = libMesh::invalid_uint;
451 73800 : Real best_val = 0;
452 73800 : bool point_found = false;
453 164000 : for (auto & group : incoming_vals_ids)
454 : {
455 : // Skip this proc if the elem wasn't in it's bounding boxes.
456 90200 : std::pair<unsigned int, dof_id_type> key(i_to, point_ids[offset]);
457 90200 : if (point_index_map[group.first].find(key) == point_index_map[group.first].end())
458 17553 : continue;
459 :
460 74347 : 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 74347 : if (_current_direction == FROM_MULTIAPP)
465 : {
466 940 : 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 74347 : if (group.second[i_pt].first == OutOfMeshValue)
472 1700 : continue;
473 :
474 72647 : best_val = group.second[i_pt].first;
475 72647 : point_found = true;
476 : }
477 :
478 73800 : 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 73800 : dof_id_type dof = elem->dof_number(sys_num, var_num, offset);
484 73800 : solution->set(dof, best_val);
485 : } // point
486 81 : } // element
487 81 : }
488 736 : solution->close();
489 736 : to_sys->update();
490 : }
491 736 : }
|
