TimeDiscretization.h

Go to the documentation of this file.

 
#pragma once
 
#include <vector>
 
#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
#include "Types.h"
 
namespace NumLib
{
 
class TimeDiscretization
{
public:
    TimeDiscretization() = default;
 
    virtual void setInitialState(const double t0) = 0;
 
    double computeRelativeChangeFromPreviousTimestep(
        GlobalVector const& x,
        GlobalVector const& x_old,
        MathLib::VecNormType norm_type);
 
    virtual void nextTimestep(const double t, const double delta_t) = 0;
 
    virtual double getCurrentTime() const = 0;
 
    virtual double getCurrentTimeIncrement() const = 0;
 
    void getXdot(GlobalVector const& x_at_new_timestep,
                 GlobalVector const& x_old,
                 GlobalVector& xdot) const
    {
        namespace LinAlg = MathLib::LinAlg;
 
        auto const dxdot_dx = getNewXWeight();
 
        // xdot = dxdot_dx * x_at_new_timestep - x_old
        getWeightedOldX(xdot, x_old);
        LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
    }
 
    virtual double getNewXWeight() const = 0;
 
    virtual void getWeightedOldX(
        GlobalVector& y, GlobalVector const& x_old) const = 0;  // = x_old
 
    virtual ~TimeDiscretization() = default;
 
protected:
    std::unique_ptr<GlobalVector> _dx;  
};
 
class BackwardEuler final : public TimeDiscretization
{
public:
    void setInitialState(const double t0) override { _t = t0; }
 
    void nextTimestep(const double t, const double delta_t) override
    {
        _t = t;
        _delta_t = delta_t;
    }
 
    double getCurrentTime() const override { return _t; }
    double getCurrentTimeIncrement() const override { return _delta_t; }
    double getNewXWeight() const override { return 1.0 / _delta_t; }
    void getWeightedOldX(GlobalVector& y,
                         GlobalVector const& x_old) const override
    {
        namespace LinAlg = MathLib::LinAlg;
 
        // y = x_old / delta_t
        LinAlg::copy(x_old, y);
        LinAlg::scale(y, 1.0 / _delta_t);
    }
 
private:
    double _t = std::numeric_limits<double>::quiet_NaN();  
    double _delta_t =
        std::numeric_limits<double>::quiet_NaN();  
};
 
}  // namespace NumLib