TwoFieldsODE.h

#ifndef gridripper_phys_minkowski_ymh_moncrief_TwoFieldsODE_h
#define gridripper_phys_minkowski_ymh_moncrief_TwoFieldsODE_h

#include <gridripper/Parameters.h>
#include <gridripper/amr1d/FieldArray.h>
#include <gridripper/odesolver/ODE.h>
#include <gridripper/lang/IllegalArgumentException.h>

namespace gridripper { namespace phys { namespace minkowski { namespace ymh {
namespace moncrief {

using namespace gridripper;
using namespace gridripper::odesolver;
using namespace gridripper::amr1d;
using namespace std;

class TwoFieldsODE: public ODE
{
private:
    // Field indices:
    enum {I_h  = 0,
          I_w  = 1,
          I_hR = 2,
          I_wR = 3,
          NUM_COMPS = 4};

    // Order of fields in the derivative matrix: h_R, w_R, h, w
    static const int MATORD[4];

    static const string COMPNAMES[NUM_COMPS];

    GReal_t kappa;
    GReal_t lambda;
    GReal_t wwR0;
    GReal_t hhR0;
    GReal_t hhRR0;
    GReal_t hh1;
    GReal_t hhR1;

    tvalarray<GReal_t> df;

public:
    TwoFieldsODE(const Parameters* p) throw(IllegalArgumentException&);

    FieldArray getInitialValuesForShooting() const;

    tvalarray<GReal_t> getInitialValuesForRelaxation(int nx) const;

    tvalarray<GReal_t> getInitialX() const {
        return tvalarray<GReal_t>(0.0, 1);
    }

    GReal_t getMinX() const {
        return 0;
    }

    GReal_t getMaxX() const {
        return 1;
    }

    void getPhysicalField(GReal_t R, GReal_t* f) const {
        GReal_t hh = f[I_h];
        GReal_t hhR = f[I_hR];
        GReal_t ww = f[I_w];
        GReal_t wwR = f[I_wR];
        f[I_h] = R*(hh - 2/kappa);
        f[I_hR] = R*hhR + hh - 2/kappa;
        f[I_w] = R*ww + 1;
        f[I_wR] = R*wwR + ww;
    }

    void initPDE(GReal_t R, const GReal_t* fode, GReal_t* fpde);

    void constrain(tvalarray<GReal_t>& F);

    void eval(const GReal_t* F, int offset, GReal_t R, GReal_t* dF);

    int getMatrixFieldOrder(int i) const {
        return MATORD[i];
    }

    bool dRHS(const GReal_t* F, int offset, GReal_t R, GReal_t* drhs);

private:

    static GReal_t interpolate_f(GReal_t f0, GReal_t fx0,
                               GReal_t f1, GReal_t fx1, GReal_t x);

    static GReal_t interpolate_fx(GReal_t f0, GReal_t fx0,
                                GReal_t f1, GReal_t fx1, GReal_t x);
};

} } } } } // namespace gridripper::phys::minkowski::ymh::moncrief;

#endif /* gridripper_phys_minkowski_ymh_moncrief_TwoFieldsODE_h */