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

#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/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>
#  ifdef DEAL_II_WITH_PETSC
#    include <deal.II/lac/petsc_block_vector.h>
#    include <deal.II/lac/petsc_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 <boost/signals2.hpp>

#  include <sundials/sundials_math.h>
#  include <sundials/sundials_types.h>

#  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
  };

} // namespace SUNDIALS

DEAL_II_NAMESPACE_CLOSE
#endif


#endif