grid_dist_amr.hpp

/*
 * grid_amr_dist.hpp
 *
 *  Created on: Sep 21, 2017
 *      Author: i-bird
 */

#ifndef AMR_GRID_AMR_DIST_HPP_
#define AMR_GRID_AMR_DIST_HPP_

#define OPENFPM_DATA_ENABLE_IO_MODULE

#include "Grid/grid_dist_id.hpp"
#include "Amr/grid_dist_amr_key_iterator.hpp"

#ifdef __NVCC__
#include "SparseGridGpu/SparseGridGpu.hpp"
#endif

#define AMR_IMPL_TRIVIAL 1
#define AMR_IMPL_PATCHES 2
#define AMR_IMPL_OPENVDB 3

template<typename Decomposition, typename garray>
class Decomposition_encap
{
    Decomposition & dec;
    garray & gd_array;

public:

    Decomposition_encap(Decomposition & dec, garray & gd_array)
    :dec(dec),gd_array(gd_array)
    {}

    Decomposition & internal_dec() const
    {
        return dec;
    }

    void decompose()
    {
        dec.decompose();

        for(size_t i = 0 ; i < gd_array.size() ; i++)
        {
            Ghost<Decomposition::dims,typename Decomposition::stype> gold = gd_array.get(i).getDecomposition().getGhost();
            gd_array.get(i).getDecomposition() = dec.duplicate(gold);
        }
    }

    void refine(size_t ts)
    {
        dec.refine();

        for(size_t i = 0 ; i < gd_array.size() ; i++)
        {
            Ghost<Decomposition::dims,typename Decomposition::stype> gold = gd_array.get(i).getDecomposition().getGhost();
            gd_array.get(i).getDecomposition() = dec.duplicate(gold);
        }
    }

    void redecompose(size_t ts)
    {
        dec.redecompose();

        for(size_t i = 0 ; i < gd_array.size() ; i++)
        {
            Ghost<Decomposition::dims,typename Decomposition::stype> gold = gd_array.get(i).getDecomposition().getGhost();
            gd_array.get(i).getDecomposition() = dec.duplicate(gold);
        }
    }

    auto getDistribution() -> decltype(dec.getDistribution())
    {
        return dec.getDistribution();
    }

    Decomposition_encap<Decomposition,garray> operator=(const Decomposition_encap<Decomposition,garray> & de) const
    {
        for(size_t i = 0 ; i < gd_array.size() ; i++)
        {gd_array.get(i).getDecomposition() = de.gd_array.get(i).getDecomposition();}

        return *this;
    }

    bool write(std::string output) const
    {
        return dec.write(output);
    }
};

template<unsigned int dim,
         typename St,
         typename T,
         unsigned int impl=AMR_IMPL_TRIVIAL ,
         typename Decomposition = CartDecomposition<dim,St>,
         typename Memory=HeapMemory ,
         typename device_grid=grid_cpu<dim,T> >
class grid_dist_amr
{

};


template<unsigned int dim,
         typename St,
         typename T,
         typename Decomposition,
         typename Memory,
         typename device_grid >
class grid_dist_amr<dim,St,T,AMR_IMPL_TRIVIAL,Decomposition,Memory,device_grid>
{
    Box<dim,St> domain;

    Ghost<dim,long int> g_int;

    periodicity<dim> bc;

    //
    openfpm::vector<grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>,
                                 HeapMemory,
                                 memory_traits_lin,
                                 openfpm::grow_policy_identity,STD_VECTOR> gd_array;

    typedef decltype(device_grid::type_of_subiterator()) device_sub_it;

    typedef decltype(device_grid::type_of_iterator()) device_it;

    openfpm::vector<grid_dist_iterator<dim,device_grid,device_sub_it,FREE>> git;

    openfpm::vector<grid_dist_iterator<dim,device_grid,device_it,FIXED>> git_g;

    openfpm::vector<grid_dist_iterator_sub<dim,device_grid>> git_sub;

    openfpm::vector<openfpm::vector<offset_mv<dim>>> mv_off;

    T bck;

    void initialize_other(size_t n_lvl, size_t (& g_sz_lvl)[dim])
    {
        for (size_t i = 0; i < n_lvl - 1 ; i++)
        {
            for (size_t j = 0 ; j < dim ; j++)
            {
                if (bc.bc[j] == NON_PERIODIC)
                {g_sz_lvl[j] = (g_sz_lvl[j]-1)*2 + 1;}
                else
                {g_sz_lvl[j] = g_sz_lvl[j]*2;}
            }

            gd_array.add(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>(gd_array.get(0).getDecomposition(),g_sz_lvl,g_int));
            gd_array.last().setBackgroundValue(bck);

            gd_array.last().getDecomposition().free_geo_cell();
            gd_array.last().getDecomposition().getDistribution().destroy_internal_graph();
            gd_array.last().getDecomposition().free_fines();
        }

        recalculate_mvoff();
    }

public:


    grid_dist_amr(const Box<dim,St> & domain, const Ghost<dim,long int> & g)
    :domain(domain),g_int(g)
    {
        // set boundary consitions to non periodic

        for (size_t i = 0; i < dim ; i++)
        {bc.bc[i] = NON_PERIODIC;}
    }

    grid_dist_amr(const Box<dim,St> & domain, const Ghost<dim,long int> & g, periodicity<dim> & bc)
    :domain(domain),g_int(g),bc(bc)
    {
    }

    void initLevels(const Decomposition & dec, size_t n_lvl,const size_t (& g_sz)[dim])
    {
        size_t g_sz_lvl[dim];

        for (size_t i = 0; i < dim ; i++)
        {g_sz_lvl[i] = g_sz[i];}

        // Add the coarse level
        gd_array.add(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>(dec,g_sz,g_int));
        gd_array.last().setBackgroundValue(bck);

        initialize_other(n_lvl,g_sz_lvl);
    }

    template<typename TT> void initLevels(const Decomposition_encap<Decomposition,TT> & dec, size_t n_lvl,const size_t (& g_sz)[dim])
    {
        initLevels(dec.internal_dec(),n_lvl,g_sz);
    }

    void recalculate_mvoff()
    {
        // Here we calculate the offset to move one level up and one level down
        // in global coordinated moving one level up is multiply the coordinates by 2
        // and moving one level down is dividing by 2. In local coordinates is the same
        // with the exception that because of the decomposition you can have an offset
        // look at the picture below
        //
        // (-1)  (0)
        //  * |   *     *    coarse level
        //  * |*  *  *  *    finer level
        //    |(0)(1)
        //
        //  Line of the decomposition
        //
        // The coarse level point 0 in local coordinates converted to the finer level is not
        // just 2*0 = 0 but is 2*(0) + 1 so a formula like 2*x+offset is required. here we calculate
        // these offset. In the case of moving from finer to coarse is the same the formula is
        // Integer_round(x+1)/2 - 1
        //
        mv_off.resize(gd_array.size());

        for (size_t i = 1 ; i < gd_array.size() ; i++)
        {
            auto & g_box_c = gd_array.get(i-1).getLocalGridsInfo();
            auto & g_box_f = gd_array.get(i).getLocalGridsInfo();

#ifdef SE_CLASS1

            if (g_box_c.size() != g_box_f.size())
            {
                std::cerr << __FILE__ << ":" << __LINE__ << " error it seem that the AMR construction between level " <<
                        i << " and " << i-1 << " is inconsistent" << std::endl;
            }

#endif

            mv_off.get(i-1).resize(g_box_f.size());
            mv_off.get(i).resize(g_box_f.size());

            for (size_t j = 0 ; j < g_box_f.size() ; j++)
            {
                for (size_t s = 0 ; s < dim ; s++)
                {
                    size_t d_orig_c = g_box_c.get(j).origin.get(s);
                    size_t d_orig_f = g_box_f.get(j).origin.get(s);

                    mv_off.get(i-1).get(j).dw.get(s) = d_orig_c*2 - d_orig_f;
                    mv_off.get(i).get(j).up.get(s) = d_orig_c*2 - d_orig_f;
                }
            }
        }
    }

    void initLevels(size_t n_lvl,const size_t (& g_sz)[dim], size_t opt = 0)
    {
        size_t g_sz_lvl[dim];

        for (size_t i = 0; i < dim ; i++)
        {g_sz_lvl[i] = g_sz[i];}

        // Add the coarse level
        gd_array.add(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>(g_sz,domain,g_int,bc,opt));

        initialize_other(n_lvl,g_sz_lvl);
    }

    template <typename Model>inline void addComputationCosts(Model md=Model(), size_t ts = 1)
    {
        gd_array.get(0).addComputationCosts(md,ts);
    }

    Decomposition_encap<Decomposition,decltype(gd_array)> getDecomposition()
    {
        Decomposition_encap<Decomposition,decltype(gd_array)> tmp(gd_array.get(0).getDecomposition(),gd_array);

        return tmp;
    }

    grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & getLevel(size_t lvl)
    {
        return gd_array.get(lvl);
    }


    grid_dist_amr_key_iterator<dim,device_grid,
                               decltype(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>::type_of_subiterator()),
                               decltype(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>::type_of_subiterator()) >
    getDomainIteratorCells()
    {
        git_sub.clear();

        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            grid_key_dx<dim> start;
            grid_key_dx<dim> stop;

            for (size_t j = 0 ; j < dim ; j++)
            {
                start.set_d(j,0);
                if (bc.bc[j] == NON_PERIODIC)
                {stop.set_d(j,getGridInfoVoid(i).size(j) - 2);}
                else
                {stop.set_d(j,getGridInfoVoid(i).size(j) - 1);}
            }

            git_sub.add(gd_array.get(i).getSubDomainIterator(start,stop));
        }

        return grid_dist_amr_key_iterator<dim,device_grid,
                                            decltype(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>::type_of_subiterator()),
                                            decltype(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>::type_of_subiterator())>(git_sub);
    }

    grid_dist_iterator_sub<dim,device_grid> getDomainIteratorCells(size_t lvl)
    {
        grid_key_dx<dim> start;
        grid_key_dx<dim> stop;

        for (size_t j = 0 ; j < dim ; j++)
        {
            start.set_d(j,0);
            if (bc.bc[j] == NON_PERIODIC)
            {stop.set_d(j,getGridInfoVoid(lvl).size(j) - 2);}
            else
            {stop.set_d(j,getGridInfoVoid(lvl).size(j) - 1);}
        }

        return gd_array.get(lvl).getSubDomainIterator(start,stop);
    }

    auto getGridGhostIterator(size_t lvl) -> decltype(gd_array.get(lvl).getGridGhostIterator(grid_key_dx<dim>(),grid_key_dx<dim>()))
    {
        grid_key_dx<dim> key_start;
        grid_key_dx<dim> key_stop;

        for (size_t i = 0 ; i < dim ; i++)
        {
            key_start.set_d(i,g_int.getLow(i));
            key_stop.set_d(i,g_int.getHigh(i) + getGridInfoVoid(lvl).size(i) -1);
        }

        return gd_array.get(lvl).getGridGhostIterator(key_start,key_stop);
    }

    auto getGridIterator(size_t lvl) -> decltype(gd_array.get(lvl).getGridIterator())
    {
        return gd_array.get(lvl).getGridIterator();
    }

    auto getGridIterator(size_t lvl, grid_key_dx<dim> & start, grid_key_dx<dim> & stop) -> decltype(gd_array.get(lvl).getGridIterator(start,stop))
    {
        return gd_array.get(lvl).getGridIterator(start,stop);
    }

#ifdef __NVCC__

    auto getGridIteratorGPU(size_t lvl) -> decltype(gd_array.get(lvl).getGridIteratorGPU())
    {
        return gd_array.get(lvl).getGridIteratorGPU();
    }

#endif

    auto getGridIteratorCells(size_t lvl) -> decltype(gd_array.get(lvl).getGridIterator())
    {
        grid_key_dx<dim> start;
        grid_key_dx<dim> stop;

        for (size_t j = 0 ; j < dim ; j++)
        {
            start.set_d(j,0);
            if (bc.bc[j] == NON_PERIODIC)
            {stop.set_d(j,getGridInfoVoid(lvl).size(j) - 2);}
            else
            {stop.set_d(j,getGridInfoVoid(lvl).size(j) - 1);}
        }

        return gd_array.get(lvl).getGridIterator(start,stop);
    }


    grid_dist_iterator<dim,device_grid,decltype(device_grid::type_of_subiterator()),FREE>
    getDomainIterator(size_t lvl) const
    {
        return gd_array.get(lvl).getDomainIterator();
    }

    grid_dist_iterator<dim,device_grid,
    decltype(device_grid::type_of_iterator()),
    FIXED>
    getDomainGhostIterator(size_t lvl) const
    {
            return gd_array.get(lvl).getDomainGhostIterator();
    }

    grid_dist_amr_key_iterator<dim,device_grid, decltype(device_grid::type_of_subiterator())>
    getDomainIterator()
    {
        git.clear();

        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            git.add(gd_array.get(i).getDomainIterator());
        }

        return grid_dist_amr_key_iterator<dim,device_grid,decltype(device_grid::type_of_subiterator())>(git);
    }

    grid_dist_amr_key_iterator<dim,device_grid, decltype(device_grid::type_of_iterator()),
                               grid_dist_iterator<dim,device_grid,decltype(device_grid::type_of_iterator()),FIXED>>
    getDomainGhostIterator()
    {
            git_g.clear();

            for (size_t i = 0 ; i < gd_array.size() ; i++)
            {
                    git_g.add(gd_array.get(i).getDomainGhostIterator());
            }

            return grid_dist_amr_key_iterator<dim,device_grid,decltype(device_grid::type_of_iterator()),
                                              grid_dist_iterator<dim,device_grid,decltype(device_grid::type_of_iterator()),FIXED>>(git_g);
    }

    template <unsigned int p>inline auto get(const grid_dist_amr_key<dim> & v1) const -> decltype(gd_array.get(v1.getLvl()).template get<p>(v1.getKey()))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(v1.getLvl()).template get<p>(v1.getKey());
    }

    template <unsigned int p>inline auto get(const grid_dist_amr_key<dim> & v1) -> decltype(gd_array.get(v1.getLvl()).template get<p>(v1.getKey()))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(v1.getLvl()).template get<p>(v1.getKey());
    }


    template <unsigned int p>inline auto get(size_t lvl, const grid_dist_key_dx<dim> & v1) const -> decltype(gd_array.get(lvl).template get<p>(v1))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(lvl).template get<p>(v1);
    }

    template <unsigned int p>inline auto get(size_t lvl, const grid_dist_key_dx<dim> & v1) -> decltype(gd_array.get(lvl).template get<p>(v1))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(lvl).template get<p>(v1);
    }



    template <unsigned int p>
    inline auto insert(const grid_dist_amr_key<dim> & v1)
    -> decltype(gd_array.get(v1.getLvl()).template insert<p>(v1.getKey()))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(v1.getLvl()).template insert<p>(v1.getKey());
    }



    template <unsigned int p>inline auto insert(size_t lvl, const grid_dist_key_dx<dim> & v1)
    -> decltype(gd_array.get(lvl).template insert<p>(v1))
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(lvl).template insert<p>(v1);
    }


    grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & getDistGrid(size_t lvl)
    {
        return gd_array.get(lvl);
    }


    inline void remove(const grid_dist_amr_key<dim> & v1)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(v1.getLvl()).remove(v1.getKey());
    }

    void remove(size_t lvl, const grid_dist_key_dx<dim> & v1)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return gd_array.get(lvl).remove(v1);
    }

    void construct_level_connections()
    {
        for (int lvl = 0 ; lvl < gd_array.size() ; lvl++)
        {
            if (lvl == 0)
            {
                gd_array.get(lvl).construct_link_dw(gd_array.get(lvl+1),mv_off.get(lvl));
            }
            else if (lvl == gd_array.size() - 1)
            {gd_array.get(lvl).construct_link_up(gd_array.get(lvl-1),mv_off.get(lvl));}
            else
            {
                gd_array.get(lvl).construct_link_dw(gd_array.get(lvl+1),mv_off.get(lvl));
                gd_array.get(lvl).construct_link_up(gd_array.get(lvl-1),mv_off.get(lvl));
            }
        }
    }

    template<typename stencil_type>
    void tagBoundaries()
    {
        for (int lvl = 0 ; lvl < gd_array.size() ; lvl++)
        {
            gd_array.get(lvl).template tagBoundaries<stencil_type>();
        }
    }


    template<int... prp> void ghost_get(size_t opt = 0)
    {
        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            gd_array.get(i).template ghost_get<prp...>(opt);
        }
    }

    void map(size_t opt = 0)
    {
        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            gd_array.get(i).map();
        }

        recalculate_mvoff();
    }

    template<template<typename,typename> class op,int... prp> void ghost_put()
    {
        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            gd_array.get(i).template ghost_put<op,prp...>();
        }
    }

    size_t size_inserted(size_t lvl)
    {
        return gd_array.get(lvl).size_local_inserted();
    }

    void setBackgroundValue(T & bv)
    {
        for (size_t i = 0 ; i < getNLvl() ; i++)
        {gd_array.get(i).setBackgroundValue(bv);}

        meta_copy<T>::meta_copy_(bv,bck);
    }

    void clear()
    {
        for (size_t i = 0 ; i < getNLvl() ; i++)
        {gd_array.get(i).clear();}
    }

    const grid_sm<dim,void> & getGridInfoVoid(size_t lvl) const
    {
        return gd_array.get(lvl).getGridInfoVoid();
    }

    size_t getNLvl()
    {
        return gd_array.size();
    }

    void moveLvlDw(grid_dist_amr_key<dim> & key)
    {
#ifdef SE_CLASS1

        if (key.getLvl() >= getNLvl() - 1)
        {std::cerr << __FILE__ << ":" << __LINE__ << " error: we are already at the last level, we cannot go one level down" << std::endl;}

#endif

        auto & key_ref = key.getKeyRef().getKeyRef();
        size_t lvl = key.getLvl();

        for (size_t i = 0 ; i < dim ; i++)
        {
            key_ref.set_d(i,(key_ref.get(i) << 1) + mv_off.get(key.getLvl()).get(key.getKeyRef().getSub()).dw.get(i) );
        }

        key.setLvl(lvl+1);
    }

    grid_dist_key_dx<dim> moveDw(int lvl, const grid_dist_key_dx<dim> & key)
    {
#ifdef SE_CLASS1

        if (lvl >= getNLvl() - 1)
        {std::cerr << __FILE__ << ":" << __LINE__ << " error: we are already at the last level, we cannot go one level down" << std::endl;}

#endif

        grid_dist_key_dx<dim> out;

        for (size_t i = 0 ; i < dim ; i++)
        {
            out.getKeyRef().set_d(i,(key.getKeyRef().get(i) << 1) + mv_off.get(lvl).get(key.getSub()).dw.get(i) );
        }

        out.setSub(key.getSub());

        return out;
    }

    inline grid_dist_amr_key<dim> getAMRKey(size_t lvl, grid_dist_key_dx<dim> key)
    {
        return grid_dist_amr_key<dim>(lvl,key);
    }

    void moveLvlUp(grid_dist_amr_key<dim> & key)
    {
#ifdef SE_CLASS1

        if (key.getLvl() == 0)
        {std::cerr << __FILE__ << ":" << __LINE__ << " error: we are already at the top level, we cannot go one level up" << std::endl;}

#endif

        auto & key_ref = key.getKeyRef().getKeyRef();
        size_t lvl = key.getLvl();

        for (size_t i = 0 ; i < dim ; i++)
        {
            key_ref.set_d(i,(key_ref.get(i) - mv_off.get(key.getLvl()).get(key.getKeyRef().getSub()).up.get(i)) >> 1);
        }

        key.setLvl(lvl-1);
    }

    grid_dist_key_dx<dim> moveUp(int lvl, const grid_dist_key_dx<dim> & key)
    {
#ifdef SE_CLASS1

        if (lvl == 0)
        {std::cerr << __FILE__ << ":" << __LINE__ << " error: we are already at the top level, we cannot go one level up" << std::endl;}

#endif

        grid_dist_key_dx<dim> out;

        for (size_t i = 0 ; i < dim ; i++)
        {
            out.getKeyRef().set_d(i,(key.getKeyRef().get(i) - mv_off.get(lvl).get(key.getSub()).up.get(i)) >> 1);
        }

        out.setSub(key.getSub());

        return out;
    }

    grid_key_dx<dim> getGKey(const grid_dist_amr_key<dim> & v1)
    {
        return gd_array.get(v1.getLvl()).getGKey(v1.getKey());
    }

    grid_key_dx<dim> getGKey(int lvl, const grid_dist_key_dx<dim> & v1)
    {
        return gd_array.get(lvl).getGKey(v1);
    }

    Point<dim,St> getPos(const grid_dist_amr_key<dim> & v1)
    {
        return gd_array.get(v1.getLvl()).getPos(v1.getKey());
    }

    Point<dim,St> getPos(int lvl, const grid_dist_key_dx<dim> & v1)
    {
        return gd_array.get(lvl).getPos(v1);
    }

    Point<dim,St> getSpacing(size_t lvl)
    {
        return gd_array.get(lvl).getSpacing();
    }

    /* \brief Check if a point exist
     *
     * \param v1 point to checl
     * 
     */
    bool existPoint(const grid_dist_amr_key<dim> & v1)
    {
        return gd_array.get(v1.getLvl()).existPoint(v1.getKey());
    }

    /* \brief Check if a point exist
     *
     * \param v1 point to checl
     * 
     */
    bool existPoint(int lvl, const grid_dist_key_dx<dim> & v1)
    {
        return gd_array.get(lvl).existPoint(v1);
    }


    bool write(std::string output, size_t opt = VTK_WRITER | FORMAT_ASCII )
    {
        bool ret = true;

        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            ret &= gd_array.get(i).write(output + "_" + std::to_string(i),opt);
        }

        return ret;
    }

#ifdef __NVCC__

    template<unsigned int ... prp> void deviceToHost()
    {
        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            gd_array.get(i).template deviceToHost<prp ...>();
        }
    }

    template<unsigned int ... prp> void hostToDevice()
    {
        for (size_t i = 0 ; i < gd_array.size() ; i++)
        {
            gd_array.get(i).template hostToDevice<prp ...>();
        }
    }

#endif
};

template<unsigned int dim, typename St, typename T>
using sgrid_dist_amr = grid_dist_amr<dim,St,T,AMR_IMPL_TRIVIAL,CartDecomposition<dim,St>,HeapMemory,sgrid_cpu<dim,T,HeapMemory>>;

#ifdef __NVCC__

template<unsigned int dim, typename St, typename T, unsigned int blockEdgeSize = 8>
using sgrid_dist_amr_gpu = grid_dist_amr<dim,St,T,AMR_IMPL_TRIVIAL,CartDecomposition<dim,St,CudaMemory,memory_traits_inte>,CudaMemory,SparseGridGpu<dim,T,blockEdgeSize,IntPow<blockEdgeSize,dim>::value >>;

#endif

#endif /* AMR_GRID_AMR_DIST_HPP_ */