MFEM  v4.0 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
28
29 public:
31
32  /// Associate a TimeDependentOperator with the ODE solver.
33  /** This method has to be called:
34  - Before the first call to Step().
35  - When the dimensions of the associated TimeDependentOperator change.
36  - When a time stepping sequence has to be restarted.
37  - To change the associated TimeDependentOperator. */
38  virtual void Init(TimeDependentOperator &f);
39
40  /** @brief Perform a time step from time @a t [in] to time @a t [out] based
41  on the requested step size @a dt [in]. */
42  /** @param[in,out] x Approximate solution.
43  @param[in,out] t Time associated with the approximate solution @a x.
44  @param[in,out] dt Time step size.
45
46  The following rules describe the common behavior of the method:
47  - The input @a x [in] is the approximate solution for the input time
48  @a t [in].
49  - The input @a dt [in] is the desired time step size, defining the desired
50  target time: t [target] = @a t [in] + @a dt [in].
51  - The output @a x [out] is the approximate solution for the output time
52  @a t [out].
53  - The output @a dt [out] is the last time step taken by the method which
54  may be smaller or larger than the input @a dt [in] value, e.g. because
55  of time step control.
56  - The method may perform more than one time step internally; in this case
57  @a dt [out] is the last internal time step size.
58  - The output value of @a t [out] may be smaller or larger than
59  t [target], however, it is not smaller than @a t [in] + @a dt [out], if
60  at least one internal time step was performed.
61  - The value @a x [out] may be obtained by interpolation using internally
62  stored data.
63  - In some cases, the contents of @a x [in] may not be used, e.g. when
64  @a x [out] from a previous Step() call was obtained by interpolation.
65  - In consecutive calls to this method, the output @a t [out] of one
66  Step() call has to be the same as the input @a t [in] to the next
67  Step() call.
68  - If the previous rule has to be broken, e.g. to restart a time stepping
69  sequence, then the ODE solver must be re-initialized by calling Init()
70  between the two Step() calls. */
71  virtual void Step(Vector &x, double &t, double &dt) = 0;
72
73  /// Perform time integration from time @a t [in] to time @a tf [in].
74  /** @param[in,out] x Approximate solution.
75  @param[in,out] t Time associated with the approximate solution @a x.
76  @param[in,out] dt Time step size.
77  @param[in] tf Requested final time.
78
79  The default implementation makes consecutive calls to Step() until
80  reaching @a tf.
81  The following rules describe the common behavior of the method:
82  - The input @a x [in] is the approximate solution for the input time
83  @a t [in].
84  - The input @a dt [in] is the initial time step size.
85  - The output @a dt [out] is the last time step taken by the method which
86  may be smaller or larger than the input @a dt [in] value, e.g. because
87  of time step control.
88  - The output value of @a t [out] is not smaller than @a tf [in]. */
89  virtual void Run(Vector &x, double &t, double &dt, double tf)
90  {
91  while (t < tf) { Step(x, t, dt); }
92  }
93
94  virtual ~ODESolver() { }
95 };
96
97
98 /// The classical forward Euler method
100 {
101 private:
102  Vector dxdt;
103
104 public:
105  virtual void Init(TimeDependentOperator &_f);
106
107  virtual void Step(Vector &x, double &t, double &dt);
108 };
109
110
111 /** A family of explicit second-order RK2 methods. Some choices for the
112  parameter 'a' are:
113  a = 1/2 - the midpoint method
114  a = 1 - Heun's method
115  a = 2/3 - default, has minimal truncation error. */
116 class RK2Solver : public ODESolver
117 {
118 private:
119  double a;
120  Vector dxdt, x1;
121
122 public:
123  RK2Solver(const double _a = 2./3.) : a(_a) { }
124
125  virtual void Init(TimeDependentOperator &_f);
126
127  virtual void Step(Vector &x, double &t, double &dt);
128 };
129
130
131 /// Third-order, strong stability preserving (SSP) Runge-Kutta method
132 class RK3SSPSolver : public ODESolver
133 {
134 private:
135  Vector y, k;
136
137 public:
138  virtual void Init(TimeDependentOperator &_f);
139
140  virtual void Step(Vector &x, double &t, double &dt);
141 };
142
143
144 /// The classical explicit forth-order Runge-Kutta method, RK4
145 class RK4Solver : public ODESolver
146 {
147 private:
148  Vector y, k, z;
149
150 public:
151  virtual void Init(TimeDependentOperator &_f);
152
153  virtual void Step(Vector &x, double &t, double &dt);
154 };
155
156
157 /** An explicit Runge-Kutta method corresponding to a general Butcher tableau
158  +--------+----------------------+
159  | c | a |
160  | c | a a |
161  | ... | ... |
162  | c[s-2] | ... a[s(s-1)/2-1] |
163  +--------+----------------------+
164  | | b b ... b[s-1] |
165  +--------+----------------------+ */
167 {
168 private:
169  int s;
170  const double *a, *b, *c;
171  Vector y, *k;
172
173 public:
174  ExplicitRKSolver(int _s, const double *_a, const double *_b,
175  const double *_c);
176
177  virtual void Init(TimeDependentOperator &_f);
178
179  virtual void Step(Vector &x, double &t, double &dt);
180
181  virtual ~ExplicitRKSolver();
182 };
183
184
185 /** An 8-stage, 6th order RK method. From Verner's "efficient" 9-stage 6(5)
186  pair. */
188 {
189 private:
190  static const double a, b, c;
191
192 public:
193  RK6Solver() : ExplicitRKSolver(8, a, b, c) { }
194 };
195
196
197 /** A 12-stage, 8th order RK method. From Verner's "efficient" 13-stage 8(7)
198  pair. */
200 {
201 private:
202  static const double a, b, c;
203
204 public:
205  RK8Solver() : ExplicitRKSolver(12, a, b, c) { }
206 };
207
208
209 /// Backward Euler ODE solver. L-stable.
211 {
212 protected:
214
215 public:
216  virtual void Init(TimeDependentOperator &_f);
217
218  virtual void Step(Vector &x, double &t, double &dt);
219 };
220
221
222 /// Implicit midpoint method. A-stable, not L-stable.
224 {
225 protected:
227
228 public:
229  virtual void Init(TimeDependentOperator &_f);
230
231  virtual void Step(Vector &x, double &t, double &dt);
232 };
233
234
235 /** Two stage, singly diagonal implicit Runge-Kutta (SDIRK) methods;
236  the choices for gamma_opt are:
237  0 - 3rd order method, not A-stable
238  1 - 3rd order method, A-stable, not L-stable (default)
239  2 - 2nd order method, L-stable
240  3 - 2nd order method, L-stable (has solves outside [t,t+dt]). */
241 class SDIRK23Solver : public ODESolver
242 {
243 protected:
244  double gamma;
246
247 public:
248  SDIRK23Solver(int gamma_opt = 1);
249
250  virtual void Init(TimeDependentOperator &_f);
251
252  virtual void Step(Vector &x, double &t, double &dt);
253 };
254
255
256 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
257  order 4. A-stable, not L-stable. */
258 class SDIRK34Solver : public ODESolver
259 {
260 protected:
261  Vector k, y, z;
262
263 public:
264  virtual void Init(TimeDependentOperator &_f);
265
266  virtual void Step(Vector &x, double &t, double &dt);
267 };
268
269
270 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
271  order 3. L-stable. */
272 class SDIRK33Solver : public ODESolver
273 {
274 protected:
276
277 public:
278  virtual void Init(TimeDependentOperator &_f);
279
280  virtual void Step(Vector &x, double &t, double &dt);
281 };
282
283
284 /// Generalized-alpha ODE solver from "A generalized-α method for integrating
285 /// the filtered Navier–Stokes equations with a stabilized finite element
286 /// method" by K.E. Jansen, C.H. Whiting and G.M. Hulbert.
288 {
289 protected:
292  bool first;
293
294  void SetRhoInf(double rho_inf);
295  void PrintProperties(std::ostream &out = mfem::out);
296 public:
297
298  GeneralizedAlphaSolver(double rho = 1.0) { SetRhoInf(rho); };
299
300  virtual void Init(TimeDependentOperator &_f);
301
302  virtual void Step(Vector &x, double &t, double &dt);
303 };
304
305
306 /// The SIASolver class is based on the Symplectic Integration Algorithm
307 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
308 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,
309 /// Vol. 92, pages 230-256 (1991).
310
311 /** The Symplectic Integration Algorithm (SIA) is designed for systems of first
312  order ODEs derived from a Hamiltonian.
313  H(q,p,t) = T(p) + V(q,t)
314  Which leads to the equations:
315  dq/dt = dT/dp
316  dp/dt = -dV/dq
317  In the integrator the operators P and F are defined to be:
318  P = dT/dp
319  F = -dV/dq
320  */
322 {
323 public:
324  SIASolver() : F_(NULL), P_(NULL) {}
325
326  virtual void Init(Operator &P, TimeDependentOperator & F);
327
328  virtual void Step(Vector &q, Vector &p, double &t, double &dt) = 0;
329
330  virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
331  {
332  while (t < tf) { Step(q, p, t, dt); }
333  }
334
335  virtual ~SIASolver() {}
336
337 protected:
338  TimeDependentOperator * F_; // p_{i+1} = p_{i} + dt F(q_{i})
339  Operator * P_; // q_{i+1} = q_{i} + dt P(p_{i+1})
340
341  mutable Vector dp_;
342  mutable Vector dq_;
343 };
344
345 // First Order Symplectic Integration Algorithm
346 class SIA1Solver : public SIASolver
347 {
348 public:
350  void Step(Vector &q, Vector &p, double &t, double &dt);
351 };
352
353 // Second Order Symplectic Integration Algorithm
354 class SIA2Solver : public SIASolver
355 {
356 public:
358  void Step(Vector &q, Vector &p, double &t, double &dt);
359 };
360
361 // Variable order Symplectic Integration Algorithm (orders 1-4)
362 class SIAVSolver : public SIASolver
363 {
364 public:
365  SIAVSolver(int order);
366  void Step(Vector &q, Vector &p, double &t, double &dt);
367
368 private:
369  int order_;
370
371  Array<double> a_;
372  Array<double> b_;
373 };
374
375 }
376
377 #endif
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:515
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:571
void SetRhoInf(double rho_inf)
Definition: ode.cpp:505
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:580
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:397
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:109
Base abstract class for time dependent operators.
Definition: operator.hpp:162
MemoryType mem_type
Definition: ode.hpp:27
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]...
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
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:159
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:210
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:468
virtual ~ODESolver()
Definition: ode.hpp:94
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:432
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:362
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
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:47
Vector dq_
Definition: ode.hpp:342
virtual ~SIASolver()
Definition: ode.hpp:335
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
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:368
ExplicitRKSolver(int _s, const double *_a, const double *_b, const double *_c)
Definition: ode.cpp:138
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:353
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:347
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:495
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:424
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:145
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
TimeDependentOperator * F_
Definition: ode.hpp:338
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:27
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:132
Vector dp_
Definition: ode.hpp:341
Operator * P_
Definition: ode.hpp:339
RK2Solver(const double _a=2./3.)
Definition: ode.hpp:123
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:223
Host memory; using new[] and delete[].
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:593
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:650
SIAVSolver(int order)
Definition: ode.cpp:608
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:542
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:76
virtual void Step(Vector &q, Vector &p, double &t, double &dt)=0
GeneralizedAlphaSolver(double rho=1.0)
Definition: ode.hpp:298
virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
Definition: ode.hpp:330
Vector data type.
Definition: vector.hpp:48
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:461
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:89
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:404
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
The classical forward Euler method.
Definition: ode.hpp:99
Abstract operator.
Definition: operator.hpp:21
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:30
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:377
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148