SubProblem.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "Problem.h"
#include "DiracKernelInfo.h"
#include "GeometricSearchData.h"
#include "MooseTypes.h"
#include "VectorTag.h"
#include "MooseError.h"
#include "FunctorMaterialProperty.h"
#include "RawValueFunctor.h"
#include "ADWrapperFunctor.h"

#include "libmesh/coupling_matrix.h"
#include "libmesh/parameters.h"

#include <memory>
#include <unordered_map>
#include <map>
#include <vector>

namespace libMesh
{
template <typename>
class VectorValue;
typedef VectorValue<Real> RealVectorValue;
class GhostingFunctor;
}

class MooseMesh;
class Factory;
class Assembly;
class MooseVariableFieldBase;
class MooseVariableScalar;
template <typename>
class MooseVariableFE;
typedef MooseVariableFE<Real> MooseVariable;
typedef MooseVariableFE<RealVectorValue> VectorMooseVariable;
typedef MooseVariableFE<RealEigenVector> ArrayMooseVariable;
class RestartableDataValue;
class SystemBase;
class LineSearch;
class FaceInfo;
class MooseObjectName;
namespace Moose
{
class FunctorEnvelopeBase;
}

// libMesh forward declarations
namespace libMesh
{
class EquationSystems;
class DofMap;
class CouplingMatrix;
template <typename T>
class SparseMatrix;
template <typename T>
class NumericVector;
class System;
} // namespace libMesh

using libMesh::CouplingMatrix;
using libMesh::EquationSystems;

class SubProblem : public Problem
{
public:
  static InputParameters validParams();

  SubProblem(const InputParameters & parameters);
  virtual ~SubProblem();

  virtual libMesh::EquationSystems & es() = 0;
  virtual MooseMesh & mesh() = 0;
  virtual const MooseMesh & mesh() const = 0;
  virtual const MooseMesh & mesh(bool use_displaced) const = 0;

  virtual bool checkNonlocalCouplingRequirement() const = 0;

  virtual bool solverSystemConverged(const unsigned int sys_num) { return converged(sys_num); }

  virtual bool nlConverged(const unsigned int nl_sys_num);

  virtual bool converged(const unsigned int sys_num) { return solverSystemConverged(sys_num); }

  virtual unsigned int nlSysNum(const NonlinearSystemName & nl_sys_name) const = 0;

  virtual unsigned int linearSysNum(const LinearSystemName & linear_sys_name) const = 0;

  virtual unsigned int solverSysNum(const SolverSystemName & solver_sys_name) const = 0;

  virtual void onTimestepBegin() = 0;
  virtual void onTimestepEnd() = 0;

  virtual bool isTransient() const = 0;

  virtual void needFV() = 0;

  virtual bool haveFV() const = 0;

  bool defaultGhosting() { return _default_ghosting; }

  virtual TagID addVectorTag(const TagName & tag_name,
                             const Moose::VectorTagType type = Moose::VECTOR_TAG_RESIDUAL);

  void addNotZeroedVectorTag(const TagID tag);

  bool vectorTagNotZeroed(const TagID tag) const;

  virtual const VectorTag & getVectorTag(const TagID tag_id) const;
  std::vector<VectorTag> getVectorTags(const std::set<TagID> & tag_ids) const;

  virtual TagID getVectorTagID(const TagName & tag_name) const;

  virtual TagName vectorTagName(const TagID tag) const;

  virtual const std::vector<VectorTag> &
  getVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const;

  virtual bool vectorTagExists(const TagID tag_id) const { return tag_id < _vector_tags.size(); }

  virtual bool vectorTagExists(const TagName & tag_name) const;

  virtual unsigned int numVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const;

  virtual Moose::VectorTagType vectorTagType(const TagID tag_id) const;

  virtual TagID addMatrixTag(TagName tag_name);

  virtual TagID getMatrixTagID(const TagName & tag_name) const;

  virtual TagName matrixTagName(TagID tag);

  virtual bool matrixTagExists(const TagName & tag_name) const;

  virtual bool matrixTagExists(TagID tag_id) const;

  virtual unsigned int numMatrixTags() const { return _matrix_tag_name_to_tag_id.size(); }

  virtual std::map<TagName, TagID> & getMatrixTags() { return _matrix_tag_name_to_tag_id; }

  virtual bool hasVariable(const std::string & var_name) const = 0;

  virtual bool hasLinearVariable(const std::string & var_name) const;

  virtual bool hasAuxiliaryVariable(const std::string & var_name) const;

  virtual const MooseVariableFieldBase & getVariable(
      const THREAD_ID tid,
      const std::string & var_name,
      Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
      Moose::VarFieldType expected_var_field_type = Moose::VarFieldType::VAR_FIELD_ANY) const = 0;
  virtual MooseVariableFieldBase &
  getVariable(const THREAD_ID tid,
              const std::string & var_name,
              Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
              Moose::VarFieldType expected_var_field_type = Moose::VarFieldType::VAR_FIELD_ANY)
  {
    return const_cast<MooseVariableFieldBase &>(const_cast<const SubProblem *>(this)->getVariable(
        tid, var_name, expected_var_type, expected_var_field_type));
  }

  virtual MooseVariable & getStandardVariable(const THREAD_ID tid,
                                              const std::string & var_name) = 0;

  virtual MooseVariableFieldBase & getActualFieldVariable(const THREAD_ID tid,
                                                          const std::string & var_name) = 0;

  virtual VectorMooseVariable & getVectorVariable(const THREAD_ID tid,
                                                  const std::string & var_name) = 0;

  virtual ArrayMooseVariable & getArrayVariable(const THREAD_ID tid,
                                                const std::string & var_name) = 0;

  virtual bool hasScalarVariable(const std::string & var_name) const = 0;

  virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid,
                                                  const std::string & var_name) = 0;

  virtual libMesh::System & getSystem(const std::string & var_name) = 0;

  virtual void
  setActiveElementalMooseVariables(const std::set<MooseVariableFieldBase *> & moose_vars,
                                   const THREAD_ID tid);

  virtual const std::set<MooseVariableFieldBase *> &
  getActiveElementalMooseVariables(const THREAD_ID tid) const;

  virtual bool hasActiveElementalMooseVariables(const THREAD_ID tid) const;

  virtual void clearActiveElementalMooseVariables(const THREAD_ID tid);

  virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) = 0;
  virtual const Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) const = 0;

  virtual const SystemBase & systemBaseNonlinear(const unsigned int sys_num) const = 0;
  virtual SystemBase & systemBaseNonlinear(const unsigned int sys_num) = 0;
  virtual const SystemBase & systemBaseLinear(const unsigned int sys_num) const = 0;
  virtual SystemBase & systemBaseLinear(const unsigned int sys_num) = 0;
  virtual const SystemBase & systemBaseSolver(const unsigned int sys_num) const = 0;
  virtual SystemBase & systemBaseSolver(const unsigned int sys_num) = 0;
  virtual const SystemBase & systemBaseAuxiliary() const = 0;
  virtual SystemBase & systemBaseAuxiliary() = 0;

  virtual void prepareShapes(unsigned int var, const THREAD_ID tid) = 0;
  virtual void prepareFaceShapes(unsigned int var, const THREAD_ID tid) = 0;
  virtual void prepareNeighborShapes(unsigned int var, const THREAD_ID tid) = 0;
  Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) const;

  unsigned int getAxisymmetricRadialCoord() const;

  virtual DiracKernelInfo & diracKernelInfo();
  virtual Real finalNonlinearResidual(const unsigned int nl_sys_num) const;
  virtual unsigned int nNonlinearIterations(const unsigned int nl_sys_num) const;
  virtual unsigned int nLinearIterations(const unsigned int nl_sys_num) const;

  virtual void addResidual(const THREAD_ID tid) = 0;
  virtual void addResidualNeighbor(const THREAD_ID tid) = 0;
  virtual void addResidualLower(const THREAD_ID tid) = 0;

  virtual void cacheResidual(const THREAD_ID tid);
  virtual void cacheResidualNeighbor(const THREAD_ID tid);
  virtual void addCachedResidual(const THREAD_ID tid);

  virtual void setResidual(libMesh::NumericVector<libMesh::Number> & residual,
                           const THREAD_ID tid) = 0;
  virtual void setResidualNeighbor(libMesh::NumericVector<libMesh::Number> & residual,
                                   const THREAD_ID tid) = 0;

  virtual void addJacobian(const THREAD_ID tid) = 0;
  virtual void addJacobianNeighbor(const THREAD_ID tid) = 0;
  virtual void addJacobianNeighborLowerD(const THREAD_ID tid) = 0;
  virtual void addJacobianLowerD(const THREAD_ID tid) = 0;
  virtual void addJacobianNeighbor(libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                   unsigned int ivar,
                                   unsigned int jvar,
                                   const libMesh::DofMap & dof_map,
                                   std::vector<dof_id_type> & dof_indices,
                                   std::vector<dof_id_type> & neighbor_dof_indices,
                                   const std::set<TagID> & tags,
                                   const THREAD_ID tid) = 0;

  virtual void cacheJacobian(const THREAD_ID tid);
  virtual void cacheJacobianNeighbor(const THREAD_ID tid);
  virtual void addCachedJacobian(const THREAD_ID tid);

  virtual void prepare(const Elem * elem, const THREAD_ID tid) = 0;
  virtual void prepareFace(const Elem * elem, const THREAD_ID tid) = 0;
  virtual void prepare(const Elem * elem,
                       unsigned int ivar,
                       unsigned int jvar,
                       const std::vector<dof_id_type> & dof_indices,
                       const THREAD_ID tid) = 0;
  virtual void setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid) = 0;
  virtual void
  setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
  virtual void prepareAssembly(const THREAD_ID tid) = 0;

  virtual void reinitElem(const Elem * elem, const THREAD_ID tid) = 0;
  virtual void reinitElemPhys(const Elem * elem,
                              const std::vector<Point> & phys_points_in_elem,
                              const THREAD_ID tid) = 0;
  virtual void reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
  virtual void reinitLowerDElem(const Elem * lower_d_elem,
                                const THREAD_ID tid,
                                const std::vector<Point> * const pts = nullptr,
                                const std::vector<Real> * const weights = nullptr);
  virtual void reinitNode(const Node * node, const THREAD_ID tid) = 0;
  virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid) = 0;
  virtual void reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) = 0;
  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) = 0;
  virtual void reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  unsigned int neighbor_side,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) = 0;
  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) = 0;
  virtual void
  reinitElemNeighborAndLowerD(const Elem * elem, unsigned int side, const THREAD_ID tid) = 0;
  virtual void reinitScalars(const THREAD_ID tid,
                             bool reinit_for_derivative_reordering = false) = 0;
  virtual void reinitOffDiagScalars(const THREAD_ID tid) = 0;

  virtual void setCurrentBoundaryID(BoundaryID bid, const THREAD_ID tid);

  virtual void reinitElemFaceRef(const Elem * elem,
                                 unsigned int side,
                                 Real tolerance,
                                 const std::vector<Point> * const pts,
                                 const std::vector<Real> * const weights = nullptr,
                                 const THREAD_ID tid = 0);

  virtual void reinitNeighborFaceRef(const Elem * neighbor_elem,
                                     unsigned int neighbor_side,
                                     Real tolerance,
                                     const std::vector<Point> * const pts,
                                     const std::vector<Real> * const weights = nullptr,
                                     const THREAD_ID tid = 0);

  void reinitNeighborLowerDElem(const Elem * elem, const THREAD_ID tid = 0);

  void reinitMortarElem(const Elem * elem, const THREAD_ID tid = 0);

  virtual bool reinitDirac(const Elem * elem, const THREAD_ID tid) = 0;
  virtual void getDiracElements(std::set<const Elem *> & elems) = 0;
  virtual void clearDiracInfo() = 0;

  // Geom Search
  virtual void
  updateGeomSearch(GeometricSearchData::GeometricSearchType type = GeometricSearchData::ALL) = 0;

  virtual GeometricSearchData & geomSearchData() = 0;

  virtual void storeSubdomainMatPropName(SubdomainID block_id, const std::string & name);

  virtual void storeBoundaryMatPropName(BoundaryID boundary_id, const std::string & name);

  virtual void storeSubdomainZeroMatProp(SubdomainID block_id, const MaterialPropertyName & name);

  virtual void storeBoundaryZeroMatProp(BoundaryID boundary_id, const MaterialPropertyName & name);

  virtual void storeSubdomainDelayedCheckMatProp(const std::string & requestor,
                                                 SubdomainID block_id,
                                                 const std::string & name);

  virtual void storeBoundaryDelayedCheckMatProp(const std::string & requestor,
                                                BoundaryID boundary_id,
                                                const std::string & name);

  virtual void checkBlockMatProps();

  virtual void checkBoundaryMatProps();

  virtual void markMatPropRequested(const std::string &);

  virtual bool isMatPropRequested(const std::string & prop_name) const;

  void addConsumedPropertyName(const MooseObjectName & obj_name, const std::string & prop_name);

  const std::map<MooseObjectName, std::set<std::string>> & getConsumedPropertyMap() const;

  virtual void addGhostedElem(dof_id_type elem_id) = 0;

  virtual void addGhostedBoundary(BoundaryID boundary_id) = 0;

  virtual void ghostGhostedBoundaries() = 0;

  virtual std::set<SubdomainID> getMaterialPropertyBlocks(const std::string & prop_name);

  virtual std::vector<SubdomainName> getMaterialPropertyBlockNames(const std::string & prop_name);

  virtual bool hasBlockMaterialProperty(SubdomainID block_id, const std::string & prop_name);

  virtual std::set<BoundaryID> getMaterialPropertyBoundaryIDs(const std::string & prop_name);

  virtual std::vector<BoundaryName> getMaterialPropertyBoundaryNames(const std::string & prop_name);

  virtual bool hasBoundaryMaterialProperty(BoundaryID boundary_id, const std::string & prop_name);

  virtual bool computingPreSMOResidual(const unsigned int nl_sys_num) const = 0;

  virtual std::set<dof_id_type> & ghostedElems() { return _ghosted_elems; }

  std::map<std::string, std::vector<dof_id_type>> _var_dof_map;

  virtual const libMesh::CouplingMatrix & nonlocalCouplingMatrix(const unsigned i) const = 0;

  const bool & currentlyComputingJacobian() const { return _currently_computing_jacobian; };

  void setCurrentlyComputingJacobian(const bool currently_computing_jacobian)
  {
    _currently_computing_jacobian = currently_computing_jacobian;
  }

  const bool & currentlyComputingResidualAndJacobian() const;

  void setCurrentlyComputingResidualAndJacobian(bool currently_computing_residual_and_jacobian);

  bool computingNonlinearResid() const { return _computing_nonlinear_residual; }

  virtual void computingNonlinearResid(const bool computing_nonlinear_residual)
  {
    _computing_nonlinear_residual = computing_nonlinear_residual;
  }

  const bool & currentlyComputingResidual() const { return _currently_computing_residual; }

  virtual void setCurrentlyComputingResidual(const bool currently_computing_residual)
  {
    _currently_computing_residual = currently_computing_residual;
  }

  virtual bool safeAccessTaggedMatrices() const { return _safe_access_tagged_matrices; }

  virtual bool safeAccessTaggedVectors() const { return _safe_access_tagged_vectors; }

  virtual void clearActiveFEVariableCoupleableMatrixTags(const THREAD_ID tid);

  virtual void clearActiveFEVariableCoupleableVectorTags(const THREAD_ID tid);

  virtual void setActiveFEVariableCoupleableVectorTags(std::set<TagID> & vtags,
                                                       const THREAD_ID tid);

  virtual void setActiveFEVariableCoupleableMatrixTags(std::set<TagID> & mtags,
                                                       const THREAD_ID tid);

  virtual void clearActiveScalarVariableCoupleableMatrixTags(const THREAD_ID tid);

  virtual void clearActiveScalarVariableCoupleableVectorTags(const THREAD_ID tid);

  virtual void setActiveScalarVariableCoupleableVectorTags(std::set<TagID> & vtags,
                                                           const THREAD_ID tid);

  virtual void setActiveScalarVariableCoupleableMatrixTags(std::set<TagID> & mtags,
                                                           const THREAD_ID tid);

  const std::set<TagID> & getActiveScalarVariableCoupleableVectorTags(const THREAD_ID tid) const;

  const std::set<TagID> & getActiveScalarVariableCoupleableMatrixTags(const THREAD_ID tid) const;

  const std::set<TagID> & getActiveFEVariableCoupleableVectorTags(const THREAD_ID tid) const;

  const std::set<TagID> & getActiveFEVariableCoupleableMatrixTags(const THREAD_ID tid) const;

  virtual void haveADObjects(bool have_ad_objects) { _have_ad_objects = have_ad_objects; }
  bool haveADObjects() const { return _have_ad_objects; }

  virtual LineSearch * getLineSearch() = 0;

  virtual const libMesh::CouplingMatrix * couplingMatrix(const unsigned int nl_sys_num) const = 0;

private:
  void cloneAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf, bool to_mesh = true);

  void cloneCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf, bool to_mesh = true);

public:
  void addAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf, bool to_mesh = true);

  void addCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf, bool to_mesh = true);

  void removeAlgebraicGhostingFunctor(libMesh::GhostingFunctor & algebraic_gf);

  void removeCouplingGhostingFunctor(libMesh::GhostingFunctor & coupling_gf);

  virtual void automaticScaling(bool automatic_scaling);

  bool automaticScaling() const;

  void hasScalingVector(const unsigned int nl_sys_num);

  virtual bool haveDisplaced() const = 0;

  virtual bool computingScalingJacobian() const = 0;

  virtual bool computingScalingResidual() const = 0;

  void clearAllDofIndices();

  template <typename T>
  const Moose::Functor<T> & getFunctor(const std::string & name,
                                       const THREAD_ID tid,
                                       const std::string & requestor_name,
                                       bool requestor_is_ad);

  bool hasFunctor(const std::string & name, const THREAD_ID tid) const;

  template <typename T>
  bool hasFunctorWithType(const std::string & name, const THREAD_ID tid) const;

  template <typename T>
  void
  addFunctor(const std::string & name, const Moose::FunctorBase<T> & functor, const THREAD_ID tid);

  template <typename T, typename PolymorphicLambda>
  const Moose::FunctorBase<T> &
  addPiecewiseByBlockLambdaFunctor(const std::string & name,
                                   PolymorphicLambda my_lammy,
                                   const std::set<ExecFlagType> & clearance_schedule,
                                   const MooseMesh & mesh,
                                   const std::set<SubdomainID> & block_ids,
                                   const THREAD_ID tid);

  virtual void initialSetup();
  virtual void timestepSetup();
  virtual void customSetup(const ExecFlagType & exec_type);
  virtual void residualSetup();
  virtual void jacobianSetup();

  void setFunctorOutput(bool set_output) { _output_functors = set_output; }

  virtual std::size_t numNonlinearSystems() const = 0;

  virtual unsigned int currentNlSysNum() const = 0;

  virtual std::size_t numLinearSystems() const = 0;

  virtual std::size_t numSolverSystems() const = 0;

  virtual unsigned int currentLinearSysNum() const = 0;

  template <typename T>
  void registerUnfilledFunctorRequest(T * functor_interface,
                                      const std::string & functor_name,
                                      const THREAD_ID tid);

  virtual const std::vector<VectorTag> & currentResidualVectorTags() const = 0;

  static void selectVectorTagsFromSystem(const SystemBase & system,
                                         const std::vector<VectorTag> & input_vector_tags,
                                         std::set<TagID> & selected_tags);

  static void selectMatrixTagsFromSystem(const SystemBase & system,
                                         const std::map<TagName, TagID> & input_matrix_tags,
                                         std::set<TagID> & selected_tags);

  void reinitFVFace(const THREAD_ID tid, const FaceInfo & fi);

  virtual bool hasNonlocalCoupling() const = 0;

  void preparePRefinement();

  [[nodiscard]] bool doingPRefinement() const;

  [[nodiscard]] bool havePRefinement() const { return _have_p_refinement; }

  virtual void setCurrentLowerDElem(const Elem * const lower_d_elem, const THREAD_ID tid);

protected:
  template <typename T>
  MooseVariableFieldBase & getVariableHelper(const THREAD_ID tid,
                                             const std::string & var_name,
                                             Moose::VarKindType expected_var_type,
                                             Moose::VarFieldType expected_var_field_type,
                                             const std::vector<T> & nls,
                                             const SystemBase & aux) const;

  bool verifyVectorTags() const;

  void markFamilyPRefinement(const InputParameters & params);

  std::map<TagName, TagID> _matrix_tag_name_to_tag_id;

  std::map<TagID, TagName> _matrix_tag_id_to_tag_name;

  Factory & _factory;

  DiracKernelInfo _dirac_kernel_info;

  std::map<SubdomainID, std::set<std::string>> _map_block_material_props;

  std::map<BoundaryID, std::set<std::string>> _map_boundary_material_props;

  std::map<SubdomainID, std::set<MaterialPropertyName>> _zero_block_material_props;
  std::map<BoundaryID, std::set<MaterialPropertyName>> _zero_boundary_material_props;

  std::set<std::string> _material_property_requested;


  std::map<SubdomainID, std::multimap<std::string, std::string>> _map_block_material_props_check;
  std::map<BoundaryID, std::multimap<std::string, std::string>> _map_boundary_material_props_check;

  std::vector<std::set<MooseVariableFieldBase *>> _active_elemental_moose_variables;

  /* This needs to remain <unsigned int> for threading purposes */
  std::vector<unsigned int> _has_active_elemental_moose_variables;

  std::vector<std::set<TagID>> _active_fe_var_coupleable_matrix_tags;

  std::vector<std::set<TagID>> _active_fe_var_coupleable_vector_tags;

  std::vector<std::set<TagID>> _active_sc_var_coupleable_matrix_tags;

  std::vector<std::set<TagID>> _active_sc_var_coupleable_vector_tags;

  bool _default_ghosting;

  std::set<dof_id_type> _ghosted_elems;

  bool _currently_computing_jacobian;

  bool _currently_computing_residual_and_jacobian;

  bool _computing_nonlinear_residual;

  bool _currently_computing_residual;

  bool _safe_access_tagged_matrices;

  bool _safe_access_tagged_vectors;

  bool _have_ad_objects;

  std::unordered_set<TagID> _not_zeroed_tagged_vectors;

private:
  virtual std::pair<bool, unsigned int>
  determineSolverSystem(const std::string & var_name, bool error_if_not_found = false) const = 0;

  enum class TrueFunctorIs
  {
    UNSET,
    NONAD,
    AD
  };

  std::vector<std::multimap<std::string,
                            std::tuple<TrueFunctorIs,
                                       std::unique_ptr<Moose::FunctorEnvelopeBase>,
                                       std::unique_ptr<Moose::FunctorEnvelopeBase>>>>
      _functors;

  std::vector<std::map<std::string, std::unique_ptr<Moose::FunctorAbstract>>> _pbblf_functors;

  void showFunctors() const;

  void showFunctorRequestors() const;

  std::map<std::string, std::set<std::string>> _functor_to_requestors;

  std::vector<std::multimap<std::string, std::pair<bool, bool>>> _functor_to_request_info;

  bool _output_functors;

  std::vector<VectorTag> _vector_tags;

  std::vector<std::vector<VectorTag>> _typed_vector_tags;

  std::map<TagName, TagID> _vector_tags_name_map;

  std::string restrictionSubdomainCheckName(SubdomainID check_id);
  std::string restrictionBoundaryCheckName(BoundaryID check_id);

  // Contains properties consumed by objects, see addConsumedPropertyName
  std::map<MooseObjectName, std::set<std::string>> _consumed_material_properties;

  std::unordered_map<libMesh::GhostingFunctor *,
                     std::vector<std::shared_ptr<libMesh::GhostingFunctor>>>
      _root_alg_gf_to_sys_clones;

  std::unordered_map<libMesh::GhostingFunctor *,
                     std::vector<std::shared_ptr<libMesh::GhostingFunctor>>>
      _root_coupling_gf_to_sys_clones;

  bool _have_p_refinement;

  std::unordered_map<FEFamily, bool> _family_for_p_refinement;
  static const std::unordered_set<FEFamily> _default_families_without_p_refinement;

  friend class Restartable;
};

template <typename T>
const Moose::Functor<T> &
SubProblem::getFunctor(const std::string & name,
                       const THREAD_ID tid,
                       const std::string & requestor_name,
                       const bool requestor_is_ad)
{
  mooseAssert(tid < _functors.size(), "Too large a thread ID");

  // Log the requestor
  _functor_to_requestors["wraps_" + name].insert(requestor_name);

  constexpr bool requested_functor_is_ad =
      !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;

  auto & functor_to_request_info = _functor_to_request_info[tid];

  // Get the requested functor if we already have it
  auto & functors = _functors[tid];
  if (auto find_ret = functors.find("wraps_" + name); find_ret != functors.end())
  {
    if (functors.count("wraps_" + name) > 1)
      mooseError("Attempted to get a functor with the name '",
                 name,
                 "' but multiple (" + std::to_string(functors.count("wraps_" + name)) +
                     ") functors match. Make sure that you do not have functor material "
                     "properties, functions, postprocessors or variables with the same names.");

    auto & [true_functor_is, non_ad_functor, ad_functor] = find_ret->second;
    auto & functor_wrapper = requested_functor_is_ad ? *ad_functor : *non_ad_functor;

    auto * const functor = dynamic_cast<Moose::Functor<T> *>(&functor_wrapper);
    if (!functor)
      mooseError("A call to SubProblem::getFunctor requested a functor named '",
                 name,
                 "' that returns the type: '",
                 libMesh::demangle(typeid(T).name()),
                 "'. However, that functor already exists and returns a different type: '",
                 functor_wrapper.returnType(),
                 "'");

    if (functor->template wrapsType<Moose::NullFunctor<T>>())
      // Store for future checking when the actual functor gets added
      functor_to_request_info.emplace(name,
                                      std::make_pair(requested_functor_is_ad, requestor_is_ad));
    else
    {
      // We already have the actual functor
      if (true_functor_is == SubProblem::TrueFunctorIs::UNSET)
        mooseError("We already have the functor; it should not be unset");

      // Check for whether this is a valid request
      // We allow auxiliary variables and linear variables to be retrieved as non AD
      if (!requested_functor_is_ad && requestor_is_ad &&
          true_functor_is == SubProblem::TrueFunctorIs::AD &&
          !(hasAuxiliaryVariable(name) || hasLinearVariable(name)))
        mooseError("The AD object '",
                   requestor_name,
                   "' is requesting the functor '",
                   name,
                   "' as a non-AD functor even though it is truly an AD functor, which is not "
                   "allowed, since this may unintentionally drop derivatives.");
    }

    return *functor;
  }

  // We don't have the functor yet but we could have it in the future. We'll create null functors
  // for now
  functor_to_request_info.emplace(name, std::make_pair(requested_functor_is_ad, requestor_is_ad));
  if constexpr (requested_functor_is_ad)
  {
    typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
    typedef T ADType;

    auto emplace_ret =
        functors.emplace("wraps_" + name,
                         std::make_tuple(SubProblem::TrueFunctorIs::UNSET,
                                         std::make_unique<Moose::Functor<NonADType>>(
                                             std::make_unique<Moose::NullFunctor<NonADType>>()),
                                         std::make_unique<Moose::Functor<ADType>>(
                                             std::make_unique<Moose::NullFunctor<ADType>>())));

    return static_cast<Moose::Functor<T> &>(*(requested_functor_is_ad
                                                  ? std::get<2>(emplace_ret->second)
                                                  : std::get<1>(emplace_ret->second)));
  }
  else
  {
    typedef T NonADType;
    typedef typename Moose::ADType<T>::type ADType;

    auto emplace_ret =
        functors.emplace("wraps_" + name,
                         std::make_tuple(SubProblem::TrueFunctorIs::UNSET,
                                         std::make_unique<Moose::Functor<NonADType>>(
                                             std::make_unique<Moose::NullFunctor<NonADType>>()),
                                         std::make_unique<Moose::Functor<ADType>>(
                                             std::make_unique<Moose::NullFunctor<ADType>>())));

    return static_cast<Moose::Functor<T> &>(*(requested_functor_is_ad
                                                  ? std::get<2>(emplace_ret->second)
                                                  : std::get<1>(emplace_ret->second)));
  }
}

template <typename T>
bool
SubProblem::hasFunctorWithType(const std::string & name, const THREAD_ID tid) const
{
  mooseAssert(tid < _functors.size(), "Too large a thread ID");
  auto & functors = _functors[tid];

  const auto & it = functors.find("wraps_" + name);
  constexpr bool requested_functor_is_ad =
      !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;

  if (it == functors.end())
    return false;
  else
    return dynamic_cast<Moose::Functor<T> *>(
        requested_functor_is_ad ? std::get<2>(it->second).get() : std::get<1>(it->second).get());
}

template <typename T, typename PolymorphicLambda>
const Moose::FunctorBase<T> &
SubProblem::addPiecewiseByBlockLambdaFunctor(const std::string & name,
                                             PolymorphicLambda my_lammy,
                                             const std::set<ExecFlagType> & clearance_schedule,
                                             const MooseMesh & mesh,
                                             const std::set<SubdomainID> & block_ids,
                                             const THREAD_ID tid)
{
  auto & pbblf_functors = _pbblf_functors[tid];

  auto [it, first_time_added] =
      pbblf_functors.emplace(name,
                             std::make_unique<PiecewiseByBlockLambdaFunctor<T>>(
                                 name, my_lammy, clearance_schedule, mesh, block_ids));

  auto * functor = dynamic_cast<PiecewiseByBlockLambdaFunctor<T> *>(it->second.get());
  if (!functor)
  {
    if (first_time_added)
      mooseError("This should be impossible. If this was the first time we added the functor, then "
                 "the dynamic cast absolutely should have succeeded");
    else
      mooseError("Attempted to add a lambda functor with the name '",
                 name,
                 "' but another lambda functor of that name returns a different type");
  }

  if (first_time_added)
    addFunctor(name, *functor, tid);
  else
    // The functor already exists
    functor->setFunctor(mesh, block_ids, my_lammy);

  return *functor;
}

template <typename T>
void
SubProblem::addFunctor(const std::string & name,
                       const Moose::FunctorBase<T> & functor,
                       const THREAD_ID tid)
{
  constexpr bool added_functor_is_ad =
      !std::is_same<T, typename MetaPhysicL::RawType<T>::value_type>::value;

  mooseAssert(tid < _functors.size(), "Too large a thread ID");

  auto & functor_to_request_info = _functor_to_request_info[tid];
  auto & functors = _functors[tid];
  auto it = functors.find("wraps_" + name);
  if (it != functors.end())
  {
    // We have this functor already. If it's a null functor, we want to replace it with the valid
    // functor we have now. If it's not then we'll add a new entry into the multimap and then we'll
    // error later if a user requests a functor because their request is ambiguous. This is the
    // reason that the functors container is a multimap: for nice error messages
    auto * const existing_wrapper_base =
        added_functor_is_ad ? std::get<2>(it->second).get() : std::get<1>(it->second).get();
    auto * const existing_wrapper = dynamic_cast<Moose::Functor<T> *>(existing_wrapper_base);
    if (existing_wrapper && existing_wrapper->template wrapsType<Moose::NullFunctor<T>>())
    {
      // Sanity check
      auto [request_info_it, request_info_end_it] = functor_to_request_info.equal_range(name);
      if (request_info_it == request_info_end_it)
        mooseError("We are wrapping a NullFunctor but we don't have any unfilled functor request "
                   "info. This doesn't make sense.");

      // Check for valid requests
      while (request_info_it != request_info_end_it)
      {
        auto & [requested_functor_is_ad, requestor_is_ad] = request_info_it->second;
        if (!requested_functor_is_ad && requestor_is_ad && added_functor_is_ad)
          mooseError("We are requesting a non-AD functor '" + name +
                     "' from an AD object, but the true functor is AD. This means we could be "
                     "dropping important derivatives. We will not allow this");
        // We're going to eventually check whether we've fulfilled all functor requests and our
        // check will be that the multimap is empty. This request is fulfilled, so erase it from the
        // map now
        request_info_it = functor_to_request_info.erase(request_info_it);
      }

      // Ok we didn't have the functor before, so we will add it now
      std::get<0>(it->second) =
          added_functor_is_ad ? SubProblem::TrueFunctorIs::AD : SubProblem::TrueFunctorIs::NONAD;
      existing_wrapper->assign(functor);
      // Finally we create the non-AD or AD complement of the just added functor
      if constexpr (added_functor_is_ad)
      {
        typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
        auto * const existing_non_ad_wrapper_base = std::get<1>(it->second).get();
        auto * const existing_non_ad_wrapper =
            dynamic_cast<Moose::Functor<NonADType> *>(existing_non_ad_wrapper_base);
        mooseAssert(existing_non_ad_wrapper->template wrapsType<Moose::NullFunctor<NonADType>>(),
                    "Both members of pair should have been wrapping a NullFunctor");
        existing_non_ad_wrapper->assign(
            std::make_unique<Moose::RawValueFunctor<NonADType>>(functor));
      }
      else
      {
        typedef typename Moose::ADType<T>::type ADType;
        auto * const existing_ad_wrapper_base = std::get<2>(it->second).get();
        auto * const existing_ad_wrapper =
            dynamic_cast<Moose::Functor<ADType> *>(existing_ad_wrapper_base);
        mooseAssert(existing_ad_wrapper->template wrapsType<Moose::NullFunctor<ADType>>(),
                    "Both members of pair should have been wrapping a NullFunctor");
        existing_ad_wrapper->assign(std::make_unique<Moose::ADWrapperFunctor<ADType>>(functor));
      }
      return;
    }
    else if (!existing_wrapper)
    {
      // Functor was emplaced but the cast failed. This could be a double definition with
      // different types, or it could be a request with one type then a definition with another
      // type. Either way it is going to error later, but it is cleaner to catch it now
      mooseError("Functor '",
                 name,
                 "' is being added with return type '",
                 MooseUtils::prettyCppType<T>(),
                 "' but it has already been defined or requested with return type '",
                 existing_wrapper_base->returnType(),
                 "'.");
    }
  }

  // We are a new functor, create the opposite ADType one and store it with other functors
  if constexpr (added_functor_is_ad)
  {
    typedef typename MetaPhysicL::RawType<T>::value_type NonADType;
    auto new_non_ad_wrapper = std::make_unique<Moose::Functor<NonADType>>(
        std::make_unique<Moose::RawValueFunctor<NonADType>>(functor));
    auto new_ad_wrapper = std::make_unique<Moose::Functor<T>>(functor);
    _functors[tid].emplace("wraps_" + name,
                           std::make_tuple(SubProblem::TrueFunctorIs::AD,
                                           std::move(new_non_ad_wrapper),
                                           std::move(new_ad_wrapper)));
  }
  else
  {
    typedef typename Moose::ADType<T>::type ADType;
    auto new_non_ad_wrapper = std::make_unique<Moose::Functor<T>>((functor));
    auto new_ad_wrapper = std::make_unique<Moose::Functor<ADType>>(
        std::make_unique<Moose::ADWrapperFunctor<ADType>>(functor));
    _functors[tid].emplace("wraps_" + name,
                           std::make_tuple(SubProblem::TrueFunctorIs::NONAD,
                                           std::move(new_non_ad_wrapper),
                                           std::move(new_ad_wrapper)));
  }
}

inline const bool &
SubProblem::currentlyComputingResidualAndJacobian() const
{
  return _currently_computing_residual_and_jacobian;
}

inline void
SubProblem::setCurrentlyComputingResidualAndJacobian(
    const bool currently_computing_residual_and_jacobian)
{
  _currently_computing_residual_and_jacobian = currently_computing_residual_and_jacobian;
}

namespace Moose
{
void initial_condition(libMesh::EquationSystems & es, const std::string & system_name);
} // namespace Moose