MFEM  v4.1.0 Finite element discretization library
ode.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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);
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);
60
61  f->SetTime(t + a*dt);
62  f->Mult(x, dxdt);
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)
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)
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
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
122
123  f->SetTime(t + dt/2);
124  f->Mult(y, k); // k2
127
128  f->Mult(y, k); // k3
131
132  f->SetTime(t + dt);
133  f->Mult(y, k); // k4
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  {
174  for (int j = 1; j < i; j++)
175  {
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  {
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  s = 0;
350  smax = std::min(_s,5);
351  a = _a;
352  k = new Vector[5];
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  ODESolver::Init(_f);
371  RKsolver->Init(_f);
372  idx.SetSize(smax);
373  for (int i = 0; i < smax; i++)
374  {
375  idx[i] = (smax-i)%smax;
376  k[i].SetSize(f->Width());
377  }
378  s = 0;
379 }
380
381 void AdamsBashforthSolver::Step(Vector &x, double &t, double &dt)
382 {
383  s++;
384  s = std::min(s, smax);
385  if (s == smax)
386  {
387  f->SetTime(t);
388  f->Mult(x, k[idx[0]]);
389  for (int i = 0; i < s; i++)
390  {
392  }
393  }
394  else
395  {
396  f->Mult(x,k[idx[0]]);
397  RKsolver->Step(x,t,dt);
398  }
399  t += dt;
400
401  // Shift the index
402  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
403 }
404
405 const double AB1Solver::a[] =
406 {1.0};
407 const double AB2Solver::a[] =
408 {1.5,-0.5};
409 const double AB3Solver::a[] =
410 {23.0/12.0,-4.0/3.0, 5.0/12.0};
411 const double AB4Solver::a[] =
412 {55.0/24.0,-59.0/24.0, 37.0/24.0,-9.0/24.0};
413 const double AB5Solver::a[] =
414 {1901.0/720.0,-2774.0/720.0, 2616.0/720.0,-1274.0/720.0, 251.0/720.0};
415
417 {
418  s = 0;
419  smax = std::min(_s+1,5);
420  a = _a;
421  k = new Vector[5];
422
423  if (smax <= 3)
424  {
425  RKsolver = new SDIRK23Solver();
426  }
427  else
428  {
429  RKsolver = new SDIRK34Solver();
430  }
431 }
432
434 {
435  ODESolver::Init(_f);
436  RKsolver->Init(_f);
437  int n = f->Width();
438  idx.SetSize(smax);
439  for (int i = 0; i < smax; i++)
440  {
441  idx[i] = (smax-i)%smax;
442  k[i].SetSize(n);
443  }
444  s = 0;
445 }
446
447 void AdamsMoultonSolver::Step(Vector &x, double &t, double &dt)
448 {
449  if ((s == 0)&&(smax>1))
450  {
451  f->Mult(x,k[idx[1]]);
452  }
453  s++;
454  s = std::min(s, smax);
455
456  if (s >= smax-1)
457  {
458  f->SetTime(t);
459  for (int i = 1; i < smax; i++)
460  {
462  }
463  f->ImplicitSolve(a[0]*dt, x, k[idx[0]]);
465  }
466  else
467  {
468  RKsolver->Step(x,t,dt);
469  f->Mult(x,k[idx[0]]);
470  }
471  t += dt;
472
473  // Shift the index
474  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
475 }
476
477 const double AM0Solver::a[] =
478 {1.0};
479 const double AM1Solver::a[] =
480 {0.5, 0.5};
481 const double AM2Solver::a[] =
482 {5.0/12.0, 2.0/3.0, -1.0/12.0};
483 const double AM3Solver::a[] =
484 {3.0/8.0, 19.0/24.0,-5.0/24.0, 1.0/24.0};
485 const double AM4Solver::a[] =
486 {251.0/720.0,646.0/720.0,-264.0/720.0, 106.0/720.0, -19.0/720.0};
487
488
490 {
491  ODESolver::Init(_f);
492  k.SetSize(f->Width(), mem_type);
493 }
494
495 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
496 {
497  f->SetTime(t + dt);
498  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
500  t += dt;
501 }
502
503
505 {
506  ODESolver::Init(_f);
507  k.SetSize(f->Width(), mem_type);
508 }
509
510 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
511 {
512  f->SetTime(t + dt/2);
513  f->ImplicitSolve(dt/2, x, k);
515  t += dt;
516 }
517
518
520 {
521  if (gamma_opt == 0)
522  {
523  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
524  }
525  else if (gamma_opt == 2)
526  {
527  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
528  }
529  else if (gamma_opt == 3)
530  {
531  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
532  }
533  else
534  {
535  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
536  }
537 }
538
540 {
541  ODESolver::Init(_f);
542  k.SetSize(f->Width(), mem_type);
543  y.SetSize(f->Width(), mem_type);
544 }
545
546 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
547 {
548  // with a = gamma:
549  // a | a
550  // 1-a | 1-2a a
551  // ------+-----------
552  // | 1/2 1/2
553  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
554  f->SetTime(t + gamma*dt);
555  f->ImplicitSolve(gamma*dt, x, k);
556  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
558
559  f->SetTime(t + (1.-gamma)*dt);
560  f->ImplicitSolve(gamma*dt, y, k);
562  t += dt;
563 }
564
565
567 {
568  ODESolver::Init(_f);
569  k.SetSize(f->Width(), mem_type);
570  y.SetSize(f->Width(), mem_type);
571  z.SetSize(f->Width(), mem_type);
572 }
573
574 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
575 {
576  // a | a
577  // 1/2 | 1/2-a a
578  // 1-a | 2a 1-4a a
579  // ------+--------------------
580  // | b 1-2b b
581  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
582  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
583  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
584
585  f->SetTime(t + a*dt);
586  f->ImplicitSolve(a*dt, x, k);
590
591  f->SetTime(t + dt/2);
592  f->ImplicitSolve(a*dt, y, k);
595
596  f->SetTime(t + (1.-a)*dt);
597  f->ImplicitSolve(a*dt, z, k);
599  t += dt;
600 }
601
602
604 {
605  ODESolver::Init(_f);
606  k.SetSize(f->Width(), mem_type);
607  y.SetSize(f->Width(), mem_type);
608 }
609
610 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
611 {
612  // a | a
613  // c | c-a a
614  // 1 | b 1-a-b a
615  // -----+----------------
616  // | b 1-a-b a
617  const double a = 0.435866521508458999416019;
618  const double b = 1.20849664917601007033648;
619  const double c = 0.717933260754229499708010;
620
621  f->SetTime(t + a*dt);
622  f->ImplicitSolve(a*dt, x, k);
625
626  f->SetTime(t + c*dt);
627  f->ImplicitSolve(a*dt, y, k);
629
630  f->SetTime(t + dt);
631  f->ImplicitSolve(a*dt, x, k);
633  t += dt;
634 }
635
636
638 {
639  ODESolver::Init(_f);
640  k.SetSize(f->Width(), mem_type);
641  y.SetSize(f->Width(), mem_type);
642  xdot.SetSize(f->Width(), mem_type);
643  xdot = 0.0;
644  first = true;
645 }
646
648 {
649  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
650  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
651
652  // According to Jansen
653  alpha_m = 0.5*(3.0 - rho_inf)/(1.0 + rho_inf);
654  alpha_f = 1.0/(1.0 + rho_inf);
655  gamma = 0.5 + alpha_m - alpha_f;
656 }
657
659 {
660  out << "Generalized alpha time integrator:" << std::endl;
661  out << "alpha_m = " << alpha_m << std::endl;
662  out << "alpha_f = " << alpha_f << std::endl;
663  out << "gamma = " << gamma << std::endl;
664
665  if (gamma == 0.5 + alpha_m - alpha_f)
666  {
667  out<<"Second order"<<" and ";
668  }
669  else
670  {
671  out<<"First order"<<" and ";
672  }
673
674  if ((alpha_m >= alpha_f)&&(alpha_f >= 0.5))
675  {
676  out<<"Stable"<<std::endl;
677  }
678  else
679  {
680  out<<"Unstable"<<std::endl;
681  }
682 }
683
684 // This routine assumes xdot is initialized.
685 void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
686 {
687  if (first)
688  {
689  f->Mult(x,xdot);
690  first = false;
691  }
692
693  // Set y = x + alpha_f*(1.0 - (gamma/alpha_m))*dt*xdot
694  add(x, alpha_f*(1.0 - (gamma/alpha_m))*dt, xdot, y);
695
696  // Solve k = f(y + dt_eff*k)
697  double dt_eff = (gamma*alpha_f/alpha_m)*dt;
698  f->SetTime(t + alpha_f*dt);
699  f->ImplicitSolve(dt_eff, y, k);
700
701  // Update x and xdot
704
705  xdot *= (1.0-(1.0/alpha_m));
707
708  t += dt;
709 }
710
711
712 void
714 {
715  P_ = &P; F_ = &F;
716
717  dp_.SetSize(F_->Height());
718  dq_.SetSize(P_->Height());
719 }
720
721 void
722 SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
723 {
724  F_->SetTime(t);
725  F_->Mult(q,dp_);
727
728  P_->Mult(p,dq_);
730
731  t += dt;
732 }
733
734 void
735 SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
736 {
737  P_->Mult(p,dq_);
739
740  F_->SetTime(t+0.5*dt);
741  F_->Mult(q,dp_);
743
744  P_->Mult(p,dq_);
746
747  t += dt;
748 }
749
751  : order_(order)
752 {
753  a_.SetSize(order);
754  b_.SetSize(order);
755
756  switch (order_)
757  {
758  case 1:
759  a_[0] = 1.0;
760  b_[0] = 1.0;
761  break;
762  case 2:
763  a_[0] = 0.5;
764  a_[1] = 0.5;
765  b_[0] = 0.0;
766  b_[1] = 1.0;
767  break;
768  case 3:
769  a_[0] = 2.0/3.0;
770  a_[1] = -2.0/3.0;
771  a_[2] = 1.0;
772  b_[0] = 7.0/24.0;
773  b_[1] = 0.75;
774  b_[2] = -1.0/24.0;
775  break;
776  case 4:
777  a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
778  a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
779  a_[2] = a_[1];
780  a_[3] = a_[0];
781  b_[0] = 0.0;
782  b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
783  b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
784  b_[3] = b_[1];
785  break;
786  default:
787  MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
788  };
789 }
790
791 void
792 SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
793 {
794  for (int i=0; i<order_; i++)
795  {
796  if ( b_[i] != 0.0 )
797  {
798  F_->SetTime(t);
799  if ( F_->isExplicit() )
800  {
801  F_->Mult(q, dp_);
802  }
803  else
804  {
805  F_->ImplicitSolve(b_[i] * dt, q, dp_);
806  }
808  }
809
810  P_->Mult(p, dq_);
812
813  t += a_[i] * dt;
814  }
815 }
816
818 {
819  this->f = &f;
821 }
822
824 {
826  d2xdt2.SetSize(f->Width());
827  d2xdt2 = 0.0;
828  first = true;
829 }
830
832 {
833  out << "Newmark time integrator:" << std::endl;
834  out << "beta = " << beta << std::endl;
835  out << "gamma = " << gamma << std::endl;
836
837  if (gamma == 0.5)
838  {
839  out<<"Second order"<<" and ";
840  }
841  else
842  {
843  out<<"First order"<<" and ";
844  }
845
846  if ((gamma >= 0.5) && (beta >= (gamma + 0.5)*(gamma + 0.5)/4))
847  {
848  out<<"A-Stable"<<std::endl;
849  }
850  else if ((gamma >= 0.5) && (beta >= 0.5*gamma))
851  {
852  out<<"Conditionally stable"<<std::endl;
853  }
854  else
855  {
856  out<<"Unstable"<<std::endl;
857  }
858 }
859
860 void NewmarkSolver::Step(Vector &x, Vector &dxdt, double &t, double &dt)
861 {
862  double fac0 = 0.5 - beta;
863  double fac2 = 1.0 - gamma;
864  double fac3 = beta;
865  double fac4 = gamma;
866
867  // In the first pass compute d2xdt2 directy from operator.
868  if (first)
869  {
870  f->Mult(x, dxdt, d2xdt2);
871  first = false;
872  }
873  f->SetTime(t + dt);
874
878
879  f->SetTime(t + dt);
880  f->ImplicitSolve(fac3*dt*dt, fac4*dt, x, dxdt, d2xdt2);
881
884  t += dt;
885 }
886
888 {
890  xa.SetSize(f->Width());
891  va.SetSize(f->Width());
892  aa.SetSize(f->Width());
893  d2xdt2.SetSize(f->Width());
894  d2xdt2 = 0.0;
895  first = true;
896 }
897
899 {
900  out << "Generalized alpha time integrator:" << std::endl;
901  out << "alpha_m = " << alpha_m << std::endl;
902  out << "alpha_f = " << alpha_f << std::endl;
903  out << "beta = " << beta << std::endl;
904  out << "gamma = " << gamma << std::endl;
905
906  if (gamma == 0.5 + alpha_m - alpha_f)
907  {
908  out<<"Second order"<<" and ";
909  }
910  else
911  {
912  out<<"First order"<<" and ";
913  }
914
915  if ((alpha_m >= alpha_f)&&
916  (alpha_f >= 0.5) &&
917  (beta >= 0.25 + 0.5*(alpha_m - alpha_f)))
918  {
919  out<<"Stable"<<std::endl;
920  }
921  else
922  {
923  out<<"Unstable"<<std::endl;
924  }
925 }
926
928  double &t, double &dt)
929 {
930  double fac0 = (0.5 - (beta/alpha_m));
931  double fac1 = alpha_f;
932  double fac2 = alpha_f*(1.0 - (gamma/alpha_m));
933  double fac3 = beta*alpha_f/alpha_m;
934  double fac4 = gamma*alpha_f/alpha_m;
935  double fac5 = alpha_m;
936
937  // In the first pass compute d2xdt2 directy from operator.
938  if (first)
939  {
940  f->Mult(x, dxdt, d2xdt2);
941  first = false;
942  }
943
944
945  // Predict alpha levels
949
950  // Solve alpha levels
951  f->SetTime(t + dt);
952  f->ImplicitSolve(fac3*dt*dt, fac4*dt, xa, va, aa);
953
954  // Correct alpha levels
957
958  // Extrapolate
959  x *= 1.0 - 1.0/fac1;
961
962  dxdt *= 1.0 - 1.0/fac1;
964
965  d2xdt2 *= 1.0 - 1.0/fac5;
967
968  t += dt;
969 }
970
971 }
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:284
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:658
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:713
void SetRhoInf(double rho_inf)
Definition: ode.cpp:647
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:722
virtual void Step(Vector &x, Vector &dxdt, 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:927
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:539
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
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
Base abstract class for first order time dependent operators.
Definition: operator.hpp:259
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:230
MemoryType mem_type
Definition: ode.hpp:27
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
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
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:303
virtual void Step(Vector &x, Vector &dxdt, 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:860
virtual void Init(SecondOrderTimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:887
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:447
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
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:610
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:238
virtual void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:831
Definition: ode.cpp:416
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:574
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:504
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:77
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:368
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
double b
Definition: lissajous.cpp:42
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:492
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:817
virtual void Init(SecondOrderTimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:823
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:306
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:510
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:495
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:637
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:566
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:488
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:433
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:132
Base abstract class for second order time dependent operators.
Definition: operator.hpp:451
virtual void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:898
Vector dp_
Definition: ode.hpp:491
Operator * P_
Definition: ode.hpp:489
double a
Definition: lissajous.cpp:41
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:735
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:792
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:49
SIAVSolver(int order)
Definition: ode.cpp:750
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: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:685
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
Definition: ode.cpp:347
Vector data type.
Definition: vector.hpp:48
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:603
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:546
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:66
Abstract operator.
Definition: operator.hpp:24
virtual void ImplicitSolve(const double dt0, const double dt1, const Vector &x, const Vector &dxdt, Vector &k)
Solve the equation: k = f(x + 1/2 dt0^2 k, dxdt + dt1 k, t), for the unknown k at the current time t...
Definition: operator.cpp:291
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:519
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:381
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:532
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148