Preconditioner reuse in petsc_nonlinear_solver

This alters the logical used in PetscNonlinearSolver to:

1. Not free the SNES instance between solves
2. Use SNESReset by default in between solves
3. Add two new options "reuse preconditioner" and "reuse preconditioner maximum iterations" to support preconditioner reuse

include/solvers/nonlinear_solver.h

@@ -357,6 +357,16 @@ class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
    */
   bool converged;
 
+  /**
+   * Whether we should reuse the linear preconditioner
+   */
+  bool reuse_preconditioner;
+
+  /**
+   * Number of linear iterations to retain the preconditioner
+   */
+  unsigned int reuse_preconditioner_max_its;
+
   /**
    * Set the solver configuration object.
    */
@@ -423,6 +433,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);
@@ -104,7 +105,7 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    * May assign a name to the solver in some implementations
    */
   virtual void init (const char * name = nullptr) override;
-
+  
   /**
    * \returns The raw PETSc snes context pointer.
    */
@@ -261,6 +262,16 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    */
   bool _computing_base_vector;
 
+  /**
+    *  Whether we have a live snes instance
+    */
+  bool _snes_allocated;
+
+  /**
+    * 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 +280,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/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);
+    if (ierr != 0) return ierr;
+
+    PetscInt niter;
+    ierr = KSPGetIterationNumber(ksp, &niter);
+    if (ierr != 0) return 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);
+      if (ierr !=0) return ierr;
+    }
+    return 0;
+  }
 
   //-------------------------------------------------------------------
   // this function is called by PETSc at the end of each nonlinear step
@@ -615,15 +642,24 @@ 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),
+  _snes_allocated(false),
+  _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_allocated) 
+    {
+    _snes.destroy();
+    _snes_allocated = false;
+    }
+}
 
 
 template <typename T>
@@ -632,9 +668,16 @@ void PetscNonlinearSolver<T>::clear ()
   if (this->initialized())
     {
       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 (!(this->reuse_preconditioner)) 
+        {
+        PetscErrorCode ierr=0;
+        ierr = SNESReset(_snes);
+        LIBMESH_CHKERR(ierr);
+        }
 
       // Reset the nonlinear iteration counter.  This information is only relevant
       // *during* the solve().  After the solve is completed it should return to
@@ -654,9 +697,17 @@ void PetscNonlinearSolver<T>::init (const char * name)
       this->_is_initialized = true;
 
       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_allocated) {
+        ierr = SNESCreate(this->comm().get(), _snes.get());
+        LIBMESH_CHKERR(ierr);
+        _snes_allocated = true;
+      }
+
+      // I believe all of the following can be safely repeated
+      // even on an old snes instance from the last solve
 
       if (name)
         {
@@ -731,7 +782,6 @@ void PetscNonlinearSolver<T>::init (const char * name)
 }
 
 
-
 template <typename T>
 SNES PetscNonlinearSolver<T>::snes()
 {
@@ -815,6 +865,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 +876,38 @@ 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 ((this->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
+      ierr = SNESMonitorSet(_snes, &libmesh_petsc_recalculate_monitor, 
+                            (void*)
+                            &(this->reuse_preconditioner_max_its),
+                            NULL);
+      LIBMESH_CHKERR(ierr);
+    }
+  else if (!(this->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);
+        }
+    }
+
   ierr = SNESSetFunction (_snes, r->vec(), libmesh_petsc_snes_residual, this);
   LIBMESH_CHKERR(ierr);
 
@@ -997,7 +1080,8 @@ 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.

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->reuse_preconditioner =  reuse_preconditioner;
+  nonlinear_solver->reuse_preconditioner_max_its = reuse_preconditioner_max_its; 
 
   if (diff_solver.get())
     {