time_stepping.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 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_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() override = 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) override;

    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) override;

    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) override;

    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 override;

  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) override;

    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) override;

    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 override;

  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() override
    {
      free_memory();
    }

    void
    free_memory();

    void
    initialize(const runge_kutta_method method) override;

    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) override;

    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 override;

  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;
  };
} // namespace TimeStepping

DEAL_II_NAMESPACE_CLOSE

#endif