PDE.h

#ifndef gridripper_amr1d_PDE_h
#define gridripper_amr1d_PDE_h

#include <gridripper/config.h>
#include <gridripper/tvalarray.h>
#include <gridripper/tvector.h>
#include "FieldWrapper.h"
#include "GridFunction.h"
#include <map>
#include <string>
#include <gridripper/math/Integrator1D.h>

namespace gridripper { namespace amr1d {

class Grad;
class Grid;

using namespace std;
using namespace gridripper::math;

class PDE
{
public:
    enum {LEFT = 1, RIGHT = 2};

    static const double EPSILON;
//
// Abstract classes and arrays
//

    class DensityQuantity
    {
    private:
        const PDE* that;
    public:
        DensityQuantity(const PDE* eq): that(eq) { }
        virtual ~DensityQuantity() { }
        virtual GReal_t density(GReal_t t, GReal_t x,
                                const tvalarray<GReal_t>& f) const =0;
        virtual GReal_t flow(GReal_t t, GReal_t x,
                             const tvalarray<GReal_t>& f) const =0;
        virtual DensityQuantity* cloneDensityQuantity() const =0;
        const PDE& getPDE() const {
            return *that;
        }
    };

    class PhysicalQuantity
    {
    private:
        const PDE* that;
        tvalarray<string> nameAndFormats;
    public:
        PhysicalQuantity(const PDE* eq, const string& name);

        PhysicalQuantity(const PDE* eq, const string& name,
                         const string& format, const string& headerFormat);

        virtual ~PhysicalQuantity() { }
        virtual GReal_t eval(const Grid& g, GReal_t dt)
            throw(IllegalArgumentException&) =0;
        string getQuantityName() const {
            return nameAndFormats[0];
        }
        string getFormat() const {
            return nameAndFormats[1];
        }
        string getHeaderFormat() const {
            return nameAndFormats[2];
        }
        virtual PhysicalQuantity* clonePhysicalQuantity() const =0;
        const PDE& getPDE() const {
            return *that;
        }
    };

    class PhysicalQuantityArray: public tvalarray<PhysicalQuantity*>
    {
    public:
        PhysicalQuantityArray();
        PhysicalQuantityArray(int n);
        PhysicalQuantityArray(const PhysicalQuantityArray& o);
        PhysicalQuantityArray& operator =(const PhysicalQuantityArray& q);
        ~PhysicalQuantityArray();
    };

    class MarkerQuantity
    {
    private:
        const PDE* that;
    public:
        MarkerQuantity(const PDE* eq): that(eq) { }
        virtual ~MarkerQuantity() { }
        virtual string getQuantityName() const =0;
        virtual GReal_t eval(GReal_t t) =0;
        virtual MarkerQuantity* cloneMarkerQuantity() const =0;
        const PDE& getPDE() const {
            return *that;
        }
    };

    class MarkerQuantityArray: public tvalarray<MarkerQuantity*>
    {
    public:
        MarkerQuantityArray();
        MarkerQuantityArray(int n);
        MarkerQuantityArray(const MarkerQuantityArray& o);
        MarkerQuantityArray& operator =(const MarkerQuantityArray& q);
        ~MarkerQuantityArray();
    };

private:

    int numComponents;

    int numIndependentComponents;

    tvalarray<string> coordinateLabels;

    tvalarray<string> timeFormats;

    string* componentNames;

    string* physicalComponentNames;

    const string* constraintNames;

    map<string, int> componentIndices;

    map<string, int> gridFunctionIndices;

    tvector<GridFunction*> gridFunctions;

    tvalarray<GReal_t> defaultSigmaArray;

    const DensityQuantity* theEnergy;

    bool energySet;

protected:
    PDE(int nc, int nindep, const string* compnames,
        const string* constrnames, const string* coordlabels);

public:
    virtual ~PDE();

    virtual GReal_t getMinX() const =0;

    virtual GReal_t getMaxX() const =0;

    virtual void evalConstrainedComponents(Grid& g, int i,
                                          FieldWrapper& f) =0;

    virtual void setNonevolvingComponents(Grid& g, Grad& d) {
    }

    virtual int getNumConstraintEqs() const {
        return 0;
    }

    virtual void evalConstraintEqs(const Grid& g, int i, const GReal_t* F,
                          GReal_t* lhs, GReal_t* rhs) const;

    int doEval(GReal_t* dF, GReal_t t,
               const GReal_t* F, GReal_t* tmp,
               const Grid& g, int i, Grad& d, GReal_t dt);
//
// Overridable methods
//

    virtual int eval(GReal_t* dF, GReal_t t, GReal_t x,
                     const GReal_t* F, const GReal_t* Fx,
                     const Grid& g, int i, Grad& d, GReal_t dt) =0;

    virtual bool isScalable() const {
        return false;
    }

    virtual GReal_t calcScaleFactor(const Grid* g) const {
        return 1.0;
    }

    virtual FieldWrapper* createField(int level) const =0;

    virtual PhysicalQuantityArray getPhysicalQuantities() const;

    virtual MarkerQuantityArray getMarkerQuantities() const;

    virtual GReal_t getPhysicalMinX() const {
        return getMinX();
    }

    virtual GReal_t getPhysicalMaxX() const {
        return getMaxX();
    }

    virtual bool isXPeriodic() const {
        return false;
    }

    virtual bool canBeExtended(int where) const;

    virtual void mirrorLeft(GReal_t x, FieldWrapper& f) const {
    }

    virtual void mirrorRight(GReal_t x, FieldWrapper& f) const;

    virtual void getExtendedField(const Grid& g, int i, FieldWrapper& f) const;

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

    virtual void getExtendedRefinedField(const Grid& g, int i, FieldWrapper& f)
        const;

    virtual string getPhysicalComponentName(int k) const {
        return physicalComponentNames[k];
    }

    virtual void getFieldWrapper(GReal_t t, GReal_t x, GReal_t* f) const {
    }

    int getNumGridFunctions() const {
        return gridFunctions.size();
    }

    const GridFunction& getGridFunction(int i) const {
        return *gridFunctions[i];
    }

    virtual GReal_t calcBestDeltaT(const Grid& g, int i, const FieldWrapper& F,
                                   GReal_t dt0, GReal_t dt) const {
        return dt;
    }

    virtual const double* getComponentSigmaArray() const {
        return (const double*)defaultSigmaArray;
    }

//
// Non-overridable methods
//
public:

    string getTimeLabel() const {
        return coordinateLabels[0];
    }

    string getTimeFormat() const {
        return timeFormats[0];
    }

    string getTimeLabelFormat() const {
        return timeFormats[1];
    }

    string getXcoordLabel() const {
        return coordinateLabels[1];
    }

    int getNumComponents() const {
        return numComponents;
    }

    int getNumIndependentComponents() const {
        return numIndependentComponents;
    }

    string getComponentName(int k) const {
        return componentNames[k];
    }

    int getComponentIndex(const string& name) const;

    int getGridFunctionIndex(const string& name) const;

    string getConstraintName(int k) const {
        return constraintNames[k];
    }

    double constraintError(Grid& g, int i);

    const DensityQuantity& energy() const {
        return *theEnergy;
    }

protected:
    void setEnergy(const DensityQuantity& e);

    void addGridFunction(GridFunction* f);

//
// Implementations
//
public:

    class ZeroDensity: public PDE::DensityQuantity {
    public:
        ZeroDensity(const PDE* pde): PDE::DensityQuantity(pde) { }
        GReal_t density(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F)
                const {
            return 0;
        }
        GReal_t flow(GReal_t T, GReal_t R, const tvalarray<GReal_t>& F) const {
            return 0;
        }
        PDE::DensityQuantity* cloneDensityQuantity() const {
            return new ZeroDensity(&getPDE());
        }
    };

    class FiniteIntegral: public PhysicalQuantity
    {
    private:
        const DensityQuantity* q;
    public:
        FiniteIntegral(const PDE* eq, const string& name,
                       const DensityQuantity& q);

        FiniteIntegral(const PDE* eq, const string& name,
                       const string& format, const string& headerFormat,
                       const DensityQuantity& q);

        virtual ~FiniteIntegral();

        GReal_t integrate(const Grid& g, GReal_t xmin, GReal_t xmax);
    protected:
        const DensityQuantity& getDensityQuantity() const {
            return *q;
        }
    };

    class TotalIntegral: public FiniteIntegral
    {
    public:
        TotalIntegral(const PDE* eq, const string& name,
                      const DensityQuantity& q);

        TotalIntegral(const PDE* eq, const string& name,
                      const string& format, const string& headerFormat,
                      const DensityQuantity& q);

        virtual PhysicalQuantity* clonePhysicalQuantity() const;

        virtual GReal_t eval(const Grid& g, GReal_t dt)
                throw(IllegalArgumentException&) {
            return integrate(g, getPDE().getPhysicalMinX(),
                             getPDE().getPhysicalMaxX());
        }
    };

    class FiniteX0Integral: public FiniteIntegral
    {
    protected:
        GReal_t xminvalue;
        GReal_t x0value;
    public:
        FiniteX0Integral(const PDE* eq, const string& name,
                         const DensityQuantity& q,
                         GReal_t xmin, GReal_t x0);

        FiniteX0Integral(const PDE* eq, const string& name,
                         const string& format, const string& headerFormat,
                         const DensityQuantity& q,
                         GReal_t xmin, GReal_t x0);

        virtual PhysicalQuantity* clonePhysicalQuantity() const;

        GReal_t eval(const Grid& g, GReal_t dt)
                throw(IllegalArgumentException&) {
            return integrate(g, xminvalue, x0value);
        }
    };

    class TotalX0Integral;

    class TimeIntegral: public PhysicalQuantity
    {
    private:
        PhysicalQuantity* q;
    protected:
        IntegratorIncremental1D* integrator;

        TimeIntegral(const PDE* eq, const string& name,
                     const string& format, const string& headerFormat,
                     PhysicalQuantity& q, const IntegratorIncremental1D& integ);
    public:
        TimeIntegral(const PDE* eq, const string& name, PhysicalQuantity& q);

        virtual ~TimeIntegral();
        virtual PhysicalQuantity* clonePhysicalQuantity() const;
        virtual GReal_t eval(const Grid& g, GReal_t dt)
            throw(IllegalArgumentException&);
    };

    class X0Flow: public PhysicalQuantity
    {
    private:
        const DensityQuantity* q;
    protected:
        GReal_t x0value;
        IntegratorIncremental1D* integrator;

        X0Flow(const PDE* eq, const string& name,
               const string& format, const string& headerFormat,
               const DensityQuantity& q,
               GReal_t x0, const IntegratorIncremental1D& integ);
    public:
        X0Flow(const PDE* eq, const string& name, const DensityQuantity& q,
               GReal_t x0);

        virtual ~X0Flow();
        PhysicalQuantity* clonePhysicalQuantity() const;
        GReal_t eval(const Grid& g, GReal_t dt)
            throw(IllegalArgumentException&);
        const IntegratorIncremental1D& getIntegrator() const {
            return *integrator;
        }
        friend class TotalX0Integral; // must be friend to be able to access
                                      // the protected constructor
    };

    class TotalX0Integral: public FiniteX0Integral
    {
    protected:
        X0Flow flow1;
        X0Flow flow2;

        TotalX0Integral(const PDE* eq, const string& name,
                        const string& format, const string& headerFormat,
                        const DensityQuantity& q,
                        GReal_t xmin, GReal_t x0,
                        const X0Flow& f1,
                        const X0Flow& f2);
    public:
        TotalX0Integral(const PDE* eq, const string& name,
                        const DensityQuantity& q, GReal_t xmin, GReal_t x0);
        virtual PhysicalQuantity* clonePhysicalQuantity() const;
        GReal_t integrate(const Grid& g) {
            return FiniteIntegral::integrate(g, xminvalue, x0value);
        }
        GReal_t eval(const Grid& g, GReal_t dt)
            throw(IllegalArgumentException&);
    };

private:
    GReal_t flow(const DensityQuantity& q, const Grid& g, double t, double x0,
                 tvalarray<GReal_t>& F) const;

    static tvalarray<string> splitNameFormats(const string& name);

    class ComponentGridFunction: public GridFunction {
    private:
        int componentIndex;
    public:
        ComponentGridFunction(const PDE& pde, int k):
                GridFunction(pde.getComponentName(k)) {
            componentIndex = k;
        }
        GReal_t eval(GReal_t t, GReal_t x, const tvalarray<GReal_t>& F) const {
            return F[componentIndex];
        }
    };
    class EnergyGridFunction: public GridFunction {
    private:
        const DensityQuantity* energy;
    public:
        EnergyGridFunction(const DensityQuantity* e): GridFunction("energy") {
            energy = e;
        }
        GReal_t eval(GReal_t t, GReal_t x, const tvalarray<GReal_t>& F) const {
            return energy->density(t, x, F);
        };
    };
};

} } // namespace gridripper::amr1d

#endif /* gridripper_amr1d_PDE_h */