odestepper.hpp

Go to the documentation of this file.

#ifndef __FDL_ODESTEPPER_HPP
#define __FDL_ODESTEPPER_HPP

#include <string>
#include <sstream>
#include "macbase.hpp"

namespace fdl {

template <typename scalar_t>
class ODEStepper
{
protected:
        stepper_T rule_; 
        scalar_t tol_; 
        scalar_t err_; 
        unsigned npts_; 
        std::vector<scalar_t> xtmp, k1, k2, k3, k4;
        
public:
        ODEStepper(unsigned npts =1);
        ~ODEStepper() {}
        
        void setRule(stepper_T rule) {rule_ = rule;}
        std::string getRule();
        void setTolerance(scalar_t tol) {tol_ = tol;}
        scalar_t getTolerance() const {return tol_;}
        scalar_t getError() const {return err_;}
        
        void resize(unsigned npts);
        
        void step(MACBase<scalar_t>* macGrid, 
                                          unsigned npts, scalar_t& dt,
                                          const scalar_t x0[], scalar_t x1[]);
};

template <typename scalar_t>
ODEStepper<scalar_t>::ODEStepper(unsigned npts)
: rule_(RK2), tol_(1), xtmp(), k1(), k2(), k3(), k4()
{
        resize(npts);
}

template <typename scalar_t>
void ODEStepper<scalar_t>::resize(unsigned npts)
{
        xtmp.resize(npts * 3);
        k1.resize(npts * 3);
        k2.resize(npts * 3);
        k3.resize(npts * 3);
        k4.resize(npts * 3);
        npts_ = npts;
}

template <typename scalar_t>
std::string ODEStepper<scalar_t>::getRule()
{
        std::ostringstream oss;
        switch (rule_) {
                case ARK23:
                        oss << "Adaptive Runge-Kutta (2, 3)";
                        break;
                case RK2:
                default:
                        oss << "Runge-Kutta order 2";
                        break;
        }
        return oss.str();
}

template <typename scalar_t>
void ODEStepper<scalar_t>::step(MACBase<scalar_t>* macGrid, 
                                                                                                unsigned npts, scalar_t& dt,
                                                                                                const scalar_t x0[], scalar_t x1[])
{
        if (npts > npts_) {
                DEV() << "Resizing ODEStepper scratch space";
                resize(npts);
        }
        
        /* We use the standard notation for Butcher tableaux; see
         <http://www.scholarpedia.org/article/Runge-Kutta_methods>.
         The c_j coefficients are not needed because the RHS is "stationary". 
         The coefficients e_j are for the error term. 
         */
        scalar_t a21, a31, a32, a41, a42, a43;
        scalar_t b1, b2, b3, b4, e1, e2, e3, e4;
        scalar_t delta;
        
        size_t j, dim = npts * 3;
        
        switch (rule_) 
        {                               
                case ARK23: // Bogacki-Shampine method 
                        a21 = scalar_t(1)/2;
                        a31 = scalar_t(0); a32 = scalar_t(3)/4;
                        a41 = scalar_t(2)/9, a42 = scalar_t(1)/3, a43 = scalar_t(4)/9;
                        b1 = a41, b2 = a42, b3 = a43, b4 = scalar_t(0);
                        e1 = b1 - scalar_t(7)/24, 
                        e2 = b2 - scalar_t(1)/4, 
                        e3 = b3 - scalar_t(1)/3,
                        e4 = b4 - scalar_t(1)/8;
                        
                        do {
                                macGrid->interpField(npts, &x0[0], &k1[0]);
                                
                                for (j = 0; j < dim; ++j)
                                        xtmp[j] = x0[j] + dt * a21 * k1[j];
                                // xtmp[j] = x0[j] - dt * a21 * k1[j];
                                macGrid->interpField(npts, &xtmp[0], &k2[0]);
                                
                                for (j = 0; j < dim; ++j)
                                        xtmp[j] = x0[j] + dt * (a31 * k1[j] + a32 * k2[j]);
                                // xtmp[j] = x0[j] - dt * (a31 * k1[j] + a32 * k2[j]);
                                macGrid->interpField(npts, &xtmp[0], &k3[0]);
                                
                                for (j = 0; j < dim; ++j)
                                        xtmp[j] = x0[j] + dt * (a41 * k1[j] + a42 * k2[j] + a43 * k3[j]);
                                //      xtmp[j] = x0[j] - dt * (a41 * k1[j] + a42 * k2[j] + a43 * k3[j]);
                                macGrid->interpField(npts, &xtmp[0], &k4[0]);
                                
                                err_ = scalar_t(0);
                                for (j = 0; j < dim; ++j) {
                                        x1[j] = xtmp[j];
                                        // compute sup-norm |err| on the fly...
                                        delta = dt * fabs(e1 * k1[j] + e2 * k2[j] + e3 * k3[j] + e4 * k4[j]);
                                        if (delta > err_) 
                                                err_ = delta;
                                }
                                if (err_ >= tol_) {
                                        dt /= 2;
                                        DEV() << "Time step adaptively adjusted: dt = " << dt;
                                }
                        } while (err_ >= tol_);
                        break;
                        
                case RK2:
                default: // non-adaptive Runge-Kutta 2. 
                        a21 = scalar_t(1)/2;
                        macGrid->interpField(npts, &x0[0], &k1[0]);
                        
                        for (j = 0; j < dim; ++j) 
                                xtmp[j] = x0[j] + dt * a21 * k1[j];
                        //      xtmp[j] = x0[j] - dt * a21 * k1[j];
                        macGrid->interpField(npts, &xtmp[0], &k2[0]);
                        
                        for (j = 0; j < dim; ++j)
                                x1[j] = x0[j] + dt * k2[j];
                        // x1[j] = x0[j] - dt * k2[j];
                        break;
        }
}       
        
} // namespace fdl

#endif // __FDL_ODESTEPPER_HPP