Reference documentation for deal.II version 8.4.1
time_stepping.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2016 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__time_stepping_h
17 #define dealii__time_stepping_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/std_cxx11/function.h>
22 
23 #include <vector>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
34 namespace TimeStepping
35 {
55  enum runge_kutta_method { FORWARD_EULER, RK_THIRD_ORDER, RK_CLASSIC_FOURTH_ORDER,
56  BACKWARD_EULER, IMPLICIT_MIDPOINT, CRANK_NICOLSON,
57  SDIRK_TWO_STAGES, HEUN_EULER, BOGACKI_SHAMPINE, DOPRI,
58  FEHLBERG, CASH_KARP
59  };
60 
61 
62 
69  enum embedded_runge_kutta_time_step { DELTA_T, MIN_DELTA_T, MAX_DELTA_T };
70 
71 
72 
77  template <typename VectorType>
79  {
80  public:
84  virtual ~TimeStepping() {}
85 
96  virtual double evolve_one_time_step
97  (std::vector<std_cxx11::function<VectorType (const double, const VectorType &)> > &F,
98  std::vector<std_cxx11::function<VectorType (const double, const double, const VectorType &)> > &J_inverse,
99  double t,
100  double delta_t,
101  VectorType &y) = 0;
102 
106  struct Status {};
107 
111  virtual const Status &get_status() const = 0;
112  };
113 
114 
115 
122  template <typename VectorType>
123  class RungeKutta : public TimeStepping<VectorType>
124  {
125  public:
129  virtual ~RungeKutta() {}
130 
134  virtual void initialize(runge_kutta_method method) = 0;
146  double evolve_one_time_step
147  (std::vector<std_cxx11::function<VectorType (const double, const VectorType &)> > &F,
148  std::vector<std_cxx11::function<VectorType (const double, const double, const VectorType &)> > &J_inverse,
149  double t,
150  double delta_t,
151  VectorType &y);
152 
164  virtual double evolve_one_time_step
165  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
166  std_cxx11::function<VectorType (const double, const double, const VectorType &)> id_minus_tau_J_inverse,
167  double t,
168  double delta_t,
169  VectorType &y) = 0;
170 
171  protected:
175  unsigned int n_stages;
176 
180  std::vector<double> b;
181 
185  std::vector<double> c;
186 
190  std::vector<std::vector<double> > a;
191  };
192 
193 
194 
199  template <typename VectorType>
200  class ExplicitRungeKutta : public RungeKutta<VectorType>
201  {
202  public:
204 
210 
215 
219  void initialize(runge_kutta_method method);
220 
232  double evolve_one_time_step
233  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
234  std_cxx11::function<VectorType (const double, const double, const VectorType &)> id_minus_tau_J_inverse,
235  double t,
236  double delta_t,
237  VectorType &y);
238 
246  double evolve_one_time_step
247  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
248  double t,
249  double delta_t,
250  VectorType &y);
251 
255  struct Status : public TimeStepping<VectorType>::Status
256  {
257  runge_kutta_method method;
258  };
259 
263  const Status &get_status() const;
264 
265  private:
269  void compute_stages
270  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
271  const double t,
272  const double delta_t,
273  const VectorType &y,
274  std::vector<VectorType> &f_stages) const;
275 
280  };
281 
282 
283 
288  template <typename VectorType>
289  class ImplicitRungeKutta : public RungeKutta<VectorType>
290  {
291  public:
293 
300 
306  ImplicitRungeKutta(runge_kutta_method method, unsigned int max_it=100, double tolerance=1e-6);
307 
311  void initialize(runge_kutta_method method);
312 
324  double evolve_one_time_step
325  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
326  std_cxx11::function<VectorType (const double, const double, const VectorType &)> id_minus_tau_J_inverse,
327  double t,
328  double delta_t,
329  VectorType &y);
330 
335  void set_newton_solver_parameters(unsigned int max_it, double tolerance);
336 
341  struct Status : public TimeStepping<VectorType>::Status
342  {
343  runge_kutta_method method;
344  unsigned int n_iterations;
345  double norm_residual;
346  };
347 
351  const Status &get_status() const;
352 
353  private:
357  void compute_stages
358  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
359  std_cxx11::function<VectorType (const double, const double, const VectorType &)> id_minus_tau_J_inverse,
360  double t,
361  double delta_t,
362  VectorType &y,
363  std::vector<VectorType> &f_stages);
364 
368  void newton_solve(std_cxx11::function<void (const VectorType &,VectorType &)> get_residual,
369  std_cxx11::function<VectorType (const VectorType &)> id_minus_tau_J_inverse,
370  VectorType &y);
371 
375  void compute_residual(std_cxx11::function<VectorType (const double, const VectorType &)> f,
376  double t,
377  double delta_t,
378  const VectorType &old_y,
379  const VectorType &y,
380  VectorType &tendency,
381  VectorType &residual) const;
382 
389 
393  unsigned int max_it;
394 
398  double tolerance;
399 
404  };
405 
406 
407 
412  template <typename VectorType>
413  class EmbeddedExplicitRungeKutta : public RungeKutta<VectorType>
414  {
415  public:
417 
424 
430  double coarsen_param = 1.2,
431  double refine_param = 0.8,
432  double min_delta = 1e-14,
433  double max_delta = 1e100,
434  double refine_tol = 1e-8,
435  double coarsen_tol = 1e-12);
436 
441  {
442  free_memory();
443  }
444 
448  void free_memory();
449 
453  void initialize(runge_kutta_method method);
454 
466  double evolve_one_time_step
467  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
468  std_cxx11::function<VectorType (const double, const double, const VectorType &)> id_minus_tau_J_inverse,
469  double t,
470  double delta_t,
471  VectorType &y);
472 
480  double evolve_one_time_step
481  (std_cxx11::function<VectorType (const double, const VectorType &)> f,
482  double t,
483  double delta_t,
484  VectorType &y);
485 
490  double refine_param,
491  double min_delta,
492  double max_delta,
493  double refine_tol,
494  double coarsen_tol);
495 
502  struct Status : public TimeStepping<VectorType>::Status
503  {
504  runge_kutta_method method;
505  embedded_runge_kutta_time_step exit_delta_t;
506  unsigned int n_iterations;
507  double delta_t_guess;
508  double error_norm;
509  };
510 
514  const Status &get_status() const;
515 
516  private:
520  void compute_stages(std_cxx11::function<VectorType (const double, const VectorType &)> f,
521  const double t,
522  const double delta_t,
523  const VectorType &y,
524  std::vector<VectorType> &f_stages);
525 
531 
536  double refine_param;
537 
541  double min_delta_t;
542 
546  double max_delta_t;
547 
552  double refine_tol;
553 
558  double coarsen_tol;
559 
565 
569  std::vector<double> b1;
570 
574  std::vector<double> b2;
575 
580  VectorType *last_stage;
581 
586  };
587 }
588 
589 DEAL_II_NAMESPACE_CLOSE
590 
591 #endif
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
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)
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 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)
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)
virtual void initialize(runge_kutta_method method)=0
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 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)
void set_time_adaptation_parameters(double coarsen_param, double refine_param, double min_delta, double max_delta, double refine_tol, double coarsen_tol)
const Status & get_status() const
void initialize(runge_kutta_method method)
std::vector< double > b
const Status & get_status() const
std::vector< double > c
std::vector< std::vector< double > > a
void initialize(runge_kutta_method method)
void set_newton_solver_parameters(unsigned int max_it, double tolerance)
embedded_runge_kutta_time_step
Definition: time_stepping.h:69
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)
void initialize(runge_kutta_method method)
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