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