HiggsEq.h

#ifndef gridripper_phys_gr_fixmp_higgsmp_HiggsEq_h
#define gridripper_phys_gr_fixmp_higgsmp_HiggsEq_h


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


namespace gridripper { namespace phys { namespace gr { namespace fixmp {
namespace higgsmp {


using namespace gridripper;
using namespace gridripper::metric;
using namespace gridripper::amr1d;
using namespace std;


class HiggsEq: public PDE
{
public:
    static const GReal_t EPSILON=1e-8;

    /* Field component indices. */
    enum { ind_f=0, ind_ftau=1, ind_frho=2 };

    const int maxOrder;
    const int arraySize;

    const int I_f;      // <- ind_f*arraySize
    const int I_ftau;   // <- ind_ftau*arraySize
    const int I_frho;   // <- ind_frho*arraySize
    static const int NUM_INDEP_COMPS=2;
    static const int NUM_CONSTRAINTS=1;
    static const int NUM_COMPS=NUM_INDEP_COMPS+NUM_CONSTRAINTS;

    static const string COMP_NAMES[NUM_COMPS];
    static const string CONSTRAINT_NAMES[NUM_CONSTRAINTS];
    static const string COORD_LABELS[2];

private:

    tvalarray<GReal_t> df;

    Metric *metric;
    mutable MetricValues g;

    GReal_t lambda;    // <- self coupling constant
    GReal_t phi0;      // <- equilibrium field value
    // (mass=2*sqrt(lambda)*phi0)
    bool selfInteraction;
    // The field is f=(phi-phi0)*r in the simulation if selfinteration 
    // is on, otherwise f=phi*r.

    GReal_t extensionRatio;
    GReal_t EmeasurementRatio;

    tvector<GReal_t> param;

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

public:
    HiggsEq(const Parameters* p, const int maxorder, const int arraysize) throw(IllegalArgumentException&);
    ~HiggsEq();

    GReal_t getMinX() const;
    GReal_t getMaxX() const;
    GReal_t getPhysicalMinX() const;
    GReal_t getPhysicalMaxX() const;
    bool isXPeriodic() const;

    int eval(GReal_t* dF, GReal_t tau, GReal_t rho,
             const GReal_t* F, const GReal_t* F_rho,
             const Grid& G, int i, Grad& d, GReal_t dtau);

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

    FieldWrapper* createField(int level) const;

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

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

    bool canBeExtended(int where) const;

    void mirrorLeft(GReal_t rho, FieldWrapper& f) const;

    void mirrorRight(GReal_t rho, FieldWrapper& f) const;

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

    PhysicalQuantityArray getPhysicalQuantities() const;

    class Eextract : public MarkerQuantity
    {
        private:
            GReal_t position;
        public:
            Eextract(const PDE* eq, GReal_t Position)
             : MarkerQuantity(eq), position(Position)
            {  }
            ~Eextract() {  }
            std::string getQuantityName() const { return "Eextract"; }
            GReal_t eval(GReal_t tau);
            MarkerQuantity* cloneMarkerQuantity() const
            { return new Eextract(&getPDE(), position); }
    };

    class Lightray : public MarkerQuantity
    {
        private:
            Metric* metric;
            GReal_t tau0;
            GReal_t rho0;
            GReal_t kappa;
            int direction;
        public:
            Lightray(const PDE* eq, Metric* g, GReal_t Tau0, GReal_t Rho0, GReal_t Kappa, int Direction)
             : MarkerQuantity(eq), metric(g), tau0(Tau0), rho0(Rho0), kappa(Kappa), direction(Direction)
            {  }
            ~Lightray() {  }
            std::string getQuantityName() const { return "Lightray"; }
            GReal_t eval(GReal_t tau);
            MarkerQuantity* cloneMarkerQuantity() const
            { return new Lightray(&getPDE(), metric, tau0, rho0, kappa, direction); }
    };

    class Horizon : public MarkerQuantity
    {
        private:
            Metric* metric;
            GReal_t tau0;
            GReal_t rho0;
            GReal_t kappa;
            GReal_t Mass;
            int direction;
        public:
            Horizon(const PDE* eq, Metric* g, GReal_t Tau0, GReal_t Rho0, GReal_t Kappa, GReal_t mass, int Direction)
             : MarkerQuantity(eq), metric(g), tau0(Tau0), rho0(Rho0), kappa(Kappa), Mass(mass), direction(Direction)
            {  }
            ~Horizon() {  }
            std::string getQuantityName() const { return "Horizon"; }
            GReal_t eval(GReal_t tau);
            MarkerQuantity* cloneMarkerQuantity() const
            { return new Horizon(&getPDE(), metric, tau0, rho0, kappa, Mass, direction); }
    };

    MarkerQuantityArray getMarkerQuantities() const;

};


} } } } } // namespace gridripper::phys::gr::fixmp::higgsmp


#endif /* gridripper_phys_gr_fixmp_higgsmp_HiggsEq_h */