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 "RadialAverage.h"
11 : #include "ThreadedRadialAverageLoop.h"
12 : #include "FEProblemBase.h"
13 :
14 : #include "libmesh/nanoflann.hpp"
15 : #include "libmesh/parallel_algebra.h"
16 : #include "libmesh/mesh_tools.h"
17 : #include "libmesh/int_range.h"
18 :
19 : #include <list>
20 : #include <iterator>
21 : #include <algorithm>
22 :
23 : // Make newer nanoflann API compatible with older nanoflann versions
24 : #if NANOFLANN_VERSION < 0x150
25 : namespace nanoflann
26 : {
27 : typedef SearchParams SearchParameters;
28 :
29 : template <typename T, typename U>
30 : using ResultItem = std::pair<T, U>;
31 : }
32 : #endif
33 :
34 : registerMooseObject("MooseApp", RadialAverage);
35 :
36 : // specialization for PointListAdaptor<RadialAverage::QPData>
37 : template <>
38 : const Point &
39 28469778 : PointListAdaptor<RadialAverage::QPData>::getPoint(const RadialAverage::QPData & item) const
40 : {
41 28469778 : return item._q_point;
42 : }
43 :
44 : InputParameters
45 14328 : RadialAverage::validParams()
46 : {
47 14328 : InputParameters params = ElementUserObject::validParams();
48 14328 : params.addClassDescription("Perform a radial average of a material property");
49 14328 : params.addRequiredParam<MaterialPropertyName>("prop_name",
50 : "Name of the material property to average");
51 14328 : MooseEnum weights_type("constant linear cosine", "linear");
52 14328 : params.addRequiredParam<MooseEnum>("weights", weights_type, "Distance based weight function");
53 14328 : params.addRequiredParam<Real>("radius", "Cut-off radius for the averaging");
54 42984 : params.addRangeCheckedParam<Real>(
55 : "padding",
56 28656 : 0.0,
57 : "padding >= 0",
58 : "Padding for communication. This gets added to the radius when determining which QPs to send "
59 : "to other processors. Increase this gradually if inconsistent parallel results occur.");
60 :
61 : // we run this object once at the beginning of the timestep by default
62 14328 : params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
63 :
64 : // make sure we always have geometric point neighbors ghosted
65 14328 : params.addRelationshipManager("ElementPointNeighborLayers",
66 : Moose::RelationshipManagerType::GEOMETRIC |
67 : Moose::RelationshipManagerType::ALGEBRAIC);
68 14328 : params.addParamNamesToGroup("padding", "Advanced");
69 28656 : return params;
70 14328 : }
71 :
72 33 : RadialAverage::RadialAverage(const InputParameters & parameters)
73 : : ElementUserObject(parameters),
74 33 : _weights_type(getParam<MooseEnum>("weights").getEnum<WeightsType>()),
75 33 : _prop(getMaterialProperty<Real>("prop_name")),
76 33 : _radius(getParam<Real>("radius")),
77 33 : _padding(getParam<Real>("padding")),
78 33 : _update_communication_lists(false),
79 33 : _my_pid(processor_id()),
80 33 : _perf_meshchanged(registerTimedSection("meshChanged", 3)),
81 33 : _perf_updatelists(registerTimedSection("updateCommunicationLists", 3)),
82 99 : _perf_finalize(registerTimedSection("finalize", 2))
83 : {
84 33 : }
85 :
86 : void
87 33 : RadialAverage::initialSetup()
88 : {
89 : // set up processor boundary node list
90 33 : meshChanged();
91 33 : }
92 :
93 : void
94 132 : RadialAverage::initialize()
95 : {
96 132 : _qp_data.clear();
97 132 : }
98 :
99 : void
100 8400 : RadialAverage::execute()
101 : {
102 8400 : auto id = _current_elem->id();
103 :
104 : // collect all QP data
105 84000 : for (const auto qp : make_range(_qrule->n_points()))
106 75600 : _qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _prop[qp]);
107 :
108 : // make sure the result map entry for the current element is sized correctly
109 8400 : auto i = _average.find(id);
110 8400 : if (i == _average.end())
111 2100 : _average.insert(std::make_pair(id, std::vector<Real>(_qrule->n_points())));
112 : else
113 6300 : i->second.resize(_qrule->n_points());
114 8400 : }
115 :
116 : void
117 120 : RadialAverage::finalize()
118 : {
119 120 : TIME_SECTION(_perf_finalize);
120 :
121 : // the first chunk of data is always the local data - remember its size
122 120 : unsigned int local_size = _qp_data.size();
123 :
124 : // communicate the qp data list if n_proc > 1
125 120 : if (_app.n_processors() > 1)
126 : {
127 : // !!!!!!!!!!!
128 : // !!CAREFUL!! Is it guaranteed that _qp_data is in the same order if the mesh has not changed?
129 : // According to @friedmud it is not guaranteed if threads are used
130 : // !!!!!!!!!!!
131 :
132 : // update after mesh changes and/or if a displaced problem exists
133 126 : if (_update_communication_lists || _fe_problem.getDisplacedProblem() ||
134 54 : libMesh::n_threads() > 1)
135 18 : updateCommunicationLists();
136 :
137 : // data structures for sparse point to point communication
138 72 : std::vector<std::vector<QPData>> send(_candidate_procs.size());
139 72 : std::vector<Parallel::Request> send_requests(_candidate_procs.size());
140 72 : Parallel::MessageTag send_tag = _communicator.get_unique_tag(4711);
141 72 : std::vector<QPData> receive;
142 :
143 72 : const auto item_type = TIMPI::StandardType<QPData>(&(_qp_data[0]));
144 :
145 : // fill buffer and send structures
146 144 : for (const auto i : index_range(_candidate_procs))
147 : {
148 72 : const auto pid = _candidate_procs[i];
149 72 : const auto & list = _communication_lists[pid];
150 :
151 : // fill send buffer for transfer to pid
152 72 : send[i].reserve(list.size());
153 9768 : for (const auto & item : list)
154 9696 : send[i].push_back(_qp_data[item]);
155 :
156 : // issue non-blocking send
157 72 : _communicator.send(pid, send[i], send_requests[i], send_tag);
158 : }
159 :
160 : // receive messages - we assume that we receive as many messages as we send!
161 : // bounding box overlapp is transitive, but data exhange between overlapping procs could still
162 : // be unidirectional!
163 144 : for (const auto i : index_range(_candidate_procs))
164 : {
165 72 : libmesh_ignore(i);
166 :
167 : // inspect incoming message
168 72 : Parallel::Status status(_communicator.probe(Parallel::any_source, send_tag));
169 72 : const auto source_pid = TIMPI::cast_int<processor_id_type>(status.source());
170 72 : const auto message_size = status.size(item_type);
171 :
172 : // resize receive buffer accordingly and receive data
173 72 : receive.resize(message_size);
174 72 : _communicator.receive(source_pid, receive, send_tag);
175 :
176 : // append communicated data
177 72 : _qp_data.insert(_qp_data.end(), receive.begin(), receive.end());
178 : }
179 :
180 : // wait until all send requests are at least buffered and we can destroy
181 : // the send buffers by going out of scope
182 72 : Parallel::wait(send_requests);
183 72 : }
184 :
185 : // build KD-Tree using data we just received
186 120 : const unsigned int max_leaf_size = 20; // slightly affects runtime
187 120 : auto point_list = PointListAdaptor<QPData>(_qp_data.begin(), _qp_data.end());
188 120 : _kd_tree = std::make_unique<KDTreeType>(
189 120 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
190 :
191 : mooseAssert(_kd_tree != nullptr, "KDTree was not properly initialized.");
192 120 : _kd_tree->buildIndex();
193 :
194 : // build thread loop functor
195 120 : ThreadedRadialAverageLoop rgcl(*this);
196 :
197 : // run threads
198 120 : auto local_range_begin = _qp_data.begin();
199 120 : auto local_range_end = local_range_begin;
200 120 : std::advance(local_range_end, local_size);
201 120 : Threads::parallel_reduce(QPDataRange(local_range_begin, local_range_end), rgcl);
202 120 : }
203 :
204 : void
205 12 : RadialAverage::threadJoin(const UserObject & y)
206 : {
207 12 : const auto & uo = static_cast<const RadialAverage &>(y);
208 12 : _qp_data.insert(_qp_data.begin(), uo._qp_data.begin(), uo._qp_data.end());
209 12 : _average.insert(uo._average.begin(), uo._average.end());
210 12 : }
211 :
212 : void
213 33 : RadialAverage::meshChanged()
214 : {
215 33 : TIME_SECTION(_perf_meshchanged);
216 :
217 : // get underlying libMesh mesh
218 33 : auto & mesh = _mesh.getMesh();
219 :
220 : // Build a new node to element map
221 33 : _nodes_to_elem_map.clear();
222 33 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), _nodes_to_elem_map);
223 :
224 : // clear procesor boundary nodes set
225 33 : _boundary_nodes.clear();
226 :
227 : //
228 : // iterate over active local elements and store all processor boundary node locations
229 : //
230 33 : const auto end = mesh.active_local_elements_end();
231 4833 : for (auto it = mesh.active_local_elements_begin(); it != end; ++it)
232 : // find faces at processor boundaries
233 12000 : for (const auto s : make_range((*it)->n_sides()))
234 : {
235 9600 : const auto * neighbor = (*it)->neighbor_ptr(s);
236 9600 : if (neighbor && neighbor->processor_id() != _my_pid)
237 : // add all nodes on the processor boundary
238 1020 : for (const auto n : make_range((*it)->n_nodes()))
239 816 : if ((*it)->is_node_on_side(n, s))
240 408 : _boundary_nodes.insert((*it)->node_ref(n));
241 :
242 : // request communication list update
243 9600 : _update_communication_lists = true;
244 33 : }
245 33 : }
246 :
247 : void
248 18 : RadialAverage::updateCommunicationLists()
249 : {
250 18 : TIME_SECTION(_perf_updatelists);
251 :
252 : // clear communication lists
253 18 : _communication_lists.clear();
254 18 : _communication_lists.resize(n_processors());
255 :
256 : // build KD-Tree using local qpoint data
257 18 : const unsigned int max_leaf_size = 20; // slightly affects runtime
258 18 : auto point_list = PointListAdaptor<QPData>(_qp_data.begin(), _qp_data.end());
259 : auto kd_tree = std::make_unique<KDTreeType>(
260 18 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
261 : mooseAssert(kd_tree != nullptr, "KDTree was not properly initialized.");
262 18 : kd_tree->buildIndex();
263 :
264 18 : std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
265 18 : nanoflann::SearchParameters search_params;
266 :
267 : // iterate over all boundary nodes and collect all boundary-near data points
268 18 : _boundary_data_indices.clear();
269 240 : for (const auto & bn : _boundary_nodes)
270 : {
271 222 : ret_matches.clear();
272 444 : kd_tree->radiusSearch(
273 222 : &(bn(0)), Utility::pow<2>(_radius + _padding), ret_matches, search_params);
274 6234 : for (auto & match : ret_matches)
275 6012 : _boundary_data_indices.insert(match.first);
276 : }
277 :
278 : // gather all processor bounding boxes (communicate as pairs)
279 18 : std::vector<std::pair<Point, Point>> pps(n_processors());
280 18 : const auto mybb = _mesh.getInflatedProcessorBoundingBox(0);
281 18 : std::pair<Point, Point> mypp = mybb;
282 18 : _communicator.allgather(mypp, pps);
283 :
284 : // inflate all processor bounding boxes by radius (no padding)
285 18 : const auto rpoint = Point(1, 1, 1) * _radius;
286 18 : std::vector<BoundingBox> bbs;
287 54 : for (const auto & pp : pps)
288 36 : bbs.emplace_back(pp.first - rpoint, pp.second + rpoint);
289 :
290 : // get candidate processors (overlapping bounding boxes)
291 18 : _candidate_procs.clear();
292 54 : for (const auto pid : index_range(bbs))
293 36 : if (pid != _my_pid && bbs[pid].intersects(mypp))
294 18 : _candidate_procs.push_back(pid);
295 :
296 : // go over all boundary data items and send them to the proc they overlap with
297 2442 : for (const auto i : _boundary_data_indices)
298 4848 : for (const auto pid : _candidate_procs)
299 2424 : if (bbs[pid].contains_point(_qp_data[i]._q_point))
300 2424 : _communication_lists[pid].insert(i);
301 :
302 : // done
303 18 : _update_communication_lists = false;
304 18 : }
305 :
306 : namespace TIMPI
307 : {
308 :
309 216 : StandardType<RadialAverage::QPData>::StandardType(const RadialAverage::QPData * example)
310 : {
311 : // We need an example for MPI_Address to use
312 216 : static const RadialAverage::QPData p;
313 216 : if (!example)
314 0 : example = &p;
315 :
316 : #ifdef LIBMESH_HAVE_MPI
317 :
318 : // Get the sub-data-types, and make sure they live long enough
319 : // to construct the derived type
320 216 : StandardType<Point> d1(&example->_q_point);
321 216 : StandardType<dof_id_type> d2(&example->_elem_id);
322 216 : StandardType<short> d3(&example->_qp);
323 216 : StandardType<Real> d4(&example->_volume);
324 216 : StandardType<Real> d5(&example->_value);
325 :
326 : MPI_Datatype types[] = {
327 216 : (data_type)d1, (data_type)d2, (data_type)d3, (data_type)d4, (data_type)d5};
328 216 : int blocklengths[] = {1, 1, 1, 1, 1};
329 : MPI_Aint displs[5], start;
330 :
331 216 : libmesh_call_mpi(MPI_Get_address(const_cast<RadialAverage::QPData *>(example), &start));
332 216 : libmesh_call_mpi(MPI_Get_address(const_cast<Point *>(&example->_q_point), &displs[0]));
333 216 : libmesh_call_mpi(MPI_Get_address(const_cast<dof_id_type *>(&example->_elem_id), &displs[1]));
334 216 : libmesh_call_mpi(MPI_Get_address(const_cast<short *>(&example->_qp), &displs[2]));
335 216 : libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_volume), &displs[3]));
336 216 : libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_value), &displs[4]));
337 :
338 1296 : for (std::size_t i = 0; i < 5; ++i)
339 1080 : displs[i] -= start;
340 :
341 : // create a prototype structure
342 : MPI_Datatype tmptype;
343 216 : libmesh_call_mpi(MPI_Type_create_struct(5, blocklengths, displs, types, &tmptype));
344 216 : libmesh_call_mpi(MPI_Type_commit(&tmptype));
345 :
346 : // resize the structure type to account for padding, if any
347 216 : libmesh_call_mpi(MPI_Type_create_resized(tmptype, 0, sizeof(RadialAverage::QPData), &_datatype));
348 216 : libmesh_call_mpi(MPI_Type_free(&tmptype));
349 :
350 216 : this->commit();
351 :
352 : #endif // LIBMESH_HAVE_MPI
353 216 : }
354 :
355 0 : StandardType<RadialAverage::QPData>::StandardType(const StandardType<RadialAverage::QPData> & t)
356 0 : : DataType(t._datatype)
357 : {
358 : #ifdef LIBMESH_HAVE_MPI
359 0 : libmesh_call_mpi(MPI_Type_dup(t._datatype, &_datatype));
360 : #endif
361 0 : }
362 :
363 : } // namespace TIMPI
|
