newmarkbeta.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_TIMESTEPPING_NEWMARK_BETA_HH
#define DUMUX_TIMESTEPPING_NEWMARK_BETA_HH
 
#include <dumux/common/parameters.hh>
 
namespace Dumux::Experimental {
 
template<class Scalar, class Solution>
class NewmarkBeta
{
    using Block = typename Solution::block_type;
public:
    NewmarkBeta(const Scalar beta, const Scalar gamma)
    : beta_(beta)
    , gamma_(gamma)
    {}
 
    NewmarkBeta()
    : NewmarkBeta(
        getParam<Scalar>("NewmarkBeta.Beta", 0.25),
        getParam<Scalar>("NewmarkBeta.Gamma", 0.5)
    ) {}
 
    void initialize(const Solution& x)
    {
        x_ = x;
 
        // make sure the size is correct, but the values are zero
        v_ = x; v_ = 0.0;
        a_ = x; a_ = 0.0;
    }
 
    void initialize(const Solution& x, const Solution& v, const Solution& a)
    {
        x_ = x;
        v_ = v;
        a_ = a;
    }
 
    void update(const Scalar dt, const Solution& x)
    {
        for (std::size_t i = 0; i < x_.size(); ++i)
        {
            const auto xNew = x[i];
            const auto aNew = acceleration(i, dt, xNew);
            const auto vNew = velocity(i, dt, xNew, aNew);
 
            x_[i] = xNew;
            v_[i] = vNew;
            a_[i] = aNew;
        }
    }
 
    Block acceleration(const std::size_t index, const Scalar dt, const Block& x) const
    {
        const auto& x0 = x_[index];
        const auto& v0 = v_[index];
        const auto& a0 = a_[index];
 
        return (x - x0 - dt*v0)/(beta_*dt*dt) + (beta_ - 0.5)/beta_*a0;
    }
 
    Block velocity(const std::size_t index, const Scalar dt, const Block& x, const Block& a) const
    {
        const auto& v0 = v_[index];
        const auto& a0 = a_[index];
 
        return v0 + dt*((1.0-gamma_)*a0 + gamma_*a);
    }
 
    Block velocity(const std::size_t index, const Scalar dt, const Block& x) const
    {
        const auto a = acceleration(index, dt, x);
        return velocity(index, dt, x, a);
    }
 
    Scalar position(const std::size_t index) const
    { return x_[index]; }
 
    Scalar velocity(const std::size_t index) const
    { return v_[index]; }
 
    Scalar acceleration(const std::size_t index) const
    { return a_[index]; }
 
private:
    Scalar beta_, gamma_; 
    Solution x_, v_, a_; 
};
 
} // end namespace Dumux::Experimental
 
#endif