kinsol.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    Copyright (C) 2017 - 2021 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.
//
//---------------------------------------------------------------
 
#ifndef dealii_sundials_kinsol_h
#define dealii_sundials_kinsol_h
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_SUNDIALS
 
 
#  include <deal.II/base/conditional_ostream.h>
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/logstream.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/base/parameter_handler.h>
 
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_memory.h>
 
#  include <boost/signals2.hpp>
 
#  include <kinsol/kinsol.h>
#  if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
#    include <kinsol/kinsol_impl.h>
#  endif
#  include <nvector/nvector_serial.h>
#  include <sundials/sundials_math.h>
#  include <sundials/sundials_types.h>
 
#  include <memory>
 
 
DEAL_II_NAMESPACE_OPEN
 
// Shorthand notation for KINSOL error codes.
#  define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code))
 
namespace SUNDIALS
{
  template <typename VectorType = Vector<double>>
  class KINSOL
  {
  public:
    class AdditionalData
    {
    public:
      enum SolutionStrategy
      {
        newton = KIN_NONE,
        linesearch = KIN_LINESEARCH,
        fixed_point = KIN_FP,
        picard = KIN_PICARD,
      };
 
      AdditionalData(const SolutionStrategy &strategy = linesearch,
                     const unsigned int maximum_non_linear_iterations = 200,
                     const double       function_tolerance            = 0.0,
                     const double       step_tolerance                = 0.0,
                     const bool         no_init_setup                 = false,
                     const unsigned int maximum_setup_calls           = 0,
                     const double       maximum_newton_step           = 0.0,
                     const double       dq_relative_error             = 0.0,
                     const unsigned int maximum_beta_failures         = 0,
                     const unsigned int anderson_subspace_size        = 0);
 
      void
      add_parameters(ParameterHandler &prm);
 
      SolutionStrategy strategy;
 
      unsigned int maximum_non_linear_iterations;
 
      double function_tolerance;
 
      double step_tolerance;
 
      bool no_init_setup;
 
      unsigned int maximum_setup_calls;
 
      double maximum_newton_step;
 
      double dq_relative_error;
 
      unsigned int maximum_beta_failures;
 
      unsigned int anderson_subspace_size;
    };
 
    KINSOL(const AdditionalData &data     = AdditionalData(),
           const MPI_Comm &      mpi_comm = MPI_COMM_WORLD);
 
    ~KINSOL();
 
    unsigned int
    solve(VectorType &initial_guess_and_solution);
 
    std::function<void(VectorType &)> reinit_vector;
 
    std::function<int(const VectorType &src, VectorType &dst)> residual;
 
    std::function<int(const VectorType &src, VectorType &dst)>
      iteration_function;
 
    std::function<int(const VectorType &current_u, const VectorType &current_f)>
      setup_jacobian;
 
    DEAL_II_DEPRECATED_EARLY
    std::function<int(const VectorType &ycur,
                      const VectorType &fcur,
                      const VectorType &rhs,
                      VectorType &      dst)>
      solve_jacobian_system;
 
    std::function<
      int(const VectorType &rhs, VectorType &dst, const double tolerance)>
      solve_with_jacobian;
 
    std::function<VectorType &()> get_solution_scaling;
 
    std::function<VectorType &()> get_function_scaling;
 
    DeclException1(ExcKINSOLError,
                   int,
                   << "One of the SUNDIALS KINSOL internal functions "
                   << "returned a negative error code: " << arg1
                   << ". Please consult SUNDIALS manual.");
 
 
  private:
    DeclException1(ExcFunctionNotProvided,
                   std::string,
                   << "Please provide an implementation for the function \""
                   << arg1 << "\"");
 
    void
    set_functions_to_trigger_an_assert();
 
    AdditionalData data;
 
    void *kinsol_mem;
 
    N_Vector solution;
 
    N_Vector u_scale;
 
    N_Vector f_scale;
 
    MPI_Comm communicator;
 
    GrowingVectorMemory<VectorType> mem;
  };
 
} // namespace SUNDIALS
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif
 
#endif