MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ode.cpp
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.googlecode.com.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "operator.hpp"
13 #include "ode.hpp"
14 
15 namespace mfem
16 {
17 
19 {
20  ODESolver::Init(_f);
21  dxdt.SetSize(f->Width());
22 }
23 
24 void ForwardEulerSolver::Step(Vector &x, double &t, double &dt)
25 {
26  f->SetTime(t);
27  f->Mult(x, dxdt);
28  x.Add(dt, dxdt);
29  t += dt;
30 }
31 
32 
34 {
35  ODESolver::Init(_f);
36  int n = f->Width();
37  dxdt.SetSize(n);
38  x1.SetSize(n);
39 }
40 
41 void RK2Solver::Step(Vector &x, double &t, double &dt)
42 {
43  // 0 |
44  // a | a
45  // ---+--------
46  // | 1-b b b = 1/(2a)
47 
48  const double b = 0.5/a;
49 
50  f->SetTime(t);
51  f->Mult(x, dxdt);
52  add(x, (1. - b)*dt, dxdt, x1);
53  x.Add(a*dt, dxdt);
54 
55  f->SetTime(t + a*dt);
56  f->Mult(x, dxdt);
57  add(x1, b*dt, dxdt, x);
58  t += dt;
59 }
60 
61 
63 {
64  ODESolver::Init(_f);
65  int n = f->Width();
66  y.SetSize(n);
67  k.SetSize(n);
68 }
69 
70 void RK3SSPSolver::Step(Vector &x, double &t, double &dt)
71 {
72  // x0 = x, t0 = t, k0 = dt*f(t0, x0)
73  f->SetTime(t);
74  f->Mult(x, k);
75 
76  // x1 = x + k0, t1 = t + dt, k1 = dt*f(t1, x1)
77  add(x, dt, k, y);
78  f->SetTime(t + dt);
79  f->Mult(y, k);
80 
81  // x2 = 3/4*x + 1/4*(x1 + k1), t2 = t + 1/2*dt, k2 = dt*f(t2, x2)
82  y.Add(dt, k);
83  add(3./4, x, 1./4, y, y);
84  f->SetTime(t + dt/2);
85  f->Mult(y, k);
86 
87  // x3 = 1/3*x + 2/3*(x2 + k2), t3 = t + dt
88  y.Add(dt, k);
89  add(1./3, x, 2./3, y, x);
90  t += dt;
91 }
92 
93 
95 {
96  ODESolver::Init(_f);
97  int n = f->Width();
98  y.SetSize(n);
99  k.SetSize(n);
100  z.SetSize(n);
101 }
102 
103 void RK4Solver::Step(Vector &x, double &t, double &dt)
104 {
105  // 0 |
106  // 1/2 | 1/2
107  // 1/2 | 0 1/2
108  // 1 | 0 0 1
109  // -----+-------------------
110  // | 1/6 1/3 1/3 1/6
111 
112  f->SetTime(t);
113  f->Mult(x, k); // k1
114  add(x, dt/2, k, y);
115  add(x, dt/6, k, z);
116 
117  f->SetTime(t + dt/2);
118  f->Mult(y, k); // k2
119  add(x, dt/2, k, y);
120  z.Add(dt/3, k);
121 
122  f->Mult(y, k); // k3
123  add(x, dt, k, y);
124  z.Add(dt/3, k);
125 
126  f->SetTime(t + dt);
127  f->Mult(y, k); // k4
128  add(z, dt/6, k, x);
129  t += dt;
130 }
131 
132 ExplicitRKSolver::ExplicitRKSolver(int _s, const double *_a, const double *_b,
133  const double *_c)
134 {
135  s = _s;
136  a = _a;
137  b = _b;
138  c = _c;
139  k = new Vector[s];
140 }
141 
143 {
144  ODESolver::Init(_f);
145  int n = f->Width();
146  y.SetSize(n);
147  for (int i = 0; i < s; i++)
148  k[i].SetSize(n);
149 }
150 
151 void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
152 {
153  // 0 |
154  // c[0] | a[0]
155  // c[1] | a[1] a[2]
156  // ... | ...
157  // c[s-2] | ... a[s(s-1)/2-1]
158  // --------+---------------------
159  // | b[0] b[1] ... b[s-1]
160 
161  f->SetTime(t);
162  f->Mult(x, k[0]);
163  for (int l = 0, i = 1; i < s; i++)
164  {
165  add(x, a[l++]*dt, k[0], y);
166  for (int j = 1; j < i; j++)
167  y.Add(a[l++]*dt, k[j]);
168 
169  f->SetTime(t + c[i-1]*dt);
170  f->Mult(y, k[i]);
171  }
172  for (int i = 0; i < s; i++)
173  x.Add(b[i]*dt, k[i]);
174  t += dt;
175 }
176 
178 {
179  delete [] k;
180 }
181 
182 const double RK6Solver::a[] = {
183  .6e-1,
184  .1923996296296296296296296296296296296296e-1,
185  .7669337037037037037037037037037037037037e-1,
186  .35975e-1,
187  0.,
188  .107925,
189  1.318683415233148260919747276431735612861,
190  0.,
191  -5.042058063628562225427761634715637693344,
192  4.220674648395413964508014358283902080483,
193  -41.87259166432751461803757780644346812905,
194  0.,
195  159.4325621631374917700365669070346830453,
196  -122.1192135650100309202516203389242140663,
197  5.531743066200053768252631238332999150076,
198  -54.43015693531650433250642051294142461271,
199  0.,
200  207.0672513650184644273657173866509835987,
201  -158.6108137845899991828742424365058599469,
202  6.991816585950242321992597280791793907096,
203  -.1859723106220323397765171799549294623692e-1,
204  -54.66374178728197680241215648050386959351,
205  0.,
206  207.9528062553893734515824816699834244238,
207  -159.2889574744995071508959805871426654216,
208  7.018743740796944434698170760964252490817,
209  -.1833878590504572306472782005141738268361e-1,
210  -.5119484997882099077875432497245168395840e-3
211 };
212 const double RK6Solver::b[] = {
213  .3438957868357036009278820124728322386520e-1,
214  0.,
215  0.,
216  .2582624555633503404659558098586120858767,
217  .4209371189673537150642551514069801967032,
218  4.405396469669310170148836816197095664891,
219  -176.4831190242986576151740942499002125029,
220  172.3641334014150730294022582711902413315
221 };
222 const double RK6Solver::c[] = {
223  .6e-1,
224  .9593333333333333333333333333333333333333e-1,
225  .1439,
226  .4973,
227  .9725,
228  .9995,
229  1.,
230 };
231 
232 const double RK8Solver::a[] = {
233  .5e-1,
234  -.69931640625e-2,
235  .1135556640625,
236  .399609375e-1,
237  0.,
238  .1198828125,
239  .3613975628004575124052940721184028345129,
240  0.,
241  -1.341524066700492771819987788202715834917,
242  1.370126503900035259414693716084313000404,
243  .490472027972027972027972027972027972028e-1,
244  0.,
245  0.,
246  .2350972042214404739862988335493427143122,
247  .180855592981356728810903963653454488485,
248  .6169289044289044289044289044289044289044e-1,
249  0.,
250  0.,
251  .1123656831464027662262557035130015442303,
252  -.3885046071451366767049048108111244567456e-1,
253  .1979188712522045855379188712522045855379e-1,
254  -1.767630240222326875735597119572145586714,
255  0.,
256  0.,
257  -62.5,
258  -6.061889377376669100821361459659331999758,
259  5.650823198222763138561298030600840174201,
260  65.62169641937623283799566054863063741227,
261  -1.180945066554970799825116282628297957882,
262  0.,
263  0.,
264  -41.50473441114320841606641502701994225874,
265  -4.434438319103725011225169229846100211776,
266  4.260408188586133024812193710744693240761,
267  43.75364022446171584987676829438379303004,
268  .787142548991231068744647504422630755086e-2,
269  -1.281405999441488405459510291182054246266,
270  0.,
271  0.,
272  -45.04713996013986630220754257136007322267,
273  -4.731362069449576477311464265491282810943,
274  4.514967016593807841185851584597240996214,
275  47.44909557172985134869022392235929015114,
276  .1059228297111661135687393955516542875228e-1,
277  -.5746842263844616254432318478286296232021e-2,
278  -1.724470134262485191756709817484481861731,
279  0.,
280  0.,
281  -60.92349008483054016518434619253765246063,
282  -5.95151837622239245520283276706185486829,
283  5.556523730698456235979791650843592496839,
284  63.98301198033305336837536378635995939281,
285  .1464202825041496159275921391759452676003e-1,
286  .6460408772358203603621865144977650714892e-1,
287  -.7930323169008878984024452548693373291447e-1,
288  -3.301622667747079016353994789790983625569,
289  0.,
290  0.,
291  -118.011272359752508566692330395789886851,
292  -10.14142238845611248642783916034510897595,
293  9.139311332232057923544012273556827000619,
294  123.3759428284042683684847180986501894364,
295  4.623244378874580474839807625067630924792,
296  -3.383277738068201923652550971536811240814,
297  4.527592100324618189451265339351129035325,
298  -5.828495485811622963193088019162985703755
299 };
300 const double RK8Solver::b[] = {
301  .4427989419007951074716746668098518862111e-1,
302  0.,
303  0.,
304  0.,
305  0.,
306  .3541049391724448744815552028733568354121,
307  .2479692154956437828667629415370663023884,
308  -15.69420203883808405099207034271191213468,
309  25.08406496555856261343930031237186278518,
310  -31.73836778626027646833156112007297739997,
311  22.93828327398878395231483560344797018313,
312  -.2361324633071542145259900641263517600737
313 };
314 const double RK8Solver::c[] = {
315  .5e-1,
316  .1065625,
317  .15984375,
318  .39,
319  .465,
320  .155,
321  .943,
322  .901802041735856958259707940678372149956,
323  .909,
324  .94,
325  1.,
326 };
327 
328 
330 {
331  ODESolver::Init(_f);
332  k.SetSize(f->Width());
333 }
334 
335 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
336 {
337  f->SetTime(t + dt);
338  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
339  x.Add(dt, k);
340  t += dt;
341 }
342 
343 
345 {
346  ODESolver::Init(_f);
347  k.SetSize(f->Width());
348 }
349 
350 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
351 {
352  f->SetTime(t + dt/2);
353  f->ImplicitSolve(dt/2, x, k);
354  x.Add(dt, k);
355  t += dt;
356 }
357 
358 
360 {
361  if (gamma_opt == 0)
362  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
363  else if (gamma_opt == 2)
364  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
365  else if (gamma_opt == 3)
366  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
367  else
368  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
369 }
370 
372 {
373  ODESolver::Init(_f);
374  k.SetSize(f->Width());
375  y.SetSize(f->Width());
376 }
377 
378 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
379 {
380  // with a = gamma:
381  // a | a
382  // 1-a | 1-2a a
383  // ------+-----------
384  // | 1/2 1/2
385  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
386  f->SetTime(t + gamma*dt);
387  f->ImplicitSolve(gamma*dt, x, k);
388  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
389  x.Add(dt/2, k);
390 
391  f->SetTime(t + (1.-gamma)*dt);
392  f->ImplicitSolve(gamma*dt, y, k);
393  x.Add(dt/2, k);
394  t += dt;
395 }
396 
397 
399 {
400  ODESolver::Init(_f);
401  k.SetSize(f->Width());
402  y.SetSize(f->Width());
403  z.SetSize(f->Width());
404 }
405 
406 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
407 {
408  // a | a
409  // 1/2 | 1/2-a a
410  // 1-a | 2a 1-4a a
411  // ------+--------------------
412  // | b 1-2b b
413  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
414  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
415  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
416 
417  f->SetTime(t + a*dt);
418  f->ImplicitSolve(a*dt, x, k);
419  add(x, (0.5-a)*dt, k, y);
420  add(x, (2.*a)*dt, k, z);
421  x.Add(b*dt, k);
422 
423  f->SetTime(t + dt/2);
424  f->ImplicitSolve(a*dt, y, k);
425  z.Add((1.-4.*a)*dt, k);
426  x.Add((1.-2.*b)*dt, k);
427 
428  f->SetTime(t + (1.-a)*dt);
429  f->ImplicitSolve(a*dt, z, k);
430  x.Add(b*dt, k);
431  t += dt;
432 }
433 
434 
436 {
437  ODESolver::Init(_f);
438  k.SetSize(f->Width());
439  y.SetSize(f->Width());
440 }
441 
442 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
443 {
444  // a | a
445  // c | c-a a
446  // 1 | b 1-a-b a
447  // -----+----------------
448  // | b 1-a-b a
449  const double a = 0.435866521508458999416019;
450  const double b = 1.20849664917601007033648;
451  const double c = 0.717933260754229499708010;
452 
453  f->SetTime(t + a*dt);
454  f->ImplicitSolve(a*dt, x, k);
455  add(x, (c-a)*dt, k, y);
456  x.Add(b*dt, k);
457 
458  f->SetTime(t + c*dt);
459  f->ImplicitSolve(a*dt, y, k);
460  x.Add((1.-a-b)*dt, k);
461 
462  f->SetTime(t + dt);
463  f->ImplicitSolve(a*dt, x, k);
464  x.Add(a*dt, k);
465  t += dt;
466 }
467 
468 }
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:33
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:371
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:103
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:248
Base abstract class for time dependent operators: (x,t) -&gt; f(x,t)
Definition: operator.hpp:68
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
virtual void SetTime(const double _t)
Definition: operator.hpp:86
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:151
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application.
virtual void Init(TimeDependentOperator &_f)
Definition: ode.hpp:30
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:442
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:227
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:406
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:344
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:94
virtual ~ExplicitRKSolver()
Definition: ode.cpp:177
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:41
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:62
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Definition: operator.hpp:92
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:350
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)
Definition: ode.cpp:335
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:329
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:398
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:18
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:168
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:70
Vector data type.
Definition: vector.hpp:29
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:435
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:378
TimeDependentOperator * f
Definition: ode.hpp:25
virtual void Step(Vector &x, double &t, double &dt)
Definition: ode.cpp:24
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:359
virtual void Init(TimeDependentOperator &_f)
Definition: ode.cpp:142