MFEM  v3.4 Finite element discretization library
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.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 #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);
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);
54
55  f->SetTime(t + a*dt);
56  f->Mult(x, dxdt);
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)
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)
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
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
116
117  f->SetTime(t + dt/2);
118  f->Mult(y, k); // k2
121
122  f->Mult(y, k); // k3
125
126  f->SetTime(t + dt);
127  f->Mult(y, k); // k4
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  {
149  k[i].SetSize(n);
150  }
151 }
152
153 void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
154 {
155  // 0 |
156  // c[0] | a[0]
157  // c[1] | a[1] a[2]
158  // ... | ...
159  // c[s-2] | ... a[s(s-1)/2-1]
160  // --------+---------------------
161  // | b[0] b[1] ... b[s-1]
162
163  f->SetTime(t);
164  f->Mult(x, k[0]);
165  for (int l = 0, i = 1; i < s; i++)
166  {
168  for (int j = 1; j < i; j++)
169  {
171  }
172
173  f->SetTime(t + c[i-1]*dt);
174  f->Mult(y, k[i]);
175  }
176  for (int i = 0; i < s; i++)
177  {
179  }
180  t += dt;
181 }
182
184 {
185  delete [] k;
186 }
187
188 const double RK6Solver::a[] =
189 {
190  .6e-1,
191  .1923996296296296296296296296296296296296e-1,
192  .7669337037037037037037037037037037037037e-1,
193  .35975e-1,
194  0.,
195  .107925,
196  1.318683415233148260919747276431735612861,
197  0.,
198  -5.042058063628562225427761634715637693344,
199  4.220674648395413964508014358283902080483,
200  -41.87259166432751461803757780644346812905,
201  0.,
202  159.4325621631374917700365669070346830453,
203  -122.1192135650100309202516203389242140663,
204  5.531743066200053768252631238332999150076,
205  -54.43015693531650433250642051294142461271,
206  0.,
207  207.0672513650184644273657173866509835987,
208  -158.6108137845899991828742424365058599469,
209  6.991816585950242321992597280791793907096,
210  -.1859723106220323397765171799549294623692e-1,
211  -54.66374178728197680241215648050386959351,
212  0.,
213  207.9528062553893734515824816699834244238,
214  -159.2889574744995071508959805871426654216,
215  7.018743740796944434698170760964252490817,
216  -.1833878590504572306472782005141738268361e-1,
217  -.5119484997882099077875432497245168395840e-3
218 };
219 const double RK6Solver::b[] =
220 {
221  .3438957868357036009278820124728322386520e-1,
222  0.,
223  0.,
224  .2582624555633503404659558098586120858767,
225  .4209371189673537150642551514069801967032,
226  4.405396469669310170148836816197095664891,
227  -176.4831190242986576151740942499002125029,
228  172.3641334014150730294022582711902413315
229 };
230 const double RK6Solver::c[] =
231 {
232  .6e-1,
233  .9593333333333333333333333333333333333333e-1,
234  .1439,
235  .4973,
236  .9725,
237  .9995,
238  1.,
239 };
240
241 const double RK8Solver::a[] =
242 {
243  .5e-1,
244  -.69931640625e-2,
245  .1135556640625,
246  .399609375e-1,
247  0.,
248  .1198828125,
249  .3613975628004575124052940721184028345129,
250  0.,
251  -1.341524066700492771819987788202715834917,
252  1.370126503900035259414693716084313000404,
253  .490472027972027972027972027972027972028e-1,
254  0.,
255  0.,
256  .2350972042214404739862988335493427143122,
257  .180855592981356728810903963653454488485,
258  .6169289044289044289044289044289044289044e-1,
259  0.,
260  0.,
261  .1123656831464027662262557035130015442303,
262  -.3885046071451366767049048108111244567456e-1,
263  .1979188712522045855379188712522045855379e-1,
264  -1.767630240222326875735597119572145586714,
265  0.,
266  0.,
267  -62.5,
268  -6.061889377376669100821361459659331999758,
269  5.650823198222763138561298030600840174201,
270  65.62169641937623283799566054863063741227,
271  -1.180945066554970799825116282628297957882,
272  0.,
273  0.,
274  -41.50473441114320841606641502701994225874,
275  -4.434438319103725011225169229846100211776,
276  4.260408188586133024812193710744693240761,
277  43.75364022446171584987676829438379303004,
278  .787142548991231068744647504422630755086e-2,
279  -1.281405999441488405459510291182054246266,
280  0.,
281  0.,
282  -45.04713996013986630220754257136007322267,
283  -4.731362069449576477311464265491282810943,
284  4.514967016593807841185851584597240996214,
285  47.44909557172985134869022392235929015114,
286  .1059228297111661135687393955516542875228e-1,
287  -.5746842263844616254432318478286296232021e-2,
288  -1.724470134262485191756709817484481861731,
289  0.,
290  0.,
291  -60.92349008483054016518434619253765246063,
292  -5.95151837622239245520283276706185486829,
293  5.556523730698456235979791650843592496839,
294  63.98301198033305336837536378635995939281,
295  .1464202825041496159275921391759452676003e-1,
296  .6460408772358203603621865144977650714892e-1,
297  -.7930323169008878984024452548693373291447e-1,
298  -3.301622667747079016353994789790983625569,
299  0.,
300  0.,
301  -118.011272359752508566692330395789886851,
302  -10.14142238845611248642783916034510897595,
303  9.139311332232057923544012273556827000619,
304  123.3759428284042683684847180986501894364,
305  4.623244378874580474839807625067630924792,
306  -3.383277738068201923652550971536811240814,
307  4.527592100324618189451265339351129035325,
308  -5.828495485811622963193088019162985703755
309 };
310 const double RK8Solver::b[] =
311 {
312  .4427989419007951074716746668098518862111e-1,
313  0.,
314  0.,
315  0.,
316  0.,
317  .3541049391724448744815552028733568354121,
318  .2479692154956437828667629415370663023884,
319  -15.69420203883808405099207034271191213468,
320  25.08406496555856261343930031237186278518,
321  -31.73836778626027646833156112007297739997,
322  22.93828327398878395231483560344797018313,
323  -.2361324633071542145259900641263517600737
324 };
325 const double RK8Solver::c[] =
326 {
327  .5e-1,
328  .1065625,
329  .15984375,
330  .39,
331  .465,
332  .155,
333  .943,
334  .901802041735856958259707940678372149956,
335  .909,
336  .94,
337  1.,
338 };
339
340
342 {
343  ODESolver::Init(_f);
344  k.SetSize(f->Width());
345 }
346
347 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
348 {
349  f->SetTime(t + dt);
350  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
352  t += dt;
353 }
354
355
357 {
358  ODESolver::Init(_f);
359  k.SetSize(f->Width());
360 }
361
362 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
363 {
364  f->SetTime(t + dt/2);
365  f->ImplicitSolve(dt/2, x, k);
367  t += dt;
368 }
369
370
372 {
373  if (gamma_opt == 0)
374  {
375  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
376  }
377  else if (gamma_opt == 2)
378  {
379  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
380  }
381  else if (gamma_opt == 3)
382  {
383  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
384  }
385  else
386  {
387  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
388  }
389 }
390
392 {
393  ODESolver::Init(_f);
394  k.SetSize(f->Width());
395  y.SetSize(f->Width());
396 }
397
398 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
399 {
400  // with a = gamma:
401  // a | a
402  // 1-a | 1-2a a
403  // ------+-----------
404  // | 1/2 1/2
405  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
406  f->SetTime(t + gamma*dt);
407  f->ImplicitSolve(gamma*dt, x, k);
408  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
410
411  f->SetTime(t + (1.-gamma)*dt);
412  f->ImplicitSolve(gamma*dt, y, k);
414  t += dt;
415 }
416
417
419 {
420  ODESolver::Init(_f);
421  k.SetSize(f->Width());
422  y.SetSize(f->Width());
423  z.SetSize(f->Width());
424 }
425
426 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
427 {
428  // a | a
429  // 1/2 | 1/2-a a
430  // 1-a | 2a 1-4a a
431  // ------+--------------------
432  // | b 1-2b b
433  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
434  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
435  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
436
437  f->SetTime(t + a*dt);
438  f->ImplicitSolve(a*dt, x, k);
442
443  f->SetTime(t + dt/2);
444  f->ImplicitSolve(a*dt, y, k);
447
448  f->SetTime(t + (1.-a)*dt);
449  f->ImplicitSolve(a*dt, z, k);
451  t += dt;
452 }
453
454
456 {
457  ODESolver::Init(_f);
458  k.SetSize(f->Width());
459  y.SetSize(f->Width());
460 }
461
462 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
463 {
464  // a | a
465  // c | c-a a
466  // 1 | b 1-a-b a
467  // -----+----------------
468  // | b 1-a-b a
469  const double a = 0.435866521508458999416019;
470  const double b = 1.20849664917601007033648;
471  const double c = 0.717933260754229499708010;
472
473  f->SetTime(t + a*dt);
474  f->ImplicitSolve(a*dt, x, k);
477
478  f->SetTime(t + c*dt);
479  f->ImplicitSolve(a*dt, y, k);
481
482  f->SetTime(t + dt);
483  f->ImplicitSolve(a*dt, x, k);
485  t += dt;
486 }
487
488
490 {
491  ODESolver::Init(_f);
492  k.SetSize(f->Width());
493  y.SetSize(f->Width());
494  xdot.SetSize(f->Width());
495  xdot = 0.0;
496  first = true;
497 }
498
500 {
501  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
502  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
503
504  alpha_m = 0.5*(3.0 - rho_inf)/(1.0 + rho_inf);
505  alpha_f = 1.0/(1.0 + rho_inf);
506  gamma = 0.5 + alpha_m - alpha_f;
507 }
508
510 {
511  out << "Generalized alpha time integrator:" << std::endl;
512  out << "alpha_m = " << alpha_m << std::endl;
513  out << "alpha_f = " << alpha_f << std::endl;
514  out << "gamma = " << gamma << std::endl;
515
516  if (gamma == 0.5 + alpha_m - alpha_f)
517  {
518  out<<"Second order"<<" and ";
519  }
520  else
521  {
522  out<<"First order"<<" and ";
523  }
524
525  if ((alpha_m >= alpha_f)&&(alpha_f >= 0.5))
526  {
527  out<<"Stable"<<std::endl;
528  }
529  else
530  {
531  out<<"Unstable"<<std::endl;
532  }
533 }
534
535 // This routine assumes xdot is initialized.
536 void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
537 {
538  double dt_fac1 = alpha_f*(1.0 - gamma/alpha_m);
539  double dt_fac2 = alpha_f*gamma/alpha_m;
540  double dt_fac3 = 1.0/alpha_m;
541
542  // In the first pass xdot is not yet computed. If parameter choices requires
543  // xdot midpoint rule is used instead for the first step only.
544  if (first && (dt_fac1 != 0.0))
545  {
546  dt_fac1 = 0.0;
547  dt_fac2 = 0.5;
548  dt_fac3 = 2.0;
549  first = false;
550  }
551
553  f->SetTime(t + dt_fac2*dt);
554  f->ImplicitSolve(dt_fac2*dt, y, k);
555
559
560  t += dt;
561 }
562
563
564 void
566 {
567  P_ = &P; F_ = &F;
568
569  dp_.SetSize(F_->Height());
570  dq_.SetSize(P_->Height());
571 }
572
573 void
574 SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
575 {
576  F_->SetTime(t);
577  F_->Mult(q,dp_);
579
580  P_->Mult(p,dq_);
582
583  t += dt;
584 }
585
586 void
587 SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
588 {
589  P_->Mult(p,dq_);
591
592  F_->SetTime(t+0.5*dt);
593  F_->Mult(q,dp_);
595
596  P_->Mult(p,dq_);
598
599  t += dt;
600 }
601
603  : order_(order)
604 {
605  a_.SetSize(order);
606  b_.SetSize(order);
607
608  switch (order_)
609  {
610  case 1:
611  a_[0] = 1.0;
612  b_[0] = 1.0;
613  break;
614  case 2:
615  a_[0] = 0.5;
616  a_[1] = 0.5;
617  b_[0] = 0.0;
618  b_[1] = 1.0;
619  break;
620  case 3:
621  a_[0] = 2.0/3.0;
622  a_[1] = -2.0/3.0;
623  a_[2] = 1.0;
624  b_[0] = 7.0/24.0;
625  b_[1] = 0.75;
626  b_[2] = -1.0/24.0;
627  break;
628  case 4:
629  a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
630  a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
631  a_[2] = a_[1];
632  a_[3] = a_[0];
633  b_[0] = 0.0;
634  b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
635  b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
636  b_[3] = b_[1];
637  break;
638  default:
639  MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
640  };
641 }
642
643 void
644 SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
645 {
646  for (int i=0; i<order_; i++)
647  {
648  if ( b_[i] != 0.0 )
649  {
650  F_->SetTime(t);
651  if ( F_->isExplicit() )
652  {
653  F_->Mult(q, dp_);
654  }
655  else
656  {
657  F_->ImplicitSolve(b_[i] * dt, q, dp_);
658  }
660  }
661
662  P_->Mult(p, dq_);
664
665  t += a_[i] * dt;
666  }
667 }
668
669 }
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:33
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:509
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:565
void SetRhoInf(double rho_inf)
Definition: ode.cpp:499
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:574
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
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
Base abstract class for time dependent operators.
Definition: operator.hpp:151
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:181
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
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:260
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
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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
Vector dq_
Definition: ode.hpp:344
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:62
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:184
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: operator.hpp:234
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:489
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:418
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: operator.hpp:213
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
TimeDependentOperator * F_
Definition: ode.hpp:340
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
Vector dp_
Definition: ode.hpp:343
Operator * P_
Definition: ode.hpp:341
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:201
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:587
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:644
SIAVSolver(int order)
Definition: ode.cpp:602
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:536
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:48
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:455
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
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
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: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