MFEM  v4.5.2
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 
353  if (smax <= 2)
354  {
355  RKsolver = new RK2Solver();
356  }
357  else if (smax == 3)
358  {
359  RKsolver = new RK3SSPSolver();
360  }
361  else
362  {
363  RKsolver = new RK4Solver();
364  }
365 }
366 
368 {
369  MFEM_ASSERT( (i >= 0) && ( i < s ),
370  " AdamsBashforthSolver::GetStateVector \n" <<
371  " - Tried to get non-existent state "<<i);
372 
373  state = k[idx[i]];
374 }
375 
377 {
378  MFEM_ASSERT( (i >= 0) && ( i < s ),
379  " AdamsBashforthSolver::GetStateVector \n" <<
380  " - Tried to get non-existent state "<<i);
381 
382  return k[idx[i]];
383 }
384 
385 
387 {
388  MFEM_ASSERT( (i >= 0) && ( i < smax ),
389  " AdamsBashforthSolver::SetStateVector \n" <<
390  " - Tried to set non-existent state "<<i);
391  k[idx[i]] = state;
392  s = std::max(i,s);
393 }
394 
396 {
397  ODESolver::Init(f_);
398  RKsolver->Init(f_);
399  idx.SetSize(smax);
400  for (int i = 0; i < smax; i++)
401  {
402  idx[i] = (smax-i)%smax;
403  k[i].SetSize(f->Width());
404  }
405  s = 0;
406 }
407 
408 void AdamsBashforthSolver::Step(Vector &x, double &t, double &dt)
409 {
410  s++;
411  s = std::min(s, smax);
412  if (s == smax)
413  {
414  f->SetTime(t);
415  f->Mult(x, k[idx[0]]);
416  for (int i = 0; i < s; i++)
417  {
418  x.Add(a[i]*dt, k[idx[i]]);
419  }
420  }
421  else
422  {
423  f->Mult(x,k[idx[0]]);
424  RKsolver->Step(x,t,dt);
425  }
426  t += dt;
427 
428  // Shift the index
429  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
430 }
431 
432 const double AB1Solver::a[] =
433 {1.0};
434 const double AB2Solver::a[] =
435 {1.5,-0.5};
436 const double AB3Solver::a[] =
437 {23.0/12.0,-4.0/3.0, 5.0/12.0};
438 const double AB4Solver::a[] =
439 {55.0/24.0,-59.0/24.0, 37.0/24.0,-9.0/24.0};
440 const double AB5Solver::a[] =
441 {1901.0/720.0,-2774.0/720.0, 2616.0/720.0,-1274.0/720.0, 251.0/720.0};
442 
443 AdamsMoultonSolver::AdamsMoultonSolver(int s_, const double *a_)
444 {
445  s = 0;
446  smax = std::min(s_+1,5);
447  a = a_;
448  k = new Vector[5];
449 
450  if (smax <= 3)
451  {
452  RKsolver = new SDIRK23Solver();
453  }
454  else
455  {
456  RKsolver = new SDIRK34Solver();
457  }
458 }
459 
461 {
462  MFEM_ASSERT( (i >= 0) && ( i < s ),
463  " AdamsMoultonSolver::GetStateVector \n" <<
464  " - Tried to get non-existent state "<<i);
465  return k[idx[i+1]];
466 }
467 
469 {
470  MFEM_ASSERT( (i >= 0) && ( i < s ),
471  " AdamsMoultonSolver::GetStateVector \n" <<
472  " - Tried to get non-existent state "<<i);
473  state = k[idx[i+1]];
474 }
475 
477 {
478  MFEM_ASSERT( (i >= 0) && ( i < smax ),
479  " AdamsMoultonSolver::SetStateVector \n" <<
480  " - Tried to set non-existent state "<<i);
481  k[idx[i+1]] = state;
482  s = std::max(i,s);
483 }
484 
486 {
487  ODESolver::Init(f_);
488  RKsolver->Init(f_);
489  int n = f->Width();
490  idx.SetSize(smax);
491  for (int i = 0; i < smax; i++)
492  {
493  idx[i] = (smax-i)%smax;
494  k[i].SetSize(n);
495  }
496  s = 0;
497 }
498 
499 void AdamsMoultonSolver::Step(Vector &x, double &t, double &dt)
500 {
501  if ((s == 0)&&(smax>1))
502  {
503  f->Mult(x,k[idx[1]]);
504  }
505  s++;
506  s = std::min(s, smax);
507 
508  if (s >= smax-1)
509  {
510  f->SetTime(t);
511  for (int i = 1; i < smax; i++)
512  {
513  x.Add(a[i]*dt, k[idx[i]]);
514  }
515  f->ImplicitSolve(a[0]*dt, x, k[idx[0]]);
516  x.Add(a[0]*dt, k[idx[0]]);
517  }
518  else
519  {
520  RKsolver->Step(x,t,dt);
521  f->Mult(x,k[idx[0]]);
522  }
523  t += dt;
524 
525  // Shift the index
526  for (int i = 0; i < smax; i++) { idx[i] = ++idx[i]%smax; }
527 }
528 
529 const double AM0Solver::a[] =
530 {1.0};
531 const double AM1Solver::a[] =
532 {0.5, 0.5};
533 const double AM2Solver::a[] =
534 {5.0/12.0, 2.0/3.0, -1.0/12.0};
535 const double AM3Solver::a[] =
536 {3.0/8.0, 19.0/24.0,-5.0/24.0, 1.0/24.0};
537 const double AM4Solver::a[] =
538 {251.0/720.0,646.0/720.0,-264.0/720.0, 106.0/720.0, -19.0/720.0};
539 
540 
542 {
543  ODESolver::Init(f_);
544  k.SetSize(f->Width(), mem_type);
545 }
546 
547 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
548 {
549  f->SetTime(t + dt);
550  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
551  x.Add(dt, k);
552  t += dt;
553 }
554 
555 
557 {
558  ODESolver::Init(f_);
559  k.SetSize(f->Width(), mem_type);
560 }
561 
562 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
563 {
564  f->SetTime(t + dt/2);
565  f->ImplicitSolve(dt/2, x, k);
566  x.Add(dt, k);
567  t += dt;
568 }
569 
570 
572 {
573  if (gamma_opt == 0)
574  {
575  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
576  }
577  else if (gamma_opt == 2)
578  {
579  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
580  }
581  else if (gamma_opt == 3)
582  {
583  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
584  }
585  else
586  {
587  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
588  }
589 }
590 
592 {
593  ODESolver::Init(f_);
594  k.SetSize(f->Width(), mem_type);
595  y.SetSize(f->Width(), mem_type);
596 }
597 
598 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
599 {
600  // with a = gamma:
601  // a | a
602  // 1-a | 1-2a a
603  // ------+-----------
604  // | 1/2 1/2
605  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
606  f->SetTime(t + gamma*dt);
607  f->ImplicitSolve(gamma*dt, x, k);
608  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
609  x.Add(dt/2, k);
610 
611  f->SetTime(t + (1.-gamma)*dt);
612  f->ImplicitSolve(gamma*dt, y, k);
613  x.Add(dt/2, k);
614  t += dt;
615 }
616 
617 
619 {
620  ODESolver::Init(f_);
621  k.SetSize(f->Width(), mem_type);
622  y.SetSize(f->Width(), mem_type);
623  z.SetSize(f->Width(), mem_type);
624 }
625 
626 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
627 {
628  // a | a
629  // 1/2 | 1/2-a a
630  // 1-a | 2a 1-4a a
631  // ------+--------------------
632  // | b 1-2b b
633  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
634  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
635  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
636 
637  f->SetTime(t + a*dt);
638  f->ImplicitSolve(a*dt, x, k);
639  add(x, (0.5-a)*dt, k, y);
640  add(x, (2.*a)*dt, k, z);
641  x.Add(b*dt, k);
642 
643  f->SetTime(t + dt/2);
644  f->ImplicitSolve(a*dt, y, k);
645  z.Add((1.-4.*a)*dt, k);
646  x.Add((1.-2.*b)*dt, k);
647 
648  f->SetTime(t + (1.-a)*dt);
649  f->ImplicitSolve(a*dt, z, k);
650  x.Add(b*dt, k);
651  t += dt;
652 }
653 
654 
656 {
657  ODESolver::Init(f_);
658  k.SetSize(f->Width(), mem_type);
659  y.SetSize(f->Width(), mem_type);
660 }
661 
662 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
663 {
664  // a | a
665  // c | c-a a
666  // 1 | b 1-a-b a
667  // -----+----------------
668  // | b 1-a-b a
669  const double a = 0.435866521508458999416019;
670  const double b = 1.20849664917601007033648;
671  const double c = 0.717933260754229499708010;
672 
673  f->SetTime(t + a*dt);
674  f->ImplicitSolve(a*dt, x, k);
675  add(x, (c-a)*dt, k, y);
676  x.Add(b*dt, k);
677 
678  f->SetTime(t + c*dt);
679  f->ImplicitSolve(a*dt, y, k);
680  x.Add((1.0-a-b)*dt, k);
681 
682  f->SetTime(t + dt);
683  f->ImplicitSolve(a*dt, x, k);
684  x.Add(a*dt, k);
685  t += dt;
686 }
687 
689 {
690  ODESolver::Init(f_);
691  k.SetSize(f->Width(), mem_type);
692  y.SetSize(f->Width(), mem_type);
693 }
694 
695 void TrapezoidalRuleSolver::Step(Vector &x, double &t, double &dt)
696 {
697  // 0 | 0 0
698  // 1 | 1/2 1/2
699  // ------+-----------
700  // | 1/2 1/2
701  f->SetTime(t);
702  f->Mult(x,k);
703  add(x, dt/2.0, k, y);
704  x.Add(dt/2.0, k);
705 
706  f->SetTime(t + dt);
707  f->ImplicitSolve(dt/2.0, y, k);
708  x.Add(dt/2.0, k);
709  t += dt;
710 }
711 
713 {
714  ODESolver::Init(f_);
715  k.SetSize(f->Width(), mem_type);
716  y.SetSize(f->Width(), mem_type);
717  z.SetSize(f->Width(), mem_type);
718 }
719 
720 void ESDIRK32Solver::Step(Vector &x, double &t, double &dt)
721 {
722  // 0 | 0 0 0
723  // 2a | a a 0
724  // 1 | 1-b-a b a
725  // ------+--------------------
726  // | 1-b-a b a
727  const double a = (2.0 - sqrt(2.0)) / 2.0;
728  const double b = (1.0 - 2.0*a) / (4.0*a);
729 
730  f->SetTime(t);
731  f->Mult(x,k);
732  add(x, a*dt, k, y);
733  add(x, (1.0-b-a)*dt, k, z);
734  x.Add((1.0-b-a)*dt, k);
735 
736  f->SetTime(t + (2.0*a)*dt);
737  f->ImplicitSolve(a*dt, y, k);
738  z.Add(b*dt, k);
739  x.Add(b*dt, k);
740 
741  f->SetTime(t + dt);
742  f->ImplicitSolve(a*dt, z, k);
743  x.Add(a*dt, k);
744  t += dt;
745 }
746 
748 {
749  ODESolver::Init(f_);
750  k.SetSize(f->Width(), mem_type);
751  y.SetSize(f->Width(), mem_type);
752  z.SetSize(f->Width(), mem_type);
753 }
754 
755 void ESDIRK33Solver::Step(Vector &x, double &t, double &dt)
756 {
757  // 0 | 0 0 0
758  // 2a | a a 0
759  // 1 | 1-b-a b a
760  // ------+----------------------------
761  // | 1-b_2-b_3 b_2 b_3
762  const double a = (3.0 + sqrt(3.0)) / 6.0;
763  const double b = (1.0 - 2.0*a) / (4.0*a);
764  const double b_2 = 1.0 / ( 12.0*a*(1.0 - 2.0*a) );
765  const double b_3 = (1.0 - 3.0*a) / ( 3.0*(1.0 - 2.0*a) );
766 
767  f->SetTime(t);
768  f->Mult(x,k);
769  add(x, a*dt, k, y);
770  add(x, (1.0-b-a)*dt, k, z);
771  x.Add((1.0-b_2-b_3)*dt, k);
772 
773  f->SetTime(t + (2.0*a)*dt);
774  f->ImplicitSolve(a*dt, y, k);
775  z.Add(b*dt, k);
776  x.Add(b_2*dt, k);
777 
778  f->SetTime(t + dt);
779  f->ImplicitSolve(a*dt, z, k);
780  x.Add(b_3*dt, k);
781  t += dt;
782 }
783 
785 {
786  ODESolver::Init(f_);
787  k.SetSize(f->Width(), mem_type);
788  y.SetSize(f->Width(), mem_type);
789  xdot.SetSize(f->Width(), mem_type);
790  xdot = 0.0;
791  nstate = 0;
792 }
793 
795 {
796  MFEM_ASSERT( (i == 0) && (nstate == 1),
797  "GeneralizedAlphaSolver::GetStateVector \n" <<
798  " - Tried to get non-existent state "<<i);
799  return xdot;
800 }
801 
803 {
804  MFEM_ASSERT( (i == 0) && (nstate == 1),
805  "GeneralizedAlphaSolver::GetStateVector \n" <<
806  " - Tried to get non-existent state "<<i);
807  state = xdot;
808 }
809 
811 {
812  MFEM_ASSERT( (i == 0),
813  "GeneralizedAlphaSolver::SetStateVector \n" <<
814  " - Tried to set non-existent state "<<i);
815  xdot = state;
816  nstate = 1;
817 }
818 
820 {
821  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
822  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
823 
824  // According to Jansen
825  alpha_m = 0.5*(3.0 - rho_inf)/(1.0 + rho_inf);
826  alpha_f = 1.0/(1.0 + rho_inf);
827  gamma = 0.5 + alpha_m - alpha_f;
828 }
829 
831 {
832  os << "Generalized alpha time integrator:" << std::endl;
833  os << "alpha_m = " << alpha_m << std::endl;
834  os << "alpha_f = " << alpha_f << std::endl;
835  os << "gamma = " << gamma << std::endl;
836 
837  if (gamma == 0.5 + alpha_m - alpha_f)
838  {
839  os<<"Second order"<<" and ";
840  }
841  else
842  {
843  os<<"First order"<<" and ";
844  }
845 
846  if ((alpha_m >= alpha_f)&&(alpha_f >= 0.5))
847  {
848  os<<"Stable"<<std::endl;
849  }
850  else
851  {
852  os<<"Unstable"<<std::endl;
853  }
854 }
855 
856 // This routine assumes xdot is initialized.
857 void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
858 {
859  if (nstate == 0)
860  {
861  f->Mult(x,xdot);
862  nstate = 1;
863  }
864 
865  // Set y = x + alpha_f*(1.0 - (gamma/alpha_m))*dt*xdot
866  add(x, alpha_f*(1.0 - (gamma/alpha_m))*dt, xdot, y);
867 
868  // Solve k = f(y + dt_eff*k)
869  double dt_eff = (gamma*alpha_f/alpha_m)*dt;
870  f->SetTime(t + alpha_f*dt);
871  f->ImplicitSolve(dt_eff, y, k);
872 
873  // Update x and xdot
874  x.Add((1.0 - (gamma/alpha_m))*dt, xdot);
875  x.Add( (gamma/alpha_m) *dt, k);
876 
877  xdot *= (1.0-(1.0/alpha_m));
878  xdot.Add((1.0/alpha_m),k);
879 
880  t += dt;
881 }
882 
883 
884 void
886 {
887  P_ = &P; F_ = &F;
888 
889  dp_.SetSize(F_->Height());
890  dq_.SetSize(P_->Height());
891 }
892 
893 void
894 SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
895 {
896  F_->SetTime(t);
897  F_->Mult(q,dp_);
898  p.Add(dt,dp_);
899 
900  P_->Mult(p,dq_);
901  q.Add(dt,dq_);
902 
903  t += dt;
904 }
905 
906 void
907 SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
908 {
909  P_->Mult(p,dq_);
910  q.Add(0.5*dt,dq_);
911 
912  F_->SetTime(t+0.5*dt);
913  F_->Mult(q,dp_);
914  p.Add(dt,dp_);
915 
916  P_->Mult(p,dq_);
917  q.Add(0.5*dt,dq_);
918 
919  t += dt;
920 }
921 
923  : order_(order)
924 {
925  a_.SetSize(order);
926  b_.SetSize(order);
927 
928  switch (order_)
929  {
930  case 1:
931  a_[0] = 1.0;
932  b_[0] = 1.0;
933  break;
934  case 2:
935  a_[0] = 0.5;
936  a_[1] = 0.5;
937  b_[0] = 0.0;
938  b_[1] = 1.0;
939  break;
940  case 3:
941  a_[0] = 2.0/3.0;
942  a_[1] = -2.0/3.0;
943  a_[2] = 1.0;
944  b_[0] = 7.0/24.0;
945  b_[1] = 0.75;
946  b_[2] = -1.0/24.0;
947  break;
948  case 4:
949  a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
950  a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
951  a_[2] = a_[1];
952  a_[3] = a_[0];
953  b_[0] = 0.0;
954  b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
955  b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
956  b_[3] = b_[1];
957  break;
958  default:
959  MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
960  };
961 }
962 
963 void
964 SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
965 {
966  for (int i=0; i<order_; i++)
967  {
968  if ( b_[i] != 0.0 )
969  {
970  F_->SetTime(t);
971  if ( F_->isExplicit() )
972  {
973  F_->Mult(q, dp_);
974  }
975  else
976  {
977  F_->ImplicitSolve(b_[i] * dt, q, dp_);
978  }
979  p.Add(b_[i] * dt, dp_);
980  }
981 
982  P_->Mult(p, dq_);
983  q.Add(a_[i] * dt, dq_);
984 
985  t += a_[i] * dt;
986  }
987 }
988 
990 {
991  this->f = &f_;
993 }
994 
996 {
998  d2xdt2.SetSize(f->Width());
999  d2xdt2 = 0.0;
1000  first = true;
1001 }
1002 
1003 void NewmarkSolver::PrintProperties(std::ostream &os)
1004 {
1005  os << "Newmark time integrator:" << std::endl;
1006  os << "beta = " << beta << std::endl;
1007  os << "gamma = " << gamma << std::endl;
1008 
1009  if (gamma == 0.5)
1010  {
1011  os<<"Second order"<<" and ";
1012  }
1013  else
1014  {
1015  os<<"First order"<<" and ";
1016  }
1017 
1018  if ((gamma >= 0.5) && (beta >= (gamma + 0.5)*(gamma + 0.5)/4))
1019  {
1020  os<<"A-Stable"<<std::endl;
1021  }
1022  else if ((gamma >= 0.5) && (beta >= 0.5*gamma))
1023  {
1024  os<<"Conditionally stable"<<std::endl;
1025  }
1026  else
1027  {
1028  os<<"Unstable"<<std::endl;
1029  }
1030 }
1031 
1032 void NewmarkSolver::Step(Vector &x, Vector &dxdt, double &t, double &dt)
1033 {
1034  double fac0 = 0.5 - beta;
1035  double fac2 = 1.0 - gamma;
1036  double fac3 = beta;
1037  double fac4 = gamma;
1038 
1039  // In the first pass compute d2xdt2 directly from operator.
1040  if (first)
1041  {
1042  f->Mult(x, dxdt, d2xdt2);
1043  first = false;
1044  }
1045  f->SetTime(t + dt);
1046 
1047  x.Add(dt, dxdt);
1048  x.Add(fac0*dt*dt, d2xdt2);
1049  dxdt.Add(fac2*dt, d2xdt2);
1050 
1051  f->SetTime(t + dt);
1052  f->ImplicitSolve(fac3*dt*dt, fac4*dt, x, dxdt, d2xdt2);
1053 
1054  x .Add(fac3*dt*dt, d2xdt2);
1055  dxdt.Add(fac4*dt, d2xdt2);
1056  t += dt;
1057 }
1058 
1060 {
1062  xa.SetSize(f->Width());
1063  va.SetSize(f->Width());
1064  aa.SetSize(f->Width());
1065  d2xdt2.SetSize(f->Width());
1066  d2xdt2 = 0.0;
1067  nstate = 0;
1068 }
1069 
1071 {
1072  MFEM_ASSERT( (i == 0) && (nstate == 1),
1073  "GeneralizedAlpha2Solver::GetStateVector \n" <<
1074  " - Tried to get non-existent state "<<i);
1075  return d2xdt2;
1076 }
1077 
1078 
1080 {
1081  MFEM_ASSERT( (i == 0) && (nstate == 1),
1082  "GeneralizedAlpha2Solver::GetStateVector \n" <<
1083  " - Tried to get non-existent state "<<i);
1084  state = d2xdt2;
1085 }
1086 
1088 {
1089  MFEM_ASSERT( (i == 0),
1090  "GeneralizedAlpha2Solver::SetStateVector \n" <<
1091  " - Tried to set non-existent state "<<i);
1092  d2xdt2 = state;
1093  nstate = 1;
1094 }
1095 
1097 {
1098  os << "Generalized alpha time integrator:" << std::endl;
1099  os << "alpha_m = " << alpha_m << std::endl;
1100  os << "alpha_f = " << alpha_f << std::endl;
1101  os << "beta = " << beta << std::endl;
1102  os << "gamma = " << gamma << std::endl;
1103 
1104  if (gamma == 0.5 + alpha_m - alpha_f)
1105  {
1106  os<<"Second order"<<" and ";
1107  }
1108  else
1109  {
1110  os<<"First order"<<" and ";
1111  }
1112 
1113  if ((alpha_m >= alpha_f)&&
1114  (alpha_f >= 0.5) &&
1115  (beta >= 0.25 + 0.5*(alpha_m - alpha_f)))
1116  {
1117  os<<"Stable"<<std::endl;
1118  }
1119  else
1120  {
1121  os<<"Unstable"<<std::endl;
1122  }
1123 }
1124 
1126  double &t, double &dt)
1127 {
1128  double fac0 = (0.5 - (beta/alpha_m));
1129  double fac1 = alpha_f;
1130  double fac2 = alpha_f*(1.0 - (gamma/alpha_m));
1131  double fac3 = beta*alpha_f/alpha_m;
1132  double fac4 = gamma*alpha_f/alpha_m;
1133  double fac5 = alpha_m;
1134 
1135  // In the first pass compute d2xdt2 directly from operator.
1136  if (nstate == 0)
1137  {
1138  f->Mult(x, dxdt, d2xdt2);
1139  nstate = 1;
1140  }
1141 
1142  // Predict alpha levels
1143  add(dxdt, fac0*dt, d2xdt2, va);
1144  add(x, fac1*dt, va, xa);
1145  add(dxdt, fac2*dt, d2xdt2, va);
1146 
1147  // Solve alpha levels
1148  f->SetTime(t + dt);
1149  f->ImplicitSolve(fac3*dt*dt, fac4*dt, xa, va, aa);
1150 
1151  // Correct alpha levels
1152  xa.Add(fac3*dt*dt, aa);
1153  va.Add(fac4*dt, aa);
1154 
1155  // Extrapolate
1156  x *= 1.0 - 1.0/fac1;
1157  x.Add (1.0/fac1, xa);
1158 
1159  dxdt *= 1.0 - 1.0/fac1;
1160  dxdt.Add (1.0/fac1, va);
1161 
1162  d2xdt2 *= 1.0 - 1.0/fac5;
1163  d2xdt2.Add (1.0/fac5, aa);
1164 
1165  t += dt;
1166 }
1167 
1168 }
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:598
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1059
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:720
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:408
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:830
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:885
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:562
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:688
void SetRhoInf(double rho_inf)
Definition: ode.cpp:819
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:541
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
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
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:27
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:695
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:1125
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:907
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:376
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:655
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:747
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:964
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:1032
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:784
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:1003
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:1087
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:618
Vector dq_
Definition: ode.hpp:571
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:989
void SetStateVector(int i, Vector &state) 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: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:857
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:755
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:794
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:626
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:712
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
TimeDependentOperator * F_
Definition: ode.hpp:567
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:556
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
Base abstract class for second order time dependent operators.
Definition: operator.hpp:634
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1096
Vector dp_
Definition: ode.hpp:570
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:591
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:995
AdamsBashforthSolver(int s_, const double *a_)
Definition: ode.cpp:347
Operator * P_
Definition: ode.hpp:568
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:386
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:253
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:485
const Vector & GetStateVector(int i) override
Definition: ode.cpp:1070
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:55
SIAVSolver(int order)
Definition: ode.cpp:922
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:499
const Vector & GetStateVector(int i) override
Definition: ode.cpp:460
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:662
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:547
Vector data type.
Definition: vector.hpp:60
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:395
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:26
RefCoord s[3]
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:810
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:894
Abstract operator.
Definition: operator.hpp:24
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:571
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:443
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:611