NewtonSolver.h

// Copyright (C) 2005-2021 Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#ifdef HAS_PETSC
 
#include <dolfinx/common/MPI.h>
#include <dolfinx/la/petsc.h>
#include <functional>
#include <memory>
#include <petscmat.h>
#include <petscvec.h>
#include <utility>
 
namespace dolfinx
{
 
namespace la::petsc
{
class KrylovSolver;
} // namespace la::petsc
 
namespace nls::petsc
{
class NewtonSolver
{
public:
  explicit NewtonSolver(MPI_Comm comm);
 
  NewtonSolver(NewtonSolver&& solver) = default;
 
  // Copy constructor (deleted)
  NewtonSolver(const NewtonSolver& solver) = delete;
 
  NewtonSolver& operator=(NewtonSolver&& solver) = default;
 
  // Assignment operator (deleted)
  NewtonSolver& operator=(const NewtonSolver& solver) = delete;
 
  ~NewtonSolver();
 
  void setF(std::function<void(const Vec, Vec)> F, Vec b);
 
  void setJ(std::function<void(const Vec, Mat)> J, Mat Jmat);
 
  void setP(std::function<void(const Vec, Mat)> P, Mat Pmat);
 
  const la::petsc::KrylovSolver& get_krylov_solver() const;
 
  la::petsc::KrylovSolver& get_krylov_solver();
 
  void set_form(std::function<void(Vec)> form);
 
  void set_convergence_check(
      std::function<std::pair<double, bool>(const NewtonSolver&, const Vec)> c);
 
  void set_update(
      std::function<void(const NewtonSolver& solver, const Vec, Vec)> update);
 
  std::pair<int, bool> solve(Vec x);
 
  int iteration() const;
 
  int krylov_iterations() const;
 
  double residual() const;
 
  double residual0() const;
 
  MPI_Comm comm() const;
 
  int max_it = 50;
 
  double rtol = 1e-9;
 
  double atol = 1e-10;
 
  std::string convergence_criterion = "residual";
 
  bool report = true;
 
  bool error_on_nonconvergence = true;
 
  double relaxation_parameter = 1.0;
 
private:
  // Function for computing the residual vector. The first argument is
  // the latest solution vector x and the second argument is the
  // residual vector.
  std::function<void(const Vec x, Vec b)> _fnF;
 
  // Function for computing the Jacobian matrix operator. The first
  // argument is the latest solution vector x and the second argument is
  // the matrix operator.
  std::function<void(const Vec x, Mat J)> _fnJ;
 
  // Function for computing the preconditioner matrix operator. The
  // first argument is the latest solution vector x and the second
  // argument is the matrix operator.
  std::function<void(const Vec x, Mat P)> _fnP;
 
  // Function called before the residual and Jacobian function at each
  // iteration.
  std::function<void(const Vec x)> _system;
 
  // Residual vector
  Vec _b = nullptr;
 
  // Jacobian matrix and preconditioner matrix
  Mat _matJ = nullptr, _matP = nullptr;
 
  // Function to check for convergence
  std::function<std::pair<double, bool>(const NewtonSolver& solver,
                                        const Vec r)>
      _converged;
 
  // Function to update the solution after inner solve for the Newton increment
  std::function<void(const NewtonSolver& solver, const Vec dx, Vec x)>
      _update_solution;
 
  // Accumulated number of Krylov iterations since solve began
  int _krylov_iterations;
 
  // Number of iterations
  int _iteration;
 
  // Most recent residual and initial residual
  double _residual, _residual0;
 
  // Linear solver
  la::petsc::KrylovSolver _solver;
 
  // Solution vector
  Vec _dx = nullptr;
 
  // MPI communicator
  dolfinx::MPI::Comm _comm;
};
} // namespace nls::petsc
} // namespace dolfinx
 
#endif