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