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 "MultiAppGeneralFieldKDTreeTransferBase.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 : InputParameters
26 17671 : MultiAppGeneralFieldKDTreeTransferBase::validParams()
27 : {
28 17671 : InputParameters params = MultiAppGeneralFieldTransfer::validParams();
29 :
30 53013 : params.addParam<unsigned int>("num_nearest_points",
31 35342 : 1,
32 : "Number of nearest source (from) points will be chosen to "
33 : "construct a value for the target point. All points will be "
34 : "selected from the same origin mesh!");
35 :
36 : // choose whether to include data from multiple apps when performing nearest-position/
37 : // mesh-divisions based transfers
38 35342 : params.addParam<bool>("group_subapps",
39 35342 : false,
40 : "Whether to group source locations and values from all subapps "
41 : "when working with a nearest-position or source mesh-division");
42 :
43 17671 : return params;
44 0 : }
45 :
46 4241 : MultiAppGeneralFieldKDTreeTransferBase::MultiAppGeneralFieldKDTreeTransferBase(
47 4241 : const InputParameters & parameters)
48 : : MultiAppGeneralFieldTransfer(parameters),
49 8482 : _num_nearest_points(getParam<unsigned int>("num_nearest_points")),
50 12723 : _group_subapps(getParam<bool>("group_subapps"))
51 : {
52 4241 : if (_source_app_must_contain_point && _nearest_positions_obj)
53 6 : paramError("use_nearest_position",
54 : "We do not support using both nearest positions matching and checking if target "
55 : "points are within an app domain because the KDTrees for nearest-positions matching "
56 : "are (currently) built with data from multiple applications.");
57 4505 : if (_nearest_positions_obj &&
58 5573 : (isParamValid("from_mesh_divisions") || isParamValid("to_mesh_divisions")))
59 0 : paramError("use_nearest_position", "Cannot use nearest positions with mesh divisions");
60 :
61 : // Parameter checks on grouping subapp values
62 4238 : if (_group_subapps && _from_mesh_divisions.empty() && !_nearest_positions_obj)
63 0 : paramError(
64 : "group_subapps",
65 : "This option is only available for using mesh divisions or nearest positions regions");
66 4361 : else if (_group_subapps &&
67 123 : (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
68 123 : _from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX))
69 0 : paramError("group_subapps",
70 : "Cannot group subapps when considering nearest-location data as we would lose "
71 : "track of the division index of the source locations");
72 4238 : else if (_group_subapps && _use_nearest_app)
73 0 : paramError(
74 : "group_subapps",
75 : "When using the 'nearest child application' data, the source data (positions and values) "
76 : "are grouped on a per-application basis, so it cannot be agglomerated over all child "
77 : "applications.\nNote that the option to use nearest applications for source restrictions, "
78 : "but further split each child application's domain by regions closest to each position "
79 : "(here the the child application's centroid), which could be conceived when "
80 : "'group_subapps' = false, is also not available.");
81 4238 : }
82 :
83 : void
84 4127 : MultiAppGeneralFieldKDTreeTransferBase::initialSetup()
85 : {
86 4127 : MultiAppGeneralFieldTransfer::initialSetup();
87 :
88 : // We need to improve the indexing if we are to allow this
89 4091 : if (!_from_mesh_divisions.empty())
90 444 : for (const auto mesh_div : _from_mesh_divisions)
91 270 : if (mesh_div->getNumDivisions() != _from_mesh_divisions[0]->getNumDivisions())
92 0 : paramError("from_mesh_division",
93 : "This transfer has only been implemented with a uniform number of source mesh "
94 : "divisions across all source applications");
95 4091 : }
96 :
97 : void
98 48536 : MultiAppGeneralFieldKDTreeTransferBase::prepareEvaluationOfInterpValues(
99 : const unsigned int var_index)
100 : {
101 48536 : _local_kdtrees.clear();
102 48536 : _local_points.clear();
103 48536 : _local_values.clear();
104 48536 : buildKDTrees(var_index);
105 48536 : }
106 :
107 : bool
108 40768 : MultiAppGeneralFieldKDTreeTransferBase::inBlocks(const std::set<SubdomainID> & blocks,
109 : const MooseMesh & mesh,
110 : const Elem * elem) const
111 : {
112 : // We need to override the definition of block restriction for an element
113 : // because we have to consider whether each node of an element is adjacent to a block
114 109024 : for (const auto & i_node : make_range(elem->n_nodes()))
115 : {
116 99120 : const auto & node = elem->node_ptr(i_node);
117 99120 : const auto & node_blocks = mesh.getNodeBlockIds(*node);
118 99120 : std::set<SubdomainID> u;
119 99120 : std::set_intersection(blocks.begin(),
120 : blocks.end(),
121 : node_blocks.begin(),
122 : node_blocks.end(),
123 : std::inserter(u, u.begin()));
124 99120 : if (!u.empty())
125 30864 : return true;
126 99120 : }
127 9904 : return false;
128 : }
129 :
130 : void
131 48536 : MultiAppGeneralFieldKDTreeTransferBase::computeNumSources()
132 : {
133 : // Number of source = number of KDTrees.
134 : // Using mesh divisions or nearest-positions, for every app we use 1 tree per division
135 96912 : if (!_from_mesh_divisions.empty() ||
136 48376 : (!_use_nearest_app && _nearest_positions_obj && !_group_subapps))
137 226 : _num_sources = _from_problems.size() * getNumDivisions();
138 : // If we group apps, then we only use one tree per division (nearest-position region)
139 48310 : else if (_nearest_positions_obj && _group_subapps)
140 113 : _num_sources = _nearest_positions_obj->getNumPositions(_fe_problem.getCurrentExecuteOnFlag() ==
141 : EXEC_INITIAL);
142 : // Regular case: 1 KDTree per app
143 : // Also if use_nearest_app = true, the number of problems is better than the number of positions,
144 : // because some of the positions are positions of child applications that are not local
145 : else
146 48197 : _num_sources = _from_problems.size();
147 48536 : }
148 :
149 : unsigned int
150 6539500 : MultiAppGeneralFieldKDTreeTransferBase::getAppIndex(unsigned int kdtree_index,
151 : unsigned int nested_loop_on_app_index) const
152 : {
153 : // Each app is mapped to a single KD Tree
154 6539500 : if (_use_nearest_app)
155 7720 : return kdtree_index;
156 : // We are looping over all the apps that are grouped together
157 6531780 : else if (_group_subapps)
158 25636 : return nested_loop_on_app_index;
159 : // There are num_divisions trees for each app, inner ordering is divisions, so dividing by the
160 : // number of divisions gets us the index of the application
161 : else
162 6506144 : return kdtree_index / getNumDivisions();
163 : }
164 :
165 : unsigned int
166 49730 : MultiAppGeneralFieldKDTreeTransferBase::getNumAppsPerTree() const
167 : {
168 49730 : if (_use_nearest_app)
169 66 : return 1;
170 49664 : else if (_group_subapps)
171 181 : return _from_meshes.size();
172 : else
173 49483 : return 1;
174 : }
175 :
176 : unsigned int
177 13283435 : MultiAppGeneralFieldKDTreeTransferBase::getNumDivisions() const
178 : {
179 : // This is not used currently, but conceptually it is better to only divide the domain with the
180 : // local of local applications rather than the global number of positions (# global applications
181 : // here)
182 13283435 : if (_use_nearest_app)
183 7720 : return _from_meshes.size();
184 : // Each nearest-position region is a division
185 13275715 : else if (_nearest_positions_obj && !_group_subapps)
186 46314 : return _nearest_positions_obj->getNumPositions(_fe_problem.getCurrentExecuteOnFlag() ==
187 46314 : EXEC_INITIAL);
188 : // Assume all mesh divisions (on each sub-app) has the same number of divisions. This is checked
189 13229401 : else if (!_from_mesh_divisions.empty())
190 790554 : return _from_mesh_divisions[0]->getNumDivisions();
191 : // Grouping subapps or no special mode, we do not subdivide
192 : else
193 12438847 : return 1;
194 : }
195 :
196 : Point
197 15395 : MultiAppGeneralFieldKDTreeTransferBase::getPointInSourceKDTreeFrame(unsigned int i_source,
198 : const Point & pt) const
199 : {
200 15395 : if (_nearest_positions_obj && _group_subapps)
201 : // KD-trees are built in global coords when grouping subapps with nearest-positions
202 6075 : return pt;
203 : else
204 27960 : return getPointInSourceAppFrame(pt, i_source, "KD-tree neighbor search");
205 : }
206 :
207 : bool
208 6364238 : MultiAppGeneralFieldKDTreeTransferBase::checkRestrictionsForSource(const Point & pt,
209 : const unsigned int mesh_div,
210 : const unsigned int i_from) const
211 : {
212 : // Only use the KDTree from the closest position if in "nearest-position" mode
213 6364238 : if (_nearest_positions_obj)
214 : {
215 : // See computeNumSources for the number of sources. i_from is the index in the source loop
216 : // i_from is local if looping on _from_problems as sources, positions are indexed globally
217 : // i_from is already indexing in positions if using group_subapps
218 70889 : auto position_index = i_from; // if _group_subapps
219 70889 : if (_use_nearest_app)
220 9949 : position_index = getGlobalSourceAppIndex(i_from);
221 60940 : else if (!_group_subapps)
222 19760 : position_index = i_from % getNumDivisions();
223 :
224 : // NOTE: if two positions are equi-distant to the point, this will chose one
225 : // This problem is detected if using search_value_conflicts in this call
226 70889 : if (!closestToPosition(position_index, pt))
227 39976 : return false;
228 : }
229 :
230 : // Application index depends on which source/grouping mode we are using
231 6324262 : const unsigned int app_index = getAppIndex(i_from, i_from / getNumDivisions());
232 :
233 : // Check mesh restriction before anything
234 6324262 : if (_source_app_must_contain_point)
235 : {
236 0 : if (_nearest_positions_obj)
237 0 : mooseError("Nearest-positions + source_app_must_contain_point not implemented");
238 : // Transform the point to place it in the local coordinate system
239 0 : const auto local_pt = getPointInSourceAppFrame(pt, app_index, "Source mesh containment check");
240 0 : if (!inMesh(_from_point_locators[app_index].get(), local_pt))
241 0 : return false;
242 : }
243 :
244 : // Check the mesh division. We have handled the restriction of the source locations when
245 : // building the nearest-neighbor trees. We only need to check that we meet the required
246 : // source division index.
247 6324262 : if (!_from_mesh_divisions.empty())
248 : {
249 : mooseAssert(mesh_div != MooseMeshDivision::INVALID_DIVISION_INDEX,
250 : "We should not be receiving point requests with an invalid "
251 : "source mesh division index");
252 208800 : const unsigned int kd_div_index = i_from % getNumDivisions();
253 :
254 : // If matching source mesh divisions to target apps, we check that the index of the target
255 : // application, which was passed in the point request, is equal to the current mesh division
256 208800 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
257 : mesh_div != kd_div_index)
258 25917 : return false;
259 : // If matching source mesh divisions to target mesh divisions, we check that the index of the
260 : // target mesh division, which was passed in the point request, is equal to the current mesh
261 : // division
262 182883 : else if ((_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
263 182883 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX) &&
264 : mesh_div != kd_div_index)
265 133152 : return false;
266 : }
267 :
268 : // If matching target apps to source mesh divisions, we check that the global index of the
269 : // application is equal to the target mesh division index, which was passed in the point request
270 6169165 : if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
271 3972 : mesh_div != getGlobalSourceAppIndex(app_index))
272 2916 : return false;
273 :
274 6162277 : return true;
275 : }
276 :
277 : void
278 5895129 : MultiAppGeneralFieldKDTreeTransferBase::evaluateNearestNodeFromKDTrees(
279 : const Point & pt,
280 : unsigned int source_index,
281 : std::pair<Real, Real> & outgoing_val,
282 : bool & point_found)
283 : {
284 : // First pass: find the nearest neighbor value across all local sources
285 12072089 : for (const auto i_from : make_range(_num_sources))
286 : {
287 6176960 : if (!checkRestrictionsForSource(pt, source_index, i_from))
288 176101 : continue;
289 :
290 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the searches
291 12001718 : std::vector<std::size_t> return_index(_num_nearest_points);
292 6000859 : std::vector<Real> return_dist_sqr(_num_nearest_points);
293 :
294 : // KD-Tree can be empty if no points are within block/boundary/bounding box restrictions
295 6000859 : if (_local_kdtrees[i_from]->numberCandidatePoints())
296 : {
297 5993033 : point_found = true;
298 : // Note that we do not need to use the transformed_pt (in the source app frame)
299 : // because the KDTree has been created in the reference frame
300 5993033 : _local_kdtrees[i_from]->neighborSearch(
301 : pt, _num_nearest_points, return_index, return_dist_sqr);
302 5993033 : Real val_sum = 0, dist_sum = 0;
303 11986066 : for (const auto index : return_index)
304 : {
305 5993033 : val_sum += _local_values[i_from][index];
306 5993033 : dist_sum += (_local_points[i_from][index] - pt).norm();
307 : }
308 :
309 : // If the new value found is closer than for other sources, use it
310 5993033 : const auto new_distance = dist_sum / return_dist_sqr.size();
311 5993033 : if (new_distance < outgoing_val.second)
312 5926971 : outgoing_val = {val_sum / return_index.size(), new_distance};
313 : }
314 6000859 : }
315 :
316 : // Second pass: search-value-conflict detection
317 5895129 : if (point_found && _search_value_conflicts)
318 : {
319 15032 : unsigned int num_equidistant_problems = 0;
320 :
321 49939 : for (const auto i_from : make_range(_num_sources))
322 : {
323 34907 : if (!checkRestrictionsForSource(pt, source_index, i_from))
324 18292 : continue;
325 :
326 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the searches
327 33230 : std::vector<std::size_t> return_index(_num_nearest_points + 1);
328 16615 : std::vector<Real> return_dist_sqr(_num_nearest_points + 1);
329 :
330 16615 : const auto app_index = getAppIndex(i_from, i_from / getNumDivisions());
331 16615 : const auto num_search = _num_nearest_points + 1;
332 :
333 16615 : if (_local_kdtrees[i_from]->numberCandidatePoints())
334 : {
335 15395 : _local_kdtrees[i_from]->neighborSearch(pt, num_search, return_index, return_dist_sqr);
336 15395 : auto num_found = return_dist_sqr.size();
337 :
338 : // Local coordinates only accessible when not using nearest-position,
339 : // as we did not keep the index of the source app, only the position index
340 15395 : const Point local_pt = getPointInSourceKDTreeFrame(app_index, pt);
341 :
342 20811 : if (!_nearest_positions_obj &&
343 5416 : _from_transforms[getGlobalSourceAppIndex(app_index)]->hasCoordinateSystemTypeChange())
344 0 : if (!_skip_coordinate_collapsing)
345 0 : mooseInfo("Search value conflict cannot find the origin point due to the "
346 : "non-uniqueness of the coordinate collapsing reverse mapping");
347 :
348 : // Look for too many equidistant nodes within a problem. First zip then sort by distance
349 15395 : std::vector<std::pair<Real, std::size_t>> zipped_nearest_points;
350 46040 : for (const auto i : make_range(num_found))
351 30645 : zipped_nearest_points.push_back(std::make_pair(return_dist_sqr[i], return_index[i]));
352 15395 : std::sort(zipped_nearest_points.begin(), zipped_nearest_points.end());
353 :
354 : // If two furthest are equally far from target point, then we have an indetermination in
355 : // what is sent in this communication round from this process. However, it may not
356 : // materialize to an actual conflict, as values sent from another process for the
357 : // desired target point could be closer (nearest). There is no way to know at this point
358 : // in the communication that a closer value exists somewhere else
359 30645 : if (num_found > 1 && num_found == num_search &&
360 15250 : MooseUtils::absoluteFuzzyEqual(zipped_nearest_points[num_found - 1].first,
361 15250 : zipped_nearest_points[num_found - 2].first))
362 : {
363 288 : if (_nearest_positions_obj)
364 12 : registerConflict(app_index, 0, pt, outgoing_val.second, true);
365 : else
366 276 : registerConflict(app_index, 0, local_pt, outgoing_val.second, true);
367 : }
368 :
369 : // Recompute the distance for this problem. If it matches the cached value more than
370 : // once it means multiple problems provide equidistant values for this point
371 15395 : Real dist_sum = 0;
372 30790 : for (const auto i : make_range(num_search - 1))
373 : {
374 15395 : auto index = zipped_nearest_points[i].second;
375 15395 : dist_sum += (_local_points[i_from][index] - pt).norm();
376 : }
377 :
378 : // Compare to the selected value found after looking at all the problems
379 15395 : if (MooseUtils::absoluteFuzzyEqual(dist_sum / return_dist_sqr.size(), outgoing_val.second))
380 : {
381 211 : num_equidistant_problems++;
382 211 : if (num_equidistant_problems > 1)
383 : {
384 0 : if (_nearest_positions_obj)
385 0 : registerConflict(app_index, 0, pt, outgoing_val.second, true);
386 : else
387 0 : registerConflict(app_index, 0, local_pt, outgoing_val.second, true);
388 : }
389 : }
390 15395 : }
391 16615 : }
392 : }
393 5895129 : }
|
