KerrHiggs.h

#ifndef gridripper_phys_gr_fixmp_kerrhiggs_KerrHiggs_h
#define gridripper_phys_gr_fixmp_kerrhiggs_KerrHiggs_h


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


namespace gridripper { namespace phys { namespace gr { namespace fixmp {
namespace kerrhiggs {


using namespace gridripper;
using namespace gridripper::multipole;
using namespace std;


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

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

    const int maxOrder;
    const int arraySize;

    const int I_f;      // <- ind_f*arraySize
    const int I_ft;     // <- ind_ft*arraySize
    const int I_frh;    // <- ind_frh*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;

    GReal_t M;          // <- mass
    GReal_t a;          // <- angular momentum
    GReal_t e;          // <- charge

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

    GReal_t minrh;
    GReal_t maxrh;
    GReal_t minrhIntq;
    GReal_t maxrhIntq;

    mutable GReal_t r;
    mutable GReal_t Delta;
    mutable ScalarFieldMP<GComplex_t> Sigma;
    mutable ScalarFieldMP<GComplex_t> Gamma;
    mutable GComplex_t H2S2SigmaSigma;
    mutable GComplex_t H2S2GammaGamma;
    mutable GComplex_t scale;
    mutable ScalarFieldMP<GComplex_t> scaledSigma;
    mutable ScalarFieldMP<GComplex_t> scaledGamma;
    ScalarFieldMP<GComplex_t> Y00;
    GReal_t neumannTolerance;
    mutable GReal_t drhdr;
    mutable GReal_t d2rdrh2;
    mutable ScalarFieldMP<GComplex_t> physf;
    mutable ScalarFieldMP<GComplex_t> physft;
    mutable ScalarFieldMP<GComplex_t> physfr;
    mutable ScalarFieldMP<GComplex_t> physfperSigma;

    /* Start time. */
    GReal_t t0;
    GReal_t rh0;

    static GReal_t horizon;
    static bool isMinkowski;
    static int transformedR;

    bool isSommerfeld;

    int exclude;

    int thresholdIndex;

    GReal_t Theta;
    tvector< tvector< tvector<GReal_t> > > coeffstorage1;
    tvector< tvector< tvector<GReal_t> > > coeffstorage2;

public:

    static GReal_t calc_rh(const GReal_t r);
    static GReal_t calc_r(const GReal_t rh);
    static GReal_t calc_drdrh(const GReal_t rh);
    static GReal_t calc_d2rdrh2(const GReal_t rh);

    KerrHiggs(const Parameters* p, const int maxorder, const int arraysize) throw(IllegalArgumentException&);

    ~KerrHiggs();

    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 t, GReal_t rh,
             const GReal_t* F, const GReal_t* F_rh,
             const Grid& G, int i, Grad& d, GReal_t dt);

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

    FieldWrapper* createField(int level) const;

    GReal_t energyDensity(GReal_t t, GReal_t rh,
                          const tvalarray<GReal_t>& F) const;

    GReal_t energyFlow(GReal_t t, GReal_t rh,
                       const tvalarray<GReal_t>& F) const;

    // Limited to [0,theta]x[0,2pi[.
    GReal_t energyFlowLimited(GReal_t t, GReal_t rh, GReal_t theta, 
                       tvector< tvector< tvector<GReal_t> > > *coeffstorage, 
                       const tvalarray<GReal_t>& F) const;

    GReal_t angmomDensity(GReal_t t, GReal_t rh,
                          const tvalarray<GReal_t>& F) const;

    GReal_t angmomFlow(GReal_t t, GReal_t rh,
                       const tvalarray<GReal_t>& F) const;

    // Limited to [0,theta]x[0,2pi[.
    GReal_t angmomFlowLimited(GReal_t t, GReal_t rh, GReal_t theta, 
                       tvector< tvector< tvector<GReal_t> > > *coeffstorage, 
                       const tvalarray<GReal_t>& F) const;

    bool canBeExtended(int where) const;

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

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

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

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

    class EnergyFlowLimited: public DensityQuantity {
    private:
        const KerrHiggs& pde;
        GReal_t theta;
        tvector< tvector< tvector<GReal_t> > > *coeffstorage;
    public:
        EnergyFlowLimited(const KerrHiggs* p, GReal_t Theta, tvector< tvector< tvector<GReal_t> > > *Coeffstorage): DensityQuantity(p), pde(*p), theta(Theta), coeffstorage(Coeffstorage) { }
        GReal_t density(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                        const {
            return 0.0;
        }
        GReal_t flow(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F) const {
            return pde.energyFlowLimited(t, rh, theta, coeffstorage, F);
        }
        DensityQuantity* cloneDensityQuantity() const {
            return new EnergyFlowLimited(&pde, theta, coeffstorage);
        }
        GReal_t getTheta() const { return theta; }
        void setTheta(GReal_t Theta)
        {
            if ( coeffstorage!=NULL )
            {
                if ( coeffstorage->size() )
                    throw IllegalArgumentException("Cannot change theta after initialized coefficient table!", "gridripper_phys_gr_fixmp_kerrhiggs::KerrHiggs::EnergyFlowLimited::setTheta");
            }
            theta=Theta;
        }
    };
    EnergyFlowLimited theEnergyFlowLimited1, theEnergyFlowLimited2;

    class Angmom: public DensityQuantity {
    private:
        const KerrHiggs& pde;
    public:
        Angmom(const KerrHiggs* p): DensityQuantity(p), pde(*p) { }
        GReal_t density(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                        const {
            return pde.angmomDensity(t, rh, F);
        }
        GReal_t flow(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                     const {
            return pde.angmomFlow(t, rh, F);
        }
        DensityQuantity* cloneDensityQuantity() const {
            return new Angmom(&pde);
        }
    };
    Angmom theAngmom;

    class AngmomFlowLimited: public DensityQuantity {
    private:
        const KerrHiggs& pde;
        GReal_t theta;
        tvector< tvector< tvector<GReal_t> > > *coeffstorage;
    public:
        AngmomFlowLimited(const KerrHiggs* p, GReal_t Theta, tvector< tvector< tvector<GReal_t> > > *Coeffstorage): DensityQuantity(p), pde(*p), theta(Theta), coeffstorage(Coeffstorage) { }
        GReal_t density(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                        const {
            return 0.0;
        }
        GReal_t flow(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                     const {
            return pde.angmomFlowLimited(t, rh, theta, coeffstorage, F);
        }
        DensityQuantity* cloneDensityQuantity() const {
            return new AngmomFlowLimited(&pde, theta, coeffstorage);
        }
        GReal_t getTheta() const { return theta; }
        void setTheta(GReal_t Theta)
        {
            if ( coeffstorage!=NULL )
            {
                if ( coeffstorage->size() )
                    throw IllegalArgumentException("Cannot change theta after initialized coefficient table!", "gridripper_phys_gr_fixmp_kerrhiggs::KerrHiggs::AngmomFlowLimited::setTheta");
            }
            theta=Theta;
        }
    };
    AngmomFlowLimited theAngmomFlowLimited1, theAngmomFlowLimited2;

    PhysicalQuantityArray getPhysicalQuantities() const;

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

    class Lightray : public MarkerQuantity
    {
        private:
            const GReal_t M;
            const GReal_t a;
            const GReal_t e;
            const GReal_t r0;
            const int direction;
            bool firststep;
            GReal_t r;
            GReal_t tprev;
        public:
            Lightray(const PDE* eq, GReal_t Marg, GReal_t aarg, GReal_t earg, GReal_t rh0arg, int directionarg)
             : MarkerQuantity(eq), M(Marg), a(aarg), e(earg), r0(calc_r(rh0arg)), direction(directionarg), firststep(true)
            {  }
            ~Lightray() {  }
            std::string getQuantityName() const { return "Lightray"; }
            GReal_t eval(GReal_t t);
            GReal_t rightside(GReal_t r) const;
            MarkerQuantity* cloneMarkerQuantity() const
            { return new Lightray(&getPDE(), M, a, e, calc_rh(r0), direction); }
    };

    MarkerQuantityArray getMarkerQuantities() const;

};


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


#endif /* gridripper_phys_gr_fixmp_kerrhiggs_KerrHiggs_h */