KGEq.h

#ifndef gridripper_phys_minkowski_kg_rtlh_KGEq_h
#define gridripper_phys_minkowski_kg_rtlh_KGEq_h

#include <gridripper/amr1d/PDE.h>
#include <gridripper/Parameters.h>
#include <gridripper/lang/IllegalArgumentException.h>

namespace gridripper { namespace phys { namespace minkowski { namespace kg {
namespace rtlh {

using namespace gridripper::util;
using namespace gridripper::amr1d;
using namespace std;

class KGEq: public PDE
{
public:
    enum {I_f = 0,
          I_ft = 1,
          I_fr = 2,
          NUM_INDEP_COMPS = 2,
          NUM_COMPS = NUM_INDEP_COMPS + 1};

    // Field bits:
    static const int F_f;
    static const int F_ft;
    static const int F_fr;

    static const string COMPNAMES[NUM_COMPS];

    static const string CONSTRAINT_NAMES[1];

private:
    static const double EPSILON;

    GReal_t mass;

    GReal_t massSq;

    GReal_t rmin, rmax;

    GReal_t r0;

    tvalarray<GReal_t> df;

    class Energy: public DensityQuantity {
    private:
        const KGEq& pde;
    public:
        Energy(const KGEq* p): DensityQuantity(p), pde(*p) { }
        GReal_t density(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F)
                        const;
        GReal_t flow(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F) const;
        DensityQuantity* cloneDensityQuantity() const {
            return new Energy(&pde);
        }
    } theEnergy;

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

    GReal_t getMinX() const {
        return rmin;
    }

    GReal_t getMaxX() const {
        return rmax;
    }

    GReal_t getPhysicalMaxX() const {
        return rmax;
    }

    int eval(GReal_t* dF, GReal_t t, GReal_t r,
             const GReal_t* F, const GReal_t* Fr,
             const Grid& g, int i, Grad& d, GReal_t dt);

    void evalConstrainedComponents(Grid& g, int i, FieldWrapper& v);

    FieldWrapper* createField(int level) const;

    bool canBeExtended(int where) const {
        return where == LEFT;
//          || where == RIGHT; // zeros beyond right edge
    }

    void mirrorLeft(GReal_t R, FieldWrapper& f) const {
        GReal_t* v = f.data();
        v[I_f]  = v[I_f];
        v[I_ft] = v[I_ft];
        v[I_fr] = -v[I_fr]; //<--!!!!
    }

/*    void mirrorRight(GReal_t R, FieldWrapper& f) const {
        // zeros beyond right edge
        GReal_t* v = f.data();
        v[I_f]  = -v[I_f];
        v[I_ft] = -v[I_ft];
    }*/

    bool hasExtendedRefinedFieldData(const Grid& g, int i) const;

    void getGridField(GReal_t t, GReal_t r, tvalarray<GReal_t>& f) const {
//      f[I_f] *= r;
//      f[I_ft] *= r;
//      f[I_fr] = f[I_f] + r*f[I_fr];
    }

    PhysicalQuantityArray getPhysicalQuantities() const;
};

} } } } } // namespace gridripper::phys::minkowski::kg::rtlh;

#endif /* gridripper_phys_minkowski_kg_rtlh_KGEq_h */