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 at
//    the top level of the deal.II distribution.
//
//---------------------------------------------------------------

#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/logstream.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_memory.h>

#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 <boost/signals2.hpp>
#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;
  };

}


DEAL_II_NAMESPACE_CLOSE
#endif


#endif