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 "AuxKernel.h"
11 :
12 : // local includes
13 : #include "FEProblem.h"
14 : #include "SubProblem.h"
15 : #include "AuxiliarySystem.h"
16 : #include "MooseTypes.h"
17 : #include "Assembly.h"
18 : #include "MortarNodalAuxKernel.h"
19 :
20 : #include "libmesh/numeric_vector.h"
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/quadrature.h"
23 : #include "libmesh/boundary_info.h"
24 :
25 : template <typename ComputeValueType>
26 : InputParameters
27 578902 : AuxKernelTempl<ComputeValueType>::validParams()
28 : {
29 578902 : InputParameters params = AuxKernelBase::validParams();
30 :
31 578902 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(VectorAuxKernel).name())
32 43726 : params.registerBase("VectorAuxKernel");
33 578902 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(ArrayAuxKernel).name())
34 49636 : params.registerBase("ArrayAuxKernel");
35 :
36 578902 : params.registerSystemAttributeName("AuxKernel");
37 :
38 578902 : return params;
39 0 : }
40 :
41 : template <typename ComputeValueType>
42 70559 : AuxKernelTempl<ComputeValueType>::AuxKernelTempl(const InputParameters & parameters)
43 : : AuxKernelBase(parameters),
44 : MooseVariableInterface<ComputeValueType>(
45 : this,
46 : parameters.getCheckedPointerParam<SystemBase *>("_sys")
47 282224 : ->getVariable(parameters.get<THREAD_ID>("_tid"),
48 141112 : parameters.get<AuxVariableName>("variable"))
49 70556 : .isNodal(),
50 : "variable",
51 : Moose::VarKindType::VAR_AUXILIARY,
52 : std::is_same<Real, ComputeValueType>::value
53 : ? Moose::VarFieldType::VAR_FIELD_STANDARD
54 : : (std::is_same<RealVectorValue, ComputeValueType>::value
55 : ? Moose::VarFieldType::VAR_FIELD_VECTOR
56 : : Moose::VarFieldType::VAR_FIELD_ARRAY)),
57 :
58 211668 : _var(_aux_sys.getActualFieldVariable<ComputeValueType>(
59 70556 : _tid, parameters.get<AuxVariableName>("variable"))),
60 70556 : _nodal(_var.isNodal()),
61 70556 : _u(_nodal ? _var.nodalValueArray() : _var.sln()),
62 :
63 70556 : _test(_bnd ? _var.phiFace() : _var.phi()),
64 70556 : _q_point(_bnd ? _assembly.qPointsFace() : _assembly.qPoints()),
65 70556 : _qrule(_bnd ? _assembly.qRuleFace() : _assembly.qRule()),
66 70556 : _JxW(_bnd ? _assembly.JxWFace() : _assembly.JxW()),
67 70556 : _coord(_assembly.coordTransformation()),
68 :
69 70556 : _current_elem(_assembly.elem()),
70 70556 : _current_side(_assembly.side()),
71 70556 : _current_elem_volume(_assembly.elemVolume()),
72 70556 : _current_side_volume(_assembly.sideElemVolume()),
73 :
74 70556 : _current_node(_assembly.node()),
75 70556 : _current_boundary_id(_assembly.currentBoundaryID()),
76 70556 : _solution(_aux_sys.solution()),
77 :
78 423339 : _current_lower_d_elem(_assembly.lowerDElem())
79 : {
80 :
81 70556 : if (!_bnd || _nodal)
82 : // If we're not boundary restricted then we cannot be a coincident lower-d calculation
83 69175 : _coincident_lower_d_calc = false;
84 70556 : }
85 :
86 : template <typename ComputeValueType>
87 : const VariableValue &
88 68 : AuxKernelTempl<ComputeValueType>::coupledDot(const std::string & var_name, unsigned int comp) const
89 : {
90 68 : auto var = getVar(var_name, comp);
91 68 : if (var->kind() == Moose::VAR_AUXILIARY)
92 3 : mooseError(
93 3 : name(),
94 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
95 :
96 65 : return Coupleable::coupledDot(var_name, comp);
97 : }
98 :
99 : template <typename ComputeValueType>
100 : const VariableValue &
101 0 : AuxKernelTempl<ComputeValueType>::coupledDotDu(const std::string & var_name,
102 : unsigned int comp) const
103 : {
104 0 : auto var = getVar(var_name, comp);
105 0 : if (var->kind() == Moose::VAR_AUXILIARY)
106 0 : mooseError(
107 0 : name(),
108 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
109 :
110 0 : return Coupleable::coupledDotDu(var_name, comp);
111 : }
112 :
113 : template <>
114 : void
115 173208 : AuxKernelTempl<Real>::setDofValueHelper(const Real & value)
116 : {
117 : mooseAssert(_n_shapes == 1,
118 : "Should only be calling setDofValue if there is one dof for the aux var");
119 173208 : _var.setDofValue(value, 0);
120 173208 : }
121 :
122 : template <>
123 : void
124 0 : AuxKernelTempl<RealVectorValue>::setDofValueHelper(const RealVectorValue &)
125 : {
126 0 : mooseError("Not implemented");
127 : }
128 :
129 : template <typename ComputeValueType>
130 : void
131 122516 : AuxKernelTempl<ComputeValueType>::insert()
132 : {
133 : mooseAssert(_coincident_lower_d_calc.has_value(),
134 : "We should have set _coincident_lower_d_calc by now");
135 122516 : if (_coincident_lower_d_calc.value())
136 1644 : _var.insertLower(_aux_sys.solution());
137 : else
138 120872 : _var.insert(_aux_sys.solution());
139 122516 : }
140 :
141 : template <typename ComputeValueType>
142 : void
143 46438090 : AuxKernelTempl<ComputeValueType>::compute()
144 : {
145 46438090 : precalculateValue();
146 :
147 46438090 : if (isNodal()) /* nodal variables */
148 : {
149 30971948 : if (_var.isNodalDefined())
150 : {
151 28540649 : _qp = 0;
152 28540649 : ComputeValueType value = computeValue();
153 : // update variable data, which is referenced by other kernels, so the value is up-to-date
154 28540647 : _var.setNodalValue(value);
155 : }
156 : }
157 : else /* elemental variables */
158 : {
159 1644 : _n_shapes =
160 15466142 : _coincident_lower_d_calc.value() ? _var.dofIndicesLower().size() : _var.numberOfDofs();
161 :
162 15466142 : if (_n_shapes == 1) /* p0 */
163 : {
164 14889091 : ComputeValueType value = 0;
165 88421862 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
166 73532782 : value += _JxW[_qp] * _coord[_qp] * computeValue();
167 14889080 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
168 14889080 : if (_var.isFV())
169 173208 : setDofValueHelper(value);
170 : else
171 : {
172 : // update the variable data referenced by other kernels.
173 : // Note that this will update the values at the quadrature points too
174 : // (because this is an Elemental variable)
175 14715872 : if (_coincident_lower_d_calc.value())
176 : {
177 822 : _local_sol.resize(1);
178 : if constexpr (std::is_same<Real, ComputeValueType>::value)
179 822 : _local_sol(0) = value;
180 : else
181 : mooseAssert(false, "We should not enter the single dof branch with a vector variable");
182 822 : _var.setLowerDofValues(_local_sol);
183 : }
184 : else
185 14715050 : _var.setNodalValue(value);
186 : }
187 : }
188 : else /* high-order */
189 : {
190 577051 : _local_re.resize(_n_shapes);
191 577051 : _local_re.zero();
192 577051 : _local_ke.resize(_n_shapes, _n_shapes);
193 577051 : _local_ke.zero();
194 :
195 577051 : const auto & test = _coincident_lower_d_calc.value() ? _var.phiLower() : _test;
196 :
197 : // assemble the local mass matrix and the load
198 3333027 : for (unsigned int i = 0; i < test.size(); i++)
199 25143102 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
200 : {
201 22387126 : ComputeValueType t = _JxW[_qp] * _coord[_qp] * test[i][_qp];
202 22387126 : _local_re(i) += t * computeValue();
203 174525132 : for (unsigned int j = 0; j < test.size(); j++)
204 152138010 : _local_ke(i, j) += t * test[j][_qp];
205 : }
206 : // mass matrix is always SPD but in case of boundary restricted, it will be rank deficient
207 577047 : _local_sol.resize(_n_shapes);
208 577047 : if (_bnd)
209 5094 : _local_ke.svd_solve(_local_re, _local_sol);
210 : else
211 571953 : _local_ke.cholesky_solve(_local_re, _local_sol);
212 :
213 577047 : _coincident_lower_d_calc.value() ? _var.setLowerDofValues(_local_sol)
214 576225 : : _var.setDofValues(_local_sol);
215 : }
216 : }
217 46438073 : }
218 :
219 : template <>
220 : void
221 589931 : AuxKernelTempl<RealEigenVector>::compute()
222 : {
223 589931 : precalculateValue();
224 :
225 589931 : if (isNodal()) /* nodal variables */
226 : {
227 585450 : if (_var.isNodalDefined())
228 : {
229 572650 : _qp = 0;
230 572650 : RealEigenVector value = computeValue();
231 : // update variable data, which is referenced by other kernels, so the value is up-to-date
232 572650 : _var.setNodalValue(value);
233 572650 : }
234 : }
235 : else /* elemental variables */
236 : {
237 : mooseAssert(
238 : _var.numberOfDofs() % _var.count() == 0,
239 : "The number of degrees of freedom should be cleanly divisible by the variable count");
240 4481 : _n_shapes = _var.numberOfDofs() / _var.count();
241 4481 : if (_n_shapes == 1) /* p0 */
242 : {
243 0 : RealEigenVector value = RealEigenVector::Zero(_var.count());
244 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
245 0 : value += _JxW[_qp] * _coord[_qp] * computeValue();
246 0 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
247 : // update the variable data referenced by other kernels.
248 : // Note that this will update the values at the quadrature points too
249 : // (because this is an Elemental variable)
250 0 : _var.setNodalValue(value);
251 0 : }
252 : else /* high-order */
253 : {
254 4481 : _local_re.resize(_n_shapes);
255 117955 : for (unsigned int i = 0; i < _local_re.size(); ++i)
256 113474 : _local_re(i) = RealEigenVector::Zero(_var.count());
257 4481 : _local_ke.resize(_n_shapes, _n_shapes);
258 4481 : _local_ke.zero();
259 :
260 : // assemble the local mass matrix and the load
261 117955 : for (unsigned int i = 0; i < _test.size(); i++)
262 5606470 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
263 : {
264 5492996 : Real t = _JxW[_qp] * _coord[_qp] * _test[i][_qp];
265 5492996 : _local_re(i) += t * computeValue();
266 159175180 : for (unsigned int j = 0; j < _test.size(); j++)
267 153682184 : _local_ke(i, j) += t * _test[j][_qp];
268 : }
269 :
270 : // mass matrix is always SPD
271 4481 : _local_sol.resize(_n_shapes);
272 117955 : for (unsigned int i = 0; i < _local_re.size(); ++i)
273 113474 : _local_sol(i) = RealEigenVector::Zero(_var.count());
274 4481 : DenseVector<Number> re(_n_shapes);
275 4481 : DenseVector<Number> sol(_n_shapes);
276 17443 : for (unsigned int i = 0; i < _var.count(); ++i)
277 : {
278 351910 : for (unsigned int j = 0; j < _n_shapes; ++j)
279 338948 : re(j) = _local_re(j)(i);
280 :
281 12962 : if (_bnd)
282 0 : _local_ke.svd_solve(re, sol);
283 : else
284 12962 : _local_ke.cholesky_solve(re, sol);
285 :
286 351910 : for (unsigned int j = 0; j < _n_shapes; ++j)
287 338948 : _local_sol(j)(i) = sol(j);
288 : }
289 :
290 4481 : _var.setDofValues(_local_sol);
291 4481 : }
292 : }
293 589931 : }
294 :
295 : template <typename ComputeValueType>
296 : const typename OutputTools<ComputeValueType>::VariableValue &
297 117 : AuxKernelTempl<ComputeValueType>::uOld() const
298 : {
299 117 : if (_sys.solutionStatesInitialized())
300 0 : mooseError("The solution states have already been initialized when calling ",
301 0 : type(),
302 : "::uOld().\n\n",
303 : "Make sure to call uOld() within the object constructor.");
304 :
305 117 : return _nodal ? _var.nodalValueOldArray() : _var.slnOld();
306 : }
307 :
308 : template <typename ComputeValueType>
309 : const typename OutputTools<ComputeValueType>::VariableValue &
310 13 : AuxKernelTempl<ComputeValueType>::uOlder() const
311 : {
312 13 : if (_sys.solutionStatesInitialized())
313 0 : mooseError("The solution states have already been initialized when calling ",
314 0 : type(),
315 : "::uOlder().\n\n",
316 : "Make sure to call uOlder() within the object constructor.");
317 :
318 13 : return _nodal ? _var.nodalValueOlderArray() : _var.slnOlder();
319 : }
320 :
321 : template <typename ComputeValueType>
322 : bool
323 22682 : AuxKernelTempl<ComputeValueType>::isMortar()
324 : {
325 22682 : return dynamic_cast<MortarNodalAuxKernelTempl<ComputeValueType> *>(this) != nullptr;
326 : }
327 :
328 : template <typename ComputeValueType>
329 : void
330 122524 : AuxKernelTempl<ComputeValueType>::determineWhetherCoincidentLowerDCalc()
331 : {
332 : mooseAssert(_bnd && !isNodal(),
333 : "We can never be a lower-dimensional calculation if we are not boundary restricted "
334 : "or if we are a nodal auxiliary kernel");
335 :
336 : // MOOSE maintains three copies of finite element data. One is for the current canonical element,
337 : // one is for the element neighbor, and one is for a lower-d element. These three copies have
338 : // supported all of MOOSE's finite element use cases to date, including mortar. For AuxKernel, we
339 : // do not use the neighbor data copy, but we do use the other two copies. The numberOfDofs() API
340 : // returns the number of degrees of freedom for the canonical element. If there are none, then
341 : // there must be degrees of freedom associated with the lower-d element
342 122524 : _coincident_lower_d_calc = !_var.numberOfDofs();
343 :
344 122524 : if (_coincident_lower_d_calc.value())
345 : {
346 1708 : static const std::string lower_error = "Make sure that the lower-d variable lives on a "
347 : "lower-d block that is a superset of the boundary";
348 1652 : if (!_current_lower_d_elem)
349 4 : mooseError("No lower-dimensional element. ", lower_error);
350 1648 : if (!_var.dofIndicesLower().size())
351 4 : mooseError("No degrees of freedom. ", lower_error);
352 : }
353 122516 : }
354 :
355 : // Explicitly instantiates the three versions of the AuxKernelTempl class
356 : template class AuxKernelTempl<Real>;
357 : template class AuxKernelTempl<RealVectorValue>;
358 : template class AuxKernelTempl<RealEigenVector>;
|
