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
 
#ifdef MOOSE_KOKKOS_ENABLED
#include "KokkosAssembly.h"
#include "KokkosFESystem.h"
#endif
 
// MOOSE includes
#include "SubProblem.h"
#include "GeometricSearchData.h"
#include "MeshDivision.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 "MortarInterfaceWarehouse.h"
#include "Mortar3DSubpatchPlane.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 MortarInterfaceWarehouse;
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 UserObjectBase;
class FVInterpolationMethod;
class FVFaceInterpolationMethod;
class FVAdvectedInterpolationMethod;
class AutomaticMortarGeneration;
class VectorPostprocessor;
class Convergence;
class MooseAppCoordTransform;
class MortarUserObject;
class SolutionInvalidity;
 
namespace Moose
{
class FunctionBase;
}
 
#ifdef MOOSE_KOKKOS_ENABLED
namespace Moose::Kokkos
{
class MaterialPropertyStorage;
class Function;
class UserObject;
}
#endif
 
// 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();
 
  [[nodiscard]] bool initialized() const { return _initialized; }
 
  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;
  MooseMesh & mesh(bool use_displaced);
 
  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; }
 
#ifndef NDEBUG
  virtual bool checkResidualForNans() const override { return _check_residual_for_nans; }
 
  void setCheckResidualForNans(bool check_residual_for_nans)
  {
    _check_residual_for_nans = check_residual_for_nans;
  }
#endif
 
  void trustUserCouplingMatrix();
 
  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
  couplingEntries(const THREAD_ID tid, const unsigned int nl_sys_num);
  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
  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;
 
  const RestartableEquationSystems & getRestartableEquationSystems() const;
 
  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;
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::Assembly & kokkosAssembly() { return _kokkos_assembly; }
  const Moose::Kokkos::Assembly & kokkosAssembly() const { return _kokkos_assembly; }
#endif
 
  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 prepareAssemblyNeighbor(const THREAD_ID tid);
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void initKokkos();
#endif
 
  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();
 
  void skipNextForwardSolutionCopyToOld();
 
  virtual void advanceState();
 
  virtual void restoreSolutions();
 
  virtual void saveOldSolutions();
 
  virtual void restoreOldSolutions();
 
  void needSolutionState(unsigned int oldest_needed, Moose::SolutionIterationType iteration_type);
 
  bool hasSolutionState(unsigned int state, Moose::SolutionIterationType iteration_type) const;
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosFunction(const std::string & type,
                                 const std::string & name,
                                 InputParameters & parameters);
  virtual bool hasKokkosFunction(const std::string & name) const;
  virtual Moose::Kokkos::Function getKokkosFunction(const std::string & name);
  template <typename T>
  T & getKokkosFunction(const std::string & name);
#endif
 
  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 bool hasDistribution(const std::string & name) const;
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::Array<Moose::Kokkos::System> & getKokkosSystems() { return _kokkos_systems; }
  const Moose::Kokkos::Array<Moose::Kokkos::System> & getKokkosSystems() const
  {
    return _kokkos_systems;
  }
 
  Moose::Kokkos::Array<Moose::Kokkos::FESystem> & getKokkosFESystems()
  {
    return _kokkos_fe_systems;
  }
  const Moose::Kokkos::Array<Moose::Kokkos::FESystem> & getKokkosFESystems() const
  {
    return _kokkos_fe_systems;
  }
 
  Moose::Kokkos::System & getKokkosSystem(const unsigned int sys_num);
  const Moose::Kokkos::System & getKokkosSystem(const unsigned int sys_num) const;
 
  Moose::Kokkos::FESystem & getKokkosFESystem(const unsigned int sys_num);
  const Moose::Kokkos::FESystem & getKokkosFESystem(const unsigned int sys_num) const;
#endif
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosKernel(const std::string & kernel_name,
                               const std::string & name,
                               InputParameters & parameters);
  virtual void addKokkosNodalKernel(const std::string & kernel_name,
                                    const std::string & name,
                                    InputParameters & parameters);
  virtual void addKokkosBoundaryCondition(const std::string & bc_name,
                                          const std::string & name,
                                          InputParameters & parameters);
  virtual void addKokkosLinearFVKernel(const std::string & kernel_name,
                                       const std::string & name,
                                       InputParameters & parameters);
  virtual void addKokkosLinearFVBC(const std::string & bc_name,
                                   const std::string & name,
                                   InputParameters & parameters);
#endif
 
  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 addElementalFieldVariable(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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosAuxKernel(const std::string & kernel_name,
                                  const std::string & name,
                                  InputParameters & parameters);
#endif
 
  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 std::vector<VariableName> & target_vars);
 
  // 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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosMaterial(const std::string & material_name,
                                 const std::string & name,
                                 InputParameters & parameters);
#endif
 
  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
  reinitMaterialsFaceOnBoundary(const BoundaryID boundary_id,
                                const SubdomainID blk_id,
                                const THREAD_ID tid,
                                const bool swap_stateful = true,
                                const std::deque<MaterialBase *> * const reinit_mats = nullptr);
 
  void
  reinitMaterialsNeighborOnBoundary(const BoundaryID boundary_id,
                                    const SubdomainID blk_id,
                                    const THREAD_ID tid,
                                    const bool swap_stateful = true,
                                    const std::deque<MaterialBase *> * const 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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void prepareKokkosMaterials(const std::unordered_set<unsigned int> & consumer_needed_mat_props);
  void reinitKokkosMaterials();
#endif
 
  /*
   * 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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosPostprocessor(const std::string & pp_name,
                                      const std::string & name,
                                      InputParameters & parameters);
  virtual void addKokkosVectorPostprocessor(const std::string & pp_name,
                                            const std::string & name,
                                            InputParameters & parameters);
  virtual void addKokkosReporter(const std::string & type,
                                 const std::string & name,
                                 InputParameters & parameters);
#endif
 
  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);
 
  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;
 
  bool hasUserObject(const std::string & name) const;
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosUserObject(const std::string & user_object_name,
                                   const std::string & name,
                                   InputParameters & parameters);
 
  template <class T>
  const T & getKokkosUserObject(const std::string & name) const
  {
    std::vector<T *> objs;
    theWarehouse()
        .query()
        .condition<AttribSystem>("KokkosUserObject")
        .condition<AttribName>(name)
        .queryInto(objs);
    if (objs.empty())
      mooseError("Unable to find Kokkos user object with name '" + name + "'");
    return *(objs[0]);
  }
 
  bool hasKokkosUserObject(const std::string & name) const;
#endif
 
  void checkUserObjectNameCollision(const std::string & name, const std::string & type) const;
 
  const Positions & getPositionsObject(const std::string & name) const;
 
  virtual void addFVInterpolationMethod(const std::string & method_type,
                                        const std::string & name,
                                        InputParameters & parameters);
 
  const FVInterpolationMethod & getFVInterpolationMethod(const InterpolationMethodName & name,
                                                         const THREAD_ID tid = 0) const;
 
  const FVFaceInterpolationMethod &
  getFVFaceInterpolationMethod(const InterpolationMethodName & name, const THREAD_ID tid = 0) const;
 
  const FVAdvectedInterpolationMethod &
  getFVAdvectedInterpolationMethod(const InterpolationMethodName & name,
                                   const THREAD_ID tid = 0) const;
 
  bool hasFVInterpolationMethod(const InterpolationMethodName & name) const;
 
  bool hasPostprocessorValueByName(const PostprocessorName & name) const;
 
  const Postprocessor & getPostprocessorObjectByName(const PostprocessorName & object_name,
                                                     const THREAD_ID tid = 0) 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 Mortar3DSubpatchPlane mortar_3d_subpatch_plane,
      const MooseEnum & triangulation,
      const bool triangulate_triangles,
      const Mortar3DQuadraturePointMapping mortar_3d_qp_mapping =
          Mortar3DQuadraturePointMapping::NORMAL_PROJECTION);
 
  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>, MortarInterfaceConfig> &
  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;
  }
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::MaterialPropertyStorage & getKokkosMaterialPropertyStorage()
  {
    return _kokkos_material_props;
  }
  Moose::Kokkos::MaterialPropertyStorage & getKokkosBndMaterialPropertyStorage()
  {
    return _kokkos_bnd_material_props;
  }
  Moose::Kokkos::MaterialPropertyStorage & getKokkosNeighborMaterialPropertyStorage()
  {
    return _kokkos_neighbor_material_props;
  }
#endif
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void initKokkosStatefulProps();
#endif
 
  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 needInternalNeighborSideMaterial(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 reference to function warehouse.
   */
  const MooseObjectWarehouse<Function> & getFunctionWarehouse() { return _functions; }
 
  /*
   * 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; }
 
#ifdef MOOSE_KOKKOS_ENABLED
  /*
   * Return a reference to the material warehouse of Kokkos Material objects to be computed.
   */
  const MaterialWarehouse & getKokkosMaterialsWarehouse() const { return _kokkos_materials; }
#endif
 
  std::shared_ptr<MaterialBase> getMaterial(std::string name,
                                            Moose::MaterialDataType type,
                                            const THREAD_ID tid = 0,
                                            bool no_warn = false);
 
  MaterialData & getMaterialData(Moose::MaterialDataType type,
                                 const THREAD_ID tid = 0,
                                 const MooseObject * object = nullptr) const;
 
#ifdef MOOSE_KOKKOS_ENABLED
  MaterialData & getKokkosMaterialData(Moose::MaterialDataType type,
                                       const MooseObject * object = nullptr) const;
#endif
 
  const std::set<const MooseObject *> &
  getMaterialPropertyStorageConsumers(Moose::MaterialDataType type) const;
 
#ifdef MOOSE_KOKKOS_ENABLED
  const std::set<const MooseObject *> &
  getKokkosMaterialPropertyStorageConsumers(Moose::MaterialDataType type) const;
#endif
 
  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;
 
  void needsPreviousMultiAppFixedPointIterationSolution(bool needed,
                                                        const unsigned int solver_sys_num);
 
  bool needsPreviousMultiAppFixedPointIterationSolution(const unsigned int solver_sys_num) const;
 
  void needsPreviousMultiAppFixedPointIterationAuxiliary(bool state);
 
  bool needsPreviousMultiAppFixedPointIterationAuxiliary() const;
 
  void needsPreviousMultiSystemFixedPointIterationSolution(bool needed,
                                                           const unsigned int solver_sys_num);
 
  bool needsPreviousMultiSystemFixedPointIterationSolution(const unsigned int solver_sys_num) const;
 
  void needsPreviousMultiSystemFixedPointIterationAuxiliary(bool state);
 
  bool needsPreviousMultiSystemFixedPointIterationAuxiliary() 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 MortarInterfaceWarehouse & mortarData() const { return *_mortar_data; }
  MortarInterfaceWarehouse & 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; }
 
  bool sideUOInterfaceMatPropIntegrityCheck() const
  {
    return _side_uo_interface_mat_prop_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);
 
  bool useHashTableMatrixAssembly() const { return _use_hash_table_matrix_assembly; }
 
#ifdef MOOSE_KOKKOS_ENABLED
  bool hasKokkosObjects() const { return _has_kokkos_objects; }
  bool hasKokkosResidualObjects() const { return _has_kokkos_residual_objects; }
  void addKokkosMeshInitializationHook(std::function<void()> function)
  {
    _kokkos_mesh_initialization_hooks.push_back(function);
  }
#endif
 
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 & parameters,
                                     const unsigned int nl_sys_num,
                                     const std::string & base_name,
                                     bool & reinit_displaced);
 
  void setAuxKernelParamsAndLog(const std::string & ak_name,
                                const std::string & name,
                                InputParameters & parameters,
                                const std::string & base_name);
 
  static SolverParams makeLinearSolverParams();
 
  TheWarehouse::Query getUOQuery(const std::string & system,
                                 const ExecFlagType & type,
                                 const Moose::AuxGroup & group) const;
 
  void getUOExecutionGroups(TheWarehouse::Query & query, std::set<int> & execution_groups) const;
 
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; 
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::Array<Moose::Kokkos::System> _kokkos_systems;
  Moose::Kokkos::Array<Moose::Kokkos::FESystem> _kokkos_fe_systems;
#endif
 
  std::map<std::string, unsigned int> _subspace_dim;
 
  std::vector<std::vector<std::unique_ptr<Assembly>>> _assembly;
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::Assembly _kokkos_assembly;
#endif
 
  MooseObjectWarehouse<MeshDivision> _mesh_divisions;
 
  MooseObjectWarehouse<Function> _functions;
 
#ifdef MOOSE_KOKKOS_ENABLED
  MooseObjectWarehouse<Moose::FunctionBase> _kokkos_functions;
#endif
 
  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;
 
#ifdef MOOSE_KOKKOS_ENABLED
  Moose::Kokkos::MaterialPropertyStorage & _kokkos_material_props;
  Moose::Kokkos::MaterialPropertyStorage & _kokkos_bnd_material_props;
  Moose::Kokkos::MaterialPropertyStorage & _kokkos_neighbor_material_props;
#endif
  // 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
 
#ifdef MOOSE_KOKKOS_ENABLED
  MaterialWarehouse _kokkos_materials; // Kokkos materials
#endif
 
  // 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;
 
  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, TheWarehouse::Query & query);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void computeKokkosUserObjectsInternal(const ExecFlagType & type, TheWarehouse::Query & query);
#endif
 
  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;
  std::unique_ptr<MortarInterfaceWarehouse> _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;
  std::vector<bool> _previous_multiapp_fp_nl_solution_required;
  bool _previous_multiapp_fp_aux_solution_required;
  std::vector<bool> _previous_multisystem_fp_nl_solution_required;
  bool _previous_multisystem_fp_aux_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;
 
  const bool _side_uo_interface_mat_prop_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;
 
#ifndef NDEBUG
  bool _check_residual_for_nans;
#endif
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void kokkosJoinAndFinalize(const std::vector<Moose::Kokkos::UserObject *> & userobjs);
#endif
 
  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;
 
#ifdef MOOSE_KOKKOS_ENABLED
  bool _has_kokkos_objects = false;
 
  bool _has_kokkos_residual_objects = false;
 
  std::vector<std::function<void()>> _kokkos_mesh_initialization_hooks;
#endif
 
  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();
}
 
#ifdef MOOSE_KOKKOS_ENABLED
template <typename T>
T &
FEProblemBase::getKokkosFunction(const std::string & name)
{
  if (!hasKokkosFunction(name))
  {
    // If we didn't find a function, it might be a default function, attempt to construct one now
    std::istringstream ss(name);
    Real real_value;
 
    // First see if it's just a constant. If it is, build a ConstantFunction
    if (ss >> real_value && ss.eof())
    {
      InputParameters params = _factory.getValidParams("KokkosConstantFunction");
      params.set<Real>("value") = real_value;
      addKokkosFunction("KokkosConstantFunction", ss.str(), params);
    }
 
    // Try once more
    if (!hasKokkosFunction(name))
      mooseError("Unable to find Kokkos function '" + name, "'");
  }
 
  auto * const ret = dynamic_cast<T *>(_kokkos_functions.getActiveObject(name).get());
  if (!ret)
    mooseError("No Kokkos function named '", name, "' of appropriate type");
 
  return *ret;
}
#endif