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 "Problem.h"
13 : #include "DiracKernelInfo.h"
14 : #include "GeometricSearchData.h"
15 : #include "MooseTypes.h"
16 : #include "VectorTag.h"
17 : #include "MooseError.h"
18 : #include "FunctorMaterialProperty.h"
19 : #include "RawValueFunctor.h"
20 : #include "ADWrapperFunctor.h"
21 :
22 : #include "libmesh/coupling_matrix.h"
23 : #include "libmesh/parameters.h"
24 :
25 : #include <memory>
26 : #include <unordered_map>
27 : #include <map>
28 : #include <vector>
29 :
30 : namespace libMesh
31 : {
32 : template <typename>
33 : class VectorValue;
34 : typedef VectorValue<Real> RealVectorValue;
35 : class GhostingFunctor;
36 : }
37 :
38 : class MooseMesh;
39 : class Factory;
40 : class Assembly;
41 : class MooseVariableFieldBase;
42 : class MooseVariableScalar;
43 : template <typename>
44 : class MooseVariableFE;
45 : typedef MooseVariableFE<Real> MooseVariable;
46 : typedef MooseVariableFE<RealVectorValue> VectorMooseVariable;
47 : typedef MooseVariableFE<RealEigenVector> ArrayMooseVariable;
48 : class RestartableDataValue;
49 : class SystemBase;
50 : class LineSearch;
51 : class FaceInfo;
52 : class MooseObjectName;
53 : namespace Moose
54 : {
55 : class FunctorEnvelopeBase;
56 : }
57 :
58 : // libMesh forward declarations
59 : namespace libMesh
60 : {
61 : class EquationSystems;
62 : class DofMap;
63 : class CouplingMatrix;
64 : template <typename T>
65 : class SparseMatrix;
66 : template <typename T>
67 : class NumericVector;
68 : class System;
69 : } // namespace libMesh
70 :
71 : using libMesh::CouplingMatrix;
72 : using libMesh::EquationSystems;
73 :
74 : /**
75 : * Generic class for solving transient nonlinear problems
76 : *
77 : */
78 : class SubProblem : public Problem
79 : {
80 : public:
81 : static InputParameters validParams();
82 :
83 : SubProblem(const InputParameters & parameters);
84 : virtual ~SubProblem();
85 :
86 : virtual libMesh::EquationSystems & es() = 0;
87 : virtual MooseMesh & mesh() = 0;
88 : virtual const MooseMesh & mesh() const = 0;
89 : virtual const MooseMesh & mesh(bool use_displaced) const = 0;
90 :
91 : /**
92 : * @returns whether there will be nonlocal coupling at any point in the simulation, e.g. whether
93 : * there are any active \emph or inactive nonlocal kernels or boundary conditions
94 : */
95 : virtual bool checkNonlocalCouplingRequirement() const = 0;
96 :
97 : /**
98 : * @return whether the given solver system \p sys_num is converged
99 : */
100 0 : virtual bool solverSystemConverged(const unsigned int sys_num) { return converged(sys_num); }
101 :
102 : /**
103 : * @return whether the given nonlinear system \p nl_sys_num is converged.
104 : */
105 : virtual bool nlConverged(const unsigned int nl_sys_num);
106 :
107 : /**
108 : * Eventually we want to convert this virtual over to taking a solver system number argument.
109 : * We will have to first convert apps to use solverSystemConverged, and then once that is done, we
110 : * can change this signature. Then we can go through the apps again and convert back to this
111 : * changed API
112 : */
113 336132 : virtual bool converged(const unsigned int sys_num) { return solverSystemConverged(sys_num); }
114 :
115 : /**
116 : * @return the nonlinear system number corresponding to the provided \p nl_sys_name
117 : */
118 : virtual unsigned int nlSysNum(const NonlinearSystemName & nl_sys_name) const = 0;
119 :
120 : /**
121 : * @return the linear system number corresponding to the provided \p linear_sys_name
122 : */
123 : virtual unsigned int linearSysNum(const LinearSystemName & linear_sys_name) const = 0;
124 :
125 : /**
126 : * @return the solver system number corresponding to the provided \p solver_sys_name
127 : */
128 : virtual unsigned int solverSysNum(const SolverSystemName & solver_sys_name) const = 0;
129 :
130 : virtual void onTimestepBegin() = 0;
131 : virtual void onTimestepEnd() = 0;
132 :
133 : virtual bool isTransient() const = 0;
134 :
135 : /// marks this problem as including/needing finite volume functionality.
136 : virtual void needFV() = 0;
137 :
138 : /// returns true if this problem includes/needs finite volume functionality.
139 : virtual bool haveFV() const = 0;
140 :
141 : /**
142 : * Whether or not the user has requested default ghosting ot be on.
143 : */
144 129778 : bool defaultGhosting() { return _default_ghosting; }
145 :
146 : /**
147 : * Create a Tag. Tags can be associated with Vectors and Matrices and allow objects
148 : * (such as Kernels) to arbitrarily contribute values to any set of vectors/matrics
149 : *
150 : * Note: If the tag is already present then this will simply return the TagID of that Tag, but the
151 : * type must be the same.
152 : *
153 : * @param tag_name The name of the tag to create, the TagID will get automatically generated
154 : * @param type The type of the tag
155 : */
156 : virtual TagID addVectorTag(const TagName & tag_name,
157 : const Moose::VectorTagType type = Moose::VECTOR_TAG_RESIDUAL);
158 :
159 : /**
160 : * Adds a vector tag to the list of vectors that will not be zeroed
161 : * when other tagged vectors are
162 : * @param tag the TagID of the vector that will be manually managed
163 : */
164 : void addNotZeroedVectorTag(const TagID tag);
165 :
166 : /**
167 : * Checks if a vector tag is in the list of vectors that will not be zeroed
168 : * when other tagged vectors are
169 : * @param tag the TagID of the vector that is currently being checked
170 : * @returns false if the tag is not within the set of vectors that are
171 : * intended to not be zero or if the set is empty. returns true otherwise
172 : */
173 : bool vectorTagNotZeroed(const TagID tag) const;
174 :
175 : /**
176 : * Get a VectorTag from a TagID.
177 : */
178 : virtual const VectorTag & getVectorTag(const TagID tag_id) const;
179 : std::vector<VectorTag> getVectorTags(const std::set<TagID> & tag_ids) const;
180 :
181 : /**
182 : * Get a TagID from a TagName.
183 : */
184 : virtual TagID getVectorTagID(const TagName & tag_name) const;
185 :
186 : /**
187 : * Retrieve the name associated with a TagID
188 : */
189 : virtual TagName vectorTagName(const TagID tag) const;
190 :
191 : /**
192 : * Return all vector tags, where a tag is represented by a map from name to ID. Can optionally be
193 : * limited to a vector tag type.
194 : */
195 : virtual const std::vector<VectorTag> &
196 : getVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const;
197 :
198 : /**
199 : * Check to see if a particular Tag exists
200 : */
201 2922724519 : virtual bool vectorTagExists(const TagID tag_id) const { return tag_id < _vector_tags.size(); }
202 :
203 : /**
204 : * Check to see if a particular Tag exists by using Tag name
205 : */
206 : virtual bool vectorTagExists(const TagName & tag_name) const;
207 :
208 : /**
209 : * The total number of tags, which can be limited to the tag type
210 : */
211 : virtual unsigned int numVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const;
212 :
213 : virtual Moose::VectorTagType vectorTagType(const TagID tag_id) const;
214 :
215 : /**
216 : * Create a Tag. Tags can be associated with Vectors and Matrices and allow objects
217 : * (such as Kernels) to arbitrarily contribute values to any set of vectors/matrics
218 : *
219 : * Note: If the tag is already present then this will simply return the TagID of that Tag
220 : *
221 : * @param tag_name The name of the tag to create, the TagID will get automatically generated
222 : */
223 : virtual TagID addMatrixTag(TagName tag_name);
224 :
225 : /**
226 : * Get a TagID from a TagName.
227 : */
228 : virtual TagID getMatrixTagID(const TagName & tag_name) const;
229 :
230 : /**
231 : * Retrieve the name associated with a TagID
232 : */
233 : virtual TagName matrixTagName(TagID tag);
234 :
235 : /**
236 : * Check to see if a particular Tag exists
237 : */
238 : virtual bool matrixTagExists(const TagName & tag_name) const;
239 :
240 : /**
241 : * Check to see if a particular Tag exists
242 : */
243 : virtual bool matrixTagExists(TagID tag_id) const;
244 :
245 : /**
246 : * The total number of tags
247 : */
248 8891055 : virtual unsigned int numMatrixTags() const { return _matrix_tag_name_to_tag_id.size(); }
249 :
250 : /**
251 : * Return all matrix tags in the system, where a tag is represented by a map from name to ID
252 : */
253 511308 : virtual std::map<TagName, TagID> & getMatrixTags() { return _matrix_tag_name_to_tag_id; }
254 :
255 : /// Whether or not this problem has the variable
256 : virtual bool hasVariable(const std::string & var_name) const = 0;
257 :
258 : /// Whether or not this problem has this linear variable
259 : virtual bool hasLinearVariable(const std::string & var_name) const;
260 :
261 : /// Whether or not this problem has this auxiliary variable
262 : virtual bool hasAuxiliaryVariable(const std::string & var_name) const;
263 :
264 : /**
265 : * Returns the variable reference for requested variable which must
266 : * be of the expected_var_type (Nonlinear vs. Auxiliary) and
267 : * expected_var_field_type (standard, scalar, vector). The default
268 : * values of VAR_ANY and VAR_FIELD_ANY should be used when "any"
269 : * type of variable is acceptable. Throws an error if the variable
270 : * in question is not in the expected System or of the expected
271 : * type.
272 : */
273 : virtual const MooseVariableFieldBase & getVariable(
274 : const THREAD_ID tid,
275 : const std::string & var_name,
276 : Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
277 : Moose::VarFieldType expected_var_field_type = Moose::VarFieldType::VAR_FIELD_ANY) const = 0;
278 : virtual MooseVariableFieldBase &
279 4323722 : getVariable(const THREAD_ID tid,
280 : const std::string & var_name,
281 : Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
282 : Moose::VarFieldType expected_var_field_type = Moose::VarFieldType::VAR_FIELD_ANY)
283 : {
284 4323722 : return const_cast<MooseVariableFieldBase &>(const_cast<const SubProblem *>(this)->getVariable(
285 4323716 : tid, var_name, expected_var_type, expected_var_field_type));
286 : }
287 :
288 : /// Returns the variable reference for requested MooseVariable which may be in any system
289 : virtual MooseVariable & getStandardVariable(const THREAD_ID tid,
290 : const std::string & var_name) = 0;
291 :
292 : /// Returns the variable reference for requested MooseVariableField which may be in any system
293 : virtual MooseVariableFieldBase & getActualFieldVariable(const THREAD_ID tid,
294 : const std::string & var_name) = 0;
295 :
296 : /// Returns the variable reference for requested VectorMooseVariable which may be in any system
297 : virtual VectorMooseVariable & getVectorVariable(const THREAD_ID tid,
298 : const std::string & var_name) = 0;
299 :
300 : /// Returns the variable reference for requested ArrayMooseVariable which may be in any system
301 : virtual ArrayMooseVariable & getArrayVariable(const THREAD_ID tid,
302 : const std::string & var_name) = 0;
303 :
304 : /// Returns a Boolean indicating whether any system contains a variable with the name provided
305 : virtual bool hasScalarVariable(const std::string & var_name) const = 0;
306 :
307 : /// Returns the scalar variable reference from whichever system contains it
308 : virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid,
309 : const std::string & var_name) = 0;
310 :
311 : /// Returns the equation system containing the variable provided
312 : virtual libMesh::System & getSystem(const std::string & var_name) = 0;
313 :
314 : /**
315 : * Set the MOOSE variables to be reinited on each element.
316 : * @param moose_vars A set of variables that need to be reinited each time reinit() is called.
317 : *
318 : * @param tid The thread id
319 : */
320 : virtual void
321 : setActiveElementalMooseVariables(const std::set<MooseVariableFieldBase *> & moose_vars,
322 : const THREAD_ID tid);
323 :
324 : /**
325 : * Get the MOOSE variables to be reinited on each element.
326 : *
327 : * @param tid The thread id
328 : */
329 : virtual const std::set<MooseVariableFieldBase *> &
330 : getActiveElementalMooseVariables(const THREAD_ID tid) const;
331 :
332 : /**
333 : * Whether or not a list of active elemental moose variables has been set.
334 : *
335 : * @return True if there has been a list of active elemental moose variables set, False otherwise
336 : */
337 : virtual bool hasActiveElementalMooseVariables(const THREAD_ID tid) const;
338 :
339 : /**
340 : * Clear the active elemental MooseVariableFieldBase. If there are no active variables then they
341 : * will all be reinited. Call this after finishing the computation that was using a restricted set
342 : * of MooseVariableFieldBase
343 : *
344 : * @param tid The thread id
345 : */
346 : virtual void clearActiveElementalMooseVariables(const THREAD_ID tid);
347 :
348 : virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) = 0;
349 : virtual const Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) const = 0;
350 :
351 : /**
352 : * Return the nonlinear system object as a base class reference given the system number
353 : */
354 : virtual const SystemBase & systemBaseNonlinear(const unsigned int sys_num) const = 0;
355 : virtual SystemBase & systemBaseNonlinear(const unsigned int sys_num) = 0;
356 : /**
357 : * Return the linear system object as a base class reference given the system number
358 : */
359 : virtual const SystemBase & systemBaseLinear(const unsigned int sys_num) const = 0;
360 : virtual SystemBase & systemBaseLinear(const unsigned int sys_num) = 0;
361 : /**
362 : * Return the solver system object as a base class reference given the system number
363 : */
364 : virtual const SystemBase & systemBaseSolver(const unsigned int sys_num) const = 0;
365 : virtual SystemBase & systemBaseSolver(const unsigned int sys_num) = 0;
366 : /**
367 : * Return the auxiliary system object as a base class reference
368 : */
369 : virtual const SystemBase & systemBaseAuxiliary() const = 0;
370 : virtual SystemBase & systemBaseAuxiliary() = 0;
371 :
372 : virtual void prepareShapes(unsigned int var, const THREAD_ID tid) = 0;
373 : virtual void prepareFaceShapes(unsigned int var, const THREAD_ID tid) = 0;
374 : virtual void prepareNeighborShapes(unsigned int var, const THREAD_ID tid) = 0;
375 : Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) const;
376 :
377 : /**
378 : * Returns the desired radial direction for RZ coordinate transformation
379 : * @return The coordinate direction for the radial direction
380 : */
381 : unsigned int getAxisymmetricRadialCoord() const;
382 :
383 : virtual DiracKernelInfo & diracKernelInfo();
384 : virtual Real finalNonlinearResidual(const unsigned int nl_sys_num) const;
385 : virtual unsigned int nNonlinearIterations(const unsigned int nl_sys_num) const;
386 : virtual unsigned int nLinearIterations(const unsigned int nl_sys_num) const;
387 :
388 : virtual void addResidual(const THREAD_ID tid) = 0;
389 : virtual void addResidualNeighbor(const THREAD_ID tid) = 0;
390 : virtual void addResidualLower(const THREAD_ID tid) = 0;
391 :
392 : virtual void cacheResidual(const THREAD_ID tid);
393 : virtual void cacheResidualNeighbor(const THREAD_ID tid);
394 : virtual void addCachedResidual(const THREAD_ID tid);
395 :
396 : virtual void setResidual(libMesh::NumericVector<libMesh::Number> & residual,
397 : const THREAD_ID tid) = 0;
398 : virtual void setResidualNeighbor(libMesh::NumericVector<libMesh::Number> & residual,
399 : const THREAD_ID tid) = 0;
400 :
401 : virtual void addJacobian(const THREAD_ID tid) = 0;
402 : virtual void addJacobianNeighbor(const THREAD_ID tid) = 0;
403 : virtual void addJacobianNeighborLowerD(const THREAD_ID tid) = 0;
404 : virtual void addJacobianLowerD(const THREAD_ID tid) = 0;
405 : virtual void addJacobianNeighbor(libMesh::SparseMatrix<libMesh::Number> & jacobian,
406 : unsigned int ivar,
407 : unsigned int jvar,
408 : const libMesh::DofMap & dof_map,
409 : std::vector<dof_id_type> & dof_indices,
410 : std::vector<dof_id_type> & neighbor_dof_indices,
411 : const std::set<TagID> & tags,
412 : const THREAD_ID tid) = 0;
413 :
414 : virtual void cacheJacobian(const THREAD_ID tid);
415 : virtual void cacheJacobianNeighbor(const THREAD_ID tid);
416 : virtual void addCachedJacobian(const THREAD_ID tid);
417 :
418 : virtual void prepare(const Elem * elem, const THREAD_ID tid) = 0;
419 : virtual void prepareFace(const Elem * elem, const THREAD_ID tid) = 0;
420 : virtual void prepare(const Elem * elem,
421 : unsigned int ivar,
422 : unsigned int jvar,
423 : const std::vector<dof_id_type> & dof_indices,
424 : const THREAD_ID tid) = 0;
425 : virtual void setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid) = 0;
426 : virtual void
427 : setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
428 : virtual void prepareAssembly(const THREAD_ID tid) = 0;
429 :
430 : virtual void reinitElem(const Elem * elem, const THREAD_ID tid) = 0;
431 : virtual void reinitElemPhys(const Elem * elem,
432 : const std::vector<Point> & phys_points_in_elem,
433 : const THREAD_ID tid) = 0;
434 : virtual void reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
435 : virtual void reinitLowerDElem(const Elem * lower_d_elem,
436 : const THREAD_ID tid,
437 : const std::vector<Point> * const pts = nullptr,
438 : const std::vector<Real> * const weights = nullptr);
439 : virtual void reinitNode(const Node * node, const THREAD_ID tid) = 0;
440 : virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid) = 0;
441 : virtual void reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) = 0;
442 : virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) = 0;
443 : virtual void reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
444 : virtual void reinitNeighborPhys(const Elem * neighbor,
445 : unsigned int neighbor_side,
446 : const std::vector<Point> & physical_points,
447 : const THREAD_ID tid) = 0;
448 : virtual void reinitNeighborPhys(const Elem * neighbor,
449 : const std::vector<Point> & physical_points,
450 : const THREAD_ID tid) = 0;
451 : virtual void
452 : reinitElemNeighborAndLowerD(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
453 : /**
454 : * fills the VariableValue arrays for scalar variables from the solution vector
455 : * @param tid The thread id
456 : * @param reinit_for_derivative_reordering A flag indicating whether we are reinitializing for the
457 : * purpose of re-ordering derivative information for ADNodalBCs
458 : */
459 : virtual void reinitScalars(const THREAD_ID tid,
460 : bool reinit_for_derivative_reordering = false) = 0;
461 : virtual void reinitOffDiagScalars(const THREAD_ID tid) = 0;
462 :
463 : /// sets the current boundary ID in assembly
464 : virtual void setCurrentBoundaryID(BoundaryID bid, const THREAD_ID tid);
465 :
466 : /**
467 : * reinitialize FE objects on a given element on a given side at a given set of reference
468 : * points and then compute variable data. Note that this method makes no assumptions about what's
469 : * been called beforehand, e.g. you don't have to call some prepare method before this one. This
470 : * is an all-in-one reinit
471 : */
472 : virtual void reinitElemFaceRef(const Elem * elem,
473 : unsigned int side,
474 : Real tolerance,
475 : const std::vector<Point> * const pts,
476 : const std::vector<Real> * const weights = nullptr,
477 : const THREAD_ID tid = 0);
478 :
479 : /**
480 : * reinitialize FE objects on a given neighbor element on a given side at a given set of reference
481 : * points and then compute variable data. Note that this method makes no assumptions about what's
482 : * been called beforehand, e.g. you don't have to call some prepare method before this one. This
483 : * is an all-in-one reinit
484 : */
485 : virtual void reinitNeighborFaceRef(const Elem * neighbor_elem,
486 : unsigned int neighbor_side,
487 : Real tolerance,
488 : const std::vector<Point> * const pts,
489 : const std::vector<Real> * const weights = nullptr,
490 : const THREAD_ID tid = 0);
491 :
492 : /**
493 : * reinitialize a neighboring lower dimensional element
494 : */
495 : void reinitNeighborLowerDElem(const Elem * elem, const THREAD_ID tid = 0);
496 :
497 : /**
498 : * Reinit a mortar element to obtain a valid JxW
499 : */
500 : void reinitMortarElem(const Elem * elem, const THREAD_ID tid = 0);
501 :
502 : /**
503 : * Returns true if the Problem has Dirac kernels it needs to compute on elem.
504 : */
505 : virtual bool reinitDirac(const Elem * elem, const THREAD_ID tid) = 0;
506 : /**
507 : * Fills "elems" with the elements that should be looped over for Dirac Kernels
508 : */
509 : virtual void getDiracElements(std::set<const Elem *> & elems) = 0;
510 : /**
511 : * Gets called before Dirac Kernels are asked to add the points they are supposed to be evaluated
512 : * in
513 : */
514 : virtual void clearDiracInfo() = 0;
515 :
516 : /**
517 : * update geometric search data
518 : */
519 : virtual void
520 : updateGeomSearch(GeometricSearchData::GeometricSearchType type = GeometricSearchData::ALL) = 0;
521 :
522 : /**
523 : * reinitialize this object's geometric search data, e.g. do things like clear and re-add
524 : * quadrature nodes
525 : */
526 : void reinitGeomSearch();
527 :
528 : virtual GeometricSearchData & geomSearchData() = 0;
529 :
530 : /**
531 : * Adds the given material property to a storage map based on block ids
532 : *
533 : * This is method is called from within the Material class when the property
534 : * is first registered.
535 : * @param block_id The block id for the MaterialProperty
536 : * @param name The name of the property
537 : */
538 : virtual void storeSubdomainMatPropName(SubdomainID block_id, const std::string & name);
539 :
540 : /**
541 : * Adds the given material property to a storage map based on boundary ids
542 : *
543 : * This is method is called from within the Material class when the property
544 : * is first registered.
545 : * @param boundary_id The block id for the MaterialProperty
546 : * @param name The name of the property
547 : */
548 : virtual void storeBoundaryMatPropName(BoundaryID boundary_id, const std::string & name);
549 :
550 : /**
551 : * Adds to a map based on block ids of material properties for which a zero
552 : * value can be returned. Thes properties are optional and will not trigger a
553 : * missing material property error.
554 : *
555 : * @param block_id The block id for the MaterialProperty
556 : * @param name The name of the property
557 : */
558 : virtual void storeSubdomainZeroMatProp(SubdomainID block_id, const MaterialPropertyName & name);
559 :
560 : /**
561 : * Adds to a map based on boundary ids of material properties for which a zero
562 : * value can be returned. Thes properties are optional and will not trigger a
563 : * missing material property error.
564 : *
565 : * @param boundary_id The block id for the MaterialProperty
566 : * @param name The name of the property
567 : */
568 : virtual void storeBoundaryZeroMatProp(BoundaryID boundary_id, const MaterialPropertyName & name);
569 :
570 : /**
571 : * Adds to a map based on block ids of material properties to validate
572 : *
573 : * @param block_id The block id for the MaterialProperty
574 : * @param name The name of the property
575 : */
576 : virtual void storeSubdomainDelayedCheckMatProp(const std::string & requestor,
577 : SubdomainID block_id,
578 : const std::string & name);
579 :
580 : /**
581 : * Adds to a map based on boundary ids of material properties to validate
582 : *
583 : * @param requestor The MOOSE object name requesting the material property
584 : * @param boundary_id The block id for the MaterialProperty
585 : * @param name The name of the property
586 : */
587 : virtual void storeBoundaryDelayedCheckMatProp(const std::string & requestor,
588 : BoundaryID boundary_id,
589 : const std::string & name);
590 :
591 : /**
592 : * Checks block material properties integrity
593 : *
594 : * \see FEProblemBase::checkProblemIntegrity
595 : */
596 : virtual void checkBlockMatProps();
597 :
598 : /**
599 : * Checks boundary material properties integrity
600 : *
601 : * \see FEProblemBase::checkProblemIntegrity
602 : */
603 : virtual void checkBoundaryMatProps();
604 :
605 : /**
606 : * Helper method for adding a material property name to the _material_property_requested set
607 : */
608 : virtual void markMatPropRequested(const std::string &);
609 :
610 : /**
611 : * Find out if a material property has been requested by any object
612 : */
613 : virtual bool isMatPropRequested(const std::string & prop_name) const;
614 :
615 : /**
616 : * Helper for tracking the object that is consuming a property for MaterialPropertyDebugOutput
617 : */
618 : void addConsumedPropertyName(const MooseObjectName & obj_name, const std::string & prop_name);
619 :
620 : /**
621 : * Return the map that tracks the object with consumed material properties
622 : */
623 : const std::map<MooseObjectName, std::set<std::string>> & getConsumedPropertyMap() const;
624 :
625 : /**
626 : * Will make sure that all dofs connected to elem_id are ghosted to this processor
627 : */
628 : virtual void addGhostedElem(dof_id_type elem_id) = 0;
629 :
630 : /**
631 : * Will make sure that all necessary elements from boundary_id are ghosted to this processor
632 : */
633 : virtual void addGhostedBoundary(BoundaryID boundary_id) = 0;
634 :
635 : /**
636 : * Causes the boundaries added using addGhostedBoundary to actually be ghosted.
637 : */
638 : virtual void ghostGhostedBoundaries() = 0;
639 :
640 : /**
641 : * Get a vector containing the block ids the material property is defined on.
642 : */
643 : virtual std::set<SubdomainID> getMaterialPropertyBlocks(const std::string & prop_name);
644 :
645 : /**
646 : * Get a vector of block id equivalences that the material property is defined on.
647 : */
648 : virtual std::vector<SubdomainName> getMaterialPropertyBlockNames(const std::string & prop_name);
649 :
650 : /**
651 : * Check if a material property is defined on a block.
652 : */
653 : virtual bool hasBlockMaterialProperty(SubdomainID block_id, const std::string & prop_name);
654 :
655 : /**
656 : * Get a vector containing the block ids the material property is defined on.
657 : */
658 : virtual std::set<BoundaryID> getMaterialPropertyBoundaryIDs(const std::string & prop_name);
659 :
660 : /**
661 : * Get a vector of block id equivalences that the material property is defined on.
662 : */
663 : virtual std::vector<BoundaryName> getMaterialPropertyBoundaryNames(const std::string & prop_name);
664 :
665 : /**
666 : * Check if a material property is defined on a block.
667 : */
668 : virtual bool hasBoundaryMaterialProperty(BoundaryID boundary_id, const std::string & prop_name);
669 :
670 : /**
671 : * Returns true if the problem is in the process of computing it's initial residual.
672 : * @return Whether or not the problem is currently computing the initial residual.
673 : */
674 : virtual bool computingPreSMOResidual(const unsigned int nl_sys_num) const = 0;
675 :
676 : /**
677 : * Return the list of elements that should have their DoFs ghosted to this processor.
678 : * @return The list
679 : */
680 79420 : virtual std::set<dof_id_type> & ghostedElems() { return _ghosted_elems; }
681 :
682 : std::map<std::string, std::vector<dof_id_type>> _var_dof_map;
683 :
684 : /**
685 : * @return the nonlocal coupling matrix for the i'th nonlinear system
686 : */
687 : virtual const libMesh::CouplingMatrix & nonlocalCouplingMatrix(const unsigned i) const = 0;
688 :
689 : /**
690 : * Returns true if the problem is in the process of computing the Jacobian
691 : */
692 393682837 : const bool & currentlyComputingJacobian() const { return _currently_computing_jacobian; };
693 :
694 : /**
695 : * Set whether or not the problem is in the process of computing the Jacobian
696 : */
697 3711198 : void setCurrentlyComputingJacobian(const bool currently_computing_jacobian)
698 : {
699 3711198 : _currently_computing_jacobian = currently_computing_jacobian;
700 3711198 : }
701 :
702 : /**
703 : * Returns true if the problem is in the process of computing the residual and the Jacobian
704 : */
705 : const bool & currentlyComputingResidualAndJacobian() const;
706 :
707 : /**
708 : * Set whether or not the problem is in the process of computing the Jacobian
709 : */
710 : void setCurrentlyComputingResidualAndJacobian(bool currently_computing_residual_and_jacobian);
711 :
712 : /**
713 : * Returns true if the problem is in the process of computing the nonlinear residual
714 : */
715 21734 : bool computingNonlinearResid() const { return _computing_nonlinear_residual; }
716 :
717 : /**
718 : * Set whether or not the problem is in the process of computing the nonlinear residual
719 : */
720 192292 : virtual void computingNonlinearResid(const bool computing_nonlinear_residual)
721 : {
722 192292 : _computing_nonlinear_residual = computing_nonlinear_residual;
723 192292 : }
724 :
725 : /**
726 : * Returns true if the problem is in the process of computing the residual
727 : */
728 71601 : const bool & currentlyComputingResidual() const { return _currently_computing_residual; }
729 :
730 : /**
731 : * Set whether or not the problem is in the process of computing the residual
732 : */
733 538876 : virtual void setCurrentlyComputingResidual(const bool currently_computing_residual)
734 : {
735 538876 : _currently_computing_residual = currently_computing_residual;
736 538876 : }
737 :
738 : /// Is it safe to access the tagged matrices
739 1046876715 : virtual bool safeAccessTaggedMatrices() const { return _safe_access_tagged_matrices; }
740 :
741 : /// Is it safe to access the tagged vectors
742 547996 : virtual bool safeAccessTaggedVectors() const { return _safe_access_tagged_vectors; }
743 :
744 : virtual void clearActiveFEVariableCoupleableMatrixTags(const THREAD_ID tid);
745 :
746 : virtual void clearActiveFEVariableCoupleableVectorTags(const THREAD_ID tid);
747 :
748 : virtual void setActiveFEVariableCoupleableVectorTags(std::set<TagID> & vtags,
749 : const THREAD_ID tid);
750 :
751 : virtual void setActiveFEVariableCoupleableMatrixTags(std::set<TagID> & mtags,
752 : const THREAD_ID tid);
753 :
754 : virtual void clearActiveScalarVariableCoupleableMatrixTags(const THREAD_ID tid);
755 :
756 : virtual void clearActiveScalarVariableCoupleableVectorTags(const THREAD_ID tid);
757 :
758 : virtual void setActiveScalarVariableCoupleableVectorTags(std::set<TagID> & vtags,
759 : const THREAD_ID tid);
760 :
761 : virtual void setActiveScalarVariableCoupleableMatrixTags(std::set<TagID> & mtags,
762 : const THREAD_ID tid);
763 :
764 : const std::set<TagID> & getActiveScalarVariableCoupleableVectorTags(const THREAD_ID tid) const;
765 :
766 : const std::set<TagID> & getActiveScalarVariableCoupleableMatrixTags(const THREAD_ID tid) const;
767 :
768 : const std::set<TagID> & getActiveFEVariableCoupleableVectorTags(const THREAD_ID tid) const;
769 :
770 : const std::set<TagID> & getActiveFEVariableCoupleableMatrixTags(const THREAD_ID tid) const;
771 :
772 : /**
773 : * Method for setting whether we have any ad objects
774 : */
775 511 : virtual void haveADObjects(bool have_ad_objects) { _have_ad_objects = have_ad_objects; }
776 : /**
777 : * Method for reading wehther we have any ad objects
778 : */
779 394863369 : bool haveADObjects() const { return _have_ad_objects; }
780 :
781 : virtual LineSearch * getLineSearch() = 0;
782 :
783 : /**
784 : * The coupling matrix defining what blocks exist in the preconditioning matrix
785 : */
786 : virtual const libMesh::CouplingMatrix * couplingMatrix(const unsigned int nl_sys_num) const = 0;
787 :
788 : private:
789 : /**
790 : * Creates (n_sys - 1) clones of the provided algebraic ghosting functor (corresponding to the
791 : * nonlinear system algebraic ghosting functor), initializes the clone with the appropriate
792 : * DofMap, and then adds the clone to said DofMap
793 : * @param algebraic_gf the (nonlinear system's) algebraic ghosting functor to clone
794 : * @param to_mesh whether the clone should be added to the corresponding DofMap's underlying
795 : * MeshBase (the underlying MeshBase will be the same for every system held by this object's
796 : * EquationSystems object)
797 : */
798 : void cloneAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf, bool to_mesh = true);
799 :
800 : /**
801 : * Creates (n_sys - 1) clones of the provided coupling ghosting functor (corresponding to the
802 : * nonlinear system coupling ghosting functor), initializes the clone with the appropriate
803 : * DofMap, and then adds the clone to said DofMap
804 : * @param coupling_gf the (nonlinear system's) coupling ghosting functor to clone
805 : * @param to_mesh whether the clone should be added to the corresponding DofMap's underlying
806 : * MeshBase (the underlying MeshBase will be the same for every system held by this object's
807 : * EquationSystems object)
808 : */
809 : void cloneCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf, bool to_mesh = true);
810 :
811 : public:
812 : /**
813 : * Add an algebraic ghosting functor to this problem's DofMaps
814 : */
815 : void addAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf, bool to_mesh = true);
816 :
817 : /**
818 : * Add a coupling functor to this problem's DofMaps
819 : */
820 : void addCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf, bool to_mesh = true);
821 :
822 : /**
823 : * Remove an algebraic ghosting functor from this problem's DofMaps
824 : */
825 : void removeAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf);
826 :
827 : /**
828 : * Remove a coupling ghosting functor from this problem's DofMaps
829 : */
830 : void removeCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf);
831 :
832 : /**
833 : * Automatic scaling setter
834 : * @param automatic_scaling A boolean representing whether we are performing automatic scaling
835 : */
836 : virtual void automaticScaling(bool automatic_scaling);
837 :
838 : /**
839 : * Automatic scaling getter
840 : * @return A boolean representing whether we are performing automatic scaling
841 : */
842 : bool automaticScaling() const;
843 :
844 : /**
845 : * Tells this problem that the assembly associated with the given nonlinear system number involves
846 : * a scaling vector
847 : */
848 : void hasScalingVector(const unsigned int nl_sys_num);
849 :
850 : /**
851 : * Whether we have a displaced problem in our simulation
852 : */
853 : virtual bool haveDisplaced() const = 0;
854 :
855 : /**
856 : * Getter for whether we're computing the scaling jacobian
857 : */
858 : virtual bool computingScalingJacobian() const = 0;
859 :
860 : /**
861 : * Getter for whether we're computing the scaling residual
862 : */
863 : virtual bool computingScalingResidual() const = 0;
864 :
865 : /**
866 : * Clear dof indices from variables in nl and aux systems
867 : */
868 : void clearAllDofIndices();
869 :
870 : /**
871 : * @tparam T The type that the functor will return when evaluated, e.g. \p
872 : ADReal or \p Real
873 : * @param name The name of the functor to retrieve
874 : * @param tid The thread ID that we are retrieving the functor property for
875 : * @param requestor_name The name of the object that is requesting this functor property
876 : * @param requestor_is_ad Whether the requesting object is an AD object
877 : * @return a constant reference to the functor
878 : */
879 : template <typename T>
880 : const Moose::Functor<T> & getFunctor(const std::string & name,
881 : const THREAD_ID tid,
882 : const std::string & requestor_name,
883 : bool requestor_is_ad);
884 :
885 : /**
886 : * checks whether we have a functor corresponding to \p name on the thread id \p tid
887 : */
888 : bool hasFunctor(const std::string & name, const THREAD_ID tid) const;
889 :
890 : /**
891 : * checks whether we have a functor of type T corresponding to \p name on the thread id \p tid
892 : */
893 : template <typename T>
894 : bool hasFunctorWithType(const std::string & name, const THREAD_ID tid) const;
895 :
896 : /**
897 : * add a functor to the problem functor container
898 : */
899 : template <typename T>
900 : void
901 : addFunctor(const std::string & name, const Moose::FunctorBase<T> & functor, const THREAD_ID tid);
902 :
903 : /**
904 : * Add a functor that has block-wise lambda definitions, e.g. the evaluations of the functor are
905 : * based on a user-provided lambda expression.
906 : * @param name The name of the functor to add
907 : * @param my_lammy The lambda expression that will be called when the functor is evaluated
908 : * @param clearance_schedule How often to clear functor evaluations. The default value is always,
909 : * which means that the functor will be re-evaluated every time it is called. If it is something
910 : * other than always, than cached values may be returned
911 : * @param mesh The mesh on which this functor operates
912 : * @param block_ids The blocks on which the lambda expression is defined
913 : * @param tid The thread on which the functor we are adding will run
914 : * @return The added functor
915 : */
916 : template <typename T, typename PolymorphicLambda>
917 : const Moose::FunctorBase<T> &
918 : addPiecewiseByBlockLambdaFunctor(const std::string & name,
919 : PolymorphicLambda my_lammy,
920 : const std::set<ExecFlagType> & clearance_schedule,
921 : const MooseMesh & mesh,
922 : const std::set<SubdomainID> & block_ids,
923 : const THREAD_ID tid);
924 :
925 : virtual void initialSetup();
926 : virtual void timestepSetup();
927 : virtual void customSetup(const ExecFlagType & exec_type);
928 : virtual void residualSetup();
929 : virtual void jacobianSetup();
930 :
931 : /// Setter for debug functor output
932 33 : void setFunctorOutput(bool set_output) { _show_functors = set_output; }
933 : /// Setter for debug chain control data output
934 12 : void setChainControlDataOutput(bool set_output) { _show_chain_control_data = set_output; }
935 :
936 : #ifndef NDEBUG
937 : /// Whether to check residual for NaN/Inf values
938 : virtual bool checkResidualForNans() const = 0;
939 : #endif
940 :
941 : /**
942 : * @return the number of nonlinear systems in the problem
943 : */
944 : virtual std::size_t numNonlinearSystems() const = 0;
945 :
946 : /**
947 : * @return the current nonlinear system number
948 : */
949 : virtual unsigned int currentNlSysNum() const = 0;
950 :
951 : /**
952 : * @return the number of linear systems in the problem
953 : */
954 : virtual std::size_t numLinearSystems() const = 0;
955 :
956 : /**
957 : * @return the number of solver systems in the problem
958 : */
959 : virtual std::size_t numSolverSystems() const = 0;
960 :
961 : /**
962 : * @return the current linear system number
963 : */
964 : virtual unsigned int currentLinearSysNum() const = 0;
965 :
966 : /**
967 : * Register an unfulfilled functor request
968 : */
969 : template <typename T>
970 : void registerUnfilledFunctorRequest(T * functor_interface,
971 : const std::string & functor_name,
972 : const THREAD_ID tid);
973 :
974 : /**
975 : * Return the residual vector tags we are currently computing
976 : */
977 : virtual const std::vector<VectorTag> & currentResidualVectorTags() const = 0;
978 :
979 : /**
980 : * Select the vector tags which belong to a specific system
981 : * @param system Reference to the system
982 : * @param input_vector_tags A vector of vector tags
983 : * @param selected_tags A set which gets populated by the tag-ids that belong to the system
984 : */
985 : static void selectVectorTagsFromSystem(const SystemBase & system,
986 : const std::vector<VectorTag> & input_vector_tags,
987 : std::set<TagID> & selected_tags);
988 :
989 : /**
990 : * Select the matrix tags which belong to a specific system
991 : * @param system Reference to the system
992 : * @param input_matrix_tags A map of matrix tags
993 : * @param selected_tags A set which gets populated by the tag-ids that belong to the system
994 : */
995 : static void selectMatrixTagsFromSystem(const SystemBase & system,
996 : const std::map<TagName, TagID> & input_matrix_tags,
997 : std::set<TagID> & selected_tags);
998 :
999 : /**
1000 : * reinitialize the finite volume assembly data for the provided face and thread
1001 : */
1002 : void reinitFVFace(const THREAD_ID tid, const FaceInfo & fi);
1003 :
1004 : /**
1005 : * Whether the simulation has active nonlocal coupling which should be accounted for in the
1006 : * Jacobian. For this to return true, there must be at least one active nonlocal kernel or
1007 : * boundary condition
1008 : */
1009 : virtual bool hasNonlocalCoupling() const = 0;
1010 :
1011 : /**
1012 : * Prepare \p DofMap and \p Assembly classes with our p-refinement information
1013 : */
1014 : void preparePRefinement();
1015 :
1016 : /**
1017 : * @returns whether we're doing p-refinement
1018 : */
1019 : [[nodiscard]] bool doingPRefinement() const;
1020 :
1021 : /**
1022 : * Query whether p-refinement has been requested at any point during the simulation
1023 : */
1024 276243 : [[nodiscard]] bool havePRefinement() const { return _have_p_refinement; }
1025 :
1026 : /**
1027 : * Mark a variable family for either disabling or enabling p-refinement with valid parameters of a
1028 : * variable
1029 : */
1030 : void markFamilyPRefinement(const InputParameters & params);
1031 :
1032 : /**
1033 : * Set the current lower dimensional element. This can be null
1034 : */
1035 : virtual void setCurrentLowerDElem(const Elem * const lower_d_elem, const THREAD_ID tid);
1036 :
1037 : protected:
1038 : /**
1039 : * Helper function called by getVariable that handles the logic for
1040 : * checking whether Variables of the requested type are available.
1041 : */
1042 : template <typename T>
1043 : MooseVariableFieldBase & getVariableHelper(const THREAD_ID tid,
1044 : const std::string & var_name,
1045 : Moose::VarKindType expected_var_type,
1046 : Moose::VarFieldType expected_var_field_type,
1047 : const std::vector<T> & nls,
1048 : const SystemBase & aux) const;
1049 :
1050 : /**
1051 : * Verify the integrity of _vector_tags and _typed_vector_tags
1052 : */
1053 : bool verifyVectorTags() const;
1054 :
1055 : /// The currently declared tags
1056 : std::map<TagName, TagID> _matrix_tag_name_to_tag_id;
1057 :
1058 : /// Reverse map
1059 : std::map<TagID, TagName> _matrix_tag_id_to_tag_name;
1060 :
1061 : /// The Factory for building objects
1062 : Factory & _factory;
1063 :
1064 : DiracKernelInfo _dirac_kernel_info;
1065 :
1066 : /// Map of material properties (block_id -> list of properties)
1067 : std::map<SubdomainID, std::set<std::string>> _map_block_material_props;
1068 :
1069 : /// Map for boundary material properties (boundary_id -> list of properties)
1070 : std::map<BoundaryID, std::set<std::string>> _map_boundary_material_props;
1071 :
1072 : /// Set of properties returned as zero properties
1073 : std::map<SubdomainID, std::set<MaterialPropertyName>> _zero_block_material_props;
1074 : std::map<BoundaryID, std::set<MaterialPropertyName>> _zero_boundary_material_props;
1075 :
1076 : /// set containing all material property names that have been requested by getMaterialProperty*
1077 : std::set<std::string> _material_property_requested;
1078 :
1079 : ///@{
1080 : /**
1081 : * Data structures of the requested material properties. We store them in a map
1082 : * from boundary/block id to multimap. Each of the multimaps is a list of
1083 : * requestor object names to material property names.
1084 : */
1085 : std::map<SubdomainID, std::multimap<std::string, std::string>> _map_block_material_props_check;
1086 : std::map<BoundaryID, std::multimap<std::string, std::string>> _map_boundary_material_props_check;
1087 : ///@}
1088 :
1089 : /// This is the set of MooseVariableFieldBase that will actually get reinited by a call to reinit(elem)
1090 : std::vector<std::set<MooseVariableFieldBase *>> _active_elemental_moose_variables;
1091 :
1092 : /// Whether or not there is currently a list of active elemental moose variables
1093 : /* This needs to remain <unsigned int> for threading purposes */
1094 : std::vector<unsigned int> _has_active_elemental_moose_variables;
1095 :
1096 : std::vector<std::set<TagID>> _active_fe_var_coupleable_matrix_tags;
1097 :
1098 : std::vector<std::set<TagID>> _active_fe_var_coupleable_vector_tags;
1099 :
1100 : std::vector<std::set<TagID>> _active_sc_var_coupleable_matrix_tags;
1101 :
1102 : std::vector<std::set<TagID>> _active_sc_var_coupleable_vector_tags;
1103 :
1104 : /// Whether or not to use default libMesh coupling
1105 : bool _default_ghosting;
1106 :
1107 : /// Elements that should have Dofs ghosted to the local processor
1108 : std::set<dof_id_type> _ghosted_elems;
1109 :
1110 : /// Flag to determine whether the problem is currently computing Jacobian
1111 : bool _currently_computing_jacobian;
1112 :
1113 : /// Flag to determine whether the problem is currently computing the residual and Jacobian
1114 : bool _currently_computing_residual_and_jacobian;
1115 :
1116 : /// Whether the non-linear residual is being evaluated
1117 : bool _computing_nonlinear_residual;
1118 :
1119 : /// Whether the residual is being evaluated
1120 : bool _currently_computing_residual;
1121 :
1122 : /// Is it safe to retrieve data from tagged matrices
1123 : bool _safe_access_tagged_matrices;
1124 :
1125 : /// Is it safe to retrieve data from tagged vectors
1126 : bool _safe_access_tagged_vectors;
1127 :
1128 : /// AD flag indicating whether **any** AD objects have been added
1129 : bool _have_ad_objects;
1130 :
1131 : /// the list of vector tags that will not be zeroed when all other tags are
1132 : std::unordered_set<TagID> _not_zeroed_tagged_vectors;
1133 :
1134 : private:
1135 : /**
1136 : * @return whether a given variable name is in the solver systems (reflected by the first
1137 : * member of the returned pair which is a boolean) and if so, what solver system number it is
1138 : * in (the second member of the returned pair; if the variable is not in the solver systems,
1139 : * then this will be an invalid unsigned integer)
1140 : */
1141 : virtual std::pair<bool, unsigned int>
1142 : determineSolverSystem(const std::string & var_name, bool error_if_not_found = false) const = 0;
1143 :
1144 : enum class TrueFunctorIs
1145 : {
1146 : UNSET,
1147 : NONAD,
1148 : AD
1149 : };
1150 :
1151 : /// A container holding pointers to all the functors in our problem. We hold a tuple where the
1152 : /// zeroth item in the tuple is an enumerator that describes what type of functor the "true"
1153 : /// functor is (either NONAD or AD), the first item in the tuple is the non-AD version of the
1154 : /// functor, and the second item in the tuple is the AD version of the functor
1155 : std::vector<std::multimap<std::string,
1156 : std::tuple<TrueFunctorIs,
1157 : std::unique_ptr<Moose::FunctorEnvelopeBase>,
1158 : std::unique_ptr<Moose::FunctorEnvelopeBase>>>>
1159 : _functors;
1160 :
1161 : /// Container to hold PiecewiseByBlockLambdaFunctors
1162 : std::vector<std::map<std::string, std::unique_ptr<Moose::FunctorAbstract>>> _pbblf_functors;
1163 :
1164 : /// Lists all functors in the problem
1165 : void showFunctors() const;
1166 :
1167 : /// Lists all functors and all the objects that requested them
1168 : void showFunctorRequestors() const;
1169 :
1170 : /// The requestors of functors where the key is the prop name and the value is a set of names of
1171 : /// requestors
1172 : std::map<std::string, std::set<std::string>> _functor_to_requestors;
1173 :
1174 : /// A multimap (for each thread) from unfilled functor requests to whether the requests were for
1175 : /// AD functors and whether the requestor was an AD object
1176 : std::vector<std::multimap<std::string, std::pair<bool, bool>>> _functor_to_request_info;
1177 :
1178 : /// Whether to output a list of the functors used and requested (currently only at initialSetup)
1179 : bool _show_functors;
1180 :
1181 : /// Whether to output a list of all the chain control data
1182 : bool _show_chain_control_data;
1183 :
1184 : /// The declared vector tags
1185 : std::vector<VectorTag> _vector_tags;
1186 :
1187 : /**
1188 : * The vector tags associated with each VectorTagType
1189 : * This is kept separate from _vector_tags for quick access into typed vector tags in places where
1190 : * we don't want to build a new vector every call (like in residual evaluation)
1191 : */
1192 : std::vector<std::vector<VectorTag>> _typed_vector_tags;
1193 :
1194 : /// Map of vector tag TagName to TagID
1195 : std::map<TagName, TagID> _vector_tags_name_map;
1196 :
1197 : ///@{ Helper functions for checking MaterialProperties
1198 : std::string restrictionSubdomainCheckName(SubdomainID check_id);
1199 : std::string restrictionBoundaryCheckName(BoundaryID check_id);
1200 : ///@}
1201 :
1202 : // Contains properties consumed by objects, see addConsumedPropertyName
1203 : std::map<MooseObjectName, std::set<std::string>> _consumed_material_properties;
1204 :
1205 : /// A map from a root algebraic ghosting functor, e.g. the ghosting functor passed into \p
1206 : /// removeAlgebraicGhostingFunctor, to its clones in other systems, e.g. systems other than system
1207 : /// 0
1208 : std::unordered_map<libMesh::GhostingFunctor *,
1209 : std::vector<std::shared_ptr<libMesh::GhostingFunctor>>>
1210 : _root_alg_gf_to_sys_clones;
1211 :
1212 : /// A map from a root coupling ghosting functor, e.g. the ghosting functor passed into \p
1213 : /// removeCouplingGhostingFunctor, to its clones in other systems, e.g. systems other than system
1214 : /// 0
1215 : std::unordered_map<libMesh::GhostingFunctor *,
1216 : std::vector<std::shared_ptr<libMesh::GhostingFunctor>>>
1217 : _root_coupling_gf_to_sys_clones;
1218 :
1219 : /// Whether p-refinement has been requested at any point during the simulation
1220 : bool _have_p_refinement;
1221 :
1222 : /// Indicate whether a family is disabled for p-refinement
1223 : std::unordered_map<FEFamily, bool> _family_for_p_refinement;
1224 : /// The set of variable families by default disable p-refinement
1225 : static const std::unordered_set<FEFamily> _default_families_without_p_refinement;
1226 :
1227 : friend class Restartable;
1228 : };
1229 :
1230 : template <typename T>
1231 : const Moose::Functor<T> &
1232 19085 : SubProblem::getFunctor(const std::string & name,
1233 : const THREAD_ID tid,
1234 : const std::string & requestor_name,
1235 : const bool requestor_is_ad)
1236 : {
1237 : mooseAssert(tid < _functors.size(), "Too large a thread ID");
1238 :
1239 : // Log the requestor
1240 19085 : _functor_to_requestors["wraps_" + name].insert(requestor_name);
1241 :
1242 19085 : constexpr bool requested_functor_is_ad =
1243 : !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;
1244 :
1245 19085 : auto & functor_to_request_info = _functor_to_request_info[tid];
1246 :
1247 : // Get the requested functor if we already have it
1248 19085 : auto & functors = _functors[tid];
1249 19085 : if (auto find_ret = functors.find("wraps_" + name); find_ret != functors.end())
1250 : {
1251 18950 : if (functors.count("wraps_" + name) > 1)
1252 0 : mooseError("Attempted to get a functor with the name '",
1253 : name,
1254 : "' but multiple (" + std::to_string(functors.count("wraps_" + name)) +
1255 : ") functors match. Make sure that you do not have functor material "
1256 : "properties, functions, postprocessors or variables with the same names.");
1257 :
1258 18950 : auto & [true_functor_is, non_ad_functor, ad_functor] = find_ret->second;
1259 18950 : auto & functor_wrapper = requested_functor_is_ad ? *ad_functor : *non_ad_functor;
1260 :
1261 18950 : auto * const functor = dynamic_cast<Moose::Functor<T> *>(&functor_wrapper);
1262 18950 : if (!functor)
1263 0 : mooseError("A call to SubProblem::getFunctor requested a functor named '",
1264 : name,
1265 : "' that returns the type: '",
1266 : libMesh::demangle(typeid(T).name()),
1267 : "'. However, that functor already exists and returns a different type: '",
1268 0 : functor_wrapper.returnType(),
1269 : "'");
1270 :
1271 18950 : if (functor->template wrapsType<Moose::NullFunctor<T>>())
1272 : // Store for future checking when the actual functor gets added
1273 247 : functor_to_request_info.emplace(name,
1274 494 : std::make_pair(requested_functor_is_ad, requestor_is_ad));
1275 : else
1276 : {
1277 : // We already have the actual functor
1278 18703 : if (true_functor_is == SubProblem::TrueFunctorIs::UNSET)
1279 0 : mooseError("We already have the functor; it should not be unset");
1280 :
1281 : // Check for whether this is a valid request
1282 : // We allow auxiliary variables and linear variables to be retrieved as non AD
1283 15922 : if (!requested_functor_is_ad && requestor_is_ad &&
1284 7600 : true_functor_is == SubProblem::TrueFunctorIs::AD &&
1285 3 : !(hasAuxiliaryVariable(name) || hasLinearVariable(name)))
1286 3 : mooseError("The AD object '",
1287 : requestor_name,
1288 : "' is requesting the functor '",
1289 : name,
1290 : "' as a non-AD functor even though it is truly an AD functor, which is not "
1291 : "allowed, since this may unintentionally drop derivatives.");
1292 : }
1293 :
1294 18947 : return *functor;
1295 : }
1296 :
1297 : // We don't have the functor yet but we could have it in the future. We'll create null functors
1298 : // for now
1299 135 : functor_to_request_info.emplace(name, std::make_pair(requested_functor_is_ad, requestor_is_ad));
1300 : if constexpr (requested_functor_is_ad)
1301 : {
1302 : typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
1303 : typedef T ADType;
1304 :
1305 : auto emplace_ret =
1306 129 : functors.emplace("wraps_" + name,
1307 86 : std::make_tuple(SubProblem::TrueFunctorIs::UNSET,
1308 : std::make_unique<Moose::Functor<NonADType>>(
1309 : std::make_unique<Moose::NullFunctor<NonADType>>()),
1310 : std::make_unique<Moose::Functor<ADType>>(
1311 : std::make_unique<Moose::NullFunctor<ADType>>())));
1312 :
1313 43 : return static_cast<Moose::Functor<T> &>(*(requested_functor_is_ad
1314 43 : ? std::get<2>(emplace_ret->second)
1315 43 : : std::get<1>(emplace_ret->second)));
1316 : }
1317 : else
1318 : {
1319 : typedef T NonADType;
1320 : typedef typename Moose::ADType<T>::type ADType;
1321 :
1322 : auto emplace_ret =
1323 276 : functors.emplace("wraps_" + name,
1324 184 : std::make_tuple(SubProblem::TrueFunctorIs::UNSET,
1325 : std::make_unique<Moose::Functor<NonADType>>(
1326 : std::make_unique<Moose::NullFunctor<NonADType>>()),
1327 : std::make_unique<Moose::Functor<ADType>>(
1328 : std::make_unique<Moose::NullFunctor<ADType>>())));
1329 :
1330 92 : return static_cast<Moose::Functor<T> &>(*(requested_functor_is_ad
1331 : ? std::get<2>(emplace_ret->second)
1332 184 : : std::get<1>(emplace_ret->second)));
1333 : }
1334 : }
1335 :
1336 : template <typename T>
1337 : bool
1338 192 : SubProblem::hasFunctorWithType(const std::string & name, const THREAD_ID tid) const
1339 : {
1340 : mooseAssert(tid < _functors.size(), "Too large a thread ID");
1341 192 : auto & functors = _functors[tid];
1342 :
1343 192 : const auto & it = functors.find("wraps_" + name);
1344 192 : constexpr bool requested_functor_is_ad =
1345 : !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;
1346 :
1347 192 : if (it == functors.end())
1348 0 : return false;
1349 : else
1350 192 : return dynamic_cast<Moose::Functor<T> *>(
1351 384 : requested_functor_is_ad ? std::get<2>(it->second).get() : std::get<1>(it->second).get());
1352 : }
1353 :
1354 : template <typename T, typename PolymorphicLambda>
1355 : const Moose::FunctorBase<T> &
1356 7655 : SubProblem::addPiecewiseByBlockLambdaFunctor(const std::string & name,
1357 : PolymorphicLambda my_lammy,
1358 : const std::set<ExecFlagType> & clearance_schedule,
1359 : const MooseMesh & mesh,
1360 : const std::set<SubdomainID> & block_ids,
1361 : const THREAD_ID tid)
1362 : {
1363 7655 : auto & pbblf_functors = _pbblf_functors[tid];
1364 :
1365 7655 : auto [it, first_time_added] =
1366 7655 : pbblf_functors.emplace(name,
1367 : std::make_unique<PiecewiseByBlockLambdaFunctor<T>>(
1368 : name, my_lammy, clearance_schedule, mesh, block_ids));
1369 :
1370 7655 : auto * functor = dynamic_cast<PiecewiseByBlockLambdaFunctor<T> *>(it->second.get());
1371 7655 : if (!functor)
1372 : {
1373 0 : if (first_time_added)
1374 0 : mooseError("This should be impossible. If this was the first time we added the functor, then "
1375 : "the dynamic cast absolutely should have succeeded");
1376 : else
1377 0 : mooseError("Attempted to add a lambda functor with the name '",
1378 : name,
1379 : "' but another lambda functor of that name returns a different type");
1380 : }
1381 :
1382 7655 : if (first_time_added)
1383 6828 : addFunctor(name, *functor, tid);
1384 : else
1385 : // The functor already exists
1386 827 : functor->setFunctor(mesh, block_ids, my_lammy);
1387 :
1388 7652 : return *functor;
1389 : }
1390 :
1391 : template <typename T>
1392 : void
1393 307231 : SubProblem::addFunctor(const std::string & name,
1394 : const Moose::FunctorBase<T> & functor,
1395 : const THREAD_ID tid)
1396 : {
1397 307231 : constexpr bool added_functor_is_ad =
1398 : !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;
1399 :
1400 : mooseAssert(tid < _functors.size(), "Too large a thread ID");
1401 :
1402 307231 : auto & functor_to_request_info = _functor_to_request_info[tid];
1403 307231 : auto & functors = _functors[tid];
1404 307231 : auto it = functors.find("wraps_" + name);
1405 307231 : if (it != functors.end())
1406 : {
1407 : // We have this functor already. If it's a null functor, we want to replace it with the valid
1408 : // functor we have now. If it's not then we'll add a new entry into the multimap and then we'll
1409 : // error later if a user requests a functor because their request is ambiguous. This is the
1410 : // reason that the functors container is a multimap: for nice error messages
1411 : auto * const existing_wrapper_base =
1412 6494 : added_functor_is_ad ? std::get<2>(it->second).get() : std::get<1>(it->second).get();
1413 6494 : auto * const existing_wrapper = dynamic_cast<Moose::Functor<T> *>(existing_wrapper_base);
1414 6494 : if (existing_wrapper && existing_wrapper->template wrapsType<Moose::NullFunctor<T>>())
1415 : {
1416 : // Sanity check
1417 135 : auto [request_info_it, request_info_end_it] = functor_to_request_info.equal_range(name);
1418 135 : if (request_info_it == request_info_end_it)
1419 0 : mooseError("We are wrapping a NullFunctor but we don't have any unfilled functor request "
1420 : "info. This doesn't make sense.");
1421 :
1422 : // Check for valid requests
1423 517 : while (request_info_it != request_info_end_it)
1424 : {
1425 382 : auto & [requested_functor_is_ad, requestor_is_ad] = request_info_it->second;
1426 4 : if (!requested_functor_is_ad && requestor_is_ad && added_functor_is_ad)
1427 0 : mooseError("We are requesting a non-AD functor '" + name +
1428 : "' from an AD object, but the true functor is AD. This means we could be "
1429 : "dropping important derivatives. We will not allow this");
1430 : // We're going to eventually check whether we've fulfilled all functor requests and our
1431 : // check will be that the multimap is empty. This request is fulfilled, so erase it from the
1432 : // map now
1433 382 : request_info_it = functor_to_request_info.erase(request_info_it);
1434 : }
1435 :
1436 : // Ok we didn't have the functor before, so we will add it now
1437 135 : std::get<0>(it->second) =
1438 : added_functor_is_ad ? SubProblem::TrueFunctorIs::AD : SubProblem::TrueFunctorIs::NONAD;
1439 135 : existing_wrapper->assign(functor);
1440 : // Finally we create the non-AD or AD complement of the just added functor
1441 : if constexpr (added_functor_is_ad)
1442 : {
1443 : typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
1444 4 : auto * const existing_non_ad_wrapper_base = std::get<1>(it->second).get();
1445 4 : auto * const existing_non_ad_wrapper =
1446 4 : dynamic_cast<Moose::Functor<NonADType> *>(existing_non_ad_wrapper_base);
1447 : mooseAssert(existing_non_ad_wrapper->template wrapsType<Moose::NullFunctor<NonADType>>(),
1448 : "Both members of pair should have been wrapping a NullFunctor");
1449 4 : existing_non_ad_wrapper->assign(
1450 : std::make_unique<Moose::RawValueFunctor<NonADType>>(functor));
1451 : }
1452 : else
1453 : {
1454 : typedef typename Moose::ADType<T>::type ADType;
1455 131 : auto * const existing_ad_wrapper_base = std::get<2>(it->second).get();
1456 131 : auto * const existing_ad_wrapper =
1457 131 : dynamic_cast<Moose::Functor<ADType> *>(existing_ad_wrapper_base);
1458 : mooseAssert(existing_ad_wrapper->template wrapsType<Moose::NullFunctor<ADType>>(),
1459 : "Both members of pair should have been wrapping a NullFunctor");
1460 131 : existing_ad_wrapper->assign(std::make_unique<Moose::ADWrapperFunctor<ADType>>(functor));
1461 : }
1462 135 : return;
1463 : }
1464 6359 : else if (!existing_wrapper)
1465 : {
1466 : // Functor was emplaced but the cast failed. This could be a double definition with
1467 : // different types, or it could be a request with one type then a definition with another
1468 : // type. Either way it is going to error later, but it is cleaner to catch it now
1469 0 : mooseError("Functor '",
1470 : name,
1471 : "' is being added with return type '",
1472 : MooseUtils::prettyCppType<T>(),
1473 : "' but it has already been defined or requested with return type '",
1474 0 : existing_wrapper_base->returnType(),
1475 : "'.");
1476 : }
1477 : }
1478 :
1479 : // We are a new functor, create the opposite ADType one and store it with other functors
1480 : if constexpr (added_functor_is_ad)
1481 : {
1482 : typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
1483 190119 : auto new_non_ad_wrapper = std::make_unique<Moose::Functor<NonADType>>(
1484 : std::make_unique<Moose::RawValueFunctor<NonADType>>(functor));
1485 190119 : auto new_ad_wrapper = std::make_unique<Moose::Functor<T>>(functor);
1486 380238 : _functors[tid].emplace("wraps_" + name,
1487 190119 : std::make_tuple(SubProblem::TrueFunctorIs::AD,
1488 190119 : std::move(new_non_ad_wrapper),
1489 190119 : std::move(new_ad_wrapper)));
1490 190119 : }
1491 : else
1492 : {
1493 : typedef typename Moose::ADType<T>::type ADType;
1494 116977 : auto new_non_ad_wrapper = std::make_unique<Moose::Functor<T>>((functor));
1495 116977 : auto new_ad_wrapper = std::make_unique<Moose::Functor<ADType>>(
1496 : std::make_unique<Moose::ADWrapperFunctor<ADType>>(functor));
1497 233954 : _functors[tid].emplace("wraps_" + name,
1498 116977 : std::make_tuple(SubProblem::TrueFunctorIs::NONAD,
1499 116977 : std::move(new_non_ad_wrapper),
1500 116977 : std::move(new_ad_wrapper)));
1501 116977 : }
1502 : }
1503 :
1504 : inline const bool &
1505 298161 : SubProblem::currentlyComputingResidualAndJacobian() const
1506 : {
1507 298161 : return _currently_computing_residual_and_jacobian;
1508 : }
1509 :
1510 : inline void
1511 3690051 : SubProblem::setCurrentlyComputingResidualAndJacobian(
1512 : const bool currently_computing_residual_and_jacobian)
1513 : {
1514 3690051 : _currently_computing_residual_and_jacobian = currently_computing_residual_and_jacobian;
1515 3690051 : }
1516 :
1517 : namespace Moose
1518 : {
1519 : void initial_condition(libMesh::EquationSystems & es, const std::string & system_name);
1520 : } // namespace Moose
|
