time_integration.h

#ifndef time_integration_h
#define time_integration_h
 
#include "config.h"  // gamma0
 
class timing
    {
public:
    inline timing(const double _tf, const double _dtmin, const double _dtmax)
        : tf(_tf), DTMIN(_dtmin), DTMAX(_dtmax), TAUR(100. * DTMAX), t(0)
        {
        set_dt(sqrt(DTMIN * DTMAX));
        }
 
    const double tf;
 
    const double DTMIN;  // 1e-14;
 
    const double DTMAX;
 
    inline double get_prefactor(void) const { return prefactor; }
 
    inline double get_dt() const { return dt; }
 
    inline double get_dt_phys() const { return dt / gamma0; }
 
    inline void set_dt(const double _dt)
        {
        dt = _dt;
        double t_tilde = _dt / TAUR;
        prefactor = (1. + t_tilde * abs(log( t_tilde)))*THETA*dt;
        }
 
    template <int NPI>
    Eigen::Matrix<double,NPI,1> get_alpha_eff(
            const Eigen::Ref<const Eigen::Matrix<double,NPI,1>> alpha_g,
            const Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff) const
        {
        Eigen::Matrix<double,NPI,1> a_eff = alpha_g;
        const double r(0.1);
 
        for(int npi=0;npi<NPI;npi++)
            {
            const double alpha = alpha_g(npi);
            const double M = 2. * alpha * r / dt;
            const double h = uHeff(npi);
            if (h > 0.)
                {
                if (h > M)
                    a_eff(npi) = alpha*(1.+r);
                else
                    a_eff(npi) = alpha + dt / 2. * h;
                }
            else
                {
                if (h < -M)
                    a_eff(npi) = alpha/(1.+r);
                else
                    a_eff(npi) = alpha / (1. - dt / (2. * alpha) * h);
                }
            }
        return a_eff;
        }
 
 
    inline bool is_dt_TooSmall() const { return (dt < DTMIN); }
 
    inline void inc_t() { t += dt; }
 
    inline double get_t() const { return t; }
 
    inline double get_t_phys() const { return t / gamma0; }
 
    inline double get_TAUR() const { return TAUR; }
 
    inline void set_t(const double _t) { t = _t; }
 
private:
    const double TAUR;
 
    double t;
 
    double dt;
 
    double prefactor;
    };
 
#endif