1 // ============================================================================= 2 // PROJECT CHRONO - http://projectchrono.org 3 // 4 // Copyright (c) 2014 projectchrono.org 5 // All rights reserved. 6 // 7 // Use of this source code is governed by a BSD-style license that can be found 8 // in the LICENSE file at the top level of the distribution and at 9 // http://projectchrono.org/license-chrono.txt. 10 // 11 // ============================================================================= 12 // Authors: Alessandro Tasora, Radu Serban 13 // ============================================================================= 14 15 #ifndef CHTIMESTEPPER_H 16 #define CHTIMESTEPPER_H 17 18 #include <cstdlib> 19 #include "chrono/core/ChApiCE.h" 20 #include "chrono/core/ChMath.h" 21 #include "chrono/serialization/ChArchive.h" 22 #include "chrono/timestepper/ChIntegrable.h" 23 #include "chrono/timestepper/ChState.h" 24 25 namespace chrono { 26 27 /// @addtogroup chrono_timestepper 28 /// @{ 29 30 /// Base class for timesteppers, i.e., time integrators that can advance a system state. 31 /// It operates on systems inherited from ChIntegrable. 32 class ChApi ChTimestepper { 33 public: 34 /// Available methods for time integration (time steppers). 35 enum class Type { 36 EULER_IMPLICIT_LINEARIZED = 0, 37 EULER_IMPLICIT_PROJECTED = 1, 38 EULER_IMPLICIT = 2, 39 TRAPEZOIDAL = 3, 40 TRAPEZOIDAL_LINEARIZED = 4, 41 HHT = 5, 42 HEUN = 6, 43 RUNGEKUTTA45 = 7, 44 EULER_EXPLICIT = 8, 45 LEAPFROG = 9, 46 NEWMARK = 10, 47 CUSTOM = 20 48 }; 49 50 /// Constructor 51 ChTimestepper(ChIntegrable* intgr = nullptr) integrable(intgr)52 : integrable(intgr), T(0), verbose(false), Qc_do_clamp(false), Qc_clamping(0) {} 53 54 /// Destructor ~ChTimestepper()55 virtual ~ChTimestepper() {} 56 57 /// Return type of the integration method. 58 /// Default is CUSTOM. Derived classes should override this function. GetType()59 virtual Type GetType() const { return Type::CUSTOM; } 60 61 /// Performs an integration timestep 62 virtual void Advance(const double dt ///< timestep to advance 63 ) = 0; 64 65 /// Access the lagrangian multipliers, if any. get_L()66 virtual ChVectorDynamic<>& get_L() { return L; } 67 68 /// Set the integrable object. SetIntegrable(ChIntegrable * intgr)69 virtual void SetIntegrable(ChIntegrable* intgr) { integrable = intgr; } 70 71 /// Get the integrable object. GetIntegrable()72 ChIntegrable* GetIntegrable() { return integrable; } 73 74 /// Get the current time. GetTime()75 virtual double GetTime() const { return T; } 76 77 /// Set the current time. SetTime(double mt)78 virtual void SetTime(double mt) { T = mt; } 79 80 /// Turn on/off logging of messages. SetVerbose(bool verb)81 void SetVerbose(bool verb) { verbose = verb; } 82 83 /// Turn on/off clamping on the Qcterm. SetQcDoClamp(bool dc)84 void SetQcDoClamp(bool dc) { Qc_do_clamp = dc; } 85 86 /// Turn on/off clamping on the Qcterm. SetQcClamping(double cl)87 void SetQcClamping(double cl) { Qc_clamping = cl; } 88 89 /// Method to allow serialization of transient data to archives. 90 virtual void ArchiveOUT(ChArchiveOut& archive); 91 92 /// Method to allow de-serialization of transient data from archives. 93 virtual void ArchiveIN(ChArchiveIn& archive); 94 95 protected: 96 ChIntegrable* integrable; 97 double T; 98 99 ChVectorDynamic<> L; 100 101 bool verbose; 102 103 bool Qc_do_clamp; 104 double Qc_clamping; 105 }; 106 107 /// Base class for 1st order timesteppers, that is a time integrator for a ChIntegrable. 108 class ChApi ChTimestepperIorder : public ChTimestepper { 109 protected: 110 ChState Y; 111 ChStateDelta dYdt; 112 113 public: 114 /// Constructor ChTimestepper(intgr)115 ChTimestepperIorder(ChIntegrable* intgr = nullptr) : ChTimestepper(intgr) { SetIntegrable(intgr); } 116 117 /// Destructor ~ChTimestepperIorder()118 virtual ~ChTimestepperIorder() {} 119 120 /// Access the state at current time get_Y()121 virtual ChState& get_Y() { return Y; } 122 123 /// Access the derivative of state at current time get_dYdt()124 virtual ChStateDelta& get_dYdt() { return dYdt; } 125 126 /// Set the integrable object SetIntegrable(ChIntegrable * intgr)127 virtual void SetIntegrable(ChIntegrable* intgr) { 128 ChTimestepper::SetIntegrable(intgr); 129 Y.setZero(1, intgr); 130 dYdt.setZero(1, intgr); 131 } 132 }; 133 134 /// Base class for 2nd order timesteppers, i.e., a time integrator for a ChIntegrableIIorder. 135 /// A ChIntegrableIIorder is a special subclass of integrable objects that have a state comprised 136 /// of position and velocity y={x,v}, and state derivative dy/dt={v,a}, where a=acceleration. 137 class ChApi ChTimestepperIIorder : public ChTimestepper { 138 protected: 139 ChState X; 140 ChStateDelta V; 141 ChStateDelta A; 142 143 public: 144 /// Constructor ChTimestepper(intgr)145 ChTimestepperIIorder(ChIntegrableIIorder* intgr = nullptr) : ChTimestepper(intgr) { SetIntegrable(intgr); } 146 147 /// Destructor ~ChTimestepperIIorder()148 virtual ~ChTimestepperIIorder() {} 149 150 /// Access the state, position part, at current time get_X()151 virtual ChState& get_X() { return X; } 152 153 /// Access the state, speed part, at current time get_V()154 virtual ChStateDelta& get_V() { return V; } 155 156 /// Access the acceleration, at current time get_A()157 virtual ChStateDelta& get_A() { return A; } 158 159 /// Set the integrable object SetIntegrable(ChIntegrableIIorder * intgr)160 virtual void SetIntegrable(ChIntegrableIIorder* intgr) { 161 ChTimestepper::SetIntegrable(intgr); 162 X.setZero(1, intgr); 163 V.setZero(1, intgr); 164 A.setZero(1, intgr); 165 } 166 167 private: 168 using ChTimestepper::SetIntegrable; 169 }; 170 171 /// Base class for implicit solvers (double inheritance) 172 class ChApi ChImplicitTimestepper {}; 173 174 /// Base properties for implicit solvers. 175 /// Such integrators require solution of a nonlinear problem, typically solved 176 /// using an iterative process, up to a desired tolerance. At each iteration, 177 /// a linear system must be solved. 178 class ChApi ChImplicitIterativeTimestepper : public ChImplicitTimestepper { 179 protected: 180 int maxiters; ///< maximum number of iterations 181 double reltol; ///< relative tolerance 182 double abstolS; ///< absolute tolerance (states) 183 double abstolL; ///< absolute tolerance (Lagrange multipliers) 184 185 int numiters; ///< number of iterations 186 int numsetups; ///< number of calls to the solver's Setup function 187 int numsolves; ///< number of calls to the solver's Solve function 188 189 public: ChImplicitIterativeTimestepper()190 ChImplicitIterativeTimestepper() 191 : maxiters(6), reltol(1e-4), abstolS(1e-10), abstolL(1e-10), numiters(0), numsetups(0), numsolves(0) {} ~ChImplicitIterativeTimestepper()192 virtual ~ChImplicitIterativeTimestepper() {} 193 194 /// Set the max number of iterations using the Newton Raphson procedure SetMaxiters(int iters)195 void SetMaxiters(int iters) { maxiters = iters; } 196 /// Get the max number of iterations using the Newton Raphson procedure GetMaxiters()197 double GetMaxiters() { return maxiters; } 198 199 /// Set the relative tolerance. 200 /// This tolerance is optionally used by derived classes in the Newton-Raphson 201 /// convergence test. SetRelTolerance(double rel_tol)202 void SetRelTolerance(double rel_tol) { reltol = rel_tol; } 203 204 /// Set the absolute tolerances. 205 /// These tolerances are optionally used by derived classes in the Newton-Raphson 206 /// convergence test. This version sets separate absolute tolerances for states 207 /// and Lagrange multipliers. SetAbsTolerances(double abs_tolS,double abs_tolL)208 void SetAbsTolerances(double abs_tolS, double abs_tolL) { 209 abstolS = abs_tolS; 210 abstolL = abs_tolL; 211 } 212 213 /// Set the absolute tolerances. 214 /// These tolerances are optionally used by derived classes in the Newton-Raphson 215 /// convergence test. This version sets equal absolute tolerances for states and 216 /// Lagrange multipliers. SetAbsTolerances(double abs_tol)217 void SetAbsTolerances(double abs_tol) { 218 abstolS = abs_tol; 219 abstolL = abs_tol; 220 } 221 222 /// Return the number of iterations. GetNumIterations()223 int GetNumIterations() const { return numiters; } 224 225 /// Return the number of calls to the solver's Setup function. GetNumSetupCalls()226 int GetNumSetupCalls() const { return numsetups; } 227 228 /// Return the number of calls to the solver's Solve function. GetNumSolveCalls()229 int GetNumSolveCalls() const { return numsolves; } 230 231 /// Method to allow serialization of transient data to archives. ArchiveOUT(ChArchiveOut & archive)232 virtual void ArchiveOUT(ChArchiveOut& archive) { 233 // version number 234 archive.VersionWrite(1); 235 // serialize all member data: 236 archive << CHNVP(maxiters); 237 archive << CHNVP(reltol); 238 archive << CHNVP(abstolS); 239 archive << CHNVP(abstolL); 240 } 241 242 /// Method to allow de-serialization of transient data from archives. ArchiveIN(ChArchiveIn & archive)243 virtual void ArchiveIN(ChArchiveIn& archive) { 244 // version number 245 /*int version =*/ archive.VersionRead(); 246 // stream in all member data: 247 archive >> CHNVP(maxiters); 248 archive >> CHNVP(reltol); 249 archive >> CHNVP(abstolS); 250 archive >> CHNVP(abstolL); 251 } 252 }; 253 254 /// Euler explicit timestepper. 255 /// This performs the typical y_new = y+ dy/dt * dt integration with Euler formula. 256 class ChApi ChTimestepperEulerExpl : public ChTimestepperIorder { 257 public: 258 /// Constructors (default empty) ChTimestepperIorder(intgr)259 ChTimestepperEulerExpl(ChIntegrable* intgr = nullptr) : ChTimestepperIorder(intgr) {} 260 261 /// Performs an integration timestep 262 virtual void Advance(const double dt ///< timestep to advance 263 ); 264 }; 265 266 /// Euler explicit timestepper customized for II order. 267 /// (It gives the same results of ChTimestepperEulerExpl, but performs a bit faster because it 268 /// can exploit the special structure of ChIntegrableIIorder) 269 /// This integrator implements the typical Euler scheme: 270 /// x_new = x + v * dt 271 /// v_new = v + a * dt 272 class ChApi ChTimestepperEulerExplIIorder : public ChTimestepperIIorder { 273 protected: 274 ChStateDelta Dv; 275 276 public: 277 /// Constructors (default empty) ChTimestepperIIorder(intgr)278 ChTimestepperEulerExplIIorder(ChIntegrableIIorder* intgr = nullptr) : ChTimestepperIIorder(intgr) {} 279 GetType()280 virtual Type GetType() const override { return Type::EULER_EXPLICIT; } 281 282 /// Performs an integration timestep 283 virtual void Advance(const double dt ///< timestep to advance 284 ) override; 285 }; 286 287 /// Euler semi-implicit timestepper. 288 /// This performs the typical 289 /// v_new = v + a * dt 290 /// x_new = x + v_new * dt 291 /// integration with Euler semi-implicit formula. 292 class ChApi ChTimestepperEulerSemiImplicit : public ChTimestepperIIorder { 293 public: 294 /// Constructors (default empty) ChTimestepperIIorder(intgr)295 ChTimestepperEulerSemiImplicit(ChIntegrableIIorder* intgr = nullptr) : ChTimestepperIIorder(intgr) {} 296 297 /// Performs an integration timestep 298 virtual void Advance(const double dt ///< timestep to advance 299 ); 300 }; 301 302 /// Performs a step of a 4th order explicit Runge-Kutta integration scheme. 303 class ChApi ChTimestepperRungeKuttaExpl : public ChTimestepperIorder { 304 protected: 305 ChState y_new; 306 ChStateDelta Dydt1; 307 ChStateDelta Dydt2; 308 ChStateDelta Dydt3; 309 ChStateDelta Dydt4; 310 311 public: 312 /// Constructors (default empty) ChTimestepperIorder(intgr)313 ChTimestepperRungeKuttaExpl(ChIntegrable* intgr = nullptr) : ChTimestepperIorder(intgr) {} 314 GetType()315 virtual Type GetType() const override { return Type::RUNGEKUTTA45; } 316 317 /// Performs an integration timestep 318 virtual void Advance(const double dt ///< timestep to advance 319 ) override; 320 }; 321 322 /// Performs a step of a Heun explicit integrator. It is like a 2nd Runge Kutta. 323 class ChApi ChTimestepperHeun : public ChTimestepperIorder { 324 protected: 325 ChState y_new; 326 ChStateDelta Dydt1; 327 ChStateDelta Dydt2; 328 329 public: 330 /// Constructors (default empty) ChTimestepperIorder(intgr)331 ChTimestepperHeun(ChIntegrable* intgr = nullptr) : ChTimestepperIorder(intgr) {} 332 GetType()333 virtual Type GetType() const override { return Type::HEUN; } 334 335 /// Performs an integration timestep 336 virtual void Advance(const double dt ///< timestep to advance 337 ) override; 338 }; 339 340 /// Performs a step of a Leapfrog explicit integrator. 341 /// This is a symplectic method, with 2nd order accuracy, at least when F depends on positions only. 342 /// Note: uses last step acceleration: changing or resorting the numbering of DOFs will invalidate it. 343 /// Suggestion: use the ChTimestepperEulerSemiImplicit, it gives the same accuracy with better performance. 344 class ChApi ChTimestepperLeapfrog : public ChTimestepperIIorder { 345 protected: 346 ChStateDelta Aold; 347 348 public: 349 /// Constructors (default empty) ChTimestepperIIorder(intgr)350 ChTimestepperLeapfrog(ChIntegrableIIorder* intgr = nullptr) : ChTimestepperIIorder(intgr) {} 351 GetType()352 virtual Type GetType() const override { return Type::LEAPFROG; } 353 354 /// Performs an integration timestep 355 virtual void Advance(const double dt ///< timestep to advance 356 ) override; 357 }; 358 359 /// Performs a step of Euler implicit for II order systems. 360 class ChApi ChTimestepperEulerImplicit : public ChTimestepperIIorder, public ChImplicitIterativeTimestepper { 361 protected: 362 ChStateDelta Dv; 363 ChVectorDynamic<> Dl; 364 ChState Xnew; 365 ChStateDelta Vnew; 366 ChVectorDynamic<> R; 367 ChVectorDynamic<> Qc; 368 369 public: 370 /// Constructors (default empty) 371 ChTimestepperEulerImplicit(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)372 : ChTimestepperIIorder(intgr), ChImplicitIterativeTimestepper() {} 373 GetType()374 virtual Type GetType() const override { return Type::EULER_IMPLICIT; } 375 376 /// Performs an integration timestep 377 virtual void Advance(const double dt ///< timestep to advance 378 ) override; 379 }; 380 381 /// Performs a step of Euler implicit for II order systems using the Anitescu/Stewart/Trinkle 382 /// single-iteration method, that is a bit like an implicit Euler where one performs only the 383 /// first Newton corrector iteration. 384 /// If using an underlying CCP complementarity solver, this is the typical Anitescu stabilized 385 /// timestepper for DVIs. 386 class ChApi ChTimestepperEulerImplicitLinearized : public ChTimestepperIIorder, public ChImplicitTimestepper { 387 protected: 388 ChStateDelta Vold; 389 ChVectorDynamic<> Dl; 390 ChVectorDynamic<> R; 391 ChVectorDynamic<> Qc; 392 393 public: 394 /// Constructors (default empty) 395 ChTimestepperEulerImplicitLinearized(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)396 : ChTimestepperIIorder(intgr), ChImplicitTimestepper() {} 397 GetType()398 virtual Type GetType() const override { return Type::EULER_IMPLICIT_LINEARIZED; } 399 400 /// Performs an integration timestep 401 virtual void Advance(const double dt ///< timestep to advance 402 ) override; 403 }; 404 405 /// Performs a step of Euler implicit for II order systems using a semi implicit Euler without 406 /// constraint stabilization, followed by a projection. That is: a speed problem followed by a 407 /// position problem that keeps constraint drifting 'closed' by using a projection. 408 /// If using an underlying CCP complementarity solver, this is the typical Tasora stabilized 409 /// timestepper for DVIs. 410 class ChApi ChTimestepperEulerImplicitProjected : public ChTimestepperIIorder, public ChImplicitTimestepper { 411 protected: 412 ChStateDelta Vold; 413 ChVectorDynamic<> Dl; 414 ChVectorDynamic<> R; 415 ChVectorDynamic<> Qc; 416 417 public: 418 /// Constructors (default empty) 419 ChTimestepperEulerImplicitProjected(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)420 : ChTimestepperIIorder(intgr), ChImplicitTimestepper() {} 421 GetType()422 virtual Type GetType() const override { return Type::EULER_IMPLICIT_PROJECTED; } 423 424 /// Performs an integration timestep 425 virtual void Advance(const double dt ///< timestep to advance 426 ) override; 427 }; 428 429 /// Performs a step of trapezoidal implicit for II order systems. 430 /// NOTE this is a modified version of the trapezoidal for DAE: the original derivation would lead 431 /// to a scheme that produces oscillatory reactions in constraints, so this is a modified version 432 /// that is first order in constraint reactions. Use damped HHT or damped Newmark for more advanced options. 433 class ChApi ChTimestepperTrapezoidal : public ChTimestepperIIorder, public ChImplicitIterativeTimestepper { 434 protected: 435 ChStateDelta Dv; 436 ChVectorDynamic<> Dl; 437 ChState Xnew; 438 ChStateDelta Vnew; 439 ChVectorDynamic<> R; 440 ChVectorDynamic<> Rold; 441 ChVectorDynamic<> Qc; 442 443 public: 444 /// Constructors (default empty) 445 ChTimestepperTrapezoidal(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)446 : ChTimestepperIIorder(intgr), ChImplicitIterativeTimestepper() {} 447 GetType()448 virtual Type GetType() const override { return Type::TRAPEZOIDAL; } 449 450 /// Performs an integration timestep 451 virtual void Advance(const double dt ///< timestep to advance 452 ) override; 453 }; 454 455 /// Performs a step of trapezoidal implicit linearized for II order systems. 456 class ChApi ChTimestepperTrapezoidalLinearized : public ChTimestepperIIorder, public ChImplicitIterativeTimestepper { 457 protected: 458 ChStateDelta Dv; 459 ChVectorDynamic<> Dl; 460 ChState Xnew; 461 ChStateDelta Vnew; 462 ChVectorDynamic<> R; 463 ChVectorDynamic<> Rold; 464 ChVectorDynamic<> Qc; 465 466 public: 467 /// Constructors (default empty) 468 ChTimestepperTrapezoidalLinearized(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)469 : ChTimestepperIIorder(intgr), ChImplicitIterativeTimestepper() {} 470 GetType()471 virtual Type GetType() const override { return Type::TRAPEZOIDAL_LINEARIZED; } 472 473 /// Performs an integration timestep 474 virtual void Advance(const double dt ///< timestep to advance 475 ) override; 476 }; 477 478 /// Performs a step of trapezoidal implicit linearized for II order systems. 479 ///*** SIMPLIFIED VERSION -DOES NOT WORK - PREFER ChTimestepperTrapezoidalLinearized 480 class ChApi ChTimestepperTrapezoidalLinearized2 : public ChTimestepperIIorder, public ChImplicitIterativeTimestepper { 481 protected: 482 ChStateDelta Dv; 483 ChState Xnew; 484 ChStateDelta Vnew; 485 ChVectorDynamic<> R; 486 ChVectorDynamic<> Qc; 487 488 public: 489 /// Constructors (default empty) 490 ChTimestepperTrapezoidalLinearized2(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)491 : ChTimestepperIIorder(intgr), ChImplicitIterativeTimestepper() {} 492 493 /// Performs an integration timestep 494 virtual void Advance(const double dt ///< timestep to advance 495 ); 496 }; 497 498 /// Performs a step of Newmark constrained implicit for II order DAE systems. 499 /// See Negrut et al. 2007. 500 class ChApi ChTimestepperNewmark : public ChTimestepperIIorder, public ChImplicitIterativeTimestepper { 501 private: 502 double gamma; 503 double beta; 504 ChStateDelta Da; 505 ChVectorDynamic<> Dl; 506 ChState Xnew; 507 ChStateDelta Vnew; 508 ChStateDelta Anew; 509 ChVectorDynamic<> R; 510 ChVectorDynamic<> Rold; 511 ChVectorDynamic<> Qc; 512 bool modified_Newton; 513 514 public: 515 /// Constructors (default empty) 516 ChTimestepperNewmark(ChIntegrableIIorder* intgr = nullptr) ChTimestepperIIorder(intgr)517 : ChTimestepperIIorder(intgr), ChImplicitIterativeTimestepper() { 518 SetGammaBeta(0.6, 0.3); // default values with some damping, and that works also with DAE constraints 519 modified_Newton = true; // default use modified Newton with jacobian factorization only at beginning 520 } 521 GetType()522 virtual Type GetType() const override { return Type::NEWMARK; } 523 524 /// Set the numerical damping parameter gamma and the beta parameter. 525 /// Gamma: in the [1/2, 1] interval. 526 /// For gamma = 1/2, no numerical damping 527 /// For gamma > 1/2, more damping 528 /// Beta: in the [0, 1] interval. 529 /// For beta = 1/4, gamma = 1/2 -> constant acceleration method 530 /// For beta = 1/6, gamma = 1/2 -> linear acceleration method 531 /// Method is second order accurate only for gamma = 1/2 532 void SetGammaBeta(double mgamma, double mbeta); 533 GetGamma()534 double GetGamma() { return gamma; } 535 GetBeta()536 double GetBeta() { return beta; } 537 538 /// Enable/disable modified Newton. 539 /// If enabled, the Newton matrix is evaluated, assembled, and factorized only once per step. 540 /// If disabled, the Newton matrix is evaluated at every iteration of the nonlinear solver. 541 /// Modified Newton iteration is enabled by default. SetModifiedNewton(bool val)542 void SetModifiedNewton(bool val) { modified_Newton = val; } 543 544 /// Performs an integration timestep 545 virtual void Advance(const double dt ///< timestep to advance 546 ) override; 547 548 /// Method to allow serialization of transient data to archives. 549 virtual void ArchiveOUT(ChArchiveOut& archive) override; 550 551 /// Method to allow de-serialization of transient data from archives. 552 virtual void ArchiveIN(ChArchiveIn& archive) override; 553 }; 554 555 /// @} chrono_timestepper 556 557 } // end namespace chrono 558 559 #endif 560