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 36324 : removeDuplicateEntries(const std::vector<dof_id_type> & ids)
28 : {
29 36324 : std::vector<dof_id_type> key = ids;
30 36324 : std::sort(key.begin(), key.end());
31 36324 : key.erase(std::unique(key.begin(), key.end()), key.end());
32 36324 : return key;
33 0 : }
34 :
35 : InputParameters
36 7049 : NodalPatchRecoveryBase::validParams()
37 : {
38 7049 : InputParameters params = ElementUserObject::validParams();
39 :
40 28196 : MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
41 21147 : 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 21147 : params.addRelationshipManager("ElementSideNeighborLayers",
48 : Moose::RelationshipManagerType::ALGEBRAIC,
49 8282 : [](const InputParameters &, InputParameters & rm_params)
50 : {
51 2466 : rm_params.set<bool>("use_point_neighbors") = true;
52 1233 : rm_params.set<unsigned short>("layers") = 1;
53 1233 : });
54 :
55 21147 : params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
56 :
57 14098 : return params;
58 7049 : }
59 :
60 454 : NodalPatchRecoveryBase::NodalPatchRecoveryBase(const InputParameters & parameters)
61 : : ElementUserObject(parameters),
62 454 : _qp(0),
63 454 : _patch_polynomial_order(
64 454 : static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
65 454 : _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
66 454 : _q(_multi_index.size()),
67 454 : _distributed_mesh(_mesh.isDistributedMesh()),
68 1816 : _proc_ids(n_processors())
69 : {
70 454 : std::iota(_proc_ids.begin(), _proc_ids.end(), 0);
71 454 : }
72 :
73 : Real
74 34850 : NodalPatchRecoveryBase::nodalPatchRecovery(const Point & x,
75 : const std::vector<dof_id_type> & elem_ids) const
76 : {
77 34850 : const RealEigenVector coef = getCoefficients(elem_ids); // const version
78 : // Compute the fitted nodal value
79 34850 : RealEigenVector p = evaluateBasisFunctions(x);
80 69700 : return p.dot(coef);
81 34850 : }
82 :
83 : const RealEigenVector
84 35587 : NodalPatchRecoveryBase::getCoefficients(const std::vector<dof_id_type> & elem_ids) const
85 : {
86 35587 : auto elem_ids_reduced = removeDuplicateEntries(elem_ids);
87 :
88 35587 : RealEigenVector coef = RealEigenVector::Zero(_q);
89 : // Before we go, check if we have enough sample points for solving the least square fitting
90 35587 : if (_q_point.size() * elem_ids_reduced.size() < _q)
91 0 : mooseError("There are not enough sample points to recover the nodal value, try reducing the "
92 : "polynomial order or using a higher-order quadrature scheme.");
93 :
94 : // Assemble the least squares problem over the patch
95 35587 : RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
96 35587 : RealEigenVector b = RealEigenVector::Zero(_q);
97 228699 : for (auto elem_id : elem_ids_reduced)
98 : {
99 193112 : const auto elem = _mesh.elemPtr(elem_id);
100 193112 : if (elem /*prevent segmentation fault in distributed mesh*/)
101 192791 : if (!hasBlocks(elem->subdomain_id()))
102 0 : mooseError("Element with id = ",
103 : elem_id,
104 : " is not in the block. "
105 : "Please use nodalPatchRecovery with elements in the block only.");
106 :
107 193112 : if (_Ae.find(elem_id) == _Ae.end())
108 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _Ae.");
109 193112 : if (_be.find(elem_id) == _be.end())
110 0 : mooseError("Missing entry for elem_id = ", elem_id, " in _be.");
111 :
112 193112 : A += libmesh_map_find(_Ae, elem_id);
113 193112 : b += libmesh_map_find(_be, elem_id);
114 : }
115 :
116 : // Solve the least squares fitting
117 35587 : coef = A.completeOrthogonalDecomposition().solve(b);
118 :
119 71174 : return coef;
120 35587 : }
121 :
122 : const RealEigenVector
123 737 : NodalPatchRecoveryBase::getCachedCoefficients(const std::vector<dof_id_type> & elem_ids)
124 : {
125 : // Check cache
126 737 : auto key = removeDuplicateEntries(elem_ids);
127 :
128 737 : if (key == _cached_elem_ids)
129 0 : return _cached_coef;
130 : else
131 : {
132 737 : const auto coef = getCoefficients(key); // const version
133 :
134 737 : _cached_elem_ids = key; // Update the cached element IDs
135 737 : _cached_coef = coef; // Update the cached coefficients
136 :
137 737 : return coef;
138 737 : }
139 737 : }
140 :
141 : RealEigenVector
142 319624 : NodalPatchRecoveryBase::evaluateBasisFunctions(const Point & q_point) const
143 : {
144 319624 : RealEigenVector p(_q);
145 : Real polynomial;
146 1705697 : for (unsigned int r = 0; r < _multi_index.size(); r++)
147 : {
148 1386073 : polynomial = 1.0;
149 : mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
150 5277199 : for (unsigned int c = 0; c < _multi_index[r].size(); c++)
151 5154183 : for (unsigned int p = 0; p < _multi_index[r][c]; p++)
152 1263057 : polynomial *= q_point(c);
153 1386073 : p(r) = polynomial;
154 : }
155 319624 : return p;
156 0 : }
157 :
158 : void
159 1630 : NodalPatchRecoveryBase::initialize()
160 : {
161 : // Clear cached data to ensure coefficients are correctly recomputed in the next patch recovery
162 : // iteration. _Ae and _be must also be reset, as the associated patch elements may differ in the
163 : // upcoming iteration.
164 :
165 1630 : _cached_elem_ids.clear();
166 1630 : _cached_coef = RealEigenVector::Zero(_q);
167 1630 : _Ae.clear();
168 1630 : _be.clear();
169 1630 : }
170 :
171 : void
172 45268 : NodalPatchRecoveryBase::execute()
173 : {
174 45268 : RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
175 45268 : RealEigenVector be = RealEigenVector::Zero(_q);
176 330042 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
177 : {
178 284774 : RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
179 284774 : Ae += p * p.transpose();
180 284774 : be += computeValue() * p;
181 284774 : }
182 :
183 45268 : dof_id_type elem_id = _current_elem->id();
184 :
185 45268 : _Ae[elem_id] = Ae;
186 45268 : _be[elem_id] = be;
187 45268 : }
188 :
189 : void
190 143 : NodalPatchRecoveryBase::threadJoin(const UserObject & uo)
191 : {
192 143 : const auto & npr = static_cast<const NodalPatchRecoveryBase &>(uo);
193 143 : _Ae.insert(npr._Ae.begin(), npr._Ae.end());
194 143 : _be.insert(npr._be.begin(), npr._be.end());
195 143 : }
196 :
197 : void
198 1487 : NodalPatchRecoveryBase::finalize()
199 : {
200 : // When calling nodalPatchRecovery, we may need to know _Ae and _be on algebraically ghosted
201 : // elements. However, this userobject is only run on local elements, so we need to query those
202 : // information from other processors in this finalize() method.
203 1487 : sync();
204 1487 : }
205 :
206 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
207 737 : NodalPatchRecoveryBase::gatherRequestList(const std::vector<dof_id_type> & specific_elems)
208 : {
209 737 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
210 :
211 : typedef std::pair<processor_id_type, dof_id_type> PidElemPair;
212 737 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
213 :
214 6725 : for (const auto & entry : specific_elems)
215 : {
216 5988 : const auto * elem = _mesh.elemPtr(entry);
217 5988 : if (_distributed_mesh)
218 : {
219 1048 : if (!elem)
220 321 : continue; // Prevent segmentation fault in distributed mesh
221 :
222 727 : if (hasBlocks(elem->subdomain_id()))
223 2181 : for (processor_id_type pid = 0; pid < n_processors(); ++pid)
224 1454 : if (pid != processor_id())
225 727 : push_data[pid].push_back(std::make_pair(elem->processor_id(), elem->id()));
226 : }
227 : else
228 : // Add to query_ids if this element's data is on a different processor
229 4940 : addToQuery(elem, query_ids);
230 : }
231 :
232 737 : if (_distributed_mesh)
233 : {
234 : auto push_receiver =
235 96 : [&](const processor_id_type, const std::vector<PidElemPair> & received_data)
236 : {
237 823 : for (const auto & [pid, id] : received_data)
238 727 : query_ids[pid].push_back(id);
239 96 : };
240 :
241 130 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
242 : }
243 :
244 1474 : return query_ids;
245 737 : }
246 :
247 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
248 1487 : NodalPatchRecoveryBase::gatherRequestList()
249 : {
250 1487 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
251 :
252 : typedef std::pair<processor_id_type, dof_id_type> PidElemPair;
253 1487 : std::unordered_map<processor_id_type, std::vector<PidElemPair>> push_data;
254 :
255 109087 : for (const auto & elem : _fe_problem.getEvaluableElementRange())
256 107600 : addToQuery(elem, query_ids);
257 :
258 2974 : return query_ids;
259 1487 : }
260 :
261 : void
262 1487 : NodalPatchRecoveryBase::sync()
263 : {
264 1487 : const auto query_ids = gatherRequestList();
265 1487 : syncHelper(query_ids);
266 1487 : }
267 :
268 : void
269 737 : NodalPatchRecoveryBase::sync(const std::vector<dof_id_type> & specific_elems)
270 : {
271 737 : const auto query_ids = gatherRequestList(specific_elems);
272 737 : syncHelper(query_ids);
273 737 : }
274 :
275 : void
276 2224 : NodalPatchRecoveryBase::syncHelper(
277 : const std::unordered_map<processor_id_type, std::vector<dof_id_type>> & query_ids)
278 : {
279 : typedef std::pair<RealEigenMatrix, RealEigenVector> AbPair;
280 :
281 : // Answer queries received from other processors
282 2617 : auto gather_data = [this](const processor_id_type /*pid*/,
283 : const std::vector<dof_id_type> & elem_ids,
284 : std::vector<AbPair> & ab_pairs)
285 : {
286 30220 : for (const auto & elem_id : elem_ids)
287 27603 : ab_pairs.emplace_back(libmesh_map_find(_Ae, elem_id), libmesh_map_find(_be, elem_id));
288 2617 : };
289 :
290 : // Gather answers received from other processors
291 2617 : auto act_on_data = [this](const processor_id_type /*pid*/,
292 : const std::vector<dof_id_type> & elem_ids,
293 : const std::vector<AbPair> & ab_pairs)
294 : {
295 30220 : for (const auto i : index_range(elem_ids))
296 : {
297 27603 : const auto elem_id = elem_ids[i];
298 27603 : const auto & [Ae, be] = ab_pairs[i];
299 27603 : _Ae[elem_id] = Ae;
300 27603 : _be[elem_id] = be;
301 : }
302 2617 : };
303 :
304 : // the send and receive are called inside the pull_parallel_vector_data function
305 2224 : libMesh::Parallel::pull_parallel_vector_data<AbPair>(
306 2224 : _communicator, query_ids, gather_data, act_on_data, 0);
307 2224 : }
308 :
309 : void
310 112540 : NodalPatchRecoveryBase::addToQuery(
311 : const libMesh::Elem * elem,
312 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> & query_ids) const
313 : {
314 112540 : if (hasBlocks(elem->subdomain_id()) && elem->processor_id() != processor_id())
315 26876 : query_ids[elem->processor_id()].push_back(elem->id());
316 112540 : }
|
