RK4.h

#ifndef gridripper_amr1d_gridint_RK4_h
#define gridripper_amr1d_gridint_RK4_h

#include <gridripper/amr1d/Integrator.h>
#include <gridripper/amr1d/Dissipation.h>

namespace gridripper { namespace amr1d {

class Grid;
class Grad;

namespace gridint {

using namespace std;

class RK4: public Integrator
{
private:
    Grid* tmpg1;
    Grid* tmpg2;
    Grid* tmpg3;
    Grid* resultGrid;
    Grid* richardsonResultGrid1;
    Grid* richardsonResultGrid2;
    O6Dissipation dissipationOp;

public:
    RK4(PDE& pde, const tvalarray<string>* opts, Sigma* sigma);

protected:
    RK4(const RK4& gi, Grid& g);

private:
    RK4(Grid* g, PDE& pde, Grad* d, Grad* d2, int deftype, Sigma* sigma);

public:
    ~RK4();

    int getNumVelocity() const {
        return 12;
    }

    void integrate(int ir);

    void richardson(const AMError& amerr, int ir);

    Integrator* createIntegrator(Grid& g, const Sigma* s) const;

    Integrator* cloneIntegrator(Grid& g) const;

    void getResult(int i, tvalarray<GReal_t>& w) const;

private:
    void integrate(const Grid& g0, Grid& g1, Grid& g2, Grid& g3, Grid& g4,
                   Grad& grad, Dissipation* dissip,
                   GReal_t dt, int ir);

    void step(GReal_t t, GReal_t dt, const Grid& g0, const Grid& gin,
              Grid& gout, Grad& grad, Dissipation* dissip, int ir);

    void step4(GReal_t t, GReal_t dt,
               const Grid& g0, const Grid& g1, const Grid& g2, const Grid& g3,
               Grid& g4, Grad& grad, Dissipation* dissip, int ir);
};

} } } // namespace gridripper::amr1d::gridint

#endif /* gridripper_amr1d_gridint_RK4_h */