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 "ProjectedStatefulMaterialNodalPatchRecovery.h"
11 : #include "ElementUserObject.h"
12 : #include "MaterialBase.h"
13 : #include "MathUtils.h"
14 : #include "Assembly.h"
15 :
16 : // TIMPI includes
17 : #include "timpi/communicator.h"
18 : #include "timpi/parallel_sync.h"
19 : #include "libmesh/parallel_eigen.h"
20 :
21 : registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryReal);
22 : registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryReal);
23 : registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue);
24 : registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue);
25 : registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor);
26 : registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor);
27 : registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor);
28 : registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor);
29 :
30 : InputParameters
31 114328 : ProjectedStatefulMaterialNodalPatchRecoveryBase::validParams()
32 : {
33 114328 : InputParameters params = ElementUserObject::validParams();
34 114328 : params.addClassDescription(
35 : "Prepare patches for use in nodal patch recovery based on a material property for material "
36 : "property states projected onto nodal variables.");
37 114328 : params.addRelationshipManager("ElementSideNeighborLayers",
38 : Moose::RelationshipManagerType::ALGEBRAIC,
39 0 : [](const InputParameters &, InputParameters & rm_params)
40 108 : { rm_params.set<unsigned short>("layers") = 2; });
41 :
42 114328 : return params;
43 0 : }
44 :
45 : template <typename T, bool is_ad>
46 : InputParameters
47 114220 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::validParams()
48 : {
49 114220 : InputParameters params = ProjectedStatefulMaterialNodalPatchRecoveryBase::validParams();
50 114220 : MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
51 114220 : params.addRequiredParam<MaterialPropertyName>(
52 : "property", "Name of the material property to perform nodal patch recovery on");
53 114220 : params.addRequiredParam<MooseEnum>(
54 : "patch_polynomial_order",
55 : orders,
56 : "Polynomial order used in least squares fitting of material property "
57 : "over the local patch of elements connected to a given node");
58 :
59 114220 : params.addRelationshipManager("ElementSideNeighborLayers",
60 : Moose::RelationshipManagerType::ALGEBRAIC,
61 0 : [](const InputParameters &, InputParameters & rm_params)
62 0 : { rm_params.set<unsigned short>("layers") = 2; });
63 :
64 114220 : params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
65 :
66 228440 : return params;
67 114220 : }
68 :
69 : template <typename T, bool is_ad>
70 52 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::
71 : ProjectedStatefulMaterialNodalPatchRecoveryTempl(const InputParameters & parameters)
72 : : ProjectedStatefulMaterialNodalPatchRecoveryBase(parameters),
73 52 : _qp(0),
74 104 : _n_components(Moose::SerialAccess<T>::size()),
75 52 : _patch_polynomial_order(
76 52 : static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
77 52 : _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
78 52 : _q(_multi_index.size()),
79 52 : _prop(getGenericMaterialProperty<T, is_ad>("property")),
80 104 : _current_subdomain_id(_assembly.currentSubdomainID())
81 : {
82 52 : }
83 :
84 : template <typename T, bool is_ad>
85 : void
86 52 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::initialSetup()
87 : {
88 : // get all material classes that provide properties for this object
89 52 : _required_materials = buildRequiredMaterials();
90 52 : }
91 :
92 : template <typename T, bool is_ad>
93 : Real
94 94000 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::nodalPatchRecovery(
95 : const Point & x, const std::vector<dof_id_type> & elem_ids, std::size_t component) const
96 : {
97 : // Before we go, check if we have enough sample points for solving the least square fitting
98 94000 : if (_q_point.size() * elem_ids.size() < _q)
99 0 : mooseError("There are not enough sample points to recover the nodal value, try reducing the "
100 : "polynomial order or using a higher-order quadrature scheme.");
101 :
102 : // Assemble the least squares problem over the patch
103 94000 : RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
104 94000 : RealEigenVector b = RealEigenVector::Zero(_q);
105 379760 : for (const auto elem_id : elem_ids)
106 : {
107 285760 : const auto abs = libmesh_map_find(_abs, elem_id);
108 285760 : A += abs.first;
109 285760 : b += abs.second[component];
110 : }
111 :
112 : // Solve the least squares fitting
113 94000 : const RealEigenVector coef = A.completeOrthogonalDecomposition().solve(b);
114 :
115 : // Compute the fitted nodal value
116 94000 : const RealEigenVector p = evaluateBasisFunctions(x);
117 188000 : return p.dot(coef);
118 94000 : }
119 :
120 : template <typename T, bool is_ad>
121 : RealEigenVector
122 104240 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::evaluateBasisFunctions(
123 : const Point & q_point) const
124 : {
125 104240 : RealEigenVector p(_q);
126 : Real polynomial;
127 416960 : for (const auto r : index_range(_multi_index))
128 : {
129 312720 : polynomial = 1.0;
130 : mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
131 938160 : for (const auto c : index_range(_multi_index[r]))
132 625440 : polynomial *= MathUtils::pow(q_point(c), _multi_index[r][c]);
133 312720 : p(r) = polynomial;
134 : }
135 104240 : return p;
136 0 : }
137 :
138 : template <typename T, bool is_ad>
139 : void
140 240 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::initialize()
141 : {
142 240 : _abs.clear();
143 240 : }
144 :
145 : template <typename T, bool is_ad>
146 : void
147 2560 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::execute()
148 : {
149 2560 : if (_t_step == 0)
150 1024 : for (const auto & mat : _required_materials[_current_subdomain_id])
151 512 : mat->initStatefulProperties(_qrule->size());
152 :
153 2560 : std::vector<RealEigenVector> bs(_n_components, RealEigenVector::Zero(_q));
154 2560 : RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
155 :
156 12800 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
157 : {
158 10240 : RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
159 10240 : Ae += p * p.transpose();
160 :
161 10240 : std::size_t index = 0;
162 250880 : for (const auto & v : Moose::serialAccess(_prop[_qp]))
163 240640 : bs[index++] += MetaPhysicL::raw_value(v) * p;
164 : }
165 :
166 2560 : const dof_id_type elem_id = _current_elem->id();
167 2560 : _abs[elem_id] = {Ae, bs};
168 2560 : }
169 :
170 : template <typename T, bool is_ad>
171 : void
172 20 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::threadJoin(const UserObject & uo)
173 : {
174 20 : const auto & npr =
175 : static_cast<const ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad> &>(uo);
176 20 : _abs.insert(npr._abs.begin(), npr._abs.end());
177 20 : }
178 :
179 : template <typename T, bool is_ad>
180 : void
181 220 : ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::finalize()
182 : {
183 : // When calling nodalPatchRecovery, we may need to know _Ae and _be on algebraically ghosted
184 : // elements. However, this userobject is only run on local elements, so we need to query those
185 : // information from other processors in this finalize() method.
186 :
187 : // Populate algebraically ghosted elements to query
188 220 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
189 220 : const ConstElemRange evaluable_elem_range = _fe_problem.getEvaluableElementRange();
190 3580 : for (const auto * const elem : evaluable_elem_range)
191 3360 : if (elem->processor_id() != processor_id())
192 800 : query_ids[elem->processor_id()].push_back(elem->id());
193 :
194 : // Answer queries received from other processors
195 340 : auto gather_data = [this](const processor_id_type /*pid*/,
196 : const std::vector<dof_id_type> & elem_ids,
197 : std::vector<ElementData> & abs_data)
198 : {
199 920 : for (const auto i : index_range(elem_ids))
200 800 : abs_data.emplace_back(libmesh_map_find(_abs, elem_ids[i]));
201 : };
202 :
203 : // Gather answers received from other processors
204 340 : auto act_on_data = [this](const processor_id_type /*pid*/,
205 : const std::vector<dof_id_type> & elem_ids,
206 : const std::vector<ElementData> & abs_data)
207 : {
208 920 : for (const auto i : index_range(elem_ids))
209 800 : _abs[elem_ids[i]] = abs_data[i];
210 : };
211 :
212 220 : libMesh::Parallel::pull_parallel_vector_data<ElementData>(
213 220 : _communicator, query_ids, gather_data, act_on_data, 0);
214 220 : }
215 :
216 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, false>;
217 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, true>;
218 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, false>;
219 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, true>;
220 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, false>;
221 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, true>;
222 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, false>;
223 : template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, true>;
|
