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 29244 : MultiAppGeometricInterpolationTransfer::validParams()
36 : {
37 29244 : InputParameters params = MultiAppConservativeTransfer::validParams();
38 29244 : 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 87732 : params.addParam<unsigned int>(
43 58488 : "num_points", 3, "The number of nearest points to use for interpolation.");
44 87732 : params.addParam<Real>(
45 58488 : "power", 2, "The polynomial power to use for calculation of the decay in the interpolation.");
46 :
47 29244 : MooseEnum interp_type("inverse_distance radial_basis", "inverse_distance");
48 29244 : params.addParam<MooseEnum>("interp_type", interp_type, "The algorithm to use for interpolation.");
49 :
50 87732 : params.addParam<Real>("radius",
51 58488 : -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 87732 : params.addParam<Real>(
56 : "shrink_gap_width",
57 58488 : 0,
58 : "gap width with which we want to temporarily shrink mesh in transfering solution");
59 :
60 29244 : MooseEnum shrink_type("SOURCE TARGET", "SOURCE");
61 29244 : params.addParam<MooseEnum>("shrink_mesh", shrink_type, "Which mesh we want to shrink");
62 :
63 29244 : 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 87732 : params.addParam<Real>("distance_tol",
69 58488 : 1e-10,
70 : "If the distance between two points is smaller than distance_tol, two "
71 : "points will be considered as identical");
72 :
73 58488 : return params;
74 29244 : }
75 :
76 357 : MultiAppGeometricInterpolationTransfer::MultiAppGeometricInterpolationTransfer(
77 357 : const InputParameters & parameters)
78 : : MultiAppConservativeTransfer(parameters),
79 357 : _num_points(getParam<unsigned int>("num_points")),
80 357 : _power(getParam<Real>("power")),
81 357 : _interp_type(getParam<MooseEnum>("interp_type")),
82 357 : _radius(getParam<Real>("radius")),
83 357 : _shrink_gap_width(getParam<Real>("shrink_gap_width")),
84 357 : _shrink_mesh(getParam<MooseEnum>("shrink_mesh")),
85 357 : _exclude_gap_blocks(getParam<std::vector<SubdomainName>>("exclude_gap_blocks")),
86 714 : _distance_tol(getParam<Real>("distance_tol"))
87 : {
88 : // This transfer does not work with DistributedMesh
89 357 : _fe_problem.mesh().errorIfDistributedMesh("MultiAppGeometricInterpolationTransfer");
90 :
91 357 : if (_to_var_names.size() != 1)
92 0 : paramError("variable", " Support single to-variable only ");
93 :
94 357 : if (_from_var_names.size() != 1)
95 0 : paramError("source_variable", " Support single from-variable only ");
96 357 : }
97 :
98 : void
99 40344 : 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 40344 : auto & node_to_elem = mesh.nodeToElemMap();
107 40344 : auto node_to_elem_pair = node_to_elem.find(node.id());
108 :
109 40344 : if (node_to_elem_pair == node_to_elem.end())
110 0 : mooseError("Can not find elements for node ", node.id());
111 :
112 40344 : subdomainids.clear();
113 : // Add all subdomain IDs that are attached to node
114 193944 : for (auto element : node_to_elem_pair->second)
115 : {
116 153600 : auto & elem = mesh.getMesh().elem_ref(element);
117 153600 : auto subdomain = elem.subdomain_id();
118 :
119 153600 : subdomainids.insert(subdomain);
120 : }
121 40344 : }
122 :
123 : void
124 661 : 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 661 : auto & from_moose_mesh = from_problem.mesh(_displaced_source_mesh);
131 661 : const auto & from_mesh = from_moose_mesh.getMesh();
132 :
133 : // Moose system
134 661 : const SystemBase & from_system_base = from_var.sys();
135 : // libmesh system
136 661 : const System & from_sys = from_system_base.system();
137 :
138 : // System number and var number
139 661 : auto from_sys_num = from_sys.number();
140 661 : 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 661 : const auto & fe_type = from_sys.variable_type(from_var_num);
144 661 : bool from_is_constant = fe_type.order == CONSTANT;
145 661 : bool from_is_nodal = fe_type.family == LAGRANGE;
146 :
147 : // Currently, for an elemental variable, we support the constant and first order
148 661 : 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 661 : std::vector<Point> & src_pts(idi->get_source_points());
154 661 : std::vector<Number> & src_vals(idi->get_source_vals());
155 :
156 : // How much we want to translate mesh if users ask
157 661 : std::unordered_map<dof_id_type, Point> from_tranforms;
158 661 : std::set<subdomain_id_type> exclude_block_ids;
159 661 : if (_shrink_gap_width > 0 && _shrink_mesh == "source")
160 : {
161 20 : computeTransformation(from_moose_mesh, from_tranforms);
162 20 : auto exclude_subdomainids = from_moose_mesh.getSubdomainIDs(_exclude_gap_blocks);
163 20 : exclude_block_ids.insert(exclude_subdomainids.begin(), exclude_subdomainids.end());
164 20 : }
165 :
166 : // The solution from the system with which the from_var is associated
167 661 : const NumericVector<Number> & from_solution = *from_sys.solution;
168 :
169 661 : std::set<subdomain_id_type> subdomainids;
170 661 : std::vector<subdomain_id_type> include_block_ids;
171 661 : if (from_is_nodal)
172 : {
173 160607 : for (const auto * const from_node : from_mesh.local_node_ptr_range())
174 : {
175 : // Assuming LAGRANGE!
176 80016 : if (from_node->n_comp(from_sys_num, from_var_num) == 0)
177 9184 : continue;
178 :
179 80016 : Point translate(0);
180 :
181 80016 : if (from_tranforms.size() > 0)
182 : {
183 26896 : 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 26896 : include_block_ids.clear();
187 26896 : include_block_ids.resize(std::max(subdomainids.size(), exclude_block_ids.size()));
188 26896 : 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 26896 : include_block_ids.resize(it - include_block_ids.begin());
195 : // Node is not excluded
196 26896 : if (include_block_ids.size())
197 17712 : translate = from_tranforms[*include_block_ids.begin()];
198 : else
199 9184 : continue;
200 : }
201 :
202 : // Push value and point to KDTree
203 70832 : dof_id_type from_dof = from_node->dof_number(from_sys_num, from_var_num, 0);
204 70832 : src_vals.push_back(from_solution(from_dof));
205 70832 : src_pts.push_back(from_app_transform(*from_node) + translate);
206 575 : }
207 : }
208 : else
209 : {
210 86 : std::vector<Point> points;
211 86 : for (const auto * const from_elem :
212 15372 : 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 7600 : if (from_elem->n_dofs(from_sys_num, from_var_num) < 1)
216 0 : continue;
217 :
218 7600 : points.clear();
219 7600 : if (from_is_constant)
220 7600 : 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 7600 : unsigned int n_comp = from_elem->n_comp(from_sys_num, from_var_num);
226 7600 : auto n_points = points.size();
227 : // We assume each point corresponds to one component of elemental variable
228 7600 : 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 7600 : unsigned int offset = 0;
235 :
236 7600 : Point translate(0);
237 :
238 7600 : 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 15200 : for (const auto & point : points)
253 : {
254 7600 : dof_id_type from_dof = from_elem->dof_number(from_sys_num, from_var_num, offset++);
255 7600 : src_vals.push_back(from_solution(from_dof));
256 7600 : src_pts.push_back(from_app_transform(point) + translate);
257 : }
258 86 : }
259 86 : }
260 661 : }
261 :
262 : void
263 504 : 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 504 : SystemBase & to_system_base = to_var.sys();
272 : // libmesh system
273 504 : System & to_sys = to_system_base.system();
274 :
275 : // System number and var number
276 504 : auto to_sys_num = to_sys.number();
277 504 : auto to_var_num = to_sys.variable_number(to_var.name());
278 :
279 504 : const MooseMesh & to_moose_mesh = to_problem.mesh(_displaced_target_mesh);
280 504 : const MeshBase & to_mesh = to_moose_mesh.getMesh();
281 :
282 : // Compute transform info
283 504 : std::unordered_map<dof_id_type, Point> to_tranforms;
284 504 : std::set<subdomain_id_type> exclude_block_ids;
285 504 : if (_shrink_gap_width > 0 && _shrink_mesh == "target")
286 : {
287 20 : computeTransformation(to_moose_mesh, to_tranforms);
288 20 : auto exclude_subdomainids = to_moose_mesh.getSubdomainIDs(_exclude_gap_blocks);
289 20 : exclude_block_ids.insert(exclude_subdomainids.begin(), exclude_subdomainids.end());
290 20 : }
291 :
292 504 : const auto & to_fe_type = to_sys.variable_type(to_var_num);
293 504 : bool to_is_constant = to_fe_type.order == CONSTANT;
294 504 : bool to_is_nodal = to_fe_type.family == LAGRANGE;
295 :
296 504 : 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 504 : std::set<subdomain_id_type> subdomainids;
300 504 : std::vector<subdomain_id_type> include_block_ids;
301 504 : std::vector<Point> pts;
302 504 : std::vector<Number> vals;
303 504 : if (to_is_nodal)
304 : {
305 110958 : for (const auto * const node : to_mesh.local_node_ptr_range())
306 : {
307 55286 : if (node->n_dofs(to_sys_num, to_var_num) <= 0) // If this variable has dofs at this node
308 4592 : continue;
309 :
310 55286 : Point translate(0);
311 55286 : if (to_tranforms.size() > 0)
312 : {
313 13448 : 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 13448 : include_block_ids.clear();
317 13448 : include_block_ids.resize(std::max(subdomainids.size(), exclude_block_ids.size()));
318 13448 : 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 13448 : include_block_ids.resize(it - include_block_ids.begin());
324 13448 : if (include_block_ids.size())
325 8856 : translate = to_tranforms[*include_block_ids.begin()];
326 : else
327 4592 : continue;
328 : }
329 :
330 50694 : pts.clear();
331 50694 : pts.push_back(to_app_transform(*node) + translate);
332 50694 : vals.resize(1);
333 :
334 152082 : idi->interpolate_field_data({_to_var_name}, pts, vals);
335 50694 : dof_id_type dof = node->dof_number(to_sys_num, to_var_num, 0);
336 50694 : to_solution.set(dof, vals.front());
337 386 : }
338 : }
339 : else // Elemental
340 : {
341 118 : std::vector<Point> points;
342 118 : for (const auto * const elem :
343 69836 : 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 34800 : if (elem->n_dofs(to_sys_num, to_var_num) < 1)
347 11840 : continue;
348 :
349 28080 : points.clear();
350 28080 : if (to_is_constant)
351 28080 : 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 28080 : auto n_points = points.size();
357 28080 : 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 28080 : 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 28080 : Point translate(0);
366 :
367 28080 : if (to_tranforms.size() > 0)
368 : {
369 12800 : auto subdomain = elem->subdomain_id();
370 :
371 12800 : if (subdomain == Moose::INVALID_BLOCK_ID)
372 0 : mooseError("subdomain id does not make sense", subdomain);
373 :
374 12800 : if (exclude_block_ids.find(subdomain) == exclude_block_ids.end())
375 7680 : translate = to_tranforms[subdomain];
376 : else
377 5120 : continue;
378 : }
379 :
380 22960 : unsigned int offset = 0;
381 45920 : for (const auto & point : points)
382 : {
383 22960 : pts.clear();
384 22960 : pts.push_back(to_app_transform(point) + translate);
385 22960 : vals.resize(1);
386 :
387 68880 : idi->interpolate_field_data({_to_var_name}, pts, vals);
388 22960 : dof_id_type dof = elem->dof_number(to_sys_num, to_var_num, offset++);
389 22960 : to_solution.set(dof, vals.front());
390 : } // point
391 118 : } // auto elem
392 118 : } // else
393 :
394 504 : to_solution.close();
395 504 : to_sys.update();
396 74158 : }
397 :
398 : void
399 504 : MultiAppGeometricInterpolationTransfer::execute()
400 : {
401 504 : TIME_SECTION("MultiAppGeometricInterpolationTransfer::execute()",
402 : 5,
403 : "Transferring variables based on node interpolation");
404 :
405 : const FEProblemBase & fe_problem =
406 504 : hasFromMultiApp() ? getFromMultiApp()->problemBase() : getToMultiApp()->problemBase();
407 504 : std::unique_ptr<InverseDistanceInterpolation<LIBMESH_DIM>> idi;
408 504 : switch (_interp_type)
409 : {
410 452 : case 0:
411 904 : idi = std::make_unique<InverseDistanceInterpolation<LIBMESH_DIM>>(
412 904 : fe_problem.comm(), _num_points, _power);
413 452 : break;
414 52 : case 1:
415 52 : idi = std::make_unique<RadialBasisInterpolation<LIBMESH_DIM>>(fe_problem.comm(), _radius);
416 52 : break;
417 0 : default:
418 0 : mooseError("Unknown interpolation type!");
419 : }
420 :
421 1512 : idi->set_field_variables({_to_var_name});
422 :
423 504 : switch (_current_direction)
424 : {
425 157 : case TO_MULTIAPP:
426 : {
427 157 : FEProblemBase & from_problem = getToMultiApp()->problemBase();
428 157 : 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 157 : 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 157 : idi->prepare_for_use();
436 :
437 314 : for (unsigned int i = 0; i < getToMultiApp()->numGlobalApps(); i++)
438 : {
439 157 : if (getToMultiApp()->hasLocalApp(i))
440 : {
441 157 : auto & to_problem = getToMultiApp()->appProblemBase(i);
442 157 : Moose::ScopedCommSwapper swapper(to_problem.comm().get());
443 157 : auto & to_var = to_problem.getVariable(0,
444 : _to_var_name,
445 : Moose::VarKindType::VAR_ANY,
446 : Moose::VarFieldType::VAR_FIELD_STANDARD);
447 :
448 157 : auto & to_solution = getToMultiApp()->appTransferVector(i, _to_var_name);
449 :
450 157 : interpolateTargetPoints(to_problem, to_var, to_solution, *_to_transforms[i], idi);
451 157 : }
452 : }
453 :
454 157 : break;
455 : }
456 :
457 347 : case FROM_MULTIAPP:
458 : {
459 1001 : for (unsigned int i = 0; i < getFromMultiApp()->numGlobalApps(); i++)
460 : {
461 654 : if (getFromMultiApp()->hasLocalApp(i))
462 : {
463 504 : auto & from_problem = getFromMultiApp()->appProblemBase(i);
464 504 : Moose::ScopedCommSwapper swapper(from_problem.comm().get());
465 504 : 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 504 : fillSourceInterpolationPoints(from_problem, from_var, *_from_transforms[i], idi);
471 504 : }
472 : }
473 :
474 347 : idi->prepare_for_use();
475 :
476 347 : FEProblemBase & to_problem = getFromMultiApp()->problemBase();
477 347 : MooseVariableFieldBase & to_var = to_problem.getVariable(
478 : 0, _to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
479 :
480 347 : auto & to_solution = *to_var.sys().system().solution;
481 :
482 : mooseAssert(_to_transforms.size() == 1, "This should be size 1");
483 347 : interpolateTargetPoints(to_problem, to_var, to_solution, *_to_transforms[0], idi);
484 :
485 347 : break;
486 : }
487 0 : default:
488 : {
489 0 : mooseError("Unsupported transfer direction ", _current_direction);
490 : break;
491 : }
492 : }
493 1008 : }
494 :
495 : void
496 40 : MultiAppGeometricInterpolationTransfer::computeTransformation(
497 : const MooseMesh & mesh, std::unordered_map<dof_id_type, Point> & transformation)
498 : {
499 40 : auto & libmesh_mesh = mesh.getMesh();
500 :
501 40 : auto & subdomainids = mesh.meshSubdomains();
502 :
503 40 : subdomain_id_type max_subdomain_id = 0;
504 :
505 : // max_subdomain_id will be used to represent the center of the entire domain
506 240 : for (auto subdomain_id : subdomainids)
507 : {
508 200 : max_subdomain_id = max_subdomain_id > subdomain_id ? max_subdomain_id : subdomain_id;
509 : }
510 :
511 40 : max_subdomain_id += 1;
512 :
513 40 : std::unordered_map<dof_id_type, Point> subdomain_centers;
514 40 : std::unordered_map<dof_id_type, dof_id_type> nelems;
515 :
516 40 : for (auto & elem :
517 68080 : as_range(libmesh_mesh.local_elements_begin(), libmesh_mesh.local_elements_end()))
518 : {
519 : // Compute center of the entire domain
520 68000 : subdomain_centers[max_subdomain_id] += elem->vertex_average();
521 68000 : nelems[max_subdomain_id] += 1;
522 :
523 68000 : auto subdomain = elem->subdomain_id();
524 :
525 68000 : if (subdomain == Moose::INVALID_BLOCK_ID)
526 0 : mooseError("block is invalid");
527 :
528 : // Centers for subdomains
529 68000 : subdomain_centers[subdomain] += elem->vertex_average();
530 :
531 68000 : nelems[subdomain] += 1;
532 40 : }
533 :
534 40 : comm().sum(subdomain_centers);
535 :
536 40 : comm().sum(nelems);
537 :
538 40 : subdomain_centers[max_subdomain_id] /= nelems[max_subdomain_id];
539 :
540 240 : for (auto subdomain_id : subdomainids)
541 : {
542 200 : 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 40 : transformation.clear();
548 240 : for (auto subdomain_id : subdomainids)
549 : {
550 200 : transformation[subdomain_id] =
551 200 : subdomain_centers[max_subdomain_id] - subdomain_centers[subdomain_id];
552 :
553 200 : 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 200 : if (norm > _distance_tol)
558 160 : transformation[subdomain_id] /= norm;
559 :
560 200 : transformation[subdomain_id] *= _shrink_gap_width;
561 : }
562 40 : }
|
