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 2166942 : AuxKernelTempl<ComputeValueType>::validParams()
28 : {
29 2166942 : InputParameters params = MooseObject::validParams();
30 2166942 : params += BlockRestrictable::validParams();
31 2166942 : params += BoundaryRestrictable::validParams();
32 2166942 : params += RandomInterface::validParams();
33 2166942 : params += MeshChangedInterface::validParams();
34 2166942 : params += MaterialPropertyInterface::validParams();
35 2166942 : params += FunctorInterface::validParams();
36 :
37 : // Add the SetupInterface parameter 'execute_on' with 'linear' and 'timestep_end'
38 2166942 : params += SetupInterface::validParams();
39 2166942 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
40 2166942 : exec_enum.addAvailableFlags(EXEC_PRE_DISPLACE);
41 6500826 : exec_enum = {EXEC_LINEAR, EXEC_TIMESTEP_END};
42 2166942 : params.setDocString("execute_on", exec_enum.getDocString());
43 :
44 2166942 : params.addRequiredParam<AuxVariableName>("variable",
45 : "The name of the variable that this object applies to");
46 :
47 6500826 : params.addParam<bool>("use_displaced_mesh",
48 4333884 : false,
49 : "Whether or not this object should use the "
50 : "displaced mesh for computation. Note that "
51 : "in the case this is true but no "
52 : "displacements are provided in the Mesh block "
53 : "the undisplaced mesh will still be used.");
54 2166942 : params.addParamNamesToGroup("use_displaced_mesh", "Advanced");
55 6500826 : params.addParam<bool>("check_boundary_restricted",
56 4333884 : true,
57 : "Whether to check for multiple element sides on the boundary "
58 : "in the case of a boundary restricted, element aux variable. "
59 : "Setting this to false will allow contribution to a single element's "
60 : "elemental value(s) from multiple boundary sides on the same element "
61 : "(example: when the restricted boundary exists on two or more sides "
62 : "of an element, such as at a corner of a mesh");
63 :
64 : // This flag is set to true if the AuxKernelTempl is being used on a boundary
65 2166942 : params.addPrivateParam<bool>("_on_boundary", false);
66 :
67 2166942 : params.addRelationshipManager("GhostLowerDElems",
68 : Moose::RelationshipManagerType::GEOMETRIC |
69 : Moose::RelationshipManagerType::ALGEBRAIC);
70 :
71 2166942 : params.declareControllable("enable"); // allows Control to enable/disable this type of object
72 2166942 : params.registerBase("AuxKernel");
73 2166942 : params.registerSystemAttributeName("AuxKernel");
74 :
75 2166942 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(VectorAuxKernel).name())
76 100280 : params.registerBase("VectorAuxKernel");
77 2166942 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(ArrayAuxKernel).name())
78 114466 : params.registerBase("ArrayAuxKernel");
79 2166942 : return params;
80 2166942 : }
81 :
82 : template <typename ComputeValueType>
83 73515 : AuxKernelTempl<ComputeValueType>::AuxKernelTempl(const InputParameters & parameters)
84 : : MooseObject(parameters),
85 : MooseVariableInterface<ComputeValueType>(
86 : this,
87 : parameters.getCheckedPointerParam<SystemBase *>("_sys")
88 220541 : ->getVariable(parameters.get<THREAD_ID>("_tid"),
89 147026 : parameters.get<AuxVariableName>("variable"))
90 73511 : .isNodal(),
91 : "variable",
92 : Moose::VarKindType::VAR_AUXILIARY,
93 : std::is_same<Real, ComputeValueType>::value
94 : ? Moose::VarFieldType::VAR_FIELD_STANDARD
95 : : (std::is_same<RealVectorValue, ComputeValueType>::value
96 : ? Moose::VarFieldType::VAR_FIELD_VECTOR
97 : : Moose::VarFieldType::VAR_FIELD_ARRAY)),
98 : BlockRestrictable(this),
99 73511 : BoundaryRestrictable(this, mooseVariableBase()->isNodal()),
100 : SetupInterface(this),
101 73511 : CoupleableMooseVariableDependencyIntermediateInterface(this, mooseVariableBase()->isNodal()),
102 : FunctionInterface(this),
103 : UserObjectInterface(this),
104 : TransientInterface(this),
105 : MaterialPropertyInterface(this, blockIDs(), boundaryIDs()),
106 : PostprocessorInterface(this),
107 : DependencyResolverInterface(),
108 : RandomInterface(parameters,
109 147022 : *parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base"),
110 73511 : parameters.get<THREAD_ID>("_tid"),
111 73511 : mooseVariableBase()->isNodal()),
112 : GeometricSearchInterface(this),
113 : Restartable(this, "AuxKernels"),
114 : MeshChangedInterface(parameters),
115 : VectorPostprocessorInterface(this),
116 : ElementIDInterface(this),
117 : NonADFunctorInterface(this),
118 73511 : _check_boundary_restricted(getParam<bool>("check_boundary_restricted")),
119 73511 : _subproblem(*getCheckedPointerParam<SubProblem *>("_subproblem")),
120 73511 : _sys(*getCheckedPointerParam<SystemBase *>("_sys")),
121 73511 : _nl_sys(*getCheckedPointerParam<SystemBase *>("_nl_sys")),
122 73511 : _aux_sys(static_cast<AuxiliarySystem &>(_sys)),
123 73511 : _tid(parameters.get<THREAD_ID>("_tid")),
124 220533 : _var(_aux_sys.getActualFieldVariable<ComputeValueType>(
125 73511 : _tid, parameters.get<AuxVariableName>("variable"))),
126 73511 : _nodal(_var.isNodal()),
127 73511 : _u(_nodal ? _var.nodalValueArray() : _var.sln()),
128 :
129 73511 : _assembly(_subproblem.assembly(_tid, 0)),
130 73511 : _bnd(boundaryRestricted()),
131 73511 : _mesh(_subproblem.mesh()),
132 :
133 73511 : _test(_bnd ? _var.phiFace() : _var.phi()),
134 73511 : _q_point(_bnd ? _assembly.qPointsFace() : _assembly.qPoints()),
135 73511 : _qrule(_bnd ? _assembly.qRuleFace() : _assembly.qRule()),
136 73511 : _JxW(_bnd ? _assembly.JxWFace() : _assembly.JxW()),
137 73511 : _coord(_assembly.coordTransformation()),
138 :
139 73511 : _current_elem(_assembly.elem()),
140 73511 : _current_side(_assembly.side()),
141 73511 : _current_elem_volume(_assembly.elemVolume()),
142 73511 : _current_side_volume(_assembly.sideElemVolume()),
143 :
144 73511 : _current_node(_assembly.node()),
145 73511 : _current_boundary_id(_assembly.currentBoundaryID()),
146 73511 : _solution(_aux_sys.solution()),
147 :
148 73511 : _current_lower_d_elem(_assembly.lowerDElem()),
149 735114 : _coincident_lower_d_calc(_bnd && !isNodal() && _var.isLowerD())
150 : {
151 73511 : addMooseVariableDependency(&_var);
152 73511 : _supplied_vars.insert(parameters.get<AuxVariableName>("variable"));
153 :
154 73511 : if (_bnd && !isNodal() && !_coincident_lower_d_calc && _check_boundary_restricted)
155 : {
156 : // when the variable is elemental and this aux kernel operates on boundaries,
157 : // we need to check that no elements are visited more than once through visiting
158 : // all the sides on the boundaries
159 1310 : auto boundaries = _mesh.getMesh().get_boundary_info().build_side_list();
160 1310 : std::set<dof_id_type> elements;
161 194805 : for (const auto & t : boundaries)
162 : {
163 193499 : if (hasBoundary(std::get<2>(t)))
164 : {
165 72459 : const auto eid = std::get<0>(t);
166 72459 : const auto stat = elements.insert(eid);
167 72459 : if (!stat.second) // already existed in the set
168 4 : mooseError(
169 : "Boundary restricted auxiliary kernel '",
170 4 : name(),
171 : "' has element (id=",
172 : eid,
173 : ") connected with more than one boundary sides.\nTo skip this error check, "
174 : "set 'check_boundary_restricted = false'.\nRefer to the AuxKernel "
175 : "documentation on boundary restricted aux kernels for understanding this error.");
176 : }
177 : }
178 1306 : }
179 :
180 : // Check for supported variable types
181 : // Any 'nodal' family that actually has DoFs outside of nodes, or gradient dofs at nodes is
182 : // not properly set by AuxKernelTempl::compute
183 : // NOTE: We could add a few exceptions, lower order from certain unsupported families and on
184 : // certain element types only have value-DoFs on nodes
185 73507 : const auto type = _var.feType();
186 73507 : if (_var.isNodal() && !((type.family == LAGRANGE) || (type.order <= FIRST)))
187 0 : paramError("variable",
188 0 : "Variable family " + Moose::stringify(type.family) + " is not supported at order " +
189 : Moose::stringify(type.order) + " by the AuxKernel system.");
190 73507 : }
191 :
192 : template <typename ComputeValueType>
193 : const std::set<std::string> &
194 239398 : AuxKernelTempl<ComputeValueType>::getRequestedItems()
195 : {
196 239398 : return _depend_vars;
197 : }
198 :
199 : template <typename ComputeValueType>
200 : const std::set<std::string> &
201 239398 : AuxKernelTempl<ComputeValueType>::getSuppliedItems()
202 : {
203 239398 : return _supplied_vars;
204 : }
205 :
206 : template <typename ComputeValueType>
207 : void
208 5398 : AuxKernelTempl<ComputeValueType>::addUserObjectDependencyHelper(const UserObject & uo) const
209 : {
210 5398 : _depend_uo.insert(uo.name());
211 5478 : for (const auto & indirect_dependent : uo.getDependObjects())
212 80 : _depend_uo.insert(indirect_dependent);
213 5398 : }
214 :
215 : template <typename ComputeValueType>
216 : void
217 223 : AuxKernelTempl<ComputeValueType>::addPostprocessorDependencyHelper(
218 : const PostprocessorName & name) const
219 : {
220 223 : getUserObjectBaseByName(name); // getting the UO will call addUserObjectDependencyHelper()
221 223 : }
222 :
223 : template <typename ComputeValueType>
224 : void
225 250 : AuxKernelTempl<ComputeValueType>::addVectorPostprocessorDependencyHelper(
226 : const VectorPostprocessorName & name) const
227 : {
228 250 : getUserObjectBaseByName(name); // getting the UO will call addUserObjectDependencyHelper()
229 250 : }
230 :
231 : template <typename ComputeValueType>
232 : void
233 8936 : AuxKernelTempl<ComputeValueType>::coupledCallback(const std::string & var_name, bool is_old) const
234 : {
235 8936 : if (!is_old)
236 : {
237 8730 : const auto & var_names = getParam<std::vector<VariableName>>(var_name);
238 8730 : _depend_vars.insert(var_names.begin(), var_names.end());
239 : }
240 8936 : }
241 :
242 : template <typename ComputeValueType>
243 : const VariableValue &
244 74 : AuxKernelTempl<ComputeValueType>::coupledDot(const std::string & var_name, unsigned int comp) const
245 : {
246 74 : auto var = getVar(var_name, comp);
247 74 : if (var->kind() == Moose::VAR_AUXILIARY)
248 8 : mooseError(
249 4 : name(),
250 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
251 :
252 70 : return Coupleable::coupledDot(var_name, comp);
253 : }
254 :
255 : template <typename ComputeValueType>
256 : const VariableValue &
257 0 : AuxKernelTempl<ComputeValueType>::coupledDotDu(const std::string & var_name,
258 : unsigned int comp) const
259 : {
260 0 : auto var = getVar(var_name, comp);
261 0 : if (var->kind() == Moose::VAR_AUXILIARY)
262 0 : mooseError(
263 0 : name(),
264 : ": Unable to couple time derivative of an auxiliary variable into the auxiliary system.");
265 :
266 0 : return Coupleable::coupledDotDu(var_name, comp);
267 : }
268 :
269 : template <>
270 : void
271 194724 : AuxKernelTempl<Real>::setDofValueHelper(const Real & value)
272 : {
273 : mooseAssert(_n_local_dofs == 1,
274 : "Should only be calling setDofValue if there is one dof for the aux var");
275 194724 : _var.setDofValue(value, 0);
276 194724 : }
277 :
278 : template <>
279 : void
280 0 : AuxKernelTempl<RealVectorValue>::setDofValueHelper(const RealVectorValue &)
281 : {
282 0 : mooseError("Not implemented");
283 : }
284 :
285 : template <typename ComputeValueType>
286 : void
287 132721 : AuxKernelTempl<ComputeValueType>::insert()
288 : {
289 132721 : if (_coincident_lower_d_calc)
290 1848 : _var.insertLower(_aux_sys.solution());
291 : else
292 130873 : _var.insert(_aux_sys.solution());
293 132721 : }
294 :
295 : template <typename ComputeValueType>
296 : void
297 51806531 : AuxKernelTempl<ComputeValueType>::compute()
298 : {
299 51806531 : precalculateValue();
300 :
301 51806531 : if (isNodal()) /* nodal variables */
302 : {
303 : mooseAssert(!_coincident_lower_d_calc,
304 : "Nodal evaluations are point evaluations. We don't have to concern ourselves with "
305 : "coincidence of lower-d blocks and higher-d faces because they share nodes");
306 34592287 : if (_var.isNodalDefined())
307 : {
308 31880530 : _qp = 0;
309 31880530 : ComputeValueType value = computeValue();
310 : // update variable data, which is referenced by other kernels, so the value is up-to-date
311 31880526 : _var.setNodalValue(value);
312 : }
313 : }
314 : else /* elemental variables */
315 : {
316 17214244 : _n_local_dofs = _coincident_lower_d_calc ? _var.dofIndicesLower().size() : _var.numberOfDofs();
317 :
318 17214244 : if (_coincident_lower_d_calc)
319 : {
320 1854 : static const std::string lower_error = "Make sure that the lower-d variable lives on a "
321 : "lower-d block that is a superset of the boundary";
322 1854 : if (!_current_lower_d_elem)
323 3 : mooseError("No lower-dimensional element. ", lower_error);
324 1851 : if (!_n_local_dofs)
325 3 : mooseError("No degrees of freedom. ", lower_error);
326 : }
327 :
328 17214238 : if (_n_local_dofs == 1) /* p0 */
329 : {
330 16569577 : ComputeValueType value = 0;
331 98391266 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
332 81821701 : value += _JxW[_qp] * _coord[_qp] * computeValue();
333 16569565 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
334 16569565 : if (_var.isFV())
335 194724 : setDofValueHelper(value);
336 : else
337 : {
338 : // update the variable data referenced by other kernels.
339 : // Note that this will update the values at the quadrature points too
340 : // (because this is an Elemental variable)
341 16374841 : if (_coincident_lower_d_calc)
342 : {
343 924 : _local_sol.resize(1);
344 : if constexpr (std::is_same<Real, ComputeValueType>::value)
345 924 : _local_sol(0) = value;
346 : else
347 : mooseAssert(false, "We should not enter the single dof branch with a vector variable");
348 924 : _var.setLowerDofValues(_local_sol);
349 : }
350 : else
351 16373917 : _var.setNodalValue(value);
352 : }
353 : }
354 : else /* high-order */
355 : {
356 644661 : _local_re.resize(_n_local_dofs);
357 644661 : _local_re.zero();
358 644661 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
359 644661 : _local_ke.zero();
360 :
361 644661 : const auto & test = _coincident_lower_d_calc ? _var.phiLower() : _test;
362 :
363 : // assemble the local mass matrix and the load
364 3726389 : for (unsigned int i = 0; i < test.size(); i++)
365 28038661 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
366 : {
367 24956933 : ComputeValueType t = _JxW[_qp] * _coord[_qp] * test[i][_qp];
368 24956933 : _local_re(i) += t * computeValue();
369 195059096 : for (unsigned int j = 0; j < test.size(); j++)
370 170102166 : _local_ke(i, j) += t * test[j][_qp];
371 : }
372 : // mass matrix is always SPD but in case of boundary restricted, it will be rank deficient
373 644658 : _local_sol.resize(_n_local_dofs);
374 644658 : if (_bnd)
375 5676 : _local_ke.svd_solve(_local_re, _local_sol);
376 : else
377 638982 : _local_ke.cholesky_solve(_local_re, _local_sol);
378 :
379 644658 : _coincident_lower_d_calc ? _var.setLowerDofValues(_local_sol) : _var.setDofValues(_local_sol);
380 : }
381 : }
382 51806506 : }
383 :
384 : template <>
385 : void
386 662905 : AuxKernelTempl<RealEigenVector>::compute()
387 : {
388 662905 : precalculateValue();
389 :
390 662905 : if (isNodal()) /* nodal variables */
391 : {
392 657865 : if (_var.isNodalDefined())
393 : {
394 643465 : _qp = 0;
395 643465 : RealEigenVector value = computeValue();
396 : // update variable data, which is referenced by other kernels, so the value is up-to-date
397 643465 : _var.setNodalValue(value);
398 643465 : }
399 : }
400 : else /* elemental variables */
401 : {
402 5040 : _n_local_dofs = _var.numberOfDofs();
403 5040 : if (_n_local_dofs == 1) /* p0 */
404 : {
405 0 : RealEigenVector value = RealEigenVector::Zero(_var.count());
406 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
407 0 : value += _JxW[_qp] * _coord[_qp] * computeValue();
408 0 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
409 : // update the variable data referenced by other kernels.
410 : // Note that this will update the values at the quadrature points too
411 : // (because this is an Elemental variable)
412 0 : _var.setNodalValue(value);
413 0 : }
414 : else /* high-order */
415 : {
416 5040 : _local_re.resize(_n_local_dofs);
417 132696 : for (unsigned int i = 0; i < _local_re.size(); ++i)
418 127656 : _local_re(i) = RealEigenVector::Zero(_var.count());
419 5040 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
420 5040 : _local_ke.zero();
421 :
422 : // assemble the local mass matrix and the load
423 132696 : for (unsigned int i = 0; i < _test.size(); i++)
424 6307272 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
425 : {
426 6179616 : Real t = _JxW[_qp] * _coord[_qp] * _test[i][_qp];
427 6179616 : _local_re(i) += t * computeValue();
428 179072064 : for (unsigned int j = 0; j < _test.size(); j++)
429 172892448 : _local_ke(i, j) += t * _test[j][_qp];
430 : }
431 :
432 : // mass matrix is always SPD
433 5040 : _local_sol.resize(_n_local_dofs);
434 132696 : for (unsigned int i = 0; i < _local_re.size(); ++i)
435 127656 : _local_sol(i) = RealEigenVector::Zero(_var.count());
436 5040 : DenseVector<Number> re(_n_local_dofs);
437 5040 : DenseVector<Number> sol(_n_local_dofs);
438 19620 : for (unsigned int i = 0; i < _var.count(); ++i)
439 : {
440 395892 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
441 381312 : re(j) = _local_re(j)(i);
442 :
443 14580 : if (_bnd)
444 0 : _local_ke.svd_solve(re, sol);
445 : else
446 14580 : _local_ke.cholesky_solve(re, sol);
447 :
448 395892 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
449 381312 : _local_sol(j)(i) = sol(j);
450 : }
451 :
452 5040 : _var.setDofValues(_local_sol);
453 5040 : }
454 : }
455 662905 : }
456 :
457 : template <typename ComputeValueType>
458 : const typename OutputTools<ComputeValueType>::VariableValue &
459 126 : AuxKernelTempl<ComputeValueType>::uOld() const
460 : {
461 126 : if (_sys.solutionStatesInitialized())
462 0 : mooseError("The solution states have already been initialized when calling ",
463 0 : type(),
464 : "::uOld().\n\n",
465 : "Make sure to call uOld() within the object constructor.");
466 :
467 126 : return _nodal ? _var.nodalValueOldArray() : _var.slnOld();
468 : }
469 :
470 : template <typename ComputeValueType>
471 : const typename OutputTools<ComputeValueType>::VariableValue &
472 14 : AuxKernelTempl<ComputeValueType>::uOlder() const
473 : {
474 14 : if (_sys.solutionStatesInitialized())
475 0 : mooseError("The solution states have already been initialized when calling ",
476 0 : type(),
477 : "::uOlder().\n\n",
478 : "Make sure to call uOlder() within the object constructor.");
479 :
480 14 : return _nodal ? _var.nodalValueOlderArray() : _var.slnOlder();
481 : }
482 :
483 : template <typename ComputeValueType>
484 : bool
485 23621 : AuxKernelTempl<ComputeValueType>::isMortar()
486 : {
487 23621 : return dynamic_cast<MortarNodalAuxKernelTempl<ComputeValueType> *>(this) != nullptr;
488 : }
489 :
490 : // Explicitly instantiates the three versions of the AuxKernelTempl class
491 : template class AuxKernelTempl<Real>;
492 : template class AuxKernelTempl<RealVectorValue>;
493 : template class AuxKernelTempl<RealEigenVector>;
|
