time_stepping.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 2017 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_time_stepping_h
#define dealii_time_stepping_h


#include <deal.II/base/config.h>
#include <deal.II/base/signaling_nan.h>

#include <functional>
#include <vector>

DEAL_II_NAMESPACE_OPEN

namespace TimeStepping
{
  enum runge_kutta_method { FORWARD_EULER, RK_THIRD_ORDER, RK_CLASSIC_FOURTH_ORDER,
                            BACKWARD_EULER, IMPLICIT_MIDPOINT, CRANK_NICOLSON,
                            SDIRK_TWO_STAGES, HEUN_EULER, BOGACKI_SHAMPINE, DOPRI,
                            FEHLBERG, CASH_KARP,
                            invalid
                          };



  enum embedded_runge_kutta_time_step { DELTA_T, MIN_DELTA_T, MAX_DELTA_T };



  template <typename VectorType>
  class TimeStepping
  {
  public:
    virtual ~TimeStepping() = default;

    virtual double evolve_one_time_step
    (std::vector<std::function<VectorType (const double, const VectorType &)> >               &F,
     std::vector<std::function<VectorType (const double, const double, const VectorType &)> > &J_inverse,
     double                                                                                   t,
     double                                                                                   delta_t,
     VectorType                                                                               &y) = 0;

    struct Status {};

    virtual const Status &get_status() const = 0;
  };



  template <typename VectorType>
  class RungeKutta : public TimeStepping<VectorType>
  {
  public:
    virtual ~RungeKutta() = default;

    virtual void initialize(const runge_kutta_method method) = 0;

    double evolve_one_time_step
    (std::vector<std::function<VectorType (const double, const VectorType &)> >               &F,
     std::vector<std::function<VectorType (const double, const double, const VectorType &)> > &J_inverse,
     double                                                                                   t,
     double                                                                                   delta_t,
     VectorType                                                                               &y);

    virtual double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)>               &f,
     const std::function<VectorType (const double, const double, const VectorType &)> &id_minus_tau_J_inverse,
     double                                                                     t,
     double                                                                     delta_t,
     VectorType                                                                 &y) = 0;

  protected:
    unsigned int n_stages;

    std::vector<double> b;

    std::vector<double> c;

    std::vector<std::vector<double> > a;
  };



  template <typename VectorType>
  class ExplicitRungeKutta : public RungeKutta<VectorType>
  {
  public:
    using RungeKutta<VectorType>::evolve_one_time_step;

    ExplicitRungeKutta() = default;

    ExplicitRungeKutta(const runge_kutta_method method);

    void initialize(const runge_kutta_method method);

    double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)>                &f,
     const std::function<VectorType (const double, const double, const VectorType &)>  &id_minus_tau_J_inverse,
     double                                                                      t,
     double                                                                      delta_t,
     VectorType                                                                  &y);

    double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)> &f,
     double                                                       t,
     double                                                       delta_t,
     VectorType                                                   &y);

    struct Status : public TimeStepping<VectorType>::Status
    {
      Status ()
        :
        method (invalid)
      {}

      runge_kutta_method method;
    };

    const Status &get_status() const;

  private:
    void compute_stages
    (const std::function<VectorType (const double, const VectorType &)> &f,
     const double                                                 t,
     const double                                                 delta_t,
     const VectorType                                             &y,
     std::vector<VectorType>                                      &f_stages) const;

    Status status;
  };



  template <typename VectorType>
  class ImplicitRungeKutta : public RungeKutta<VectorType>
  {
  public:
    using RungeKutta<VectorType>::evolve_one_time_step;

    ImplicitRungeKutta() = default;

    ImplicitRungeKutta(const runge_kutta_method method,
                       const unsigned int       max_it    = 100,
                       const double             tolerance = 1e-6);

    void initialize(const runge_kutta_method method);

    double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)>               &f,
     const std::function<VectorType (const double, const double, const VectorType &)> &id_minus_tau_J_inverse,
     double                                                                           t,
     double                                                                           delta_t,
     VectorType                                                                       &y);

    void set_newton_solver_parameters(const unsigned int max_it,
                                      const double       tolerance);

    struct Status : public TimeStepping<VectorType>::Status
    {
      Status ()
        :
        method (invalid),
        n_iterations (numbers::invalid_unsigned_int),
        norm_residual (numbers::signaling_nan<double>())
      {}

      runge_kutta_method method;
      unsigned int       n_iterations;
      double             norm_residual;
    };

    const Status &get_status() const;

  private:
    void compute_stages
    (const std::function<VectorType (const double, const VectorType &)>               &f,
     const std::function<VectorType (const double, const double, const VectorType &)> &id_minus_tau_J_inverse,
     double                                                                           t,
     double                                                                           delta_t,
     VectorType                                                                       &y,
     std::vector<VectorType> &f_stages);

    void newton_solve(const std::function<void (const VectorType &,VectorType &)> &get_residual,
                      const std::function<VectorType (const VectorType &)>        &id_minus_tau_J_inverse,
                      VectorType                                                  &y);

    void compute_residual(const std::function<VectorType (const double, const VectorType &)> &f,
                          double                                                             t,
                          double                                                             delta_t,
                          const VectorType                                                   &new_y,
                          const VectorType                                                   &y,
                          VectorType                                                         &tendency,
                          VectorType                                                         &residual) const;

    bool skip_linear_combi;

    unsigned int max_it;

    double tolerance;

    Status status;
  };



  template <typename VectorType>
  class EmbeddedExplicitRungeKutta : public RungeKutta<VectorType>
  {
  public:
    using RungeKutta<VectorType>::evolve_one_time_step;

    EmbeddedExplicitRungeKutta() = default;

    EmbeddedExplicitRungeKutta(const runge_kutta_method method,
                               const double coarsen_param = 1.2,
                               const double refine_param = 0.8,
                               const double min_delta = 1e-14,
                               const double max_delta = 1e100,
                               const double refine_tol = 1e-8,
                               const double coarsen_tol = 1e-12);

    ~EmbeddedExplicitRungeKutta()
    {
      free_memory();
    }

    void free_memory();

    void initialize(const runge_kutta_method method);

    double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)>               &f,
     const std::function<VectorType (const double, const double, const VectorType &)> &id_minus_tau_J_inverse,
     double                                                                     t,
     double                                                                     delta_t,
     VectorType &y);

    double evolve_one_time_step
    (const std::function<VectorType (const double, const VectorType &)> &f,
     double                                                       t,
     double                                                       delta_t,
     VectorType &y);

    void set_time_adaptation_parameters(const double coarsen_param,
                                        const double refine_param,
                                        const double min_delta,
                                        const double max_delta,
                                        const double refine_tol,
                                        const double coarsen_tol);

    struct Status : public TimeStepping<VectorType>::Status
    {
      runge_kutta_method method;
      embedded_runge_kutta_time_step exit_delta_t;
      unsigned int n_iterations;
      double delta_t_guess;
      double error_norm;
    };

    const Status &get_status() const;

  private:
    void compute_stages(const std::function<VectorType (const double, const VectorType &)> &f,
                        const double                                                 t,
                        const double                                                 delta_t,
                        const VectorType                                             &y,
                        std::vector<VectorType>                                      &f_stages);

    double coarsen_param;

    double refine_param;

    double min_delta_t;

    double max_delta_t;

    double refine_tol;

    double coarsen_tol;

    bool last_same_as_first;

    std::vector<double> b1;

    std::vector<double> b2;

    VectorType *last_stage;

    Status status;
  };
}

DEAL_II_NAMESPACE_CLOSE

#endif