FEProblemBase.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* 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 "MortarData.h"
#include "PostprocessorData.h"
#include "VectorPostprocessorData.h"
#include "Adaptivity.h"
#include "InitialConditionWarehouse.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 "libmesh/enum_quadrature_type.h"
#include "libmesh/equation_systems.h"

#include <unordered_map>

// Forward declarations
class AuxiliarySystem;
class DisplacedProblem;
class FEProblemBase;
class MooseMesh;
class NonlinearSystemBase;
class NonlinearSystem;
class RandomInterface;
class RandomData;
class MeshChangedInterface;
class MultiMooseEnum;
class MaterialPropertyStorage;
class MaterialData;
class MooseEnum;
class Resurrector;
class Assembly;
class JacobianBlock;
class Control;
class MultiApp;
class TransientMultiApp;
class ScalarInitialCondition;
class Indicator;
class InternalSideIndicator;
class Marker;
class Material;
class Transfer;
class XFEMInterface;
class SideUserObject;
class NodalUserObject;
class ElementUserObject;
class InternalSideUserObject;
class InterfaceUserObject;
class GeneralUserObject;
class Function;
class Distribution;
class Sampler;
class KernelBase;
class IntegratedBCBase;
class LineSearch;
class UserObject;
class AutomaticMortarGeneration;

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

template <>
InputParameters validParams<FEProblemBase>();

enum class MooseNonlinearConvergenceReason
{
  ITERATING = 0,
  CONVERGED_FNORM_ABS = 2,
  CONVERGED_FNORM_RELATIVE = 3,
  CONVERGED_SNORM_RELATIVE = 4,
  DIVERGED_FUNCTION_COUNT = -2,
  DIVERGED_FNORM_NAN = -4,
  DIVERGED_LINE_SEARCH = -6
};

// The idea with these enums is to abstract the reasons for
// convergence/divergence, i.e. they could be used with linear algebra
// packages other than PETSc.  They were directly inspired by PETSc,
// though.  This enum could also be combined with the
// MooseNonlinearConvergenceReason enum but there might be some
// confusion (?)
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:
  FEProblemBase(const InputParameters & parameters);
  virtual ~FEProblemBase();

  virtual EquationSystems & es() override { return _eq; }
  virtual MooseMesh & mesh() override { return _mesh; }
  virtual const MooseMesh & mesh() const override { return _mesh; }

  virtual Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) override;
  virtual 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() { return _coupling; }

  void setCouplingMatrix(std::unique_ptr<CouplingMatrix> cm);

  // DEPRECATED METHOD
  void setCouplingMatrix(CouplingMatrix * cm);

  const CouplingMatrix * couplingMatrix() { return _cm.get(); }

  void setNonlocalCouplingMatrix();

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

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

  virtual MooseNonlinearConvergenceReason
  checkNonlinearConvergence(std::string & msg,
                            const PetscInt it,
                            const Real xnorm,
                            const Real snorm,
                            const Real fnorm,
                            const Real rtol,
                            const Real stol,
                            const Real abstol,
                            const PetscInt nfuncs,
                            const PetscInt max_funcs,
                            const PetscBool force_iteration,
                            const Real initial_residual_before_preset_bcs,
                            const Real div_threshold);

  virtual MooseLinearConvergenceReason checkLinearConvergence(std::string & msg,
                                                              const PetscInt n,
                                                              const Real rnorm,
                                                              const Real rtol,
                                                              const Real atol,
                                                              const Real dtol,
                                                              const PetscInt maxits);

  virtual bool hasVariable(const std::string & var_name) const override;
  virtual MooseVariableFEBase & getVariable(
      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) override;
  virtual MooseVariable & getStandardVariable(THREAD_ID tid, const std::string & var_name) override;
  virtual VectorMooseVariable & getVectorVariable(THREAD_ID tid,
                                                  const std::string & var_name) override;

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

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

  virtual void clearActiveElementalMooseVariables(THREAD_ID tid) override;

  virtual void clearActiveFEVariableCoupleableMatrixTags(THREAD_ID tid) override;

  virtual void clearActiveFEVariableCoupleableVectorTags(THREAD_ID tid) override;

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

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

  virtual void clearActiveScalarVariableCoupleableMatrixTags(THREAD_ID tid) override;

  virtual void clearActiveScalarVariableCoupleableVectorTags(THREAD_ID tid) override;

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

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

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

  virtual void clearActiveMaterialProperties(THREAD_ID tid) override;

  virtual void createQRules(QuadratureType type,
                            Order order,
                            Order volume_order = INVALID_ORDER,
                            Order face_order = INVALID_ORDER);

  unsigned int getMaxQps() const;

  unsigned int getMaxShapeFunctions() const;

  Order getMaxScalarOrder() const;

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

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

  Assembly & assembly(THREAD_ID tid) override
  {
    mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
    return *_assembly[tid];
  }
  const Assembly & assembly(THREAD_ID tid) const override
  {
    mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
    return *_assembly[tid];
  }

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

  virtual void addExtraVectors();
  virtual void initialSetup();
  virtual void timestepSetup();

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

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

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

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

  virtual void newAssemblyArray(NonlinearSystemBase & nl);
  virtual void initNullSpaceVectors(const InputParameters & parameters, NonlinearSystemBase & nl);

  virtual void init() override;
  virtual void solve() override;

  const ConstElemRange & getEvaluableElementRange();

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

  virtual bool hasException() { return _has_exception; }

  virtual void checkExceptionAndStopSolve();

  virtual bool converged() override;
  virtual unsigned int nNonlinearIterations() const override;
  virtual unsigned int nLinearIterations() const override;
  virtual Real finalNonlinearResidual() const override;
  virtual bool computingInitialResidual() const override;

  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; }

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

  virtual void outputStep(ExecFlagType type);

  virtual void postExecute();


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

  void forceOutput();

  void initPetscOutput();

#ifdef LIBMESH_HAVE_PETSC

  Moose::PetscSupport::PetscOptions & getPetscOptions() { return _petsc_options; }
#endif // LIBMESH_HAVE_PETSC

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

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

  virtual void lineSearch();

  std::shared_ptr<LineSearch> getLineSearch() { return _line_search; }

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

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

  // NL /////
  NonlinearSystemBase & getNonlinearSystemBase() { return *_nl; }
  const NonlinearSystemBase & getNonlinearSystemBase() const { return *_nl; }

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

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

  virtual NonlinearSystem & getNonlinearSystem();

  virtual void addVariable(const std::string & var_name,
                           const FEType & type,
                           Real scale_factor,
                           const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void addScalarVariable(const std::string & var_name,
                                 Order order,
                                 Real scale_factor = 1.,
                                 const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void
  addKernel(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_name,
                              const FEType & type,
                              const std::set<SubdomainID> * const active_subdomains = NULL);
  virtual void addAuxScalarVariable(const std::string & var_name,
                                    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);

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

  void projectSolution();

  // Materials /////
  virtual void addMaterial(const std::string & kernel_name,
                           const std::string & name,
                           InputParameters parameters);
  virtual void addADResidualMaterial(const std::string & kernel_name,
                                     const std::string & name,
                                     InputParameters parameters);
  virtual void addADJacobianMaterial(const std::string & kernel_name,
                                     const std::string & name,
                                     InputParameters parameters);
  virtual void addMaterialHelper(std::vector<MaterialWarehouse *> warehouse,
                                 const std::string & kernel_name,
                                 const std::string & name,
                                 InputParameters parameters);

  virtual void prepareMaterials(SubdomainID blk_id, THREAD_ID tid);

  virtual void reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful = true);
  virtual void reinitMaterialsFace(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful = true);
  virtual void
  reinitMaterialsNeighbor(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful = true);
  virtual void
  reinitMaterialsBoundary(BoundaryID boundary_id, THREAD_ID tid, bool swap_stateful = true);
  /*
   * Swap back underlying data storing stateful material properties
   */
  virtual void swapBackMaterials(THREAD_ID tid);
  virtual void swapBackMaterialsFace(THREAD_ID tid);
  virtual void swapBackMaterialsNeighbor(THREAD_ID tid);

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

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

  void initPostprocessorData(const std::string & name);

  void initVectorPostprocessorData(const std::string & name);

  // UserObjects /////
  virtual void
  addUserObject(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<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;

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

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

  PostprocessorValue & getPostprocessorValue(const PostprocessorName & name);

  PostprocessorValue & getPostprocessorValueOld(const std::string & name);

  PostprocessorValue & getPostprocessorValueOlder(const std::string & name);


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

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

  VectorPostprocessorValue & getVectorPostprocessorValue(const VectorPostprocessorName & name,
                                                         const std::string & vector_name);

  VectorPostprocessorValue & getVectorPostprocessorValueOld(const std::string & name,
                                                            const std::string & vector_name);

  VectorPostprocessorValue & getVectorPostprocessorValue(const VectorPostprocessorName & name,
                                                         const std::string & vector_name,
                                                         bool needs_broadcast);

  VectorPostprocessorValue & getVectorPostprocessorValueOld(const std::string & name,
                                                            const std::string & vector_name,
                                                            bool needs_broadcast);

  ScatterVectorPostprocessorValue &
  getScatterVectorPostprocessorValue(const VectorPostprocessorName & vpp_name,
                                     const std::string & vector_name);

  ScatterVectorPostprocessorValue &
  getScatterVectorPostprocessorValueOld(const VectorPostprocessorName & vpp_name,
                                        const std::string & vector_name);

  VectorPostprocessorValue & declareVectorPostprocessorVector(const VectorPostprocessorName & name,
                                                              const std::string & vector_name,
                                                              bool contains_complete_history,
                                                              bool is_broadcast);

  bool vectorPostprocessorHasVectors(const std::string & vpp_name)
  {
    return _vpps_data.hasVectors(vpp_name);
  }

  const std::vector<std::pair<std::string, VectorPostprocessorData::VectorPostprocessorState>> &
  getVectorPostprocessorVectors(const std::string & vpp_name);

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

  bool hasDampers() { return _has_dampers; }

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

  // Markers //////
  virtual void
  addMarker(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,
                                                      MultiAppTransfer::DIRECTION direction) const;

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

  void execMultiAppTransfers(ExecFlagType type, MultiAppTransfer::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);

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

  void computeTransientImplicitResidual(Real time,
                                        const NumericVector<Number> & u,
                                        const NumericVector<Number> & udot,
                                        const NumericVector<Number> & udotdot,
                                        NumericVector<Number> & residual);

  void computeTransientImplicitJacobian(Real time,
                                        const NumericVector<Number> & u,
                                        const NumericVector<Number> & udot,
                                        const NumericVector<Number> & udotdot,
                                        Real duDotDu_shift,
                                        Real duDotDotDu_shift,
                                        SparseMatrix<Number> & jacobian);


  virtual Real computeResidualL2Norm();

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

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

  virtual void computeJacobianSys(NonlinearImplicitSystem & sys,
                                  const NumericVector<Number> & soln,
                                  SparseMatrix<Number> & jacobian);
  virtual void computeJacobian(const NumericVector<Number> & soln, SparseMatrix<Number> & jacobian);

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

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

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

  virtual void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks);

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

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

  virtual bool shouldUpdateSolution();

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

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

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

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

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

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

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

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

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

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

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

  // Displaced problem /////
  virtual void addDisplacedProblem(std::shared_ptr<DisplacedProblem> 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> & master_slave_boundary_pair,
                        const std::pair<SubdomainID, SubdomainID> & master_slave_subdomain_pair,
                        bool on_displaced,
                        bool periodic);

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

  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 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<InternalSideIndicator> & getInternalSideIndicatorWarehouse()
  {
    return _internal_side_indicators;
  }
  const MooseObjectWarehouse<Marker> & getMarkerWarehouse() { return _markers; }

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

  SolverParams & solverParams();

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

  void notifyWhenMeshChanges(MeshChangedInterface * mci);

  virtual void checkProblemIntegrity();

  void serializeSolution();

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

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

  void setKernelCoverageCheck(bool flag) { _kernel_coverage_check = flag; }

  void setMaterialCoverageCheck(bool flag) { _material_coverage_check = flag; }

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

  void parentOutputPositionChanged();


  bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid);
  bool needSubdomainMaterialOnSide(SubdomainID subdomain_id, 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 & getResidualMaterialsWarehouse() const { return _residual_materials; }
  const MaterialWarehouse & getJacobianMaterialsWarehouse() const { return _jacobian_materials; }
  const MaterialWarehouse & getDiscreteMaterialWarehouse() const { return _discrete_materials; }

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

  /*
   * Return a pointer to the MaterialData
   */
  std::shared_ptr<MaterialData> getMaterialData(Moose::MaterialDataType type, THREAD_ID tid = 0);

  bool errorOnJacobianNonzeroReallocation() const
  {
    return _error_on_jacobian_nonzero_reallocation;
  }

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

  bool ignoreZerosInJacobian() const { return _ignore_zeros_in_jacobian; }

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

  bool hasTimeIntegrator() const { return _has_time_integrator; }


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

  virtual void execute(const ExecFlagType & exec_type);

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

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

  virtual void computeAuxiliaryKernels(const ExecFlagType & type);

  void needsPreviousNewtonIteration(bool state);

  bool needsPreviousNewtonIteration() const;

  bool skipAdditionalRestartData() const { return _skip_additional_restart_data; }


  std::vector<Real> _real_zero;
  std::vector<VariableValue> _scalar_zero;
  std::vector<VariableValue> _zero;
  std::vector<MooseArray<DualReal>> _ad_zero;
  std::vector<VariableGradient> _grad_zero;
  std::vector<MooseArray<DualRealVectorValue>> _ad_grad_zero;
  std::vector<VariableSecond> _second_zero;
  std::vector<MooseArray<DualRealTensorValue>> _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; }

  const VectorPostprocessorData & getVectorPostprocessorData() const;

  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; }

  bool isSNESMFReuseBaseSetbyUser() { return _snesmf_reuse_base_set_by_user; }

  void usingADMatProps(bool using_ad_mat_props) { _using_ad_mat_props = using_ad_mat_props; }

  bool usingADMatProps() const { return _using_ad_mat_props; }

  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; }

protected:
  void createTagVectors();

  MooseMesh & _mesh;
  EquationSystems _eq;
  bool _initialized;

  std::set<TagID> _fe_vector_tags;

  std::set<TagID> _fe_matrix_tags;

  bool _solve;

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

  std::shared_ptr<NonlinearSystemBase> _nl;
  std::shared_ptr<AuxiliarySystem> _aux;

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

  // Dimension of the subspace spanned by the vectors with a given prefix
  std::map<std::string, unsigned int> _subspace_dim;

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

  MooseObjectWarehouse<Function> _functions;

  MooseObjectWarehouseBase<Distribution> _distributions;

  ExecuteMooseObjectWarehouse<Sampler> _samplers;

  MooseObjectWarehouse<KernelBase> _nonlocal_kernels;

  MooseObjectWarehouse<IntegratedBCBase> _nonlocal_integrated_bcs;

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

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

  std::vector<std::shared_ptr<MaterialData>> _material_data;
  std::vector<std::shared_ptr<MaterialData>> _bnd_material_data;
  std::vector<std::shared_ptr<MaterialData>> _neighbor_material_data;

  // Material Warehouses
  MaterialWarehouse _residual_materials; // Residual materials. This is the union of traditional
                                         // materials and the residual copy of an ADMaterial
  MaterialWarehouse _jacobian_materials; // Jacobian materials. This is the union of traditional
                                         // materials and the Jacobian copy of an ADMaterial
  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<InternalSideIndicator> _internal_side_indicators;

  // Marker Warehouse
  MooseObjectWarehouse<Marker> _markers;

  // postprocessors
  PostprocessorData _pps_data;

  // VectorPostprocessors
  VectorPostprocessorData _vpps_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;

  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<MeshChangedInterface *> _notify_when_mesh_changes;

  void meshChangedHelper(bool intermediate_change = false);

  bool duplicateVariableCheck(const std::string & var_name, const FEType & type, bool is_aux);

  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<Material>>> & materials_map);

  void checkCoordinateSystems();

  void reinitBecauseOfGhostingOrNewGeomObjects();

#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 _input_file_saved;

  bool _has_dampers;

  bool _has_constraints;

  bool _snesmf_reuse_base;

  bool _snesmf_reuse_base_set_by_user;

  bool _has_initialized_stateful;

  std::unique_ptr<Resurrector> _resurrector;

  bool _const_jacobian;

  bool _has_jacobian;

  bool _needs_old_newton_iter;

  bool _has_nonlocal_coupling;
  bool _calculate_jacobian_in_uo;

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

  SolverParams _solver_params;

  bool _kernel_coverage_check;

  bool _material_coverage_check;

  unsigned int _max_qps;

  unsigned int _max_shape_funcs;

  Order _max_scalar_order;

  bool _has_time_integrator;

  bool _has_exception;

  bool _parallel_barrier_messaging;

  std::string _exception_message;

  ExecFlagType _current_execute_on_flag;

  ExecuteMooseObjectWarehouse<Control> _control_warehouse;

#ifdef LIBMESH_HAVE_PETSC
  Moose::PetscSupport::PetscOptions _petsc_options;
#endif // LIBMESH_HAVE_PETSC

  std::shared_ptr<LineSearch> _line_search;

  std::unique_ptr<ConstElemRange> _evaluable_local_elem_range;

  bool _using_ad_mat_props;

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

  bool _error_on_jacobian_nonzero_reallocation;
  bool _ignore_zeros_in_jacobian;
  const bool _force_restart;
  const bool _skip_additional_restart_data;
  const bool _skip_nl_system_check;
  bool _fail_next_linear_convergence_check;

  bool _started_initial_setup;

  bool _has_internal_edge_residual_objects;

  const PerfID _initial_setup_timer;
  const PerfID _project_solution_timer;
  const PerfID _compute_indicators_timer;
  const PerfID _compute_markers_timer;
  const PerfID _compute_user_objects_timer;
  const PerfID _execute_controls_timer;
  const PerfID _execute_samplers_timer;
  const PerfID _update_active_objects_timer;
  const PerfID _reinit_because_of_ghosting_or_new_geom_objects_timer;
  const PerfID _exec_multi_app_transfers_timer;
  const PerfID _init_timer;
  const PerfID _eq_init_timer;
  const PerfID _solve_timer;
  const PerfID _check_exception_and_stop_solve_timer;
  const PerfID _advance_state_timer;
  const PerfID _restore_solutions_timer;
  const PerfID _save_old_solutions_timer;
  const PerfID _restore_old_solutions_timer;
  const PerfID _output_step_timer;
  const PerfID _on_timestep_begin_timer;
  const PerfID _compute_residual_l2_norm_timer;
  const PerfID _compute_residual_sys_timer;
  const PerfID _compute_residual_internal_timer;
  const PerfID _compute_residual_type_timer;
  const PerfID _compute_transient_implicit_residual_timer;
  const PerfID _compute_residual_tags_timer;
  const PerfID _compute_jacobian_internal_timer;
  const PerfID _compute_jacobian_tags_timer;
  const PerfID _compute_jacobian_blocks_timer;
  const PerfID _compute_bounds_timer;
  const PerfID _compute_post_check_timer;
  const PerfID _compute_damping_timer;
  const PerfID _possibly_rebuild_geom_search_patches_timer;
  const PerfID _initial_adapt_mesh_timer;
  const PerfID _adapt_mesh_timer;
  const PerfID _update_mesh_xfem_timer;
  const PerfID _mesh_changed_timer;
  const PerfID _mesh_changed_helper_timer;
  const PerfID _check_problem_integrity_timer;
  const PerfID _serialize_solution_timer;
  const PerfID _check_nonlinear_convergence_timer;
  const PerfID _check_linear_convergence_timer;
  const PerfID _update_geometric_search_timer;
  const PerfID _exec_multi_apps_timer;
  const PerfID _backup_multi_apps_timer;

  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 RestartableDataIO;
  friend class Restartable;
  friend class DisplacedProblem;
};

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