arkode.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_arkode_h
#define dealii_sundials_arkode_h

#include <deal.II/base/config.h>
#include <deal.II/base/mpi.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>
#ifdef DEAL_II_WITH_PETSC
#  include <deal.II/lac/petsc_parallel_block_vector.h>
#  include <deal.II/lac/petsc_parallel_vector.h>
#endif
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_memory.h>


#include <arkode/arkode.h>
#include <arkode/arkode_impl.h>
#include <nvector/nvector_serial.h>
#ifdef DEAL_II_WITH_MPI
#  include <nvector/nvector_parallel.h>
#endif
#include <sundials/sundials_math.h>
#include <sundials/sundials_types.h>

#include <boost/signals2.hpp>
#include <memory>


DEAL_II_NAMESPACE_OPEN

// Shorthand notation for ARKODE error codes.
#define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))

namespace SUNDIALS
{
  template<typename VectorType=Vector<double> >
  class ARKode
  {
  public:

    class AdditionalData
    {
    public:
      AdditionalData(
        // Initial parameters
        const double &initial_time = 0.0,
        const double &final_time = 1.0,
        const double &initial_step_size = 1e-2,
        const double &output_period = 1e-1,
        // Running parameters
        const double &minimum_step_size = 1e-6,
        const unsigned int &maximum_order = 5,
        const unsigned int &maximum_non_linear_iterations = 10,
        const bool implicit_function_is_linear = false,
        const bool implicit_function_is_time_independent = false,
        // Error parameters
        const double &absolute_tolerance = 1e-6,
        const double &relative_tolerance = 1e-5) :
        initial_time(initial_time),
        final_time(final_time),
        initial_step_size(initial_step_size),
        minimum_step_size(minimum_step_size),
        absolute_tolerance(absolute_tolerance),
        relative_tolerance(relative_tolerance),
        maximum_order(maximum_order),
        output_period(output_period),
        maximum_non_linear_iterations(maximum_non_linear_iterations),
        implicit_function_is_linear(implicit_function_is_linear),
        implicit_function_is_time_independent(implicit_function_is_time_independent)
      {}

      void add_parameters(ParameterHandler &prm)
      {
        prm.add_parameter("Initial time", initial_time);
        prm.add_parameter("Final time", final_time);
        prm.add_parameter("Time interval between each output", output_period);

        prm.enter_subsection("Running parameters");
        prm.add_parameter("Initial step size",initial_step_size);
        prm.add_parameter("Minimum step size", minimum_step_size);
        prm.add_parameter("Maximum order of ARK", maximum_order);
        prm.add_parameter("Maximum number of nonlinear iterations", maximum_non_linear_iterations);
        prm.add_parameter("Implicit function is linear", implicit_function_is_linear);
        prm.add_parameter("Implicit function is time independent", implicit_function_is_time_independent);
        prm.leave_subsection();

        prm.enter_subsection("Error control");
        prm.add_parameter("Absolute error tolerance", absolute_tolerance);
        prm.add_parameter("Relative error tolerance", relative_tolerance);
        prm.leave_subsection();
      }

      double initial_time;

      double final_time;

      double initial_step_size;

      double minimum_step_size;

      double absolute_tolerance;

      double relative_tolerance;

      unsigned int maximum_order;

      double output_period;

      unsigned int maximum_non_linear_iterations;

      bool implicit_function_is_linear;

      bool implicit_function_is_time_independent;
    };

    ARKode(const AdditionalData &data=AdditionalData(),
           const MPI_Comm mpi_comm = MPI_COMM_WORLD);

    ~ARKode();

    unsigned int solve_ode(VectorType &solution);

    void reset(const double &t,
               const double &h,
               const VectorType &y);

    std::function<void(VectorType &)> reinit_vector;

    std::function<int(const double t,
                      const VectorType &y,
                      VectorType &explicit_f)> explicit_function;

    std::function<int(const double t,
                      const VectorType &y,
                      VectorType &res)> implicit_function;

    std::function<int(const int convfail,
                      const double t,
                      const double gamma,
                      const VectorType &ypred,
                      const VectorType &fpred,
                      bool &j_is_current)> setup_jacobian;

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


    std::function<int(const double t)> setup_mass;

    std::function<int(const VectorType &rhs, VectorType &dst)> solve_mass_system;

    std::function<void(const double t,
                       const VectorType &sol,
                       const unsigned int step_number)> output_step;

    std::function<bool (const double t,
                        VectorType &sol)> solver_should_restart;

    std::function<VectorType&()> get_local_tolerances;

    DeclException1(ExcARKodeError, int, << "One of the SUNDIALS ARKode 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 *arkode_mem;

    N_Vector yy;

    N_Vector abs_tolls;

    MPI_Comm communicator;

    GrowingVectorMemory<VectorType> mem;

#ifdef DEAL_II_WITH_PETSC
#ifdef PETSC_USE_COMPLEX
    static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
                  "Sundials does not support complex scalar types, "
                  "but PETSc is configured to use a complex scalar type!");

    static_assert(!std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
                  "Sundials does not support complex scalar types, "
                  "but PETSc is configured to use a complex scalar type!");
#endif // PETSC_USE_COMPLEX
#endif // DEAL_II_WITH_PETSC
  };

}

DEAL_II_NAMESPACE_CLOSE
#endif


#endif