NonlinearSystemBase.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 "SolverSystem.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"
#include "libmesh/linear_solver.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 DirichletBCBase;
class ADDirichletBCBase;
class DGKernelBase;
class InterfaceKernelBase;
class ScalarKernelBase;
class DiracKernelBase;
class NodalKernelBase;
class Split;
class KernelBase;
class HDGKernel;
class BoundaryCondition;
class ResidualObject;
class PenetrationInfo;

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

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

  virtual void preInit() override;

  bool computedScalingJacobian() const { return _computed_scaling; }

  virtual void turnOffJacobian();

  virtual void solve() override = 0;

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

  virtual SNES getSNES() = 0;

  virtual unsigned int getCurrentNonlinearIterationNumber() = 0;

  bool computingPreSMOResidual() { return _computing_pre_smo_residual; }

  // Setup Functions ////
  virtual void initialSetup() override;
  virtual void timestepSetup() override;
  virtual void customSetup(const ExecFlagType & exec_type) override;
  virtual void residualSetup() override;
  virtual void jacobianSetup() override;

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

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

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

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

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

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

  MooseObjectWarehouseBase<Split> & getSplits() { return _splits; }

  bool shouldEvaluatePreSMOResidual() const;

  void setPreSMOResidual(bool use) { _use_pre_smo_residual = use; }

  const bool & usePreSMOResidual() const { return _use_pre_smo_residual; }

  Real referenceResidual() const;

  Real preSMOResidual() const;

  Real initialResidual() const;

  void setInitialResidual(Real r);

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

  void setInitialSolution();

  void setConstraintSecondaryValues(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 computeResidualAndJacobianTags(const std::set<TagID> & vector_tags,
                                      const std::set<TagID> & matrix_tags);

  void computeResidualAndJacobianInternal(const std::set<TagID> & vector_tags,
                                          const std::set<TagID> & matrix_tags);

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

  void addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data);

  void constraintJacobians(const SparseMatrix<Number> & jacobian_to_view, bool displaced);

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

  bool computeScaling();

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

  void computeJacobian(libMesh::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 onTimestepBegin();

  using SystemBase::subdomainSetup;
  virtual void subdomainSetup(SubdomainID subdomain, THREAD_ID tid);

  void overwriteNodeFace(NumericVector<Number> & soln);

  void updateActive(THREAD_ID tid);

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

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

  NumericVector<Number> & getResidualTimeVector();

  NumericVector<Number> & getResidualNonTimeVector();

  NumericVector<Number> & residualVector(TagID tag);

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

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

  virtual void augmentSparsity(libMesh::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;
  }

  virtual void attachPreconditioner(libMesh::Preconditioner<Number> * preconditioner) = 0;

  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;
  virtual bool containsTimeKernel() override;
  virtual std::vector<std::string> timeKernelVariableNames() override;

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

  bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;

  bool needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;

  bool needSubdomainMaterialOnSide(SubdomainID subdomain_id, THREAD_ID tid) const;

  bool doingDG() const;


  MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() { return _kernels; }
  const MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() const { return _kernels; }
  MooseObjectTagWarehouse<DGKernelBase> & getDGKernelWarehouse() { return _dg_kernels; }
  MooseObjectTagWarehouse<InterfaceKernelBase> & getInterfaceKernelWarehouse()
  {
    return _interface_kernels;
  }
  MooseObjectTagWarehouse<DiracKernelBase> & getDiracKernelWarehouse() { return _dirac_kernels; }
  MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() { return _integrated_bcs; }
  const MooseObjectTagWarehouse<ScalarKernelBase> & getScalarKernelWarehouse() const
  {
    return _scalar_kernels;
  }
  const MooseObjectTagWarehouse<NodalKernelBase> & getNodalKernelWarehouse() const
  {
    return _nodal_kernels;
  }
  MooseObjectTagWarehouse<HDGKernel> & getHDGKernelWarehouse() { return _hybridized_kernels; }
  const MooseObjectWarehouse<ElementDamper> & getElementDamperWarehouse() const
  {
    return _element_dampers;
  }
  const MooseObjectWarehouse<NodalDamper> & getNodalDamperWarehouse() const
  {
    return _nodal_dampers;
  }
  const ConstraintWarehouse & getConstraintWarehouse() const { return _constraints; }

  const MooseObjectTagWarehouse<NodalBCBase> & getNodalBCWarehouse() const { return _nodal_bcs; }

  const MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() const
  {
    return _integrated_bcs;
  }


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

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

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

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

  TagID timeVectorTag() const override { return _Re_time_tag; }
  TagID nonTimeVectorTag() const override { return _Re_non_time_tag; }
  TagID residualVectorTag() const override { return _Re_tag; }
  TagID systemMatrixTag() const override { return _Ke_system_tag; }

  virtual void residualAndJacobianTogether() = 0;

  bool computeScalingOnce() const { return _compute_scaling_once; }
  void computeScalingOnce(bool compute_scaling_once)
  {
    _compute_scaling_once = compute_scaling_once;
  }

  void autoScalingParam(Real resid_vs_jac_scaling_param)
  {
    _resid_vs_jac_scaling_param = resid_vs_jac_scaling_param;
  }

  void scalingGroupVariables(const std::vector<std::vector<std::string>> & scaling_group_variables)
  {
    _scaling_group_variables = scaling_group_variables;
  }

  void
  ignoreVariablesForAutoscaling(const std::vector<std::string> & ignore_variables_for_autoscaling)
  {
    _ignore_variables_for_autoscaling = ignore_variables_for_autoscaling;
  }

  bool offDiagonalsInAutoScaling() const { return _off_diagonals_in_auto_scaling; }
  void offDiagonalsInAutoScaling(bool off_diagonals_in_auto_scaling)
  {
    _off_diagonals_in_auto_scaling = off_diagonals_in_auto_scaling;
  }

  libMesh::System & _sys;
  // FIXME: make these protected and create getters/setters
  Real _last_nl_rnorm;
  std::vector<unsigned int> _current_l_its;
  unsigned int _current_nl_its;

  void setupDM();

  using SystemBase::reinitNodeFace;

  virtual void potentiallySetupFiniteDifferencing() {}

  void destroyColoring();

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

  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 mortarConstraints(Moose::ComputeType compute_type,
                         const std::set<TagID> & vector_tags,
                         const std::set<TagID> & matrix_tags);

  virtual void computeScalingJacobian() = 0;

  virtual void computeScalingResidual() = 0;

  void assembleScalingVector();

  virtual void postAddResidualObject(ResidualObject &) {}

  void reinitNodeFace(const Node & secondary_node,
                      const BoundaryID secondary_boundary,
                      const PenetrationInfo & info,
                      const bool displaced);

  bool preSolve();

  NumericVector<Number> * _residual_ghosted;

  std::unique_ptr<NumericVector<Number>> _residual_copy;

  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<HDGKernel> _hybridized_kernels;
  MooseObjectTagWarehouse<ScalarKernelBase> _scalar_kernels;
  MooseObjectTagWarehouse<DGKernelBase> _dg_kernels;
  MooseObjectTagWarehouse<InterfaceKernelBase> _interface_kernels;


  MooseObjectTagWarehouse<IntegratedBCBase> _integrated_bcs;
  MooseObjectTagWarehouse<NodalBCBase> _nodal_bcs;
  MooseObjectWarehouse<DirichletBCBase> _preset_nodal_bcs;
  MooseObjectWarehouse<ADDirichletBCBase> _ad_preset_nodal_bcs;

  MooseObjectTagWarehouse<DiracKernelBase> _dirac_kernels;

  MooseObjectWarehouse<ElementDamper> _element_dampers;

  MooseObjectWarehouse<NodalDamper> _nodal_dampers;

  MooseObjectWarehouse<GeneralDamper> _general_dampers;

  MooseObjectTagWarehouse<NodalKernelBase> _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;

  bool _use_finite_differenced_preconditioner;

  MatFDColoring _fdcoloring;

  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_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_pre_smo_residual;

  Real _pre_smo_residual;
  Real _initial_residual;
  bool _use_pre_smo_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);

  bool _computed_scaling;

  bool _compute_scaling_once;

  Real _resid_vs_jac_scaling_param;

  std::vector<std::vector<std::string>> _scaling_group_variables;

  std::vector<bool> _variable_autoscaled;

  std::vector<std::string> _ignore_variables_for_autoscaling;

  bool _off_diagonals_in_auto_scaling;

  std::unique_ptr<libMesh::DiagonalMatrix<Number>> _scaling_matrix;

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

  void setupScalingData();

  std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
      _undisplaced_mortar_functors;

  std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
      _displaced_mortar_functors;

  std::vector<NumericVector<Number> *> _solution_state;

  bool _auto_scaling_initd;

  std::unordered_map<unsigned int, unsigned int> _var_to_group_var;

  std::size_t _num_scaling_groups;
};