Reference documentation for deal.II version 9.1.0-pre
time_stepping.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2017 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 {
57  enum runge_kutta_method
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 
82  enum embedded_runge_kutta_time_step
83  {
84  DELTA_T,
85  MIN_DELTA_T,
86  MAX_DELTA_T
87  };
88 
89 
90 
95  template <typename VectorType>
96  class TimeStepping
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
172  evolve_one_time_step(
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
193  evolve_one_time_step(
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:
234  using RungeKutta<VectorType>::evolve_one_time_step;
235 
241  ExplicitRungeKutta() = default;
242 
246  ExplicitRungeKutta(const runge_kutta_method method);
247 
251  void
252  initialize(const runge_kutta_method method) override;
253 
265  double
266  evolve_one_time_step(
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
283  evolve_one_time_step(
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
312  compute_stages(
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 
322  Status status;
323  };
324 
325 
326 
331  template <typename VectorType>
332  class ImplicitRungeKutta : public RungeKutta<VectorType>
333  {
334  public:
335  using RungeKutta<VectorType>::evolve_one_time_step;
336 
342  ImplicitRungeKutta() = default;
343 
349  ImplicitRungeKutta(const runge_kutta_method method,
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
371  evolve_one_time_step(
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
416  compute_stages(
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
440  compute_residual(
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 
454  bool skip_linear_combi;
455 
459  unsigned int max_it;
460 
464  double tolerance;
465 
469  Status status;
470  };
471 
472 
473 
478  template <typename VectorType>
479  class EmbeddedExplicitRungeKutta : public RungeKutta<VectorType>
480  {
481  public:
482  using RungeKutta<VectorType>::evolve_one_time_step;
483 
489  EmbeddedExplicitRungeKutta() = default;
490 
495  EmbeddedExplicitRungeKutta(const runge_kutta_method method,
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 
506  ~EmbeddedExplicitRungeKutta() override
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
535  evolve_one_time_step(
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
552  evolve_one_time_step(
553  const std::function<VectorType(const double, const VectorType &)> &f,
554  double t,
555  double delta_t,
556  VectorType &y);
557 
561  void
562  set_time_adaptation_parameters(const double coarsen_param,
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
595  compute_stages(
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 
606  double coarsen_param;
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 
640  bool last_same_as_first;
641 
645  std::vector<double> b1;
646 
650  std::vector<double> b2;
651 
656  VectorType *last_stage;
657 
661  Status status;
662  };
663 } // namespace TimeStepping
664 
665 DEAL_II_NAMESPACE_CLOSE
666 
667 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
TimeStepping::RungeKutta::b
std::vector< double > b
Definition: time_stepping.h:211
TimeStepping::RungeKutta::c
std::vector< double > c
Definition: time_stepping.h:216
TimeStepping::RungeKutta::a
std::vector< std::vector< double > > a
Definition: time_stepping.h:221
TimeStepping::embedded_runge_kutta_time_step
embedded_runge_kutta_time_step
Definition: time_stepping.h:82
Generated by   doxygen 1.8.11