kinsol.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2017 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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/mpi_stub.h>
#  include <deal.II/base/parameter_handler.h>
 
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_memory.h>
 
#  include <deal.II/sundials/sundials_types.h>
 
#  include <boost/signals2.hpp>
 
#  include <kinsol/kinsol.h>
#  include <nvector/nvector_serial.h>
#  include <sundials/sundials_math.h>
 
#  include <exception>
#  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,
      };
 
 
      enum OrthogonalizationStrategy
      {
        modified_gram_schmidt = 0,
 
        inverse_compact = 1,
 
        classical_gram_schmidt = 2,
 
        delayed_classical_gram_schmidt = 3,
      };
 
      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,
                     const OrthogonalizationStrategy
                       anderson_qr_orthogonalization = modified_gram_schmidt);
 
      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;
 
      OrthogonalizationStrategy anderson_qr_orthogonalization;
    };
 
    KINSOL(const AdditionalData &data = AdditionalData());
 
    KINSOL(const AdditionalData &data, const MPI_Comm mpi_comm);
 
    ~KINSOL();
 
    unsigned int
    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 &src, VectorType &dst)>
      iteration_function;
 
    std::function<void(const VectorType &current_u,
                       const VectorType &current_f)>
      setup_jacobian;
 
    std::function<
      void(const VectorType &rhs, VectorType &dst, const double tolerance)>
      solve_with_jacobian;
 
    std::function<VectorType &()> get_solution_scaling;
 
    std::function<VectorType &()> get_function_scaling;
 
 
    std::function<void(void *kinsol_mem)> custom_setup;
 
    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;
 
    MPI_Comm mpi_communicator;
 
    void *kinsol_mem;
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
    SUNContext kinsol_ctx;
#  endif
 
 
    GrowingVectorMemory<VectorType> mem;
 
    mutable std::exception_ptr pending_exception;
  };
 
} // namespace SUNDIALS
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif
 
#endif