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 "NodalPatchRecoveryBase.h"
11 : #include "MathUtils.h"
12 :
13 : #include <Eigen/Dense>
14 :
15 : // TIMPI includes
16 : #include "timpi/communicator.h"
17 : #include "timpi/parallel_sync.h"
18 : #include "libmesh/parallel_eigen.h"
19 : #include <iomanip>
20 :
21 : // Remove duplicates and sort entries from a vector of element IDs.
22 : // This is necessary because user input may contain repeated element IDs,
23 : // which is problematic when using these IDs as keys.
24 : // In Patch Recovery, we also want to avoid duplicate elements contributing
25 : // multiple times to the Ae and be.
26 : static std::vector<dof_id_type>
27 7096 : removeDuplicateEntries(const std::vector<dof_id_type> & ids)
28 : {
29 7096 : std::vector<dof_id_type> key = ids;
30 7096 : std::sort(key.begin(), key.end());
31 7096 : key.erase(std::unique(key.begin(), key.end()), key.end());
32 7096 : return key;
33 0 : }
34 :
35 : InputParameters
36 29267 : NodalPatchRecoveryBase::validParams()
37 : {
38 29267 : InputParameters params = ElementUserObject::validParams();
39 :
40 117068 : MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
41 87801 : params.addRequiredParam<MooseEnum>(
42 : "patch_polynomial_order",
43 : orders,
44 : "Polynomial order used in least squares fitting of material property "
45 : "over the local patch of elements connected to a given node");
46 :
47 87801 : params.addRelationshipManager("ElementSideNeighborLayers",
48 : Moose::RelationshipManagerType::ALGEBRAIC,
49 30314 : [](const InputParameters &, InputParameters & rm_params)
50 2094 : { rm_params.set<unsigned short>("layers") = 2; });
51 :
52 87801 : params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
53 :
54 58534 : return params;
55 29267 : }
56 :
57 380 : NodalPatchRecoveryBase::NodalPatchRecoveryBase(const InputParameters & parameters)
58 : : ElementUserObject(parameters),
59 380 : _qp(0),
60 380 : _patch_polynomial_order(
61 380 : static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
62 380 : _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
63 380 : _q(_multi_index.size()),
64 380 : _distributed_mesh(_mesh.isDistributedMesh()),
65 1520 : _proc_ids(n_processors())
66 : {
67 380 : std::iota(_proc_ids.begin(), _proc_ids.end(), 0);
68 380 : }
69 :
70 : Real
71 5566 : NodalPatchRecoveryBase::nodalPatchRecovery(const Point & x,
72 : const std::vector<dof_id_type> & elem_ids) const
73 : {
74 5566 : const RealEigenVector coef = getCoefficients(elem_ids); // const version
75 : // Compute the fitted nodal value
76 5566 : RealEigenVector p = evaluateBasisFunctions(x);
77 11132 : return p.dot(coef);
78 5566 : }
79 :
80 : const RealEigenVector
81 6331 : NodalPatchRecoveryBase::getCoefficients(const std::vector<dof_id_type> & elem_ids) const
82 : {
83 6331 : auto elem_ids_reduced = removeDuplicateEntries(elem_ids);
84 :
85 6331 : RealEigenVector coef = RealEigenVector::Zero(_q);
86 : // Before we go, check if we have enough sample points for solving the least square fitting
87 6331 : if (_q_point.size() * elem_ids_reduced.size() < _q)
88 0 : mooseError("There are not enough sample points to recover the nodal value, try reducing the "
89 : "polynomial order or using a higher-order quadrature scheme.");
90 :
91 : // Assemble the least squares problem over the patch
92 6331 : RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
93 6331 : RealEigenVector b = RealEigenVector::Zero(_q);
94 31647 : for (auto elem_id : elem_ids_reduced)
95 : {
96 25316 : const auto elem = _mesh.elemPtr(elem_id);
97 25316 : if (elem /*prevent segmentation fault in distributed mesh*/)
98 25003 : if (!hasBlocks(elem->subdomain_id()))
99 0 : mooseError("Element with id = ",
100 : elem_id,
101 : " is not in the block. "
102 : "Please use nodalPatchRecovery with elements in the block only.");
103 :
104 25316 : if (_Ae.find(elem_id) == _Ae.end())
105 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _Ae.");
106 25316 : if (_be.find(elem_id) == _be.end())
107 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _be.");
108 :
109 25316 : A += libmesh_map_find(_Ae, elem_id);
110 25316 : b += libmesh_map_find(_be, elem_id);
111 : }
112 :
113 : // Solve the least squares fitting
114 6331 : coef = A.completeOrthogonalDecomposition().solve(b);
115 :
116 12662 : return coef;
117 6331 : }
118 :
119 : const RealEigenVector
120 765 : NodalPatchRecoveryBase::getCachedCoefficients(const std::vector<dof_id_type> & elem_ids)
121 : {
122 : // Check cache
123 765 : auto key = removeDuplicateEntries(elem_ids);
124 :
125 765 : if (key == _cached_elem_ids)
126 0 : return _cached_coef;
127 : else
128 : {
129 765 : const auto coef = getCoefficients(key); // const version
130 :
131 765 : _cached_elem_ids = key; // Update the cached element IDs
132 765 : _cached_coef = coef; // Update the cached coefficients
133 :
134 765 : return coef;
135 765 : }
136 765 : }
137 :
138 : RealEigenVector
139 133421 : NodalPatchRecoveryBase::evaluateBasisFunctions(const Point & q_point) const
140 : {
141 133421 : RealEigenVector p(_q);
142 : Real polynomial;
143 789940 : for (unsigned int r = 0; r < _multi_index.size(); r++)
144 : {
145 656519 : polynomial = 1.0;
146 : mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
147 2335637 : for (unsigned int c = 0; c < _multi_index[r].size(); c++)
148 2421864 : for (unsigned int p = 0; p < _multi_index[r][c]; p++)
149 742746 : polynomial *= q_point(c);
150 656519 : p(r) = polynomial;
151 : }
152 133421 : return p;
153 0 : }
154 :
155 : void
156 1312 : NodalPatchRecoveryBase::initialize()
157 : {
158 : // Clear cached data to ensure coefficients are correctly recomputed in the next patch recovery
159 : // iteration. _Ae and _be must also be reset, as the associated patch elements may differ in the
160 : // upcoming iteration.
161 :
162 1312 : _cached_elem_ids.clear();
163 1312 : _cached_coef = RealEigenVector::Zero(_q);
164 1312 : _Ae.clear();
165 1312 : _be.clear();
166 1312 : }
167 :
168 : void
169 26449 : NodalPatchRecoveryBase::execute()
170 : {
171 26449 : RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
172 26449 : RealEigenVector be = RealEigenVector::Zero(_q);
173 154304 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
174 : {
175 127855 : RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
176 127855 : Ae += p * p.transpose();
177 127855 : be += computeValue() * p;
178 127855 : }
179 :
180 26449 : dof_id_type elem_id = _current_elem->id();
181 :
182 26449 : _Ae[elem_id] = Ae;
183 26449 : _be[elem_id] = be;
184 26449 : }
185 :
186 : void
187 99 : NodalPatchRecoveryBase::threadJoin(const UserObject & uo)
188 : {
189 99 : const auto & npr = static_cast<const NodalPatchRecoveryBase &>(uo);
190 99 : _Ae.insert(npr._Ae.begin(), npr._Ae.end());
191 99 : _be.insert(npr._be.begin(), npr._be.end());
192 99 : }
193 :
194 : void
195 1213 : NodalPatchRecoveryBase::finalize()
196 : {
197 : // When calling nodalPatchRecovery, we may need to know _Ae and _be on algebraically ghosted
198 : // elements. However, this userobject is only run on local elements, so we need to query those
199 : // information from other processors in this finalize() method.
200 1213 : sync();
201 1213 : }
202 :
203 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
204 765 : NodalPatchRecoveryBase::gatherRequestList(const std::vector<dof_id_type> & specific_elems)
205 : {
206 765 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
207 :
208 : typedef std::pair<processor_id_type, dof_id_type> PidElemPair;
209 765 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
210 :
211 7129 : for (const auto & entry : specific_elems)
212 : {
213 6364 : const auto * elem = _mesh.elemPtr(entry);
214 6364 : if (_distributed_mesh)
215 : {
216 1024 : if (!elem)
217 313 : continue; // Prevent segmentation fault in distributed mesh
218 :
219 711 : if (hasBlocks(elem->subdomain_id()))
220 2133 : for (processor_id_type pid = 0; pid < n_processors(); ++pid)
221 1422 : if (pid != processor_id())
222 711 : push_data[pid].push_back(std::make_pair(elem->processor_id(), elem->id()));
223 : }
224 : else
225 : // Add to query_ids if this element's data is on a different processor
226 5340 : addToQuery(elem, query_ids);
227 : }
228 :
229 765 : if (_distributed_mesh)
230 : {
231 : auto push_receiver =
232 92 : [&](const processor_id_type, const std::vector<PidElemPair> & received_data)
233 : {
234 803 : for (const auto & [pid, id] : received_data)
235 711 : query_ids[pid].push_back(id);
236 92 : };
237 :
238 124 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
239 : }
240 :
241 1530 : return query_ids;
242 765 : }
243 :
244 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
245 1213 : NodalPatchRecoveryBase::gatherRequestList()
246 : {
247 1213 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
248 :
249 : typedef std::pair<processor_id_type, dof_id_type> PidElemPair;
250 1213 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
251 :
252 71816 : for (const auto & elem : _fe_problem.getEvaluableElementRange())
253 70603 : addToQuery(elem, query_ids);
254 :
255 2426 : return query_ids;
256 1213 : }
257 :
258 : void
259 1213 : NodalPatchRecoveryBase::sync()
260 : {
261 1213 : const auto query_ids = gatherRequestList();
262 1213 : syncHelper(query_ids);
263 1213 : }
264 :
265 : void
266 765 : NodalPatchRecoveryBase::sync(const std::vector<dof_id_type> & specific_elems)
267 : {
268 765 : const auto query_ids = gatherRequestList(specific_elems);
269 765 : syncHelper(query_ids);
270 765 : }
271 :
272 : void
273 1978 : NodalPatchRecoveryBase::syncHelper(
274 : const std::unordered_map<processor_id_type, std::vector<dof_id_type>> & query_ids)
275 : {
276 : typedef std::pair<RealEigenMatrix, RealEigenVector> AbPair;
277 :
278 : // Answer queries received from other processors
279 815 : auto gather_data = [this](const processor_id_type /*pid*/,
280 : const std::vector<dof_id_type> & elem_ids,
281 : std::vector<AbPair> & ab_pairs)
282 : {
283 5990 : for (const auto & elem_id : elem_ids)
284 5175 : ab_pairs.emplace_back(libmesh_map_find(_Ae, elem_id), libmesh_map_find(_be, elem_id));
285 815 : };
286 :
287 : // Gather answers received from other processors
288 815 : auto act_on_data = [this](const processor_id_type /*pid*/,
289 : const std::vector<dof_id_type> & elem_ids,
290 : const std::vector<AbPair> & ab_pairs)
291 : {
292 5990 : for (const auto i : index_range(elem_ids))
293 : {
294 5175 : const auto elem_id = elem_ids[i];
295 5175 : const auto & [Ae, be] = ab_pairs[i];
296 5175 : _Ae[elem_id] = Ae;
297 5175 : _be[elem_id] = be;
298 : }
299 815 : };
300 :
301 : // the send and receive are called inside the pull_parallel_vector_data function
302 1978 : libMesh::Parallel::pull_parallel_vector_data<AbPair>(
303 1978 : _communicator, query_ids, gather_data, act_on_data, 0);
304 1978 : }
305 :
306 : void
307 75943 : NodalPatchRecoveryBase::addToQuery(
308 : const libMesh::Elem * elem,
309 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> & query_ids) const
310 : {
311 75943 : if (hasBlocks(elem->subdomain_id()) && elem->processor_id() != processor_id())
312 4464 : query_ids[elem->processor_id()].push_back(elem->id());
313 75943 : }
|
