| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Experimental | ||
| 10 | * \brief Newmark-beta time integration scheme | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_TIMESTEPPING_NEWMARK_BETA_HH | ||
| 14 | #define DUMUX_TIMESTEPPING_NEWMARK_BETA_HH | ||
| 15 | |||
| 16 | #include <dumux/common/parameters.hh> | ||
| 17 | |||
| 18 | namespace Dumux::Experimental { | ||
| 19 | |||
| 20 | /*! | ||
| 21 | * \brief Newmark-beta time integration scheme | ||
| 22 | * \tparam Scalar The scalar type used for the solution | ||
| 23 | * \tparam Solution The solution type used for the variables | ||
| 24 | * \note Newmark (1959) "A method of computation for structural dynamics" | ||
| 25 | * Journal of the Engineering Mechanics Division, | ||
| 26 | * https://doi.org/10.1061/JMCEA3.0000098 | ||
| 27 | * \note This is typically used for PDEs of the form M*a + C*v + K*x = f | ||
| 28 | * with M the mass matrix, C the damping matrix, K the stiffness matrix, | ||
| 29 | * x the displacement, v the velocity and a the acceleration. | ||
| 30 | * An example is dynamic mechanics. | ||
| 31 | */ | ||
| 32 | template<class Scalar, class Solution> | ||
| 33 | class NewmarkBeta | ||
| 34 | { | ||
| 35 | using Block = typename Solution::block_type; | ||
| 36 | public: | ||
| 37 | //! Construct the Newmark-beta time integration scheme | ||
| 38 | //! with the given β and γ parameters | ||
| 39 | 2 | NewmarkBeta(const Scalar beta, const Scalar gamma) | |
| 40 | 2 | : beta_(beta) | |
| 41 | 2 | , gamma_(gamma) | |
| 42 | {} | ||
| 43 | |||
| 44 | //! Construct the Newmark-beta time integration scheme | ||
| 45 | //! where parameters are read from the parameter file | ||
| 46 | 2 | NewmarkBeta() | |
| 47 | : NewmarkBeta( | ||
| 48 | 4 | getParam<Scalar>("NewmarkBeta.Beta", 0.25), | |
| 49 | 4 | getParam<Scalar>("NewmarkBeta.Gamma", 0.5) | |
| 50 | 2 | ) {} | |
| 51 | |||
| 52 | //! Initialize the time integration scheme with the current solution | ||
| 53 | //! while setting the velocity and acceleration to zero | ||
| 54 | 1 | void initialize(const Solution& x) | |
| 55 | { | ||
| 56 | 1 | x_ = x; | |
| 57 | |||
| 58 | // make sure the size is correct, but the values are zero | ||
| 59 | 1 | v_ = x; v_ = 0.0; | |
| 60 | 1 | a_ = x; a_ = 0.0; | |
| 61 | 1 | } | |
| 62 | |||
| 63 | //! Initialize the time integration scheme with the current solution | ||
| 64 | 1 | void initialize(const Solution& x, const Solution& v, const Solution& a) | |
| 65 | { | ||
| 66 | 1 | x_ = x; | |
| 67 | 1 | v_ = v; | |
| 68 | 1 | a_ = a; | |
| 69 | 1 | } | |
| 70 | |||
| 71 | //! Update with new position x | ||
| 72 | 602 | void update(const Scalar dt, const Solution& x) | |
| 73 | { | ||
| 74 |
2/2✓ Branch 0 taken 109202 times.
✓ Branch 1 taken 602 times.
|
109804 | for (std::size_t i = 0; i < x_.size(); ++i) |
| 75 | { | ||
| 76 | 109202 | const auto xNew = x[i]; | |
| 77 | 109202 | const auto aNew = acceleration(i, dt, xNew); | |
| 78 | 109202 | const auto vNew = velocity(i, dt, xNew, aNew); | |
| 79 | |||
| 80 | 109202 | x_[i] = xNew; | |
| 81 | 109202 | v_[i] = vNew; | |
| 82 | 109202 | a_[i] = aNew; | |
| 83 | } | ||
| 84 | 602 | } | |
| 85 | |||
| 86 | //! new a in terms of the old x, v, a, the new x and the time step size | ||
| 87 | 11582997 | Block acceleration(const std::size_t index, const Scalar dt, const Block& x) const | |
| 88 | { | ||
| 89 | 11582997 | const auto& x0 = x_[index]; | |
| 90 | 11582997 | const auto& v0 = v_[index]; | |
| 91 | 11582997 | const auto& a0 = a_[index]; | |
| 92 | |||
| 93 | 69470357 | return (x - x0 - dt*v0)/(beta_*dt*dt) + (beta_ - 0.5)/beta_*a0; | |
| 94 | } | ||
| 95 | |||
| 96 | //! new v in terms of the old v, a, the new x, a and the time step size | ||
| 97 | 109202 | Block velocity(const std::size_t index, const Scalar dt, const Block& x, const Block& a) const | |
| 98 | { | ||
| 99 | 109202 | const auto& v0 = v_[index]; | |
| 100 | 109202 | const auto& a0 = a_[index]; | |
| 101 | |||
| 102 | 544002 | return v0 + dt*((1.0-gamma_)*a0 + gamma_*a); | |
| 103 | } | ||
| 104 | |||
| 105 | //! new v in terms of the old v, a, the new x and the time step size | ||
| 106 | Block velocity(const std::size_t index, const Scalar dt, const Block& x) const | ||
| 107 | { | ||
| 108 | const auto a = acceleration(index, dt, x); | ||
| 109 | return velocity(index, dt, x, a); | ||
| 110 | } | ||
| 111 | |||
| 112 | //! current position | ||
| 113 | 502 | Scalar position(const std::size_t index) const | |
| 114 |
1/2✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
|
502 | { return x_[index]; } |
| 115 | |||
| 116 | //! current velocity | ||
| 117 | 502 | Scalar velocity(const std::size_t index) const | |
| 118 |
1/2✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
|
502 | { return v_[index]; } |
| 119 | |||
| 120 | //! current acceleration | ||
| 121 | 502 | Scalar acceleration(const std::size_t index) const | |
| 122 |
1/2✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
|
502 | { return a_[index]; } |
| 123 | |||
| 124 | private: | ||
| 125 | Scalar beta_, gamma_; //!< Newmark-beta parameters | ||
| 126 | Solution x_, v_, a_; //!< current position, velocity and acceleration | ||
| 127 | }; | ||
| 128 | |||
| 129 | } // end namespace Dumux::Experimental | ||
| 130 | |||
| 131 | #endif | ||
| 132 |