Bdf.h

Go to the documentation of this file.

 
#ifndef PBAT_SIM_INTEGRATION_BDF_H
#define PBAT_SIM_INTEGRATION_BDF_H
 
#include "pbat/Aliases.h"
#include "pbat/common/Concepts.h"
#include "pbat/common/ConstexprFor.h"
#include "pbat/common/Modulo.h"
#include "pbat/io/Archive.h"
 
#include <cassert>
#include <exception>
#include <tuple>
 
namespace pbat::sim::integration {

template <class TScalar = Scalar, class TIndex = Index>
class Bdf
{
  public:
    using ScalarType = TScalar; 
    using IndexType  = TIndex;  

    Bdf(int step = 1, int order = 2);
 
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>
        xt; 
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>
        xtilde;   
    IndexType ti; 
    ScalarType h; 

    [[maybe_unused]] auto Order() const { return mOrder; }
    [[maybe_unused]] auto Step() const { return mStep; }
    [[maybe_unused]] auto Dimensions() const { return xt.rows(); }
    [[maybe_unused]] auto TimeStep() const { return h; }
    [[maybe_unused]] auto Inertia(int o) const { return xtilde.col(o); }
    auto State(int k, int o = 0) const -> decltype(xt.col(0));
    auto State(int k, int o = 0) -> decltype(xt.col(0));
    auto CurrentState(int o = 0) const -> decltype(xt.col(0));
    auto CurrentState(int o = 0) -> decltype(xt.col(0));
    [[maybe_unused]] auto Alpha() const { return mAlpha(Eigen::seqN(0, mStep)); }
    [[maybe_unused]] auto Beta() const { return mBeta; }
    [[maybe_unused]] auto BetaTilde() const { return mBeta * h; }

    void SetOrder(int order);
    void SetStep(int step);
    void SetTimeStep(ScalarType dt);
    void ConstructEquations();
    template <class TDerivedX>
    void SetInitialConditions(Eigen::DenseBase<TDerivedX> const& x0);
    template <class... TDerivedX>
    void SetInitialConditions(Eigen::DenseBase<TDerivedX> const&... x0);
    template <class TDerivedX>
    void Step(Eigen::DenseBase<TDerivedX> const& x);
    template <class... TDerivedX>
    void Step(Eigen::DenseBase<TDerivedX> const&... xs);
    [[maybe_unused]] void Tick() { ++ti; }
    void Serialize(io::Archive& archive) const;
    void Deserialize(io::Archive const& archive);
 
  private:
    int mOrder; 
    int mStep;  
    Eigen::Vector<ScalarType, 6>
        mAlpha; 
    ScalarType
        mBeta; 
};
 
template <class TScalar, class TIndex>
Bdf<TScalar, TIndex>::Bdf(int step, int order)
    : xt(), xtilde(), ti(0), h(Scalar(0.01)), mOrder(), mStep(), mAlpha(), mBeta()
{
    SetStep(step);
    SetOrder(order);
}
 
template <class TScalar, class TIndex>
auto Bdf<TScalar, TIndex>::State(int k, int o) const -> decltype(xt.col(0))
{
    if (k < 0 || k > mStep)
    {
        throw std::out_of_range("0 <= k <= s");
    }
    if (o < 0 || o >= mOrder)
    {
        throw std::out_of_range("0 <= o < order");
    }
    auto kt = common::Modulo(ti /*- mStep*/ + k, mStep);
    return xt.col(o * mStep + kt);
}
 
template <class TScalar, class TIndex>
auto Bdf<TScalar, TIndex>::State(int k, int o) -> decltype(xt.col(0))
{
    if (k < 0 || k > mStep)
    {
        throw std::out_of_range("0 <= k <= s");
    }
    if (o < 0 || o >= mOrder)
    {
        throw std::out_of_range("0 <= o < order");
    }
    auto kt = common::Modulo(ti /*- mStep*/ + k, mStep);
    return xt.col(o * mStep + kt);
}
 
template <class TScalar, class TIndex>
auto Bdf<TScalar, TIndex>::CurrentState(int o) const -> decltype(xt.col(0))
{
    return State(mStep - 1, o);
}
 
template <class TScalar, class TIndex>
auto Bdf<TScalar, TIndex>::CurrentState(int o) -> decltype(xt.col(0))
{
    return State(mStep - 1, o);
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::SetOrder(int order)
{
    if (order <= 0)
    {
        throw std::invalid_argument("order > 0");
    }
    mOrder = order;
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::SetStep(int step)
{
    if (step < 1 || step > 6)
    {
        throw std::invalid_argument("0 < s < 7.");
    }
    mStep = step;
    mAlpha.setZero();
    switch (mStep)
    {
        case 1:
            mAlpha(0) = Scalar(-1);
            mBeta     = Scalar(1);
            break;
        case 2:
            mAlpha(0) = Scalar(1) / 3;
            mAlpha(1) = Scalar(-4) / 3;
            mBeta     = Scalar(2) / 3;
            break;
        case 3:
            mAlpha(0) = Scalar(-2) / 11;
            mAlpha(1) = Scalar(9) / 11;
            mAlpha(2) = Scalar(-18) / 11;
            mBeta     = Scalar(6) / 11;
            break;
        case 4:
            mAlpha(0) = Scalar(3) / 25;
            mAlpha(1) = Scalar(-16) / 25;
            mAlpha(2) = Scalar(36) / 25;
            mAlpha(3) = Scalar(-48) / 25;
            mBeta     = Scalar(12) / 25;
            break;
        case 5:
            mAlpha(0) = Scalar(-12) / 137;
            mAlpha(1) = Scalar(75) / 137;
            mAlpha(2) = Scalar(-200) / 137;
            mAlpha(3) = Scalar(300) / 137;
            mAlpha(4) = Scalar(-300) / 137;
            mBeta     = Scalar(60) / 137;
            break;
        case 6:
            mAlpha(0) = Scalar(10) / 147;
            mAlpha(1) = Scalar(-72) / 147;
            mAlpha(2) = Scalar(225) / 147;
            mAlpha(3) = Scalar(-400) / 147;
            mAlpha(4) = Scalar(450) / 147;
            mAlpha(5) = Scalar(-360) / 147;
            mBeta     = Scalar(60) / 147;
            break;
    }
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::SetTimeStep(ScalarType dt)
{
    if (dt <= 0)
    {
        throw std::invalid_argument("dt > 0");
    }
    h = dt;
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::ConstructEquations()
{
    xtilde.setZero();
    for (auto o = 0; o < mOrder; ++o)
        for (auto k = 0; k < mStep; ++k)
            xtilde.col(o) += mAlpha(k) * State(k, o);
}
 
template <class TScalar, class TIndex>
template <class TDerivedX>
inline void Bdf<TScalar, TIndex>::SetInitialConditions(Eigen::DenseBase<TDerivedX> const& x0)
{
    ti     = 0;
    auto n = x0.rows();
    xt.resize(n, mStep * mOrder);
    xtilde.resize(n, mOrder);
    for (auto o = 0; o < mOrder; ++o)
        xt.middleCols(o * mStep, mStep).colwise() = x0.col(o);
}
 
template <class TScalar, class TIndex>
template <class... TDerivedX>
inline void Bdf<TScalar, TIndex>::SetInitialConditions(Eigen::DenseBase<TDerivedX> const&... x0)
{
    auto constexpr nDerivs = sizeof...(TDerivedX);
    assert(nDerivs == mOrder);
    std::tuple<decltype(x0)...> tup{x0...};
    ti     = 0;
    auto n = std::get<0>(tup).rows();
    xt.resize(n, mStep * mOrder);
    xtilde.resize(n, mOrder);
#include "pbat/warning/Push.h"
#include "pbat/warning/SignConversion.h"
    common::ForRange<0, nDerivs>(
        [&]<auto o>() { xt.middleCols(o * mStep, mStep).colwise() = std::get<o>(tup); });
#include "pbat/warning/Pop.h"
}
 
template <class TScalar, class TIndex>
template <class TDerivedX>
inline void Bdf<TScalar, TIndex>::Step(Eigen::DenseBase<TDerivedX> const& x)
{
    Tick();
    for (auto o = 0; o < mOrder; ++o)
        CurrentState(o) = x.col(o);
}
 
template <class TScalar, class TIndex>
template <class... TDerivedX>
inline void Bdf<TScalar, TIndex>::Step(Eigen::DenseBase<TDerivedX> const&... x)
{
    auto constexpr nDerivs = sizeof...(TDerivedX);
    assert(nDerivs == mOrder);
    Tick();
    std::tuple<decltype(x)...> tup{x...};
    common::ForRange<0, nDerivs>([&]<auto o>() { CurrentState(o) = std::get<o>(tup); });
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::Serialize(io::Archive& archive) const
{
    io::Archive bdfArchive = archive["pbat.sim.integration.Bdf"];
    bdfArchive.WriteData("xt", xt);
    bdfArchive.WriteData("xtilde", xtilde);
    bdfArchive.WriteMetaData("h", h);
    bdfArchive.WriteMetaData("ti", ti);
    bdfArchive.WriteMetaData("order", mOrder);
    bdfArchive.WriteMetaData("step", mStep);
    bdfArchive.WriteMetaData("alpha", mAlpha);
    bdfArchive.WriteMetaData("beta", mBeta);
}
 
template <class TScalar, class TIndex>
void Bdf<TScalar, TIndex>::Deserialize(io::Archive const& archive)
{
    io::Archive const group = archive["pbat.sim.integration.Bdf"];
    xt                      = group.ReadData<MatrixX>("xt");
    xtilde                  = group.ReadData<MatrixX>("xtilde");
    h                       = group.ReadMetaData<Scalar>("h");
    ti                      = group.ReadMetaData<Index>("ti");
    mOrder                  = group.ReadMetaData<int>("order");
    mStep                   = group.ReadMetaData<int>("step");
    mAlpha                  = group.ReadMetaData<Vector<6>>("alpha");
    mBeta                   = group.ReadMetaData<Scalar>("beta");
}
 
} // namespace pbat::sim::integration
 
#endif // PBAT_SIM_INTEGRATION_BDF_H