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 "MultiAppGeometricInterpolationTransfer.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 "MultiApp.h"
19 : #include "MooseAppCoordTransform.h"
20 :
21 : #include "libmesh/parallel_algebra.h"
22 : #include "libmesh/meshfree_interpolation.h"
23 : #include "libmesh/system.h"
24 : #include "libmesh/radial_basis_interpolation.h"
25 :
26 : using namespace libMesh;
27 :
28 : registerMooseObject("MooseApp", MultiAppGeometricInterpolationTransfer);
29 : registerMooseObjectRenamed("MooseApp",
30 : MultiAppInterpolationTransfer,
31 : "12/31/2023 24:00",
32 : MultiAppGeometricInterpolationTransfer);
33 :
34 : InputParameters
35 29186 : MultiAppGeometricInterpolationTransfer::validParams()
36 : {
37 29186 : InputParameters params = MultiAppConservativeTransfer::validParams();
38 29186 : params.addClassDescription(
39 : "Transfers the value to the target domain from a combination/interpolation of the values on "
40 : "the nearest nodes in the source domain, using coefficients based on the distance to each "
41 : "node.");
42 87558 : params.addParam<unsigned int>(
43 58372 : "num_points", 3, "The number of nearest points to use for interpolation.");
44 87558 : params.addParam<Real>(
45 58372 : "power", 2, "The polynomial power to use for calculation of the decay in the interpolation.");
46 :
47 29186 : MooseEnum interp_type("inverse_distance radial_basis", "inverse_distance");
48 29186 : params.addParam<MooseEnum>("interp_type", interp_type, "The algorithm to use for interpolation.");
49 :
50 87558 : params.addParam<Real>("radius",
51 58372 : -1,
52 : "Radius to use for radial_basis interpolation. If negative "
53 : "then the radius is taken as the max distance between "
54 : "points.");
55 87558 : params.addParam<Real>(
56 : "shrink_gap_width",
57 58372 : 0,
58 : "gap width with which we want to temporarily shrink mesh in transfering solution");
59 :
60 29186 : MooseEnum shrink_type("SOURCE TARGET", "SOURCE");
61 29186 : params.addParam<MooseEnum>("shrink_mesh", shrink_type, "Which mesh we want to shrink");
62 :
63 29186 : params.addParam<std::vector<SubdomainName>>(
64 : "exclude_gap_blocks",
65 : {},
66 : "Gap subdomains we want to exclude when constructing/using virtually translated points");
67 :
68 87558 : params.addParam<Real>("distance_tol",
69 58372 : 1e-10,
70 : "If the distance between two points is smaller than distance_tol, two "
71 : "points will be considered as identical");
72 :
73 58372 : return params;
74 29186 : }
75 :
76 328 : MultiAppGeometricInterpolationTransfer::MultiAppGeometricInterpolationTransfer(
77 328 : const InputParameters & parameters)
78 : : MultiAppConservativeTransfer(parameters),
79 328 : _num_points(getParam<unsigned int>("num_points")),
80 328 : _power(getParam<Real>("power")),
81 328 : _interp_type(getParam<MooseEnum>("interp_type")),
82 328 : _radius(getParam<Real>("radius")),
83 328 : _shrink_gap_width(getParam<Real>("shrink_gap_width")),
84 328 : _shrink_mesh(getParam<MooseEnum>("shrink_mesh")),
85 328 : _exclude_gap_blocks(getParam<std::vector<SubdomainName>>("exclude_gap_blocks")),
86 656 : _distance_tol(getParam<Real>("distance_tol"))
87 : {
88 : // This transfer does not work with DistributedMesh
89 328 : _fe_problem.mesh().errorIfDistributedMesh("MultiAppGeometricInterpolationTransfer");
90 :
91 328 : if (_to_var_names.size() != 1)
92 0 : paramError("variable", " Support single to-variable only ");
93 :
94 328 : if (_from_var_names.size() != 1)
95 0 : paramError("source_variable", " Support single from-variable only ");
96 328 : }
97 :
98 : void
99 35301 : MultiAppGeometricInterpolationTransfer::subdomainIDsNode(MooseMesh & mesh,
100 : const Node & node,
101 : std::set<subdomain_id_type> & subdomainids)
102 : {
103 : // We need this map to figure out to which subdomains a given mesh point is attached
104 : // We can not make mesh const here because we may need to create a node-to-elems map
105 : // if it does not exists
106 35301 : auto & node_to_elem = mesh.nodeToElemMap();
107 35301 : auto node_to_elem_pair = node_to_elem.find(node.id());
108 :
109 35301 : if (node_to_elem_pair == node_to_elem.end())
110 0 : mooseError("Can not find elements for node ", node.id());
111 :
112 35301 : subdomainids.clear();
113 : // Add all subdomain IDs that are attached to node
114 169701 : for (auto element : node_to_elem_pair->second)
115 : {
116 134400 : auto & elem = mesh.getMesh().elem_ref(element);
117 134400 : auto subdomain = elem.subdomain_id();
118 :
119 134400 : subdomainids.insert(subdomain);
120 : }
121 35301 : }
122 :
123 : void
124 563 : MultiAppGeometricInterpolationTransfer::fillSourceInterpolationPoints(
125 : FEProblemBase & from_problem,
126 : const MooseVariableFieldBase & from_var,
127 : const MultiAppCoordTransform & from_app_transform,
128 : std::unique_ptr<InverseDistanceInterpolation<LIBMESH_DIM>> & idi)
129 : {
130 563 : auto & from_moose_mesh = from_problem.mesh(_displaced_source_mesh);
131 563 : const auto & from_mesh = from_moose_mesh.getMesh();
132 :
133 : // Moose system
134 563 : const SystemBase & from_system_base = from_var.sys();
135 : // libmesh system
136 563 : const System & from_sys = from_system_base.system();
137 :
138 : // System number and var number
139 563 : auto from_sys_num = from_sys.number();
140 563 : auto from_var_num = from_sys.variable_number(from_var.name());
141 :
142 : // Check FE type so we can figure out how to sample points
143 563 : const auto & fe_type = from_sys.variable_type(from_var_num);
144 563 : bool from_is_constant = fe_type.order == CONSTANT;
145 563 : bool from_is_nodal = fe_type.family == LAGRANGE;
146 :
147 : // Currently, for an elemental variable, we support the constant and first order
148 563 : if (fe_type.order > FIRST && !from_is_nodal)
149 0 : mooseError("We don't currently support second order or higher elemental variable ");
150 :
151 : // Containers for points and values
152 : // We later will push data into these containers
153 563 : std::vector<Point> & src_pts(idi->get_source_points());
154 563 : std::vector<Number> & src_vals(idi->get_source_vals());
155 :
156 : // How much we want to translate mesh if users ask
157 563 : std::unordered_map<dof_id_type, Point> from_tranforms;
158 563 : std::set<subdomain_id_type> exclude_block_ids;
159 563 : if (_shrink_gap_width > 0 && _shrink_mesh == "source")
160 : {
161 18 : computeTransformation(from_moose_mesh, from_tranforms);
162 18 : auto exclude_subdomainids = from_moose_mesh.getSubdomainIDs(_exclude_gap_blocks);
163 18 : exclude_block_ids.insert(exclude_subdomainids.begin(), exclude_subdomainids.end());
164 18 : }
165 :
166 : // The solution from the system with which the from_var is associated
167 563 : const NumericVector<Number> & from_solution = *from_sys.solution;
168 :
169 563 : std::set<subdomain_id_type> subdomainids;
170 563 : std::vector<subdomain_id_type> include_block_ids;
171 563 : if (from_is_nodal)
172 : {
173 136639 : for (const auto * const from_node : from_mesh.local_node_ptr_range())
174 : {
175 : // Assuming LAGRANGE!
176 68074 : if (from_node->n_comp(from_sys_num, from_var_num) == 0)
177 8036 : continue;
178 :
179 68074 : Point translate(0);
180 :
181 68074 : if (from_tranforms.size() > 0)
182 : {
183 23534 : subdomainIDsNode(const_cast<MooseMesh &>(from_moose_mesh), *from_node, subdomainids);
184 : // Check if node is excluded
185 : // Node will be excluded if it is in the interior of excluded subdomains
186 23534 : include_block_ids.clear();
187 23534 : include_block_ids.resize(std::max(subdomainids.size(), exclude_block_ids.size()));
188 23534 : auto it = std::set_difference(subdomainids.begin(),
189 : subdomainids.end(),
190 : exclude_block_ids.begin(),
191 : exclude_block_ids.end(),
192 : include_block_ids.begin());
193 :
194 23534 : include_block_ids.resize(it - include_block_ids.begin());
195 : // Node is not excluded
196 23534 : if (include_block_ids.size())
197 15498 : translate = from_tranforms[*include_block_ids.begin()];
198 : else
199 8036 : continue;
200 : }
201 :
202 : // Push value and point to KDTree
203 60038 : dof_id_type from_dof = from_node->dof_number(from_sys_num, from_var_num, 0);
204 60038 : src_vals.push_back(from_solution(from_dof));
205 60038 : src_pts.push_back(from_app_transform(*from_node) + translate);
206 491 : }
207 : }
208 : else
209 : {
210 72 : std::vector<Point> points;
211 72 : for (const auto * const from_elem :
212 12544 : as_range(from_mesh.local_elements_begin(), from_mesh.local_elements_end()))
213 : {
214 : // Skip this element if the variable has no dofs at it.
215 6200 : if (from_elem->n_dofs(from_sys_num, from_var_num) < 1)
216 0 : continue;
217 :
218 6200 : points.clear();
219 6200 : if (from_is_constant)
220 6200 : points.push_back(from_elem->vertex_average());
221 : else
222 0 : for (const auto & node : from_elem->node_ref_range())
223 0 : points.push_back(node);
224 :
225 6200 : unsigned int n_comp = from_elem->n_comp(from_sys_num, from_var_num);
226 6200 : auto n_points = points.size();
227 : // We assume each point corresponds to one component of elemental variable
228 6200 : if (n_points != n_comp)
229 0 : mooseError(" Number of points ",
230 : n_points,
231 : " does not equal to number of variable components ",
232 : n_comp);
233 :
234 6200 : unsigned int offset = 0;
235 :
236 6200 : Point translate(0);
237 :
238 6200 : if (from_tranforms.size() > 0)
239 : {
240 0 : auto subdomain = from_elem->subdomain_id();
241 :
242 0 : if (subdomain == Moose::INVALID_BLOCK_ID)
243 0 : mooseError("subdomain id does not make sense", subdomain);
244 :
245 : // subdomain is not excluded
246 0 : if (exclude_block_ids.find(subdomain) == exclude_block_ids.end())
247 0 : translate = from_tranforms[subdomain];
248 : else
249 0 : continue;
250 : }
251 :
252 12400 : for (const auto & point : points)
253 : {
254 6200 : dof_id_type from_dof = from_elem->dof_number(from_sys_num, from_var_num, offset++);
255 6200 : src_vals.push_back(from_solution(from_dof));
256 6200 : src_pts.push_back(from_app_transform(point) + translate);
257 : }
258 72 : }
259 72 : }
260 563 : }
261 :
262 : void
263 442 : MultiAppGeometricInterpolationTransfer::interpolateTargetPoints(
264 : FEProblemBase & to_problem,
265 : MooseVariableFieldBase & to_var,
266 : NumericVector<Real> & to_solution,
267 : const MultiAppCoordTransform & to_app_transform,
268 : const std::unique_ptr<InverseDistanceInterpolation<LIBMESH_DIM>> & idi)
269 : {
270 : // Moose system
271 442 : SystemBase & to_system_base = to_var.sys();
272 : // libmesh system
273 442 : System & to_sys = to_system_base.system();
274 :
275 : // System number and var number
276 442 : auto to_sys_num = to_sys.number();
277 442 : auto to_var_num = to_sys.variable_number(to_var.name());
278 :
279 442 : const MooseMesh & to_moose_mesh = to_problem.mesh(_displaced_target_mesh);
280 442 : const MeshBase & to_mesh = to_moose_mesh.getMesh();
281 :
282 : // Compute transform info
283 442 : std::unordered_map<dof_id_type, Point> to_tranforms;
284 442 : std::set<subdomain_id_type> exclude_block_ids;
285 442 : if (_shrink_gap_width > 0 && _shrink_mesh == "target")
286 : {
287 18 : computeTransformation(to_moose_mesh, to_tranforms);
288 18 : auto exclude_subdomainids = to_moose_mesh.getSubdomainIDs(_exclude_gap_blocks);
289 18 : exclude_block_ids.insert(exclude_subdomainids.begin(), exclude_subdomainids.end());
290 18 : }
291 :
292 442 : const auto & to_fe_type = to_sys.variable_type(to_var_num);
293 442 : bool to_is_constant = to_fe_type.order == CONSTANT;
294 442 : bool to_is_nodal = to_fe_type.family == LAGRANGE;
295 :
296 442 : if (to_fe_type.order > FIRST && !to_is_nodal)
297 0 : mooseError("We don't currently support second order or higher elemental variable ");
298 :
299 442 : std::set<subdomain_id_type> subdomainids;
300 442 : std::vector<subdomain_id_type> include_block_ids;
301 442 : std::vector<Point> pts;
302 442 : std::vector<Number> vals;
303 442 : if (to_is_nodal)
304 : {
305 94852 : for (const auto * const node : to_mesh.local_node_ptr_range())
306 : {
307 47256 : if (node->n_dofs(to_sys_num, to_var_num) <= 0) // If this variable has dofs at this node
308 4018 : continue;
309 :
310 47256 : Point translate(0);
311 47256 : if (to_tranforms.size() > 0)
312 : {
313 11767 : subdomainIDsNode(const_cast<MooseMesh &>(to_moose_mesh), *node, subdomainids);
314 : // Check if node is excluded
315 : // Node will be excluded if it is in the interior of excluded subdomains
316 11767 : include_block_ids.clear();
317 11767 : include_block_ids.resize(std::max(subdomainids.size(), exclude_block_ids.size()));
318 11767 : auto it = std::set_difference(subdomainids.begin(),
319 : subdomainids.end(),
320 : exclude_block_ids.begin(),
321 : exclude_block_ids.end(),
322 : include_block_ids.begin());
323 11767 : include_block_ids.resize(it - include_block_ids.begin());
324 11767 : if (include_block_ids.size())
325 7749 : translate = to_tranforms[*include_block_ids.begin()];
326 : else
327 4018 : continue;
328 : }
329 :
330 43238 : pts.clear();
331 43238 : pts.push_back(to_app_transform(*node) + translate);
332 43238 : vals.resize(1);
333 :
334 129714 : idi->interpolate_field_data({_to_var_name}, pts, vals);
335 43238 : dof_id_type dof = node->dof_number(to_sys_num, to_var_num, 0);
336 43238 : to_solution.set(dof, vals.front());
337 340 : }
338 : }
339 : else // Elemental
340 : {
341 102 : std::vector<Point> points;
342 102 : for (const auto * const elem :
343 60204 : as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
344 : {
345 : // Skip this element if the variable has no dofs at it.
346 30000 : if (elem->n_dofs(to_sys_num, to_var_num) < 1)
347 10360 : continue;
348 :
349 24120 : points.clear();
350 24120 : if (to_is_constant)
351 24120 : points.push_back(elem->vertex_average());
352 : else
353 0 : for (const auto & node : elem->node_ref_range())
354 0 : points.push_back(node);
355 :
356 24120 : auto n_points = points.size();
357 24120 : unsigned int n_comp = elem->n_comp(to_sys_num, to_var_num);
358 : // We assume each point corresponds to one component of elemental variable
359 24120 : if (n_points != n_comp)
360 0 : mooseError(" Number of points ",
361 : n_points,
362 : " does not equal to number of variable components ",
363 : n_comp);
364 :
365 24120 : Point translate(0);
366 :
367 24120 : if (to_tranforms.size() > 0)
368 : {
369 11200 : auto subdomain = elem->subdomain_id();
370 :
371 11200 : if (subdomain == Moose::INVALID_BLOCK_ID)
372 0 : mooseError("subdomain id does not make sense", subdomain);
373 :
374 11200 : if (exclude_block_ids.find(subdomain) == exclude_block_ids.end())
375 6720 : translate = to_tranforms[subdomain];
376 : else
377 4480 : continue;
378 : }
379 :
380 19640 : unsigned int offset = 0;
381 39280 : for (const auto & point : points)
382 : {
383 19640 : pts.clear();
384 19640 : pts.push_back(to_app_transform(point) + translate);
385 19640 : vals.resize(1);
386 :
387 58920 : idi->interpolate_field_data({_to_var_name}, pts, vals);
388 19640 : dof_id_type dof = elem->dof_number(to_sys_num, to_var_num, offset++);
389 19640 : to_solution.set(dof, vals.front());
390 : } // point
391 102 : } // auto elem
392 102 : } // else
393 :
394 442 : to_solution.close();
395 442 : to_sys.update();
396 63320 : }
397 :
398 : void
399 442 : MultiAppGeometricInterpolationTransfer::execute()
400 : {
401 442 : TIME_SECTION("MultiAppGeometricInterpolationTransfer::execute()",
402 : 5,
403 : "Transferring variables based on node interpolation");
404 :
405 : const FEProblemBase & fe_problem =
406 442 : hasFromMultiApp() ? getFromMultiApp()->problemBase() : getToMultiApp()->problemBase();
407 442 : std::unique_ptr<InverseDistanceInterpolation<LIBMESH_DIM>> idi;
408 442 : switch (_interp_type)
409 : {
410 398 : case 0:
411 796 : idi = std::make_unique<InverseDistanceInterpolation<LIBMESH_DIM>>(
412 796 : fe_problem.comm(), _num_points, _power);
413 398 : break;
414 44 : case 1:
415 44 : idi = std::make_unique<RadialBasisInterpolation<LIBMESH_DIM>>(fe_problem.comm(), _radius);
416 44 : break;
417 0 : default:
418 0 : mooseError("Unknown interpolation type!");
419 : }
420 :
421 1326 : idi->set_field_variables({_to_var_name});
422 :
423 442 : switch (_current_direction)
424 : {
425 135 : case TO_MULTIAPP:
426 : {
427 135 : FEProblemBase & from_problem = getToMultiApp()->problemBase();
428 135 : const auto & from_var = from_problem.getVariable(
429 : 0, _from_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
430 :
431 : mooseAssert(_from_transforms.size() == 1, "This should be size 1");
432 135 : fillSourceInterpolationPoints(from_problem, from_var, *_from_transforms[0], idi);
433 :
434 : // We have only set local values - prepare for use by gathering remote gata
435 135 : idi->prepare_for_use();
436 :
437 270 : for (unsigned int i = 0; i < getToMultiApp()->numGlobalApps(); i++)
438 : {
439 135 : if (getToMultiApp()->hasLocalApp(i))
440 : {
441 135 : auto & to_problem = getToMultiApp()->appProblemBase(i);
442 135 : Moose::ScopedCommSwapper swapper(to_problem.comm().get());
443 135 : auto & to_var = to_problem.getVariable(0,
444 : _to_var_name,
445 : Moose::VarKindType::VAR_ANY,
446 : Moose::VarFieldType::VAR_FIELD_STANDARD);
447 :
448 135 : auto & to_solution = getToMultiApp()->appTransferVector(i, _to_var_name);
449 :
450 135 : interpolateTargetPoints(to_problem, to_var, to_solution, *_to_transforms[i], idi);
451 135 : }
452 : }
453 :
454 135 : break;
455 : }
456 :
457 307 : case FROM_MULTIAPP:
458 : {
459 885 : for (unsigned int i = 0; i < getFromMultiApp()->numGlobalApps(); i++)
460 : {
461 578 : if (getFromMultiApp()->hasLocalApp(i))
462 : {
463 428 : auto & from_problem = getFromMultiApp()->appProblemBase(i);
464 428 : Moose::ScopedCommSwapper swapper(from_problem.comm().get());
465 428 : const auto & from_var = from_problem.getVariable(0,
466 : _from_var_name,
467 : Moose::VarKindType::VAR_ANY,
468 : Moose::VarFieldType::VAR_FIELD_STANDARD);
469 :
470 428 : fillSourceInterpolationPoints(from_problem, from_var, *_from_transforms[i], idi);
471 428 : }
472 : }
473 :
474 307 : idi->prepare_for_use();
475 :
476 307 : FEProblemBase & to_problem = getFromMultiApp()->problemBase();
477 307 : MooseVariableFieldBase & to_var = to_problem.getVariable(
478 : 0, _to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
479 :
480 307 : auto & to_solution = *to_var.sys().system().solution;
481 :
482 : mooseAssert(_to_transforms.size() == 1, "This should be size 1");
483 307 : interpolateTargetPoints(to_problem, to_var, to_solution, *_to_transforms[0], idi);
484 :
485 307 : break;
486 : }
487 0 : default:
488 : {
489 0 : mooseError("Unsupported transfer direction ", _current_direction);
490 : break;
491 : }
492 : }
493 884 : }
494 :
495 : void
496 36 : MultiAppGeometricInterpolationTransfer::computeTransformation(
497 : const MooseMesh & mesh, std::unordered_map<dof_id_type, Point> & transformation)
498 : {
499 36 : auto & libmesh_mesh = mesh.getMesh();
500 :
501 36 : auto & subdomainids = mesh.meshSubdomains();
502 :
503 36 : subdomain_id_type max_subdomain_id = 0;
504 :
505 : // max_subdomain_id will be used to represent the center of the entire domain
506 216 : for (auto subdomain_id : subdomainids)
507 : {
508 180 : max_subdomain_id = max_subdomain_id > subdomain_id ? max_subdomain_id : subdomain_id;
509 : }
510 :
511 36 : max_subdomain_id += 1;
512 :
513 36 : std::unordered_map<dof_id_type, Point> subdomain_centers;
514 36 : std::unordered_map<dof_id_type, dof_id_type> nelems;
515 :
516 36 : for (auto & elem :
517 59572 : as_range(libmesh_mesh.local_elements_begin(), libmesh_mesh.local_elements_end()))
518 : {
519 : // Compute center of the entire domain
520 59500 : subdomain_centers[max_subdomain_id] += elem->vertex_average();
521 59500 : nelems[max_subdomain_id] += 1;
522 :
523 59500 : auto subdomain = elem->subdomain_id();
524 :
525 59500 : if (subdomain == Moose::INVALID_BLOCK_ID)
526 0 : mooseError("block is invalid");
527 :
528 : // Centers for subdomains
529 59500 : subdomain_centers[subdomain] += elem->vertex_average();
530 :
531 59500 : nelems[subdomain] += 1;
532 36 : }
533 :
534 36 : comm().sum(subdomain_centers);
535 :
536 36 : comm().sum(nelems);
537 :
538 36 : subdomain_centers[max_subdomain_id] /= nelems[max_subdomain_id];
539 :
540 216 : for (auto subdomain_id : subdomainids)
541 : {
542 180 : subdomain_centers[subdomain_id] /= nelems[subdomain_id];
543 : }
544 :
545 : // Compute unit vectors representing directions in which we want to shrink mesh
546 : // The unit vectors is scaled by 'shrink_gap_width'
547 36 : transformation.clear();
548 216 : for (auto subdomain_id : subdomainids)
549 : {
550 180 : transformation[subdomain_id] =
551 180 : subdomain_centers[max_subdomain_id] - subdomain_centers[subdomain_id];
552 :
553 180 : auto norm = transformation[subdomain_id].norm();
554 :
555 : // The current subdomain is the center of the entire domain,
556 : // then we do not move this subdomain
557 180 : if (norm > _distance_tol)
558 144 : transformation[subdomain_id] /= norm;
559 :
560 180 : transformation[subdomain_id] *= _shrink_gap_width;
561 : }
562 36 : }
|
