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 2186384 : AuxKernelTempl<ComputeValueType>::validParams()
28 : {
29 2186384 : InputParameters params = AuxKernelBase::validParams();
30 :
31 2186384 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(VectorAuxKernel).name())
32 200582 : params.registerBase("VectorAuxKernel");
33 2186384 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(ArrayAuxKernel).name())
34 228952 : params.registerBase("ArrayAuxKernel");
35 :
36 2186384 : return params;
37 0 : }
38 :
39 : template <typename ComputeValueType>
40 76108 : AuxKernelTempl<ComputeValueType>::AuxKernelTempl(const InputParameters & parameters)
41 : : AuxKernelBase(parameters),
42 : MooseVariableInterface<ComputeValueType>(
43 : this,
44 : parameters.getCheckedPointerParam<SystemBase *>("_sys")
45 304400 : ->getVariable(parameters.get<THREAD_ID>("_tid"),
46 152200 : parameters.get<AuxVariableName>("variable"))
47 76100 : .isNodal(),
48 : "variable",
49 : Moose::VarKindType::VAR_AUXILIARY,
50 : std::is_same<Real, ComputeValueType>::value
51 : ? Moose::VarFieldType::VAR_FIELD_STANDARD
52 : : (std::is_same<RealVectorValue, ComputeValueType>::value
53 : ? Moose::VarFieldType::VAR_FIELD_VECTOR
54 : : Moose::VarFieldType::VAR_FIELD_ARRAY)),
55 :
56 228300 : _var(_aux_sys.getActualFieldVariable<ComputeValueType>(
57 76100 : _tid, parameters.get<AuxVariableName>("variable"))),
58 76100 : _nodal(_var.isNodal()),
59 76100 : _u(_nodal ? _var.nodalValueArray() : _var.sln()),
60 :
61 76100 : _test(_bnd ? _var.phiFace() : _var.phi()),
62 76100 : _q_point(_bnd ? _assembly.qPointsFace() : _assembly.qPoints()),
63 76100 : _qrule(_bnd ? _assembly.qRuleFace() : _assembly.qRule()),
64 76100 : _JxW(_bnd ? _assembly.JxWFace() : _assembly.JxW()),
65 76100 : _coord(_assembly.coordTransformation()),
66 :
67 76100 : _current_elem(_assembly.elem()),
68 76100 : _current_side(_assembly.side()),
69 76100 : _current_elem_volume(_assembly.elemVolume()),
70 76100 : _current_side_volume(_assembly.sideElemVolume()),
71 :
72 76100 : _current_node(_assembly.node()),
73 76100 : _current_boundary_id(_assembly.currentBoundaryID()),
74 76100 : _solution(_aux_sys.solution()),
75 :
76 456608 : _current_lower_d_elem(_assembly.lowerDElem())
77 : {
78 76100 : }
79 :
80 : template <typename ComputeValueType>
81 : const VariableValue &
82 74 : AuxKernelTempl<ComputeValueType>::coupledDot(const std::string & var_name, unsigned int comp) const
83 : {
84 74 : auto var = getVar(var_name, comp);
85 74 : if (var->kind() == Moose::VAR_AUXILIARY)
86 4 : mooseError(
87 4 : name(),
88 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
89 :
90 70 : return Coupleable::coupledDot(var_name, comp);
91 : }
92 :
93 : template <typename ComputeValueType>
94 : const VariableValue &
95 0 : AuxKernelTempl<ComputeValueType>::coupledDotDu(const std::string & var_name,
96 : unsigned int comp) const
97 : {
98 0 : auto var = getVar(var_name, comp);
99 0 : if (var->kind() == Moose::VAR_AUXILIARY)
100 0 : mooseError(
101 0 : name(),
102 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
103 :
104 0 : return Coupleable::coupledDotDu(var_name, comp);
105 : }
106 :
107 : template <>
108 : void
109 194724 : AuxKernelTempl<Real>::setDofValueHelper(const Real & value)
110 : {
111 : mooseAssert(_n_local_dofs == 1,
112 : "Should only be calling setDofValue if there is one dof for the aux var");
113 194724 : _var.setDofValue(value, 0);
114 194724 : }
115 :
116 : template <>
117 : void
118 0 : AuxKernelTempl<RealVectorValue>::setDofValueHelper(const RealVectorValue &)
119 : {
120 0 : mooseError("Not implemented");
121 : }
122 :
123 : template <typename ComputeValueType>
124 : void
125 137941 : AuxKernelTempl<ComputeValueType>::insert()
126 : {
127 137941 : if (_coincident_lower_d_calc)
128 1852 : _var.insertLower(_aux_sys.solution());
129 : else
130 136089 : _var.insert(_aux_sys.solution());
131 137941 : }
132 :
133 : template <typename ComputeValueType>
134 : void
135 52299322 : AuxKernelTempl<ComputeValueType>::compute()
136 : {
137 52299322 : precalculateValue();
138 :
139 52299322 : if (isNodal()) /* nodal variables */
140 : {
141 : mooseAssert(!_coincident_lower_d_calc,
142 : "Nodal evaluations are point evaluations. We don't have to concern ourselves with "
143 : "coincidence of lower-d blocks and higher-d faces because they share nodes");
144 34699511 : if (_var.isNodalDefined())
145 : {
146 31981481 : _qp = 0;
147 31981481 : ComputeValueType value = computeValue();
148 : // update variable data, which is referenced by other kernels, so the value is up-to-date
149 31981477 : _var.setNodalValue(value);
150 : }
151 : }
152 : else /* elemental variables */
153 : {
154 17599811 : _n_local_dofs = _coincident_lower_d_calc ? _var.dofIndicesLower().size() : _var.numberOfDofs();
155 :
156 17599811 : if (_coincident_lower_d_calc)
157 : {
158 1926 : static const std::string lower_error = "Make sure that the lower-d variable lives on a "
159 : "lower-d block that is a superset of the boundary";
160 1862 : if (!_current_lower_d_elem)
161 5 : mooseError("No lower-dimensional element. ", lower_error);
162 1857 : if (!_n_local_dofs)
163 5 : mooseError("No degrees of freedom. ", lower_error);
164 : }
165 :
166 17599801 : if (_n_local_dofs == 1) /* p0 */
167 : {
168 16955136 : ComputeValueType value = 0;
169 100310607 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
170 83355485 : value += _JxW[_qp] * _coord[_qp] * computeValue();
171 16955122 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
172 16955122 : if (_var.isFV())
173 194724 : setDofValueHelper(value);
174 : else
175 : {
176 : // update the variable data referenced by other kernels.
177 : // Note that this will update the values at the quadrature points too
178 : // (because this is an Elemental variable)
179 16760398 : if (_coincident_lower_d_calc)
180 : {
181 926 : _local_sol.resize(1);
182 : if constexpr (std::is_same<Real, ComputeValueType>::value)
183 926 : _local_sol(0) = value;
184 : else
185 : mooseAssert(false, "We should not enter the single dof branch with a vector variable");
186 926 : _var.setLowerDofValues(_local_sol);
187 : }
188 : else
189 16759472 : _var.setNodalValue(value);
190 : }
191 : }
192 : else /* high-order */
193 : {
194 644665 : _local_re.resize(_n_local_dofs);
195 644665 : _local_re.zero();
196 644665 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
197 644665 : _local_ke.zero();
198 :
199 644665 : const auto & test = _coincident_lower_d_calc ? _var.phiLower() : _test;
200 :
201 : // assemble the local mass matrix and the load
202 3726397 : for (unsigned int i = 0; i < test.size(); i++)
203 28038675 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
204 : {
205 24956943 : ComputeValueType t = _JxW[_qp] * _coord[_qp] * test[i][_qp];
206 24956943 : _local_re(i) += t * computeValue();
207 195059120 : for (unsigned int j = 0; j < test.size(); j++)
208 170102182 : _local_ke(i, j) += t * test[j][_qp];
209 : }
210 : // mass matrix is always SPD but in case of boundary restricted, it will be rank deficient
211 644660 : _local_sol.resize(_n_local_dofs);
212 644660 : if (_bnd)
213 5678 : _local_ke.svd_solve(_local_re, _local_sol);
214 : else
215 638982 : _local_ke.cholesky_solve(_local_re, _local_sol);
216 :
217 644660 : _coincident_lower_d_calc ? _var.setLowerDofValues(_local_sol) : _var.setDofValues(_local_sol);
218 : }
219 : }
220 52299289 : }
221 :
222 : template <>
223 : void
224 662905 : AuxKernelTempl<RealEigenVector>::compute()
225 : {
226 662905 : precalculateValue();
227 :
228 662905 : if (isNodal()) /* nodal variables */
229 : {
230 657865 : if (_var.isNodalDefined())
231 : {
232 643465 : _qp = 0;
233 643465 : RealEigenVector value = computeValue();
234 : // update variable data, which is referenced by other kernels, so the value is up-to-date
235 643465 : _var.setNodalValue(value);
236 643465 : }
237 : }
238 : else /* elemental variables */
239 : {
240 5040 : _n_local_dofs = _var.numberOfDofs();
241 5040 : if (_n_local_dofs == 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 5040 : _local_re.resize(_n_local_dofs);
255 132696 : for (unsigned int i = 0; i < _local_re.size(); ++i)
256 127656 : _local_re(i) = RealEigenVector::Zero(_var.count());
257 5040 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
258 5040 : _local_ke.zero();
259 :
260 : // assemble the local mass matrix and the load
261 132696 : for (unsigned int i = 0; i < _test.size(); i++)
262 6307272 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
263 : {
264 6179616 : Real t = _JxW[_qp] * _coord[_qp] * _test[i][_qp];
265 6179616 : _local_re(i) += t * computeValue();
266 179072064 : for (unsigned int j = 0; j < _test.size(); j++)
267 172892448 : _local_ke(i, j) += t * _test[j][_qp];
268 : }
269 :
270 : // mass matrix is always SPD
271 5040 : _local_sol.resize(_n_local_dofs);
272 132696 : for (unsigned int i = 0; i < _local_re.size(); ++i)
273 127656 : _local_sol(i) = RealEigenVector::Zero(_var.count());
274 5040 : DenseVector<Number> re(_n_local_dofs);
275 5040 : DenseVector<Number> sol(_n_local_dofs);
276 19620 : for (unsigned int i = 0; i < _var.count(); ++i)
277 : {
278 395892 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
279 381312 : re(j) = _local_re(j)(i);
280 :
281 14580 : if (_bnd)
282 0 : _local_ke.svd_solve(re, sol);
283 : else
284 14580 : _local_ke.cholesky_solve(re, sol);
285 :
286 395892 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
287 381312 : _local_sol(j)(i) = sol(j);
288 : }
289 :
290 5040 : _var.setDofValues(_local_sol);
291 5040 : }
292 : }
293 662905 : }
294 :
295 : template <typename ComputeValueType>
296 : const typename OutputTools<ComputeValueType>::VariableValue &
297 126 : AuxKernelTempl<ComputeValueType>::uOld() const
298 : {
299 126 : 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 126 : return _nodal ? _var.nodalValueOldArray() : _var.slnOld();
306 : }
307 :
308 : template <typename ComputeValueType>
309 : const typename OutputTools<ComputeValueType>::VariableValue &
310 14 : AuxKernelTempl<ComputeValueType>::uOlder() const
311 : {
312 14 : 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 14 : return _nodal ? _var.nodalValueOlderArray() : _var.slnOlder();
319 : }
320 :
321 : template <typename ComputeValueType>
322 : bool
323 24115 : AuxKernelTempl<ComputeValueType>::isMortar()
324 : {
325 24115 : return dynamic_cast<MortarNodalAuxKernelTempl<ComputeValueType> *>(this) != nullptr;
326 : }
327 :
328 : // Explicitly instantiates the three versions of the AuxKernelTempl class
329 : template class AuxKernelTempl<Real>;
330 : template class AuxKernelTempl<RealVectorValue>;
331 : template class AuxKernelTempl<RealEigenVector>;
|
