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>

#include <time.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;

    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;

    tvector< FieldComponents<GReal_t> > fieldat;

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

    GReal_t lambda;    // <- self coupling constant
    GReal_t Phi0;      // <- equilibrium field value
    // (field 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 d2rhdr2;
    mutable GReal_t dphihdr;
    mutable GReal_t d2phihdr2;
    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;

    bool isSommerfeld;

    int exclude;

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

public:

    class CoordTransf
    {
        private:
            const bool isMinkowski;
            const bool isSchwarzschild;
            const bool isRotating;
            const GReal_t horizono;
            const GReal_t horizoni;
            const GReal_t surfacegravityo;
            const GReal_t surfacegravityi;
            const GReal_t specangmom;
            const int transformedR;
            const int transformedPhi;
            const int rstarTableDiscretization;
            const GReal_t newtonraphsonTolerance;
            const int newtonraphsonMaxIteration;
            tvector<GReal_t> invrstarTable;
            const GReal_t rstarMin;
            const GReal_t rstarMax;
            const GReal_t rstarStep;
        public:
            long double calc_rstar(const long double r);
            GReal_t calc_invrstar_newtonraphson(const GReal_t rstar) throw(IllegalArgumentException&);
            GReal_t calc_invrstar(const GReal_t rstar);
            GReal_t calc_phitilde(const GReal_t r, const GReal_t phi);
            GReal_t calc_invphitilde(const GReal_t r, const GReal_t phitilde);
        public:
            CoordTransf(const GReal_t MArg, const GReal_t aArg, const GReal_t eArg, 
                        const int transformedRArg, const int transformedPhiArg, 
                        const GReal_t minrhArg, const GReal_t maxrhArg, 
                        const int rstarTableDiscretizationArg, 
                        const GReal_t newtonraphsonToleranceArg, 
                        const int newtonraphsonMaxIterationArg);
            ~CoordTransf() {  }
            GReal_t calc_rh(const GReal_t r);
            GReal_t calc_r(const GReal_t rh);
            GReal_t calc_drhdr(const GReal_t r);
            GReal_t calc_d2rhdr2(const GReal_t r);
            GReal_t calc_phih(const GReal_t r, const GReal_t phi);
            GReal_t calc_phi(const GReal_t r, const GReal_t phih);
            GReal_t calc_dphihdr(const GReal_t r);
            GReal_t calc_d2phihdr2(const GReal_t r);
    };
    CoordTransf* coordTransf;
    static inline ScalarFieldMP<GComplex_t> partial_phih(const ScalarFieldMP<GComplex_t>& field)
    {
        return partial_phi(field);
    }
    static inline ScalarFieldMP<GComplex_t> Sqr_partial_phih(const ScalarFieldMP<GComplex_t>& field)
    {
        return Sqr_partial_phi(field);
    }

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

    void setNonevolvingComponents(Grid& G, Grad& d);

    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 PDE* eq;
            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* eqarg, GReal_t Marg, GReal_t aarg, GReal_t earg, GReal_t rh0arg, int directionarg)
             : MarkerQuantity(eqarg), eq(eqarg), M(Marg), a(aarg), e(earg), r0(((KerrHiggs*)eq)->coordTransf->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, ((KerrHiggs*)eq)->coordTransf->calc_rh(r0), direction); }
    };

    MarkerQuantityArray getMarkerQuantities() const;

private:

    static tvector<string> componentNames(int maxOrder);
    static tvector<string> constraintNames(int maxOrder);

};


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


#endif /* gridripper_phys_gr_fixmp_kerrhiggs_KerrHiggs_h */