Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
time_stepping.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2018 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.md at
12 // the top level directory of deal.II.
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 
22 #include <deal.II/base/signaling_nan.h>
23 
24 #include <functional>
25 #include <vector>
26 
27 DEAL_II_NAMESPACE_OPEN
28 
36 namespace TimeStepping
37 {
58  {
59  FORWARD_EULER,
60  RK_THIRD_ORDER,
61  RK_CLASSIC_FOURTH_ORDER,
62  BACKWARD_EULER,
63  IMPLICIT_MIDPOINT,
64  CRANK_NICOLSON,
65  SDIRK_TWO_STAGES,
66  HEUN_EULER,
67  BOGACKI_SHAMPINE,
68  DOPRI,
69  FEHLBERG,
70  CASH_KARP,
71  invalid
72  };
73 
74 
75 
83  {
84  DELTA_T,
85  MIN_DELTA_T,
86  MAX_DELTA_T
87  };
88 
89 
90 
95  template <typename VectorType>
97  {
98  public:
102  virtual ~TimeStepping() = default;
103 
114  virtual double
115  evolve_one_time_step(
116  std::vector<std::function<VectorType(const double, const VectorType &)>>
117  & F,
118  std::vector<std::function<
119  VectorType(const double, const double, const VectorType &)>> &J_inverse,
120  double t,
121  double delta_t,
122  VectorType & y) = 0;
123 
127  struct Status
128  {};
129 
133  virtual const Status &
134  get_status() const = 0;
135  };
136 
137 
138 
145  template <typename VectorType>
146  class RungeKutta : public TimeStepping<VectorType>
147  {
148  public:
152  virtual ~RungeKutta() override = default;
153 
157  virtual void
158  initialize(const runge_kutta_method method) = 0;
159 
171  double
173  std::vector<std::function<VectorType(const double, const VectorType &)>>
174  & F,
175  std::vector<std::function<
176  VectorType(const double, const double, const VectorType &)>> &J_inverse,
177  double t,
178  double delta_t,
179  VectorType &y) override;
180 
192  virtual double
194  const std::function<VectorType(const double, const VectorType &)> &f,
195  const std::function<
196  VectorType(const double, const double, const VectorType &)>
197  & id_minus_tau_J_inverse,
198  double t,
199  double delta_t,
200  VectorType &y) = 0;
201 
202  protected:
206  unsigned int n_stages;
207 
211  std::vector<double> b;
212 
216  std::vector<double> c;
217 
221  std::vector<std::vector<double>> a;
222  };
223 
224 
225 
230  template <typename VectorType>
231  class ExplicitRungeKutta : public RungeKutta<VectorType>
232  {
233  public:
235 
241  ExplicitRungeKutta() = default;
242 
247 
251  void
252  initialize(const runge_kutta_method method) override;
253 
265  double
267  const std::function<VectorType(const double, const VectorType &)> &f,
268  const std::function<
269  VectorType(const double, const double, const VectorType &)>
270  & id_minus_tau_J_inverse,
271  double t,
272  double delta_t,
273  VectorType &y) override;
274 
282  double
284  const std::function<VectorType(const double, const VectorType &)> &f,
285  double t,
286  double delta_t,
287  VectorType &y);
288 
292  struct Status : public TimeStepping<VectorType>::Status
293  {
294  Status()
295  : method(invalid)
296  {}
297 
298  runge_kutta_method method;
299  };
300 
304  const Status &
305  get_status() const override;
306 
307  private:
311  void
313  const std::function<VectorType(const double, const VectorType &)> &f,
314  const double t,
315  const double delta_t,
316  const VectorType & y,
317  std::vector<VectorType> &f_stages) const;
318 
323  };
324 
325 
326 
331  template <typename VectorType>
332  class ImplicitRungeKutta : public RungeKutta<VectorType>
333  {
334  public:
336 
342  ImplicitRungeKutta() = default;
343 
350  const unsigned int max_it = 100,
351  const double tolerance = 1e-6);
352 
356  void
357  initialize(const runge_kutta_method method) override;
358 
370  double
372  const std::function<VectorType(const double, const VectorType &)> &f,
373  const std::function<
374  VectorType(const double, const double, const VectorType &)>
375  & id_minus_tau_J_inverse,
376  double t,
377  double delta_t,
378  VectorType &y) override;
379 
384  void
385  set_newton_solver_parameters(const unsigned int max_it,
386  const double tolerance);
387 
392  struct Status : public TimeStepping<VectorType>::Status
393  {
394  Status()
395  : method(invalid)
396  , n_iterations(numbers::invalid_unsigned_int)
397  , norm_residual(numbers::signaling_nan<double>())
398  {}
399 
400  runge_kutta_method method;
401  unsigned int n_iterations;
402  double norm_residual;
403  };
404 
408  const Status &
409  get_status() const override;
410 
411  private:
415  void
417  const std::function<VectorType(const double, const VectorType &)> &f,
418  const std::function<
419  VectorType(const double, const double, const VectorType &)>
420  & id_minus_tau_J_inverse,
421  double t,
422  double delta_t,
423  VectorType & y,
424  std::vector<VectorType> &f_stages);
425 
429  void
430  newton_solve(
431  const std::function<void(const VectorType &, VectorType &)> &get_residual,
432  const std::function<VectorType(const VectorType &)>
433  & id_minus_tau_J_inverse,
434  VectorType &y);
435 
439  void
441  const std::function<VectorType(const double, const VectorType &)> &f,
442  double t,
443  double delta_t,
444  const VectorType &new_y,
445  const VectorType &y,
446  VectorType & tendency,
447  VectorType & residual) const;
448 
455 
459  unsigned int max_it;
460 
464  double tolerance;
465 
470  };
471 
472 
473 
478  template <typename VectorType>
479  class EmbeddedExplicitRungeKutta : public RungeKutta<VectorType>
480  {
481  public:
483 
489  EmbeddedExplicitRungeKutta() = default;
490 
496  const double coarsen_param = 1.2,
497  const double refine_param = 0.8,
498  const double min_delta = 1e-14,
499  const double max_delta = 1e100,
500  const double refine_tol = 1e-8,
501  const double coarsen_tol = 1e-12);
502 
507  {
508  free_memory();
509  }
510 
514  void
515  free_memory();
516 
520  void
521  initialize(const runge_kutta_method method) override;
522 
534  double
536  const std::function<VectorType(const double, const VectorType &)> &f,
537  const std::function<
538  VectorType(const double, const double, const VectorType &)>
539  & id_minus_tau_J_inverse,
540  double t,
541  double delta_t,
542  VectorType &y) override;
543 
551  double
553  const std::function<VectorType(const double, const VectorType &)> &f,
554  double t,
555  double delta_t,
556  VectorType &y);
557 
561  void
563  const double refine_param,
564  const double min_delta,
565  const double max_delta,
566  const double refine_tol,
567  const double coarsen_tol);
568 
575  struct Status : public TimeStepping<VectorType>::Status
576  {
577  runge_kutta_method method;
578  embedded_runge_kutta_time_step exit_delta_t;
579  unsigned int n_iterations;
580  double delta_t_guess;
581  double error_norm;
582  };
583 
587  const Status &
588  get_status() const override;
589 
590  private:
594  void
596  const std::function<VectorType(const double, const VectorType &)> &f,
597  const double t,
598  const double delta_t,
599  const VectorType & y,
600  std::vector<VectorType> &f_stages);
601 
607 
612  double refine_param;
613 
617  double min_delta_t;
618 
622  double max_delta_t;
623 
628  double refine_tol;
629 
634  double coarsen_tol;
635 
641 
645  std::vector<double> b1;
646 
650  std::vector<double> b2;
651 
656  VectorType *last_stage;
657 
662  };
663 } // namespace TimeStepping
664 
665 DEAL_II_NAMESPACE_CLOSE
666 
667 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void set_newton_solver_parameters(const unsigned int max_it, const double tolerance)
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
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
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
void initialize(const runge_kutta_method method) override
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)
virtual ~RungeKutta() override=default
const Status & get_status() const override
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_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 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, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
std::vector< double > b
std::vector< double > c
std::vector< std::vector< double > > a
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
const Status & get_status() const override
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)
virtual void initialize(const runge_kutta_method method)=0
embedded_runge_kutta_time_step
Definition: time_stepping.h:82
const Status & get_status() const override
void initialize(const runge_kutta_method method) override
void initialize(const runge_kutta_method method) override