MFEM  v4.6.0
Finite element discretization library
ode.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "operator.hpp"
13 #include "ode.hpp"
14 
15 namespace mfem
16 {
17 
19 {
20  this->f = &f_;
22 }
23 
25 {
26  ODESolver::Init(f_);
27  dxdt.SetSize(f->Width(), mem_type);
28 }
29 
30 void ForwardEulerSolver::Step(Vector &x, double &t, double &dt)
31 {
32  f->SetTime(t);
33  f->Mult(x, dxdt);
34  x.Add(dt, dxdt);
35  t += dt;
36 }
37 
38 
40 {
41  ODESolver::Init(f_);
42  int n = f->Width();
43  dxdt.SetSize(n, mem_type);
44  x1.SetSize(n, mem_type);
45 }
46 
47 void RK2Solver::Step(Vector &x, double &t, double &dt)
48 {
49  // 0 |
50  // a | a
51  // ---+--------
52  // | 1-b b b = 1/(2a)
53 
54  const double b = 0.5/a;
55 
56  f->SetTime(t);
57  f->Mult(x, dxdt);
58  add(x, (1. - b)*dt, dxdt, x1);
59  x.Add(a*dt, dxdt);
60 
61  f->SetTime(t + a*dt);
62  f->Mult(x, dxdt);
63  add(x1, b*dt, dxdt, x);
64  t += dt;
65 }
66 
67 
69 {
70  ODESolver::Init(f_);
71  int n = f->Width();
72  y.SetSize(n, mem_type);
73  k.SetSize(n, mem_type);
74 }
75 
76 void RK3SSPSolver::Step(Vector &x, double &t, double &dt)
77 {
78  // x0 = x, t0 = t, k0 = dt*f(t0, x0)
79  f->SetTime(t);
80  f->Mult(x, k);
81 
82  // x1 = x + k0, t1 = t + dt, k1 = dt*f(t1, x1)
83  add(x, dt, k, y);
84  f->SetTime(t + dt);
85  f->Mult(y, k);
86 
87  // x2 = 3/4*x + 1/4*(x1 + k1), t2 = t + 1/2*dt, k2 = dt*f(t2, x2)
88  y.Add(dt, k);
89  add(3./4, x, 1./4, y, y);
90  f->SetTime(t + dt/2);
91  f->Mult(y, k);
92 
93  // x3 = 1/3*x + 2/3*(x2 + k2), t3 = t + dt
94  y.Add(dt, k);
95  add(1./3, x, 2./3, y, x);
96  t += dt;
97 }
98 
99 
101 {
102  ODESolver::Init(f_);
103  int n = f->Width();
104  y.SetSize(n, mem_type);
105  k.SetSize(n, mem_type);
106  z.SetSize(n, mem_type);
107 }
108 
109 void RK4Solver::Step(Vector &x, double &t, double &dt)
110 {
111  // 0 |
112  // 1/2 | 1/2
113  // 1/2 | 0 1/2
114  // 1 | 0 0 1
115  // -----+-------------------
116  // | 1/6 1/3 1/3 1/6
117 
118  f->SetTime(t);
119  f->Mult(x, k); // k1
120  add(x, dt/2, k, y);
121  add(x, dt/6, k, z);
122 
123  f->SetTime(t + dt/2);
124  f->Mult(y, k); // k2
125  add(x, dt/2, k, y);
126  z.Add(dt/3, k);
127 
128  f->Mult(y, k); // k3
129  add(x, dt, k, y);
130  z.Add(dt/3, k);
131 
132  f->SetTime(t + dt);
133  f->Mult(y, k); // k4
134  add(z, dt/6, k, x);
135  t += dt;
136 }
137 
138 ExplicitRKSolver::ExplicitRKSolver(int s_, const double *a_, const double *b_,
139  const double *c_)
140 {
141  s = s_;
142  a = a_;
143  b = b_;
144  c = c_;
145  k = new Vector[s];
146 }
147 
149 {
150  ODESolver::Init(f_);
151  int n = f->Width();
152  y.SetSize(n, mem_type);
153  for (int i = 0; i < s; i++)
154  {
155  k[i].SetSize(n, mem_type);
156  }
157 }
158 
159 void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
160 {
161  // 0 |
162  // c[0] | a[0]
163  // c[1] | a[1] a[2]
164  // ... | ...
165  // c[s-2] | ... a[s(s-1)/2-1]
166  // --------+---------------------
167  // | b[0] b[1] ... b[s-1]
168 
169  f->SetTime(t);
170  f->Mult(x, k[0]);
171  for (int l = 0, i = 1; i < s; i++)
172  {
173  add(x, a[l++]*dt, k[0], y);
174  for (int j = 1; j < i; j++)
175  {
176  y.Add(a[l++]*dt, k[j]);
177  }
178 
179  f->SetTime(t + c[i-1]*dt);
180  f->Mult(y, k[i]);
181  }
182  for (int i = 0; i < s; i++)
183  {
184  x.Add(b[i]*dt, k[i]);
185  }
186  t += dt;
187 }
188 
190 {
191  delete [] k;
192 }
193 
194 const double RK6Solver::a[] =
195 {
196  .6e-1,
197  .1923996296296296296296296296296296296296e-1,
198  .7669337037037037037037037037037037037037e-1,
199  .35975e-1,
200  0.,
201  .107925,
202  1.318683415233148260919747276431735612861,
203  0.,
204  -5.042058063628562225427761634715637693344,
205  4.220674648395413964508014358283902080483,
206  -41.87259166432751461803757780644346812905,
207  0.,
208  159.4325621631374917700365669070346830453,
209  -122.1192135650100309202516203389242140663,
210  5.531743066200053768252631238332999150076,
211  -54.43015693531650433250642051294142461271,
212  0.,
213  207.0672513650184644273657173866509835987,
214  -158.6108137845899991828742424365058599469,
215  6.991816585950242321992597280791793907096,
216  -.1859723106220323397765171799549294623692e-1,
217  -54.66374178728197680241215648050386959351,
218  0.,
219  207.9528062553893734515824816699834244238,
220  -159.2889574744995071508959805871426654216,
221  7.018743740796944434698170760964252490817,
222  -.1833878590504572306472782005141738268361e-1,
223  -.5119484997882099077875432497245168395840e-3
224 };
225 const double RK6Solver::b[] =
226 {
227  .3438957868357036009278820124728322386520e-1,
228  0.,
229  0.,
230  .2582624555633503404659558098586120858767,
231  .4209371189673537150642551514069801967032,
232  4.405396469669310170148836816197095664891,
233  -176.4831190242986576151740942499002125029,
234  172.3641334014150730294022582711902413315
235 };
236 const double RK6Solver::c[] =
237 {
238  .6e-1,
239  .9593333333333333333333333333333333333333e-1,
240  .1439,
241  .4973,
242  .9725,
243  .9995,
244  1.,
245 };
246 
247 const double RK8Solver::a[] =
248 {
249  .5e-1,
250  -.69931640625e-2,
251  .1135556640625,
252  .399609375e-1,
253  0.,
254  .1198828125,
255  .3613975628004575124052940721184028345129,
256  0.,
257  -1.341524066700492771819987788202715834917,
258  1.370126503900035259414693716084313000404,
259  .490472027972027972027972027972027972028e-1,
260  0.,
261  0.,
262  .2350972042214404739862988335493427143122,
263  .180855592981356728810903963653454488485,
264  .6169289044289044289044289044289044289044e-1,
265  0.,
266  0.,
267  .1123656831464027662262557035130015442303,
268  -.3885046071451366767049048108111244567456e-1,
269  .1979188712522045855379188712522045855379e-1,
270  -1.767630240222326875735597119572145586714,
271  0.,
272  0.,
273  -62.5,
274  -6.061889377376669100821361459659331999758,
275  5.650823198222763138561298030600840174201,
276  65.62169641937623283799566054863063741227,
277  -1.180945066554970799825116282628297957882,
278  0.,
279  0.,
280  -41.50473441114320841606641502701994225874,
281  -4.434438319103725011225169229846100211776,
282  4.260408188586133024812193710744693240761,
283  43.75364022446171584987676829438379303004,
284  .787142548991231068744647504422630755086e-2,
285  -1.281405999441488405459510291182054246266,
286  0.,
287  0.,
288  -45.04713996013986630220754257136007322267,
289  -4.731362069449576477311464265491282810943,
290  4.514967016593807841185851584597240996214,
291  47.44909557172985134869022392235929015114,
292  .1059228297111661135687393955516542875228e-1,
293  -.5746842263844616254432318478286296232021e-2,
294  -1.724470134262485191756709817484481861731,
295  0.,
296  0.,
297  -60.92349008483054016518434619253765246063,
298  -5.95151837622239245520283276706185486829,
299  5.556523730698456235979791650843592496839,
300  63.98301198033305336837536378635995939281,
301  .1464202825041496159275921391759452676003e-1,
302  .6460408772358203603621865144977650714892e-1,
303  -.7930323169008878984024452548693373291447e-1,
304  -3.301622667747079016353994789790983625569,
305  0.,
306  0.,
307  -118.011272359752508566692330395789886851,
308  -10.14142238845611248642783916034510897595,
309  9.139311332232057923544012273556827000619,
310  123.3759428284042683684847180986501894364,
311  4.623244378874580474839807625067630924792,
312  -3.383277738068201923652550971536811240814,
313  4.527592100324618189451265339351129035325,
314  -5.828495485811622963193088019162985703755
315 };
316 const double RK8Solver::b[] =
317 {
318  .4427989419007951074716746668098518862111e-1,
319  0.,
320  0.,
321  0.,
322  0.,
323  .3541049391724448744815552028733568354121,
324  .2479692154956437828667629415370663023884,
325  -15.69420203883808405099207034271191213468,
326  25.08406496555856261343930031237186278518,
327  -31.73836778626027646833156112007297739997,
328  22.93828327398878395231483560344797018313,
329  -.2361324633071542145259900641263517600737
330 };
331 const double RK8Solver::c[] =
332 {
333  .5e-1,
334  .1065625,
335  .15984375,
336  .39,
337  .465,
338  .155,
339  .943,
340  .901802041735856958259707940678372149956,
341  .909,
342  .94,
343  1.,
344 };
345 
346 
348 {
349  smax = std::min(s_,5);
350  a = a_;
351  k = new Vector[5];
352  dt_ = -1.0;
353 
354  if (smax <= 2)
355  {
356  RKsolver = new RK2Solver();
357  }
358  else if (smax == 3)
359  {
360  RKsolver = new RK3SSPSolver();
361  }
362  else
363  {
364  RKsolver = new RK4Solver();
365  }
366 }
367 
369 {
370  MFEM_ASSERT( (i >= 0) && ( i < s ),
371  " AdamsBashforthSolver::GetStateVector \n" <<
372  " - Tried to get non-existent state "<<i);
373 
374  state = k[idx[i]];
375 }
376 
378 {
379  MFEM_ASSERT( (i >= 0) && ( i < s ),
380  " AdamsBashforthSolver::GetStateVector \n" <<
381  " - Tried to get non-existent state "<<i);
382 
383  return k[idx[i]];
384 }
385 
386 
388 {
389  MFEM_ASSERT( (i >= 0) && ( i < smax ),
390  " AdamsBashforthSolver::SetStateVector \n" <<
391  " - Tried to set non-existent state "<<i);
392  k[idx[i]] = state;
393  s = std::max(i,s);
394 }
395 
397 {
398  ODESolver::Init(f_);
399  RKsolver->Init(f_);
400  idx.SetSize(smax);
401  for (int i = 0; i < smax; i++)
402  {
403  idx[i] = (smax-i)%smax;
404  k[i].SetSize(f->Width());
405  }
406  s = 0;
407 }
408 
409 void AdamsBashforthSolver::Step(Vector &x, double &t, double &dt)
410 {
411  if ( (dt_ > 0.0) && (fabs(dt-dt_) >10*std::numeric_limits<double>::epsilon()))
412  {
413  s = 0;
414  dt_ = dt;
415 
416  if (print())
417  {
418  mfem::out << "WARNING:" << std::endl;
419  mfem::out << " - Time step changed" << std::endl;
420  mfem::out << " - Purging Adams-Bashforth history" << std::endl;
421  mfem::out << " - Will run Runge-Kutta to rebuild history" << std::endl;
422  }
423  }
424 
425  s++;
426  s = std::min(s, smax);
427  if (s == smax)
428  {
429  f->SetTime(t);
430  f->Mult(x, k[idx[0]]);
431  for (int i = 0; i < s; i++)
432  {
433  x.Add(a[i]*dt, k[idx[i]]);
434  }
435  t += dt;
436  }
437  else
438  {
439  f->Mult(x,k[idx[0]]);
440  RKsolver->Step(x,t,dt);
441  }
442 
443  // Shift the index
444  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
445 }
446 
447 const double AB1Solver::a[] =
448 {1.0};
449 const double AB2Solver::a[] =
450 {1.5,-0.5};
451 const double AB3Solver::a[] =
452 {23.0/12.0,-4.0/3.0, 5.0/12.0};
453 const double AB4Solver::a[] =
454 {55.0/24.0,-59.0/24.0, 37.0/24.0,-9.0/24.0};
455 const double AB5Solver::a[] =
456 {1901.0/720.0,-2774.0/720.0, 2616.0/720.0,-1274.0/720.0, 251.0/720.0};
457 
458 AdamsMoultonSolver::AdamsMoultonSolver(int s_, const double *a_)
459 {
460  s = 0;
461  smax = std::min(s_+1,5);
462  a = a_;
463  k = new Vector[5];
464  dt_ = -1.0;
465 
466  if (smax <= 3)
467  {
468  RKsolver = new SDIRK23Solver();
469  }
470  else
471  {
472  RKsolver = new SDIRK34Solver();
473  }
474 }
475 
477 {
478  MFEM_ASSERT( (i >= 0) && ( i < s ),
479  " AdamsMoultonSolver::GetStateVector \n" <<
480  " - Tried to get non-existent state "<<i);
481  return k[idx[i+1]];
482 }
483 
485 {
486  MFEM_ASSERT( (i >= 0) && ( i < s ),
487  " AdamsMoultonSolver::GetStateVector \n" <<
488  " - Tried to get non-existent state "<<i);
489  state = k[idx[i+1]];
490 }
491 
493 {
494  MFEM_ASSERT( (i >= 0) && ( i < smax ),
495  " AdamsMoultonSolver::SetStateVector \n" <<
496  " - Tried to set non-existent state "<<i);
497  k[idx[i+1]] = state;
498  s = std::max(i,s);
499 }
500 
502 {
503  ODESolver::Init(f_);
504  RKsolver->Init(f_);
505  int n = f->Width();
506  idx.SetSize(smax);
507  for (int i = 0; i < smax; i++)
508  {
509  idx[i] = (smax-i)%smax;
510  k[i].SetSize(n);
511  }
512  s = 0;
513 }
514 
515 void AdamsMoultonSolver::Step(Vector &x, double &t, double &dt)
516 {
517  if ( (dt_ > 0.0) && (fabs(dt-dt_) >10*std::numeric_limits<double>::epsilon()))
518  {
519  s = 0;
520  dt_ = dt;
521 
522  if (print())
523  {
524  mfem::out << "WARNING:" << std::endl;
525  mfem::out << " - Time step changed" << std::endl;
526  mfem::out << " - Purging Adams-Moulton history" << std::endl;
527  mfem::out << " - Will run Runge-Kutta to rebuild history" << std::endl;
528  }
529  }
530 
531  if ((s == 0)&&(smax>1))
532  {
533  f->Mult(x,k[idx[1]]);
534  }
535  s++;
536  s = std::min(s, smax);
537 
538  if (s >= smax-1)
539  {
540  f->SetTime(t);
541  for (int i = 1; i < smax; i++)
542  {
543  x.Add(a[i]*dt, k[idx[i]]);
544  }
545  f->ImplicitSolve(a[0]*dt, x, k[idx[0]]);
546  x.Add(a[0]*dt, k[idx[0]]);
547  t += dt;
548  }
549  else
550  {
551  RKsolver->Step(x,t,dt);
552  f->Mult(x,k[idx[0]]);
553  }
554 
555 
556  // Shift the index
557  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
558 }
559 
560 const double AM0Solver::a[] =
561 {1.0};
562 const double AM1Solver::a[] =
563 {0.5, 0.5};
564 const double AM2Solver::a[] =
565 {5.0/12.0, 2.0/3.0, -1.0/12.0};
566 const double AM3Solver::a[] =
567 {3.0/8.0, 19.0/24.0,-5.0/24.0, 1.0/24.0};
568 const double AM4Solver::a[] =
569 {251.0/720.0,646.0/720.0,-264.0/720.0, 106.0/720.0, -19.0/720.0};
570 
571 
573 {
574  ODESolver::Init(f_);
575  k.SetSize(f->Width(), mem_type);
576 }
577 
578 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
579 {
580  f->SetTime(t + dt);
581  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
582  x.Add(dt, k);
583  t += dt;
584 }
585 
586 
588 {
589  ODESolver::Init(f_);
590  k.SetSize(f->Width(), mem_type);
591 }
592 
593 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
594 {
595  f->SetTime(t + dt/2);
596  f->ImplicitSolve(dt/2, x, k);
597  x.Add(dt, k);
598  t += dt;
599 }
600 
601 
603 {
604  if (gamma_opt == 0)
605  {
606  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
607  }
608  else if (gamma_opt == 2)
609  {
610  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
611  }
612  else if (gamma_opt == 3)
613  {
614  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
615  }
616  else
617  {
618  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
619  }
620 }
621 
623 {
624  ODESolver::Init(f_);
625  k.SetSize(f->Width(), mem_type);
626  y.SetSize(f->Width(), mem_type);
627 }
628 
629 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
630 {
631  // with a = gamma:
632  // a | a
633  // 1-a | 1-2a a
634  // ------+-----------
635  // | 1/2 1/2
636  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
637  f->SetTime(t + gamma*dt);
638  f->ImplicitSolve(gamma*dt, x, k);
639  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
640  x.Add(dt/2, k);
641 
642  f->SetTime(t + (1.-gamma)*dt);
643  f->ImplicitSolve(gamma*dt, y, k);
644  x.Add(dt/2, k);
645  t += dt;
646 }
647 
648 
650 {
651  ODESolver::Init(f_);
652  k.SetSize(f->Width(), mem_type);
653  y.SetSize(f->Width(), mem_type);
654  z.SetSize(f->Width(), mem_type);
655 }
656 
657 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
658 {
659  // a | a
660  // 1/2 | 1/2-a a
661  // 1-a | 2a 1-4a a
662  // ------+--------------------
663  // | b 1-2b b
664  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
665  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
666  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
667 
668  f->SetTime(t + a*dt);
669  f->ImplicitSolve(a*dt, x, k);
670  add(x, (0.5-a)*dt, k, y);
671  add(x, (2.*a)*dt, k, z);
672  x.Add(b*dt, k);
673 
674  f->SetTime(t + dt/2);
675  f->ImplicitSolve(a*dt, y, k);
676  z.Add((1.-4.*a)*dt, k);
677  x.Add((1.-2.*b)*dt, k);
678 
679  f->SetTime(t + (1.-a)*dt);
680  f->ImplicitSolve(a*dt, z, k);
681  x.Add(b*dt, k);
682  t += dt;
683 }
684 
685 
687 {
688  ODESolver::Init(f_);
689  k.SetSize(f->Width(), mem_type);
690  y.SetSize(f->Width(), mem_type);
691 }
692 
693 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
694 {
695  // a | a
696  // c | c-a a
697  // 1 | b 1-a-b a
698  // -----+----------------
699  // | b 1-a-b a
700  const double a = 0.435866521508458999416019;
701  const double b = 1.20849664917601007033648;
702  const double c = 0.717933260754229499708010;
703 
704  f->SetTime(t + a*dt);
705  f->ImplicitSolve(a*dt, x, k);
706  add(x, (c-a)*dt, k, y);
707  x.Add(b*dt, k);
708 
709  f->SetTime(t + c*dt);
710  f->ImplicitSolve(a*dt, y, k);
711  x.Add((1.0-a-b)*dt, k);
712 
713  f->SetTime(t + dt);
714  f->ImplicitSolve(a*dt, x, k);
715  x.Add(a*dt, k);
716  t += dt;
717 }
718 
720 {
721  ODESolver::Init(f_);
722  k.SetSize(f->Width(), mem_type);
723  y.SetSize(f->Width(), mem_type);
724 }
725 
726 void TrapezoidalRuleSolver::Step(Vector &x, double &t, double &dt)
727 {
728  // 0 | 0 0
729  // 1 | 1/2 1/2
730  // ------+-----------
731  // | 1/2 1/2
732  f->SetTime(t);
733  f->Mult(x,k);
734  add(x, dt/2.0, k, y);
735  x.Add(dt/2.0, k);
736 
737  f->SetTime(t + dt);
738  f->ImplicitSolve(dt/2.0, y, k);
739  x.Add(dt/2.0, k);
740  t += dt;
741 }
742 
744 {
745  ODESolver::Init(f_);
746  k.SetSize(f->Width(), mem_type);
747  y.SetSize(f->Width(), mem_type);
748  z.SetSize(f->Width(), mem_type);
749 }
750 
751 void ESDIRK32Solver::Step(Vector &x, double &t, double &dt)
752 {
753  // 0 | 0 0 0
754  // 2a | a a 0
755  // 1 | 1-b-a b a
756  // ------+--------------------
757  // | 1-b-a b a
758  const double a = (2.0 - sqrt(2.0)) / 2.0;
759  const double b = (1.0 - 2.0*a) / (4.0*a);
760 
761  f->SetTime(t);
762  f->Mult(x,k);
763  add(x, a*dt, k, y);
764  add(x, (1.0-b-a)*dt, k, z);
765  x.Add((1.0-b-a)*dt, k);
766 
767  f->SetTime(t + (2.0*a)*dt);
768  f->ImplicitSolve(a*dt, y, k);
769  z.Add(b*dt, k);
770  x.Add(b*dt, k);
771 
772  f->SetTime(t + dt);
773  f->ImplicitSolve(a*dt, z, k);
774  x.Add(a*dt, k);
775  t += dt;
776 }
777 
779 {
780  ODESolver::Init(f_);
781  k.SetSize(f->Width(), mem_type);
782  y.SetSize(f->Width(), mem_type);
783  z.SetSize(f->Width(), mem_type);
784 }
785 
786 void ESDIRK33Solver::Step(Vector &x, double &t, double &dt)
787 {
788  // 0 | 0 0 0
789  // 2a | a a 0
790  // 1 | 1-b-a b a
791  // ------+----------------------------
792  // | 1-b_2-b_3 b_2 b_3
793  const double a = (3.0 + sqrt(3.0)) / 6.0;
794  const double b = (1.0 - 2.0*a) / (4.0*a);
795  const double b_2 = 1.0 / ( 12.0*a*(1.0 - 2.0*a) );
796  const double b_3 = (1.0 - 3.0*a) / ( 3.0*(1.0 - 2.0*a) );
797 
798  f->SetTime(t);
799  f->Mult(x,k);
800  add(x, a*dt, k, y);
801  add(x, (1.0-b-a)*dt, k, z);
802  x.Add((1.0-b_2-b_3)*dt, k);
803 
804  f->SetTime(t + (2.0*a)*dt);
805  f->ImplicitSolve(a*dt, y, k);
806  z.Add(b*dt, k);
807  x.Add(b_2*dt, k);
808 
809  f->SetTime(t + dt);
810  f->ImplicitSolve(a*dt, z, k);
811  x.Add(b_3*dt, k);
812  t += dt;
813 }
814 
816 {
817  ODESolver::Init(f_);
818  k.SetSize(f->Width(), mem_type);
819  y.SetSize(f->Width(), mem_type);
820  xdot.SetSize(f->Width(), mem_type);
821  xdot = 0.0;
822  nstate = 0;
823 }
824 
826 {
827  MFEM_ASSERT( (i == 0) && (nstate == 1),
828  "GeneralizedAlphaSolver::GetStateVector \n" <<
829  " - Tried to get non-existent state "<<i);
830  return xdot;
831 }
832 
834 {
835  MFEM_ASSERT( (i == 0) && (nstate == 1),
836  "GeneralizedAlphaSolver::GetStateVector \n" <<
837  " - Tried to get non-existent state "<<i);
838  state = xdot;
839 }
840 
842 {
843  MFEM_ASSERT( (i == 0),
844  "GeneralizedAlphaSolver::SetStateVector \n" <<
845  " - Tried to set non-existent state "<<i);
846  xdot = state;
847  nstate = 1;
848 }
849 
851 {
852  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
853  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
854 
855  // According to Jansen
856  alpha_m = 0.5*(3.0 - rho_inf)/(1.0 + rho_inf);
857  alpha_f = 1.0/(1.0 + rho_inf);
858  gamma = 0.5 + alpha_m - alpha_f;
859 }
860 
862 {
863  os << "Generalized alpha time integrator:" << std::endl;
864  os << "alpha_m = " << alpha_m << std::endl;
865  os << "alpha_f = " << alpha_f << std::endl;
866  os << "gamma = " << gamma << std::endl;
867 
868  if (gamma == 0.5 + alpha_m - alpha_f)
869  {
870  os<<"Second order"<<" and ";
871  }
872  else
873  {
874  os<<"First order"<<" and ";
875  }
876 
877  if ((alpha_m >= alpha_f)&&(alpha_f >= 0.5))
878  {
879  os<<"Stable"<<std::endl;
880  }
881  else
882  {
883  os<<"Unstable"<<std::endl;
884  }
885 }
886 
887 // This routine assumes xdot is initialized.
888 void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
889 {
890  if (nstate == 0)
891  {
892  f->Mult(x,xdot);
893  nstate = 1;
894  }
895 
896  // Set y = x + alpha_f*(1.0 - (gamma/alpha_m))*dt*xdot
897  add(x, alpha_f*(1.0 - (gamma/alpha_m))*dt, xdot, y);
898 
899  // Solve k = f(y + dt_eff*k)
900  double dt_eff = (gamma*alpha_f/alpha_m)*dt;
901  f->SetTime(t + alpha_f*dt);
902  f->ImplicitSolve(dt_eff, y, k);
903 
904  // Update x and xdot
905  x.Add((1.0 - (gamma/alpha_m))*dt, xdot);
906  x.Add( (gamma/alpha_m) *dt, k);
907 
908  xdot *= (1.0-(1.0/alpha_m));
909  xdot.Add((1.0/alpha_m),k);
910 
911  t += dt;
912 }
913 
914 
915 void
917 {
918  P_ = &P; F_ = &F;
919 
920  dp_.SetSize(F_->Height());
921  dq_.SetSize(P_->Height());
922 }
923 
924 void
925 SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
926 {
927  F_->SetTime(t);
928  F_->Mult(q,dp_);
929  p.Add(dt,dp_);
930 
931  P_->Mult(p,dq_);
932  q.Add(dt,dq_);
933 
934  t += dt;
935 }
936 
937 void
938 SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
939 {
940  P_->Mult(p,dq_);
941  q.Add(0.5*dt,dq_);
942 
943  F_->SetTime(t+0.5*dt);
944  F_->Mult(q,dp_);
945  p.Add(dt,dp_);
946 
947  P_->Mult(p,dq_);
948  q.Add(0.5*dt,dq_);
949 
950  t += dt;
951 }
952 
954  : order_(order)
955 {
956  a_.SetSize(order);
957  b_.SetSize(order);
958 
959  switch (order_)
960  {
961  case 1:
962  a_[0] = 1.0;
963  b_[0] = 1.0;
964  break;
965  case 2:
966  a_[0] = 0.5;
967  a_[1] = 0.5;
968  b_[0] = 0.0;
969  b_[1] = 1.0;
970  break;
971  case 3:
972  a_[0] = 2.0/3.0;
973  a_[1] = -2.0/3.0;
974  a_[2] = 1.0;
975  b_[0] = 7.0/24.0;
976  b_[1] = 0.75;
977  b_[2] = -1.0/24.0;
978  break;
979  case 4:
980  a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
981  a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
982  a_[2] = a_[1];
983  a_[3] = a_[0];
984  b_[0] = 0.0;
985  b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
986  b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
987  b_[3] = b_[1];
988  break;
989  default:
990  MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
991  };
992 }
993 
994 void
995 SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
996 {
997  for (int i=0; i<order_; i++)
998  {
999  if ( b_[i] != 0.0 )
1000  {
1001  F_->SetTime(t);
1002  if ( F_->isExplicit() )
1003  {
1004  F_->Mult(q, dp_);
1005  }
1006  else
1007  {
1008  F_->ImplicitSolve(b_[i] * dt, q, dp_);
1009  }
1010  p.Add(b_[i] * dt, dp_);
1011  }
1012 
1013  P_->Mult(p, dq_);
1014  q.Add(a_[i] * dt, dq_);
1015 
1016  t += a_[i] * dt;
1017  }
1018 }
1019 
1021 {
1022  this->f = &f_;
1024 }
1025 
1027 {
1029  d2xdt2.SetSize(f->Width());
1030  d2xdt2 = 0.0;
1031  first = true;
1032 }
1033 
1034 void NewmarkSolver::PrintProperties(std::ostream &os)
1035 {
1036  os << "Newmark time integrator:" << std::endl;
1037  os << "beta = " << beta << std::endl;
1038  os << "gamma = " << gamma << std::endl;
1039 
1040  if (gamma == 0.5)
1041  {
1042  os<<"Second order"<<" and ";
1043  }
1044  else
1045  {
1046  os<<"First order"<<" and ";
1047  }
1048 
1049  if ((gamma >= 0.5) && (beta >= (gamma + 0.5)*(gamma + 0.5)/4))
1050  {
1051  os<<"A-Stable"<<std::endl;
1052  }
1053  else if ((gamma >= 0.5) && (beta >= 0.5*gamma))
1054  {
1055  os<<"Conditionally stable"<<std::endl;
1056  }
1057  else
1058  {
1059  os<<"Unstable"<<std::endl;
1060  }
1061 }
1062 
1063 void NewmarkSolver::Step(Vector &x, Vector &dxdt, double &t, double &dt)
1064 {
1065  double fac0 = 0.5 - beta;
1066  double fac2 = 1.0 - gamma;
1067  double fac3 = beta;
1068  double fac4 = gamma;
1069 
1070  // In the first pass compute d2xdt2 directly from operator.
1071  if (first)
1072  {
1073  f->Mult(x, dxdt, d2xdt2);
1074  first = false;
1075  }
1076  f->SetTime(t + dt);
1077 
1078  x.Add(dt, dxdt);
1079  x.Add(fac0*dt*dt, d2xdt2);
1080  dxdt.Add(fac2*dt, d2xdt2);
1081 
1082  f->SetTime(t + dt);
1083  f->ImplicitSolve(fac3*dt*dt, fac4*dt, x, dxdt, d2xdt2);
1084 
1085  x .Add(fac3*dt*dt, d2xdt2);
1086  dxdt.Add(fac4*dt, d2xdt2);
1087  t += dt;
1088 }
1089 
1091 {
1093  xa.SetSize(f->Width());
1094  va.SetSize(f->Width());
1095  aa.SetSize(f->Width());
1096  d2xdt2.SetSize(f->Width());
1097  d2xdt2 = 0.0;
1098  nstate = 0;
1099 }
1100 
1102 {
1103  MFEM_ASSERT( (i == 0) && (nstate == 1),
1104  "GeneralizedAlpha2Solver::GetStateVector \n" <<
1105  " - Tried to get non-existent state "<<i);
1106  return d2xdt2;
1107 }
1108 
1109 
1111 {
1112  MFEM_ASSERT( (i == 0) && (nstate == 1),
1113  "GeneralizedAlpha2Solver::GetStateVector \n" <<
1114  " - Tried to get non-existent state "<<i);
1115  state = d2xdt2;
1116 }
1117 
1119 {
1120  MFEM_ASSERT( (i == 0),
1121  "GeneralizedAlpha2Solver::SetStateVector \n" <<
1122  " - Tried to set non-existent state "<<i);
1123  d2xdt2 = state;
1124  nstate = 1;
1125 }
1126 
1128 {
1129  os << "Generalized alpha time integrator:" << std::endl;
1130  os << "alpha_m = " << alpha_m << std::endl;
1131  os << "alpha_f = " << alpha_f << std::endl;
1132  os << "beta = " << beta << std::endl;
1133  os << "gamma = " << gamma << std::endl;
1134 
1135  if (gamma == 0.5 + alpha_m - alpha_f)
1136  {
1137  os<<"Second order"<<" and ";
1138  }
1139  else
1140  {
1141  os<<"First order"<<" and ";
1142  }
1143 
1144  if ((alpha_m >= alpha_f)&&
1145  (alpha_f >= 0.5) &&
1146  (beta >= 0.25 + 0.5*(alpha_m - alpha_f)))
1147  {
1148  os<<"Stable"<<std::endl;
1149  }
1150  else
1151  {
1152  os<<"Unstable"<<std::endl;
1153  }
1154 }
1155 
1157  double &t, double &dt)
1158 {
1159  double fac0 = (0.5 - (beta/alpha_m));
1160  double fac1 = alpha_f;
1161  double fac2 = alpha_f*(1.0 - (gamma/alpha_m));
1162  double fac3 = beta*alpha_f/alpha_m;
1163  double fac4 = gamma*alpha_f/alpha_m;
1164  double fac5 = alpha_m;
1165 
1166  // In the first pass compute d2xdt2 directly from operator.
1167  if (nstate == 0)
1168  {
1169  f->Mult(x, dxdt, d2xdt2);
1170  nstate = 1;
1171  }
1172 
1173  // Predict alpha levels
1174  add(dxdt, fac0*dt, d2xdt2, va);
1175  add(x, fac1*dt, va, xa);
1176  add(dxdt, fac2*dt, d2xdt2, va);
1177 
1178  // Solve alpha levels
1179  f->SetTime(t + dt);
1180  f->ImplicitSolve(fac3*dt*dt, fac4*dt, xa, va, aa);
1181 
1182  // Correct alpha levels
1183  xa.Add(fac3*dt*dt, aa);
1184  va.Add(fac4*dt, aa);
1185 
1186  // Extrapolate
1187  x *= 1.0 - 1.0/fac1;
1188  x.Add (1.0/fac1, xa);
1189 
1190  dxdt *= 1.0 - 1.0/fac1;
1191  dxdt.Add (1.0/fac1, va);
1192 
1193  d2xdt2 *= 1.0 - 1.0/fac5;
1194  d2xdt2.Add (1.0/fac5, aa);
1195 
1196  t += dt;
1197 }
1198 
1199 }
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:629
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1090
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:751
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:409
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:861
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:916
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:47
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:593
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:719
void SetRhoInf(double rho_inf)
Definition: ode.cpp:850
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:572
virtual void ImplicitSolve(const double fac0, const double fac1, const Vector &x, const Vector &dxdt, Vector &k)
Solve the equation: k = f(x + fac0 k, dxdt + fac1 k, t), for the unknown k at the current time t...
Definition: operator.cpp:359
double epsilon
Definition: ex25.cpp:141
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
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.cpp:298
MemoryType mem_type
Definition: ode.hpp:28
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:726
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
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]...
void Step(Vector &x, Vector &dxdt, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:1156
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:938
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
ExplicitRKSolver(int s_, const double *a_, const double *b_, const double *c_)
Definition: ode.cpp:138
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:360
const Vector & GetStateVector(int i) override
Definition: ode.cpp:377
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:686
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:778
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:995
void Step(Vector &x, Vector &dxdt, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:1063
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:815
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.cpp:293
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1034
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:1118
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
double b
Definition: lissajous.cpp:42
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:649
Vector dq_
Definition: ode.hpp:591
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1020
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:492
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:159
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:888
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:786
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:30
const Vector & GetStateVector(int i) override
Definition: ode.cpp:825
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:657
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:743
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
TimeDependentOperator * F_
Definition: ode.hpp:587
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:587
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:66
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:150
Base abstract class for second order time dependent operators.
Definition: operator.hpp:634
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1127
Vector dp_
Definition: ode.hpp:590
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:622
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1026
AdamsBashforthSolver(int s_, const double *a_)
Definition: ode.cpp:347
Operator * P_
Definition: ode.hpp:588
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:387
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const
Perform the action of the operator: y = k = f(x,@ dxdt, t), where k solves the algebraic equation F(x...
Definition: operator.cpp:352
double a
Definition: lissajous.cpp:41
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:501
const Vector & GetStateVector(int i) override
Definition: ode.cpp:1101
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:363
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:60
SIAVSolver(int order)
Definition: ode.cpp:953
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:515
const Vector & GetStateVector(int i) override
Definition: ode.cpp:476
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:693
RefCoord t[3]
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:578
Vector data type.
Definition: vector.hpp:58
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:396
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:109
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:27
RefCoord s[3]
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:841
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:925
Abstract operator.
Definition: operator.hpp:24
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:602
void Step(Vector &x, double &t, double &dt) override
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 Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
AdamsMoultonSolver(int s_, const double *a_)
Definition: ode.cpp:458
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:631