time_stepping.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2014 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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,
    SSP_THIRD_ORDER,
    RK_CLASSIC_FOURTH_ORDER,
    LOW_STORAGE_RK_STAGE3_ORDER3,
    LOW_STORAGE_RK_STAGE5_ORDER4,
    LOW_STORAGE_RK_STAGE7_ORDER4,
    LOW_STORAGE_RK_STAGE9_ORDER5,
    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 LowStorageRungeKutta : public RungeKutta<VectorType>
  {
  public:
    using RungeKutta<VectorType>::evolve_one_time_step;
 
    LowStorageRungeKutta() = default;
 
    LowStorageRungeKutta(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 &solution,
      VectorType &vec_ri,
      VectorType &vec_ki);
 
    void
    get_coefficients(std::vector<double> &a,
                     std::vector<double> &b,
                     std::vector<double> &c) const;
 
    struct Status : public TimeStepping<VectorType>::Status
    {
      Status()
        : method(invalid)
      {}
 
      runge_kutta_method method;
    };
 
    const Status &
    get_status() const override;
 
  private:
    void
    compute_one_stage(
      const std::function<VectorType(const double, const VectorType &)> &f,
      const double                                                       t,
      const double      factor_solution,
      const double      factor_ai,
      const VectorType &current_ri,
      VectorType       &vec_ki,
      VectorType       &solution,
      VectorType       &next_ri) 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;
 
    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
    {
      Status()
        : method(invalid)
      {}
 
      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 = false;
 
    std::vector<double> b1;
 
    std::vector<double> b2;
 
    VectorType *last_stage = nullptr;
 
    Status status;
  };
} // namespace TimeStepping
 
DEAL_II_NAMESPACE_CLOSE
 
#endif