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 : #pragma once
11 :
12 : #include "DenseMatrix.h"
13 : #include "MooseTypes.h"
14 : #include "MultiMooseEnum.h"
15 : #include "Assembly.h"
16 : #include "MooseVariableFE.h"
17 : #include "SystemBase.h"
18 :
19 : #include "libmesh/dense_vector.h"
20 : #include "metaphysicl/raw_type.h"
21 :
22 : #include <vector>
23 :
24 : // Forward declarations
25 : class InputParameters;
26 : class MooseObject;
27 : class SubProblem;
28 : class Assembly;
29 :
30 : template <typename T>
31 : InputParameters validParams();
32 :
33 : /**
34 : * Utility structure for packaging up all of the residual object's information needed to add into
35 : * the system residual and Jacobian in the context of automatic differentiation. The necessary
36 : * information includes the vector of residuals and Jacobians represented by dual numbers, the
37 : * vector of degrees of freedom (this should be the same length as the vector of dual numbers), and
38 : * a scaling factor that will multiply the dual numbers before their addition into the global
39 : * residual and Jacobian
40 : */
41 : struct ADResidualsPacket
42 : {
43 : const DenseVector<ADReal> & residuals;
44 : const std::vector<dof_id_type> & dof_indices;
45 : const Real scaling_factor;
46 : };
47 :
48 : class TaggingInterface
49 : {
50 : public:
51 : TaggingInterface(const MooseObject * moose_object);
52 : virtual ~TaggingInterface();
53 :
54 : static InputParameters validParams();
55 :
56 : /**
57 : * Class that is used as a parameter to some of our vector tag APIs that allows only
58 : * blessed framework classes to call them
59 : */
60 : class VectorTagsKey
61 : {
62 : friend class AttribVectorTags;
63 : friend class NonlinearEigenSystem;
64 : friend class LinearSystemContributionObject;
65 : template <typename>
66 : friend class MooseObjectTagWarehouse;
67 :
68 275588 : VectorTagsKey() {}
69 : VectorTagsKey(const VectorTagsKey &) {}
70 : };
71 :
72 : /**
73 : * Class that is used as a parameter to some of our matrix tag APIs that allows only
74 : * blessed framework classes to call them
75 : */
76 : class MatrixTagsKey
77 : {
78 : friend class AttribMatrixTags;
79 : friend class NonlinearEigenSystem;
80 : friend class LinearSystemContributionObject;
81 : template <typename>
82 : friend class MooseObjectTagWarehouse;
83 :
84 237879 : MatrixTagsKey() {}
85 : MatrixTagsKey(const MatrixTagsKey &) {}
86 : };
87 :
88 : /**
89 : * Enumerate whether a (residual) vector tag is to be of a non-reference or reference tag type
90 : */
91 : enum class ResidualTagType
92 : {
93 : NonReference,
94 : Reference
95 : };
96 :
97 : void useVectorTag(const TagName & tag_name, VectorTagsKey);
98 :
99 : void useMatrixTag(const TagName & tag_name, MatrixTagsKey);
100 :
101 : void useVectorTag(TagID tag_id, VectorTagsKey);
102 :
103 : void useMatrixTag(TagID tag_id, MatrixTagsKey);
104 :
105 : bool isVectorTagged() { return _vector_tags.size() > 0; }
106 :
107 : bool isMatrixTagged() { return _matrix_tags.size() > 0; }
108 :
109 488748355 : bool hasVectorTags() const { return !_vector_tags.empty(); }
110 :
111 272255 : const std::set<TagID> & getVectorTags(VectorTagsKey) const { return _vector_tags; }
112 :
113 232486 : const std::set<TagID> & getMatrixTags(MatrixTagsKey) const { return _matrix_tags; }
114 :
115 : protected:
116 : /**
117 : * Prepare data for computing element residual according to active tags.
118 : * Residual blocks for different tags will be extracted from Assembly.
119 : * A local residual will be zeroed. It should be called
120 : * right before the local element vector is computed.
121 : */
122 : void prepareVectorTag(Assembly & assembly, unsigned int ivar);
123 :
124 : /**
125 : * Prepare vector tags in a reference residual problem context
126 : * @param Assembly The assembly object that we obtain the local residual blocks from
127 : * @param ivar The variable which we are retrieving the local residual blocks for
128 : * @param ref_problem A pointer to a reference residual problem. This can be a nullptr
129 : * @param tag_type What type of tags to prepare
130 : */
131 : void prepareVectorTag(Assembly & assembly, unsigned int ivar, ResidualTagType tag_type);
132 :
133 : /**
134 : * Prepare data for computing element residual the according to active tags
135 : * for DG and interface kernels.
136 : * Residual blocks for different tags will be extracted from Assembly.
137 : * A local residual will be zeroed. It should be called
138 : * right before the local element vector is computed.
139 : */
140 : void prepareVectorTagNeighbor(Assembly & assembly, unsigned int ivar);
141 :
142 : /**
143 : * Prepare data for computing the residual according to active tags for mortar constraints.
144 : * Residual blocks for different tags will be extracted from Assembly. A local residual will be
145 : * zeroed. It should be called right before the local element vector is computed.
146 : */
147 : void prepareVectorTagLower(Assembly & assembly, unsigned int ivar);
148 :
149 : /**
150 : * Prepare data for computing element jacobian according to the active tags.
151 : * Jacobian blocks for different tags will be extracted from Assembly.
152 : * A local Jacobian will be zeroed. It should be called
153 : * right before the local element matrix is computed.
154 : */
155 : void prepareMatrixTag(Assembly & assembly, unsigned int ivar, unsigned int jvar);
156 :
157 : void prepareMatrixTag(Assembly & assembly,
158 : unsigned int ivar,
159 : unsigned int jvar,
160 : DenseMatrix<Number> & k) const;
161 :
162 : /**
163 : * Prepare data for computing nonlocal element jacobian according to the active tags.
164 : * Jacobian blocks for different tags will be extracted from Assembly.
165 : * A nonlocal Jacobian will be zeroed. It should be called
166 : * right before the nonlocal element matrix is computed.
167 : */
168 : void prepareMatrixTagNonlocal(Assembly & assembly, unsigned int ivar, unsigned int jvar);
169 :
170 : /**
171 : * Prepare data for computing element jacobian according to the active tags
172 : * for DG and interface kernels.
173 : * Jacobian blocks for different tags will be extracted from Assembly.
174 : * A local Jacobian will be zeroed. It should be called
175 : * right before the local element matrix is computed.
176 : */
177 : void prepareMatrixTagNeighbor(Assembly & assembly,
178 : unsigned int ivar,
179 : unsigned int jvar,
180 : Moose::DGJacobianType type);
181 :
182 : void prepareMatrixTagNeighbor(Assembly & assembly,
183 : unsigned int ivar,
184 : unsigned int jvar,
185 : Moose::DGJacobianType type,
186 : DenseMatrix<Number> & k) const;
187 :
188 : /**
189 : * Prepare data for computing the jacobian according to the active tags for mortar. Jacobian
190 : * blocks for different tags will be extracted from Assembly. A local Jacobian will be zeroed. It
191 : * should be called right before the local element matrix is computed.
192 : */
193 : void prepareMatrixTagLower(Assembly & assembly,
194 : unsigned int ivar,
195 : unsigned int jvar,
196 : Moose::ConstraintJacobianType type);
197 :
198 : /**
199 : * Local residual blocks will be appended by adding the current local kernel residual.
200 : * It should be called after the local element vector has been computed.
201 : */
202 : void accumulateTaggedLocalResidual();
203 :
204 : /**
205 : * Local residual blocks will assigned as the current local kernel residual.
206 : * It should be called after the local element vector has been computed.
207 : */
208 : void assignTaggedLocalResidual();
209 :
210 : /**
211 : * Local Jacobian blocks will be appended by adding the current local kernel Jacobian.
212 : * It should be called after the local element matrix has been computed.
213 : */
214 : void accumulateTaggedLocalMatrix();
215 :
216 : void accumulateTaggedLocalMatrix(Assembly & assembly,
217 : unsigned int ivar,
218 : unsigned int jvar,
219 : const DenseMatrix<Number> & k);
220 :
221 : void accumulateTaggedLocalMatrix(Assembly & assembly,
222 : unsigned int ivar,
223 : unsigned int jvar,
224 : Moose::DGJacobianType type,
225 : const DenseMatrix<Number> & k);
226 :
227 : /**
228 : * Nonlocal Jacobian blocks will be appended by adding the current nonlocal kernel Jacobian.
229 : * It should be called after the nonlocal element matrix has been computed.
230 : */
231 : void accumulateTaggedNonlocalMatrix();
232 :
233 : /**
234 : * Local Jacobian blocks will assigned as the current local kernel Jacobian.
235 : * It should be called after the local element matrix has been computed.
236 : */
237 : void assignTaggedLocalMatrix();
238 :
239 : /**
240 : * Add the provided incoming residuals corresponding to the provided dof indices
241 : */
242 : template <typename Residuals, typename Indices>
243 : void addResiduals(Assembly & assembly,
244 : const Residuals & residuals,
245 : const Indices & dof_indices,
246 : Real scaling_factor);
247 :
248 : /**
249 : * Add the provided incoming residuals corresponding to the provided dof indices
250 : */
251 : template <typename T, typename Indices>
252 : void addResiduals(Assembly & assembly,
253 : const DenseVector<T> & residuals,
254 : const Indices & dof_indices,
255 : Real scaling_factor);
256 :
257 : /**
258 : * Add the provided incoming residuals and derivatives for the Jacobian, corresponding to the
259 : * provided dof indices
260 : */
261 : template <typename Residuals, typename Indices>
262 : void addResidualsAndJacobian(Assembly & assembly,
263 : const Residuals & residuals,
264 : const Indices & dof_indices,
265 : Real scaling_factor);
266 :
267 : /**
268 : * Add the provided residual derivatives into the Jacobian for the provided dof indices
269 : */
270 : template <typename Residuals, typename Indices>
271 : void addJacobian(Assembly & assembly,
272 : const Residuals & residuals,
273 : const Indices & dof_indices,
274 : Real scaling_factor);
275 :
276 : /**
277 : * Add the provided incoming residuals corresponding to the provided dof indices
278 : */
279 : void addResiduals(Assembly & assembly, const ADResidualsPacket & packet);
280 :
281 : /**
282 : * Add the provided incoming residuals and derivatives for the Jacobian, corresponding to the
283 : * provided dof indices
284 : */
285 : void addResidualsAndJacobian(Assembly & assembly, const ADResidualsPacket & packet);
286 :
287 : /**
288 : * Add the provided residual derivatives into the Jacobian for the provided dof indices
289 : */
290 : void addJacobian(Assembly & assembly, const ADResidualsPacket & packet);
291 :
292 : /**
293 : * Add the provided incoming residuals corresponding to the provided dof indices. This API should
294 : * only be used if the caller knows that no libMesh-level constraints (hanging nodes or periodic
295 : * boundary conditions) apply to the provided dof indices
296 : */
297 : template <typename Residuals, typename Indices>
298 : void addResidualsWithoutConstraints(Assembly & assembly,
299 : const Residuals & residuals,
300 : const Indices & dof_indices,
301 : Real scaling_factor);
302 :
303 : /**
304 : * Add the provided incoming residuals and derivatives for the Jacobian, corresponding to the
305 : * provided dof indices. This API should only be used if the caller knows that no libMesh-level
306 : * constraints (hanging nodes or periodic boundary conditions) apply to the provided dof indices
307 : */
308 : template <typename Residuals, typename Indices>
309 : void addResidualsAndJacobianWithoutConstraints(Assembly & assembly,
310 : const Residuals & residuals,
311 : const Indices & dof_indices,
312 : Real scaling_factor);
313 :
314 : /**
315 : * Add the provided residual derivatives into the Jacobian for the provided dof indices. This API
316 : * should only be used if the caller knows that no libMesh-level constraints (hanging nodes or
317 : * periodic boundary conditions) apply to the provided dof indices
318 : */
319 : template <typename Residuals, typename Indices>
320 : void addJacobianWithoutConstraints(Assembly & assembly,
321 : const Residuals & residuals,
322 : const Indices & dof_indices,
323 : Real scaling_factor);
324 :
325 : /**
326 : * Add into a single Jacobian element
327 : */
328 : void addJacobianElement(Assembly & assembly,
329 : Real value,
330 : dof_id_type row_index,
331 : dof_id_type column_index,
332 : Real scaling_factor);
333 :
334 : /**
335 : * Add a local Jacobian matrix
336 : */
337 : void addJacobian(Assembly & assembly,
338 : DenseMatrix<Real> & local_k,
339 : const std::vector<dof_id_type> & row_indices,
340 : const std::vector<dof_id_type> & column_indices,
341 : Real scaling_factor);
342 :
343 : /**
344 : * Set residual using the variables' insertion API
345 : */
346 : template <typename T>
347 : void setResidual(SystemBase & sys, const T & residual, MooseVariableFE<T> & var);
348 :
349 : /**
350 : * Set residual at a specified degree of freedom index
351 : */
352 : void setResidual(SystemBase & sys, Real residual, dof_id_type dof_index);
353 :
354 : /**
355 : * Set residuals using the provided functor
356 : */
357 : template <typename SetResidualFunctor>
358 : void setResidual(SystemBase & sys, SetResidualFunctor set_residual_functor);
359 :
360 : /// SubProblem that contains tag info
361 : SubProblem & _subproblem;
362 :
363 : /// Holds local residual entries as they are accumulated by this Kernel
364 : DenseVector<Number> _local_re;
365 :
366 : /// Holds local Jacobian entries as they are accumulated by this Kernel
367 : DenseMatrix<Number> _local_ke;
368 :
369 : /// Holds nonlocal Jacobian entries as they are accumulated by this Kernel
370 : DenseMatrix<Number> _nonlocal_ke;
371 :
372 : private:
373 : /**
374 : * Prepare data for computing element residual according to the specified tags
375 : * Residual blocks for different tags will be extracted from Assembly.
376 : * A local residual will be zeroed. It should be called
377 : * right before the local element vector is computed.
378 : */
379 : void prepareVectorTagInternal(Assembly & assembly,
380 : unsigned int ivar,
381 : const std::set<TagID> & vector_tags,
382 : const std::set<TagID> & absolute_value_vector_tags);
383 :
384 : /// The vector tag ids this Kernel will contribute to
385 : std::set<TagID> _vector_tags;
386 :
387 : /// The absolute value residual tag ids
388 : std::set<TagID> _abs_vector_tags;
389 :
390 : /// The matrices this Kernel will contribute to
391 : std::set<TagID> _matrix_tags;
392 :
393 : /// A set to hold vector tags excluding the reference residual tag. If there is no reference
394 : /// residual problem, this container is the same as \p _vector_tags;
395 : std::set<TagID> _non_ref_vector_tags;
396 :
397 : /// A set to hold absolute value vector tags excluding the reference residual tag. If there is no
398 : /// reference residual problem, this container is the same as \p _abs_vector_tags;
399 : std::set<TagID> _non_ref_abs_vector_tags;
400 :
401 : /// A set of either size 1 or 0. If we have a reference residual problem and \p _vector_tags holds
402 : /// the reference vector tag, then this set holds the reference vector tags, otherwise it holds
403 : /// nothing
404 : std::set<TagID> _ref_vector_tags;
405 :
406 : /// A set of either size 1 or 0. If we have a reference residual problem and \p _abs_vector_tags
407 : /// holds the reference vector tag, then this set holds the reference vector tags, otherwise it
408 : /// holds nothing
409 : std::set<TagID> _ref_abs_vector_tags;
410 :
411 : /// Moose objct this tag works on
412 : const MooseObject & _moose_object;
413 :
414 : /// Parameters from moose object
415 : const InputParameters & _tag_params;
416 :
417 : /// Residual blocks Vectors For each Tag
418 : std::vector<DenseVector<Number> *> _re_blocks;
419 :
420 : /// Residual blocks for absolute value residual tags
421 : std::vector<DenseVector<Number> *> _absre_blocks;
422 :
423 : /// Kernel blocks Vectors For each Tag
424 : std::vector<DenseMatrix<Number> *> _ke_blocks;
425 :
426 : /// A container to hold absolute values of residuals passed into \p addResiduals. We maintain
427 : /// this data member to avoid constant dynamic heap allocations
428 : std::vector<Real> _absolute_residuals;
429 :
430 : friend class NonlinearSystemBase;
431 : };
432 :
433 : #define usingTaggingInterfaceMembers \
434 : using TaggingInterface::_subproblem; \
435 : using TaggingInterface::accumulateTaggedLocalResidual; \
436 : using TaggingInterface::accumulateTaggedLocalMatrix; \
437 : using TaggingInterface::prepareVectorTag; \
438 : using TaggingInterface::prepareMatrixTag; \
439 : using TaggingInterface::prepareVectorTagNeighbor; \
440 : using TaggingInterface::_local_re; \
441 : using TaggingInterface::prepareVectorTagLower; \
442 : using TaggingInterface::prepareMatrixTagNeighbor; \
443 : using TaggingInterface::prepareMatrixTagLower; \
444 : using TaggingInterface::_local_ke
445 :
446 : template <typename Residuals, typename Indices>
447 : void
448 42083804 : TaggingInterface::addResiduals(Assembly & assembly,
449 : const Residuals & residuals,
450 : const Indices & dof_indices,
451 : const Real scaling_factor)
452 : {
453 42083804 : assembly.cacheResiduals(
454 42083804 : residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _vector_tags);
455 42083804 : if (!_abs_vector_tags.empty())
456 : {
457 16380 : _absolute_residuals.resize(residuals.size());
458 49140 : for (const auto i : index_range(residuals))
459 32760 : _absolute_residuals[i] = std::abs(MetaPhysicL::raw_value(residuals[i]));
460 :
461 16380 : assembly.cacheResiduals(_absolute_residuals,
462 : dof_indices,
463 : scaling_factor,
464 16380 : Assembly::LocalDataKey{},
465 16380 : _abs_vector_tags);
466 : }
467 42083804 : }
468 :
469 : template <typename T, typename Indices>
470 : void
471 11605612 : TaggingInterface::addResiduals(Assembly & assembly,
472 : const DenseVector<T> & residuals,
473 : const Indices & dof_indices,
474 : const Real scaling_factor)
475 : {
476 11605612 : addResiduals(assembly, residuals.get_values(), dof_indices, scaling_factor);
477 11605612 : }
478 :
479 : template <typename Residuals, typename Indices>
480 : void
481 174821 : TaggingInterface::addResidualsWithoutConstraints(Assembly & assembly,
482 : const Residuals & residuals,
483 : const Indices & dof_indices,
484 : const Real scaling_factor)
485 : {
486 174821 : assembly.cacheResidualsWithoutConstraints(
487 174821 : residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _vector_tags);
488 174821 : if (!_abs_vector_tags.empty())
489 : {
490 216 : _absolute_residuals.resize(residuals.size());
491 936 : for (const auto i : index_range(residuals))
492 720 : _absolute_residuals[i] = std::abs(MetaPhysicL::raw_value(residuals[i]));
493 :
494 216 : assembly.cacheResidualsWithoutConstraints(_absolute_residuals,
495 : dof_indices,
496 : scaling_factor,
497 216 : Assembly::LocalDataKey{},
498 216 : _abs_vector_tags);
499 : }
500 174821 : }
501 :
502 : template <typename Residuals, typename Indices>
503 : void
504 15023201 : TaggingInterface::addResidualsAndJacobian(Assembly & assembly,
505 : const Residuals & residuals,
506 : const Indices & dof_indices,
507 : Real scaling_factor)
508 : {
509 15023201 : addResiduals(assembly, residuals, dof_indices, scaling_factor);
510 15023201 : addJacobian(assembly, residuals, dof_indices, scaling_factor);
511 15023201 : }
512 :
513 : template <typename Residuals, typename Indices>
514 : void
515 25369675 : TaggingInterface::addJacobian(Assembly & assembly,
516 : const Residuals & residuals,
517 : const Indices & dof_indices,
518 : Real scaling_factor)
519 : {
520 25369675 : assembly.cacheJacobian(
521 25369675 : residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
522 25369675 : }
523 :
524 : template <typename Residuals, typename Indices>
525 : void
526 174821 : TaggingInterface::addResidualsAndJacobianWithoutConstraints(Assembly & assembly,
527 : const Residuals & residuals,
528 : const Indices & dof_indices,
529 : Real scaling_factor)
530 : {
531 174821 : addResidualsWithoutConstraints(assembly, residuals, dof_indices, scaling_factor);
532 174821 : addJacobianWithoutConstraints(assembly, residuals, dof_indices, scaling_factor);
533 174821 : }
534 :
535 : template <typename Residuals, typename Indices>
536 : void
537 174821 : TaggingInterface::addJacobianWithoutConstraints(Assembly & assembly,
538 : const Residuals & residuals,
539 : const Indices & dof_indices,
540 : Real scaling_factor)
541 : {
542 174821 : assembly.cacheJacobianWithoutConstraints(
543 174821 : residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
544 174821 : }
545 :
546 : inline void
547 9973523 : TaggingInterface::addJacobianElement(Assembly & assembly,
548 : const Real value,
549 : const dof_id_type row_index,
550 : const dof_id_type column_index,
551 : const Real scaling_factor)
552 : {
553 9973523 : assembly.cacheJacobian(
554 9973523 : row_index, column_index, value * scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
555 9973523 : }
556 :
557 : inline void
558 6719819 : TaggingInterface::addJacobian(Assembly & assembly,
559 : DenseMatrix<Real> & local_k,
560 : const std::vector<dof_id_type> & row_indices,
561 : const std::vector<dof_id_type> & column_indices,
562 : const Real scaling_factor)
563 : {
564 13439638 : for (const auto matrix_tag : _matrix_tags)
565 6719819 : assembly.cacheJacobianBlock(
566 13439638 : local_k, row_indices, column_indices, scaling_factor, Assembly::LocalDataKey{}, matrix_tag);
567 6719819 : }
568 :
569 : template <typename T>
570 : void
571 61094840 : TaggingInterface::setResidual(SystemBase & sys, const T & residual, MooseVariableFE<T> & var)
572 : {
573 123030931 : for (const auto tag_id : _vector_tags)
574 61936091 : if (sys.hasVector(tag_id))
575 61142698 : var.insertNodalValue(sys.getVector(tag_id), residual);
576 61094840 : }
577 :
578 : inline void
579 918313 : TaggingInterface::setResidual(SystemBase & sys, const Real residual, const dof_id_type dof_index)
580 : {
581 1836626 : for (const auto tag_id : _vector_tags)
582 918313 : if (sys.hasVector(tag_id))
583 918313 : sys.getVector(tag_id).set(dof_index, residual);
584 918313 : }
585 :
586 : template <typename SetResidualFunctor>
587 : void
588 : TaggingInterface::setResidual(SystemBase & sys, const SetResidualFunctor set_residual_functor)
589 : {
590 : for (const auto tag_id : _vector_tags)
591 : if (sys.hasVector(tag_id))
592 : set_residual_functor(sys.getVector(tag_id));
593 : }
594 :
595 : inline void
596 6835416 : TaggingInterface::addResiduals(Assembly & assembly, const ADResidualsPacket & packet)
597 : {
598 6835416 : addResiduals(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor);
599 6835416 : }
600 :
601 : inline void
602 0 : TaggingInterface::addResidualsAndJacobian(Assembly & assembly, const ADResidualsPacket & packet)
603 : {
604 0 : addResidualsAndJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor);
605 0 : }
606 :
607 : inline void
608 3361708 : TaggingInterface::addJacobian(Assembly & assembly, const ADResidualsPacket & packet)
609 : {
610 3361708 : addJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor);
611 3361708 : }
|
