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;
class FieldSplitPreconditionerBase;
class Convergence;
 
// libMesh forward declarations
namespace libMesh
{
template <typename T>
class NumericVector;
template <typename T>
class SparseMatrix;
template <typename T>
class DiagonalMatrix;
class DofMapBase;
} // 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;
  void reinitMortarFunctors();
 
  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;
 
  bool haveFiniteDifferencedPreconditioner() const
  {
    return _use_finite_differenced_preconditioner;
  }
  bool haveFieldSplitPreconditioner() const { return _fsp; }
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  virtual void addKokkosKernel(const std::string & kernel_name,
                               const std::string & name,
                               InputParameters & parameters);
 
  virtual void addKokkosNodalKernel(const std::string & kernel_name,
                                    const std::string & name,
                                    InputParameters & parameters);
 
  void addKokkosBoundaryCondition(const std::string & bc_name,
                                  const std::string & name,
                                  InputParameters & parameters);
#endif
 
  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();
 
#ifdef MOOSE_KOKKOS_ENABLED
  void setKokkosInitialSolution();
#endif
 
  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);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void computeKokkosResidualAndJacobian(const std::set<TagID> & vector_tags,
                                        const std::set<TagID> & matrix_tags);
#endif
 
  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 useFieldSplitPreconditioner(FieldSplitPreconditionerBase * fsp) { _fsp = fsp; }
 
  FieldSplitPreconditionerBase & getFieldSplitPreconditioner();
 
  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;
 
  virtual unsigned int nNonlinearIterations() const { return _n_iters; }
 
  virtual unsigned int nLinearIterations() const { return _n_linear_iters; }
 
  unsigned int nResidualEvaluations() const { return _n_residual_evaluations; }
 
  virtual 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 needInternalNeighborSideMaterial(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;
  }
 
#ifdef MOOSE_KOKKOS_ENABLED
  MooseObjectTagWarehouse<ResidualObject> & getKokkosKernelWarehouse() { return _kokkos_kernels; }
  MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalKernelWarehouse()
  {
    return _kokkos_nodal_kernels;
  }
  MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalBCWarehouse()
  {
    return _kokkos_nodal_bcs;
  }
  MooseObjectTagWarehouse<ResidualObject> & getKokkosIntegratedBCWarehouse()
  {
    return _kokkos_integrated_bcs;
  }
#endif
 
 
  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();
 
  void setConvergenceName(const ConvergenceName & convergence_name)
  {
    _convergence_name = convergence_name;
  }
 
  Convergence & convergence();
 
protected:
  void computeResidualInternal(const std::set<TagID> & tags);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void computeKokkosResidual(const std::set<TagID> & tags);
  void computeKokkosNodalBCsResidual(const std::set<TagID> & tags);
#endif
 
  void computeNodalBCsResidual(NumericVector<Number> & residual);
 
  void computeNodalBCsResidual(NumericVector<Number> & residual, const std::set<TagID> & tags);
 
  void computeNodalBCsResidual(const std::set<TagID> & tags);
 
  void computeNodalBCsJacobian(const std::set<TagID> & tags);
 
  void computeNodalBCsResidualAndJacobian(const std::set<TagID> & vector_tags,
                                          const std::set<TagID> & matrix_tags);
 
  void computeJacobianInternal(const std::set<TagID> & tags);
 
#ifdef MOOSE_KOKKOS_ENABLED
  void computeKokkosJacobian(const std::set<TagID> & tags);
#endif
 
  void computeDiracContributions(const std::set<TagID> & vector_tags,
                                 const std::set<TagID> & matrix_tags,
                                 Moose::ComputeType compute_type);
 
  void computeScalarKernelsJacobians(const std::set<TagID> & tags);
 
  void enforceNodalConstraintsResidual(NumericVector<Number> & residual);
 
  bool enforceNodalConstraintsJacobian(const SparseMatrix<Number> & jacobian);
 
  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;
 
#ifdef MOOSE_KOKKOS_ENABLED
  MooseObjectTagWarehouse<ResidualObject> _kokkos_kernels;
  MooseObjectTagWarehouse<ResidualObject> _kokkos_integrated_bcs;
  MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_bcs;
  MooseObjectWarehouse<ResidualObject> _kokkos_preset_nodal_bcs;
  MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_kernels;
#endif
 
  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;
 
  FieldSplitPreconditionerBase * _fsp;
 
  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:
  std::vector<SetupInterface *> getFVSetupObjects(THREAD_ID tid);
 
  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;
 
  ConvergenceName _convergence_name;
};