block_time_step.h

Go to the documentation of this file.

#pragma once
 
#include "Common/Float.h"
 
namespace H4{
    // forth order Time step calculator class
    class TimeStep4th{
    public:
        Float eta_4th;  
        Float eta_2nd;  
        Float acc0_offset_sq; 
 
        TimeStep4th(): eta_4th(Float(-1.0)), eta_2nd(Float(-1.0)), acc0_offset_sq(Float(-1.0)) {}
 
 
        bool checkParams() {
            ASSERT(eta_4th>0.0);
            ASSERT(eta_2nd>0.0);
            ASSERT(acc0_offset_sq>=0.0);
            return true;
        }        
 
 
        void calcAcc0OffsetSq(const Float _mass, const Float _r, const Float _G) {
            acc0_offset_sq = _G * _mass / (_r*_r);
            acc0_offset_sq = acc0_offset_sq*acc0_offset_sq;
        }
 
 
 
        Float calcDt2nd(const Float* acc0, 
                        const Float* acc1) const {
            const Float s0 = acc0[0] * acc0[0] + acc0[1] * acc0[1] + acc0[2] * acc0[2] + acc0_offset_sq;
            const Float s1 = acc1[0] * acc1[0] + acc1[1] * acc1[1] + acc1[2] * acc1[2];
            if(s0 == acc0_offset_sq || s1 == 0.0)
                return NUMERIC_FLOAT_MAX;
            else 
                return eta_2nd * sqrt( s0 / s1 );
        }
 
 
        Float calcDt4th(const Float* acc0,
                        const Float* acc1,
                        const Float* acc2,
                        const Float* acc3) const {
            const Float s0 = acc0[0] * acc0[0] + acc0[1] * acc0[1] + acc0[2] * acc0[2] + acc0_offset_sq;
            const Float s1 = acc1[0] * acc1[0] + acc1[1] * acc1[1] + acc1[2] * acc1[2];
            const Float s2 = acc2[0] * acc2[0] + acc2[1] * acc2[1] + acc2[2] * acc2[2];
            const Float s3 = acc3[0] * acc3[0] + acc3[1] * acc3[1] + acc3[2] * acc3[2];
            if(s0 == acc0_offset_sq || s1 == 0.0) 
                return NUMERIC_FLOAT_MAX;
            else 
                return eta_4th * sqrt( (sqrt(s0*s2) + s1) / (sqrt(s1*s3) + s2) );
        }
 
        void print(std::ostream & _fout) const{
            _fout<<"eta_4th : "<<eta_4th<<std::endl
                 <<"eta_2nd : "<<eta_2nd<<std::endl;
        }
 
 
        void printColumnTitle(std::ostream & _fout, const int _width=20) {
            _fout<<std::setw(_width)<<"Eta(4th)"
                 <<std::setw(_width)<<"Eta(2nd)";
        }
 
 
        void printColumn(std::ostream & _fout, const int _width=20){
            _fout<<std::setw(_width)<<eta_4th
                 <<std::setw(_width)<<eta_2nd;
        }
 
 
        void writeBinary(FILE *_fp) const {
            fwrite(this, sizeof(*this),1,_fp);
        }
 
 
        void readBinary(FILE *_fin) {
            size_t rcount = fread(this, sizeof(*this),1,_fin);
            if (rcount<1) {
                std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
                abort();
            }
        }
    };
 
    // Block time step class
    class BlockTimeStep4th: public TimeStep4th{
    private:
        Float dt_max_; // maximum time step
        Float dt_min_; // minimum time step
    public:
        BlockTimeStep4th(): TimeStep4th(), dt_max_(-1.0), dt_min_(-1.0) {}
 
 
        void setDtRange(const Float _dt_max, const int _pow_index_min) {
            dt_max_ = _dt_max;
            dt_min_ = _dt_max;
            for (int i=0; i<_pow_index_min; i++) {
                dt_min_ *= Float(0.5);
            }
        }
 
 
        Float getDtMax() const {
            return dt_max_;
        }
 
 
        Float getDtMin() const {
            return dt_min_;
        }
 
 
        Float calcNextDtLimit(const Float _time) {
            ASSERT(dt_max_>dt_min_);
            ASSERT(dt_min_>0.0);
            // for first step, the maximum time step is OK
            if(_time==0.0) return dt_max_;
            else {
                // the binary tree for current time position in block step 
                unsigned long long int bitmap = to_double(_time/dt_min_);
                //#ifdef __GNUC__ 
                //        PS::S64 dts = __builtin_ctz(bitmap) ;
                //        PS::U64 c = (1<<dts);
                //#else
 
                // block step multiply factor 
                unsigned long long int c=1;
                // find the last zero in the binary tree to obtain the current block step level
                while((bitmap&1)==0) {
                    bitmap = (bitmap>>1);
                    c = (c<<1);
                }
                //#endif
            
                // return the maximum step allown
                return std::min(c*dt_min_,dt_max_);
            }
        }
 
 
        Float calcBlockDt2nd(const Float* _acc0, 
                             const Float* _acc1,
                             const Float _dt_limit) const{
            ASSERT(dt_max_>dt_min_);
            ASSERT(_dt_limit<=dt_max_);
            ASSERT(_dt_limit>=dt_min_);
 
            const Float dt_ref = TimeStep4th::calcDt2nd(_acc0, _acc1);
            Float dt = _dt_limit;
            while(dt > dt_ref) dt *= 0.5;
 
            if(dt<dt_min_) {
                std::cerr<<"Error: time step size too small: ("<<dt<<") < dt_min ("<<dt_min_<<")!"<<std::endl;
                DATADUMP("hard_dump");
                abort();
            }
            else return dt;
        }
 
 
        Float calcBlockDt4th(const Float* acc0, 
                             const Float* acc1,
                             const Float* acc2,
                             const Float* acc3,
                             const Float _dt_limit) const {
            ASSERT(dt_max_>dt_min_);
            ASSERT(_dt_limit<=dt_max_);
            ASSERT(_dt_limit>=dt_min_);
 
            const Float dt_ref = TimeStep4th::calcDt4th(acc0, acc1, acc2, acc3);
            Float dt = _dt_limit;
            while(dt > dt_ref) dt *= 0.5;
 
            if(dt<dt_min_) {
                std::cerr<<"Error: time step size too small: ("<<dt<<") < dt_min ("<<dt_min_<<")!"<<std::endl;
                DATADUMP("hard_dump");
                abort();
            }
            else return dt;
        }
 
        void print(std::ostream & _fout) const{
            TimeStep4th::print(_fout);
            _fout<<"dt_max: "<<dt_max_<<std::endl
                 <<"dt_min: "<<dt_min_<<std::endl;
        }
 
 
        void printColumnTitle(std::ostream & _fout, const int _width=20) {
            TimeStep4th::printColumnTitle(_fout, _width);
            _fout<<std::setw(_width)<<"Dt_max"
                 <<std::setw(_width)<<"Dt_min";
        }
 
 
        void printColumn(std::ostream & _fout, const int _width=20){
            TimeStep4th::printColumn(_fout, _width);
            _fout<<std::setw(_width)<<dt_max_
                 <<std::setw(_width)<<dt_min_;
        }
 
 
        void writeBinary(FILE *_fp) const {
            fwrite(this, sizeof(*this),1,_fp);
        }
 
 
        void readBinary(FILE *_fin) {
            size_t rcount = fread(this, sizeof(*this),1,_fin);
            if (rcount<1) {
                std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
                abort();
            }
        }
    };
 
}