MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ode.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
15namespace mfem
16{
17
19{
20 this->f = &f_;
22}
23
25{
27 dxdt.SetSize(f->Width(), mem_type);
28}
29
31{
32 f->SetTime(t);
33 f->Mult(x, dxdt);
34 x.Add(dt, dxdt);
35 t += dt;
36}
37
38
40{
42 int n = f->Width();
43 dxdt.SetSize(n, mem_type);
44 x1.SetSize(n, mem_type);
45}
46
48{
49 // 0 |
50 // a | a
51 // ---+--------
52 // | 1-b b b = 1/(2a)
53
54 const real_t 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{
71 int n = f->Width();
72 y.SetSize(n, mem_type);
73 k.SetSize(n, mem_type);
74}
75
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
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
139 const real_t *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
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
194const real_t 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};
225const real_t RK6Solver::b[] =
226{
227 .3438957868357036009278820124728322386520e-1,
228 0.,
229 0.,
230 .2582624555633503404659558098586120858767,
231 .4209371189673537150642551514069801967032,
232 4.405396469669310170148836816197095664891,
233 -176.4831190242986576151740942499002125029,
234 172.3641334014150730294022582711902413315
235};
236const real_t RK6Solver::c[] =
237{
238 .6e-1,
239 .9593333333333333333333333333333333333333e-1,
240 .1439,
241 .4973,
242 .9725,
243 .9995,
244 1.,
245};
246
247const real_t 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};
316const real_t 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};
331const real_t 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
410{
411 if ( (dt_ > 0.0) && (fabs(dt-dt_) >10*std::numeric_limits<real_t>::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
447const real_t AB1Solver::a[] =
448{1.0};
449const real_t AB2Solver::a[] =
450{1.5,-0.5};
451const real_t AB3Solver::a[] =
452{23.0/12.0,-4.0/3.0, 5.0/12.0};
453const real_t AB4Solver::a[] =
454{55.0/24.0,-59.0/24.0, 37.0/24.0,-9.0/24.0};
455const real_t AB5Solver::a[] =
456{1901.0/720.0,-2774.0/720.0, 2616.0/720.0,-1274.0/720.0, 251.0/720.0};
457
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
516{
517 if ( (dt_ > 0.0) && (fabs(dt-dt_) >10*std::numeric_limits<real_t>::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
560const real_t AM0Solver::a[] =
561{1.0};
562const real_t AM1Solver::a[] =
563{0.5, 0.5};
564const real_t AM2Solver::a[] =
565{5.0/12.0, 2.0/3.0, -1.0/12.0};
566const real_t AM3Solver::a[] =
567{3.0/8.0, 19.0/24.0,-5.0/24.0, 1.0/24.0};
568const real_t 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
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
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
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
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 real_t a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
666 const real_t 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
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 real_t a = 0.435866521508458999416019;
701 const real_t b = 1.20849664917601007033648;
702 const real_t 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
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
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 real_t a = (2.0 - sqrt(2.0)) / 2.0;
759 const real_t 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
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 real_t a = (3.0 + sqrt(3.0)) / 6.0;
794 const real_t b = (1.0 - 2.0*a) / (4.0*a);
795 const real_t b_2 = 1.0 / ( 12.0*a*(1.0 - 2.0*a) );
796 const real_t 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);
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.
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 real_t 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
915void
917{
918 P_ = &P; F_ = &F;
919
920 dp_.SetSize(F_->Height());
921 dq_.SetSize(P_->Height());
922}
923
924void
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
937void
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
994void
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
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
1064{
1065 real_t fac0 = 0.5 - beta;
1066 real_t fac2 = 1.0 - gamma;
1067 real_t fac3 = beta;
1068 real_t 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 real_t &t, real_t &dt)
1158{
1159 real_t fac0 = (0.5 - (beta/alpha_m));
1160 real_t fac1 = alpha_f;
1161 real_t fac2 = alpha_f*(1.0 - (gamma/alpha_m));
1162 real_t fac3 = beta*alpha_f/alpha_m;
1163 real_t fac4 = gamma*alpha_f/alpha_m;
1164 real_t 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 Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:396
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:387
AdamsBashforthSolver(int s_, const real_t *a_)
Definition: ode.cpp:347
const Vector & GetStateVector(int i) override
Definition: ode.cpp:377
void Step(Vector &x, real_t &t, real_t &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
const Vector & GetStateVector(int i) override
Definition: ode.cpp:476
void Step(Vector &x, real_t &t, real_t &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
AdamsMoultonSolver(int s_, const real_t *a_)
Definition: ode.cpp:458
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:501
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:492
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:697
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:572
void Step(Vector &x, real_t &t, real_t &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
virtual void Step(Vector &x, real_t &t, real_t &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
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:743
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:778
virtual void Step(Vector &x, real_t &t, real_t &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, real_t &t, real_t &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
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
ExplicitRKSolver(int s_, const real_t *a_, const real_t *b_, const real_t *c_)
Definition: ode.cpp:138
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148
void Step(Vector &x, real_t &t, real_t &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
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1090
void Step(Vector &x, Vector &dxdt, real_t &t, real_t &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 SetStateVector(int i, Vector &state) override
Definition: ode.cpp:1118
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1127
const Vector & GetStateVector(int i) override
Definition: ode.cpp:1101
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:815
void SetRhoInf(real_t rho_inf)
Definition: ode.cpp:850
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:841
void Step(Vector &x, real_t &t, real_t &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
const Vector & GetStateVector(int i) override
Definition: ode.cpp:825
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:861
void Step(Vector &x, real_t &t, real_t &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
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:587
void Step(Vector &x, Vector &dxdt, real_t &t, real_t &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 PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1034
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1026
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:27
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
MemoryType mem_type
Definition: ode.hpp:28
Abstract operator.
Definition: operator.hpp:25
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
void Step(Vector &x, real_t &t, real_t &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
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:151
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
void Step(Vector &x, real_t &t, real_t &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
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:164
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
void Step(Vector &x, real_t &t, real_t &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
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:602
void Step(Vector &x, real_t &t, real_t &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(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:622
void Step(Vector &x, real_t &t, real_t &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
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:686
void Step(Vector &x, real_t &t, real_t &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
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:649
void Step(Vector &q, Vector &p, real_t &t, real_t &dt) override
Definition: ode.cpp:925
void Step(Vector &q, Vector &p, real_t &t, real_t &dt) override
Definition: ode.cpp:938
TimeDependentOperator * F_
Definition: ode.hpp:587
Vector dq_
Definition: ode.hpp:591
Vector dp_
Definition: ode.hpp:590
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:916
Operator * P_
Definition: ode.hpp:588
void Step(Vector &q, Vector &p, real_t &t, real_t &dt) override
Definition: ode.cpp:995
SIAVSolver(int order)
Definition: ode.cpp:953
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:631
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1020
Base abstract class for second order time dependent operators.
Definition: operator.hpp:635
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
virtual void ImplicitSolve(const real_t fac0, const real_t 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
Base abstract class for first order time dependent operators.
Definition: operator.hpp:317
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:363
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
virtual void SetTime(const real_t t_)
Set the current time.
Definition: operator.hpp:360
virtual void ImplicitSolve(const real_t 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
virtual void Step(Vector &x, real_t &t, real_t &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
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:719
Vector data type.
Definition: vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:538
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:247
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
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 add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:316
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:60
float real_t
Definition: config.hpp:43
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53
RefCoord t[3]
RefCoord s[3]