Proto_Stencil.H

#ifndef _PROTO_STENCIL_H_
#define _PROTO_STENCIL_H_
#include "Proto_BoxData.H"
#include "Proto_DisjointBoxLayout.H"

#include <vector>
#include <tuple>
#include <iostream>
#include <iomanip> //for pretty printing
#include <set>

namespace Proto {

//=======================================================================================
// SHIFT ||
//=======++

template <typename T> class Stencil; //forward declaration



class Shift {
public:



    Shift() : m_shift(Point::Zeros()){};


    explicit Shift(const Point& a_pt){m_shift = a_pt;};


    template<typename... vals>
    inline explicit Shift(vals... args) : m_shift(args...) {}




    inline static Shift Basis(int a_dir, int a_scale=1)
    {return Shift(Point::Basis(a_dir, a_scale));};


    inline static Shift Zeros(){return Shift(Point::Zeros());};


    inline static Shift Ones(int a_scale=1){return Shift(Point::Ones(a_scale));};

    inline Point& shift(){return m_shift;}


    template<typename T>
    inline Stencil<T> operator*(T a_coef) const
    {
        return Stencil<T>(*this,a_coef);
    }


    inline Shift operator*(const Shift& a_shift) const
    {
        return Shift(m_shift + a_shift.m_shift);
    }

    inline int& operator[](int a_dir)
    {
      PROTO_ASSERT(((a_dir >= 0) && (a_dir < DIM)),
        "Shift::operator[](int a_dir) invalid for a_dir = %i. a_dir must be in [0,DIM=%i)",
        a_dir, DIM);
        return m_shift[a_dir];
    }

private:
    Point m_shift;
};

//=======================================================================================
// LAZYSTENCIL ||
//=============++

// forward declaration
  template <typename T, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    class BoxData;


  template <typename T, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
class LazyStencil {
public:
    inline LazyStencil() {};
    inline LazyStencil(const Stencil<T>* a_stencil,
                       const BoxData<T,C,MEMTYPE,D,E>* a_src,
                       Box a_box, T a_scale);

    inline void apply(BoxData<T,C,MEMTYPE,D,E>& a_dest, bool a_overwrite);

    inline unsigned long long int size(){return m_stencil.size();}

    Box m_range;

    std::vector<const Stencil<T>*> m_stencil;
    std::vector<BoxData<T, C,MEMTYPE, D, E>*> m_src;
    std::vector<Box> m_box;
    std::vector<T> m_scale;
};

//=======================================================================================
// STENCIL ||
//=========++


template <typename T>
class Stencil {

public:

struct coeff_holder
  {
    T elements[64];
  };

  ~Stencil()
  {
//    using std::cout;
//    using std::endl;
//    cout << "in stencil destructor" << endl;
    if(m_isClosed)
    {
#ifdef PROTO_CUDA
      protoFreeGPU(d_coeff);
      protoFreeGPU(d_offset);
#endif
      m_isClosed = false;
    }
  }
  //Set up the coefficients on the device.
  //Must be called before computations can be performed on the device
  inline void closeForDevice();
  template <typename TT, unsigned int C, MemType MEMTYPE , unsigned char D, unsigned char E>
  friend class LazyStencil;

    template <typename TT, unsigned int C , MemType MEMTYPE, unsigned char D, unsigned char E>
    friend class BoxData;

    template <typename TT, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    friend BoxData<TT,C,MEMTYPE,D,E>& operator|=(BoxData<TT,C,MEMTYPE,D,E>& a_dest,
                         LazyStencil<TT,C,MEMTYPE,D,E>&& a_op);

  template <typename TT, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    friend BoxData<TT,C,MEMTYPE,D,E>& operator+=(BoxData<TT,C,MEMTYPE,D,E>& a_dest,
                                                         LazyStencil<TT,C,MEMTYPE,D,E>&& a_op);



    Stencil();


    Stencil(Shift   a_shift,
            T       a_coef,
            Point   a_destRefratio  = Point::Ones(),
            Point   a_destShift     = Point::Zeros(),
            Point   a_srcRefratio   = Point::Ones());





    Stencil<T> operator*(const Stencil<T>& a_stencil) const;


    Stencil<T> operator*(const T a_coef) const;


    void operator*=(const Stencil<T>& a_stencil);

    void operator*=(const T a_coef);


    Stencil<T> operator+(const Stencil<T>& a_stencil) const;


    Stencil<T> operator-(const Stencil<T>& a_stencil) const;


    void operator+=(const Stencil<T>& a_stencil);


    void operator-=(const Stencil<T>& a_stencil);




    inline bool operator==(Stencil<T>& a_stencil) const;


    inline bool operator!=(Stencil<T>& a_stencil) const {return !(*this == a_stencil);}


    inline const std::vector<T>& coefs() const {return m_coefs;};


    inline const T* devCoefs() const {return d_coeff;};


    inline const Point* devOffsets() const {return d_offset;};


    inline const std::vector<Point>& offsets() const {return m_offsets;};


    inline unsigned long long int size() const {return m_coefs.size();}


    inline Box span() const {return m_span;};


    inline Point spanPoint() const;


    inline Point ghost() const;


    inline unsigned long long int numFlops(const Box& a_box) const;



    inline Point& srcRatio(){return m_srcRefratio;};

    inline const Point& srcRatio() const {return m_srcRefratio;};


    inline Point& destRatio(){return m_destRefratio;};

    inline const Point& destRatio() const {return m_destRefratio;};


    inline Point& destShift(){return m_destShift;};


    inline const Point& destShift() const {return m_destShift;};





  template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
  inline LazyStencil <T,C,MEMTYPE,D,E>
    operator()(const BoxData<T,C,MEMTYPE,D,E>&  a_src, T a_scale = 1) const;


  template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    inline LazyStencil <T,C,MEMTYPE,D,E>
    operator()(const BoxData<T,C,MEMTYPE,D,E>&  a_src, Box a_box, T a_scale = 1) const;




    inline void invert(int a_dir);


    inline void transpose(unsigned char a, unsigned char b);


    inline Box indexRange(Box a_domain) const;


    inline Box indexDomain(Box a_range) const;


    inline Box range(Box a_domain) const;


    inline Box domain(Box a_range) const;


    T diagonalValue() const;


    inline void print() const;

    // End of Stencil Operations Doxygen Module




    static Stencil<T> Derivative(int a_n, int a_dir = 0, int a_order = 2);


    static Stencil<T> Laplacian();
    #if DIM == 2
    static Stencil<T> Laplacian_9();
    #elif DIM == 3
    static Stencil<T> Laplacian_19();
    static Stencil<T> Laplacian_27();
    #endif

    static Stencil<T> LaplacianFace(int a_dir, int a_order = 2);


    static Stencil<T> CellToEdge(int a_dir, int a_order = 4);


    static Stencil<T> DiffCellToEdge(int a_dir, int a_order = 4);


    static Stencil<T> EdgeToCell(int a_dir, int a_order = 4);


    static Stencil<T> CellToEdgeL(int a_dir, int a_order = 5);


    static Stencil<T> CellToEdgeH(int a_dir, int a_order = 5);


    static Stencil<T> AvgDown(int a_refRatio = 2);


    static Stencil<T> AvgDownEdge(int a_dir, int a_refRatio = 2);


    static Stencil<T> FluxDivergence(int a_dir);

    // End of Stencil Operations Doxygen Module

  template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    void apply(const BoxData<T,C,MEMTYPE,D,E>&   a_src,
               BoxData<T,C,MEMTYPE,D,E>&         a_dest,
               const Box&                a_bx,
               bool                      a_initToZero = false,
               const T                   a_scale = 1) const;




    template<unsigned int C, unsigned char D, unsigned char E, MemType MEMTYPE>
    void protoApply(  const BoxData<T,C,MEMTYPE,D,E>&  a_src,
                      BoxData<T,C,MEMTYPE,D,E>&  a_dst,
                      const Box&               a_box,
                      bool               a_initToZero = false,
                      T                  a_scale = 1);


  template<unsigned int C, unsigned char D, unsigned char E>   
  void hostApply(  const BoxData<T,C,MemType::HOST,D,E>&  a_src,
                   BoxData<T,C,MemType::HOST,D,E>&  a_dst,
                   const Box&               a_box,
                   bool               a_initToZero = false,
                   T                  a_scale = 1) const;

 
   inline bool getIsClosed() {return m_isClosed;}

private:

    void augment(T a_coef, Point a_offset);

    std::vector<T>      m_coefs;        
    std::vector<Point>  m_offsets;      
    Point               m_srcRefratio;  
    Point               m_destRefratio; 
    Point               m_destShift;    
    Box                  m_span;         

    bool m_isClosed;
  //this copies to the device
  T*    d_coeff;
  Point* d_offset;
  Stencil<T>::coeff_holder c_coeff;

};




template <typename T>
inline Stencil<T> operator*(T a_coef, Shift a_shift)
{
    return Stencil<T>(a_shift, a_coef);
}


template <typename T>
inline Stencil<T> operator*(T a_coef, const Stencil<T> a_stencil)
{
    return a_stencil*a_coef;
}


template <typename T, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
inline BoxData<T,C,MEMTYPE,D,E>& operator|=(BoxData<T,C,MEMTYPE,D,E>& a_dest, LazyStencil<T,C,MEMTYPE,D,E>&& a_op);


template <class T, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
BoxData<T,C,MEMTYPE,D,E>& operator+=(BoxData<T,C,MEMTYPE,D,E>& a_dest, LazyStencil<T,C,MEMTYPE,D,E>&& a_op);


//=======================================================================================
// InterpStencil ||
//===============++


template <class T>
class InterpStencil
{
public:
    inline InterpStencil()
    {
        m_closed = false;
    }


    inline InterpStencil(int a_ratio)
    {
        define(Point::Ones(a_ratio));
    }


    inline InterpStencil(Point a_ratio)
    {
        define(a_ratio);
    }



    inline void define(Point a_ratio)
    {
        for (int ii = 0; ii < DIM; ii++)
        {
            PROTO_ASSERT(a_ratio[ii] > 0,
            "InterpStencil(Point ratio) invalid. All ratios must be 1 or greater");
        }
        m_k = Box(a_ratio - Point::Ones());
        m_s.resize(m_k.size());
        m_closed = false;
    }

//  inline void close()
//    {
//        for (auto iter = m_k.begin(); iter != m_k.end(); ++iter)
//        {
//            Stencil<T>& S = (*this)(*iter);
//            S.destRatio() = m_k.high() + Point::Ones();
//            S.destShift() = *iter;
//        }
//        m_closed = true;
//    }


    inline Stencil<T>& operator()(Point a_p)
    {
        PROTO_ASSERT(!m_closed,
        "Components of InterpStencil are read-only once closed. Use InterpStencil::get instead.");
        PROTO_ASSERT(m_k.contains(a_p),
        "InterpStencil::operator()(Point p) invalid.\
         p is not a member of the InterStencil kernel");

        return m_s[m_k.index(a_p)];
    }


    inline const Stencil<T>& get(Point a_p) const
    {
        PROTO_ASSERT(m_k.contains(a_p),
        "InterpStencil::operator()(Point p) invalid.\
         p is not a member of the InterStencil kernel");

        return m_s[m_k.index(a_p)];
    }


    inline Box span()
    {
        Box span;
        for (int ii = 0; ii < m_s.size(); ii++)
        {
            span = span & m_s[ii].span().low();
            span = span & m_s[ii].span().high();
        }
        return span;
    }

    inline Point spanPoint()
    {
        Box spanBox = span();
        int v[DIM];
        for (int ii = 0; ii < DIM; ii++)
        {
            v[ii] = std::max(std::abs(spanBox.low()[ii]), std::abs(spanBox.high()[ii]));
        }
        return Point(v);
    }


  template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    inline LazyStencil<T,C,MEMTYPE,D,E> operator()(const BoxData<T,C,MEMTYPE,D,E>& a_src,
                           T                       a_scale = 1) const;

  template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
    inline LazyStencil<T,C,MEMTYPE,D,E> operator()(const BoxData<T,C,MEMTYPE,D,E>& a_src,
                           Box                      a_box,
                           T                       a_scale = 1) const;

    static inline InterpStencil<T> PiecewiseConstant(Point a_ratio);

    static inline InterpStencil<T> PiecewiseLinear(Point a_ratio);

    static inline InterpStencil<T> Quadratic(int a_ratio);


    static inline InterpStencil<T> Build(int a_shiftMax, Box a_shiftKernel,
        int a_order, int a_refRatio);


    static inline InterpStencil<T> Build(
            std::vector<Point>& a_shifts,
            int a_maxOrder,
            int a_refRatio);

    static inline InterpStencil<T> Build(
            std::vector<Point>& a_shifts,
            const std::vector<Point>& a_orders,
            int a_refRatio);


    inline Box kernel() const {return m_k;}

    inline Point ratio() const {return m_r;}

    inline typename std::vector<Stencil<T>>::iterator begin(){return m_s.begin();}

    inline typename std::vector<Stencil<T>>::iterator end(){return m_s.end();}


    inline bool empty(){return (m_s.size() <= 0);}


    inline unsigned long long int size() const {return m_s.size();}
private:
    Point m_r;                      
    Box m_k;                         
    std::vector<Stencil<T>> m_s;    
    bool m_closed = false;          
};

#include "implem/Proto_StencilImplem.H"
#include "implem/Proto_StencilDefs.H"
} //end Proto namespace
#endif