time_stepping.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 2016 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/std_cxx11/function.h>

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



  enum embedded_runge_kutta_time_step { DELTA_T, MIN_DELTA_T, MAX_DELTA_T };



  template <typename VectorType>
  class TimeStepping
  {
  public:
    virtual ~TimeStepping() {}

    virtual double evolve_one_time_step
    (std::vector<std_cxx11::function<VectorType (const double, const VectorType &)> >               &F,
     std::vector<std_cxx11::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() {}

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

    virtual double evolve_one_time_step
    (std_cxx11::function<VectorType (const double, const VectorType &)>               f,
     std_cxx11::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() {}

    ExplicitRungeKutta(runge_kutta_method method);

    void initialize(runge_kutta_method method);

    double evolve_one_time_step
    (std_cxx11::function<VectorType (const double, const VectorType &)>               f,
     std_cxx11::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
    (std_cxx11::function<VectorType (const double, const VectorType &)> f,
     double                                                             t,
     double                                                             delta_t,
     VectorType                                                         &y);

    struct Status : public TimeStepping<VectorType>::Status
    {
      runge_kutta_method method;
    };

    const Status &get_status() const;

  private:
    void compute_stages
    (std_cxx11::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() {}

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

    void initialize(runge_kutta_method method);

    double evolve_one_time_step
    (std_cxx11::function<VectorType (const double, const VectorType &)>               f,
     std_cxx11::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(unsigned int max_it, double tolerance);

    struct Status : public TimeStepping<VectorType>::Status
    {
      runge_kutta_method method;
      unsigned int       n_iterations;
      double             norm_residual;
    };

    const Status &get_status() const;

  private:
    void compute_stages
    (std_cxx11::function<VectorType (const double, const VectorType &)>               f,
     std_cxx11::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(std_cxx11::function<void (const VectorType &,VectorType &)> get_residual,
                      std_cxx11::function<VectorType (const VectorType &)>        id_minus_tau_J_inverse,
                      VectorType                                                  &y);

    void compute_residual(std_cxx11::function<VectorType (const double, const VectorType &)> f,
                          double                                                             t,
                          double                                                             delta_t,
                          const VectorType                                                   &old_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() {}

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

    ~EmbeddedExplicitRungeKutta()
    {
      free_memory();
    }

    void free_memory();

    void initialize(runge_kutta_method method);

    double evolve_one_time_step
    (std_cxx11::function<VectorType (const double, const VectorType &)>               f,
     std_cxx11::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
    (std_cxx11::function<VectorType (const double, const VectorType &)> f,
     double                                                             t,
     double                                                             delta_t,
     VectorType &y);

    void set_time_adaptation_parameters(double coarsen_param,
                                        double refine_param,
                                        double min_delta,
                                        double max_delta,
                                        double refine_tol,
                                        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(std_cxx11::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