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