Nonlinear solver generalizations

include/numerics/petsc_matrix_shell_matrix.h

@@ -58,6 +58,27 @@ class PetscMatrixShellMatrix : public PetscMatrixBase<T>
 
   virtual bool require_sparsity_pattern() const override { return false; }
 
+  virtual void zero() override;
+  virtual std::unique_ptr<SparseMatrix<T>> zero_clone() const override;
+  virtual std::unique_ptr<SparseMatrix<T>> clone() const override;
+  virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override;
+  virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override;
+  virtual void add_matrix(const DenseMatrix<T> & dm,
+                          const std::vector<numeric_index_type> & rows,
+                          const std::vector<numeric_index_type> & cols) override;
+  virtual void add_matrix(const DenseMatrix<T> & dm,
+                          const std::vector<numeric_index_type> & dof_indices) override;
+  virtual void add(const T a, const SparseMatrix<T> & X) override;
+  virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override;
+  virtual Real l1_norm() const override;
+  virtual Real linfty_norm() const override;
+  virtual void print_personal(std::ostream & os = libMesh::out) const override;
+  virtual void get_diagonal(NumericVector<T> & dest) const override;
+  virtual void get_transpose(SparseMatrix<T> & dest) const override;
+  virtual void get_row(numeric_index_type i,
+                       std::vector<numeric_index_type> & indices,
+                       std::vector<T> & values) const override;
+
 private:
   // Make this private because we mark as initialized after we've done our initialization, and we
   // don't want derived classes to mistakenly register their data as initialized (or not)

include/numerics/petsc_mffd_matrix.h

@@ -48,7 +48,15 @@ class PetscMFFDMatrix : public PetscMatrixBase<T>
 
   explicit PetscMFFDMatrix(const Parallel::Communicator & comm_in);
 
-  PetscMFFDMatrix & operator=(Mat m);
+  /**
+   * Adopt an existing, externally-owned Mat, without destroying it when this
+   * object goes out of scope. \p set_context controls whether we attach a
+   * context pointer to \p m allowing \p get_context() to recover this object
+   * from the Mat later; skip this when this wrapper is short-lived (e.g. a
+   * function-local variable) so we don't leave a dangling context on \p m
+   * after we're destroyed.
+   */
+  void assign(Mat m, bool set_context);
 
   virtual void init(const numeric_index_type,
                     const numeric_index_type,
@@ -101,11 +109,14 @@ PetscMFFDMatrix<T>::PetscMFFDMatrix(const Parallel::Communicator & comm_in)
 }
 
 template <typename T>
-PetscMFFDMatrix<T> &
-PetscMFFDMatrix<T>::operator=(Mat m)
+void
+PetscMFFDMatrix<T>::assign(Mat m, bool set_context)
 {
   this->_mat = m;
-  return *this;
+  this->_is_initialized = true;
+  this->_destroy_mat_on_exit = false;
+  if (set_context)
+    this->set_context();
 }
 
 template <typename T>

include/solvers/nonlinear_solver.h

@@ -102,6 +102,23 @@ class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
                                                const double,      // Stopping tolerance
                                                const unsigned int) = 0; // N. Iterations
 
+  /**
+   * Solves the nonlinear system using \p jac_in as the actual Jacobian operator and \p pre_in as
+   * the preconditioning matrix -- which may be the same object (the common case) or genuinely
+   * distinct (e.g. a matrix-free operator paired with an assembled preconditioning matrix). The
+   * default implementation ignores \p jac_in and simply forwards to the single-matrix solve(),
+   * for solver backends that do not distinguish between the two.
+   */
+  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> & /* jac_in */,
+                                               SparseMatrix<T> & pre_in,
+                                               NumericVector<T> & x_in,
+                                               NumericVector<T> & r_in,
+                                               const double tol,
+                                               const unsigned int m_its)
+  {
+    return this->solve(pre_in, x_in, r_in, tol, m_its);
+  }
+
   /**
    * Prints a useful message about why the latest nonlinear solve
    * con(di)verged.

include/solvers/petsc_nonlinear_solver.h

@@ -107,11 +107,23 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
   SNES snes(const char * name = nullptr);
 
   /**
-   * Call the Petsc solver.  It calls the method below, using the
-   * same matrix for the system and preconditioner matrices.
+   * Call the Petsc solver, using the same matrix for the system and preconditioner matrices. Calls
+   * the two-matrix overload below with \p pre_in for both.
    */
   virtual std::pair<unsigned int, Real>
-  solve (SparseMatrix<T> &,                     // System Jacobian Matrix
+  solve (SparseMatrix<T> & pre_in,               // System Preconditioning Matrix
+         NumericVector<T> &,                    // Solution vector
+         NumericVector<T> &,                    // Residual vector
+         const double,                         // Stopping tolerance
+         const unsigned int) override; // N. Iterations
+
+  /**
+   * Call the Petsc solver, using \p jac_in as the actual SNES Jacobian operator (Amat) and
+   * \p pre_in as the preconditioning matrix (Pmat).
+   */
+  virtual std::pair<unsigned int, Real>
+  solve (SparseMatrix<T> & jac_in,               // Jacobian operator matrix (Amat)
+         SparseMatrix<T> & pre_in,               // Preconditioning matrix (Pmat)
          NumericVector<T> &,                    // Solution vector
          NumericVector<T> &,                    // Residual vector
          const double,                         // Stopping tolerance
@@ -291,9 +303,6 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
   PetscDMWrapper _dm_wrapper;
 #endif
 
-  /// Wrapper for matrix-free finite-difference Jacobians
-  PetscMFFDMatrix<Number> _mffd_jac;
-
 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);

include/systems/nonlinear_implicit_system.h

@@ -257,10 +257,20 @@ class NonlinearImplicitSystem : public ImplicitSystem
   virtual void reinit () override;
 
   /**
-   * Assembles & solves the nonlinear system R(x) = 0.
+   * Assembles & solves the nonlinear system R(x) = 0. If a matrix has been registered via
+   * set_operator_matrix(), it is used as the actual Jacobian operator (Amat) while \p matrix
+   * remains the preconditioning matrix (Pmat); otherwise \p matrix is used for both, as usual.
    */
   virtual void solve () override;
 
+  /**
+   * Set a matrix to use as the actual Jacobian operator (Amat) on the next solve(), distinct from
+   * \p matrix, which continues to be used as the preconditioning matrix (Pmat). Pass nullptr (the
+   * default) to restore the ordinary behavior of using \p matrix for both. Only solve() itself
+   * reads this value back, so there is no public getter.
+   */
+  void set_operator_matrix(SparseMatrix<Number> * mat) { _operator_matrix = mat; }
+
   /**
    * \returns An integer corresponding to the upper iteration count
    * limit and a Real corresponding to the convergence tolerance to
@@ -335,6 +345,12 @@ class NonlinearImplicitSystem : public ImplicitSystem
    * The final residual for the nonlinear system R(x)
    */
   Real _final_nonlinear_residual;
+
+  /**
+   * An optional matrix to use as the actual Jacobian operator (Amat), distinct from \p matrix
+   * (used as the preconditioning matrix, Pmat). See set_operator_matrix().
+   */
+  SparseMatrix<Number> * _operator_matrix;
 };
 
 } // namespace libMesh

src/numerics/petsc_matrix_shell_matrix.C

@@ -54,6 +54,116 @@ PetscMatrixShellMatrix<T>::init(ParallelType libmesh_dbg_var(type))
   this->set_context();
 }
 
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::zero()
+{
+  libmesh_error();
+}
+
+template <typename T>
+std::unique_ptr<SparseMatrix<T>>
+PetscMatrixShellMatrix<T>::zero_clone() const
+{
+  libmesh_error();
+}
+
+template <typename T>
+std::unique_ptr<SparseMatrix<T>>
+PetscMatrixShellMatrix<T>::clone() const
+{
+  libmesh_not_implemented();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::set(const numeric_index_type, const numeric_index_type, const T)
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::add(const numeric_index_type, const numeric_index_type, const T)
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
+                                      const std::vector<numeric_index_type> &,
+                                      const std::vector<numeric_index_type> &)
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
+                                      const std::vector<numeric_index_type> &)
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::add(const T, const SparseMatrix<T> &)
+{
+  libmesh_error();
+}
+
+template <typename T>
+T
+PetscMatrixShellMatrix<T>::operator()(const numeric_index_type, const numeric_index_type) const
+{
+  libmesh_error();
+}
+
+template <typename T>
+Real
+PetscMatrixShellMatrix<T>::l1_norm() const
+{
+  libmesh_error();
+}
+
+template <typename T>
+Real
+PetscMatrixShellMatrix<T>::linfty_norm() const
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::print_personal(std::ostream &) const
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::get_diagonal(NumericVector<T> &) const
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::get_transpose(SparseMatrix<T> &) const
+{
+  libmesh_error();
+}
+
+template <typename T>
+void
+PetscMatrixShellMatrix<T>::get_row(numeric_index_type,
+                                   std::vector<numeric_index_type> &,
+                                   std::vector<T> &) const
+{
+  libmesh_error();
+}
+
 template class LIBMESH_EXPORT PetscMatrixShellMatrix<Number>;
 
 } // namespace libMesh

src/solvers/petsc_nonlinear_solver.C

@@ -450,7 +450,9 @@ extern "C"
       {
         libmesh_assert(!Jac);
         Jac = &mffd_jac;
-        mffd_jac = jac;
+        // mffd_jac is function-local, so don't attach a context to jac here -- it would
+        // dangle once mffd_jac is destroyed at the end of this call.
+        mffd_jac.assign(jac, /*set_context=*/false);
       }
 
     // We already computed the Jacobian during the residual evaluation
@@ -696,8 +698,7 @@ PetscNonlinearSolver<T>::PetscNonlinearSolver (sys_type & system_in) :
   _default_monitor(true),
   _snesmf_reuse_base(true),
   _computing_base_vector(true),
-  _setup_reuse(false),
-  _mffd_jac(this->_communicator)
+  _setup_reuse(false)
 {
 }
 
@@ -899,6 +900,18 @@ PetscNonlinearSolver<T>::build_mat_null_space(NonlinearImplicitSystem::ComputeVe
 template <typename T>
 std::pair<unsigned int, Real>
 PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditioning Matrix
+                                NumericVector<T> & x_in,    // Solution vector
+                                NumericVector<T> & r_in,    // Residual vector
+                                const double        tol,     // Stopping tolerance
+                                const unsigned int  m_its)
+{
+  return this->solve(pre_in, pre_in, x_in, r_in, tol, m_its);
+}
+
+template <typename T>
+std::pair<unsigned int, Real>
+PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  jac_in,  // Jacobian operator matrix (Amat)
+                                SparseMatrix<T> &  pre_in,  // Preconditioning matrix (Pmat)
                                 NumericVector<T> & x_in,    // Solution vector
                                 NumericVector<T> & r_in,    // Residual vector
                                 const double,              // Stopping tolerance
@@ -910,6 +923,7 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditi
   this->init ();
 
   // Make sure the data passed in are really of Petsc types
+  PetscMatrixBase<T> * jac = cast_ptr<PetscMatrixBase<T> *>(&jac_in);
   PetscMatrixBase<T> * pre = cast_ptr<PetscMatrixBase<T> *>(&pre_in);
   PetscVector<T> * x   = cast_ptr<PetscVector<T> *>(&x_in);
   PetscVector<T> * r   = cast_ptr<PetscVector<T> *>(&r_in);
@@ -951,7 +965,7 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  pre_in,  // System Preconditi
   // Only set the jacobian function if we've been provided with something to call.
   // This allows a user to set their own jacobian function if they want to
   if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
-    LibmeshPetscCall(SNESSetJacobian (_snes, pre->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));
+    LibmeshPetscCall(SNESSetJacobian (_snes, jac->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));
 
   // Have the Krylov subspace method use our good initial guess rather than 0
   KSP ksp;

src/systems/nonlinear_implicit_system.C

@@ -38,7 +38,8 @@ NonlinearImplicitSystem::NonlinearImplicitSystem (EquationSystems & es,
   nonlinear_solver          (NonlinearSolver<Number>::build(*this)),
   diff_solver               (),
   _n_nonlinear_iterations   (0),
-  _final_nonlinear_residual (1.e20)
+  _final_nonlinear_residual (1.e20),
+  _operator_matrix          (nullptr)
 {
   // Set default parameters
   // These were chosen to match the Petsc defaults
@@ -223,10 +224,20 @@ void NonlinearImplicitSystem::solve ()
       // Solve the nonlinear system.
       // Store the number of nonlinear iterations required to
       // solve and the final residual.
-      std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
-        nonlinear_solver->solve (*matrix, *solution, *rhs,
-                                 nonlinear_solver->relative_residual_tolerance,
-                                 nonlinear_solver->max_linear_iterations);
+      //
+      // If a distinct Jacobian operator matrix has been registered (see
+      // set_operator_matrix()), use it as Amat while *matrix remains the preconditioning matrix
+      // (Pmat); otherwise use the ordinary single-matrix solve.
+      if (_operator_matrix)
+        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
+          nonlinear_solver->solve (*_operator_matrix, *matrix, *solution, *rhs,
+                                   nonlinear_solver->relative_residual_tolerance,
+                                   nonlinear_solver->max_linear_iterations);
+      else
+        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
+          nonlinear_solver->solve (*matrix, *solution, *rhs,
+                                   nonlinear_solver->relative_residual_tolerance,
+                                   nonlinear_solver->max_linear_iterations);
     }
 
   // Update the system after the solve