nonlinear.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
 
#include <deal.II/lac/petsc_snes.h>
 
#include <deal.II/sundials/kinsol.h>
 
#include <deal.II/trilinos/nox.h>
 
 
DEAL_II_NAMESPACE_OPEN
 
template <typename VectorType = Vector<double>>
class NonlinearSolverSelector
{
public:
  class AdditionalData
  {
  public:
    enum SolutionStrategy
    {
      newton,
      linesearch,
      picard,
    };
 
    enum SolverType
    {
      automatic,
      kinsol,
      nox,
      /*
       * Use the PETSc SNES solver.
       */
      petsc_snes
    };
 
    AdditionalData(const SolverType &      solver_type = automatic,
                   const SolutionStrategy &strategy    = linesearch,
                   const unsigned int      maximum_non_linear_iterations = 200,
                   const double            function_tolerance            = 1e-8,
                   const double            relative_tolerance            = 1e-5,
                   const double            step_tolerance                = 0.0,
                   const unsigned int      anderson_subspace_size        = 0);
 
    SolverType solver_type;
 
    SolutionStrategy strategy;
 
    unsigned int maximum_non_linear_iterations;
 
    double function_tolerance;
 
    const double relative_tolerance;
 
    double step_tolerance;
 
    unsigned int anderson_subspace_size;
  };
 
  NonlinearSolverSelector();
 
  NonlinearSolverSelector(const AdditionalData &additional_data);
 
  NonlinearSolverSelector(const AdditionalData &additional_data,
                          const MPI_Comm &      mpi_communicator);
 
  void
  select(const typename AdditionalData::SolverType &type);
 
  void
  set_data(const AdditionalData &additional_data);
 
#ifdef DEAL_II_WITH_TRILINOS
  void
  set_data(
    const typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
      &                                         additional_data,
    const Teuchos::RCP<Teuchos::ParameterList> &parameters =
      Teuchos::rcp(new Teuchos::ParameterList));
#endif
 
#ifdef DEAL_II_WITH_SUNDIALS
  void
  set_data(const typename SUNDIALS::KINSOL<VectorType>::AdditionalData
             &additional_data);
#endif
 
#ifdef DEAL_II_WITH_PETSC
  void
  set_data(const typename PETScWrappers::NonlinearSolverData &additional_data);
#endif
 
  void
  solve(VectorType &initial_guess_and_solution);
 
  std::function<void(VectorType &)> reinit_vector;
 
  std::function<void(const VectorType &src, VectorType &dst)> residual;
 
  std::function<void(const VectorType &current_u)> setup_jacobian;
 
  std::function<
    void(const VectorType &rhs, VectorType &dst, const double tolerance)>
    solve_with_jacobian;
 
private:
  AdditionalData additional_data;
 
private:
  MPI_Comm mpi_communicator;
 
#ifdef DEAL_II_WITH_SUNDIALS
  typename SUNDIALS::KINSOL<VectorType>::AdditionalData additional_data_kinsol;
#endif
 
#ifdef DEAL_II_WITH_TRILINOS
  typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
                                       additional_data_nox;
  Teuchos::RCP<Teuchos::ParameterList> parameters_nox =
    Teuchos::rcp(new Teuchos::ParameterList);
#endif
 
#ifdef DEAL_II_WITH_PETSC
  typename PETScWrappers::NonlinearSolverData additional_data_petsc_snes;
#endif
 
  void
  solve_with_petsc(VectorType &initial_guess_and_solution);
};
 
 
 
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::set_data(
  const AdditionalData &additional_data)
{
#ifdef DEAL_II_WITH_SUNDIALS
  // These if statements pass on the strategy to the other nonlinear solvers
  if (additional_data.strategy ==
      NonlinearSolverSelector<VectorType>::AdditionalData::linesearch)
    additional_data_kinsol.strategy =
      SUNDIALS::KINSOL<VectorType>::AdditionalData::linesearch;
  else if (additional_data.strategy ==
           NonlinearSolverSelector<VectorType>::AdditionalData::newton)
    additional_data_kinsol.strategy =
      SUNDIALS::KINSOL<VectorType>::AdditionalData::newton;
  else if (additional_data.strategy ==
           NonlinearSolverSelector<VectorType>::AdditionalData::picard)
    additional_data_kinsol.strategy =
      SUNDIALS::KINSOL<VectorType>::AdditionalData::picard;
 
  // Setting data points in the KINSOL class from the NonlinearSolverSelector
  // class
  additional_data_kinsol.maximum_non_linear_iterations =
    additional_data.maximum_non_linear_iterations;
  additional_data_kinsol.function_tolerance =
    additional_data.function_tolerance;
  additional_data_kinsol.step_tolerance = additional_data.step_tolerance;
  additional_data_kinsol.anderson_subspace_size =
    additional_data.anderson_subspace_size;
#endif
 
// Do the same thing we did above but with NOX
#ifdef DEAL_II_WITH_TRILINOS
  parameters_nox->set("Nonlinear Solver", "Line Search Based");
  Teuchos::ParameterList &Line_Search = parameters_nox->sublist("Line Search");
  Line_Search.set("Method", "Full Step");
 
  additional_data_nox.max_iter = additional_data.maximum_non_linear_iterations;
  additional_data_nox.abs_tol  = additional_data.function_tolerance;
  additional_data_nox.rel_tol  = additional_data.relative_tolerance;
#endif
 
#ifdef DEAL_II_WITH_PETSC
  additional_data_petsc_snes.options_prefix = "";
 
  if (additional_data.anderson_subspace_size > 0 &&
      additional_data.strategy ==
        NonlinearSolverSelector<VectorType>::AdditionalData::picard)
    {
      additional_data_petsc_snes.snes_type = "anderson";
      // TODO additional_data.anderson_subspace_size;
    }
  else if (additional_data.strategy ==
           NonlinearSolverSelector<VectorType>::AdditionalData::linesearch)
    {
      additional_data_petsc_snes.snes_type            = "newtonls";
      additional_data_petsc_snes.snes_linesearch_type = "bt";
    }
  else if (additional_data.strategy ==
           NonlinearSolverSelector<VectorType>::AdditionalData::newton)
    {
      additional_data_petsc_snes.snes_linesearch_type = "newtonls";
      additional_data_petsc_snes.snes_linesearch_type = "basic";
    }
  else if (additional_data.strategy ==
           NonlinearSolverSelector<VectorType>::AdditionalData::picard)
    additional_data_petsc_snes.snes_type = "nrichardson";
 
  additional_data_petsc_snes.absolute_tolerance =
    additional_data.function_tolerance;
  additional_data_petsc_snes.relative_tolerance =
    additional_data.relative_tolerance;
  additional_data_petsc_snes.step_tolerance = additional_data.step_tolerance;
  additional_data_petsc_snes.maximum_non_linear_iterations =
    additional_data.maximum_non_linear_iterations;
  additional_data_petsc_snes.max_n_function_evaluations = -1;
 
#endif
}
 
 
 
template <typename VectorType>
NonlinearSolverSelector<VectorType>::NonlinearSolverSelector() = default;
 
 
 
template <typename VectorType>
NonlinearSolverSelector<VectorType>::NonlinearSolverSelector(
  const AdditionalData &additional_data)
  : additional_data(additional_data)
{
  set_data(additional_data);
}
 
 
 
template <typename VectorType>
NonlinearSolverSelector<VectorType>::NonlinearSolverSelector(
  const AdditionalData &additional_data,
  const MPI_Comm &      mpi_communicator)
  : additional_data(additional_data)
  , mpi_communicator(mpi_communicator)
{
  set_data(additional_data);
}
 
 
 
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::select(
  const typename AdditionalData::SolverType &type)
{
  additional_data.solver_type = type;
}
 
 
 
template <typename VectorType>
NonlinearSolverSelector<VectorType>::AdditionalData::AdditionalData(
  const SolverType &      solver_type,
  const SolutionStrategy &strategy,
  const unsigned int      maximum_non_linear_iterations,
  const double            function_tolerance,
  const double            relative_tolerance,
  const double            step_tolerance,
  const unsigned int      anderson_subspace_size)
  : solver_type(solver_type)
  , strategy(strategy)
  , maximum_non_linear_iterations(maximum_non_linear_iterations)
  , function_tolerance(function_tolerance)
  , relative_tolerance(relative_tolerance)
  , step_tolerance(step_tolerance)
  , anderson_subspace_size(anderson_subspace_size)
{}
 
 
 
#ifdef DEAL_II_WITH_TRILINOS
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::set_data(
  const typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
    &                                         additional_data,
  const Teuchos::RCP<Teuchos::ParameterList> &parameters)
{
  additional_data_nox = additional_data;
  parameters_nox      = parameters;
}
#endif
 
 
 
#ifdef DEAL_II_WITH_SUNDIALS
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::set_data(
  const typename SUNDIALS::KINSOL<VectorType>::AdditionalData &additional_data)
{
  additional_data_kinsol = additional_data;
}
#endif
 
 
 
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::solve_with_petsc(
  VectorType & /*initial_guess_and_solution*/)
{
  AssertThrow(false,
              ExcMessage("PETSc SNES requires you to use PETSc vectors."));
}
 
 
 
#ifdef DEAL_II_WITH_PETSC
template <>
void
NonlinearSolverSelector<PETScWrappers::MPI::Vector>::solve_with_petsc(
  PETScWrappers::MPI::Vector &initial_guess_and_solution)
{
  PETScWrappers::NonlinearSolver<PETScWrappers::MPI::Vector> nonlinear_solver(
    additional_data_petsc_snes, mpi_communicator);
 
  nonlinear_solver.residual = residual;
 
  nonlinear_solver.setup_jacobian = setup_jacobian;
 
  nonlinear_solver.solve_with_jacobian =
    [&](const PETScWrappers::MPI::Vector &src,
        PETScWrappers::MPI::Vector &      dst) {
      // PETSc does not gives a tolerance, so we have to choose something
      // reasonable to provide to the user:
      const double tolerance = 1e-6;
      solve_with_jacobian(src, dst, tolerance);
    };
 
  nonlinear_solver.solve(initial_guess_and_solution);
}
#endif
 
 
 
template <typename VectorType>
void
NonlinearSolverSelector<VectorType>::solve(
  VectorType &initial_guess_and_solution)
{
  // The "auto" solver_type will default to kinsol, however if KINSOL is not
  // available then we will use NOX.
  if (additional_data.solver_type == AdditionalData::SolverType::automatic)
    {
#ifdef DEAL_II_WITH_PETSC
      additional_data.solver_type = AdditionalData::SolverType::petsc_snes;
#endif
#ifdef DEAL_II_WITH_TRILINOS
      additional_data.solver_type = AdditionalData::SolverType::nox;
#endif
#ifdef DEAL_II_WITH_SUNDIALS
      additional_data.solver_type = AdditionalData::SolverType::kinsol;
#endif
 
      // If "auto" is still the solver type we cannot solve the problem
      if (additional_data.solver_type == AdditionalData::SolverType::automatic)
        AssertThrow(false, ExcMessage("No valid solver type."));
    }
 
  if (additional_data.solver_type == AdditionalData::SolverType::kinsol)
    {
#ifdef DEAL_II_WITH_SUNDIALS
      SUNDIALS::KINSOL<VectorType> nonlinear_solver(additional_data_kinsol,
                                                    mpi_communicator);
 
      nonlinear_solver.reinit_vector = reinit_vector;
      nonlinear_solver.residual      = residual;
 
      // We cannot simply set these two functions equal to each other
      // because they have a different number of inputs.
      nonlinear_solver.setup_jacobian = [&](const VectorType &current_u,
                                            const VectorType & /*current_f*/) {
        return NonlinearSolverSelector<VectorType>::setup_jacobian(current_u);
      };
 
      nonlinear_solver.solve_with_jacobian = solve_with_jacobian;
 
      nonlinear_solver.solve(initial_guess_and_solution);
#else
      AssertThrow(
        false, ExcMessage("You do not have SUNDIALS configured with deal.II!"));
#endif
    }
  else if (additional_data.solver_type == AdditionalData::SolverType::nox)
    {
#ifdef DEAL_II_WITH_TRILINOS
 
      TrilinosWrappers::NOXSolver<VectorType> nonlinear_solver(
        additional_data_nox, parameters_nox);
 
      // Do the same thing for NOX that we did with KINSOL.
      nonlinear_solver.residual = residual;
 
      // setup_jacobian for NOX has the same number of arguments for the same
      // function in NonlinearSolverSelector.
      nonlinear_solver.setup_jacobian = setup_jacobian;
 
      nonlinear_solver.solve_with_jacobian = solve_with_jacobian;
 
      nonlinear_solver.solve(initial_guess_and_solution);
#else
      AssertThrow(
        false, ExcMessage("You do not have Trilinos configured with deal.II"));
#endif
    }
  else if (additional_data.solver_type ==
           AdditionalData::SolverType::petsc_snes)
    {
      // Forward to internal function specializations, which throws for
      // non-supported vector types:
      solve_with_petsc(initial_guess_and_solution);
    }
  else
    {
      const std::string solvers =
#ifdef DEAL_II_WITH_SUNDIALS
        "kinsol\n"
#endif
#ifdef DEAL_II_WITH_TRILINOS
        "NOX\n"
#endif
#ifdef DEAL_II_WITH_PETSC
        "SNES\n"
#endif
        ;
 
      AssertThrow(false,
                  ExcMessage(
                    "Invalid nonlinear solver specified. "
                    "The solvers available in your installation are:\n" +
                    solvers));
    }
}
 
DEAL_II_NAMESPACE_CLOSE