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
 
#include <dolfinx/common/MPI.h>
#include <dolfinx/la/PETScKrylovSolver.h>
#include <functional>
#include <memory>
#include <petscmat.h>
#include <petscvec.h>
#include <utility>
 
namespace dolfinx
{
 
namespace la
{
class PETScKrylovSolver;
} // namespace la
 
namespace nls
{
 
 
class NewtonSolver
{
public:
  explicit NewtonSolver(MPI_Comm comm);
 
  // Move constructor (deleted)
  NewtonSolver(NewtonSolver&& solver) = delete;
 
  // Copy constructor (deleted)
  NewtonSolver(const NewtonSolver& solver) = delete;
 
  // Assignment operator (deleted)
  NewtonSolver& operator=(const NewtonSolver& solver) = delete;
 
  // Move assignment constructor (deleted)
  NewtonSolver& operator=(const NewtonSolver&& solver) = delete;
 
  ~NewtonSolver();
 
  void setF(const std::function<void(const Vec, Vec)>& F, Vec b);
 
  void setJ(const std::function<void(const Vec, Mat)>& J, Mat Jmat);
 
  void setP(const std::function<void(const Vec, Mat)>& P, Mat Pmat);
 
  const la::PETScKrylovSolver& get_krylov_solver() const;
 
  la::PETScKrylovSolver& get_krylov_solver();
 
  void set_form(const std::function<void(Vec x)>& form);
 
  void
  set_convergence_check(const std::function<std::pair<double, bool>(
                            const nls::NewtonSolver& solver, const Vec r)>& c);
 
  void set_update(const std::function<void(const nls::NewtonSolver& solver,
                                           const Vec dx, Vec x)>& update);
 
  std::pair<int, bool> solve(Vec x);
 
  int iteration() const;
 
  int krylov_iterations() const;
 
  double residual() const;
 
  double residual0() const;
 
  MPI_Comm mpi_comm() const;
 
  int max_it = 50;
 
  double rtol = 1e-9;
 
  double atol = 1e-10;
 
  // FIXME: change to string to enum
  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 nls::NewtonSolver& solver,
                                        const Vec r)>
      _converged;
 
  // Function to update the solution once convergence is reached
  std::function<void(const nls::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::PETScKrylovSolver _solver;
 
  // Solution vector
  Vec _dx = nullptr;
 
  // MPI communicator
  dolfinx::MPI::Comm _mpi_comm;
};
} // namespace nls
} // namespace dolfinx