TwoFields.h

#ifndef gridripper_phys_minkowski_ymh_moncrief_TwoFields_h
#define gridripper_phys_minkowski_ymh_moncrief_TwoFields_h

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

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

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

class TwoFields: public PDE
{
public:
    // Field indices:
    enum {I_h  = 0,
          I_w  = 1,
          I_hT = 2,
          I_wT = 3,
          NUM_INDEP_COMPS = 4,
          NUM_COMPS = NUM_INDEP_COMPS + 2,
          I_hR = NUM_INDEP_COMPS,
          I_wR = NUM_INDEP_COMPS + 1};

private:
    // Field bits:
    static const int F_h;
    static const int F_w;
    static const int F_hT;
    static const int F_wT;
    static const int F_ALL_INDEP;
    static const int F_hR;
    static const int F_wR;
    static const int F_ALL_DEP;
    static const int F_ALL;

    static const string COMPNAMES[NUM_COMPS];

    static const string CONSTRAINT_NAMES[2];

    GReal_t kappa;

    GReal_t kappaSq;

    GReal_t lambda;

    GReal_t Rmax;

    GReal_t R0;

    tvalarray<GReal_t> df;

    class Energy: public DensityQuantity {
    private:
        const TwoFields& pde;
    public:
        Energy(const TwoFields* 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:
    TwoFields(const Parameters* p) throw(IllegalArgumentException&);

    GReal_t getMinX() const {
        return 0;
    }

    GReal_t getMaxX() const {
        return Rmax;
    }

    GReal_t getPhysicalMaxX() const {
        return 1;
    }

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

    void mirrorLeft(GReal_t R, FieldWrapper& f) const {
        GReal_t* v = f.data();
        GReal_t Rsq = R*R;
        GReal_t Omega = kappa*(1-Rsq)/2;
        GReal_t Omegasq = Omega*Omega;

        v[I_h]  = v[I_h] + 2*R/Omega; // h(-R) = h(R)+2R/Omega H_inf
        v[I_hT] = v[I_hT]; // h_T(-R) = h_T(R)
        v[I_hR] = -v[I_hR] - kappa*(1+Rsq)/Omegasq; // h_R(-R) = -h_R(R) - ...

        v[I_w]  = v[I_w]; // w(-R) = w(R)
        v[I_wT] = v[I_wT]; // w_T(-R) = w_T(R)
        v[I_wR] = -v[I_wR]; // w_R(-R) = -w_R(R)
    }

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

    PhysicalQuantityArray getPhysicalQuantities() const;

private:
    GReal_t energyMomentumT00(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F)
        const;

    GReal_t energyMomentumT01(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F)
        const;
};

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

#endif /* gridripper_phys_minkowski_ymh_moncrief_TwoFields_h */