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

  virtual void copyOldSolutions();
  virtual void copyPreviousNonlinearSolutions();
  virtual void copyPreviousFixedPointSolutions();
  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();

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

  const std::vector<std::unique_ptr<NumericVector<Number>>> & gradientContainer() const
  {
    return _raw_grad_container;
  }

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

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;

  std::vector<std::unique_ptr<NumericVector<Number>>> _raw_grad_container;

private:
  TagName oldSolutionStateVectorName(const unsigned int,
                                     Moose::SolutionIterationType iteration_type) const;

  std::array<std::vector<NumericVector<Number> *>, 3> _solution_states;
  std::vector<NumericVector<Number> *> _saved_solution_states;
};

inline bool
SystemBase::hasSolutionState(const unsigned int state,
                             const Moose::SolutionIterationType iteration_type) const
{
  return _solution_states[static_cast<unsigned short>(iteration_type)].size() > state;
}

#define PARALLEL_TRY

#define PARALLEL_CATCH _fe_problem.checkExceptionAndStopSolve();