MFEM  v3.3 Finite element discretization library
ode.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
11
12 #ifndef MFEM_ODE
13 #define MFEM_ODE
14
15 #include "../config/config.hpp"
16 #include "operator.hpp"
17
18 namespace mfem
19 {
20
21 /// Abstract class for solving systems of ODEs: dx/dt = f(x,t)
22 class ODESolver
23 {
24 protected:
25  /// Pointer to the associated TimeDependentOperator.
26  TimeDependentOperator *f; // f(.,t) : R^n --> R^n
27
28 public:
29  ODESolver() : f(NULL) { }
30
31  /// Associate a TimeDependentOperator with the ODE solver.
32  /** This method has to be called:
33  - Before the first call to Step().
34  - When the dimensions of the associated TimeDependentOperator change.
35  - When a time stepping sequence has to be restarted.
36  - To change the associated TimeDependentOperator. */
37  virtual void Init(TimeDependentOperator &f)
38  {
39  this->f = &f;
40  }
41
42  /** @brief Perform a time step from time @a t [in] to time @a t [out] based
43  on the requested step size @a dt [in]. */
44  /** @param[in,out] x Approximate solution.
45  @param[in,out] t Time associated with the approximate solution @a x.
46  @param[in,out] dt Time step size.
47
48  The following rules describe the common behavior of the method:
49  - The input @a x [in] is the approximate solution for the input time
50  @a t [in].
51  - The input @a dt [in] is the desired time step size, defining the desired
52  target time: t [target] = @a t [in] + @a dt [in].
53  - The output @a x [out] is the approximate solution for the output time
54  @a t [out].
55  - The output @a dt [out] is the last time step taken by the method which
56  may be smaller or larger than the input @a dt [in] value, e.g. because
57  of time step control.
58  - The method may perform more than one time step internally; in this case
59  @a dt [out] is the last internal time step size.
60  - The output value of @a t [out] may be smaller or larger than
61  t [target], however, it is not smaller than @a t [in] + @a dt [out], if
62  at least one internal time step was performed.
63  - The value @a x [out] may be obtained by interpolation using internally
64  stored data.
65  - In some cases, the contents of @a x [in] may not be used, e.g. when
66  @a x [out] from a previous Step() call was obtained by interpolation.
67  - In consecutive calls to this method, the output @a t [out] of one
68  Step() call has to be the same as the input @a t [in] to the next
69  Step() call.
70  - If the previous rule has to be broken, e.g. to restart a time stepping
71  sequence, then the ODE solver must be re-initialized by calling Init()
72  between the two Step() calls. */
73  virtual void Step(Vector &x, double &t, double &dt) = 0;
74
75  /// Perform time integration from time @a t [in] to time @a tf [in].
76  /** @param[in,out] x Approximate solution.
77  @param[in,out] t Time associated with the approximate solution @a x.
78  @param[in,out] dt Time step size.
79  @param[in] tf Requested final time.
80
81  The default implementation makes consecutive calls to Step() until
82  reaching @a tf.
83  The following rules describe the common behavior of the method:
84  - The input @a x [in] is the approximate solution for the input time
85  @a t [in].
86  - The input @a dt [in] is the initial time step size.
87  - The output @a dt [out] is the last time step taken by the method which
88  may be smaller or larger than the input @a dt [in] value, e.g. because
89  of time step control.
90  - The output value of @a t [out] is not smaller than @a tf [in]. */
91  virtual void Run(Vector &x, double &t, double &dt, double tf)
92  {
93  while (t < tf) { Step(x, t, dt); }
94  }
95
96  virtual ~ODESolver() { }
97 };
98
99
100 /// The classical forward Euler method
102 {
103 private:
104  Vector dxdt;
105
106 public:
107  virtual void Init(TimeDependentOperator &_f);
108
109  virtual void Step(Vector &x, double &t, double &dt);
110 };
111
112
113 /** A family of explicit second-order RK2 methods. Some choices for the
114  parameter 'a' are:
115  a = 1/2 - the midpoint method
116  a = 1 - Heun's method
117  a = 2/3 - default, has minimal truncation error. */
118 class RK2Solver : public ODESolver
119 {
120 private:
121  double a;
122  Vector dxdt, x1;
123
124 public:
125  RK2Solver(const double _a = 2./3.) : a(_a) { }
126
127  virtual void Init(TimeDependentOperator &_f);
128
129  virtual void Step(Vector &x, double &t, double &dt);
130 };
131
132
133 /// Third-order, strong stability preserving (SSP) Runge-Kutta method
134 class RK3SSPSolver : public ODESolver
135 {
136 private:
137  Vector y, k;
138
139 public:
140  virtual void Init(TimeDependentOperator &_f);
141
142  virtual void Step(Vector &x, double &t, double &dt);
143 };
144
145
146 /// The classical explicit forth-order Runge-Kutta method, RK4
147 class RK4Solver : public ODESolver
148 {
149 private:
150  Vector y, k, z;
151
152 public:
153  virtual void Init(TimeDependentOperator &_f);
154
155  virtual void Step(Vector &x, double &t, double &dt);
156 };
157
158
159 /** An explicit Runge-Kutta method corresponding to a general Butcher tableau
160  +--------+----------------------+
161  | c[0] | a[0] |
162  | c[1] | a[1] a[2] |
163  | ... | ... |
164  | c[s-2] | ... a[s(s-1)/2-1] |
165  +--------+----------------------+
166  | | b[0] b[1] ... b[s-1] |
167  +--------+----------------------+ */
169 {
170 private:
171  int s;
172  const double *a, *b, *c;
173  Vector y, *k;
174
175 public:
176  ExplicitRKSolver(int _s, const double *_a, const double *_b,
177  const double *_c);
178
179  virtual void Init(TimeDependentOperator &_f);
180
181  virtual void Step(Vector &x, double &t, double &dt);
182
183  virtual ~ExplicitRKSolver();
184 };
185
186
187 /** An 8-stage, 6th order RK method. From Verner's "efficient" 9-stage 6(5)
188  pair. */
190 {
191 private:
192  static const double a[28], b[8], c[7];
193
194 public:
195  RK6Solver() : ExplicitRKSolver(8, a, b, c) { }
196 };
197
198
199 /** A 12-stage, 8th order RK method. From Verner's "efficient" 13-stage 8(7)
200  pair. */
202 {
203 private:
204  static const double a[66], b[12], c[11];
205
206 public:
207  RK8Solver() : ExplicitRKSolver(12, a, b, c) { }
208 };
209
210
211 /// Backward Euler ODE solver. L-stable.
213 {
214 protected:
216
217 public:
218  virtual void Init(TimeDependentOperator &_f);
219
220  virtual void Step(Vector &x, double &t, double &dt);
221 };
222
223
224 /// Implicit midpoint method. A-stable, not L-stable.
226 {
227 protected:
229
230 public:
231  virtual void Init(TimeDependentOperator &_f);
232
233  virtual void Step(Vector &x, double &t, double &dt);
234 };
235
236
237 /** Two stage, singly diagonal implicit Runge-Kutta (SDIRK) methods;
238  the choices for gamma_opt are:
239  0 - 3rd order method, not A-stable
240  1 - 3rd order method, A-stable, not L-stable (default)
241  2 - 2nd order method, L-stable
242  3 - 2nd order method, L-stable (has solves outside [t,t+dt]). */
243 class SDIRK23Solver : public ODESolver
244 {
245 protected:
246  double gamma;
248
249 public:
250  SDIRK23Solver(int gamma_opt = 1);
251
252  virtual void Init(TimeDependentOperator &_f);
253
254  virtual void Step(Vector &x, double &t, double &dt);
255 };
256
257
258 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
259  order 4. A-stable, not L-stable. */
260 class SDIRK34Solver : public ODESolver
261 {
262 protected:
263  Vector k, y, z;
264
265 public:
266  virtual void Init(TimeDependentOperator &_f);
267
268  virtual void Step(Vector &x, double &t, double &dt);
269 };
270
271
272 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
273  order 3. L-stable. */
274 class SDIRK33Solver : public ODESolver
275 {
276 protected:
278
279 public:
280  virtual void Init(TimeDependentOperator &_f);
281
282  virtual void Step(Vector &x, double &t, double &dt);
283 };
284
285 }
286
287 #endif
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:33
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:391
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:103
Base abstract class for time dependent operators.
Definition: operator.hpp:140
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:153
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:212
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:462
virtual ~ODESolver()
Definition: ode.hpp:96
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:426
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:356
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:94
virtual ~ExplicitRKSolver()
Definition: ode.cpp:183
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:41
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:62
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:362
ExplicitRKSolver(int _s, const double *_a, const double *_b, const double *_c)
Definition: ode.cpp:132
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:347
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:341
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:418
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:134
RK2Solver(const double _a=2./3.)
Definition: ode.hpp:125
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:225
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:70
Vector data type.
Definition: vector.hpp:36
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:455
virtual void Run(Vector &x, double &t, double &dt, double tf)
Perform time integration from time t [in] to time tf [in].
Definition: ode.hpp:91
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:398
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
The classical forward Euler method.
Definition: ode.hpp:101
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:24
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:371
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:142