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