SystemBase.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include <map>
#include <set>
#include <vector>
 
#include "DataIO.h"
#include "MooseTypes.h"
#include "VariableWarehouse.h"
#include "InputParameters.h"
#include "MooseVariableBase.h"
#include "ConsoleStreamInterface.h"
// libMesh
#include "libmesh/exodusII_io.h"
#include "libmesh/parallel_object.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/sparse_matrix.h"
 
// Forward declarations
class Factory;
class MooseApp;
class MooseVariableFieldBase;
template <typename>
class MooseVariableFE;
typedef MooseVariableFE<Real> MooseVariable;
typedef MooseVariableFE<VectorValue<Real>> VectorMooseVariable;
class MooseMesh;
class SubProblem;
class SystemBase;
class TimeIntegrator;
class InputParameters;
class FEProblemBase;
 
// libMesh forward declarations
namespace libMesh
{
class System;
class DofMap;
class FEType;
}
 
void extraSendList(std::vector<dof_id_type> & send_list, void * context);
 
void extraSparsity(libMesh::SparsityPattern::Graph & sparsity,
                   std::vector<dof_id_type> & n_nz,
                   std::vector<dof_id_type> & n_oz,
                   void * context);
 
struct VarCopyInfo
{
  VarCopyInfo(const std::string & dest_name,
              const std::string & source_name,
              const std::string & timestep)
    : _dest_name(dest_name), _source_name(source_name), _timestep(timestep)
  {
  }
 
  std::string _dest_name;
  std::string _source_name;
  std::string _timestep;
};
 
class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
 
{
public:
  SystemBase(SubProblem & subproblem,
             FEProblemBase & fe_problem,
             const std::string & name,
             Moose::VarKindType var_kind);
  virtual ~SystemBase() {}
 
  unsigned int number() const;
  MooseMesh & mesh() { return _mesh; }
  const MooseMesh & mesh() const { return _mesh; }
  SubProblem & subproblem() { return _subproblem; }
  const SubProblem & subproblem() const { return _subproblem; }
  FEProblemBase & feProblem() { return _fe_problem; }
  const FEProblemBase & feProblem() const { return _fe_problem; }
 
  void applyScalingFactors(const std::vector<Real> & inverse_scaling_factors);
 
  bool computingScalingJacobian() const;
 
  bool automaticScaling() const { return _automatic_scaling; }
 
  void automaticScaling(bool automatic_scaling) { _automatic_scaling = automatic_scaling; }
 
  void setVerboseFlag(const bool & verbose) { _verbose = verbose; }
 
  virtual libMesh::DofMap & dofMap();
 
  virtual const libMesh::DofMap & dofMap() const;
 
  virtual libMesh::System & system() = 0;
  virtual const libMesh::System & system() const = 0;
 
  virtual void preInit() {}
 
  /*
   * This is called after the libMesh system has been init'd. This can be used to initialize MOOSE
   * system data that relies on the libMesh system data already being initialized
   */
  virtual void postInit() {}
 
  virtual void reinit() {}
 
  virtual void initializeObjects() {}
 
  void update();
 
  virtual void solve();
 
  void copyOldSolutions();
  void copyPreviousSolutions(const Moose::SolutionIterationType iteration_type);
 
  virtual void restoreSolutions();
 
  virtual const NumericVector<Number> * const & currentSolution() const = 0;
 
  NumericVector<Number> & solution() { return solutionState(0); }
  NumericVector<Number> & solutionOld() { return solutionState(1); }
  NumericVector<Number> & solutionOlder() { return solutionState(2); }
  const NumericVector<Number> & solution() const { return solutionState(0); }
  const NumericVector<Number> & solutionOld() const { return solutionState(1); }
  const NumericVector<Number> & solutionOlder() const { return solutionState(2); }
 
  virtual const NumericVector<Number> * solutionPreviousNewton() const;
  virtual NumericVector<Number> * solutionPreviousNewton();
 
  virtual void initSolutionState();
 
  const std::vector<NumericVector<Number> *> &
  getSolutionStates(const Moose::SolutionIterationType iteration_type) const;
 
  std::size_t getNumSolutionStates(const Moose::SolutionIterationType iteration_type) const;
 
  virtual NumericVector<Number> &
  solutionState(const unsigned int state,
                Moose::SolutionIterationType iteration_type = Moose::SolutionIterationType::Time);
 
  virtual const NumericVector<Number> & solutionState(
      const unsigned int state,
      Moose::SolutionIterationType iteration_type = Moose::SolutionIterationType::Time) const;
 
  libMesh::ParallelType
  solutionStateParallelType(const unsigned int state,
                            const Moose::SolutionIterationType iteration_type) const;
 
  virtual void needSolutionState(
      const unsigned int state,
      Moose::SolutionIterationType iteration_type = Moose::SolutionIterationType::Time,
      libMesh::ParallelType parallel_type = GHOSTED);
 
  virtual bool hasSolutionState(
      const unsigned int state,
      Moose::SolutionIterationType iteration_type = Moose::SolutionIterationType::Time) const;
 
  virtual void addDotVectors();
 
  virtual std::vector<Number> & duDotDus() { return _du_dot_du; }
  virtual Number & duDotDotDu() { return _du_dotdot_du; }
  virtual const Number & duDotDu(unsigned int var_num = 0) const;
  virtual const Number & duDotDotDu() const { return _du_dotdot_du; }
 
  virtual NumericVector<Number> * solutionUDot() { return _u_dot; }
  virtual NumericVector<Number> * solutionUDotDot() { return _u_dotdot; }
  virtual NumericVector<Number> * solutionUDotOld() { return _u_dot_old; }
  virtual NumericVector<Number> * solutionUDotDotOld() { return _u_dotdot_old; }
  virtual const NumericVector<Number> * solutionUDot() const { return _u_dot; }
  virtual const NumericVector<Number> * solutionUDotDot() const { return _u_dotdot; }
  virtual const NumericVector<Number> * solutionUDotOld() const { return _u_dot_old; }
  virtual const NumericVector<Number> * solutionUDotDotOld() const { return _u_dotdot_old; }
 
  virtual void saveOldSolutions();
  virtual void restoreOldSolutions();
 
  bool hasVector(const std::string & tag_name) const;
 
  virtual bool hasVector(TagID tag_id) const
  {
    return tag_id < _tagged_vectors.size() && _tagged_vectors[tag_id];
  }
 
  virtual TagID timeVectorTag() const { mooseError("Not implemented yet"); }
 
  virtual TagID systemMatrixTag() const { mooseError("Not implemented yet"); }
 
  /*
   * Return TagID for nontime
   */
  virtual TagID nonTimeVectorTag() const { mooseError("Not implemented yet"); }
 
  /*
   * Return TagID for nontime
   */
  virtual TagID residualVectorTag() const { mooseError("Not implemented yet"); }
 
  virtual std::set<TagID> defaultVectorTags() const
  {
    return {timeVectorTag(), nonTimeVectorTag(), residualVectorTag()};
  }
  virtual std::set<TagID> defaultMatrixTags() const { return {systemMatrixTag()}; }
 
  virtual NumericVector<Number> & getVector(const std::string & name);
  virtual const NumericVector<Number> & getVector(const std::string & name) const;
 
  virtual NumericVector<Number> & getVector(TagID tag);
  virtual const NumericVector<Number> & getVector(TagID tag) const;
 
  virtual void associateVectorToTag(NumericVector<Number> & vec, TagID tag);
 
  virtual void disassociateVectorFromTag(NumericVector<Number> & vec, TagID tag);
 
  virtual void disassociateVectorFromTag(TagID tag);
 
  virtual void disassociateDefaultVectorTags();
 
  virtual bool hasMatrix(TagID tag) const
  {
    return tag < _tagged_matrices.size() && _tagged_matrices[tag];
  }
 
  virtual libMesh::SparseMatrix<Number> & getMatrix(TagID tag);
 
  virtual const libMesh::SparseMatrix<Number> & getMatrix(TagID tag) const;
 
  virtual void activateAllMatrixTags();
 
  virtual bool matrixTagActive(TagID tag) const;
 
  virtual void deactivateAllMatrixTags();
 
  void closeTaggedMatrices(const std::set<TagID> & tags);
 
  void flushTaggedMatrices(const std::set<TagID> & tags);
 
  virtual void associateMatrixToTag(libMesh::SparseMatrix<Number> & matrix, TagID tag);
 
  virtual void disassociateMatrixFromTag(libMesh::SparseMatrix<Number> & matrix, TagID tag);
 
  virtual void disassociateMatrixFromTag(TagID tag);
 
  virtual void disassociateDefaultMatrixTags();
 
  virtual NumericVector<Number> & serializedSolution();
 
  virtual NumericVector<Number> & residualCopy()
  {
    mooseError("This system does not support getting a copy of the residual");
  }
  virtual NumericVector<Number> & residualGhosted()
  {
    mooseError("This system does not support getting a ghosted copy of the residual");
  }
 
  virtual void augmentSendList(std::vector<dof_id_type> & send_list);
 
  virtual void augmentSparsity(libMesh::SparsityPattern::Graph & sparsity,
                               std::vector<dof_id_type> & n_nz,
                               std::vector<dof_id_type> & n_oz) = 0;
 
  virtual void addVariable(const std::string & var_type,
                           const std::string & var_name,
                           InputParameters & parameters);
 
  virtual bool isArrayVariable(const std::string & var_name) const;
 
 
  virtual bool hasVariable(const std::string & var_name) const;
  virtual bool hasScalarVariable(const std::string & var_name) const;
 
  virtual bool isScalarVariable(unsigned int var_name) const;
 
  MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string & var_name) const;
 
  MooseVariableFieldBase & getVariable(THREAD_ID tid, unsigned int var_number) const;
 
  template <typename T>
  MooseVariableFE<T> & getFieldVariable(THREAD_ID tid, const std::string & var_name);
 
  template <typename T>
  MooseVariableField<T> & getActualFieldVariable(THREAD_ID tid, const std::string & var_name);
 
  template <typename T>
  MooseVariableFE<T> & getFieldVariable(THREAD_ID tid, unsigned int var_number);
 
  template <typename T>
  MooseVariableField<T> & getActualFieldVariable(THREAD_ID tid, unsigned int var_number);
 
  template <typename T>
  MooseVariableFV<T> & getFVVariable(THREAD_ID tid, const std::string & var_name);
 
  virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid,
                                                  const std::string & var_name) const;
 
  virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, unsigned int var_number) const;
 
  virtual const std::set<SubdomainID> * getVariableBlocks(unsigned int var_number);
 
  virtual unsigned int nVariables() const;
 
  unsigned int nFieldVariables() const;
 
  unsigned int nFVVariables() const;
 
  std::size_t getMaxVarNDofsPerElem() const { return _max_var_n_dofs_per_elem; }
 
  std::size_t getMaxVarNDofsPerNode() const { return _max_var_n_dofs_per_node; }
 
  void assignMaxVarNDofsPerElem(std::size_t max_dofs) { _max_var_n_dofs_per_elem = max_dofs; }
 
  void assignMaxVarNDofsPerNode(std::size_t max_dofs) { _max_var_n_dofs_per_node = max_dofs; }
 
  virtual void addVariableToZeroOnResidual(std::string var_name);
 
  virtual void addVariableToZeroOnJacobian(std::string var_name);
 
  virtual void zeroVariables(std::vector<std::string> & vars_to_be_zeroed);
 
  virtual void zeroVariablesForResidual();
 
  virtual void zeroVariablesForJacobian();
 
  virtual libMesh::Order getMinQuadratureOrder();
 
  virtual void prepare(THREAD_ID tid);
 
  virtual void prepareFace(THREAD_ID tid, bool resize_data);
 
  virtual void prepareNeighbor(THREAD_ID tid);
 
  virtual void prepareLowerD(THREAD_ID tid);
 
  virtual void reinitElem(const Elem * elem, THREAD_ID tid);
 
  virtual void reinitElemFace(const Elem * elem, unsigned int side, THREAD_ID tid);
 
  virtual void reinitNeighborFace(const Elem * elem, unsigned int side, THREAD_ID tid);
 
  virtual void reinitNeighbor(const Elem * elem, THREAD_ID tid);
 
  virtual void reinitLowerD(THREAD_ID tid);
 
  virtual void reinitNode(const Node * node, THREAD_ID tid);
 
  virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, THREAD_ID tid);
 
  virtual void reinitNodes(const std::vector<dof_id_type> & nodes, THREAD_ID tid);
 
  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, THREAD_ID tid);
 
  virtual void reinitScalars(THREAD_ID tid, bool reinit_for_derivative_reordering = false);
 
  virtual void addVariableToCopy(const std::string & dest_name,
                                 const std::string & source_name,
                                 const std::string & timestep);
 
  const std::vector<MooseVariableFieldBase *> & getVariables(THREAD_ID tid)
  {
    return _vars[tid].fieldVariables();
  }
 
  const VariableWarehouse & variableWarehouse(THREAD_ID tid = 0) const { return _vars[tid]; }
 
  const std::vector<MooseVariableScalar *> & getScalarVariables(THREAD_ID tid)
  {
    return _vars[tid].scalars();
  }
 
  const std::set<SubdomainID> & getSubdomainsForVar(unsigned int var_number) const
  {
    return _var_map.at(var_number);
  }
 
  const std::set<SubdomainID> & getSubdomainsForVar(const std::string & var_name) const;
 
  void removeVector(const std::string & name);
 
  NumericVector<Number> &
  addVector(const std::string & vector_name, const bool project, const libMesh::ParallelType type);
 
  NumericVector<Number> &
  addVector(TagID tag, const bool project, const libMesh::ParallelType type);
 
  void closeTaggedVector(const TagID tag);
  void closeTaggedVectors(const std::set<TagID> & tags);
 
  void zeroTaggedVector(const TagID tag);
  void zeroTaggedVectors(const std::set<TagID> & tags);
 
  void removeVector(TagID tag_id);
 
  void setVariableGlobalDoFs(const std::string & var_name);
 
  const std::vector<dof_id_type> & getVariableGlobalDoFs() { return _var_all_dof_indices; }
 
  libMesh::SparseMatrix<Number> & addMatrix(TagID tag);
 
  void removeMatrix(TagID tag);
 
  virtual const std::string & name() const;
 
  const std::vector<VariableName> & getVariableNames() const { return _vars[0].names(); }
 
  void getStandardFieldVariableNames(std::vector<VariableName> & std_field_variables) const;
 
  unsigned int getMaxVariableNumber() const { return _max_var_number; }
 
  virtual void computeVariables(const NumericVector<Number> & /*soln*/) {}
 
  void copyVars(libMesh::ExodusII_IO & io);
 
  virtual void copySolutionsBackwards();
 
  void addTimeIntegrator(const std::string & type,
                         const std::string & name,
                         InputParameters & parameters);
 
  bool hasVarCopy() const { return _var_to_copy.size() > 0; }
 
  void addScalingVector();
 
  bool solutionStatesInitialized() const { return _solution_states_initialized; }
 
  virtual void initialSetup();
  virtual void timestepSetup();
  virtual void customSetup(const ExecFlagType & exec_type);
  virtual void subdomainSetup();
  virtual void residualSetup();
  virtual void jacobianSetup();
 
  void clearAllDofIndices();
 
  void setActiveVariableCoupleableVectorTags(const std::set<TagID> & vtags, THREAD_ID tid);
 
  void setActiveScalarVariableCoupleableVectorTags(const std::set<TagID> & vtags, THREAD_ID tid);
 
  Moose::VarKindType varKind() const { return _var_kind; }
 
  virtual void compute(ExecFlagType type) = 0;
 
  void copyTimeIntegrators(const SystemBase & other_sys);
 
  const TimeIntegrator & getTimeIntegrator(const unsigned int var_num) const;
 
  const TimeIntegrator * queryTimeIntegrator(const unsigned int var_num) const;
 
  const std::vector<std::shared_ptr<TimeIntegrator>> & getTimeIntegrators();
 
  std::string prefix() const;
 
  void sizeVariableMatrixData();
 
  void skipNextSolutionToOldCopy() { _skip_next_solution_to_old_copy = true; }
 
protected:
  virtual NumericVector<Number> & solutionInternal() const = 0;
 
  SubProblem & _subproblem;
 
  FEProblemBase & _fe_problem;
 
  MooseApp & _app;
  Factory & _factory;
 
  MooseMesh & _mesh;
  std::string _name;
 
  std::vector<VariableWarehouse> _vars;
  std::map<unsigned int, std::set<SubdomainID>> _var_map;
  unsigned int _max_var_number;
 
  std::vector<std::string> _vars_to_be_zeroed_on_residual;
  std::vector<std::string> _vars_to_be_zeroed_on_jacobian;
 
  NumericVector<Number> * _u_dot;
  NumericVector<Number> * _u_dotdot;
 
  NumericVector<Number> * _u_dot_old;
  NumericVector<Number> * _u_dotdot_old;
 
  std::vector<Real> _du_dot_du;
  Real _du_dotdot_du;
 
  std::vector<NumericVector<Number> *> _tagged_vectors;
  std::vector<libMesh::SparseMatrix<Number> *> _tagged_matrices;
  std::unordered_map<TagID, libMesh::SparseMatrix<Number> *> _active_tagged_matrices;
  std::vector<bool> _matrix_tag_active_flags;
 
  // Used for saving old solutions so that they wont be accidentally changed
  NumericVector<Real> * _saved_old;
  NumericVector<Real> * _saved_older;
 
  // Used for saving old u_dot and u_dotdot so that they wont be accidentally changed
  NumericVector<Real> * _saved_dot_old;
  NumericVector<Real> * _saved_dotdot_old;
 
  Moose::VarKindType _var_kind;
 
  std::vector<VarCopyInfo> _var_to_copy;
 
  size_t _max_var_n_dofs_per_elem;
 
  size_t _max_var_n_dofs_per_node;
 
  std::vector<std::shared_ptr<TimeIntegrator>> _time_integrators;
 
  std::vector<std::vector<MooseVariableFieldBase *>> _numbered_vars;
 
  bool _automatic_scaling;
 
  bool _verbose;
 
  bool _solution_states_initialized;
 
  std::vector<dof_id_type> _var_all_dof_indices;
 
  std::unique_ptr<NumericVector<Number>> _serialized_solution;
 
private:
  TagName oldSolutionStateVectorName(const unsigned int,
                                     Moose::SolutionIterationType iteration_type) const;
 
  std::vector<NumericVector<Number> *> &
  getSolutionStates(const Moose::SolutionIterationType iteration_type)
  {
    return const_cast<std::vector<NumericVector<Number> *> &>(
        static_cast<const SystemBase *>(this)->getSolutionStates(iteration_type));
  }
 
  std::array<std::vector<NumericVector<Number> *>,
             static_cast<size_t>(Moose::SolutionIterationType::Count)>
      _solution_states;
  std::vector<NumericVector<Number> *> _saved_solution_states;
  bool _skip_next_solution_to_old_copy;
};
 
inline const std::vector<NumericVector<Number> *> &
SystemBase::getSolutionStates(const Moose::SolutionIterationType iteration_type) const
{
  const auto iteration_type_index = static_cast<std::size_t>(iteration_type);
  mooseAssert(iteration_type_index < static_cast<std::size_t>(Moose::SolutionIterationType::Count),
              "Invalid solution iteration type");
  mooseAssert(iteration_type_index < _solution_states.size(), "_solution_states sized incorrectly");
  return _solution_states[iteration_type_index];
}
 
inline std::size_t
SystemBase::getNumSolutionStates(const Moose::SolutionIterationType iteration_type) const
{
  return getSolutionStates(iteration_type).size();
}
 
inline bool
SystemBase::hasSolutionState(const unsigned int state,
                             const Moose::SolutionIterationType iteration_type) const
{
  return getNumSolutionStates(iteration_type) > state;
}
 
#define PARALLEL_TRY
 
#define PARALLEL_CATCH _fe_problem.checkExceptionAndStopSolve();