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 "MultiAppGeneralFieldNearestLocationTransfer.h"
11 :
12 : // MOOSE includes
13 : #include "FEProblem.h"
14 : #include "MooseMesh.h"
15 : #include "MooseTypes.h"
16 : #include "MooseVariableFE.h"
17 : #include "SystemBase.h"
18 : #include "Positions.h"
19 : #include "MooseAppCoordTransform.h"
20 :
21 : #include "libmesh/system.h"
22 :
23 : using namespace libMesh;
24 :
25 : registerMooseObject("MooseApp", MultiAppGeneralFieldNearestLocationTransfer);
26 : registerMooseObjectRenamed("MooseApp",
27 : MultiAppGeneralFieldNearestNodeTransfer,
28 : "12/31/2024 24:00",
29 : MultiAppGeneralFieldNearestLocationTransfer);
30 :
31 : InputParameters
32 12852 : MultiAppGeneralFieldNearestLocationTransfer::validParams()
33 : {
34 12852 : InputParameters params = MultiAppGeneralFieldKDTreeTransferBase::validParams();
35 12852 : params.addClassDescription(
36 : "Transfers field data at the MultiApp position by finding the value at the nearest "
37 : "neighbor(s) in the origin application.");
38 :
39 : // Nearest node is historically more an extrapolation transfer
40 38556 : params.set<Real>("extrapolation_constant") = GeneralFieldTransfer::OutOfMeshValue;
41 25704 : params.suppressParameter<Real>("extrapolation_constant");
42 : // We dont keep track of both point distance to app and to nearest node, so we cannot guarantee
43 : // that the nearest app (among the local apps, not globally) will hold the nearest location.
44 : // However, if the user knows this is true, we can use this heuristic to reduce the number of apps
45 : // that are requested to provide a candidate value. If the user is wrong, then the nearest
46 : // location is used, which can be from the non-nearest app.
47 64260 : params.renameParam("use_nearest_app", "assume_nearest_app_holds_nearest_location", "");
48 :
49 : // the default of node/centroid switching based on the variable is causing lots of mistakes and
50 : // bad results
51 : std::vector<MooseEnum> source_types = {
52 38556 : MooseEnum("nodes centroids variable_default", "variable_default")};
53 38556 : params.addParam<std::vector<MooseEnum>>(
54 : "source_type", source_types, "Where to get the source values from for each source variable");
55 :
56 25704 : return params;
57 64260 : }
58 :
59 3362 : MultiAppGeneralFieldNearestLocationTransfer::MultiAppGeneralFieldNearestLocationTransfer(
60 3362 : const InputParameters & parameters)
61 3362 : : MultiAppGeneralFieldKDTreeTransferBase(parameters)
62 : {
63 3359 : }
64 :
65 : void
66 3248 : MultiAppGeneralFieldNearestLocationTransfer::initialSetup()
67 : {
68 3248 : MultiAppGeneralFieldKDTreeTransferBase::initialSetup();
69 :
70 : // Handle the source types ahead of time
71 6424 : const auto & source_types = getParam<std::vector<MooseEnum>>("source_type");
72 3212 : _source_is_nodes.resize(_from_var_names.size());
73 3212 : _use_zero_dof_for_value.resize(_from_var_names.size());
74 3212 : if (source_types.size() != _from_var_names.size())
75 0 : mooseError("Not enough source types specified for this number of variables. Source types must "
76 : "be specified for transfers with multiple variables");
77 6466 : for (const auto var_index : index_range(_from_var_names))
78 : {
79 : // Local app does not own any of the source problems
80 3260 : if (_from_problems.empty())
81 3 : break;
82 :
83 : // Get some info on the source variable
84 3257 : FEProblemBase & from_problem = *_from_problems[0];
85 : MooseVariableFieldBase & from_var =
86 3257 : from_problem.getVariable(0,
87 3257 : _from_var_names[var_index],
88 : Moose::VarKindType::VAR_ANY,
89 : Moose::VarFieldType::VAR_FIELD_ANY);
90 3257 : System & from_sys = from_var.sys().system();
91 3257 : const auto & fe_type = from_sys.variable_type(from_var.number());
92 :
93 : // Select where to get the origin values
94 3257 : if (source_types[var_index] == "nodes")
95 48 : _source_is_nodes[var_index] = true;
96 3209 : else if (source_types[var_index] == "centroids")
97 48 : _source_is_nodes[var_index] = false;
98 : else
99 : {
100 : // "Nodal" variables are continuous for sure at nodes
101 3161 : if (from_var.isNodal())
102 2357 : _source_is_nodes[var_index] = true;
103 : // for everyone else, we know they are continuous at centroids
104 : else
105 804 : _source_is_nodes[var_index] = false;
106 : }
107 :
108 : // Some variables can be sampled directly at their 0 dofs
109 : // - lagrange at nodes on a first order mesh
110 : // Higher order mesh is also fine incidentally because on the 'other' (mid-face for example)
111 : // nodes, the 1st order lagrange variable has 0 dofs, and the second order lagrange
112 : // use the mid-edge nodes as lagrange points (and 0 dofs on extra ones).
113 : // Third order is different (except on edge4, but not used much)
114 : // - anything constant and elemental obviously has the 0-dof value at the centroid (or
115 : // vertex-average). However, higher order elemental, even monomial, do not hold the centroid
116 : // value at dof index 0. For example: pyramid has dof 0 at the center of the base, prism has dof
117 : // 0 on an edge etc
118 4109 : if ((_source_is_nodes[var_index] && fe_type.family == LAGRANGE && fe_type.order <= SECOND) ||
119 4109 : (!_source_is_nodes[var_index] && fe_type.order == CONSTANT))
120 3209 : _use_zero_dof_for_value[var_index] = true;
121 : else
122 48 : _use_zero_dof_for_value[var_index] = false;
123 :
124 : // Check with the source variable that it is not discontinuous at the source
125 3257 : if (_source_is_nodes[var_index])
126 4810 : if (from_var.getContinuity() == DISCONTINUOUS ||
127 2405 : from_var.getContinuity() == SIDE_DISCONTINUOUS)
128 0 : mooseError("Source variable cannot be sampled at nodes as it is discontinuous");
129 :
130 : // Local app does not own any of the target problems
131 3257 : if (_to_problems.empty())
132 3 : break;
133 :
134 : // Check with the target variable that we are not doing awful projections
135 : MooseVariableFieldBase & to_var =
136 6508 : _to_problems[0]->getVariable(0,
137 3254 : _to_var_names[var_index],
138 : Moose::VarKindType::VAR_ANY,
139 : Moose::VarFieldType::VAR_FIELD_ANY);
140 3254 : System & to_sys = to_var.sys().system();
141 3254 : const auto & to_fe_type = to_sys.variable_type(to_var.number());
142 3254 : if (_source_is_nodes[var_index])
143 : {
144 2402 : if (to_fe_type.order == CONSTANT)
145 108 : mooseWarning(
146 : "Transfer is projecting from nearest-nodes to centroids. This is likely causing "
147 : "floating point indetermination in the results because multiple nodes are 'nearest' to "
148 : "a centroid. Please consider using a ProjectionAux to build an elemental source "
149 : "variable (for example constant monomial) before transferring");
150 : }
151 852 : else if (to_var.isNodal())
152 72 : mooseWarning(
153 : "Transfer is projecting from nearest-centroids to nodes. This is likely causing "
154 : "floating point indetermination in the results because multiple centroids are "
155 : "'nearest' to a node. Please consider using a ProjectionAux to build a nodal source "
156 : "variable (for example linear Lagrange) before transferring");
157 : }
158 :
159 : // We need to improve the indexing if we are to allow this
160 3212 : if (!_from_mesh_divisions.empty())
161 384 : for (const auto mesh_div : _from_mesh_divisions)
162 240 : if (mesh_div->getNumDivisions() != _from_mesh_divisions[0]->getNumDivisions())
163 0 : paramError("from_mesh_division",
164 : "This transfer has only been implemented with a uniform number of source mesh "
165 : "divisions across all source applications");
166 3212 : }
167 :
168 : void
169 47732 : MultiAppGeneralFieldNearestLocationTransfer::buildKDTrees(const unsigned int var_index)
170 : {
171 47732 : computeNumSources();
172 47732 : const auto num_apps_per_tree = getNumAppsPerTree();
173 47732 : _local_kdtrees.resize(_num_sources);
174 47732 : _local_points.resize(_num_sources);
175 47732 : _local_values.resize(_num_sources);
176 47732 : unsigned int max_leaf_size = 0;
177 :
178 : // Construct a local KDTree for each source. A source can be a single app or multiple apps
179 : // combined (option for nearest-position / mesh-divisions)
180 98983 : for (const auto i_source : make_range(_num_sources))
181 : {
182 : // Nest a loop on apps in case multiple apps contribute to the same KD-Tree source
183 102542 : for (const auto app_i : make_range(num_apps_per_tree))
184 : {
185 : // Get the current app index
186 51291 : const auto i_from = getAppIndex(i_source, app_i);
187 : // Current position index, if using nearest positions (not used for use_nearest_app)
188 51291 : const auto i_pos = _group_subapps ? i_source : (i_source % getNumDivisions());
189 :
190 : // Get access to the variable and some variable information
191 51291 : FEProblemBase & from_problem = *_from_problems[i_from];
192 51291 : auto & from_mesh = from_problem.mesh(_displaced_source_mesh);
193 : MooseVariableFieldBase & from_var =
194 51291 : from_problem.getVariable(0,
195 51291 : _from_var_names[var_index],
196 : Moose::VarKindType::VAR_ANY,
197 : Moose::VarFieldType::VAR_FIELD_ANY);
198 51291 : System & from_sys = from_var.sys().system();
199 51291 : const unsigned int from_var_num = from_sys.variable_number(getFromVarName(var_index));
200 :
201 : // Build KDTree from nodes and nodal values
202 51291 : if (_source_is_nodes[var_index])
203 : {
204 4387929 : for (const auto & node : from_mesh.getMesh().local_node_ptr_range())
205 : {
206 : // This notably excludes first order Lagrange variables on the mid-nodes
207 : // e.g. the mid-nodes are not added to the KD tree
208 4339143 : if (node->n_dofs(from_sys.number(), from_var_num) < 1)
209 277120 : continue;
210 :
211 4216351 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, from_mesh, node))
212 5152 : continue;
213 :
214 4211199 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, node))
215 34642 : continue;
216 :
217 : // Handle the various source mesh divisions behaviors
218 : // NOTE: This could be more efficient, as instead of rejecting points in the
219 : // wrong division, we could just be adding them to the tree for the right division
220 4176557 : if (!_from_mesh_divisions.empty())
221 : {
222 87080 : const auto tree_division_index = i_source % getNumDivisions();
223 87080 : const auto node_div_index = _from_mesh_divisions[i_from]->divisionIndex(*node);
224 :
225 : // Spatial restriction is always active
226 87080 : if (node_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
227 56552 : continue;
228 : // We fill one tree per division index for matching subapp index or division index. We
229 : // only accept source data from the division index
230 30528 : else if ((_from_mesh_division_behavior ==
231 3328 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
232 3328 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
233 3328 : _from_mesh_division_behavior ==
234 33856 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
235 : tree_division_index != node_div_index)
236 26880 : continue;
237 : }
238 :
239 : // Transformed node is in the reference space, as is the _nearest_positions_obj
240 4093125 : const auto transformed_node = (*_from_transforms[getGlobalSourceAppIndex(i_from)])(*node);
241 :
242 : // Only add to the KDTree nodes that are closest to the 'position'
243 : // When querying values at a target point, the KDTree associated to the closest
244 : // position to the target point is queried
245 : // We do not need to check the positions when using nearest app as we will assume
246 : // (somewhat incorrectly) that all the points in each subapp are closer to that subapp
247 : // than to any other
248 4137421 : if (!_use_nearest_app && _nearest_positions_obj &&
249 44296 : !closestToPosition(i_pos, transformed_node))
250 31102 : continue;
251 :
252 4062023 : _local_points[i_source].push_back(transformed_node);
253 4062023 : const auto dof = node->dof_number(from_sys.number(), from_var_num, 0);
254 4062023 : _local_values[i_source].push_back((*from_sys.solution)(dof));
255 4062023 : if (!_use_zero_dof_for_value[var_index])
256 0 : flagSolutionWarning(
257 : "Nearest-location is not implemented for this source variable type on "
258 : "this mesh. Returning value at dof 0");
259 48786 : }
260 : }
261 : // Build KDTree from centroids and value at centroids
262 : else
263 : {
264 5010 : for (auto & elem : as_range(from_mesh.getMesh().local_elements_begin(),
265 168627 : from_mesh.getMesh().local_elements_end()))
266 : {
267 161112 : if (elem->n_dofs(from_sys.number(), from_var_num) < 1)
268 89928 : continue;
269 :
270 160776 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, elem))
271 3560 : continue;
272 :
273 157216 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, elem))
274 2048 : continue;
275 :
276 : // Handle the various source mesh divisions behaviors
277 155168 : if (!_from_mesh_divisions.empty())
278 : {
279 67424 : const auto tree_division_index = i_source % getNumDivisions();
280 67424 : const auto elem_div_index = _from_mesh_divisions[i_from]->divisionIndex(*elem);
281 :
282 : // Spatial restriction is always active
283 67424 : if (elem_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
284 39744 : continue;
285 : // We fill one tree per division index for matching subapp index or division index. We
286 : // only accept source data from the division index
287 27680 : else if ((_from_mesh_division_behavior ==
288 2880 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
289 2880 : _from_mesh_division_behavior ==
290 30560 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
291 : tree_division_index != elem_div_index)
292 24240 : continue;
293 : }
294 :
295 91184 : const auto vertex_average = elem->vertex_average();
296 : // Transformed element vertex average is in the reference space, as is the
297 : // _nearest_positions_obj
298 : const auto transformed_vertex_average =
299 91184 : (*_from_transforms[getGlobalSourceAppIndex(i_from)])(vertex_average);
300 :
301 : // We do not need to check the positions when using nearest app as we will assume
302 : // (somewhat incorrectly) that all the points in each subapp are closer to that subapp
303 119984 : if (!_use_nearest_app && _nearest_positions_obj &&
304 28800 : !closestToPosition(i_pos, transformed_vertex_average))
305 20000 : continue;
306 :
307 71184 : _local_points[i_source].push_back(transformed_vertex_average);
308 71184 : if (_use_zero_dof_for_value[var_index])
309 : {
310 67744 : auto dof = elem->dof_number(from_sys.number(), from_var_num, 0);
311 67744 : _local_values[i_source].push_back((*from_sys.solution)(dof));
312 : }
313 : else
314 3440 : _local_values[i_source].push_back(
315 3440 : from_sys.point_value(from_var_num, vertex_average, elem));
316 2505 : }
317 : }
318 51291 : max_leaf_size = std::max(max_leaf_size, from_mesh.getMaxLeafSize());
319 : }
320 :
321 : // Make a KDTree from the accumulated points data
322 : std::shared_ptr<KDTree> _kd_tree =
323 51251 : std::make_shared<KDTree>(_local_points[i_source], max_leaf_size);
324 51251 : _local_kdtrees[i_source] = _kd_tree;
325 51251 : }
326 47732 : }
327 :
328 : void
329 74070 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValues(
330 : const unsigned int /*var_index*/,
331 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
332 : std::vector<std::pair<Real, Real>> & outgoing_vals)
333 : {
334 74070 : evaluateInterpValuesNearestNode(incoming_points, outgoing_vals);
335 74070 : }
336 :
337 : void
338 74070 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValuesNearestNode(
339 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
340 : std::vector<std::pair<Real, Real>> & outgoing_vals)
341 : {
342 74070 : dof_id_type i_pt = 0;
343 5964634 : for (const auto & [pt, mesh_div] : incoming_points)
344 : {
345 5890564 : outgoing_vals[i_pt].second = std::numeric_limits<Real>::max();
346 5890564 : bool point_found = false;
347 5890564 : evaluateNearestNodeFromKDTrees(pt, mesh_div, outgoing_vals[i_pt], point_found);
348 5890564 : if (!point_found)
349 4713 : outgoing_vals[i_pt] = {GeneralFieldTransfer::OutOfMeshValue,
350 9426 : GeneralFieldTransfer::OutOfMeshValue};
351 5890564 : i_pt++;
352 : }
353 74070 : }
|
