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 6435 : removeDuplicateEntries(const std::vector<dof_id_type> & ids)
28 : {
29 6435 : std::vector<dof_id_type> key = ids;
30 6435 : std::sort(key.begin(), key.end());
31 6435 : key.erase(std::unique(key.begin(), key.end()), key.end());
32 6435 : return key;
33 0 : }
34 :
35 : InputParameters
36 6841 : NodalPatchRecoveryBase::validParams()
37 : {
38 6841 : InputParameters params = ElementUserObject::validParams();
39 :
40 27364 : MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
41 20523 : 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 20523 : params.addRelationshipManager("ElementSideNeighborLayers",
48 : Moose::RelationshipManagerType::ALGEBRAIC,
49 7858 : [](const InputParameters &, InputParameters & rm_params)
50 2034 : { rm_params.set<unsigned short>("layers") = 2; });
51 :
52 20523 : params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
53 :
54 13682 : return params;
55 6841 : }
56 :
57 374 : NodalPatchRecoveryBase::NodalPatchRecoveryBase(const InputParameters & parameters)
58 : : ElementUserObject(parameters),
59 374 : _qp(0),
60 374 : _patch_polynomial_order(
61 374 : static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
62 374 : _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
63 374 : _q(_multi_index.size()),
64 374 : _distributed_mesh(_mesh.isDistributedMesh()),
65 1496 : _proc_ids(n_processors())
66 : {
67 374 : std::iota(_proc_ids.begin(), _proc_ids.end(), 0);
68 374 : }
69 :
70 : Real
71 4961 : NodalPatchRecoveryBase::nodalPatchRecovery(const Point & x,
72 : const std::vector<dof_id_type> & elem_ids) const
73 : {
74 4961 : const RealEigenVector coef = getCoefficients(elem_ids); // const version
75 : // Compute the fitted nodal value
76 4961 : RealEigenVector p = evaluateBasisFunctions(x);
77 9922 : return p.dot(coef);
78 4961 : }
79 :
80 : const RealEigenVector
81 5698 : NodalPatchRecoveryBase::getCoefficients(const std::vector<dof_id_type> & elem_ids) const
82 : {
83 5698 : auto elem_ids_reduced = removeDuplicateEntries(elem_ids);
84 :
85 5698 : RealEigenVector coef = RealEigenVector::Zero(_q);
86 : // Before we go, check if we have enough sample points for solving the least square fitting
87 5698 : 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 5698 : RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
93 5698 : RealEigenVector b = RealEigenVector::Zero(_q);
94 28578 : for (auto elem_id : elem_ids_reduced)
95 : {
96 22880 : const auto elem = _mesh.elemPtr(elem_id);
97 22880 : if (elem /*prevent segmentation fault in distributed mesh*/)
98 22559 : 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 22880 : if (_Ae.find(elem_id) == _Ae.end())
105 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _Ae.");
106 22880 : if (_be.find(elem_id) == _be.end())
107 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _be.");
108 :
109 22880 : A += libmesh_map_find(_Ae, elem_id);
110 22880 : b += libmesh_map_find(_be, elem_id);
111 : }
112 :
113 : // Solve the least squares fitting
114 5698 : coef = A.completeOrthogonalDecomposition().solve(b);
115 :
116 11396 : return coef;
117 5698 : }
118 :
119 : const RealEigenVector
120 737 : NodalPatchRecoveryBase::getCachedCoefficients(const std::vector<dof_id_type> & elem_ids)
121 : {
122 : // Check cache
123 737 : auto key = removeDuplicateEntries(elem_ids);
124 :
125 737 : if (key == _cached_elem_ids)
126 0 : return _cached_coef;
127 : else
128 : {
129 737 : const auto coef = getCoefficients(key); // const version
130 :
131 737 : _cached_elem_ids = key; // Update the cached element IDs
132 737 : _cached_coef = coef; // Update the cached coefficients
133 :
134 737 : return coef;
135 737 : }
136 737 : }
137 :
138 : RealEigenVector
139 121799 : NodalPatchRecoveryBase::evaluateBasisFunctions(const Point & q_point) const
140 : {
141 121799 : RealEigenVector p(_q);
142 : Real polynomial;
143 716572 : for (unsigned int r = 0; r < _multi_index.size(); r++)
144 : {
145 594773 : polynomial = 1.0;
146 : mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
147 2111999 : for (unsigned int c = 0; c < _multi_index[r].size(); c++)
148 2186808 : for (unsigned int p = 0; p < _multi_index[r][c]; p++)
149 669582 : polynomial *= q_point(c);
150 594773 : p(r) = polynomial;
151 : }
152 121799 : return p;
153 0 : }
154 :
155 : void
156 1262 : 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 1262 : _cached_elem_ids.clear();
163 1262 : _cached_coef = RealEigenVector::Zero(_q);
164 1262 : _Ae.clear();
165 1262 : _be.clear();
166 1262 : }
167 :
168 : void
169 24276 : NodalPatchRecoveryBase::execute()
170 : {
171 24276 : RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
172 24276 : RealEigenVector be = RealEigenVector::Zero(_q);
173 141114 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
174 : {
175 116838 : RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
176 116838 : Ae += p * p.transpose();
177 116838 : be += computeValue() * p;
178 116838 : }
179 :
180 24276 : dof_id_type elem_id = _current_elem->id();
181 :
182 24276 : _Ae[elem_id] = Ae;
183 24276 : _be[elem_id] = be;
184 24276 : }
185 :
186 : void
187 103 : NodalPatchRecoveryBase::threadJoin(const UserObject & uo)
188 : {
189 103 : const auto & npr = static_cast<const NodalPatchRecoveryBase &>(uo);
190 103 : _Ae.insert(npr._Ae.begin(), npr._Ae.end());
191 103 : _be.insert(npr._be.begin(), npr._be.end());
192 103 : }
193 :
194 : void
195 1159 : 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 1159 : sync();
201 1159 : }
202 :
203 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
204 737 : NodalPatchRecoveryBase::gatherRequestList(const std::vector<dof_id_type> & specific_elems)
205 : {
206 737 : 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 737 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
210 :
211 6725 : for (const auto & entry : specific_elems)
212 : {
213 5988 : const auto * elem = _mesh.elemPtr(entry);
214 5988 : if (_distributed_mesh)
215 : {
216 1048 : if (!elem)
217 321 : continue; // Prevent segmentation fault in distributed mesh
218 :
219 727 : if (hasBlocks(elem->subdomain_id()))
220 2181 : for (processor_id_type pid = 0; pid < n_processors(); ++pid)
221 1454 : if (pid != processor_id())
222 727 : 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 4940 : addToQuery(elem, query_ids);
227 : }
228 :
229 737 : if (_distributed_mesh)
230 : {
231 : auto push_receiver =
232 96 : [&](const processor_id_type, const std::vector<PidElemPair> & received_data)
233 : {
234 823 : for (const auto & [pid, id] : received_data)
235 727 : query_ids[pid].push_back(id);
236 96 : };
237 :
238 130 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
239 : }
240 :
241 1474 : return query_ids;
242 737 : }
243 :
244 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
245 1159 : NodalPatchRecoveryBase::gatherRequestList()
246 : {
247 1159 : 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 1159 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
251 :
252 67133 : for (const auto & elem : _fe_problem.getEvaluableElementRange())
253 65974 : addToQuery(elem, query_ids);
254 :
255 2318 : return query_ids;
256 1159 : }
257 :
258 : void
259 1159 : NodalPatchRecoveryBase::sync()
260 : {
261 1159 : const auto query_ids = gatherRequestList();
262 1159 : syncHelper(query_ids);
263 1159 : }
264 :
265 : void
266 737 : NodalPatchRecoveryBase::sync(const std::vector<dof_id_type> & specific_elems)
267 : {
268 737 : const auto query_ids = gatherRequestList(specific_elems);
269 737 : syncHelper(query_ids);
270 737 : }
271 :
272 : void
273 1896 : 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 849 : 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 6178 : for (const auto & elem_id : elem_ids)
284 5329 : ab_pairs.emplace_back(libmesh_map_find(_Ae, elem_id), libmesh_map_find(_be, elem_id));
285 849 : };
286 :
287 : // Gather answers received from other processors
288 849 : 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 6178 : for (const auto i : index_range(elem_ids))
293 : {
294 5329 : const auto elem_id = elem_ids[i];
295 5329 : const auto & [Ae, be] = ab_pairs[i];
296 5329 : _Ae[elem_id] = Ae;
297 5329 : _be[elem_id] = be;
298 : }
299 849 : };
300 :
301 : // the send and receive are called inside the pull_parallel_vector_data function
302 1896 : libMesh::Parallel::pull_parallel_vector_data<AbPair>(
303 1896 : _communicator, query_ids, gather_data, act_on_data, 0);
304 1896 : }
305 :
306 : void
307 70914 : NodalPatchRecoveryBase::addToQuery(
308 : const libMesh::Elem * elem,
309 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> & query_ids) const
310 : {
311 70914 : if (hasBlocks(elem->subdomain_id()) && elem->processor_id() != processor_id())
312 4602 : query_ids[elem->processor_id()].push_back(elem->id());
313 70914 : }
|
