kinsol.h

//-----------------------------------------------------------
//
//    Copyright (C) 2017 - 2018 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>
#  include <kinsol/kinsol_impl.h>
#  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(
        // Global parameters
        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)
        : strategy(strategy)
        , maximum_non_linear_iterations(maximum_non_linear_iterations)
        , function_tolerance(function_tolerance)
        , step_tolerance(step_tolerance)
        , no_init_setup(no_init_setup)
        , maximum_setup_calls(maximum_setup_calls)
        , maximum_newton_step(maximum_newton_step)
        , dq_relative_error(dq_relative_error)
        , maximum_beta_failures(maximum_beta_failures)
        , anderson_subspace_size(anderson_subspace_size)
      {}

      void
      add_parameters(ParameterHandler &prm)
      {
        static std::string strategy_str("newton");
        prm.add_parameter("Solution strategy",
                          strategy_str,
                          "Choose among newton|linesearch|fixed_point|picard",
                          Patterns::Selection(
                            "newton|linesearch|fixed_point|picard"));
        prm.add_action("Solution strategy", [&](const std::string &value) {
          if (value == "newton")
            strategy = newton;
          else if (value == "linesearch")
            strategy = linesearch;
          else if (value == "fixed_point")
            strategy = fixed_point;
          else if (value == "picard")
            strategy = picard;
          else
            Assert(false, ExcInternalError());
        });
        prm.add_parameter("Maximum number of nonlinear iterations",
                          maximum_non_linear_iterations);
        prm.add_parameter("Function norm stopping tolerance",
                          function_tolerance);
        prm.add_parameter("Scaled step stopping tolerance", step_tolerance);

        prm.enter_subsection("Newton parameters");
        prm.add_parameter("No initial matrix setup", no_init_setup);
        prm.add_parameter("Maximum iterations without matrix setup",
                          maximum_setup_calls);
        prm.add_parameter("Maximum allowable scaled length of the Newton step",
                          maximum_newton_step);
        prm.add_parameter("Relative error for different quotient computation",
                          dq_relative_error);
        prm.leave_subsection();

        prm.enter_subsection("Linesearch parameters");
        prm.add_parameter("Maximum number of beta-condition failures",
                          maximum_beta_failures);
        prm.leave_subsection();


        prm.enter_subsection("Fixed point and Picard parameters");
        prm.add_parameter("Anderson acceleration subspace size",
                          anderson_subspace_size);
        prm.leave_subsection();
      }

      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;

    std::function<int(const VectorType &ycur,
                      const VectorType &fcur,
                      const VectorType &rhs,
                      VectorType &      dst)>
      solve_jacobian_system;

    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