MFEM  v4.4.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-2022, 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  delete [] k;
253  }
254 };
255 
256 
257 /** A 1-stage, 1st order AB method. */
259 {
260 private:
261  static const double a[1];
262 
263 public:
265 };
266 
267 /** A 2-stage, 2st order AB method. */
269 {
270 private:
271  static const double a[2];
272 
273 public:
275 };
276 
277 /** A 3-stage, 3st order AB method. */
279 {
280 private:
281  static const double a[3];
282 
283 public:
285 };
286 
287 /** A 4-stage, 4st order AB method. */
289 {
290 private:
291  static const double a[4];
292 
293 public:
295 };
296 
297 /** A 5-stage, 5st order AB method. */
299 {
300 private:
301  static const double a[5];
302 
303 public:
305 };
306 
307 
308 /** An implicit Adams-Moulton method. */
310 {
311 private:
312  int s, smax;
313  const double *a;
314  Vector *k;
315  Array<int> idx;
316  ODESolver *RKsolver;
317 
318 public:
319  AdamsMoultonSolver(int s_, const double *a_);
320 
321  void Init(TimeDependentOperator &f_) override;
322 
323  void Step(Vector &x, double &t, double &dt) override;
324 
325  int GetMaxStateSize() override { return smax-1; };
326  int GetStateSize() override { return s-1; };
327  const Vector &GetStateVector(int i) override;
328  void GetStateVector(int i, Vector &state) override;
329  void SetStateVector(int i, Vector &state) override;
330 
332  {
333  if (RKsolver) { delete RKsolver; }
334  delete [] k;
335  };
336 };
337 
338 /** A 0-stage, 1st order AM method. */
340 {
341 private:
342  static const double a[1];
343 
344 public:
346 };
347 
348 
349 /** A 1-stage, 2st order AM method. */
351 {
352 private:
353  static const double a[2];
354 
355 public:
357 };
358 
359 /** A 2-stage, 3st order AM method. */
361 {
362 private:
363  static const double a[3];
364 
365 public:
367 };
368 
369 /** A 3-stage, 4st order AM method. */
371 {
372 private:
373  static const double a[4];
374 
375 public:
377 };
378 
379 /** A 4-stage, 5st order AM method. */
381 {
382 private:
383  static const double a[5];
384 
385 public:
387 };
388 
389 
390 /// Backward Euler ODE solver. L-stable.
392 {
393 protected:
395 
396 public:
397  void Init(TimeDependentOperator &f_) override;
398 
399  void Step(Vector &x, double &t, double &dt) override;
400 };
401 
402 
403 /// Implicit midpoint method. A-stable, not L-stable.
405 {
406 protected:
408 
409 public:
410  void Init(TimeDependentOperator &f_) override;
411 
412  void Step(Vector &x, double &t, double &dt) override;
413 };
414 
415 
416 /** Two stage, singly diagonal implicit Runge-Kutta (SDIRK) methods;
417  the choices for gamma_opt are:
418  0 - 3rd order method, not A-stable
419  1 - 3rd order method, A-stable, not L-stable (default)
420  2 - 2nd order method, L-stable
421  3 - 2nd order method, L-stable (has solves outside [t,t+dt]). */
422 class SDIRK23Solver : public ODESolver
423 {
424 protected:
425  double gamma;
427 
428 public:
429  SDIRK23Solver(int gamma_opt = 1);
430 
431  void Init(TimeDependentOperator &f_) override;
432 
433  void Step(Vector &x, double &t, double &dt) override;
434 };
435 
436 
437 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
438  order 4. A-stable, not L-stable. */
439 class SDIRK34Solver : public ODESolver
440 {
441 protected:
442  Vector k, y, z;
443 
444 public:
445  void Init(TimeDependentOperator &f_) override;
446 
447  void Step(Vector &x, double &t, double &dt) override;
448 };
449 
450 
451 /** Three stage, singly diagonal implicit Runge-Kutta (SDIRK) method of
452  order 3. L-stable. */
453 class SDIRK33Solver : public ODESolver
454 {
455 protected:
457 
458 public:
459  void Init(TimeDependentOperator &f_) override;
460 
461  void Step(Vector &x, double &t, double &dt) override;
462 };
463 
464 
465 /** Two stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
466  of order 2. A-stable. */
468 {
469 protected:
471 
472 public:
473  virtual void Init(TimeDependentOperator &f_);
474 
475  virtual void Step(Vector &x, double &t, double &dt);
476 };
477 
478 
479 /** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
480  of order 2. L-stable. */
481 class ESDIRK32Solver : public ODESolver
482 {
483 protected:
484  Vector k, y, z;
485 
486 public:
487  virtual void Init(TimeDependentOperator &f_);
488 
489  virtual void Step(Vector &x, double &t, double &dt);
490 };
491 
492 
493 /** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
494  of order 3. A-stable. */
495 class ESDIRK33Solver : public ODESolver
496 {
497 protected:
498  Vector k, y, z;
499 
500 public:
501  virtual void Init(TimeDependentOperator &f_);
502 
503  virtual void Step(Vector &x, double &t, double &dt);
504 };
505 
506 
507 /// Generalized-alpha ODE solver from "A generalized-α method for integrating
508 /// the filtered Navier-Stokes equations with a stabilized finite element
509 /// method" by K.E. Jansen, C.H. Whiting and G.M. Hulbert.
511 {
512 protected:
513  mutable Vector xdot,k,y;
515  int nstate;
516 
517  void SetRhoInf(double rho_inf);
518  void PrintProperties(std::ostream &out = mfem::out);
519 public:
520 
521  GeneralizedAlphaSolver(double rho = 1.0) { SetRhoInf(rho); };
522 
523  void Init(TimeDependentOperator &f_) override;
524 
525  void Step(Vector &x, double &t, double &dt) override;
526 
527  int GetMaxStateSize() override { return 1; };
528  int GetStateSize() override { return nstate; };
529  const Vector &GetStateVector(int i) override;
530  void GetStateVector(int i, Vector &state) override;
531  void SetStateVector(int i, Vector &state) override;
532 };
533 
534 
535 /// The SIASolver class is based on the Symplectic Integration Algorithm
536 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
537 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,
538 /// Vol. 92, pages 230-256 (1991).
539 
540 /** The Symplectic Integration Algorithm (SIA) is designed for systems of first
541  order ODEs derived from a Hamiltonian.
542  H(q,p,t) = T(p) + V(q,t)
543  Which leads to the equations:
544  dq/dt = dT/dp
545  dp/dt = -dV/dq
546  In the integrator the operators P and F are defined to be:
547  P = dT/dp
548  F = -dV/dq
549  */
551 {
552 public:
553  SIASolver() : F_(NULL), P_(NULL) {}
554 
555  virtual void Init(Operator &P, TimeDependentOperator & F);
556 
557  virtual void Step(Vector &q, Vector &p, double &t, double &dt) = 0;
558 
559  virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
560  {
561  while (t < tf) { Step(q, p, t, dt); }
562  }
563 
564  virtual ~SIASolver() {}
565 
566 protected:
567  TimeDependentOperator * F_; // p_{i+1} = p_{i} + dt F(q_{i})
568  Operator * P_; // q_{i+1} = q_{i} + dt P(p_{i+1})
569 
570  mutable Vector dp_;
571  mutable Vector dq_;
572 };
573 
574 /// First Order Symplectic Integration Algorithm
575 class SIA1Solver : public SIASolver
576 {
577 public:
579  void Step(Vector &q, Vector &p, double &t, double &dt) override;
580 };
581 
582 /// Second Order Symplectic Integration Algorithm
583 class SIA2Solver : public SIASolver
584 {
585 public:
587  void Step(Vector &q, Vector &p, double &t, double &dt) override;
588 };
589 
590 /// Variable order Symplectic Integration Algorithm (orders 1-4)
591 class SIAVSolver : public SIASolver
592 {
593 public:
594  SIAVSolver(int order);
595  void Step(Vector &q, Vector &p, double &t, double &dt) override;
596 
597 private:
598  int order_;
599 
600  Array<double> a_;
601  Array<double> b_;
602 };
603 
604 
605 
606 /// Abstract class for solving systems of ODEs: d2x/dt2 = f(x,dx/dt,t)
608 {
609 protected:
610  /// Pointer to the associated TimeDependentOperator.
611  SecondOrderTimeDependentOperator *f; // f(.,.,t) : R^n x R^n --> R^n
613 
614 public:
616 
617  /// Associate a TimeDependentOperator with the ODE solver.
618  /** This method has to be called:
619  - Before the first call to Step().
620  - When the dimensions of the associated TimeDependentOperator change.
621  - When a time stepping sequence has to be restarted.
622  - To change the associated TimeDependentOperator. */
623  virtual void Init(SecondOrderTimeDependentOperator &f);
624 
625  /** @brief Perform a time step from time @a t [in] to time @a t [out] based
626  on the requested step size @a dt [in]. */
627  /** @param[in,out] x Approximate solution.
628  @param[in,out] dxdt Approximate rate.
629  @param[in,out] t Time associated with the
630  approximate solution @a x and rate @ dxdt
631  @param[in,out] dt Time step size.
632 
633  The following rules describe the common behavior of the method:
634  - The input @a x [in] is the approximate solution for the input time
635  @a t [in].
636  - The input @a dxdt [in] is the approximate rate for the input time
637  @a t [in].
638  - The input @a dt [in] is the desired time step size, defining the desired
639  target time: t [target] = @a t [in] + @a dt [in].
640  - The output @a x [out] is the approximate solution for the output time
641  @a t [out].
642  - The output @a dxdt [out] is the approximate rate for the output time
643  @a t [out].
644  - The output @a dt [out] is the last time step taken by the method which
645  may be smaller or larger than the input @a dt [in] value, e.g. because
646  of time step control.
647  - The method may perform more than one time step internally; in this case
648  @a dt [out] is the last internal time step size.
649  - The output value of @a t [out] may be smaller or larger than
650  t [target], however, it is not smaller than @a t [in] + @a dt [out], if
651  at least one internal time step was performed.
652  - The value @a x [out] may be obtained by interpolation using internally
653  stored data.
654  - In some cases, the contents of @a x [in] may not be used, e.g. when
655  @a x [out] from a previous Step() call was obtained by interpolation.
656  - In consecutive calls to this method, the output @a t [out] of one
657  Step() call has to be the same as the input @a t [in] to the next
658  Step() call.
659  - If the previous rule has to be broken, e.g. to restart a time stepping
660  sequence, then the ODE solver must be re-initialized by calling Init()
661  between the two Step() calls. */
662  virtual void Step(Vector &x, Vector &dxdt, double &t, double &dt) = 0;
663 
664  /// Perform time integration from time @a t [in] to time @a tf [in].
665  /** @param[in,out] x Approximate solution.
666  @param[in,out] dxdt Approximate rate.
667  @param[in,out] t Time associated with the approximate solution @a x.
668  @param[in,out] dt Time step size.
669  @param[in] tf Requested final time.
670 
671  The default implementation makes consecutive calls to Step() until
672  reaching @a tf.
673  The following rules describe the common behavior of the method:
674  - The input @a x [in] is the approximate solution for the input time
675  @a t [in].
676  - The input @a dxdt [in] is the approximate rate for the input time
677  @a t [in].
678  - The input @a dt [in] is the initial time step size.
679  - The output @a dt [out] is the last time step taken by the method which
680  may be smaller or larger than the input @a dt [in] value, e.g. because
681  of time step control.
682  - The output value of @a t [out] is not smaller than @a tf [in]. */
683  virtual void Run(Vector &x, Vector &dxdt, double &t, double &dt, double tf)
684  {
685  while (t < tf) { Step(x, dxdt, t, dt); }
686  }
687 
688  /// Function for getting and setting the state vectors
689  virtual int GetMaxStateSize() { return 0; };
690  virtual int GetStateSize() { return 0; }
691  virtual const Vector &GetStateVector(int i)
692  {
693  mfem_error("ODESolver has no state vectors");
694  Vector *s = NULL; return *s; // Make some compiler happy
695  }
696  virtual void GetStateVector(int i, Vector &state)
697  {
698  mfem_error("ODESolver has no state vectors");
699  }
700  virtual void SetStateVector(int i, Vector &state)
701  {
702  mfem_error("ODESolver has no state vectors");
703  }
704 
705  virtual ~SecondOrderODESolver() { }
706 };
707 
708 /// The classical newmark method.
709 /// Newmark, N. M. (1959) A method of computation for structural dynamics.
710 /// Journal of Engineering Mechanics, ASCE, 85 (EM3) 67-94.
712 {
713 private:
714  Vector d2xdt2;
715 
716  double beta, gamma;
717  bool first;
718 
719 public:
720  NewmarkSolver(double beta_ = 0.25, double gamma_ = 0.5) { beta = beta_; gamma = gamma_; };
721 
722  void PrintProperties(std::ostream &out = mfem::out);
723 
724  void Init(SecondOrderTimeDependentOperator &f_) override;
725 
726  void Step(Vector &x, Vector &dxdt, double &t, double &dt) override;
727 };
728 
730 {
731 public:
733 };
734 
736 {
737 public:
739 };
740 
742 {
743 public:
744  FoxGoodwinSolver() : NewmarkSolver(1.0/12.0, 0.5) { };
745 };
746 
747 /// Generalized-alpha ODE solver
748 /// A Time Integration Algorithm for Structural Dynamics With Improved
749 /// Numerical Dissipation: The Generalized-α Method
750 /// J.Chung and G.M. Hulbert, J. Appl. Mech 60(2), 371-375, 1993
751 /// https://doi.org/10.1115/1.2900803
752 /// rho_inf in [0,1]
754 {
755 protected:
758  int nstate;
759 
760 public:
761  GeneralizedAlpha2Solver(double rho_inf = 1.0)
762  {
763  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
764  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
765 
766  alpha_m = (2.0 - rho_inf)/(1.0 + rho_inf);
767  alpha_f = 1.0/(1.0 + rho_inf);
768  beta = 0.25*pow(1.0 + alpha_m - alpha_f,2);
769  gamma = 0.5 + alpha_m - alpha_f;
770  };
771 
772  void PrintProperties(std::ostream &out = mfem::out);
773 
774  void Init(SecondOrderTimeDependentOperator &f_) override;
775 
776  void Step(Vector &x, Vector &dxdt, double &t, double &dt) override;
777 
778  int GetMaxStateSize() override { return 1; };
779  int GetStateSize() override { return nstate; };
780  const Vector &GetStateVector(int i) override;
781  void GetStateVector(int i, Vector &state) override;
782  void SetStateVector(int i, Vector &state) override;
783 };
784 
785 /// The classical midpoint method.
787 {
788 public:
790  {
791  alpha_m = 0.5;
792  alpha_f = 0.5;
793  beta = 0.25;
794  gamma = 0.5;
795  };
796 };
797 
798 /// HHT-alpha ODE solver
799 /// Improved numerical dissipation for time integration algorithms
800 /// in structural dynamics
801 /// H.M. Hilber, T.J.R. Hughes and R.L. Taylor 1977
802 /// https://doi.org/10.1002/eqe.4290050306
803 /// alpha in [2/3,1] --> Defined differently than in paper.
805 {
806 public:
807  HHTAlphaSolver(double alpha = 1.0)
808  {
809  alpha = (alpha > 1.0) ? 1.0 : alpha;
810  alpha = (alpha < 2.0/3.0) ? 2.0/3.0 : alpha;
811 
812  alpha_m = 1.0;
813  alpha_f = alpha;
814  beta = (2-alpha)*(2-alpha)/4;
815  gamma = 0.5 + alpha_m - alpha_f;
816  };
817 
818 };
819 
820 /// WBZ-alpha ODE solver
821 /// An alpha modification of Newmark's method
822 /// W.L. Wood, M. Bossak and O.C. Zienkiewicz 1980
823 /// https://doi.org/10.1002/nme.1620151011
824 /// rho_inf in [0,1]
826 {
827 public:
828  WBZAlphaSolver(double rho_inf = 1.0)
829  {
830  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
831  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
832 
833  alpha_f = 1.0;
834  alpha_m = 2.0/(1.0 + rho_inf);
835  beta = 0.25*pow(1.0 + alpha_m - alpha_f,2);
836  gamma = 0.5 + alpha_m - alpha_f;
837  };
838 
839 };
840 
841 }
842 
843 #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:705
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:325
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:688
int GetStateSize() override
Definition: ode.hpp:779
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:285
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:528
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:607
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:391
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:747
Second Order Symplectic Integration Algorithm.
Definition: ode.hpp:583
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:761
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:786
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
virtual const Vector & GetStateVector(int i)
Definition: ode.hpp:691
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:618
Vector dq_
Definition: ode.hpp:571
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:989
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:476
virtual ~SIASolver()
Definition: ode.hpp:564
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:575
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:828
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
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
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:689
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:567
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:603
NewmarkSolver(double beta_=0.25, double gamma_=0.5)
Definition: ode.hpp:720
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:1096
Vector dp_
Definition: ode.hpp:570
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:591
void Init(SecondOrderTimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:995
AdamsBashforthSolver(int s_, const double *a_)
Definition: ode.cpp:347
Operator * P_
Definition: ode.hpp:568
void SetStateVector(int i, Vector &state) override
Definition: ode.cpp:386
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:404
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:683
SIAVSolver(int order)
Definition: ode.cpp:922
virtual int GetStateSize()
Definition: ode.hpp:690
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:521
virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
Definition: ode.hpp:559
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:778
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:696
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:700
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:807
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int GetStateSize() override
Definition: ode.hpp:326
AdamsMoultonSolver(int s_, const double *a_)
Definition: ode.cpp:443
int GetMaxStateSize() override
Function for getting and setting the state vectors.
Definition: ode.hpp:527
Variable order Symplectic Integration Algorithm (orders 1-4)
Definition: ode.hpp:591
SecondOrderTimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:611