MFEM  v4.6.0
Finite element discretization library
ode.hpp
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 #ifndef MFEM_ODE
13 #define MFEM_ODE
14 
15 #include "../config/config.hpp"
16 #include "operator.hpp"
17 #include "../general/communication.hpp"
18 
19 namespace mfem
20 {
21 
22 /// Abstract class for solving systems of ODEs: dx/dt = f(x,t)
23 class ODESolver
24 {
25 protected:
26  /// Pointer to the associated TimeDependentOperator.
27  TimeDependentOperator *f; // f(.,t) : R^n --> R^n
29 
30 public:
32 
33  /// Associate a TimeDependentOperator with the ODE solver.
34  /** This method has to be called:
35  - Before the first call to Step().
36  - When the dimensions of the associated TimeDependentOperator change.
37  - When a time stepping sequence has to be restarted.
38  - To change the associated TimeDependentOperator. */
39  virtual void Init(TimeDependentOperator &f_);
40 
41  /** @brief Perform a time step from time @a t [in] to time @a t [out] based
42  on the requested step size @a dt [in]. */
43  /** @param[in,out] x Approximate solution.
44  @param[in,out] t Time associated with the approximate solution @a x.
45  @param[in,out] dt Time step size.
46 
47  The following rules describe the common behavior of the method:
48  - The input @a x [in] is the approximate solution for the input time
49  @a t [in].
50  - The input @a dt [in] is the desired time step size, defining the desired
51  target time: t [target] = @a t [in] + @a dt [in].
52  - The output @a x [out] is the approximate solution for the output time
53  @a t [out].
54  - The output @a dt [out] is the last time step taken by the method which
55  may be smaller or larger than the input @a dt [in] value, e.g. because
56  of time step control.
57  - The method may perform more than one time step internally; in this case
58  @a dt [out] is the last internal time step size.
59  - The output value of @a t [out] may be smaller or larger than
60  t [target], however, it is not smaller than @a t [in] + @a dt [out], if
61  at least one internal time step was performed.
62  - The value @a x [out] may be obtained by interpolation using internally
63  stored data.
64  - In some cases, the contents of @a x [in] may not be used, e.g. when
65  @a x [out] from a previous Step() call was obtained by interpolation.
66  - In consecutive calls to this method, the output @a t [out] of one
67  Step() call has to be the same as the input @a t [in] to the next
68  Step() call.
69  - If the previous rule has to be broken, e.g. to restart a time stepping
70  sequence, then the ODE solver must be re-initialized by calling Init()
71  between the two Step() calls. */
72  virtual void Step(Vector &x, double &t, double &dt) = 0;
73 
74  /// Perform time integration from time @a t [in] to time @a tf [in].
75  /** @param[in,out] x Approximate solution.
76  @param[in,out] t Time associated with the approximate solution @a x.
77  @param[in,out] dt Time step size.
78  @param[in] tf Requested final time.
79 
80  The default implementation makes consecutive calls to Step() until
81  reaching @a tf.
82  The following rules describe the common behavior of the method:
83  - The input @a x [in] is the approximate solution for the input time
84  @a t [in].
85  - The input @a dt [in] is the initial time step size.
86  - The output @a dt [out] is the last time step taken by the method which
87  may be smaller or larger than the input @a dt [in] value, e.g. because
88  of time step control.
89  - The output value of @a t [out] is not smaller than @a tf [in]. */
90  virtual void Run(Vector &x, double &t, double &dt, double tf)
91  {
92  while (t < tf) { Step(x, t, dt); }
93  }
94 
95  /// Function for getting and setting the state vectors
96  virtual int GetMaxStateSize() { return 0; }
97  virtual int GetStateSize() { return 0; }
98  virtual const Vector &GetStateVector(int i)
99  {
100  mfem_error("ODESolver has no state vectors");
101  Vector *s = NULL; return *s; // Make some compiler happy
102  }
103  virtual void GetStateVector(int i, Vector &state)
104  {
105  mfem_error("ODESolver has no state vectors");
106  }
107  virtual void SetStateVector(int i, Vector &state)
108  {
109  mfem_error("ODESolver has no state vectors");
110  }
111 
112  virtual ~ODESolver() { }
113 };
114 
115 
116 /// The classical forward Euler method
118 {
119 private:
120  Vector dxdt;
121 
122 public:
123  void Init(TimeDependentOperator &f_) override;
124 
125  void Step(Vector &x, double &t, double &dt) override;
126 };
127 
128 
129 /** A family of explicit second-order RK2 methods. Some choices for the
130  parameter 'a' are:
131  a = 1/2 - the midpoint method
132  a = 1 - Heun's method
133  a = 2/3 - default, has minimal truncation error. */
134 class RK2Solver : public ODESolver
135 {
136 private:
137  double a;
138  Vector dxdt, x1;
139 
140 public:
141  RK2Solver(const double a_ = 2./3.) : a(a_) { }
142 
143  void Init(TimeDependentOperator &f_) override;
144 
145  void Step(Vector &x, double &t, double &dt) override;
146 };
147 
148 
149 /// Third-order, strong stability preserving (SSP) Runge-Kutta method
150 class RK3SSPSolver : public ODESolver
151 {
152 private:
153  Vector y, k;
154 
155 public:
156  void Init(TimeDependentOperator &f_) override;
157 
158  void Step(Vector &x, double &t, double &dt) override;
159 };
160 
161 
162 /// The classical explicit forth-order Runge-Kutta method, RK4
163 class RK4Solver : public ODESolver
164 {
165 private:
166  Vector y, k, z;
167 
168 public:
169  void Init(TimeDependentOperator &f_) override;
170 
171  void Step(Vector &x, double &t, double &dt) override;
172 };
173 
174 
175 /** An explicit Runge-Kutta method corresponding to a general Butcher tableau
176  +--------+----------------------+
177  | c[0] | a[0] |
178  | c[1] | a[1] a[2] |
179  | ... | ... |
180  | c[s-2] | ... a[s(s-1)/2-1] |
181  +--------+----------------------+
182  | | b[0] b[1] ... b[s-1] |
183  +--------+----------------------+ */
185 {
186 private:
187  int s;
188  const double *a, *b, *c;
189  Vector y, *k;
190 
191 public:
192  ExplicitRKSolver(int s_, const double *a_, const double *b_,
193  const double *c_);
194 
195  void Init(TimeDependentOperator &f_) override;
196 
197  void Step(Vector &x, double &t, double &dt) override;
198 
199  virtual ~ExplicitRKSolver();
200 };
201 
202 
203 /** An 8-stage, 6th order RK method. From Verner's "efficient" 9-stage 6(5)
204  pair. */
206 {
207 private:
208  static MFEM_EXPORT const double a[28], b[8], c[7];
209 
210 public:
211  RK6Solver() : ExplicitRKSolver(8, a, b, c) { }
212 };
213 
214 
215 /** A 12-stage, 8th order RK method. From Verner's "efficient" 13-stage 8(7)
216  pair. */
218 {
219 private:
220  static const double a[66], b[12], c[11];
221 
222 public:
223  RK8Solver() : ExplicitRKSolver(12, a, b, c) { }
224 };
225 
226 
227 /** An explicit Adams-Bashforth method. */
229 {
230 private:
231  int s, smax;
232  const double *a;
233  Vector *k;
234  Array<int> idx;
235  ODESolver *RKsolver;
236  double dt_;
237 
238  inline bool print()
239  {
240 #ifdef MFEM_USE_MPI
241  return Mpi::IsInitialized() ? Mpi::Root() : true;
242 #else
243  return true;
244 #endif
245  }
246 
247 public:
248  AdamsBashforthSolver(int s_, const double *a_);
249 
250  void Init(TimeDependentOperator &f_) override;
251 
252  void Step(Vector &x, double &t, double &dt) override;
253 
254  int GetMaxStateSize() override { return smax; };
255  int GetStateSize() override { return s; };
256  const Vector &GetStateVector(int i) override;
257  void GetStateVector(int i, Vector &state) override;
258  void SetStateVector(int i, Vector &state) override;
259 
261  {
262  if (RKsolver) { delete RKsolver; }
263  delete [] k;
264  }
265 };
266 
267 /** A 1-stage, 1st order AB method. */
269 {
270 private:
271  static MFEM_EXPORT const double a[1];
272 
273 public:
275 };
276 
277 /** A 2-stage, 2nd order AB method. */
279 {
280 private:
281  static MFEM_EXPORT const double a[2];
282 
283 public:
285 };
286 
287 /** A 3-stage, 3rd order AB method. */
289 {
290 private:
291  static MFEM_EXPORT const double a[3];
292 
293 public:
295 };
296 
297 /** A 4-stage, 4th order AB method. */
299 {
300 private:
301  static MFEM_EXPORT const double a[4];
302 
303 public:
305 };
306 
307 /** A 5-stage, 5th order AB method. */
309 {
310 private:
311  static MFEM_EXPORT const double a[5];
312 
313 public:
315 };
316 
317 
318 /** An implicit Adams-Moulton method. */
320 {
321 private:
322  int s, smax;
323  const double *a;
324  Vector *k;
325  Array<int> idx;
326  ODESolver *RKsolver;
327  double dt_;
328 
329  inline bool print()
330  {
331 #ifdef MFEM_USE_MPI
332  return Mpi::IsInitialized() ? Mpi::Root() : true;
333 #else
334  return true;
335 #endif
336  }
337 
338 public:
339  AdamsMoultonSolver(int s_, const double *a_);
340 
341  void Init(TimeDependentOperator &f_) override;
342 
343  void Step(Vector &x, double &t, double &dt) override;
344 
345  int GetMaxStateSize() override { return smax-1; };
346  int GetStateSize() override { return s-1; };
347  const Vector &GetStateVector(int i) override;
348  void GetStateVector(int i, Vector &state) override;
349  void SetStateVector(int i, Vector &state) override;
350 
352  {
353  if (RKsolver) { delete RKsolver; }
354  delete [] k;
355  };
356 };
357 
358 /** A 0-stage, 1st order AM method. */
360 {
361 private:
362  static MFEM_EXPORT const double a[1];
363 
364 public:
366 };
367 
368 
369 /** A 1-stage, 2nd order AM method. */
371 {
372 private:
373  static MFEM_EXPORT const double a[2];
374 
375 public:
377 };
378 
379 /** A 2-stage, 3rd order AM method. */
381 {
382 private:
383  static MFEM_EXPORT const double a[3];
384 
385 public:
387 };
388 
389 /** A 3-stage, 4th order AM method. */
391 {
392 private:
393  static MFEM_EXPORT const double a[4];
394 
395 public:
397 };
398 
399 /** A 4-stage, 5th order AM method. */
401 {
402 private:
403  static MFEM_EXPORT const double a[5];
404 
405 public:
407 };
408 
409 
410 /// Backward Euler ODE solver. L-stable.
412 {
413 protected:
415 
416 public:
417  void Init(TimeDependentOperator &f_) override;
418 
419  void Step(Vector &x, double &t, double &dt) override;
420 };
421 
422 
423 /// Implicit midpoint method. A-stable, not L-stable.
425 {
426 protected:
428 
429 public:
430  void Init(TimeDependentOperator &f_) override;
431 
432  void Step(Vector &x, double &t, double &dt) override;
433 };
434 
435 
436 /** Two stage, singly diagonal implicit Runge-Kutta (SDIRK) methods;
437  the choices for gamma_opt are:
438  0 - 3rd order method, not A-stable
439  1 - 3rd order method, A-stable, not L-stable (default)
440  2 - 2nd order method, L-stable
441  3 - 2nd order method, L-stable (has solves outside [t,t+dt]). */
442 class SDIRK23Solver : public ODESolver
443 {
444 protected:
445  double gamma;
447 
448 public:
449  SDIRK23Solver(int gamma_opt = 1);
450 
451  void Init(TimeDependentOperator &f_) override;
452 
453  void Step(Vector &x, double &t, double &dt) override;
454 };
455 
456 
457 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
458  order 4. A-stable, not L-stable. */
459 class SDIRK34Solver : public ODESolver
460 {
461 protected:
462  Vector k, y, z;
463 
464 public:
465  void Init(TimeDependentOperator &f_) override;
466 
467  void Step(Vector &x, double &t, double &dt) override;
468 };
469 
470 
471 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
472  order 3. L-stable. */
473 class SDIRK33Solver : public ODESolver
474 {
475 protected:
477 
478 public:
479  void Init(TimeDependentOperator &f_) override;
480 
481  void Step(Vector &x, double &t, double &dt) override;
482 };
483 
484 
485 /** Two stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
486  of order 2. A-stable. */
488 {
489 protected:
491 
492 public:
493  virtual void Init(TimeDependentOperator &f_);
494 
495  virtual void Step(Vector &x, double &t, double &dt);
496 };
497 
498 
499 /** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
500  of order 2. L-stable. */
501 class ESDIRK32Solver : public ODESolver
502 {
503 protected:
504  Vector k, y, z;
505 
506 public:
507  virtual void Init(TimeDependentOperator &f_);
508 
509  virtual void Step(Vector &x, double &t, double &dt);
510 };
511 
512 
513 /** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
514  of order 3. A-stable. */
515 class ESDIRK33Solver : public ODESolver
516 {
517 protected:
518  Vector k, y, z;
519 
520 public:
521  virtual void Init(TimeDependentOperator &f_);
522 
523  virtual void Step(Vector &x, double &t, double &dt);
524 };
525 
526 
527 /// Generalized-alpha ODE solver from "A generalized-α method for integrating
528 /// the filtered Navier-Stokes equations with a stabilized finite element
529 /// method" by K.E. Jansen, C.H. Whiting and G.M. Hulbert.
531 {
532 protected:
533  mutable Vector xdot,k,y;
535  int nstate;
536 
537  void SetRhoInf(double rho_inf);
538  void PrintProperties(std::ostream &out = mfem::out);
539 public:
540 
541  GeneralizedAlphaSolver(double rho = 1.0) { SetRhoInf(rho); };
542 
543  void Init(TimeDependentOperator &f_) override;
544 
545  void Step(Vector &x, double &t, double &dt) override;
546 
547  int GetMaxStateSize() override { return 1; };
548  int GetStateSize() override { return nstate; };
549  const Vector &GetStateVector(int i) override;
550  void GetStateVector(int i, Vector &state) override;
551  void SetStateVector(int i, Vector &state) override;
552 };
553 
554 
555 /// The SIASolver class is based on the Symplectic Integration Algorithm
556 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
557 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,
558 /// Vol. 92, pages 230-256 (1991).
559 
560 /** The Symplectic Integration Algorithm (SIA) is designed for systems of first
561  order ODEs derived from a Hamiltonian.
562  H(q,p,t) = T(p) + V(q,t)
563  Which leads to the equations:
564  dq/dt = dT/dp
565  dp/dt = -dV/dq
566  In the integrator the operators P and F are defined to be:
567  P = dT/dp
568  F = -dV/dq
569  */
571 {
572 public:
573  SIASolver() : F_(NULL), P_(NULL) {}
574 
575  virtual void Init(Operator &P, TimeDependentOperator & F);
576 
577  virtual void Step(Vector &q, Vector &p, double &t, double &dt) = 0;
578 
579  virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
580  {
581  while (t < tf) { Step(q, p, t, dt); }
582  }
583 
584  virtual ~SIASolver() {}
585 
586 protected:
587  TimeDependentOperator * F_; // p_{i+1} = p_{i} + dt F(q_{i})
588  Operator * P_; // q_{i+1} = q_{i} + dt P(p_{i+1})
589 
590  mutable Vector dp_;
591  mutable Vector dq_;
592 };
593 
594 /// First Order Symplectic Integration Algorithm
595 class SIA1Solver : public SIASolver
596 {
597 public:
599  void Step(Vector &q, Vector &p, double &t, double &dt) override;
600 };
601 
602 /// Second Order Symplectic Integration Algorithm
603 class SIA2Solver : public SIASolver
604 {
605 public:
607  void Step(Vector &q, Vector &p, double &t, double &dt) override;
608 };
609 
610 /// Variable order Symplectic Integration Algorithm (orders 1-4)
611 class SIAVSolver : public SIASolver
612 {
613 public:
614  SIAVSolver(int order);
615  void Step(Vector &q, Vector &p, double &t, double &dt) override;
616 
617 private:
618  int order_;
619 
620  Array<double> a_;
621  Array<double> b_;
622 };
623 
624 
625 
626 /// Abstract class for solving systems of ODEs: d2x/dt2 = f(x,dx/dt,t)
628 {
629 protected:
630  /// Pointer to the associated TimeDependentOperator.
631  SecondOrderTimeDependentOperator *f; // f(.,.,t) : R^n x R^n --> R^n
633 
634 public:
636 
637  /// Associate a TimeDependentOperator with the ODE solver.
638  /** This method has to be called:
639  - Before the first call to Step().
640  - When the dimensions of the associated TimeDependentOperator change.
641  - When a time stepping sequence has to be restarted.
642  - To change the associated TimeDependentOperator. */
643  virtual void Init(SecondOrderTimeDependentOperator &f);
644 
645  /** @brief Perform a time step from time @a t [in] to time @a t [out] based
646  on the requested step size @a dt [in]. */
647  /** @param[in,out] x Approximate solution.
648  @param[in,out] dxdt Approximate rate.
649  @param[in,out] t Time associated with the
650  approximate solution @a x and rate @ dxdt
651  @param[in,out] dt Time step size.
652 
653  The following rules describe the common behavior of the method:
654  - The input @a x [in] is the approximate solution for the input time
655  @a t [in].
656  - The input @a dxdt [in] is the approximate rate for the input time
657  @a t [in].
658  - The input @a dt [in] is the desired time step size, defining the desired
659  target time: t [target] = @a t [in] + @a dt [in].
660  - The output @a x [out] is the approximate solution for the output time
661  @a t [out].
662  - The output @a dxdt [out] is the approximate rate for the output time
663  @a t [out].
664  - The output @a dt [out] is the last time step taken by the method which
665  may be smaller or larger than the input @a dt [in] value, e.g. because
666  of time step control.
667  - The method may perform more than one time step internally; in this case
668  @a dt [out] is the last internal time step size.
669  - The output value of @a t [out] may be smaller or larger than
670  t [target], however, it is not smaller than @a t [in] + @a dt [out], if
671  at least one internal time step was performed.
672  - The value @a x [out] may be obtained by interpolation using internally
673  stored data.
674  - In some cases, the contents of @a x [in] may not be used, e.g. when
675  @a x [out] from a previous Step() call was obtained by interpolation.
676  - In consecutive calls to this method, the output @a t [out] of one
677  Step() call has to be the same as the input @a t [in] to the next
678  Step() call.
679  - If the previous rule has to be broken, e.g. to restart a time stepping
680  sequence, then the ODE solver must be re-initialized by calling Init()
681  between the two Step() calls. */
682  virtual void Step(Vector &x, Vector &dxdt, double &t, double &dt) = 0;
683 
684  /// Perform time integration from time @a t [in] to time @a tf [in].
685  /** @param[in,out] x Approximate solution.
686  @param[in,out] dxdt Approximate rate.
687  @param[in,out] t Time associated with the approximate solution @a x.
688  @param[in,out] dt Time step size.
689  @param[in] tf Requested final time.
690 
691  The default implementation makes consecutive calls to Step() until
692  reaching @a tf.
693  The following rules describe the common behavior of the method:
694  - The input @a x [in] is the approximate solution for the input time
695  @a t [in].
696  - The input @a dxdt [in] is the approximate rate for the input time
697  @a t [in].
698  - The input @a dt [in] is the initial time step size.
699  - The output @a dt [out] is the last time step taken by the method which
700  may be smaller or larger than the input @a dt [in] value, e.g. because
701  of time step control.
702  - The output value of @a t [out] is not smaller than @a tf [in]. */
703  virtual void Run(Vector &x, Vector &dxdt, double &t, double &dt, double tf)
704  {
705  while (t < tf) { Step(x, dxdt, t, dt); }
706  }
707 
708  /// Function for getting and setting the state vectors
709  virtual int GetMaxStateSize() { return 0; };
710  virtual int GetStateSize() { return 0; }
711  virtual const Vector &GetStateVector(int i)
712  {
713  mfem_error("ODESolver has no state vectors");
714  Vector *s = NULL; return *s; // Make some compiler happy
715  }
716  virtual void GetStateVector(int i, Vector &state)
717  {
718  mfem_error("ODESolver has no state vectors");
719  }
720  virtual void SetStateVector(int i, Vector &state)
721  {
722  mfem_error("ODESolver has no state vectors");
723  }
724 
725  virtual ~SecondOrderODESolver() { }
726 };
727 
728 /// The classical newmark method.
729 /// Newmark, N. M. (1959) A method of computation for structural dynamics.
730 /// Journal of Engineering Mechanics, ASCE, 85 (EM3) 67-94.
732 {
733 private:
734  Vector d2xdt2;
735 
736  double beta, gamma;
737  bool first;
738 
739 public:
740  NewmarkSolver(double beta_ = 0.25, double gamma_ = 0.5) { beta = beta_; gamma = gamma_; };
741 
742  void PrintProperties(std::ostream &out = mfem::out);
743 
744  void Init(SecondOrderTimeDependentOperator &f_) override;
745 
746  void Step(Vector &x, Vector &dxdt, double &t, double &dt) override;
747 };
748 
750 {
751 public:
753 };
754 
756 {
757 public:
759 };
760 
762 {
763 public:
764  FoxGoodwinSolver() : NewmarkSolver(1.0/12.0, 0.5) { };
765 };
766 
767 /// Generalized-alpha ODE solver
768 /// A Time Integration Algorithm for Structural Dynamics With Improved
769 /// Numerical Dissipation: The Generalized-α Method
770 /// J.Chung and G.M. Hulbert, J. Appl. Mech 60(2), 371-375, 1993
771 /// https://doi.org/10.1115/1.2900803
772 /// rho_inf in [0,1]
774 {
775 protected:
778  int nstate;
779 
780 public:
781  GeneralizedAlpha2Solver(double rho_inf = 1.0)
782  {
783  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
784  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
785 
786  alpha_m = (2.0 - rho_inf)/(1.0 + rho_inf);
787  alpha_f = 1.0/(1.0 + rho_inf);
788  beta = 0.25*pow(1.0 + alpha_m - alpha_f,2);
789  gamma = 0.5 + alpha_m - alpha_f;
790  };
791 
792  void PrintProperties(std::ostream &out = mfem::out);
793 
794  void Init(SecondOrderTimeDependentOperator &f_) override;
795 
796  void Step(Vector &x, Vector &dxdt, double &t, double &dt) override;
797 
798  int GetMaxStateSize() override { return 1; };
799  int GetStateSize() override { return nstate; };
800  const Vector &GetStateVector(int i) override;
801  void GetStateVector(int i, Vector &state) override;
802  void SetStateVector(int i, Vector &state) override;
803 };
804 
805 /// The classical midpoint method.
807 {
808 public:
810  {
811  alpha_m = 0.5;
812  alpha_f = 0.5;
813  beta = 0.25;
814  gamma = 0.5;
815  };
816 };
817 
818 /// HHT-alpha ODE solver
819 /// Improved numerical dissipation for time integration algorithms
820 /// in structural dynamics
821 /// H.M. Hilber, T.J.R. Hughes and R.L. Taylor 1977
822 /// https://doi.org/10.1002/eqe.4290050306
823 /// alpha in [2/3,1] --> Defined differently than in paper.
825 {
826 public:
827  HHTAlphaSolver(double alpha = 1.0)
828  {
829  alpha = (alpha > 1.0) ? 1.0 : alpha;
830  alpha = (alpha < 2.0/3.0) ? 2.0/3.0 : alpha;
831 
832  alpha_m = 1.0;
833  alpha_f = alpha;
834  beta = (2-alpha)*(2-alpha)/4;
835  gamma = 0.5 + alpha_m - alpha_f;
836  };
837 
838 };
839 
840 /// WBZ-alpha ODE solver
841 /// An alpha modification of Newmark's method
842 /// W.L. Wood, M. Bossak and O.C. Zienkiewicz 1980
843 /// https://doi.org/10.1002/nme.1620151011
844 /// rho_inf in [0,1]
846 {
847 public:
848  WBZAlphaSolver(double rho_inf = 1.0)
849  {
850  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
851  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
852 
853  alpha_f = 1.0;
854  alpha_m = 2.0/(1.0 + rho_inf);
855  beta = 0.25*pow(1.0 + alpha_m - alpha_f,2);
856  gamma = 0.5 + alpha_m - alpha_f;
857  };
858 
859 };
860 
861 }
862 
863 #endif
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:629
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1090
virtual ~SecondOrderODESolver()
Definition: ode.hpp:725
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:751
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:409
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:861
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:916
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:47
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:593
int GetMaxStateSize() override
Function for getting and setting the state vectors.
Definition: ode.hpp:345
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:719
int GetStateSize() override
Definition: ode.hpp:799
void SetRhoInf(double rho_inf)
Definition: ode.cpp:850
virtual void SetStateVector(int i, Vector &state)
Definition: ode.hpp:107
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:572
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
virtual int GetMaxStateSize()
Function for getting and setting the state vectors.
Definition: ode.hpp:96
MemoryType mem_type
Definition: ode.hpp:28
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:726
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 int GetStateSize()
Definition: ode.hpp:97
void Step(Vector &x, Vector &dxdt, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:1156
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:938
int GetStateSize() override
Definition: ode.hpp:548
RK2Solver(const double a_=2./3.)
Definition: ode.hpp:141
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:23
ExplicitRKSolver(int s_, const double *a_, const double *b_, const double *c_)
Definition: ode.cpp:138
Vector beta
const Vector & GetStateVector(int i) override
Definition: ode.cpp:377
Abstract class for solving systems of ODEs: d2x/dt2 = f(x,dx/dt,t)
Definition: ode.hpp:627
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:686
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:411
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:778
Second Order Symplectic Integration Algorithm.
Definition: ode.hpp:603
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:995
virtual const Vector & GetStateVector(int i)
Definition: ode.hpp:98
void Step(Vector &x, Vector &dxdt, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:1063
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:815
int GetMaxStateSize() override
Function for getting and setting the state vectors.
Definition: ode.hpp:254
virtual ~ODESolver()
Definition: ode.hpp:112
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1034
virtual void Step(Vector &x, Vector &dxdt, 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 SetStateVector(int i, Vector &state) override
Definition: ode.cpp:1118
GeneralizedAlpha2Solver(double rho_inf=1.0)
Definition: ode.hpp:781
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
The classical midpoint method.
Definition: ode.hpp:806
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
virtual const Vector & GetStateVector(int i)
Definition: ode.hpp:711
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:649
Vector dq_
Definition: ode.hpp:591
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1020
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:492
virtual ~SIASolver()
Definition: ode.hpp:584
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:159
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:888
First Order Symplectic Integration Algorithm.
Definition: ode.hpp:595
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:786
WBZAlphaSolver(double rho_inf=1.0)
Definition: ode.hpp:848
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:30
const Vector & GetStateVector(int i) override
Definition: ode.cpp:825
static bool IsInitialized()
Return true if MPI has been initialized.
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:657
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:743
virtual int GetMaxStateSize()
Function for getting and setting the state vectors.
Definition: ode.hpp:709
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
static MemoryType GetHostMemoryType()
Get the current Host MemoryType. This is the MemoryType used by most MFEM classes when allocating mem...
Definition: device.hpp:264
TimeDependentOperator * F_
Definition: ode.hpp:587
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:587
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:150
Base abstract class for second order time dependent operators.
Definition: operator.hpp:634
NewmarkSolver(double beta_=0.25, double gamma_=0.5)
Definition: ode.hpp:740
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1127
Vector dp_
Definition: ode.hpp:590
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:622
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:1026
AdamsBashforthSolver(int s_, const double *a_)
Definition: ode.cpp:347
Operator * P_
Definition: ode.hpp:588
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:387
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:424
double a
Definition: lissajous.cpp:41
Host memory; using new[] and delete[].
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:501
const Vector & GetStateVector(int i) override
Definition: ode.cpp:1101
virtual void Run(Vector &x, Vector &dxdt, double &t, double &dt, double tf)
Perform time integration from time t [in] to time tf [in].
Definition: ode.hpp:703
SIAVSolver(int order)
Definition: ode.cpp:953
virtual int GetStateSize()
Definition: ode.hpp:710
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:515
const Vector & GetStateVector(int i) override
Definition: ode.cpp:476
virtual void Step(Vector &q, Vector &p, double &t, double &dt)=0
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:693
int GetStateSize() override
Definition: ode.hpp:255
RefCoord t[3]
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:578
GeneralizedAlphaSolver(double rho=1.0)
Definition: ode.hpp:541
virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
Definition: ode.hpp:579
const double alpha
Definition: ex15.cpp:369
virtual void GetStateVector(int i, Vector &state)
Definition: ode.hpp:103
Vector data type.
Definition: vector.hpp:58
int GetMaxStateSize() override
Function for getting and setting the state vectors.
Definition: ode.hpp:798
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:396
virtual void GetStateVector(int i, Vector &state)
Definition: ode.hpp:716
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
virtual void SetStateVector(int i, Vector &state)
Definition: ode.hpp:720
virtual void Run(Vector &x, double &t, double &dt, double tf)
Perform time integration from time t [in] to time tf [in].
Definition: ode.hpp:90
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:27
RefCoord s[3]
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:841
void Step(Vector &q, Vector &p, double &t, double &dt) override
Definition: ode.cpp:925
The classical forward Euler method.
Definition: ode.hpp:117
Abstract operator.
Definition: operator.hpp:24
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:602
void Step(Vector &x, double &t, double &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:76
HHTAlphaSolver(double alpha=1.0)
Definition: ode.hpp:827
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int GetStateSize() override
Definition: ode.hpp:346
AdamsMoultonSolver(int s_, const double *a_)
Definition: ode.cpp:458
int GetMaxStateSize() override
Function for getting and setting the state vectors.
Definition: ode.hpp:547
Variable order Symplectic Integration Algorithm (orders 1-4)
Definition: ode.hpp:611
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:631