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

#include "SystemBase.h"
#include "ConstraintWarehouse.h"
#include "MooseObjectWarehouse.h"
#include "MooseObjectTagWarehouse.h"
#include "PerfGraphInterface.h"
#include "ComputeMortarFunctor.h"
#include "MooseHashing.h"

#include "libmesh/transient_system.h"
#include "libmesh/nonlinear_implicit_system.h"

// Forward declarations
class FEProblemBase;
class MoosePreconditioner;
class JacobianBlock;
class TimeIntegrator;
class Predictor;
class ElementDamper;
class NodalDamper;
class GeneralDamper;
class GeometricSearchData;
class IntegratedBCBase;
class NodalBCBase;
class PresetNodalBC;
template <ComputeStage>
class ADPresetNodalBC;
class DGKernelBase;
class InterfaceKernelBase;
class ScalarKernel;
class DiracKernel;
class NodalKernel;
class Split;
class KernelBase;
class BoundaryCondition;

// libMesh forward declarations
namespace libMesh
{
template <typename T>
class NumericVector;
template <typename T>
class SparseMatrix;
} // namespace libMesh

class NonlinearSystemBase : public SystemBase,
                            public ConsoleStreamInterface,
                            public PerfGraphInterface
{
public:
  NonlinearSystemBase(FEProblemBase & problem, System & sys, const std::string & name);
  virtual ~NonlinearSystemBase();

  virtual void init() override;

  void turnOffJacobian();

  virtual void addExtraVectors() override;
  virtual void solve() override = 0;
  virtual void restoreSolutions() override;

  virtual void stopSolve() = 0;

  virtual NonlinearSolver<Number> * nonlinearSolver() = 0;

  virtual unsigned int getCurrentNonlinearIterationNumber() = 0;

  virtual bool computingInitialResidual() { return _computing_initial_residual; }

  // Setup Functions ////
  virtual void initialSetup();
  virtual void timestepSetup();

  virtual void setupFiniteDifferencedPreconditioner() = 0;
  void setupFieldDecomposition();

  bool haveFiniteDifferencedPreconditioner() const
  {
    return _use_finite_differenced_preconditioner;
  }
  bool haveFieldSplitPreconditioner() const { return _use_field_split_preconditioner; }

  virtual bool converged() = 0;

  void addTimeIntegrator(const std::string & type,
                         const std::string & name,
                         InputParameters parameters) override;
  using SystemBase::addTimeIntegrator;

  void addDotVectors();

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

  void addScalarKernel(const std::string & kernel_name,
                       const std::string & name,
                       InputParameters parameters);

  void addBoundaryCondition(const std::string & bc_name,
                            const std::string & name,
                            InputParameters parameters);

  void
  addConstraint(const std::string & c_name, const std::string & name, InputParameters parameters);

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

  void
  addDGKernel(std::string dg_kernel_name, const std::string & name, InputParameters parameters);

  void addInterfaceKernel(std::string interface_kernel_name,
                          const std::string & name,
                          InputParameters parameters);

  void
  addDamper(const std::string & damper_name, const std::string & name, InputParameters parameters);

  void
  addSplit(const std::string & split_name, const std::string & name, InputParameters parameters);

  std::shared_ptr<Split> getSplit(const std::string & name);

  void zeroVectorForResidual(const std::string & vector_name);

  void setInitialSolution();

  void setConstraintSlaveValues(NumericVector<Number> & solution, bool displaced);

  void constraintResiduals(NumericVector<Number> & residual, bool displaced);

  void computeResidualTag(NumericVector<Number> & residual, TagID tag_id);

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

  void computeResidual(NumericVector<Number> & residual, TagID tag_id);

  void
  findImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data,
                                       std::map<dof_id_type, std::vector<dof_id_type>> & graph);

  void addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data);

  void constraintJacobians(bool displaced);

  void setVariableGlobalDoFs(const std::string & var_name);
  const std::vector<dof_id_type> & getVariableGlobalDoFs() { return _var_all_dof_indices; }

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

  void computeJacobian(SparseMatrix<Number> & jacobian, const std::set<TagID> & tags);

  void computeJacobian(SparseMatrix<Number> & jacobian);

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

  void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks, const std::set<TagID> & tags);

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

  void computeTimeDerivatives();

  void onTimestepBegin();

  virtual void subdomainSetup(SubdomainID subdomain, THREAD_ID tid);

  virtual void setSolution(const NumericVector<Number> & soln);

  void updateActive(THREAD_ID tid);

  virtual void setSolutionUDot(const NumericVector<Number> & udot);

  virtual void setSolutionUDotDot(const NumericVector<Number> & udotdot);

  NumericVector<Number> * solutionUDot() override { return _u_dot; }
  NumericVector<Number> * solutionUDotDot() override { return _u_dotdot; }
  NumericVector<Number> * solutionUDotOld() override { return _u_dot_old; }
  NumericVector<Number> * solutionUDotDotOld() override { return _u_dotdot_old; }
  const NumericVector<Number> * solutionUDot() const override { return _u_dot; }
  const NumericVector<Number> * solutionUDotDot() const override { return _u_dotdot; }
  const NumericVector<Number> * solutionUDotOld() const override { return _u_dot_old; }
  const NumericVector<Number> * solutionUDotDotOld() const override { return _u_dotdot_old; }

  NumericVector<Number> & getResidualTimeVector();

  NumericVector<Number> & getResidualNonTimeVector();

  NumericVector<Number> & residualVector(TagID tag);

  const NumericVector<Number> * const & currentSolution() const override
  {
    return _current_solution;
  }

  virtual void serializeSolution();
  virtual NumericVector<Number> & serializedSolution() override;

  virtual NumericVector<Number> & residualCopy() override;
  virtual NumericVector<Number> & residualGhosted() override;

  virtual NumericVector<Number> & RHS() = 0;

  virtual void augmentSparsity(SparsityPattern::Graph & sparsity,
                               std::vector<dof_id_type> & n_nz,
                               std::vector<dof_id_type> & n_oz) override;

  void setPreconditioner(std::shared_ptr<MoosePreconditioner> pc);
  MoosePreconditioner const * getPreconditioner() const;

  void useFiniteDifferencedPreconditioner(bool use = true)
  {
    _use_finite_differenced_preconditioner = use;
  }

  void setDecomposition(const std::vector<std::string> & decomposition);

  void useFieldSplitPreconditioner(bool use = true) { _use_field_split_preconditioner = use; }

  void addImplicitGeometricCouplingEntriesToJacobian(bool add = true)
  {
    _add_implicit_geometric_coupling_entries_to_jacobian = add;
  }

  void assembleConstraintsSeparately(bool separately = true)
  {
    _assemble_constraints_separately = separately;
  }

  void setupDampers();
  void reinitIncrementAtQpsForDampers(THREAD_ID tid, const std::set<MooseVariable *> & damped_vars);

  void reinitIncrementAtNodeForDampers(THREAD_ID tid,
                                       const std::set<MooseVariable *> & damped_vars);

  void checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const;
  bool containsTimeKernel();

  unsigned int nNonlinearIterations() const { return _n_iters; }

  unsigned int nLinearIterations() const { return _n_linear_iters; }

  unsigned int nResidualEvaluations() const { return _n_residual_evaluations; }

  Real finalNonlinearResidual() const { return _final_residual; }

  Real nonlinearNorm() const { return _last_nl_rnorm; }

  void printAllVariableNorms(bool state) { _print_all_var_norms = state; }

  void debuggingResiduals(bool state) { _debugging_residuals = state; }

  unsigned int _num_residual_evaluations;

  void setPredictor(std::shared_ptr<Predictor> predictor);
  Predictor * getPredictor() { return _predictor.get(); }

  void setPCSide(MooseEnum pcs);

  Moose::PCSideType getPCSide() { return _pc_side; }

  void setMooseKSPNormType(MooseEnum kspnorm);

  Moose::MooseKSPNormType getMooseKSPNormType() { return _ksp_norm; }

  bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;

  bool needSubdomainMaterialOnSide(SubdomainID subdomain_id, THREAD_ID tid) const;

  bool doingDG() const;


  MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() { return _kernels; }
  MooseObjectTagWarehouse<KernelBase> & getADJacobianKernelWarehouse()
  {
    return _ad_jacobian_kernels;
  }
  MooseObjectTagWarehouse<DGKernelBase> & getDGKernelWarehouse() { return _dg_kernels; }
  MooseObjectTagWarehouse<InterfaceKernelBase> & getInterfaceKernelWarehouse()
  {
    return _interface_kernels;
  }
  MooseObjectTagWarehouse<DiracKernel> & getDiracKernelWarehouse() { return _dirac_kernels; }
  MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() { return _integrated_bcs; }
  const MooseObjectWarehouse<ElementDamper> & getElementDamperWarehouse() const
  {
    return _element_dampers;
  }
  const MooseObjectWarehouse<NodalDamper> & getNodalDamperWarehouse() const
  {
    return _nodal_dampers;
  }
  const ConstraintWarehouse & getConstraintWarehouse() const { return _constraints; };

  bool hasSaveIn() const { return _has_save_in || _has_nodalbc_save_in; }

  bool hasDiagSaveIn() const { return _has_diag_save_in || _has_nodalbc_diag_save_in; }

  NumericVector<Number> & solution() override { return *_sys.solution; }
  const NumericVector<Number> & solution() const override { return *_sys.solution; }

  virtual System & system() override { return _sys; }
  virtual const System & system() const override { return _sys; }

  NumericVector<Number> * solutionPreviousNewton() override { return _solution_previous_nl; }
  const NumericVector<Number> * solutionPreviousNewton() const override
  {
    return _solution_previous_nl;
  }

  virtual void setSolutionUDotOld(const NumericVector<Number> & u_dot_old);

  virtual void setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old);

  virtual void setPreviousNewtonSolution(const NumericVector<Number> & soln);

  virtual TagID timeVectorTag() override { return _Re_time_tag; }

  virtual TagID nonTimeVectorTag() override { return _Re_non_time_tag; }

  virtual TagID residualVectorTag() override { return _Re_tag; }

  virtual TagID systemMatrixTag() override { return _Ke_system_tag; }

public:
  FEProblemBase & _fe_problem;
  System & _sys;
  // FIXME: make these protected and create getters/setters
  Real _last_rnorm;
  Real _last_nl_rnorm;
  Real _l_abs_step_tol;
  Real _initial_residual_before_preset_bcs;
  Real _initial_residual_after_preset_bcs;
  std::vector<unsigned int> _current_l_its;
  unsigned int _current_nl_its;
  bool _compute_initial_residual_before_preset_bcs;

protected:
  void computeResidualInternal(const std::set<TagID> & tags);

  void computeNodalBCs(NumericVector<Number> & residual);

  void computeNodalBCs(NumericVector<Number> & residual, const std::set<TagID> & tags);

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

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

  void computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian);

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

  void enforceNodalConstraintsResidual(NumericVector<Number> & residual);
  void enforceNodalConstraintsJacobian();

  void mortarResidualConstraints(bool displaced);

  void mortarJacobianConstraints(bool displaced);

protected:
  const NumericVector<Number> * _current_solution;
  NumericVector<Number> * _residual_ghosted;

  NumericVector<Number> & _serialized_solution;

  NumericVector<Number> * _solution_previous_nl;

  NumericVector<Number> & _residual_copy;

  NumericVector<Number> * _u_dot;
  NumericVector<Number> * _u_dotdot;

  NumericVector<Number> * _u_dot_old;
  NumericVector<Number> * _u_dotdot_old;

  Number _du_dot_du;
  Number _du_dotdot_du;

  TagID _Re_time_tag;

  std::set<TagID> _nl_vector_tags;

  std::set<TagID> _nl_matrix_tags;

  NumericVector<Number> * _Re_time;

  TagID _Re_non_time_tag;
  NumericVector<Number> * _Re_non_time;

  TagID _Re_tag;

  TagID _Ke_non_time_tag;

  TagID _Ke_system_tag;

  MooseObjectTagWarehouse<KernelBase> _kernels;
  MooseObjectTagWarehouse<KernelBase> _ad_jacobian_kernels;
  MooseObjectTagWarehouse<ScalarKernel> _scalar_kernels;
  MooseObjectTagWarehouse<DGKernelBase> _dg_kernels;
  MooseObjectTagWarehouse<InterfaceKernelBase> _interface_kernels;


  MooseObjectTagWarehouse<IntegratedBCBase> _integrated_bcs;
  MooseObjectTagWarehouse<NodalBCBase> _nodal_bcs;
  MooseObjectWarehouse<PresetNodalBC> _preset_nodal_bcs;
  MooseObjectWarehouse<ADPresetNodalBC<RESIDUAL>> _ad_preset_nodal_bcs;

  MooseObjectTagWarehouse<DiracKernel> _dirac_kernels;

  MooseObjectWarehouse<ElementDamper> _element_dampers;

  MooseObjectWarehouse<NodalDamper> _nodal_dampers;

  MooseObjectWarehouse<GeneralDamper> _general_dampers;

  MooseObjectTagWarehouse<NodalKernel> _nodal_kernels;

  MooseObjectWarehouseBase<Split> _splits; // use base b/c there are no setup methods

  ConstraintWarehouse _constraints;

  NumericVector<Number> * _increment_vec;
  std::shared_ptr<MoosePreconditioner> _preconditioner;
  Moose::PCSideType _pc_side;
  Moose::MooseKSPNormType _ksp_norm;

  bool _use_finite_differenced_preconditioner;
#ifdef LIBMESH_HAVE_PETSC
  MatFDColoring _fdcoloring;
#endif
  bool _have_decomposition;
  std::string _decomposition_split;
  bool _use_field_split_preconditioner;

  bool _add_implicit_geometric_coupling_entries_to_jacobian;

  bool _assemble_constraints_separately;

  bool _need_serialized_solution;

  bool _need_residual_copy;
  bool _need_residual_ghosted;
  bool _debugging_residuals;

  bool _doing_dg;

  std::vector<std::string> _vecs_to_zero_for_residual;

  unsigned int _n_iters;
  unsigned int _n_linear_iters;

  unsigned int _n_residual_evaluations;

  Real _final_residual;

  std::shared_ptr<Predictor> _predictor;

  bool _computing_initial_residual;

  bool _print_all_var_norms;

  bool _has_save_in;

  bool _has_diag_save_in;

  bool _has_nodalbc_save_in;

  bool _has_nodalbc_diag_save_in;

  void getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs);

  std::vector<dof_id_type> _var_all_dof_indices;

  PerfID _compute_residual_tags_timer;
  PerfID _compute_residual_internal_timer;
  PerfID _kernels_timer;
  PerfID _scalar_kernels_timer;
  PerfID _nodal_kernels_timer;
  PerfID _nodal_kernel_bcs_timer;
  PerfID _nodal_bcs_timer;
  PerfID _compute_jacobian_tags_timer;
  PerfID _compute_jacobian_blocks_timer;
  PerfID _compute_dampers_timer;
  PerfID _compute_dirac_timer;

private:
  std::unordered_map<std::pair<BoundaryID, BoundaryID>,
                     ComputeMortarFunctor<ComputeStage::RESIDUAL>>
      _undisplaced_mortar_residual_functors;

  std::unordered_map<std::pair<BoundaryID, BoundaryID>,
                     ComputeMortarFunctor<ComputeStage::JACOBIAN>>
      _undisplaced_mortar_jacobian_functors;

  std::unordered_map<std::pair<BoundaryID, BoundaryID>,
                     ComputeMortarFunctor<ComputeStage::RESIDUAL>>
      _displaced_mortar_residual_functors;

  std::unordered_map<std::pair<BoundaryID, BoundaryID>,
                     ComputeMortarFunctor<ComputeStage::JACOBIAN>>
      _displaced_mortar_jacobian_functors;
};