Preconditioner reuse in petsc_nonlinear_solver

include/solvers/nonlinear_solver.h

@@ -354,7 +354,38 @@ class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
    */
   void set_solver_configuration(SolverConfiguration & solver_configuration);
 
+  /**
+   *  Get the reuse_preconditioner flag
+   */
+  virtual bool reuse_preconditioner() const;
+
+  /**
+   *  Set the reuse preconditioner flag
+   */
+  virtual void set_reuse_preconditioner(bool reuse);
+
+  /**
+   *  Get the reuse_preconditioner_max_its parameter
+   */
+  virtual unsigned int reuse_preconditioner_max_its() const;
+
+  /**
+   *  Set the reuse_preconditioner_max_its parameter
+   */
+  virtual void set_reuse_preconditioner_max_its(unsigned int i);
+
+
 protected:
+  /**
+   * Whether we should reuse the linear preconditioner
+   */
+  bool _reuse_preconditioner;
+
+  /**
+   * Number of linear iterations to retain the preconditioner
+   */
+  unsigned int _reuse_preconditioner_max_its;
+
   /**
    * A reference to the system we are solving.
    */
@@ -415,6 +446,8 @@ NonlinearSolver<T>::NonlinearSolver (sys_type & s) :
   initial_linear_tolerance(0),
   minimum_linear_tolerance(0),
   converged(false),
+  _reuse_preconditioner(false),
+  _reuse_preconditioner_max_its(0),
   _system(s),
   _is_initialized (false),
   _preconditioner (nullptr),

include/solvers/petsc_nonlinear_solver.h

@@ -47,6 +47,7 @@ class ResidualContext;
 // need access to these most of the time as they are used internally by this object.
 extern "C"
 {
+  PetscErrorCode libmesh_petsc_recalculate_monitor(SNES snes, PetscInt it, PetscReal norm, void* mctx);
   PetscErrorCode libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *);
   PetscErrorCode libmesh_petsc_snes_residual (SNES, Vec x, Vec r, void * ctx);
   PetscErrorCode libmesh_petsc_snes_fd_residual (SNES, Vec x, Vec r, void * ctx);
@@ -201,6 +202,21 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    */
   std::unique_ptr<ComputeLineSearchObject> linesearch_object;
 
+  /**
+    * Setup the default monitor if required
+    */
+  void setup_default_monitor();
+
+  /**
+   * Getter for preconditioner reuse
+   */
+  virtual bool reuse_preconditioner() const override;
+
+  /**
+   *  Getter for the maximum iterations flag for preconditioner reuse
+   */
+  virtual unsigned int reuse_preconditioner_max_its() const override;
+
 protected:
 
   /**
@@ -261,6 +277,11 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    */
   bool _computing_base_vector;
 
+  /**
+    * Whether we've triggered the preconditioner reuse
+    */
+  bool _setup_reuse;
+
 private:
   friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
   friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);
@@ -269,8 +290,6 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
   friend PetscErrorCode libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
 };
 
-
-
 } // namespace libMesh
 
 

src/solvers/nonlinear_solver.C

@@ -84,6 +84,29 @@ void NonlinearSolver<T>::set_solver_configuration(SolverConfiguration & solver_c
   _solver_configuration = &solver_configuration;
 }
 
+template <typename T>
+bool NonlinearSolver<T>::reuse_preconditioner() const
+{
+  libmesh_not_implemented();
+}
+
+template <typename T>
+void NonlinearSolver<T>::set_reuse_preconditioner(bool reuse)
+{
+  _reuse_preconditioner = reuse;
+}
+
+template <typename T>
+unsigned int NonlinearSolver<T>::reuse_preconditioner_max_its() const
+{
+  libmesh_not_implemented();
+}
+
+template <typename T>
+void NonlinearSolver<T>::set_reuse_preconditioner_max_its(unsigned int i)
+{
+  _reuse_preconditioner_max_its = i;
+}
 
 //------------------------------------------------------------------
 // Explicit instantiations

src/solvers/petsc_nonlinear_solver.C

@@ -105,6 +105,33 @@ libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx)
 // Give them an obscure name to avoid namespace pollution.
 extern "C"
 {
+  // -----------------------------------------------------------------
+  // this function monitors the nonlinear solve and checks to see
+  // if we want to recalculate the preconditioner.  It only gets
+  // added to the SNES instance if we're reusing the preconditioner
+  PetscErrorCode
+  libmesh_petsc_recalculate_monitor(SNES snes, PetscInt, PetscReal, void* mctx)
+  {
+    unsigned int * max_iters = (unsigned int *) mctx;
+    PetscErrorCode ierr = 0;
+
+    KSP ksp;
+    ierr = SNESGetKSP(snes, &ksp);
+    LIBMESH_CHKERR(ierr);
+
+    PetscInt niter;
+    ierr = KSPGetIterationNumber(ksp, &niter);
+    LIBMESH_CHKERR(ierr);
+
+    if (niter > *max_iters)
+    {
+      // -2 is a magic number for "recalculate next time you need it
+      // and then not again"
+      ierr = SNESSetLagPreconditioner(snes, -2);
+      LIBMESH_CHKERR(ierr);
+    }
+    return 0;
+  }
 
   //-------------------------------------------------------------------
   // this function is called by PETSc at the end of each nonlinear step
@@ -615,15 +642,22 @@ PetscNonlinearSolver<T>::PetscNonlinearSolver (sys_type & system_in) :
   _zero_out_jacobian(true),
   _default_monitor(true),
   _snesmf_reuse_base(true),
-  _computing_base_vector(true)
+  _computing_base_vector(true),
+  _setup_reuse(false)
 {
 }
 
 
 
 template <typename T>
-PetscNonlinearSolver<T>::~PetscNonlinearSolver () = default;
-
+PetscNonlinearSolver<T>::~PetscNonlinearSolver ()
+{
+  // Take care of the SNES instance, if it was allocated
+  if (_snes)
+    {
+    _snes.destroy();
+    }
+}
 
 
 template <typename T>
@@ -633,8 +667,16 @@ void PetscNonlinearSolver<T>::clear ()
     {
       this->_is_initialized = false;
 
-      // Calls custom deleter
-      _snes.destroy();
+      // If we don't need the preconditioner next time
+      // retain the original behavior of clearing the data
+      // between solves.
+      if (!(reuse_preconditioner()))
+        {
+        // SNESReset really ought to work but replacing destory() with
+        // SNESReset causes a very slight change in behavior that
+        // manifests as two failed MOOSE tests...
+        _snes.destroy();
+        }
 
       // Reset the nonlinear iteration counter.  This information is only relevant
       // *during* the solve().  After the solve is completed it should return to
@@ -655,8 +697,15 @@ void PetscNonlinearSolver<T>::init (const char * name)
 
       PetscErrorCode ierr=0;
 
-      ierr = SNESCreate(this->comm().get(), _snes.get());
-      LIBMESH_CHKERR(ierr);
+      // Make only if we don't already have a retained snes
+      // hanging around from the last solve
+      if (!_snes) {
+        ierr = SNESCreate(this->comm().get(), _snes.get());
+        LIBMESH_CHKERR(ierr);
+      }
+
+      // I believe all of the following can be safely repeated
+      // even on an old snes instance from the last solve
 
       if (name)
         {
@@ -680,12 +729,7 @@ void PetscNonlinearSolver<T>::init (const char * name)
         // SNES now owns the reference to dm.
       }
 
-      if (_default_monitor)
-        {
-          ierr = SNESMonitorSet (_snes, libmesh_petsc_snes_monitor,
-                                 this, PETSC_NULL);
-          LIBMESH_CHKERR(ierr);
-        }
+      setup_default_monitor();
 
       // If the SolverConfiguration object is provided, use it to set
       // options during solver initialization.
@@ -731,7 +775,6 @@ void PetscNonlinearSolver<T>::init (const char * name)
 }
 
 
-
 template <typename T>
 SNES PetscNonlinearSolver<T>::snes()
 {
@@ -815,6 +858,7 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditi
   LOG_SCOPE("solve()", "PetscNonlinearSolver");
   this->init ();
 
+
   // Make sure the data passed in are really of Petsc types
   PetscMatrix<T> * pre = cast_ptr<PetscMatrix<T> *>(&pre_in);
   PetscVector<T> * x   = cast_ptr<PetscVector<T> *>(&x_in);
@@ -825,6 +869,41 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditi
   // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
   Real final_residual_norm=0.;
 
+  // We don't want to do this twice because it resets
+  // SNESSetLagPreconditioner
+  if ((reuse_preconditioner()) && (!_setup_reuse))
+    {
+      _setup_reuse = true;
+      ierr = SNESSetLagPreconditionerPersists(_snes, PETSC_TRUE);
+      LIBMESH_CHKERR(ierr);
+      // According to the PETSC 3.16.5 docs -2 is a magic number which
+      // means "recalculate the next time you need it and then not again"
+      ierr = SNESSetLagPreconditioner(_snes, -2);
+      LIBMESH_CHKERR(ierr);
+      // Add in our callback which will trigger recalculating
+      // the preconditioner when we hit reuse_preconditioner_max_its
+      unsigned int max_its = reuse_preconditioner_max_its();
+      ierr = SNESMonitorSet(_snes, &libmesh_petsc_recalculate_monitor,
+                            (void*)
+                            &max_its,
+                            NULL);
+      LIBMESH_CHKERR(ierr);
+    }
+  else if (!(reuse_preconditioner()))
+    // This covers the case where it was enabled but was then disabled
+    {
+      ierr = SNESSetLagPreconditionerPersists(_snes, PETSC_FALSE);
+      LIBMESH_CHKERR(ierr);
+      if (_setup_reuse)
+        {
+          _setup_reuse = false;
+          ierr = SNESMonitorCancel(_snes);
+          LIBMESH_CHKERR(ierr);
+          // Readd default monitor
+          setup_default_monitor();
+        }
+    }
+
   ierr = SNESSetFunction (_snes, r->vec(), libmesh_petsc_snes_residual, this);
   LIBMESH_CHKERR(ierr);
 
@@ -998,6 +1077,7 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditi
   //Based on Petsc 2.3.3 documentation all diverged reasons are negative
   this->converged = (_reason >= 0);
 
+  // Reset data structure
   this->clear();
 
   // return the # of its. and the final residual norm.
@@ -1036,6 +1116,29 @@ int PetscNonlinearSolver<T>::get_total_linear_iterations()
   return _n_linear_iterations;
 }
 
+template <typename T>
+void PetscNonlinearSolver<T>::setup_default_monitor()
+{
+  if (_default_monitor)
+    {
+      PetscErrorCode ierr = SNESMonitorSet (_snes, libmesh_petsc_snes_monitor,
+                                            this, PETSC_NULL);
+      LIBMESH_CHKERR(ierr);
+    }
+}
+
+template <typename T>
+bool PetscNonlinearSolver<T>::reuse_preconditioner() const
+{
+  return this->_reuse_preconditioner;
+}
+
+template <typename T>
+unsigned int PetscNonlinearSolver<T>::reuse_preconditioner_max_its() const
+{
+  return this->_reuse_preconditioner_max_its;
+}
+
 
 //------------------------------------------------------------------
 // Explicit instantiations

src/systems/nonlinear_implicit_system.C

@@ -52,6 +52,9 @@ NonlinearImplicitSystem::NonlinearImplicitSystem (EquationSystems & es,
   es.parameters.set<Real>("nonlinear solver divergence tolerance") = 1e+4;
   es.parameters.set<Real>("nonlinear solver absolute step tolerance") = 1e-8;
   es.parameters.set<Real>("nonlinear solver relative step tolerance") = 1e-8;
+
+  es.parameters.set<bool>("reuse preconditioner") = false;
+  es.parameters.set<unsigned int>("reuse preconditioner maximum iterations") = 1;
 }
 
 
@@ -129,6 +132,11 @@ void NonlinearImplicitSystem::set_solver_parameters ()
   const double linear_min_tol =
     double(es.parameters.get<Real>("linear solver minimum tolerance"));
 
+  const bool reuse_preconditioner =
+      es.parameters.get<bool>("reuse preconditioner");
+  const unsigned int reuse_preconditioner_max_its =
+      es.parameters.get<unsigned int>("reuse preconditioner maximum iterations");
+
   // Set all the parameters on the NonlinearSolver
   nonlinear_solver->max_nonlinear_iterations = maxits;
   nonlinear_solver->max_function_evaluations = maxfuncs;
@@ -140,6 +148,8 @@ void NonlinearImplicitSystem::set_solver_parameters ()
   nonlinear_solver->max_linear_iterations = maxlinearits;
   nonlinear_solver->initial_linear_tolerance = linear_tol;
   nonlinear_solver->minimum_linear_tolerance = linear_min_tol;
+  nonlinear_solver->set_reuse_preconditioner(reuse_preconditioner);
+  nonlinear_solver->set_reuse_preconditioner_max_its(reuse_preconditioner_max_its);
 
   if (diff_solver.get())
     {