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 2157238 : AuxKernelTempl<ComputeValueType>::validParams()
28 : {
29 2157238 : InputParameters params = MooseObject::validParams();
30 2157238 : params += BlockRestrictable::validParams();
31 2157238 : params += BoundaryRestrictable::validParams();
32 2157238 : params += RandomInterface::validParams();
33 2157238 : params += MeshChangedInterface::validParams();
34 2157238 : params += MaterialPropertyInterface::validParams();
35 2157238 : params += FunctorInterface::validParams();
36 :
37 : // Add the SetupInterface parameter 'execute_on' with 'linear' and 'timestep_end'
38 2157238 : params += SetupInterface::validParams();
39 2157238 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
40 2157238 : exec_enum.addAvailableFlags(EXEC_PRE_DISPLACE);
41 6471714 : exec_enum = {EXEC_LINEAR, EXEC_TIMESTEP_END};
42 2157238 : params.setDocString("execute_on", exec_enum.getDocString());
43 :
44 2157238 : params.addRequiredParam<AuxVariableName>("variable",
45 : "The name of the variable that this object applies to");
46 :
47 6471714 : params.addParam<bool>("use_displaced_mesh",
48 4314476 : 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 2157238 : params.addParamNamesToGroup("use_displaced_mesh", "Advanced");
55 6471714 : params.addParam<bool>("check_boundary_restricted",
56 4314476 : 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 2157238 : params.addPrivateParam<bool>("_on_boundary", false);
66 :
67 2157238 : params.addRelationshipManager("GhostLowerDElems",
68 : Moose::RelationshipManagerType::GEOMETRIC |
69 : Moose::RelationshipManagerType::ALGEBRAIC);
70 :
71 2157238 : params.declareControllable("enable"); // allows Control to enable/disable this type of object
72 2157238 : params.registerBase("AuxKernel");
73 2157238 : params.registerSystemAttributeName("AuxKernel");
74 :
75 2157238 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(VectorAuxKernel).name())
76 100248 : params.registerBase("VectorAuxKernel");
77 2157238 : if (typeid(AuxKernelTempl<ComputeValueType>).name() == typeid(ArrayAuxKernel).name())
78 114442 : params.registerBase("ArrayAuxKernel");
79 2157238 : return params;
80 2157238 : }
81 :
82 : template <typename ComputeValueType>
83 68664 : AuxKernelTempl<ComputeValueType>::AuxKernelTempl(const InputParameters & parameters)
84 : : MooseObject(parameters),
85 : MooseVariableInterface<ComputeValueType>(
86 : this,
87 : parameters.getCheckedPointerParam<SystemBase *>("_sys")
88 205988 : ->getVariable(parameters.get<THREAD_ID>("_tid"),
89 137324 : parameters.get<AuxVariableName>("variable"))
90 68660 : .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 68660 : BoundaryRestrictable(this, mooseVariableBase()->isNodal()),
100 : SetupInterface(this),
101 68660 : 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 137320 : *parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base"),
110 68660 : parameters.get<THREAD_ID>("_tid"),
111 68660 : mooseVariableBase()->isNodal()),
112 : GeometricSearchInterface(this),
113 : Restartable(this, "AuxKernels"),
114 : MeshChangedInterface(parameters),
115 : VectorPostprocessorInterface(this),
116 : ElementIDInterface(this),
117 : NonADFunctorInterface(this),
118 68660 : _check_boundary_restricted(getParam<bool>("check_boundary_restricted")),
119 68660 : _subproblem(*getCheckedPointerParam<SubProblem *>("_subproblem")),
120 68660 : _sys(*getCheckedPointerParam<SystemBase *>("_sys")),
121 68660 : _nl_sys(*getCheckedPointerParam<SystemBase *>("_nl_sys")),
122 68660 : _aux_sys(static_cast<AuxiliarySystem &>(_sys)),
123 68660 : _tid(parameters.get<THREAD_ID>("_tid")),
124 205980 : _var(_aux_sys.getActualFieldVariable<ComputeValueType>(
125 68660 : _tid, parameters.get<AuxVariableName>("variable"))),
126 68660 : _nodal(_var.isNodal()),
127 68660 : _u(_nodal ? _var.nodalValueArray() : _var.sln()),
128 :
129 68660 : _assembly(_subproblem.assembly(_tid, 0)),
130 68660 : _bnd(boundaryRestricted()),
131 68660 : _mesh(_subproblem.mesh()),
132 :
133 68660 : _test(_bnd ? _var.phiFace() : _var.phi()),
134 68660 : _q_point(_bnd ? _assembly.qPointsFace() : _assembly.qPoints()),
135 68660 : _qrule(_bnd ? _assembly.qRuleFace() : _assembly.qRule()),
136 68660 : _JxW(_bnd ? _assembly.JxWFace() : _assembly.JxW()),
137 68660 : _coord(_assembly.coordTransformation()),
138 :
139 68660 : _current_elem(_assembly.elem()),
140 68660 : _current_side(_assembly.side()),
141 68660 : _current_elem_volume(_assembly.elemVolume()),
142 68660 : _current_side_volume(_assembly.sideElemVolume()),
143 :
144 68660 : _current_node(_assembly.node()),
145 68660 : _current_boundary_id(_assembly.currentBoundaryID()),
146 68660 : _solution(_aux_sys.solution()),
147 :
148 68660 : _current_lower_d_elem(_assembly.lowerDElem()),
149 686604 : _coincident_lower_d_calc(_bnd && !isNodal() && _var.isLowerD())
150 : {
151 68660 : addMooseVariableDependency(&_var);
152 68660 : _supplied_vars.insert(parameters.get<AuxVariableName>("variable"));
153 :
154 68660 : 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 1214 : auto boundaries = _mesh.getMesh().get_boundary_info().build_side_list();
160 1214 : std::set<dof_id_type> elements;
161 180823 : for (const auto & t : boundaries)
162 : {
163 179613 : if (hasBoundary(std::get<2>(t)))
164 : {
165 67219 : const auto eid = std::get<0>(t);
166 67219 : const auto stat = elements.insert(eid);
167 67219 : 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 1210 : }
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 68656 : const auto type = _var.feType();
186 68656 : 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 68656 : }
191 :
192 : template <typename ComputeValueType>
193 : const std::set<std::string> &
194 223387 : AuxKernelTempl<ComputeValueType>::getRequestedItems()
195 : {
196 223387 : return _depend_vars;
197 : }
198 :
199 : template <typename ComputeValueType>
200 : const std::set<std::string> &
201 223387 : AuxKernelTempl<ComputeValueType>::getSuppliedItems()
202 : {
203 223387 : return _supplied_vars;
204 : }
205 :
206 : template <typename ComputeValueType>
207 : void
208 4936 : AuxKernelTempl<ComputeValueType>::addUserObjectDependencyHelper(const UserObject & uo) const
209 : {
210 4936 : _depend_uo.insert(uo.name());
211 5011 : for (const auto & indirect_dependent : uo.getDependObjects())
212 75 : _depend_uo.insert(indirect_dependent);
213 4936 : }
214 :
215 : template <typename ComputeValueType>
216 : void
217 210 : AuxKernelTempl<ComputeValueType>::addPostprocessorDependencyHelper(
218 : const PostprocessorName & name) const
219 : {
220 210 : getUserObjectBaseByName(name); // getting the UO will call addUserObjectDependencyHelper()
221 210 : }
222 :
223 : template <typename ComputeValueType>
224 : void
225 227 : AuxKernelTempl<ComputeValueType>::addVectorPostprocessorDependencyHelper(
226 : const VectorPostprocessorName & name) const
227 : {
228 227 : getUserObjectBaseByName(name); // getting the UO will call addUserObjectDependencyHelper()
229 227 : }
230 :
231 : template <typename ComputeValueType>
232 : void
233 8317 : AuxKernelTempl<ComputeValueType>::coupledCallback(const std::string & var_name, bool is_old) const
234 : {
235 8317 : if (!is_old)
236 : {
237 8126 : const auto & var_names = getParam<std::vector<VariableName>>(var_name);
238 8126 : _depend_vars.insert(var_names.begin(), var_names.end());
239 : }
240 8317 : }
241 :
242 : template <typename ComputeValueType>
243 : const VariableValue &
244 69 : AuxKernelTempl<ComputeValueType>::coupledDot(const std::string & var_name, unsigned int comp) const
245 : {
246 69 : auto var = getVar(var_name, comp);
247 69 : 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 65 : 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 172480 : 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 172480 : _var.setDofValue(value, 0);
276 172480 : }
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 118000 : AuxKernelTempl<ComputeValueType>::insert()
288 : {
289 118000 : if (_coincident_lower_d_calc)
290 1644 : _var.insertLower(_aux_sys.solution());
291 : else
292 116356 : _var.insert(_aux_sys.solution());
293 118000 : }
294 :
295 : template <typename ComputeValueType>
296 : void
297 45960456 : AuxKernelTempl<ComputeValueType>::compute()
298 : {
299 45960456 : precalculateValue();
300 :
301 45960456 : 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 30708298 : if (_var.isNodalDefined())
307 : {
308 28314274 : _qp = 0;
309 28314274 : ComputeValueType value = computeValue();
310 : // update variable data, which is referenced by other kernels, so the value is up-to-date
311 28314270 : _var.setNodalValue(value);
312 : }
313 : }
314 : else /* elemental variables */
315 : {
316 15252158 : _n_local_dofs = _coincident_lower_d_calc ? _var.dofIndicesLower().size() : _var.numberOfDofs();
317 :
318 15252158 : if (_coincident_lower_d_calc)
319 : {
320 1650 : 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 1650 : if (!_current_lower_d_elem)
323 3 : mooseError("No lower-dimensional element. ", lower_error);
324 1647 : if (!_n_local_dofs)
325 3 : mooseError("No degrees of freedom. ", lower_error);
326 : }
327 :
328 15252152 : if (_n_local_dofs == 1) /* p0 */
329 : {
330 14684475 : ComputeValueType value = 0;
331 86967800 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
332 72283337 : value += _JxW[_qp] * _coord[_qp] * computeValue();
333 14684463 : value /= (_bnd ? _current_side_volume : _current_elem_volume);
334 14684463 : if (_var.isFV())
335 172480 : 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 14511983 : if (_coincident_lower_d_calc)
342 : {
343 822 : _local_sol.resize(1);
344 : if constexpr (std::is_same<Real, ComputeValueType>::value)
345 822 : _local_sol(0) = value;
346 : else
347 : mooseAssert(false, "We should not enter the single dof branch with a vector variable");
348 822 : _var.setLowerDofValues(_local_sol);
349 : }
350 : else
351 14511161 : _var.setNodalValue(value);
352 : }
353 : }
354 : else /* high-order */
355 : {
356 567677 : _local_re.resize(_n_local_dofs);
357 567677 : _local_re.zero();
358 567677 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
359 567677 : _local_ke.zero();
360 :
361 567677 : const auto & test = _coincident_lower_d_calc ? _var.phiLower() : _test;
362 :
363 : // assemble the local mass matrix and the load
364 3278953 : for (unsigned int i = 0; i < test.size(); i++)
365 24756903 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
366 : {
367 22045627 : ComputeValueType t = _JxW[_qp] * _coord[_qp] * test[i][_qp];
368 22045627 : _local_re(i) += t * computeValue();
369 172492584 : for (unsigned int j = 0; j < test.size(); j++)
370 150446960 : _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 567674 : _local_sol.resize(_n_local_dofs);
374 567674 : if (_bnd)
375 4854 : _local_ke.svd_solve(_local_re, _local_sol);
376 : else
377 562820 : _local_ke.cholesky_solve(_local_re, _local_sol);
378 :
379 567674 : _coincident_lower_d_calc ? _var.setLowerDofValues(_local_sol) : _var.setDofValues(_local_sol);
380 : }
381 : }
382 45960431 : }
383 :
384 : template <>
385 : void
386 575344 : AuxKernelTempl<RealEigenVector>::compute()
387 : {
388 575344 : precalculateValue();
389 :
390 575344 : if (isNodal()) /* nodal variables */
391 : {
392 570872 : if (_var.isNodalDefined())
393 : {
394 558072 : _qp = 0;
395 558072 : RealEigenVector value = computeValue();
396 : // update variable data, which is referenced by other kernels, so the value is up-to-date
397 558072 : _var.setNodalValue(value);
398 558072 : }
399 : }
400 : else /* elemental variables */
401 : {
402 4472 : _n_local_dofs = _var.numberOfDofs();
403 4472 : 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 4472 : _local_re.resize(_n_local_dofs);
417 117928 : for (unsigned int i = 0; i < _local_re.size(); ++i)
418 113456 : _local_re(i) = RealEigenVector::Zero(_var.count());
419 4472 : _local_ke.resize(_n_local_dofs, _n_local_dofs);
420 4472 : _local_ke.zero();
421 :
422 : // assemble the local mass matrix and the load
423 117928 : for (unsigned int i = 0; i < _test.size(); i++)
424 5606416 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
425 : {
426 5492960 : Real t = _JxW[_qp] * _coord[_qp] * _test[i][_qp];
427 5492960 : _local_re(i) += t * computeValue();
428 159175072 : for (unsigned int j = 0; j < _test.size(); j++)
429 153682112 : _local_ke(i, j) += t * _test[j][_qp];
430 : }
431 :
432 : // mass matrix is always SPD
433 4472 : _local_sol.resize(_n_local_dofs);
434 117928 : for (unsigned int i = 0; i < _local_re.size(); ++i)
435 113456 : _local_sol(i) = RealEigenVector::Zero(_var.count());
436 4472 : DenseVector<Number> re(_n_local_dofs);
437 4472 : DenseVector<Number> sol(_n_local_dofs);
438 17416 : for (unsigned int i = 0; i < _var.count(); ++i)
439 : {
440 351856 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
441 338912 : re(j) = _local_re(j)(i);
442 :
443 12944 : if (_bnd)
444 0 : _local_ke.svd_solve(re, sol);
445 : else
446 12944 : _local_ke.cholesky_solve(re, sol);
447 :
448 351856 : for (unsigned int j = 0; j < _n_local_dofs; ++j)
449 338912 : _local_sol(j)(i) = sol(j);
450 : }
451 :
452 4472 : _var.setDofValues(_local_sol);
453 4472 : }
454 : }
455 575344 : }
456 :
457 : template <typename ComputeValueType>
458 : const typename OutputTools<ComputeValueType>::VariableValue &
459 117 : AuxKernelTempl<ComputeValueType>::uOld() const
460 : {
461 117 : 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 117 : return _nodal ? _var.nodalValueOldArray() : _var.slnOld();
468 : }
469 :
470 : template <typename ComputeValueType>
471 : const typename OutputTools<ComputeValueType>::VariableValue &
472 13 : AuxKernelTempl<ComputeValueType>::uOlder() const
473 : {
474 13 : 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 13 : return _nodal ? _var.nodalValueOlderArray() : _var.slnOlder();
481 : }
482 :
483 : template <typename ComputeValueType>
484 : bool
485 21932 : AuxKernelTempl<ComputeValueType>::isMortar()
486 : {
487 21932 : 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>;
|
