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 27 MemoryType mem_type; 28 29 public: ODESolver()30 ODESolver() : f(NULL) { mem_type = Device::GetHostMemoryType(); } 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]. */ Run(Vector & x,double & t,double & dt,double tf)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 GetMaxStateSize()95 virtual int GetMaxStateSize() { return 0; } GetStateSize()96 virtual int GetStateSize() { return 0; } GetStateVector(int i)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 } GetStateVector(int i,Vector & state)102 virtual void GetStateVector(int i, Vector &state) 103 { 104 mfem_error("ODESolver has no state vectors"); 105 } SetStateVector(int i,Vector & state)106 virtual void SetStateVector(int i, Vector &state) 107 { 108 mfem_error("ODESolver has no state vectors"); 109 } 110 ~ODESolver()111 virtual ~ODESolver() { } 112 }; 113 114 115 /// The classical forward Euler method 116 class ForwardEulerSolver : public ODESolver 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: RK2Solver(const double a_=2./3.)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 +--------+----------------------+ */ 183 class ExplicitRKSolver : public ODESolver 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. */ 204 class RK6Solver : public ExplicitRKSolver 205 { 206 private: 207 static const double a[28], b[8], c[7]; 208 209 public: RK6Solver()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. */ 216 class RK8Solver : public ExplicitRKSolver 217 { 218 private: 219 static const double a[66], b[12], c[11]; 220 221 public: RK8Solver()222 RK8Solver() : ExplicitRKSolver(12, a, b, c) { } 223 }; 224 225 226 /** An explicit Adams-Bashforth method. */ 227 class AdamsBashforthSolver : public ODESolver 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 GetMaxStateSize()243 int GetMaxStateSize() override { return smax; }; GetStateSize()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 ~AdamsBashforthSolver()249 ~AdamsBashforthSolver() 250 { 251 if (RKsolver) { delete RKsolver; } 252 } 253 }; 254 255 256 /** A 1-stage, 1st order AB method. */ 257 class AB1Solver : public AdamsBashforthSolver 258 { 259 private: 260 static const double a[1]; 261 262 public: AB1Solver()263 AB1Solver() : AdamsBashforthSolver(1, a) { } 264 }; 265 266 /** A 2-stage, 2st order AB method. */ 267 class AB2Solver : public AdamsBashforthSolver 268 { 269 private: 270 static const double a[2]; 271 272 public: AB2Solver()273 AB2Solver() : AdamsBashforthSolver(2, a) { } 274 }; 275 276 /** A 3-stage, 3st order AB method. */ 277 class AB3Solver : public AdamsBashforthSolver 278 { 279 private: 280 static const double a[3]; 281 282 public: AB3Solver()283 AB3Solver() : AdamsBashforthSolver(3, a) { } 284 }; 285 286 /** A 4-stage, 4st order AB method. */ 287 class AB4Solver : public AdamsBashforthSolver 288 { 289 private: 290 static const double a[4]; 291 292 public: AB4Solver()293 AB4Solver() : AdamsBashforthSolver(4, a) { } 294 }; 295 296 /** A 5-stage, 5st order AB method. */ 297 class AB5Solver : public AdamsBashforthSolver 298 { 299 private: 300 static const double a[5]; 301 302 public: AB5Solver()303 AB5Solver() : AdamsBashforthSolver(5, a) { } 304 }; 305 306 307 /** An implicit Adams-Moulton method. */ 308 class AdamsMoultonSolver : public ODESolver 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 GetMaxStateSize()324 int GetMaxStateSize() override { return smax-1; }; GetStateSize()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 ~AdamsMoultonSolver()330 ~AdamsMoultonSolver() 331 { 332 if (RKsolver) { delete RKsolver; } 333 }; 334 }; 335 336 /** A 0-stage, 1st order AM method. */ 337 class AM0Solver : public AdamsMoultonSolver 338 { 339 private: 340 static const double a[1]; 341 342 public: AM0Solver()343 AM0Solver() : AdamsMoultonSolver(0, a) { } 344 }; 345 346 347 /** A 1-stage, 2st order AM method. */ 348 class AM1Solver : public AdamsMoultonSolver 349 { 350 private: 351 static const double a[2]; 352 353 public: AM1Solver()354 AM1Solver() : AdamsMoultonSolver(1, a) { } 355 }; 356 357 /** A 2-stage, 3st order AM method. */ 358 class AM2Solver : public AdamsMoultonSolver 359 { 360 private: 361 static const double a[3]; 362 363 public: AM2Solver()364 AM2Solver() : AdamsMoultonSolver(2, a) { } 365 }; 366 367 /** A 3-stage, 4st order AM method. */ 368 class AM3Solver : public AdamsMoultonSolver 369 { 370 private: 371 static const double a[4]; 372 373 public: AM3Solver()374 AM3Solver() : AdamsMoultonSolver(3, a) { } 375 }; 376 377 /** A 4-stage, 5st order AM method. */ 378 class AM4Solver : public AdamsMoultonSolver 379 { 380 private: 381 static const double a[5]; 382 383 public: AM4Solver()384 AM4Solver() : AdamsMoultonSolver(4, a) { } 385 }; 386 387 388 /// Backward Euler ODE solver. L-stable. 389 class BackwardEulerSolver : public ODESolver 390 { 391 protected: 392 Vector k; 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. 402 class ImplicitMidpointSolver : public ODESolver 403 { 404 protected: 405 Vector k; 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; 424 Vector k, y; 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: 454 Vector k, y; 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. */ 465 class TrapezoidalRuleSolver : public ODESolver 466 { 467 protected: 468 Vector k, y; 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. 508 class GeneralizedAlphaSolver : public ODESolver 509 { 510 protected: 511 mutable Vector xdot,k,y; 512 double alpha_f, alpha_m, gamma; 513 int nstate; 514 515 void SetRhoInf(double rho_inf); 516 void PrintProperties(std::ostream &out = mfem::out); 517 public: 518 GeneralizedAlphaSolver(double rho=1.0)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 GetMaxStateSize()525 int GetMaxStateSize() override { return 1; }; GetStateSize()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 */ 548 class SIASolver 549 { 550 public: SIASolver()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 Run(Vector & q,Vector & p,double & t,double & dt,double tf)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 ~SIASolver()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: SIA1Solver()576 SIA1Solver() {} 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: SIA2Solver()584 SIA2Solver() {} 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) 605 class SecondOrderODESolver 606 { 607 protected: 608 /// Pointer to the associated TimeDependentOperator. 609 SecondOrderTimeDependentOperator *f; // f(.,.,t) : R^n x R^n --> R^n 610 MemoryType mem_type; 611 612 public: SecondOrderODESolver()613 SecondOrderODESolver() : f(NULL) { mem_type = MemoryType::HOST; } 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]. */ Run(Vector & x,Vector & dxdt,double & t,double & dt,double tf)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 GetMaxStateSize()687 virtual int GetMaxStateSize() { return 0; }; GetStateSize()688 virtual int GetStateSize() { return 0; } GetStateVector(int i)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 } GetStateVector(int i,Vector & state)694 virtual void GetStateVector(int i, Vector &state) 695 { 696 mfem_error("ODESolver has no state vectors"); 697 } SetStateVector(int i,Vector & state)698 virtual void SetStateVector(int i, Vector &state) 699 { 700 mfem_error("ODESolver has no state vectors"); 701 } 702 ~SecondOrderODESolver()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. 709 class NewmarkSolver : public SecondOrderODESolver 710 { 711 private: 712 Vector d2xdt2; 713 714 double beta, gamma; 715 bool first; 716 717 public: NewmarkSolver(double beta_=0.25,double gamma_=0.5)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 727 class LinearAccelerationSolver : public NewmarkSolver 728 { 729 public: LinearAccelerationSolver()730 LinearAccelerationSolver() : NewmarkSolver(1.0/6.0, 0.5) { }; 731 }; 732 733 class CentralDifferenceSolver : public NewmarkSolver 734 { 735 public: CentralDifferenceSolver()736 CentralDifferenceSolver() : NewmarkSolver(0.0, 0.5) { }; 737 }; 738 739 class FoxGoodwinSolver : public NewmarkSolver 740 { 741 public: FoxGoodwinSolver()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] 751 class GeneralizedAlpha2Solver : public SecondOrderODESolver 752 { 753 protected: 754 Vector xa,va,aa,d2xdt2; 755 double alpha_f, alpha_m, beta, gamma; 756 int nstate; 757 758 public: GeneralizedAlpha2Solver(double rho_inf=1.0)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 GetMaxStateSize()776 int GetMaxStateSize() override { return 1; }; GetStateSize()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. 784 class AverageAccelerationSolver : public GeneralizedAlpha2Solver 785 { 786 public: AverageAccelerationSolver()787 AverageAccelerationSolver() 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. 802 class HHTAlphaSolver : public GeneralizedAlpha2Solver 803 { 804 public: HHTAlphaSolver(double alpha=1.0)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] 823 class WBZAlphaSolver : public GeneralizedAlpha2Solver 824 { 825 public: WBZAlphaSolver(double rho_inf=1.0)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 842