FEProblemBase.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

// MOOSE includes
#include "SubProblem.h"
#include "GeometricSearchData.h"
#include "MeshDivision.h"
#include "MortarData.h"
#include "ReporterData.h"
#include "Adaptivity.h"
#include "InitialConditionWarehouse.h"
#include "FVInitialConditionWarehouse.h"
#include "ScalarInitialConditionWarehouse.h"
#include "Restartable.h"
#include "SolverParams.h"
#include "PetscSupport.h"
#include "MooseApp.h"
#include "ExecuteMooseObjectWarehouse.h"
#include "MaterialWarehouse.h"
#include "MooseVariableFE.h"
#include "MultiAppTransfer.h"
#include "Postprocessor.h"
#include "HashMap.h"
#include "VectorPostprocessor.h"
#include "PerfGraphInterface.h"
#include "Attributes.h"
#include "MooseObjectWarehouse.h"
#include "MaterialPropertyRegistry.h"
#include "RestartableEquationSystems.h"
#include "SolutionInvalidity.h"
#include "PetscSupport.h"

#include "libmesh/enum_quadrature_type.h"
#include "libmesh/equation_systems.h"

#include <unordered_map>
#include <memory>

// Forward declarations
class AuxiliarySystem;
class DisplacedProblem;
class MooseMesh;
class NonlinearSystemBase;
class LinearSystem;
class SolverSystem;
class NonlinearSystem;
class RandomInterface;
class RandomData;
class MeshChangedInterface;
class MeshDisplacedInterface;
class MultiMooseEnum;
class MaterialPropertyStorage;
class MaterialData;
class MooseEnum;
class Assembly;
class JacobianBlock;
class Control;
class MultiApp;
class TransientMultiApp;
class ScalarInitialCondition;
class Indicator;
class InternalSideIndicatorBase;
class Marker;
class Material;
class Transfer;
class XFEMInterface;
class SideUserObject;
class NodalUserObject;
class ElementUserObject;
class InternalSideUserObject;
class InterfaceUserObject;
class GeneralUserObject;
class Positions;
class Function;
class Distribution;
class Sampler;
class KernelBase;
class IntegratedBCBase;
class LineSearch;
class UserObject;
class AutomaticMortarGeneration;
class VectorPostprocessor;
class Convergence;
class MooseAppCoordTransform;
class MortarUserObject;
class SolutionInvalidity;

// libMesh forward declarations
namespace libMesh
{
class CouplingMatrix;
class NonlinearImplicitSystem;
class LinearImplicitSystem;
} // namespace libMesh

enum class MooseLinearConvergenceReason
{
  ITERATING = 0,
  // CONVERGED_RTOL_NORMAL        =  1,
  // CONVERGED_ATOL_NORMAL        =  9,
  CONVERGED_RTOL = 2,
  CONVERGED_ATOL = 3,
  CONVERGED_ITS = 4,
  // CONVERGED_CG_NEG_CURVE       =  5,
  // CONVERGED_CG_CONSTRAINED     =  6,
  // CONVERGED_STEP_LENGTH        =  7,
  // CONVERGED_HAPPY_BREAKDOWN    =  8,
  DIVERGED_NULL = -2,
  // DIVERGED_ITS                 = -3,
  // DIVERGED_DTOL                = -4,
  // DIVERGED_BREAKDOWN           = -5,
  // DIVERGED_BREAKDOWN_BICG      = -6,
  // DIVERGED_NONSYMMETRIC        = -7,
  // DIVERGED_INDEFINITE_PC       = -8,
  DIVERGED_NANORINF = -9,
  // DIVERGED_INDEFINITE_MAT      = -10
  DIVERGED_PCSETUP_FAILED = -11
};

class FEProblemBase : public SubProblem, public Restartable
{
public:
  static InputParameters validParams();

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

  enum class CoverageCheckMode
  {
    FALSE,
    TRUE,
    OFF,
    ON,
    SKIP_LIST,
    ONLY_LIST,
  };

  virtual libMesh::EquationSystems & es() override { return _req.set().es(); }
  virtual MooseMesh & mesh() override { return _mesh; }
  virtual const MooseMesh & mesh() const override { return _mesh; }
  const MooseMesh & mesh(bool use_displaced) const override;

  void setCoordSystem(const std::vector<SubdomainName> & blocks, const MultiMooseEnum & coord_sys);
  void setAxisymmetricCoordAxis(const MooseEnum & rz_coord_axis);

  void setCoupling(Moose::CouplingType type);

  Moose::CouplingType coupling() const { return _coupling; }

  void setCouplingMatrix(std::unique_ptr<libMesh::CouplingMatrix> cm,
                         const unsigned int nl_sys_num);

  // DEPRECATED METHOD
  void setCouplingMatrix(libMesh::CouplingMatrix * cm, const unsigned int nl_sys_num);

  const libMesh::CouplingMatrix * couplingMatrix(const unsigned int nl_sys_num) const override;

  void setNonlocalCouplingMatrix();

  bool
  areCoupled(const unsigned int ivar, const unsigned int jvar, const unsigned int nl_sys_num) const;

  bool hasUOAuxStateCheck() const { return _uo_aux_state_check; }

  bool checkingUOAuxState() const { return _checking_uo_aux_state; }

  void trustUserCouplingMatrix();

  std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> &
  couplingEntries(const THREAD_ID tid, const unsigned int nl_sys_num);
  std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> &
  nonlocalCouplingEntries(const THREAD_ID tid, const unsigned int nl_sys_num);

  virtual bool hasVariable(const std::string & var_name) const override;
  // NOTE: hasAuxiliaryVariable defined in parent class
  bool hasSolverVariable(const std::string & var_name) const;
  using SubProblem::getVariable;
  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 override;
  MooseVariableFieldBase & getActualFieldVariable(const THREAD_ID tid,
                                                  const std::string & var_name) override;
  virtual MooseVariable & getStandardVariable(const THREAD_ID tid,
                                              const std::string & var_name) override;
  virtual VectorMooseVariable & getVectorVariable(const THREAD_ID tid,
                                                  const std::string & var_name) override;
  virtual ArrayMooseVariable & getArrayVariable(const THREAD_ID tid,
                                                const std::string & var_name) override;

  virtual bool hasScalarVariable(const std::string & var_name) const override;
  virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid,
                                                  const std::string & var_name) override;
  virtual libMesh::System & getSystem(const std::string & var_name) override;

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

  virtual void clearActiveElementalMooseVariables(const THREAD_ID tid) override;

  virtual void clearActiveFEVariableCoupleableMatrixTags(const THREAD_ID tid) override;

  virtual void clearActiveFEVariableCoupleableVectorTags(const THREAD_ID tid) override;

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

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

  virtual void clearActiveScalarVariableCoupleableMatrixTags(const THREAD_ID tid) override;

  virtual void clearActiveScalarVariableCoupleableVectorTags(const THREAD_ID tid) override;

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

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

  virtual void createQRules(libMesh::QuadratureType type,
                            libMesh::Order order,
                            libMesh::Order volume_order = libMesh::INVALID_ORDER,
                            libMesh::Order face_order = libMesh::INVALID_ORDER,
                            SubdomainID block = Moose::ANY_BLOCK_ID,
                            bool allow_negative_qweights = true);

  void bumpVolumeQRuleOrder(libMesh::Order order, SubdomainID block);

  void bumpAllQRuleOrder(libMesh::Order order, SubdomainID block);

  unsigned int getMaxQps() const;

  libMesh::Order getMaxScalarOrder() const;

  void checkNonlocalCoupling();
  void checkUserObjectJacobianRequirement(THREAD_ID tid);
  void setVariableAllDoFMap(const std::vector<const MooseVariableFEBase *> & moose_vars);

  const std::vector<const MooseVariableFEBase *> &
  getUserObjectJacobianVariables(const THREAD_ID tid) const
  {
    return _uo_jacobian_moose_vars[tid];
  }

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

  virtual std::vector<VariableName> getVariableNames();

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

  virtual void prepare(const Elem * elem, const THREAD_ID tid) override;
  virtual void prepareFace(const Elem * elem, const THREAD_ID tid) override;
  virtual void prepare(const Elem * elem,
                       unsigned int ivar,
                       unsigned int jvar,
                       const std::vector<dof_id_type> & dof_indices,
                       const THREAD_ID tid) override;

  virtual void setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid) override;
  virtual void
  setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void setNeighborSubdomainID(const Elem * elem, const THREAD_ID tid);
  virtual void prepareAssembly(const THREAD_ID tid) override;

  virtual void addGhostedElem(dof_id_type elem_id) override;
  virtual void addGhostedBoundary(BoundaryID boundary_id) override;
  virtual void ghostGhostedBoundaries() override;

  virtual void sizeZeroes(unsigned int size, const THREAD_ID tid);
  virtual bool reinitDirac(const Elem * elem, const THREAD_ID tid) override;

  virtual void reinitElem(const Elem * elem, const THREAD_ID tid) override;
  virtual void reinitElemPhys(const Elem * elem,
                              const std::vector<Point> & phys_points_in_elem,
                              const THREAD_ID tid) override;
  void reinitElemFace(const Elem * elem, unsigned int side, BoundaryID, const THREAD_ID tid);
  virtual void reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  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) override;
  virtual void reinitNode(const Node * node, const THREAD_ID tid) override;
  virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid) override;
  virtual void reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) override;
  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes,
                                   const THREAD_ID tid) override;
  virtual void reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  unsigned int neighbor_side,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) override;
  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) override;
  virtual void
  reinitElemNeighborAndLowerD(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void reinitScalars(const THREAD_ID tid,
                             bool reinit_for_derivative_reordering = false) override;
  virtual void reinitOffDiagScalars(const THREAD_ID tid) override;

  virtual void getDiracElements(std::set<const Elem *> & elems) override;
  virtual void clearDiracInfo() override;

  virtual void subdomainSetup(SubdomainID subdomain, const THREAD_ID tid);
  virtual void neighborSubdomainSetup(SubdomainID subdomain, const THREAD_ID tid);

  virtual void newAssemblyArray(std::vector<std::shared_ptr<SolverSystem>> & solver_systems);
  virtual void initNullSpaceVectors(const InputParameters & parameters,
                                    std::vector<std::shared_ptr<NonlinearSystemBase>> & nl);

  virtual void init() override;
  virtual void solve(const unsigned int nl_sys_num);

  virtual void solveLinearSystem(const unsigned int linear_sys_num,
                                 const Moose::PetscSupport::PetscOptions * po = nullptr);


  const libMesh::ConstElemRange & getEvaluableElementRange();
  const libMesh::ConstElemRange & getNonlinearEvaluableElementRange();


  const libMesh::ConstElemRange & getCurrentAlgebraicElementRange();
  const libMesh::ConstNodeRange & getCurrentAlgebraicNodeRange();
  const ConstBndNodeRange & getCurrentAlgebraicBndNodeRange();


  void setCurrentAlgebraicElementRange(libMesh::ConstElemRange * range);
  void setCurrentAlgebraicNodeRange(libMesh::ConstNodeRange * range);
  void setCurrentAlgebraicBndNodeRange(ConstBndNodeRange * range);

  virtual void setException(const std::string & message);

  virtual bool hasException() { return _has_exception; }

  virtual void checkExceptionAndStopSolve(bool print_message = true);

  virtual bool solverSystemConverged(const unsigned int solver_sys_num) override;
  virtual unsigned int nNonlinearIterations(const unsigned int nl_sys_num) const override;
  virtual unsigned int nLinearIterations(const unsigned int nl_sys_num) const override;
  virtual Real finalNonlinearResidual(const unsigned int nl_sys_num) const override;
  virtual bool computingPreSMOResidual(const unsigned int nl_sys_num) const override;

  virtual std::string solverTypeString(unsigned int solver_sys_num = 0);

  virtual bool startedInitialSetup() { return _started_initial_setup; }

  virtual void onTimestepBegin() override;
  virtual void onTimestepEnd() override;

  virtual Real & time() const { return _time; }
  virtual Real & timeOld() const { return _time_old; }
  virtual int & timeStep() const { return _t_step; }
  virtual Real & dt() const { return _dt; }
  virtual Real & dtOld() const { return _dt_old; }
  Real getTimeFromStateArg(const Moose::StateArg & state) const;

  virtual void transient(bool trans) { _transient = trans; }
  virtual bool isTransient() const override { return _transient; }

  virtual void addTimeIntegrator(const std::string & type,
                                 const std::string & name,
                                 InputParameters & parameters);
  virtual void
  addPredictor(const std::string & type, const std::string & name, InputParameters & parameters);

  virtual void copySolutionsBackwards();

  virtual void advanceState();

  virtual void restoreSolutions();

  virtual void saveOldSolutions();

  virtual void restoreOldSolutions();

  void needSolutionState(unsigned int oldest_needed, Moose::SolutionIterationType iteration_type);

  virtual void outputStep(ExecFlagType type);

  virtual void postExecute();


  void allowOutput(bool state);
  template <typename T>
  void allowOutput(bool state);

  void forceOutput();

  virtual void initPetscOutputAndSomeSolverSettings();

  Moose::PetscSupport::PetscOptions & getPetscOptions() { return _petsc_options; }

  void logAdd(const std::string & system,
              const std::string & name,
              const std::string & type,
              const InputParameters & params) const;

  // Function /////
  virtual void
  addFunction(const std::string & type, const std::string & name, InputParameters & parameters);
  virtual bool hasFunction(const std::string & name, const THREAD_ID tid = 0);
  virtual Function & getFunction(const std::string & name, const THREAD_ID tid = 0);

  virtual void
  addMeshDivision(const std::string & type, const std::string & name, InputParameters & params);
  MeshDivision & getMeshDivision(const std::string & name, const THREAD_ID tid = 0) const;

  virtual void
  addConvergence(const std::string & type, const std::string & name, InputParameters & parameters);
  virtual Convergence & getConvergence(const std::string & name, const THREAD_ID tid = 0) const;
  virtual const std::vector<std::shared_ptr<Convergence>> &
  getConvergenceObjects(const THREAD_ID tid = 0) const;
  virtual bool hasConvergence(const std::string & name, const THREAD_ID tid = 0) const;
  bool needToAddDefaultNonlinearConvergence() const
  {
    return _need_to_add_default_nonlinear_convergence;
  }
  bool needToAddDefaultMultiAppFixedPointConvergence() const
  {
    return _need_to_add_default_multiapp_fixed_point_convergence;
  }
  bool needToAddDefaultSteadyStateConvergence() const
  {
    return _need_to_add_default_steady_state_convergence;
  }
  void setNeedToAddDefaultNonlinearConvergence()
  {
    _need_to_add_default_nonlinear_convergence = true;
  }
  void setNeedToAddDefaultMultiAppFixedPointConvergence()
  {
    _need_to_add_default_multiapp_fixed_point_convergence = true;
  }
  void setNeedToAddDefaultSteadyStateConvergence()
  {
    _need_to_add_default_steady_state_convergence = true;
  }
  bool hasSetMultiAppFixedPointConvergenceName() const
  {
    return _multiapp_fixed_point_convergence_name.has_value();
  }
  bool hasSetSteadyStateConvergenceName() const
  {
    return _steady_state_convergence_name.has_value();
  }
  virtual void addDefaultNonlinearConvergence(const InputParameters & params);
  virtual bool onlyAllowDefaultNonlinearConvergence() const { return false; }
  void addDefaultMultiAppFixedPointConvergence(const InputParameters & params);
  void addDefaultSteadyStateConvergence(const InputParameters & params);

  virtual void addLineSearch(const InputParameters & /*parameters*/)
  {
    mooseError("Line search not implemented for this problem type yet.");
  }

  virtual void lineSearch();

  LineSearch * getLineSearch() override { return _line_search.get(); }

  virtual void
  addDistribution(const std::string & type, const std::string & name, InputParameters & parameters);
  virtual Distribution & getDistribution(const std::string & name);

  virtual void
  addSampler(const std::string & type, const std::string & name, InputParameters & parameters);
  virtual Sampler & getSampler(const std::string & name, const THREAD_ID tid = 0);

  // NL /////
  NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num);
  const NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num) const;
  void setCurrentNonlinearSystem(const unsigned int nl_sys_num);
  NonlinearSystemBase & currentNonlinearSystem();
  const NonlinearSystemBase & currentNonlinearSystem() const;

  virtual const SystemBase & systemBaseNonlinear(const unsigned int sys_num) const override;
  virtual SystemBase & systemBaseNonlinear(const unsigned int sys_num) override;

  virtual const SystemBase & systemBaseSolver(const unsigned int sys_num) const override;
  virtual SystemBase & systemBaseSolver(const unsigned int sys_num) override;

  virtual const SystemBase & systemBaseAuxiliary() const override;
  virtual SystemBase & systemBaseAuxiliary() override;

  virtual NonlinearSystem & getNonlinearSystem(const unsigned int sys_num);

  virtual const SystemBase & getSystemBase(const unsigned int sys_num) const;

  virtual SystemBase & getSystemBase(const unsigned int sys_num);

  SystemBase & getSystemBase(const std::string & sys_name);

  LinearSystem & getLinearSystem(unsigned int sys_num);

  const LinearSystem & getLinearSystem(unsigned int sys_num) const;

  SolverSystem & getSolverSystem(unsigned int sys_num);

  const SolverSystem & getSolverSystem(unsigned int sys_num) const;

  void setCurrentLinearSystem(unsigned int sys_num);

  LinearSystem & currentLinearSystem();
  const LinearSystem & currentLinearSystem() const;

  virtual const SystemBase & systemBaseLinear(unsigned int sys_num) const override;

  virtual SystemBase & systemBaseLinear(unsigned int sys_num) override;

  virtual void
  addVariable(const std::string & var_type, const std::string & var_name, InputParameters & params);

  virtual void addKernel(const std::string & kernel_name,
                         const std::string & name,
                         InputParameters & parameters);
  virtual void addHDGKernel(const std::string & kernel_name,
                            const std::string & name,
                            InputParameters & parameters);
  virtual void addNodalKernel(const std::string & kernel_name,
                              const std::string & name,
                              InputParameters & parameters);
  virtual void addScalarKernel(const std::string & kernel_name,
                               const std::string & name,
                               InputParameters & parameters);
  virtual void addBoundaryCondition(const std::string & bc_name,
                                    const std::string & name,
                                    InputParameters & parameters);
  virtual void
  addConstraint(const std::string & c_name, const std::string & name, InputParameters & parameters);

  virtual void setInputParametersFEProblem(InputParameters & parameters)
  {
    parameters.set<FEProblemBase *>("_fe_problem_base") = this;
  }

  // Aux /////

  virtual void addAuxVariable(const std::string & var_type,
                              const std::string & var_name,
                              InputParameters & params);

  virtual void addAuxVariable(const std::string & var_name,
                              const libMesh::FEType & type,
                              const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void addAuxArrayVariable(const std::string & var_name,
                                   const libMesh::FEType & type,
                                   unsigned int components,
                                   const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void addAuxScalarVariable(const std::string & var_name,
                                    libMesh::Order order,
                                    Real scale_factor = 1.,
                                    const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void addAuxKernel(const std::string & kernel_name,
                            const std::string & name,
                            InputParameters & parameters);
  virtual void addAuxScalarKernel(const std::string & kernel_name,
                                  const std::string & name,
                                  InputParameters & parameters);

  AuxiliarySystem & getAuxiliarySystem() { return *_aux; }

  // Dirac /////
  virtual void addDiracKernel(const std::string & kernel_name,
                              const std::string & name,
                              InputParameters & parameters);

  // DG /////
  virtual void addDGKernel(const std::string & kernel_name,
                           const std::string & name,
                           InputParameters & parameters);
  // FV /////
  virtual void addFVKernel(const std::string & kernel_name,
                           const std::string & name,
                           InputParameters & parameters);

  virtual void addLinearFVKernel(const std::string & kernel_name,
                                 const std::string & name,
                                 InputParameters & parameters);
  virtual void
  addFVBC(const std::string & fv_bc_name, const std::string & name, InputParameters & parameters);
  virtual void addLinearFVBC(const std::string & fv_bc_name,
                             const std::string & name,
                             InputParameters & parameters);

  virtual void addFVInterfaceKernel(const std::string & fv_ik_name,
                                    const std::string & name,
                                    InputParameters & parameters);

  // Interface /////
  virtual void addInterfaceKernel(const std::string & kernel_name,
                                  const std::string & name,
                                  InputParameters & parameters);

  // IC /////
  virtual void addInitialCondition(const std::string & ic_name,
                                   const std::string & name,
                                   InputParameters & parameters);
  virtual void addFVInitialCondition(const std::string & ic_name,
                                     const std::string & name,
                                     InputParameters & parameters);

  void projectSolution();

  unsigned short getCurrentICState();

  void projectInitialConditionOnCustomRange(
      libMesh::ConstElemRange & elem_range,
      ConstBndNodeRange & bnd_node_range,
      const std::optional<std::set<VariableName>> & target_vars = std::nullopt);

  void projectFunctionOnCustomRange(ConstElemRange & elem_range,
                                    Number (*func)(const Point &,
                                                   const libMesh::Parameters &,
                                                   const std::string &,
                                                   const std::string &),
                                    Gradient (*func_grad)(const Point &,
                                                          const libMesh::Parameters &,
                                                          const std::string &,
                                                          const std::string &),
                                    const libMesh::Parameters & params,
                                    const VariableName & target_var);

  // Materials
  virtual void addMaterial(const std::string & material_name,
                           const std::string & name,
                           InputParameters & parameters);
  virtual void addMaterialHelper(std::vector<MaterialWarehouse *> warehouse,
                                 const std::string & material_name,
                                 const std::string & name,
                                 InputParameters & parameters);
  virtual void addInterfaceMaterial(const std::string & material_name,
                                    const std::string & name,
                                    InputParameters & parameters);
  virtual void addFunctorMaterial(const std::string & functor_material_name,
                                  const std::string & name,
                                  InputParameters & parameters);

  void prepareMaterials(const std::unordered_set<unsigned int> & consumer_needed_mat_props,
                        const SubdomainID blk_id,
                        const THREAD_ID tid);

  void reinitMaterials(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful = true);

  void reinitMaterialsFace(SubdomainID blk_id,
                           const THREAD_ID tid,
                           bool swap_stateful = true,
                           const std::deque<MaterialBase *> * reinit_mats = nullptr);

  void reinitMaterialsNeighbor(SubdomainID blk_id,
                               const THREAD_ID tid,
                               bool swap_stateful = true,
                               const std::deque<MaterialBase *> * reinit_mats = nullptr);

  void reinitMaterialsBoundary(BoundaryID boundary_id,
                               const THREAD_ID tid,
                               bool swap_stateful = true,
                               const std::deque<MaterialBase *> * reinit_mats = nullptr);

  void
  reinitMaterialsInterface(BoundaryID boundary_id, const THREAD_ID tid, bool swap_stateful = true);

  /*
   * Swap back underlying data storing stateful material properties
   */
  virtual void swapBackMaterials(const THREAD_ID tid);
  virtual void swapBackMaterialsFace(const THREAD_ID tid);
  virtual void swapBackMaterialsNeighbor(const THREAD_ID tid);

  void setActiveMaterialProperties(const std::unordered_set<unsigned int> & mat_prop_ids,
                                   const THREAD_ID tid);

  bool hasActiveMaterialProperties(const THREAD_ID tid) const;

  void clearActiveMaterialProperties(const THREAD_ID tid);

  template <typename T>
  std::vector<std::shared_ptr<T>> addObject(const std::string & type,
                                            const std::string & name,
                                            InputParameters & parameters,
                                            const bool threaded = true,
                                            const std::string & var_param_name = "variable");

  // Postprocessors /////
  virtual void addPostprocessor(const std::string & pp_name,
                                const std::string & name,
                                InputParameters & parameters);

  // VectorPostprocessors /////
  virtual void addVectorPostprocessor(const std::string & pp_name,
                                      const std::string & name,
                                      InputParameters & parameters);

  virtual void
  addReporter(const std::string & type, const std::string & name, InputParameters & parameters);

  const ReporterData & getReporterData() const { return _reporter_data; }

  ReporterData & getReporterData(ReporterData::WriteKey /*key*/) { return _reporter_data; }

  // UserObjects /////
  virtual std::vector<std::shared_ptr<UserObject>> addUserObject(
      const std::string & user_object_name, const std::string & name, InputParameters & parameters);

  // TODO: delete this function after apps have been updated to not call it
  const ExecuteMooseObjectWarehouse<UserObject> & getUserObjects() const
  {
    mooseDeprecated(
        "This function is deprecated, use theWarehouse().query() to construct a query instead");
    return _all_user_objects;
  }

  template <class T>
  T & getUserObject(const std::string & name, unsigned int tid = 0) const
  {
    std::vector<T *> objs;
    theWarehouse()
        .query()
        .condition<AttribSystem>("UserObject")
        .condition<AttribThread>(tid)
        .condition<AttribName>(name)
        .queryInto(objs);
    if (objs.empty())
      mooseError("Unable to find user object with name '" + name + "'");
    return *(objs[0]);
  }
  const UserObject & getUserObjectBase(const std::string & name, const THREAD_ID tid = 0) const;

  const Positions & getPositionsObject(const std::string & name) const;

  bool hasUserObject(const std::string & name) const;

  bool hasPostprocessorValueByName(const PostprocessorName & name) const;

  const PostprocessorValue & getPostprocessorValueByName(const PostprocessorName & name,
                                                         std::size_t t_index = 0) const;

  void setPostprocessorValueByName(const PostprocessorName & name,
                                   const PostprocessorValue & value,
                                   std::size_t t_index = 0);

  bool hasPostprocessor(const std::string & name) const;

  const VectorPostprocessorValue &
  getVectorPostprocessorValueByName(const std::string & object_name,
                                    const std::string & vector_name,
                                    std::size_t t_index = 0) const;

  void setVectorPostprocessorValueByName(const std::string & object_name,
                                         const std::string & vector_name,
                                         const VectorPostprocessorValue & value,
                                         std::size_t t_index = 0);

  const VectorPostprocessor & getVectorPostprocessorObjectByName(const std::string & object_name,
                                                                 const THREAD_ID tid = 0) const;


  bool hasMultiApps() const { return _multi_apps.hasActiveObjects(); }
  bool hasMultiApps(ExecFlagType type) const;
  bool hasMultiApp(const std::string & name) const;

  // Dampers /////
  virtual void addDamper(const std::string & damper_name,
                         const std::string & name,
                         InputParameters & parameters);
  void setupDampers();

  bool hasDampers() { return _has_dampers; }

  // Indicators /////
  virtual void addIndicator(const std::string & indicator_name,
                            const std::string & name,
                            InputParameters & parameters);

  // Markers //////
  virtual void addMarker(const std::string & marker_name,
                         const std::string & name,
                         InputParameters & parameters);

  virtual void addMultiApp(const std::string & multi_app_name,
                           const std::string & name,
                           InputParameters & parameters);

  std::shared_ptr<MultiApp> getMultiApp(const std::string & multi_app_name) const;

  std::vector<std::shared_ptr<Transfer>> getTransfers(ExecFlagType type,
                                                      Transfer::DIRECTION direction) const;
  std::vector<std::shared_ptr<Transfer>> getTransfers(Transfer::DIRECTION direction) const;

  const ExecuteMooseObjectWarehouse<Transfer> &
  getMultiAppTransferWarehouse(Transfer::DIRECTION direction) const;

  void execMultiAppTransfers(ExecFlagType type, Transfer::DIRECTION direction);

  bool execMultiApps(ExecFlagType type, bool auto_advance = true);

  void finalizeMultiApps();

  void incrementMultiAppTStep(ExecFlagType type);

  void advanceMultiApps(ExecFlagType type)
  {
    mooseDeprecated("Deprecated method; use finishMultiAppStep and/or incrementMultiAppTStep "
                    "depending on your purpose");
    finishMultiAppStep(type);
  }

  void finishMultiAppStep(ExecFlagType type, bool recurse_through_multiapp_levels = false);

  void backupMultiApps(ExecFlagType type);

  void restoreMultiApps(ExecFlagType type, bool force = false);

  Real computeMultiAppsDT(ExecFlagType type);

  virtual void addTransfer(const std::string & transfer_name,
                           const std::string & name,
                           InputParameters & parameters);

  void execTransfers(ExecFlagType type);

  Real computeResidualL2Norm(NonlinearSystemBase & sys);

  Real computeResidualL2Norm(LinearSystem & sys);

  virtual Real computeResidualL2Norm();

  virtual void computeResidualSys(libMesh::NonlinearImplicitSystem & sys,
                                  const NumericVector<libMesh::Number> & soln,
                                  NumericVector<libMesh::Number> & residual);
  void computeResidual(libMesh::NonlinearImplicitSystem & sys,
                       const NumericVector<libMesh::Number> & soln,
                       NumericVector<libMesh::Number> & residual);

  virtual void computeResidual(const NumericVector<libMesh::Number> & soln,
                               NumericVector<libMesh::Number> & residual,
                               const unsigned int nl_sys_num);

  void computeResidualAndJacobian(const NumericVector<libMesh::Number> & soln,
                                  NumericVector<libMesh::Number> & residual,
                                  libMesh::SparseMatrix<libMesh::Number> & jacobian);

  virtual void computeResidualTag(const NumericVector<libMesh::Number> & soln,
                                  NumericVector<libMesh::Number> & residual,
                                  TagID tag);
  virtual void computeResidualType(const NumericVector<libMesh::Number> & soln,
                                   NumericVector<libMesh::Number> & residual,
                                   TagID tag);

  virtual void computeResidualInternal(const NumericVector<libMesh::Number> & soln,
                                       NumericVector<libMesh::Number> & residual,
                                       const std::set<TagID> & tags);
  virtual void computeResidualTags(const std::set<TagID> & tags);

  virtual void computeJacobianSys(libMesh::NonlinearImplicitSystem & sys,
                                  const NumericVector<libMesh::Number> & soln,
                                  libMesh::SparseMatrix<libMesh::Number> & jacobian);
  virtual void computeJacobian(const NumericVector<libMesh::Number> & soln,
                               libMesh::SparseMatrix<libMesh::Number> & jacobian,
                               const unsigned int nl_sys_num);

  virtual void computeJacobianTag(const NumericVector<libMesh::Number> & soln,
                                  libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                  TagID tag);

  virtual void computeJacobianInternal(const NumericVector<libMesh::Number> & soln,
                                       libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                       const std::set<TagID> & tags);

  virtual void computeJacobianTags(const std::set<TagID> & tags);

  virtual void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
                                     const unsigned int nl_sys_num);

  virtual void computeJacobianBlock(libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                    libMesh::System & precond_system,
                                    unsigned int ivar,
                                    unsigned int jvar);

  virtual void computeLinearSystemSys(libMesh::LinearImplicitSystem & sys,
                                      libMesh::SparseMatrix<libMesh::Number> & system_matrix,
                                      NumericVector<libMesh::Number> & rhs,
                                      const bool compute_gradients = true);

  void computeLinearSystemTags(const NumericVector<libMesh::Number> & soln,
                               const std::set<TagID> & vector_tags,
                               const std::set<TagID> & matrix_tags,
                               const bool compute_gradients = true);

  virtual Real computeDamping(const NumericVector<libMesh::Number> & soln,
                              const NumericVector<libMesh::Number> & update);

  virtual bool shouldUpdateSolution();

  virtual bool updateSolution(NumericVector<libMesh::Number> & vec_solution,
                              NumericVector<libMesh::Number> & ghosted_solution);

  virtual void predictorCleanup(NumericVector<libMesh::Number> & ghosted_solution);

  virtual void computeBounds(libMesh::NonlinearImplicitSystem & sys,
                             NumericVector<libMesh::Number> & lower,
                             NumericVector<libMesh::Number> & upper);
  virtual void computeNearNullSpace(libMesh::NonlinearImplicitSystem & sys,
                                    std::vector<NumericVector<libMesh::Number> *> & sp);
  virtual void computeNullSpace(libMesh::NonlinearImplicitSystem & sys,
                                std::vector<NumericVector<libMesh::Number> *> & sp);
  virtual void computeTransposeNullSpace(libMesh::NonlinearImplicitSystem & sys,
                                         std::vector<NumericVector<libMesh::Number> *> & sp);
  virtual void computePostCheck(libMesh::NonlinearImplicitSystem & sys,
                                const NumericVector<libMesh::Number> & old_soln,
                                NumericVector<libMesh::Number> & search_direction,
                                NumericVector<libMesh::Number> & new_soln,
                                bool & changed_search_direction,
                                bool & changed_new_soln);

  virtual void computeIndicatorsAndMarkers();
  virtual void computeIndicators();
  virtual void computeMarkers();

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

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

  virtual void addCachedResidualDirectly(NumericVector<libMesh::Number> & residual,
                                         const THREAD_ID tid);

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

  virtual void addJacobian(const THREAD_ID tid) override;
  virtual void addJacobianNeighbor(const THREAD_ID tid) override;
  virtual void addJacobianNeighborLowerD(const THREAD_ID tid) override;
  virtual void addJacobianLowerD(const THREAD_ID tid) override;
  virtual void addJacobianBlockTags(libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                    unsigned int ivar,
                                    unsigned int jvar,
                                    const DofMap & dof_map,
                                    std::vector<dof_id_type> & dof_indices,
                                    const std::set<TagID> & tags,
                                    const THREAD_ID tid);
  virtual void addJacobianNeighbor(libMesh::SparseMatrix<libMesh::Number> & jacobian,
                                   unsigned int ivar,
                                   unsigned int jvar,
                                   const 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) override;
  virtual void addJacobianScalar(const THREAD_ID tid = 0);
  virtual void addJacobianOffDiagScalar(unsigned int ivar, const THREAD_ID tid = 0);

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

  virtual void prepareShapes(unsigned int var, const THREAD_ID tid) override;
  virtual void prepareFaceShapes(unsigned int var, const THREAD_ID tid) override;
  virtual void prepareNeighborShapes(unsigned int var, const THREAD_ID tid) override;

  // Displaced problem /////
  virtual void addDisplacedProblem(std::shared_ptr<DisplacedProblem> displaced_problem);
  virtual std::shared_ptr<const DisplacedProblem> getDisplacedProblem() const
  {
    return _displaced_problem;
  }
  virtual std::shared_ptr<DisplacedProblem> getDisplacedProblem() { return _displaced_problem; }

  virtual void updateGeomSearch(
      GeometricSearchData::GeometricSearchType type = GeometricSearchData::ALL) override;
  virtual void updateMortarMesh();

  void createMortarInterface(
      const std::pair<BoundaryID, BoundaryID> & primary_secondary_boundary_pair,
      const std::pair<SubdomainID, SubdomainID> & primary_secondary_subdomain_pair,
      bool on_displaced,
      bool periodic,
      const bool debug,
      const bool correct_edge_dropping,
      const Real minimum_projection_angle);

  const AutomaticMortarGeneration &
  getMortarInterface(const std::pair<BoundaryID, BoundaryID> & primary_secondary_boundary_pair,
                     const std::pair<SubdomainID, SubdomainID> & primary_secondary_subdomain_pair,
                     bool on_displaced) const;

  AutomaticMortarGeneration &
  getMortarInterface(const std::pair<BoundaryID, BoundaryID> & primary_secondary_boundary_pair,
                     const std::pair<SubdomainID, SubdomainID> & primary_secondary_subdomain_pair,
                     bool on_displaced);

  const std::unordered_map<std::pair<BoundaryID, BoundaryID>, AutomaticMortarGeneration> &
  getMortarInterfaces(bool on_displaced) const;

  virtual void possiblyRebuildGeomSearchPatches();

  virtual GeometricSearchData & geomSearchData() override { return _geometric_search_data; }

  void setRestartFile(const std::string & file_name);

  const MaterialPropertyRegistry & getMaterialPropertyRegistry() const
  {
    return _material_prop_registry;
  }

  const MaterialPropertyStorage & getMaterialPropertyStorage() { return _material_props; }
  const MaterialPropertyStorage & getBndMaterialPropertyStorage() { return _bnd_material_props; }
  const MaterialPropertyStorage & getNeighborMaterialPropertyStorage()
  {
    return _neighbor_material_props;
  }

  const MooseObjectWarehouse<Indicator> & getIndicatorWarehouse() { return _indicators; }
  const MooseObjectWarehouse<InternalSideIndicatorBase> & getInternalSideIndicatorWarehouse()
  {
    return _internal_side_indicators;
  }
  const MooseObjectWarehouse<Marker> & getMarkerWarehouse() { return _markers; }

  const InitialConditionWarehouse & getInitialConditionWarehouse() const { return _ics; }

  const FVInitialConditionWarehouse & getFVInitialConditionWarehouse() const { return _fv_ics; }

  SolverParams & solverParams(unsigned int solver_sys_num = 0);

  const SolverParams & solverParams(unsigned int solver_sys_num = 0) const;

#ifdef LIBMESH_ENABLE_AMR
  // Adaptivity /////
  Adaptivity & adaptivity() { return _adaptivity; }
  virtual void initialAdaptMesh();

  virtual bool adaptMesh();

  unsigned int getNumCyclesCompleted() { return _cycles_completed; }

  bool hasInitialAdaptivity() const { return _adaptivity.getInitialSteps() > 0; }
#else

  bool hasInitialAdaptivity() const { return false; }
#endif // LIBMESH_ENABLE_AMR

  void initXFEM(std::shared_ptr<XFEMInterface> xfem);

  std::shared_ptr<XFEMInterface> getXFEM() { return _xfem; }

  bool haveXFEM() { return _xfem != nullptr; }

  virtual bool updateMeshXFEM();

  virtual void
  meshChanged(bool intermediate_change, bool contract_mesh, bool clean_refinement_flags);

  void notifyWhenMeshChanges(MeshChangedInterface * mci);

  void notifyWhenMeshDisplaces(MeshDisplacedInterface * mdi);

  void initElementStatefulProps(const libMesh::ConstElemRange & elem_range, const bool threaded);

  virtual void checkProblemIntegrity();

  void registerRandomInterface(RandomInterface & random_interface, const std::string & name);

  void setConstJacobian(bool state) { _const_jacobian = state; }

  void setKernelCoverageCheck(CoverageCheckMode mode) { _kernel_coverage_check = mode; }

  void setKernelCoverageCheck(bool flag)
  {
    _kernel_coverage_check = flag ? CoverageCheckMode::TRUE : CoverageCheckMode::FALSE;
  }

  void setMaterialCoverageCheck(CoverageCheckMode mode) { _material_coverage_check = mode; }

  void setMaterialCoverageCheck(bool flag)
  {
    _material_coverage_check = flag ? CoverageCheckMode::TRUE : CoverageCheckMode::FALSE;
  }

  void setParallelBarrierMessaging(bool flag) { _parallel_barrier_messaging = flag; }

  void setVerboseProblem(bool verbose);

  bool verboseMultiApps() const { return _verbose_multiapps; }

  void parentOutputPositionChanged();


  bool needBoundaryMaterialOnSide(BoundaryID bnd_id, const THREAD_ID tid);
  bool needInterfaceMaterialOnSide(BoundaryID bnd_id, const THREAD_ID tid);
  bool needSubdomainMaterialOnSide(SubdomainID subdomain_id, const THREAD_ID tid);

  unsigned int subspaceDim(const std::string & prefix) const
  {
    if (_subspace_dim.count(prefix))
      return _subspace_dim.find(prefix)->second;
    else
      return 0;
  }

  /*
   * Return a reference to the material warehouse of *all* Material objects.
   */
  const MaterialWarehouse & getMaterialWarehouse() const { return _all_materials; }

  /*
   * Return a reference to the material warehouse of Material objects to be computed.
   */
  const MaterialWarehouse & getRegularMaterialsWarehouse() const { return _materials; }
  const MaterialWarehouse & getDiscreteMaterialWarehouse() const { return _discrete_materials; }
  const MaterialWarehouse & getInterfaceMaterialsWarehouse() const { return _interface_materials; }

  std::shared_ptr<MaterialBase> getMaterial(std::string name,
                                            Moose::MaterialDataType type,
                                            const THREAD_ID tid = 0,
                                            bool no_warn = false);

  /*
   * @return The MaterialData for the type \p type for thread \p tid
   */
  MaterialData & getMaterialData(Moose::MaterialDataType type, const THREAD_ID tid = 0) const;

  bool restoreOriginalNonzeroPattern() const { return _restore_original_nonzero_pattern; }

  bool errorOnJacobianNonzeroReallocation() const
  {
    return _error_on_jacobian_nonzero_reallocation;
  }

  void setErrorOnJacobianNonzeroReallocation(bool state)
  {
    _error_on_jacobian_nonzero_reallocation = state;
  }

  bool preserveMatrixSparsityPattern() const { return _preserve_matrix_sparsity_pattern; };

  void setPreserveMatrixSparsityPattern(bool preserve);

  bool ignoreZerosInJacobian() const { return _ignore_zeros_in_jacobian; }

  void setIgnoreZerosInJacobian(bool state) { _ignore_zeros_in_jacobian = state; }

  bool acceptInvalidSolution() const;
  bool allowInvalidSolution() const { return _allow_invalid_solution; }

  bool showInvalidSolutionConsole() const { return _show_invalid_solution_console; }

  bool immediatelyPrintInvalidSolution() const { return _immediately_print_invalid_solution; }

  bool hasTimeIntegrator() const { return _has_time_integrator; }


  const ExecFlagType & getCurrentExecuteOnFlag() const;
  void setCurrentExecuteOnFlag(const ExecFlagType &);

  virtual void execute(const ExecFlagType & exec_type);
  virtual void executeAllObjects(const ExecFlagType & exec_type);

  virtual Executor & getExecutor(const std::string & name) { return _app.getExecutor(name); }

  virtual void computeUserObjects(const ExecFlagType & type, const Moose::AuxGroup & group);

  virtual void computeUserObjectByName(const ExecFlagType & type,
                                       const Moose::AuxGroup & group,
                                       const std::string & name);

  void needsPreviousNewtonIteration(bool state);

  bool needsPreviousNewtonIteration() const;


  std::vector<Real> _real_zero;
  std::vector<VariableValue> _scalar_zero;
  std::vector<VariableValue> _zero;
  std::vector<VariablePhiValue> _phi_zero;
  std::vector<MooseArray<ADReal>> _ad_zero;
  std::vector<VariableGradient> _grad_zero;
  std::vector<MooseArray<ADRealVectorValue>> _ad_grad_zero;
  std::vector<VariablePhiGradient> _grad_phi_zero;
  std::vector<VariableSecond> _second_zero;
  std::vector<MooseArray<ADRealTensorValue>> _ad_second_zero;
  std::vector<VariablePhiSecond> _second_phi_zero;
  std::vector<Point> _point_zero;
  std::vector<VectorVariableValue> _vector_zero;
  std::vector<VectorVariableCurl> _vector_curl_zero;

  ExecuteMooseObjectWarehouse<Control> & getControlWarehouse() { return _control_warehouse; }

  void executeControls(const ExecFlagType & exec_type);

  void executeSamplers(const ExecFlagType & exec_type);

  virtual void updateActiveObjects();

  void reportMooseObjectDependency(MooseObject * a, MooseObject * b);

  ExecuteMooseObjectWarehouse<MultiApp> & getMultiAppWarehouse() { return _multi_apps; }

  bool hasJacobian() const;

  bool constJacobian() const;

  void addOutput(const std::string &, const std::string &, InputParameters &);

  inline TheWarehouse & theWarehouse() const { return _app.theWarehouse(); }

  void setSNESMFReuseBase(bool reuse, bool set_by_user)
  {
    _snesmf_reuse_base = reuse, _snesmf_reuse_base_set_by_user = set_by_user;
  }

  bool useSNESMFReuseBase() { return _snesmf_reuse_base; }

  void skipExceptionCheck(bool skip_exception_check)
  {
    _skip_exception_check = skip_exception_check;
  }

  bool isSNESMFReuseBaseSetbyUser() { return _snesmf_reuse_base_set_by_user; }

  bool & petscOptionsInserted() { return _is_petsc_options_inserted; }

#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
  PetscOptions & petscOptionsDatabase() { return _petsc_option_data_base; }
#endif

  virtual void setUDotRequested(const bool u_dot_requested) { _u_dot_requested = u_dot_requested; }

  virtual void setUDotDotRequested(const bool u_dotdot_requested)
  {
    _u_dotdot_requested = u_dotdot_requested;
  }

  virtual void setUDotOldRequested(const bool u_dot_old_requested)
  {
    _u_dot_old_requested = u_dot_old_requested;
  }

  virtual void setUDotDotOldRequested(const bool u_dotdot_old_requested)
  {
    _u_dotdot_old_requested = u_dotdot_old_requested;
  }

  virtual bool uDotRequested() { return _u_dot_requested; }

  virtual bool uDotDotRequested() { return _u_dotdot_requested; }

  virtual bool uDotOldRequested()
  {
    if (_u_dot_old_requested && !_u_dot_requested)
      mooseError("FEProblemBase: When requesting old time derivative of solution, current time "
                 "derivative of solution should also be stored. Please set `u_dot_requested` to "
                 "true using setUDotRequested.");

    return _u_dot_old_requested;
  }

  virtual bool uDotDotOldRequested()
  {
    if (_u_dotdot_old_requested && !_u_dotdot_requested)
      mooseError("FEProblemBase: When requesting old second time derivative of solution, current "
                 "second time derivation of solution should also be stored. Please set "
                 "`u_dotdot_requested` to true using setUDotDotRequested.");
    return _u_dotdot_old_requested;
  }

  using SubProblem::haveADObjects;
  void haveADObjects(bool have_ad_objects) override;

  // Whether or not we should solve this system
  bool shouldSolve() const { return _solve; }

  const MortarData & mortarData() const { return _mortar_data; }
  MortarData & mortarData() { return _mortar_data; }

  virtual bool hasNeighborCoupling() const { return _has_internal_edge_residual_objects; }

  virtual bool hasMortarCoupling() const { return _has_mortar; }

  using SubProblem::computingNonlinearResid;
  void computingNonlinearResid(bool computing_nonlinear_residual) final;

  using SubProblem::currentlyComputingResidual;
  void setCurrentlyComputingResidual(bool currently_computing_residual) final;

  void numGridSteps(unsigned int num_grid_steps) { _num_grid_steps = num_grid_steps; }

  void uniformRefine();

  using SubProblem::automaticScaling;
  void automaticScaling(bool automatic_scaling) override;


  template <typename T>
  static void objectSetupHelper(const std::vector<T *> & objects, const ExecFlagType & exec_flag);
  template <typename T>
  static void objectExecuteHelper(const std::vector<T *> & objects);

  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) override;

  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) override;

  bool fvBCsIntegrityCheck() const { return _fv_bcs_integrity_check; }

  void fvBCsIntegrityCheck(bool fv_bcs_integrity_check);

  void getFVMatsAndDependencies(SubdomainID block_id,
                                std::vector<std::shared_ptr<MaterialBase>> & face_materials,
                                std::vector<std::shared_ptr<MaterialBase>> & neighbor_materials,
                                std::set<MooseVariableFieldBase *> & variables,
                                const THREAD_ID tid);

  void resizeMaterialData(Moose::MaterialDataType data_type, unsigned int nqp, const THREAD_ID tid);

  bool haveDisplaced() const override final { return _displaced_problem.get(); }

  bool hasLinearConvergenceObjects() const;
  void setNonlinearConvergenceNames(const std::vector<ConvergenceName> & convergence_names);
  void setLinearConvergenceNames(const std::vector<ConvergenceName> & convergence_names);
  void setMultiAppFixedPointConvergenceName(const ConvergenceName & convergence_name);
  void setSteadyStateConvergenceName(const ConvergenceName & convergence_name);

  const std::vector<ConvergenceName> & getNonlinearConvergenceNames() const;
  const std::vector<ConvergenceName> & getLinearConvergenceNames() const;
  const ConvergenceName & getMultiAppFixedPointConvergenceName() const;
  const ConvergenceName & getSteadyStateConvergenceName() const;

  void computingScalingJacobian(bool computing_scaling_jacobian)
  {
    _computing_scaling_jacobian = computing_scaling_jacobian;
  }

  bool computingScalingJacobian() const override final { return _computing_scaling_jacobian; }

  void computingScalingResidual(bool computing_scaling_residual)
  {
    _computing_scaling_residual = computing_scaling_residual;
  }

  bool computingScalingResidual() const override final { return _computing_scaling_residual; }

  MooseAppCoordTransform & coordTransform();

  virtual std::size_t numNonlinearSystems() const override { return _num_nl_sys; }

  virtual std::size_t numLinearSystems() const override { return _num_linear_sys; }

  virtual std::size_t numSolverSystems() const override { return _num_nl_sys + _num_linear_sys; }

  bool isSolverSystemNonlinear(const unsigned int sys_num) { return sys_num < _num_nl_sys; }

  virtual unsigned int currentNlSysNum() const override;

  virtual unsigned int currentLinearSysNum() const override;

  virtual unsigned int nlSysNum(const NonlinearSystemName & nl_sys_name) const override;

  unsigned int linearSysNum(const LinearSystemName & linear_sys_name) const override;

  unsigned int solverSysNum(const SolverSystemName & solver_sys_name) const override;

  unsigned int systemNumForVariable(const VariableName & variable_name) const;

  bool getFailNextNonlinearConvergenceCheck() const { return getFailNextSystemConvergenceCheck(); }
  bool getFailNextSystemConvergenceCheck() const { return _fail_next_system_convergence_check; }

  void setFailNextNonlinearConvergenceCheck() { setFailNextSystemConvergenceCheck(); }
  void setFailNextSystemConvergenceCheck() { _fail_next_system_convergence_check = true; }

  void resetFailNextNonlinearConvergenceCheck() { resetFailNextSystemConvergenceCheck(); }
  void resetFailNextSystemConvergenceCheck() { _fail_next_system_convergence_check = false; }

  /*
   * Set the status of loop order of execution printing
   * @param print_exec set of execution flags to print on
   */
  void setExecutionPrinting(const ExecFlagEnum & print_exec) { _print_execution_on = print_exec; }

  bool shouldPrintExecution(const THREAD_ID tid) const;
  void reinitMortarUserObjects(BoundaryID primary_boundary_id,
                               BoundaryID secondary_boundary_id,
                               bool displaced);

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

  class CurrentResidualVectorTagsKey
  {
    friend class CrankNicolson;
    friend class FEProblemBase;
    CurrentResidualVectorTagsKey() {}
    CurrentResidualVectorTagsKey(const CurrentResidualVectorTagsKey &) {}
  };

  void setCurrentResidualVectorTags(const std::set<TagID> & vector_tags);

  void clearCurrentResidualVectorTags();

  void clearCurrentJacobianMatrixTags() {}

  virtual void needFV() override { _have_fv = true; }
  virtual bool haveFV() const override { return _have_fv; }

  virtual bool hasNonlocalCoupling() const override { return _has_nonlocal_coupling; }

  bool identifyVariableGroupsInNL() const { return _identify_variable_groups_in_nl; }

  virtual void setCurrentLowerDElem(const Elem * const lower_d_elem, const THREAD_ID tid) override;
  virtual void setCurrentBoundaryID(BoundaryID bid, const THREAD_ID tid) override;

  const std::vector<NonlinearSystemName> & getNonlinearSystemNames() const { return _nl_sys_names; }
  const std::vector<LinearSystemName> & getLinearSystemNames() const { return _linear_sys_names; }
  const std::vector<SolverSystemName> & getSolverSystemNames() const { return _solver_sys_names; }

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

  virtual bool checkNonlocalCouplingRequirement() const override;

  virtual Moose::FEBackend feBackend() const { return Moose::FEBackend::LibMesh; }

  class CreateTaggedMatrixKey
  {
    CreateTaggedMatrixKey() {}
    CreateTaggedMatrixKey(const CreateTaggedMatrixKey &) {}

    friend class AddTaggedMatricesAction;
  };

  void createTagMatrices(CreateTaggedMatrixKey);

protected:
  virtual void meshChanged() {}

  void createTagVectors();

  void createTagSolutions();

  virtual void meshDisplaced();

  void computeSystems(const ExecFlagType & type);

  MooseMesh & _mesh;

private:
  Restartable::ManagedValue<RestartableEquationSystems> _req;

  void setResidualObjectParamsAndLog(const std::string & ro_name,
                                     const std::string & name,
                                     InputParameters & params,
                                     const unsigned int nl_sys_num,
                                     const std::string & base_name,
                                     bool & reinit_displaced);

  static SolverParams makeLinearSolverParams();

protected:
  bool _initialized;

  std::optional<std::vector<ConvergenceName>> _nonlinear_convergence_names;
  std::optional<std::vector<ConvergenceName>> _linear_convergence_names;
  std::optional<ConvergenceName> _multiapp_fixed_point_convergence_name;
  std::optional<ConvergenceName> _steady_state_convergence_name;

  std::set<TagID> _fe_vector_tags;

  std::set<TagID> _fe_matrix_tags;

  std::set<TagID> _linear_vector_tags;

  std::set<TagID> _linear_matrix_tags;

  const bool & _solve;

  bool _transient;
  Real & _time;
  Real & _time_old;
  int & _t_step;
  Real & _dt;
  Real & _dt_old;

  bool _need_to_add_default_nonlinear_convergence;
  bool _need_to_add_default_multiapp_fixed_point_convergence;
  bool _need_to_add_default_steady_state_convergence;

  const std::vector<LinearSystemName> _linear_sys_names;

  const std::size_t _num_linear_sys;

  std::vector<std::shared_ptr<LinearSystem>> _linear_systems;

  std::map<LinearSystemName, unsigned int> _linear_sys_name_to_num;

  LinearSystem * _current_linear_sys;

  const bool _using_default_nl;

  const std::vector<NonlinearSystemName> _nl_sys_names;

  const std::size_t _num_nl_sys;

  std::vector<std::shared_ptr<NonlinearSystemBase>> _nl;

  std::map<NonlinearSystemName, unsigned int> _nl_sys_name_to_num;

  NonlinearSystemBase * _current_nl_sys;

  SolverSystem * _current_solver_sys;

  std::vector<std::shared_ptr<SolverSystem>> _solver_systems;

  std::map<SolverVariableName, unsigned int> _solver_var_to_sys_num;

  std::map<SolverSystemName, unsigned int> _solver_sys_name_to_num;

  std::vector<SolverSystemName> _solver_sys_names;

  std::shared_ptr<AuxiliarySystem> _aux;

  Moose::CouplingType _coupling;                             
  std::vector<std::unique_ptr<libMesh::CouplingMatrix>> _cm; 

  std::map<std::string, unsigned int> _subspace_dim;

  std::vector<std::vector<std::unique_ptr<Assembly>>> _assembly;

  MooseObjectWarehouse<MeshDivision> _mesh_divisions;

  MooseObjectWarehouse<Function> _functions;

  MooseObjectWarehouse<Convergence> _convergences;

  MooseObjectWarehouse<KernelBase> _nonlocal_kernels;

  MooseObjectWarehouse<IntegratedBCBase> _nonlocal_integrated_bcs;

  InitialConditionWarehouse _ics;
  FVInitialConditionWarehouse _fv_ics;
  ScalarInitialConditionWarehouse _scalar_ics; // use base b/c of setup methods

  // material properties
  MaterialPropertyRegistry _material_prop_registry;
  MaterialPropertyStorage & _material_props;
  MaterialPropertyStorage & _bnd_material_props;
  MaterialPropertyStorage & _neighbor_material_props;

  // Material Warehouses
  MaterialWarehouse _materials;           // regular materials
  MaterialWarehouse _interface_materials; // interface materials
  MaterialWarehouse _discrete_materials;  // Materials that the user must compute
  MaterialWarehouse _all_materials; // All materials for error checking and MaterialData storage

  // Indicator Warehouses
  MooseObjectWarehouse<Indicator> _indicators;
  MooseObjectWarehouse<InternalSideIndicatorBase> _internal_side_indicators;

  // Marker Warehouse
  MooseObjectWarehouse<Marker> _markers;

  // Helper class to access Reporter object values
  ReporterData _reporter_data;

  // TODO: delete this after apps have been updated to not call getUserObjects
  ExecuteMooseObjectWarehouse<UserObject> _all_user_objects;

  ExecuteMooseObjectWarehouse<MultiApp> _multi_apps;

  ExecuteMooseObjectWarehouse<TransientMultiApp> _transient_multi_apps;

  ExecuteMooseObjectWarehouse<Transfer> _transfers;

  ExecuteMooseObjectWarehouse<Transfer> _to_multi_app_transfers;

  ExecuteMooseObjectWarehouse<Transfer> _from_multi_app_transfers;

  ExecuteMooseObjectWarehouse<Transfer> _between_multi_app_transfers;

  std::map<std::string, std::unique_ptr<RandomData>> _random_data_objects;

  std::vector<std::unordered_map<SubdomainID, bool>> _block_mat_side_cache;

  std::vector<std::unordered_map<BoundaryID, bool>> _bnd_mat_side_cache;

  std::vector<std::unordered_map<BoundaryID, bool>> _interface_mat_side_cache;

  std::vector<MeshChangedInterface *> _notify_when_mesh_changes;

  std::vector<MeshDisplacedInterface *> _notify_when_mesh_displaces;

  bool duplicateVariableCheck(const std::string & var_name,
                              const libMesh::FEType & type,
                              bool is_aux,
                              const std::set<SubdomainID> * const active_subdomains);

  void computeUserObjectsInternal(const ExecFlagType & type,
                                  const Moose::AuxGroup & group,
                                  TheWarehouse::Query & query);

  void checkDisplacementOrders();

  void checkUserObjects();

  void checkDependMaterialsHelper(
      const std::map<SubdomainID, std::vector<std::shared_ptr<MaterialBase>>> & materials_map);

  void checkCoordinateSystems();

  void reinitBecauseOfGhostingOrNewGeomObjects(bool mortar_changed = false);

  void addObjectParamsHelper(InputParameters & params,
                             const std::string & object_name,
                             const std::string & var_param_name = "variable");

#ifdef LIBMESH_ENABLE_AMR
  Adaptivity _adaptivity;
  unsigned int _cycles_completed;
#endif

  std::shared_ptr<XFEMInterface> _xfem;

  // Displaced mesh /////
  MooseMesh * _displaced_mesh;
  std::shared_ptr<DisplacedProblem> _displaced_problem;
  GeometricSearchData _geometric_search_data;
  MortarData _mortar_data;

  bool _reinit_displaced_elem;
  bool _reinit_displaced_face;
  bool _reinit_displaced_neighbor;

  bool _input_file_saved;

  bool _has_dampers;

  bool _has_constraints;

  bool _snesmf_reuse_base;

  bool _skip_exception_check;

  bool _snesmf_reuse_base_set_by_user;

  bool _has_initialized_stateful;

  bool _const_jacobian;

  bool _has_jacobian;

  bool _needs_old_newton_iter;

  bool _previous_nl_solution_required;

  bool _has_nonlocal_coupling;
  bool _calculate_jacobian_in_uo;

  std::vector<std::vector<const MooseVariableFEBase *>> _uo_jacobian_moose_vars;

  std::vector<unsigned char> _has_active_material_properties;

  std::vector<SolverParams> _solver_params;

  CoverageCheckMode _kernel_coverage_check;
  std::vector<SubdomainName> _kernel_coverage_blocks;

  const bool _boundary_restricted_node_integrity_check;

  const bool _boundary_restricted_elem_integrity_check;

  CoverageCheckMode _material_coverage_check;
  std::vector<SubdomainName> _material_coverage_blocks;

  bool _fv_bcs_integrity_check;

  const bool _material_dependency_check;

  const bool _uo_aux_state_check;

  unsigned int _max_qps;

  libMesh::Order _max_scalar_order;

  bool _has_time_integrator;

  bool _has_exception;

  bool _parallel_barrier_messaging;

  MooseEnum _verbose_setup;

  bool _verbose_multiapps;

  bool _verbose_restore;

  std::string _exception_message;

  ExecFlagType _current_execute_on_flag;

  ExecuteMooseObjectWarehouse<Control> _control_warehouse;

  Moose::PetscSupport::PetscOptions _petsc_options;
#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
  PetscOptions _petsc_option_data_base;
#endif

  bool _is_petsc_options_inserted;

  std::shared_ptr<LineSearch> _line_search;

  std::unique_ptr<libMesh::ConstElemRange> _evaluable_local_elem_range;
  std::unique_ptr<libMesh::ConstElemRange> _nl_evaluable_local_elem_range;
  std::unique_ptr<libMesh::ConstElemRange> _aux_evaluable_local_elem_range;

  std::unique_ptr<libMesh::ConstElemRange> _current_algebraic_elem_range;
  std::unique_ptr<libMesh::ConstNodeRange> _current_algebraic_node_range;
  std::unique_ptr<ConstBndNodeRange> _current_algebraic_bnd_node_range;

  bool _using_ad_mat_props;

  // loop state during projection of initial conditions
  unsigned short _current_ic_state;

  const bool _use_hash_table_matrix_assembly;

private:
  void handleException(const std::string & calling_method);

  std::vector<MortarUserObject *>
  getMortarUserObjects(BoundaryID primary_boundary_id,
                       BoundaryID secondary_boundary_id,
                       bool displaced,
                       const std::vector<MortarUserObject *> & mortar_uo_superset);

  std::vector<MortarUserObject *> getMortarUserObjects(BoundaryID primary_boundary_id,
                                                       BoundaryID secondary_boundary_id,
                                                       bool displaced);

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

  void checkICRestartError(const std::string & ic_name,
                           const std::string & name,
                           const VariableName & var_name);

  /*
   * Test if stateful property redistribution is expected to be
   * necessary, and set it up if so.
   */
  void addAnyRedistributers();

  void updateMaxQps();

  void joinAndFinalize(TheWarehouse::Query query, bool isgen = false);

  virtual void resetState();

  // Parameters handling Jacobian sparsity pattern behavior
  bool _error_on_jacobian_nonzero_reallocation;
  const bool _restore_original_nonzero_pattern;
  bool _ignore_zeros_in_jacobian;
  bool _preserve_matrix_sparsity_pattern;

  const bool _force_restart;
  const bool _allow_ics_during_restart;
  const bool _skip_nl_system_check;
  bool _fail_next_system_convergence_check;
  const bool _allow_invalid_solution;
  const bool _show_invalid_solution_console;
  const bool & _immediately_print_invalid_solution;

  bool _started_initial_setup;

  bool _has_internal_edge_residual_objects;

  bool _u_dot_requested;

  bool _u_dotdot_requested;

  bool _u_dot_old_requested;

  bool _u_dotdot_old_requested;

  friend class AuxiliarySystem;
  friend class NonlinearSystemBase;
  friend class MooseEigenSystem;
  friend class Resurrector;
  friend class Restartable;
  friend class DisplacedProblem;

  bool _has_mortar;

  unsigned int _num_grid_steps;

  bool _trust_user_coupling_matrix = false;

  bool _computing_scaling_jacobian = false;

  bool _computing_scaling_residual = false;

  bool _checking_uo_aux_state = false;

  ExecFlagEnum _print_execution_on;

  const bool _identify_variable_groups_in_nl;

  std::vector<VectorTag> _current_residual_vector_tags;

  bool _have_fv = false;

  const bool _regard_general_exceptions_as_errors;

  std::vector<libMesh::CouplingMatrix> _nonlocal_cm;

  bool _requires_nonlocal_coupling;

  friend void Moose::PetscSupport::setSinglePetscOption(const std::string & name,
                                                        const std::string & value,
                                                        FEProblemBase * const problem);
};

using FVProblemBase = FEProblemBase;

template <typename T>
void
FEProblemBase::allowOutput(bool state)
{
  _app.getOutputWarehouse().allowOutput<T>(state);
}

template <typename T>
void
FEProblemBase::objectSetupHelper(const std::vector<T *> & objects, const ExecFlagType & exec_flag)
{
  if (exec_flag == EXEC_INITIAL)
  {
    for (T * obj_ptr : objects)
      obj_ptr->initialSetup();
  }

  else if (exec_flag == EXEC_TIMESTEP_BEGIN)
  {
    for (const auto obj_ptr : objects)
      obj_ptr->timestepSetup();
  }
  else if (exec_flag == EXEC_SUBDOMAIN)
  {
    for (const auto obj_ptr : objects)
      obj_ptr->subdomainSetup();
  }

  else if (exec_flag == EXEC_NONLINEAR)
  {
    for (const auto obj_ptr : objects)
      obj_ptr->jacobianSetup();
  }

  else if (exec_flag == EXEC_LINEAR)
  {
    for (const auto obj_ptr : objects)
      obj_ptr->residualSetup();
  }
}

template <typename T>
void
FEProblemBase::objectExecuteHelper(const std::vector<T *> & objects)
{
  for (T * obj_ptr : objects)
    obj_ptr->execute();
}

template <typename T>
std::vector<std::shared_ptr<T>>
FEProblemBase::addObject(const std::string & type,
                         const std::string & name,
                         InputParameters & parameters,
                         const bool threaded,
                         const std::string & var_param_name)
{
  parallel_object_only();

  logAdd(MooseUtils::prettyCppType<T>(), name, type, parameters);
  // Add the _subproblem and _sys parameters depending on use_displaced_mesh
  addObjectParamsHelper(parameters, name, var_param_name);

  const auto n_threads = threaded ? libMesh::n_threads() : 1;
  std::vector<std::shared_ptr<T>> objects(n_threads);
  for (THREAD_ID tid = 0; tid < n_threads; ++tid)
  {
    std::shared_ptr<T> obj = _factory.create<T>(type, name, parameters, tid);
    theWarehouse().add(obj);
    objects[tid] = std::move(obj);
  }

  return objects;
}

inline NonlinearSystemBase &
FEProblemBase::getNonlinearSystemBase(const unsigned int sys_num)
{
  mooseAssert(sys_num < _nl.size(), "System number greater than the number of nonlinear systems");
  return *_nl[sys_num];
}

inline const NonlinearSystemBase &
FEProblemBase::getNonlinearSystemBase(const unsigned int sys_num) const
{
  mooseAssert(sys_num < _nl.size(), "System number greater than the number of nonlinear systems");
  return *_nl[sys_num];
}

inline SolverSystem &
FEProblemBase::getSolverSystem(const unsigned int sys_num)
{
  mooseAssert(sys_num < _solver_systems.size(),
              "System number greater than the number of solver systems");
  return *_solver_systems[sys_num];
}

inline const SolverSystem &
FEProblemBase::getSolverSystem(const unsigned int sys_num) const
{
  mooseAssert(sys_num < _solver_systems.size(),
              "System number greater than the number of solver systems");
  return *_solver_systems[sys_num];
}

inline NonlinearSystemBase &
FEProblemBase::currentNonlinearSystem()
{
  mooseAssert(_current_nl_sys, "The nonlinear system is not currently set");
  return *_current_nl_sys;
}

inline const NonlinearSystemBase &
FEProblemBase::currentNonlinearSystem() const
{
  mooseAssert(_current_nl_sys, "The nonlinear system is not currently set");
  return *_current_nl_sys;
}

inline LinearSystem &
FEProblemBase::getLinearSystem(const unsigned int sys_num)
{
  mooseAssert(sys_num < _linear_systems.size(),
              "System number greater than the number of linear systems");
  return *_linear_systems[sys_num];
}

inline const LinearSystem &
FEProblemBase::getLinearSystem(const unsigned int sys_num) const
{
  mooseAssert(sys_num < _linear_systems.size(),
              "System number greater than the number of linear systems");
  return *_linear_systems[sys_num];
}

inline LinearSystem &
FEProblemBase::currentLinearSystem()
{
  mooseAssert(_current_linear_sys, "The linear system is not currently set");
  return *_current_linear_sys;
}

inline const LinearSystem &
FEProblemBase::currentLinearSystem() const
{
  mooseAssert(_current_linear_sys, "The linear system is not currently set");
  return *_current_linear_sys;
}

inline Assembly &
FEProblemBase::assembly(const THREAD_ID tid, const unsigned int sys_num)
{
  mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
  mooseAssert(sys_num < _assembly[tid].size(),
              "System number larger than the assembly container size");
  return *_assembly[tid][sys_num];
}

inline const Assembly &
FEProblemBase::assembly(const THREAD_ID tid, const unsigned int sys_num) const
{
  mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
  mooseAssert(sys_num < _assembly[tid].size(),
              "System number larger than the assembly container size");
  return *_assembly[tid][sys_num];
}

inline const libMesh::CouplingMatrix *
FEProblemBase::couplingMatrix(const unsigned int i) const
{
  return _cm[i].get();
}

inline void
FEProblemBase::fvBCsIntegrityCheck(const bool fv_bcs_integrity_check)
{
  if (!_fv_bcs_integrity_check)
    // the user has requested that we don't check integrity so we will honor that
    return;

  _fv_bcs_integrity_check = fv_bcs_integrity_check;
}

inline const std::vector<VectorTag> &
FEProblemBase::currentResidualVectorTags() const
{
  return _current_residual_vector_tags;
}

inline void
FEProblemBase::setCurrentResidualVectorTags(const std::set<TagID> & vector_tags)
{
  _current_residual_vector_tags = getVectorTags(vector_tags);
}

inline void
FEProblemBase::clearCurrentResidualVectorTags()
{
  _current_residual_vector_tags.clear();
}