KGEq.h

#ifndef gridripper_phys_minkowski_kg_moncrief_KGEq_h
#define gridripper_phys_minkowski_kg_moncrief_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 moncrief {

using namespace gridripper;
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 kappa;

    GReal_t kappaSq;

    GReal_t mass;

    GReal_t massSq;

    GReal_t 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 0;
    }

    GReal_t getMaxX() const {
        return Rmax;
    }

    GReal_t getPhysicalMaxX() const {
        return 1;
    }

    GReal_t getR0() const {
        return R0;
    }

    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& fw);

    FieldWrapper* createField(int level) const;

    bool canBeExtended(int where) const {
        return where == LEFT
            || where == RIGHT; // zeros beyond right edge in
//                             // P. Csizmadia, IJMPD15 (2006) 107
    }

    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];
    }

    void mirrorRight(GReal_t R, FieldWrapper& f) const {
        // zeros beyond right edge in P. Csizmadia, IJMPD15 (2006) 107
        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 {
        if(R > EPSILON && R <= 1 - EPSILON) {
            GReal_t c = 2*R/(1 - R*R);
            f[I_f] *= c;
            if(f.size() > 1) {
                f[I_fT] *= c;
            }
        } else {
            f[I_f] = 0;
            if(I_fT < f.size()) {
                f[I_fT] = 0;
            }
        }
    }

    PhysicalQuantityArray getPhysicalQuantities() const;
};

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

#endif /* gridripper_phys_minkowski_kg_moncrief_KGEq_h */