From 6133c34ca5c5dd51a308a0fabff008319535b920 Mon Sep 17 00:00:00 2001 From: Pietro Incardona Date: Wed, 26 Jun 2019 18:26:19 +0200 Subject: [PATCH] grid_dist_id work with SparseGridGPU --- openfpm_data | 2 +- src/Amr/grid_dist_amr.hpp | 895 ++++++++++++++++ src/Amr/grid_dist_amr_key.hpp | 98 ++ src/Amr/grid_dist_amr_key_iterator.hpp | 136 +++ src/Amr/grid_dist_amr_unit_tests.cpp | 993 ++++++++++++++++++ src/Amr/tests/amr_base_unit_tests.cpp | 171 +++ src/CMakeLists.txt | 4 +- src/Decomposition/Domain_icells_cart.hpp | 8 +- .../Iterators/grid_dist_id_iterator_dec.hpp | 50 + src/Grid/cuda/grid_dist_id_kernels.cuh | 20 + src/Grid/grid_dist_id.hpp | 109 +- src/Grid/tests/grid_dist_id_dlb_unit_test.cpp | 242 +++++ .../tests/sgrid_dist_id_gpu_unit_tests.cu | 202 ++++ src/Grid/tests/sgrid_dist_id_unit_tests.cpp | 250 +++++ src/SubdomainGraphNodes.hpp | 2 +- 15 files changed, 3171 insertions(+), 11 deletions(-) create mode 100644 src/Amr/grid_dist_amr.hpp create mode 100644 src/Amr/grid_dist_amr_key.hpp create mode 100644 src/Amr/grid_dist_amr_key_iterator.hpp create mode 100644 src/Amr/grid_dist_amr_unit_tests.cpp create mode 100644 src/Amr/tests/amr_base_unit_tests.cpp create mode 100644 src/Grid/cuda/grid_dist_id_kernels.cuh create mode 100644 src/Grid/tests/grid_dist_id_dlb_unit_test.cpp create mode 100644 src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu create mode 100644 src/Grid/tests/sgrid_dist_id_unit_tests.cpp diff --git a/openfpm_data b/openfpm_data index f982d0ac..b0330228 160000 --- a/openfpm_data +++ b/openfpm_data @@ -1 +1 @@ -Subproject commit f982d0ac3f3ae0eefbd1ca9e39a700b8898f8ee9 +Subproject commit b0330228e8503af944dcc5a05fa4fb1bc3e3b594 diff --git a/src/Amr/grid_dist_amr.hpp b/src/Amr/grid_dist_amr.hpp new file mode 100644 index 00000000..76e2a2ab --- /dev/null +++ b/src/Amr/grid_dist_amr.hpp @@ -0,0 +1,895 @@ +/* + * grid_amr_dist.hpp + * + * Created on: Sep 21, 2017 + * Author: i-bird + */ + +#ifndef AMR_GRID_AMR_DIST_HPP_ +#define AMR_GRID_AMR_DIST_HPP_ + +#include "Grid/grid_dist_id.hpp" +#include "Amr/grid_dist_amr_key_iterator.hpp" + +#define AMR_IMPL_TRIVIAL 1 +#define AMR_IMPL_PATCHES 2 +#define AMR_IMPL_OPENVDB 3 + +template +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; + } + + /*! \brief Start decomposition + * + */ + void decompose() + { + dec.decompose(); + + for(size_t i = 0 ; i < gd_array.size() ; i++) + { + Ghost gold = gd_array.get(i).getDecomposition().getGhost(); + gd_array.get(i).getDecomposition() = dec.duplicate(gold); + } + } + + /*! \brief Refine the decomposition, available only for ParMetis distribution, for Metis it is a null call + * + * \param ts number of time step from the previous load balancing + * + */ + void refine(size_t ts) + { + dec.refine(); + + for(size_t i = 0 ; i < gd_array.size() ; i++) + { + Ghost gold = gd_array.get(i).getDecomposition().getGhost(); + gd_array.get(i).getDecomposition() = dec.duplicate(gold); + } + } + + /*! \brief Refine the decomposition, available only for ParMetis distribution, for Metis it is a null call + * + * \param ts number of time step from the previous load balancing + * + */ + void redecompose(size_t ts) + { + dec.redecompose(); + + for(size_t i = 0 ; i < gd_array.size() ; i++) + { + Ghost 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 operator=(const Decomposition_encap & 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,typename Memory=HeapMemory , typename device_grid=grid_cpu > +class grid_dist_amr +{ + +}; + +/*! \brief It contain the offset necessary to move to coarser and finer level grids + * + */ +template +struct offset_mv +{ + //! offset to move up on an upper grid (coarse) + Point up; + + //! offset to move on the lower grid (finer) + Point dw; +}; + +/*! \brief AMR Adaptive Multi Resolution Grid + * + * \tparam dim Dimensionality + * \tparam St type of space + * \tparam T what each point of the grid store + * \tparam Decomposition type of decomposition + * + */ +template +class grid_dist_amr +{ + //! Simulation domain + Box domain; + + //! Ghost integer + Ghost g_int; + + //! Boundary conditions of the structure + periodicity bc; + + //! array of grids + // + openfpm::vector, + HeapMemory, + typename memory_traits_lin>::type, + memory_traits_lin, + openfpm::grow_policy_identity,STD_VECTOR> gd_array; + + //! Type of structure sub-grid iterator + typedef decltype(device_grid::type_of_subiterator()) device_sub_it; + + //! Type of structure for the grid iterator + typedef decltype(device_grid::type_of_iterator()) device_it; + + //! Domain iterator for each distributed grid + openfpm::vector> git; + + //! Domain and ghost iterator for each distributed grid + openfpm::vector> git_g; + + //! Iterator for each distributed grid + openfpm::vector> git_sub; + + //! Moving offsets + openfpm::vector>> mv_off; + + //! background level + T bck; + + /*! \brief Initialize the others levels + * + * \param n_grid_dist_idlvl number of levels + * \param g_sz_lvl grid size on each level + * + */ + 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(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: + + + /*! \brief Constructor + * + * \param domain Simulation domain + * \param g ghost extension + * + */ + grid_dist_amr(const Box & domain, const Ghost & 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;} + } + + /*! \brief Constructor + * + * \param domain Simulation domain + * \param g ghost extension + * \param bc boundary conditions + * + */ + grid_dist_amr(const Box & domain, const Ghost & g, periodicity & bc) + :domain(domain),g_int(g),bc(bc) + { + } + + /*! \brief Initialize the amr grid + * + * \param dec Decomposition (this parameter is useful in case we want to constrain the AMR to an external decomposition) + * \param n_lvl maximum number of levels (0 mean no additional levels) + * \param g_sz coarsest grid size on each direction + * + */ + 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(dec,g_sz,g_int)); + gd_array.last().setBackgroundValue(bck); + + initialize_other(n_lvl,g_sz_lvl); + } + + /*! \brief Initialize the amr grid + * + * \param dec Decomposition (this parameter is useful in case we want to constrain the AMR to an external decomposition) + * \param n_lvl maximum number of levels (0 mean no additional levels) + * \param g_sz coarsest grid size on each direction + * + */ + template void initLevels(const Decomposition_encap & dec, size_t n_lvl,const size_t (& g_sz)[dim]) + { + initLevels(dec.internal_dec(),n_lvl,g_sz); + } + + /*! \brief Recalculate the offset array for the moveLvlUp and moveLvlDw + * + * + * + */ + 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; + } + } + } + } + + /*! \brief Initialize the amr grid + * + * \param n_lvl maximum number of levels (0 mean no additional levels) + * \param g_sz coarsest grid size on each direction + * \param opt options + * + */ + 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(g_sz,domain,g_int,bc,opt)); + + initialize_other(n_lvl,g_sz_lvl); + } + + /*! \brief Add the computation cost on the decomposition using a resolution function + * + * + * \param md Model to use + * \param ts It is an optional parameter approximately should be the number of ghost get between two + * rebalancing at first decomposition this number can be ignored (default = 1) because not used + * + */ + template inline void addComputationCosts(Model md=Model(), size_t ts = 1) + { + gd_array.get(0).addComputationCosts(md,ts); + } + + /*! \brief Get the object that store the information about the decomposition + * + * \return the decomposition object + * + */ + Decomposition_encap getDecomposition() + { + Decomposition_encap tmp(gd_array.get(0).getDecomposition(),gd_array); + + return tmp; + } + + /*! \brief Get the underlying grid level + * + * \param lvl level + * + * \return the grid level + * + */ + grid_dist_id & getLevel(size_t lvl) + { + return gd_array.get(lvl); + } + + + grid_dist_amr_key_iterator::type_of_subiterator()), + decltype(grid_dist_id::type_of_subiterator()) > + getDomainIteratorCells() + { + git_sub.clear(); + + for (size_t i = 0 ; i < gd_array.size() ; i++) + { + grid_key_dx start; + grid_key_dx 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::type_of_subiterator()), + decltype(grid_dist_id::type_of_subiterator())>(git_sub); + } + + grid_dist_iterator_sub getDomainIteratorCells(size_t lvl) + { + grid_key_dx start; + grid_key_dx 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); + } + + /*! \brief Get an iterator to the grid + * + * \return an iterator to the grid + * + */ + auto getGridGhostIterator(size_t lvl) -> decltype(gd_array.get(lvl).getGridGhostIterator(grid_key_dx(),grid_key_dx())) + { + grid_key_dx key_start; + grid_key_dx 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); + } + + /*! \brief Get an iterator to the grid + * + * \return an iterator to the grid + * + */ + auto getGridIterator(size_t lvl) -> decltype(gd_array.get(lvl).getGridIterator()) + { + return gd_array.get(lvl).getGridIterator(); + } + + /*! \brief Get an iterator to the grid + * + * \return an iterator to the grid + * + */ + auto getGridIteratorCells(size_t lvl) -> decltype(gd_array.get(lvl).getGridIterator()) + { + grid_key_dx start; + grid_key_dx 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); + } + + + /*! \brief return an iterator over the level lvl + * + * \param lvl level + * + * \return an iterator over the level lvl selected + * + */ + grid_dist_iterator + getDomainIterator(size_t lvl) const + { + return gd_array.get(lvl).getDomainIterator(); + } + + /*! \brief return an iterator over the level lvl + * + * \param lvl level + * + * \return an iterator over the level lvl selected + * + */ + grid_dist_iterator + getDomainGhostIterator(size_t lvl) const + { + return gd_array.get(lvl).getDomainGhostIterator(); + } + + /*! \brief Get domain iterator + * + * \return an iterator over all the grid levels + * + */ + grid_dist_amr_key_iterator + 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(git); + } + + /*! \brief Get domain iterator + * + * \return an iterator over all the grid levels + * + */ + grid_dist_amr_key_iterator> + 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>(git_g); + } + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto get(const grid_dist_amr_key & v1) const -> typename std::add_lvalue_reference(v1.getKey()))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(v1.getLvl()).template get

(v1.getKey()); + } + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto get(const grid_dist_amr_key & v1) -> typename std::add_lvalue_reference(v1.getKey()))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(v1.getLvl()).template get

(v1.getKey()); + } + + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto get(size_t lvl, const grid_dist_key_dx & v1) const -> typename std::add_lvalue_reference(v1))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(lvl).template get

(v1); + } + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto get(size_t lvl, const grid_dist_key_dx & v1) -> typename std::add_lvalue_reference(v1))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(lvl).template get

(v1); + } + + //////////////////// Insert functions + + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template + inline auto insert(const grid_dist_amr_key & v1) + -> typename std::add_lvalue_reference(v1.getKey()))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(v1.getLvl()).template insert

(v1.getKey()); + } + + + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto insert(size_t lvl, const grid_dist_key_dx & v1) + -> typename std::add_lvalue_reference(v1))>::type + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(lvl).template insert

(v1); + } + + ////////////////////////////////////// + + /*! \brief Get the internal distributed grid + * + * \param lvl level + * + * \return the internal distributed grid + * + */ + grid_dist_id & getDistGrid(size_t lvl) + { + return gd_array.get(lvl); + } + + //////////////////// Remove functions + + /*! \brief Remove a grid point (this function make sense only in case of + * sparse grid) + * + * \param v1 grid_key that identify the element in the AMR grid to eleminate + * + */ + inline void remove(const grid_dist_amr_key & v1) + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(v1.getLvl()).remove(v1.getKey()); + } + + /*! \brief Remove a grid point (this function make sense only in case of + * sparse grid) + * + * \param v1 grid_key that identify the element in the AMR grid to eleminate + * + */ + void remove(size_t lvl, const grid_dist_key_dx & v1) + { +#ifdef SE_CLASS2 + check_valid(this,8); +#endif + return gd_array.get(lvl).remove(v1); + } + + ////////////////////////////////////// + + /*! \brief It synchronize the ghost parts + * + * \tparam prp... Properties to synchronize + * + */ + template void ghost_get() + { + for (size_t i = 0 ; i < gd_array.size() ; i++) + { + gd_array.get(i).template ghost_get(); + } + } + + /*! \brief It move all the grid parts that do not belong to the local processor to the respective processor + * + */ + void map(size_t opt = 0) + { + for (size_t i = 0 ; i < gd_array.size() ; i++) + { + gd_array.get(i).map(); + } + + recalculate_mvoff(); + } + + /*! \brief Apply the ghost put + * + * \tparam prp... Properties to apply ghost put + * + */ + template class op,int... prp> void ghost_put() + { + for (size_t i = 0 ; i < gd_array.size() ; i++) + { + gd_array.get(i).template ghost_put(); + } + } + + /*! \brief Return the number of inserted points on a particular level + * + * \return the number of inserted points + * + */ + size_t size_inserted(size_t lvl) + { + return gd_array.get(lvl).size_local_inserted(); + } + + /*! \brief set the background value + * + * You can use this function make sense in case of sparse in case of dense + * it does nothing + * + */ + void setBackgroundValue(T & bv) + { + for (size_t i = 0 ; i < getNLvl() ; i++) + {gd_array.get(i).setBackgroundValue(bv);} + + meta_copy::meta_copy_(bv,bck); + } + + /*! \brief delete all the points in the grid + * + * In case of sparse grid in delete all the inserted points, in case + * of dense it does nothing + * + */ + void clear() + { + for (size_t i = 0 ; i < getNLvl() ; i++) + {gd_array.get(i).clear();} + } + + /*! \brief Get an object containing the grid informations for a specific level + * + * \param lvl level + * + * \return an information object about this grid + * + */ + const grid_sm & getGridInfoVoid(size_t lvl) const + { + return gd_array.get(lvl).getGridInfoVoid(); + } + + /*! \brief Return the maximum number of levels in the AMR struct + * + * \return the number of levels + * + */ + size_t getNLvl() + { + return gd_array.size(); + } + + /*! \brief Move down (to finer level) the key + * + * \param key multi-resolution AMR key + * + */ + void moveLvlDw(grid_dist_amr_key & 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); + } + + /*! \brief From a distributed key it return a AMR key that contain also the grid level + * + * \param lvl level + * \param key distributed key + * + */ + inline grid_dist_amr_key getAMRKey(size_t lvl, grid_dist_key_dx key) + { + return grid_dist_amr_key(lvl,key); + } + + /*! \brief Move up (to coarser level) the key + * + * \param key multi-resolution AMR key + * + */ + void moveLvlUp(grid_dist_amr_key & 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); + } + + /*! \brief Get the position on the grid in global coordinates + * + * \param v1 amr key + * + * \return the position in global coordinates + * + */ + grid_key_dx getGKey(const grid_dist_amr_key & v1) + { + return gd_array.get(v1.getLvl()).getGKey(v1.getKey()); + } + + /*! \brief return the spacing for the grid in the level lvl + * + * \param lvl level + * + * \return return the spacing + * + */ + Point getSpacing(size_t lvl) + { + return gd_array.get(lvl).getSpacing(); + } + + /*! \brief Write on vtk file + * + * \param output filename output + * + */ + 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; + } +}; + +template +using sgrid_dist_amr = grid_dist_amr,HeapMemory,sgrid_cpu>; + +#endif /* AMR_GRID_AMR_DIST_HPP_ */ diff --git a/src/Amr/grid_dist_amr_key.hpp b/src/Amr/grid_dist_amr_key.hpp new file mode 100644 index 00000000..48b76c30 --- /dev/null +++ b/src/Amr/grid_dist_amr_key.hpp @@ -0,0 +1,98 @@ +/* + * grid_dist_amr_key.hpp + * + * Created on: Sep 23, 2017 + * Author: i-bird + */ + +#ifndef SRC_AMR_GRID_DIST_AMR_KEY_HPP_ +#define SRC_AMR_GRID_DIST_AMR_KEY_HPP_ + +/*! \brief Amr grid distributed key + * + * \tparam dim dimensionality + * + */ +template +class grid_dist_amr_key +{ + //! actual level + size_t lvl; + + //! actual position in the distributed grid + grid_dist_key_dx key; + + +public: + + /*! \constructor + * + * \param lvl level + * \param key distributed grid key + * \param offsets to move between levels + * + */ + inline grid_dist_amr_key(size_t lvl, + grid_dist_key_dx key) + :lvl(lvl),key(key) + {} + + /*! \brief Return the grid key + * + * \return the distributed key + * + */ + inline const grid_dist_key_dx & getKey() const + { + return key; + } + + /*! \brief Return the grid key (as reference) + * + * \return the distributed key + * + */ + inline grid_dist_key_dx & getKeyRef() + { + return key; + } + + + /*! \brief Return the level + * + * \return the level + * + */ + inline size_t getLvl() const + { + return lvl; + } + + /*! \brief Return the level + * + * \param lvl level to set + * + */ + inline void setLvl(size_t lvl) + { + this->lvl = lvl; + } + + /*! \brief Create a new key moving the old one + * + * \param s dimension id + * \param s number of steps + * + * \return new key + * + */ + inline grid_dist_amr_key moveSpace(size_t d,size_t s) + { + return grid_dist_amr_key(lvl,key.move(d,s)); + } +}; + + + + +#endif /* SRC_AMR_GRID_DIST_AMR_KEY_HPP_ */ diff --git a/src/Amr/grid_dist_amr_key_iterator.hpp b/src/Amr/grid_dist_amr_key_iterator.hpp new file mode 100644 index 00000000..02905e6b --- /dev/null +++ b/src/Amr/grid_dist_amr_key_iterator.hpp @@ -0,0 +1,136 @@ +/* + * grid_amr_dist_key_iterator.hpp + * + * Created on: Sep 22, 2017 + * Author: i-bird + */ + +#ifndef SRC_AMR_GRID_DIST_AMR_KEY_ITERATOR_HPP_ +#define SRC_AMR_GRID_DIST_AMR_KEY_ITERATOR_HPP_ + +#include "Vector/map_vector.hpp" +#include "Grid/Iterators/grid_dist_id_iterator.hpp" +#include "grid_dist_amr_key.hpp" + +template> +class grid_dist_amr_key_iterator +{ + //! Array of grid iterators + openfpm::vector & git; + + //! actual it type + struct actual_it + { + it_type & it; + }; + + //! Actual distributed grid iterator + it_type * a_it; + + //! iterator pointer + size_t g_c; + + + + /*! \brief from g_c increment g_c until you find a valid grid + * + */ + void selectValidGrid() + { + // When the grid has size 0 potentially all the other informations are garbage + while (g_c < git.size() && git.get(g_c).isNext() == false ) g_c++; + + // get the next grid iterator + if (g_c < git.size()) + { + a_it = &git.get(g_c); + } + } + +public: + + /*! \brief Constructor + * + * \param git vector of iterator + * + */ + grid_dist_amr_key_iterator(openfpm::vector & git) + :git(git),g_c(0) + { + a_it = &git.get(0); + + selectValidGrid(); + } + + + //! Destructor + ~grid_dist_amr_key_iterator() + { + } + + + /*! \brief Get the next element + * + * \return the next grid_key + * + */ + inline grid_dist_amr_key_iterator & operator++() + { + ++(*a_it); + + // check if a_it is at the end + + if (a_it->isNext() == true) + {return *this;} + else + { + // switch to the new iterator + g_c++; + + selectValidGrid(); + } + + return *this; + } + + /*! \brief Is there a next point + * + * \return true is there is a next point + * + */ + inline bool isNext() + { + return g_c < git.size(); + } + + /*! \brief Return the actual AMR grid iterator point + * + * + */ + inline grid_dist_amr_key get() + { + return grid_dist_amr_key(g_c,a_it->get()); + } + + /*! \brief Return the actual global grid position in the AMR struct in global + * coordinates + * + * + */ + inline grid_key_dx getGKey() + { + return git.get(g_c).getGKey(a_it->get()); + } + + /*! \brief Return the level at which we are + * + * + */ + inline size_t getLvl() const + { + return g_c; + } +}; + + +#endif /* SRC_AMR_GRID_DIST_AMR_KEY_ITERATOR_HPP_ */ diff --git a/src/Amr/grid_dist_amr_unit_tests.cpp b/src/Amr/grid_dist_amr_unit_tests.cpp new file mode 100644 index 00000000..1527d83c --- /dev/null +++ b/src/Amr/grid_dist_amr_unit_tests.cpp @@ -0,0 +1,993 @@ +/* + * grid_dist_amr_dist_unit_tests.cpp + * + * Created on: Sep 21, 2017 + * Author: i-bird + */ + + +#define BOOST_TEST_DYN_LINK + +#include +#include "grid_dist_amr.hpp" + +BOOST_AUTO_TEST_SUITE( grid_dist_amr_test ) + +/*! \brief Coarsest levels of the grid + * + * \param domain Simulation domain + * \param coars_g coarsest grid resolution + * \param n_lvl number of levels + * + */ +template +void Test3D_amr_create_levels(grid_amr & amr_g, Box<3,float> & domain, size_t coars_g, size_t n_lvl) +{ + size_t g_sz[3] = {coars_g,coars_g,coars_g}; + + size_t tot_c = (coars_g - 1)*(coars_g - 1)*(coars_g - 1); + size_t correct_result = 0; + size_t correct_result_cell = 0; + size_t fact = 1; + + for (size_t i = 0 ; i < n_lvl ; i++) + { + correct_result += coars_g*coars_g*coars_g; + correct_result_cell += tot_c*fact; + coars_g = 2*(coars_g - 1) + 1; + fact *= 8; + } + + amr_g.initLevels(n_lvl,g_sz); + + + for (size_t i = 0 ; i < amr_g.getNLvl() ; i++) + { + // Fill the AMR with something + + size_t count = 0; + + auto it = amr_g.getGridIterator(i); + + while (it.isNext()) + { + auto key = it.get_dist(); + auto akey = amr_g.getAMRKey(i,key); + + amr_g.template insert<0>(akey) = 3.0; + + count++; + + ++it; + } + } + + // Iterate across all the levels initialized + auto it = amr_g.getDomainIterator(); + + size_t count = 0; + + while (it.isNext()) + { + count++; + + ++it; + } + + Vcluster<> & v_cl = create_vcluster(); + + v_cl.sum(count); + v_cl.execute(); + + BOOST_REQUIRE_EQUAL(count,correct_result); + + auto itc = amr_g.getDomainIteratorCells(); + + size_t count_c = 0; + + while (itc.isNext()) + { + count_c++; + + ++itc; + } + + v_cl.sum(count_c); + v_cl.execute(); + + auto it_level = amr_g.getDomainIteratorCells(3); + + while (it_level.isNext()) + { + auto key = it_level.get(); + + amr_g.template get<0>(3,key); + + ++it_level; + } + + BOOST_REQUIRE_EQUAL(count_c,correct_result_cell); +} + + + +template +inline bool gr_is_inside(const grid_key_dx & key, const size_t (& sz)[dim]) +{ + for (size_t i = 0 ; i < dim ; i++) + { + if (key.get(i) >= (long int)sz[i] || key.get(i) < 0) + { + return false; + } + } + + return true; +} + +template +void Test3D_amr_child_parent_get_no_periodic(grid & amr_g, Box<3,float> & domain, size_t coars_g, size_t n_lvl) +{ + const int x = 0; + const int y = 1; + const int z = 2; + + size_t g_sz[3] = {coars_g,coars_g,coars_g}; + + size_t tot = coars_g*coars_g*coars_g; + size_t correct_result = 0; + size_t fact = 1; + + for (size_t i = 0 ; i < n_lvl ; i++) + { + correct_result += tot*fact; + fact *= 8; + } + + amr_g.initLevels(n_lvl,g_sz); + + //////// Add something ///// + + for (size_t i = 0 ; i < amr_g.getNLvl() ; i++) + { + // Fill the AMR with something + + size_t count = 0; + + auto it = amr_g.getGridIterator(i); + + while (it.isNext()) + { + auto key = it.get_dist(); + auto akey = amr_g.getAMRKey(i,key); + + amr_g.template insert<0>(akey) = 3.0; + + count++; + + ++it; + } + } + + //////////////////////////// + + std::string test = amr_g.getSpacing(0).toString(); + + // Iterate across all the levels initialized + auto it = amr_g.getDomainIterator(); + + while (it.isNext()) + { + auto key = it.get(); + auto gkey = it.getGKey(); + + amr_g.template insert<0>(key) = gkey.get(0); + amr_g.template insert<1>(key) = gkey.get(1); + amr_g.template insert<2>(key) = gkey.get(2); + + ++it; + } + + amr_g.template ghost_get<0,1,2>(); + + // now we check that move space work + + auto it2 = amr_g.getDomainIterator(); + + bool match = true; + while (it2.isNext()) + { + auto key = it2.get(); + auto gkey = it2.getGKey(); + + auto key_px = key.moveSpace(x,1); + auto key_gpx = amr_g.getGKey(key_px); + + auto key_mx = key.moveSpace(x,-1); + auto key_gmx = amr_g.getGKey(key_mx); + + auto key_py = key.moveSpace(y,1); + auto key_gpy = amr_g.getGKey(key_py); + + auto key_my = key.moveSpace(y,-1); + auto key_gmy = amr_g.getGKey(key_my); + + auto key_pz = key.moveSpace(z,1); + auto key_gpz = amr_g.getGKey(key_pz); + + auto key_mz = key.moveSpace(z,-1); + auto key_gmz = amr_g.getGKey(key_mz); + + if (gr_is_inside(key_gpx,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_px) == gkey.get(0) + 1; + match &= amr_g.template get<1>(key_px) == gkey.get(1); + match &= amr_g.template get<2>(key_px) == gkey.get(2); + } + + if (gr_is_inside(key_gmx,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_mx) == gkey.get(0) - 1; + match &= amr_g.template get<1>(key_mx) == gkey.get(1); + match &= amr_g.template get<2>(key_mx) == gkey.get(2); + } + + if (gr_is_inside(key_gpy,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_py) == gkey.get(0); + match &= amr_g.template get<1>(key_py) == gkey.get(1) + 1; + match &= amr_g.template get<2>(key_py) == gkey.get(2); + } + + if (gr_is_inside(key_gmy,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_my) == gkey.get(0); + match &= amr_g.template get<1>(key_my) == gkey.get(1) - 1; + match &= amr_g.template get<2>(key_my) == gkey.get(2); + } + + if (gr_is_inside(key_gpz,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_pz) == gkey.get(0); + match &= amr_g.template get<1>(key_pz) == gkey.get(1); + match &= amr_g.template get<2>(key_pz) == gkey.get(2) + 1; + } + + if (gr_is_inside(key_gmz,amr_g.getGridInfoVoid(it2.getLvl()).getSize()) == true) + { + match &= amr_g.template get<0>(key_mz) == gkey.get(0); + match &= amr_g.template get<1>(key_mz) == gkey.get(1); + match &= amr_g.template get<2>(key_mz) == gkey.get(2) - 1; + } + + if (match == false) + { + int debug = 0; + debug++; + } + + // Test to go to all the levels down + + size_t lvl = it2.getLvl(); + + if (lvl < amr_g.getNLvl() - 1) + { + auto key_l1 = key; + amr_g.moveLvlDw(key_l1); + auto key_gl1 = amr_g.getGKey(key_l1); + + for (size_t s = 0 ; s < 3 ; s++) + { + match &= key_gl1.get(s) >> 1 == gkey.get(s); + match &= amr_g.template get<0>(key_l1) == key_gl1.get(0); + match &= amr_g.template get<1>(key_l1) == key_gl1.get(1); + match &= amr_g.template get<2>(key_l1) == key_gl1.get(2); + } + } + + if (lvl != 0) + { + auto key_l1 = key; + amr_g.moveLvlUp(key_l1); + auto key_gl1 = amr_g.getGKey(key_l1); + + for (size_t s = 0 ; s < 3 ; s++) + { + match &= gkey.get(s) >> 1 == key_gl1.get(s); + + match &= amr_g.template get<0>(key_l1) == key_gl1.get(0); + match &= amr_g.template get<1>(key_l1) == key_gl1.get(1); + match &= amr_g.template get<2>(key_l1) == key_gl1.get(2); + } + } + + if (match == false) + { + int debug = 0; + debug++; + } + + ++it2; + } + + BOOST_REQUIRE_EQUAL(match,true); +} + + +template +void Test3D_amr_child_parent_get_periodic(grid & amr_g, Box<3,float> & domain, size_t coars_g, size_t n_lvl) +{ + const int x = 0; + const int y = 1; + const int z = 2; + + size_t g_sz[3] = {coars_g,coars_g,coars_g}; + + size_t tot = coars_g*coars_g*coars_g; + size_t correct_result = 0; + size_t fact = 1; + + for (size_t i = 0 ; i < n_lvl ; i++) + { + correct_result += tot*fact; + fact *= 8; + } + + amr_g.initLevels(n_lvl,g_sz); + + //////// Add something ///// + + for (size_t i = 0 ; i < amr_g.getNLvl() ; i++) + { + // Fill the AMR with something + + size_t count = 0; + + auto it = amr_g.getGridIterator(i); + + while (it.isNext()) + { + auto key = it.get_dist(); + auto akey = amr_g.getAMRKey(i,key); + + amr_g.template insert<0>(akey) = 3.0; + + count++; + + ++it; + } + } + + //////////////////////////// + + std::string test = amr_g.getSpacing(0).toString(); + + // Iterate across all the levels initialized + auto it = amr_g.getDomainIterator(); + + while (it.isNext()) + { + auto key = it.get(); + auto gkey = it.getGKey(); + + amr_g.template insert<0>(key) = gkey.get(0); + amr_g.template insert<1>(key) = gkey.get(1); + amr_g.template insert<2>(key) = gkey.get(2); + + ++it; + } + + amr_g.template ghost_get<0,1,2>(); + + // now we check that move space work + + auto it2 = amr_g.getDomainIterator(); + + bool match = true; + while (it2.isNext()) + { + auto key = it2.get(); + auto gkey = it2.getGKey(); + + auto key_px = key.moveSpace(x,1); + auto key_mx = key.moveSpace(x,-1); + auto key_py = key.moveSpace(y,1); + auto key_my = key.moveSpace(y,-1); + auto key_pz = key.moveSpace(z,1); + auto key_mz = key.moveSpace(z,-1); + + match &= amr_g.template get<0>(key_px) == openfpm::math::positive_modulo(gkey.get(0) + 1,amr_g.getGridInfoVoid(it2.getLvl()).size(0)); + match &= amr_g.template get<1>(key_px) == gkey.get(1); + match &= amr_g.template get<2>(key_px) == gkey.get(2); + + match &= amr_g.template get<0>(key_mx) == openfpm::math::positive_modulo(gkey.get(0) - 1,amr_g.getGridInfoVoid(it2.getLvl()).size(0)); + match &= amr_g.template get<1>(key_mx) == gkey.get(1); + match &= amr_g.template get<2>(key_mx) == gkey.get(2); + + match &= amr_g.template get<0>(key_py) == gkey.get(0); + match &= amr_g.template get<1>(key_py) == openfpm::math::positive_modulo(gkey.get(1) + 1,amr_g.getGridInfoVoid(it2.getLvl()).size(1)); + match &= amr_g.template get<2>(key_py) == gkey.get(2); + + match &= amr_g.template get<0>(key_my) == gkey.get(0); + match &= amr_g.template get<1>(key_my) == openfpm::math::positive_modulo(gkey.get(1) - 1,amr_g.getGridInfoVoid(it2.getLvl()).size(1)); + match &= amr_g.template get<2>(key_my) == gkey.get(2); + + match &= amr_g.template get<0>(key_pz) == gkey.get(0); + match &= amr_g.template get<1>(key_pz) == gkey.get(1); + match &= amr_g.template get<2>(key_pz) == openfpm::math::positive_modulo(gkey.get(2) + 1,amr_g.getGridInfoVoid(it2.getLvl()).size(2)); + + match &= amr_g.template get<0>(key_mz) == gkey.get(0); + match &= amr_g.template get<1>(key_mz) == gkey.get(1); + match &= amr_g.template get<2>(key_mz) == openfpm::math::positive_modulo(gkey.get(2) - 1,amr_g.getGridInfoVoid(it2.getLvl()).size(2)); + + // Test to go to all the levels down + + size_t lvl = it2.getLvl(); + + if (lvl < amr_g.getNLvl() - 1) + { + auto key_l1 = key; + amr_g.moveLvlDw(key_l1); + auto key_gl1 = amr_g.getGKey(key_l1); + + for (size_t s = 0 ; s < 3 ; s++) + { + match &= key_gl1.get(s) >> 1 == gkey.get(s); + match &= amr_g.template get<0>(key_l1) == key_gl1.get(0); + match &= amr_g.template get<1>(key_l1) == key_gl1.get(1); + match &= amr_g.template get<2>(key_l1) == key_gl1.get(2); + } + } + + if (lvl != 0) + { + auto key_l1 = key; + amr_g.moveLvlUp(key_l1); + auto key_gl1 = amr_g.getGKey(key_l1); + + for (size_t s = 0 ; s < 3 ; s++) + { + match &= gkey.get(s) >> 1 == key_gl1.get(s); + + match &= amr_g.template get<0>(key_l1) == key_gl1.get(0); + match &= amr_g.template get<1>(key_l1) == key_gl1.get(1); + match &= amr_g.template get<2>(key_l1) == key_gl1.get(2); + } + } + + ++it2; + } + + BOOST_REQUIRE_EQUAL(match,true); +} + +template +void Test3D_amr_ghost_it(grid & amr_g, Box<3,float> & domain, size_t coars_g, size_t n_lvl) +{ + size_t g_sz[3] = {coars_g,coars_g,coars_g}; + + size_t tot = coars_g*coars_g*coars_g; + size_t correct_result = 0; + size_t fact = 1; + + for (size_t i = 0 ; i < n_lvl ; i++) + { + correct_result += tot*fact; + fact *= 8; + } + + amr_g.initLevels(n_lvl,g_sz); + + //////// Add something ///// + + for (size_t i = 0 ; i < amr_g.getNLvl() ; i++) + { + // Fill the AMR with something + + size_t count = 0; + + auto it = amr_g.getGridGhostIterator(i); + + while (it.isNext()) + { + auto key = it.get_dist(); + auto gkey = it.get(); + auto akey = amr_g.getAMRKey(i,key); + + amr_g.template insert<0>(akey) = gkey.get(0); + amr_g.template insert<1>(akey) = gkey.get(1); + amr_g.template insert<2>(akey) = gkey.get(2); + + if (gkey.get(0) == -1 || gkey.get(0) == (long int)amr_g.getGridInfoVoid(i).size(0) || + gkey.get(1) == -1 || gkey.get(1) == (long int)amr_g.getGridInfoVoid(i).size(1) || + gkey.get(2) == -1 || gkey.get(2) == (long int)amr_g.getGridInfoVoid(i).size(2)) + {count++;} + + ++it; + } + + size_t tot = (amr_g.getGridInfoVoid(i).size(0) + 2)* + (amr_g.getGridInfoVoid(i).size(1) + 2)* + (amr_g.getGridInfoVoid(i).size(2) + 2) - amr_g.getGridInfoVoid(i).size(); + + auto & v_cl = create_vcluster(); + + if (v_cl.size() == 1) + { + v_cl.sum(count); + v_cl.execute(); + + BOOST_REQUIRE_EQUAL(tot,count); + } + + bool match = true; + auto it2 = amr_g.getDomainIterator(i); + + while (it2.isNext()) + { + auto key = it2.get(); + + // move -x + + auto key_m1 = key.move(0,-1); + auto key_gm1 = it2.getGKey(key_m1); + match &= amr_g.template get<0>(i,key_m1) == key_gm1.get(0); + + // move +x + + auto key_p1 = key.move(0,1); + auto key_gp1 = it2.getGKey(key_p1); + match &= amr_g.template get<0>(i,key_p1) == key_gp1.get(0); + + // move -y + + key_m1 = key.move(1,-1); + key_gm1 = it2.getGKey(key_m1); + match &= amr_g.template get<1>(i,key_m1) == key_gm1.get(1); + + // move +y + + key_p1 = key.move(1,1); + key_gp1 = it2.getGKey(key_p1); + match &= amr_g.template get<1>(i,key_p1) == key_gp1.get(1); + + // move -z + + key_m1 = key.move(2,-1); + key_gm1 = it2.getGKey(key_m1); + match &= amr_g.template get<2>(i,key_m1) == key_gm1.get(2); + + // move +z + + key_p1 = key.move(2,1); + key_gp1 = it2.getGKey(key_p1); + match &= amr_g.template get<2>(i,key_p1) == key_gp1.get(2); + + ++it2; + } + + BOOST_REQUIRE_EQUAL(match,true); + } +} + + + +template +void Test3D_amr_domain_ghost_it(grid & amr_g, Box<3,float> & domain, size_t coars_g, size_t n_lvl) +{ + size_t g_sz[3] = {coars_g,coars_g,coars_g}; + + size_t tot = coars_g*coars_g*coars_g; + size_t correct_result = 0; + size_t fact = 1; + + for (size_t i = 0 ; i < n_lvl ; i++) + { + correct_result += tot*fact; + fact *= 8; + } + + amr_g.initLevels(n_lvl,g_sz); + + size_t total_all_level = 0; + + //////// Add something ///// + + for (size_t i = 0 ; i < amr_g.getNLvl() ; i++) + { + // Fill the AMR with something + + size_t count = 0; + + auto it = amr_g.getGridGhostIterator(i); + + while (it.isNext()) + { + auto key = it.get_dist(); + auto gkey = it.get(); + auto akey = amr_g.getAMRKey(i,key); + + amr_g.template insert<0>(akey) = gkey.get(0); + amr_g.template insert<1>(akey) = gkey.get(1); + amr_g.template insert<2>(akey) = gkey.get(2); + + count++; + + ++it; + } + + size_t tot = (amr_g.getGridInfoVoid(i).size(0) + 2)* + (amr_g.getGridInfoVoid(i).size(1) + 2)* + (amr_g.getGridInfoVoid(i).size(2) + 2); + + auto & v_cl = create_vcluster(); + + if (v_cl.size() == 1) + { + v_cl.sum(count); + v_cl.execute(); + + BOOST_REQUIRE_EQUAL(tot,count); + } + + size_t amr_cnt = 0; + bool match = true; + auto it2 = amr_g.getDomainGhostIterator(i); + + while (it2.isNext()) + { + auto key = it2.get(); + auto key_g = it2.getGKey(key); + match &= amr_g.template get<0>(i,key) == key_g.get(0); + + total_all_level++; + amr_cnt++; + + ++it2; + } + + BOOST_REQUIRE_EQUAL(amr_cnt,count); + BOOST_REQUIRE_EQUAL(match,true); + } + + // test the total iterator + + size_t gtot_count = 0; + auto tot_it = amr_g.getDomainGhostIterator(); + + while (tot_it.isNext()) + { + gtot_count++; + + ++tot_it; + } + + BOOST_REQUIRE_EQUAL(gtot_count,total_all_level); +} + +template +void Test3D_ghost_put(grid_amr & g_dist_amr, long int k) +{ + size_t sz[3] = {(size_t)k,(size_t)k,(size_t)k}; + + g_dist_amr.initLevels(4,sz); + + // Grid sm + grid_sm<3,void> info(sz); + + size_t count = 0; + + for (size_t i = 0 ; i < g_dist_amr.getNLvl() ; i++) + { + auto dom = g_dist_amr.getGridIterator(i); + + while (dom.isNext()) + { + auto key = dom.get_dist(); + + g_dist_amr.template insert<0>(i,key) = -6.0; + + // Count the points + count++; + + ++dom; + } + } + + // Set to zero the full grid + + { + auto dom = g_dist_amr.getDomainIterator(); + + while (dom.isNext()) + { + auto key = dom.get(); + + g_dist_amr.template insert<0>(key.moveSpace(0,1)) += 1.0; + g_dist_amr.template insert<0>(key.moveSpace(0,-1)) += 1.0; + g_dist_amr.template insert<0>(key.moveSpace(1,1)) += 1.0; + g_dist_amr.template insert<0>(key.moveSpace(1,-1)) += 1.0; + g_dist_amr.template insert<0>(key.moveSpace(2,1)) += 1.0; + g_dist_amr.template insert<0>(key.moveSpace(2,-1)) += 1.0; + + ++dom; + } + } + + bool correct = true; + + // Domain + Ghost iterator + auto dom_gi = g_dist_amr.getDomainIterator(); + + while (dom_gi.isNext()) + { + auto key = dom_gi.get(); + + correct &= (g_dist_amr.template get<0>(key) == 0); + + ++dom_gi; + } + + g_dist_amr.template ghost_put(); + + if (count != 0) + {BOOST_REQUIRE_EQUAL(correct, false);} + + // sync the ghosts + g_dist_amr.template ghost_get<0>(); + + correct = true; + + // Domain + Ghost iterator + auto dom_gi2 = g_dist_amr.getDomainIterator(); + + while (dom_gi2.isNext()) + { + auto key = dom_gi2.get(); + + correct &= (g_dist_amr.template get<0>(key) == 0); + + ++dom_gi2; + } + + BOOST_REQUIRE_EQUAL(correct, true); +} + +template struct Debug; + +BOOST_AUTO_TEST_CASE( grid_dist_amr_get_child_test_nop ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 16*16*16*create_vcluster().getProcessingUnits(); + k = std::pow(k, 1/3.); + + Ghost<3,long int> g(1); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g); + + Test3D_amr_child_parent_get_no_periodic(amr_g,domain3,k,4); +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_get_child_test_p ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 16*16*16*create_vcluster().getProcessingUnits(); + k = std::pow(k, 1/3.); + + periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC}; + + Ghost<3,long int> g(1); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g,bc); + + Test3D_amr_child_parent_get_periodic(amr_g,domain3,k,4); +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_test ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 16*16*16*create_vcluster().getProcessingUnits(); + k = std::pow(k, 1/3.); + + Ghost<3,float> g(0.05); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g); + + Test3D_amr_create_levels(amr_g,domain3,k,4); + + sgrid_dist_amr<3,float,aggregate> amr_g2(domain3,g); + + Test3D_amr_create_levels(amr_g2,domain3,k,4); +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_ghost_it_test ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 16*16*16*create_vcluster().getProcessingUnits(); + k = std::pow(k, 1/3.); + + Ghost<3,long int> g(1); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g); + + Test3D_amr_ghost_it(amr_g,domain3,k,4); + + sgrid_dist_amr<3,float,aggregate> amr_g2(domain3,g); + + Test3D_amr_ghost_it(amr_g2,domain3,k,4); + + for (size_t i = 0 ; i < amr_g2.getNLvl() ; i++) + {BOOST_REQUIRE(amr_g2.size_inserted(i) != 0ul);} + + amr_g2.clear(); + + for (size_t i = 0 ; i < amr_g2.getNLvl() ; i++) + {BOOST_REQUIRE_EQUAL(amr_g2.size_inserted(i),0ul);} +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_domain_ghost_it_test ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 16*16*16*create_vcluster().getProcessingUnits(); + k = std::pow(k, 1/3.); + + Ghost<3,long int> g(1); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g); + + Test3D_amr_domain_ghost_it(amr_g,domain3,k,4); + + sgrid_dist_amr<3,float,aggregate> amr_g2(domain3,g); + + Test3D_amr_domain_ghost_it(amr_g2,domain3,k,4); + + for (size_t i = 0 ; i < amr_g2.getNLvl() ; i++) + {BOOST_REQUIRE(amr_g2.size_inserted(i) != 0ul);} + + amr_g2.clear(); + + for (size_t i = 0 ; i < amr_g2.getNLvl() ; i++) + {BOOST_REQUIRE_EQUAL(amr_g2.size_inserted(i),0ul);} +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_get_child_test_low_res ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + long int k = 2; + + Ghost<3,long int> g(1); + grid_dist_amr<3,float,aggregate> amr_g(domain3,g); + + Test3D_amr_child_parent_get_no_periodic(amr_g,domain3,k,4); + + sgrid_dist_amr<3,float,aggregate> amr_g2(domain3,g); + + Test3D_amr_child_parent_get_no_periodic(amr_g2,domain3,k,4); +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_test_background_value ) +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + Ghost<3,long int> g(1); + + sgrid_dist_amr<3,float,aggregate> amr_g2(domain3,g); + + size_t g_sz[3] = {4,4,4}; + + amr_g2.initLevels(4,g_sz); + + aggregate bck; + bck.get<0>() = -57; + bck.get<1>() = -90; + bck.get<2>() = -123; + + amr_g2.setBackgroundValue(bck); + + // Get a non existent point to check that + // the background value work + + grid_dist_key_dx<3> key(0,grid_key_dx<3>({0,0,0})); + long int bck0 = amr_g2.get<0>(2,key); + BOOST_REQUIRE_EQUAL(bck0,-57); + long int bck1 = amr_g2.get<1>(2,key); + BOOST_REQUIRE_EQUAL(bck1,-90); + long int bck2 = amr_g2.get<2>(2,key); + BOOST_REQUIRE_EQUAL(bck2,-123); + + // Now we insert that point and we check the subsequent point + amr_g2.insert<0>(2,key) = 5; + + grid_dist_key_dx<3> key2(0,grid_key_dx<3>({1,0,0})); + bck0 = amr_g2.get<0>(2,key2); + BOOST_REQUIRE_EQUAL(bck0,-57); + bck1 = amr_g2.get<1>(2,key2); + BOOST_REQUIRE_EQUAL(bck1,-90); + bck2 = amr_g2.get<2>(2,key2); + BOOST_REQUIRE_EQUAL(bck2,-123); + + auto & g_dist_lvl2 = amr_g2.getDistGrid(2); + g_dist_lvl2.get_loc_grid(0).internal_clear_cache(); + + bck0 = amr_g2.get<0>(2,key2); + BOOST_REQUIRE_EQUAL(bck0,-57); + bck1 = amr_g2.get<1>(2,key2); + BOOST_REQUIRE_EQUAL(bck1,-90); + bck2 = amr_g2.get<2>(2,key2); + BOOST_REQUIRE_EQUAL(bck2,-123); + +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_get_domain_ghost_check ) +{ + // Test grid periodic + + Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0}); + + Vcluster<> & v_cl = create_vcluster(); + + if ( v_cl.getProcessingUnits() > 32 ) + {return;} + + long int k = 13; + + BOOST_TEST_CHECKPOINT( "Testing grid periodic k<=" << k ); + + // Ghost + Ghost<3,long int> g(1); + + // periodicity + periodicity<3> pr = {{PERIODIC,PERIODIC,PERIODIC}}; + + // Distributed grid with id decomposition + grid_dist_amr<3, float, aggregate> g_dist(domain,g,pr); + + Test3D_ghost_put(g_dist,k); + + // Distributed grid with id decomposition + sgrid_dist_amr<3, float, aggregate> sg_dist(domain,g,pr); + + Test3D_ghost_put(sg_dist,k); +} + +BOOST_AUTO_TEST_CASE( grid_dist_amr_ghost_put_create ) +{ + // Test grid periodic + + Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0}); + + Vcluster<> & v_cl = create_vcluster(); + + if ( v_cl.getProcessingUnits() > 32 ) + {return;} + + long int k = 13; + + BOOST_TEST_CHECKPOINT( "Testing grid periodic k<=" << k ); + + // Ghost + Ghost<3,long int> g(1); + + // periodicity + periodicity<3> pr = {{PERIODIC,PERIODIC,PERIODIC}}; + + // Distributed grid with id decomposition + grid_dist_amr<3, float, aggregate> g_dist(domain,g,pr); + + Test3D_ghost_put(g_dist,k); + + // Distributed grid with id decomposition + sgrid_dist_amr<3, float, aggregate> sg_dist(domain,g,pr); + + Test3D_ghost_put(sg_dist,k); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/src/Amr/tests/amr_base_unit_tests.cpp b/src/Amr/tests/amr_base_unit_tests.cpp new file mode 100644 index 00000000..973e6e6a --- /dev/null +++ b/src/Amr/tests/amr_base_unit_tests.cpp @@ -0,0 +1,171 @@ +/* + * amr_base_unit_test.cpp + * + * Created on: Oct 5, 2017 + * Author: i-bird + */ +#define BOOST_TEST_DYN_LINK +#include + +#include "Grid/grid_dist_id.hpp" +#include "Point_test.hpp" +#include "Grid/tests/grid_dist_id_util_tests.hpp" + +BOOST_AUTO_TEST_SUITE( amr_grid_dist_id_test ) + + +BOOST_AUTO_TEST_CASE( grid_dist_id_amr ) +{ + // Domain + Box<2,float> domain2({0.0,0.0},{1.0,1.0}); + + size_t sz[2] = {100,100}; + + // Ghost + Ghost<2,long int> g(1); + + // periodicity + periodicity<2> pr = {{PERIODIC,PERIODIC}}; + + openfpm::vector> C_draw; + C_draw.add(Box<2,long int>({20,20},{50,24})); + C_draw.add(Box<2,long int>({51,20},{60,24})); + C_draw.add(Box<2,long int>({61,20},{70,24})); + C_draw.add(Box<2,long int>({20,25},{24,66})); + C_draw.add(Box<2,long int>({15,67},{49,85})); + C_draw.add(Box<2,long int>({50,76},{70,81})); + C_draw.add(Box<2,long int>({30,25},{34,37})); + C_draw.add(Box<2,long int>({50,66},{70,70})); + + size_t volume_key = 0; + for (size_t i = 0 ; i < C_draw.size() ; i++) + { + volume_key += Box<2,long int>(C_draw.get(i)).getVolumeKey(); + } + + // Distributed grid with id decomposition + grid_dist_id<2,float,Point_test> g_dist(sz,domain2,g,pr,C_draw); + + // fill with gkey + + auto git = g_dist.getDomainIterator(); + grid_sm<2,void> gs(sz); + + size_t count = 0; + + while (git.isNext()) + { + auto key = git.get(); + auto gkey = git.getGKey(key); + + g_dist.template get<0>(key) = gs.LinId(gkey); + + count++; + + ++git; + } + + Vcluster<> & vcl = create_vcluster(); + + vcl.sum(count); + vcl.execute(); + + BOOST_REQUIRE_EQUAL(count,volume_key); + + g_dist.ghost_get<0>(); + + // Check it is correct + + bool check = true; + size_t check_count = 0; + + auto git2 = g_dist.getDomainGhostIterator(); + while (git2.isNext()) + { + auto key = git2.get(); + auto gkey = git2.getGKey(key); + + float value = g_dist.template get<0>(key); + + // check if the point is inside or outside the domain + + for (size_t k = 0; k < C_draw.size() ; k++) + { + if (Box<2,long int>(C_draw.get(k)).isInside(gkey.toPoint()) == true) + { + check &= value == gs.LinId(gkey); + + // get the gdb_ext + auto & gdb_ext = g_dist.getLocalGridsInfo(); + + for (size_t s = 0 ; s < gdb_ext.size() ; s++) + { + Box<2,long int> bx = gdb_ext.get(s).Dbox; + bx += gdb_ext.get(s).origin; + if (bx.isInside(gkey.toPoint())) + { + check_count++; + break; + } + } + break; + } + } + + ++git2; + } + + vcl.sum(check_count); + vcl.execute(); + + BOOST_REQUIRE_EQUAL(check,true); + BOOST_REQUIRE(check_count >= volume_key); +} + +BOOST_AUTO_TEST_CASE( amr_grid_dist_id_iterator_test_use_2D) +{ + // Domain + Box<2,float> domain({0.0,0.0},{1.0,1.0}); + +#ifdef TEST_COVERAGE_MODE + long int k = 256*256*create_vcluster().getProcessingUnits(); +#else + long int k = 1024*1024*create_vcluster().getProcessingUnits(); +#endif + k = std::pow(k, 1/2.); + + long int big_step = k / 30; + big_step = (big_step == 0)?1:big_step; + long int small_step = 21; + + print_test( "AMR Testing 2D full grid k<=",k); + + // 2D test + for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step ) + { + BOOST_TEST_CHECKPOINT( "AMR Testing 2D full grid k=" << k ); + + //! [Create and access a distributed grid] + + // grid size + size_t sz[2]; + sz[0] = k; + sz[1] = k; + + // periodicity + periodicity<2> pr = {{PERIODIC,PERIODIC}}; + + // Ghost + Ghost<2,long int> g(1); + + openfpm::vector> bx_def; + bx_def.add(Box<2,long int>({0,0},{k-1,k-1})); + + // Distributed grid with id decomposition + grid_dist_id<2, float, aggregate> g_dist(sz,domain,g,pr,bx_def); + + Test2D_core(g_dist,sz,k); + } +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a56b261a..0deec844 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,7 +4,9 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR) ########################### Executables if(CUDA_FOUND) - set(CUDA_SOURCES Vector/cuda/vector_dist_gpu_MP_tests.cu + set(CUDA_SOURCES + Grid/tests/sgrid_dist_id_gpu_unit_tests.cu + Vector/cuda/vector_dist_gpu_MP_tests.cu Vector/cuda/vector_dist_cuda_func_test.cu Decomposition/cuda/decomposition_cuda_tests.cu Vector/cuda/vector_dist_gpu_unit_tests.cu diff --git a/src/Decomposition/Domain_icells_cart.hpp b/src/Decomposition/Domain_icells_cart.hpp index 9055ff39..20c0a66b 100644 --- a/src/Decomposition/Domain_icells_cart.hpp +++ b/src/Decomposition/Domain_icells_cart.hpp @@ -156,7 +156,7 @@ class domain_icell_calculator openfpm::vector_sparse_gpu> vs; openfpm::vector_sparse_gpu> vsi; - vs.template getBackground<0>() = 0; + vs.template setBackground<0>(0); // insert Domain cells @@ -182,7 +182,7 @@ class domain_icell_calculator CUDA_LAUNCH((insert_icell),ite,vsi.toKernel(),cld,ite.start,p2); - vsi.template flush<>(v_cl.getmgpuContext(),flust_type::FLUSH_ON_DEVICE); + vsi.template flush<>(v_cl.getmgpuContext(),flush_type::FLUSH_ON_DEVICE); } // calculate the number of kernel launch @@ -209,8 +209,8 @@ class domain_icell_calculator CUDA_LAUNCH(insert_remove_icell,ite,vs.toKernel(),vsi.toKernel(),cld,ite.start,p2); - vs.template flush<>(v_cl.getmgpuContext(),flust_type::FLUSH_ON_DEVICE); - vsi.flush_remove(v_cl.getmgpuContext(),flust_type::FLUSH_ON_DEVICE); + vs.template flush<>(v_cl.getmgpuContext(),flush_type::FLUSH_ON_DEVICE); + vsi.flush_remove(v_cl.getmgpuContext(),flush_type::FLUSH_ON_DEVICE); } vs.swapIndexVector(icells); diff --git a/src/Grid/Iterators/grid_dist_id_iterator_dec.hpp b/src/Grid/Iterators/grid_dist_id_iterator_dec.hpp index 91fe79e7..3f5d37fc 100644 --- a/src/Grid/Iterators/grid_dist_id_iterator_dec.hpp +++ b/src/Grid/Iterators/grid_dist_id_iterator_dec.hpp @@ -155,6 +155,56 @@ class grid_dist_id_iterator_dec { } + /*! \brief Return true if we point to a valid grid + * + * \return true if valid grid + * + */ + inline bool isNextGrid() + { + return g_c < gdb_ext.size(); + } + + /*! \brief Return the index of the grid in which we are iterating + * + * + */ + inline size_t getGridId() + { + return g_c; + } + + /*! \brief next grid + * + * + */ + inline void nextGrid() + { + g_c++; + selectValidGrid(); + } + + /*! \brief Return the actual pointed grid + * + * \return the grid index + * + */ + inline Box getGridBox() + { + Box bx; + + auto start = a_it.getStart(); + auto stop = a_it.getStop(); + + for (int i = 0 ; i < Decomposition::dims ; i++) + { + bx.setHigh(i,stop.get(i)); + bx.setLow(i,start.get(i)); + } + + return bx; + } + /*! \brief Get the next element * * \return the next grid_key diff --git a/src/Grid/cuda/grid_dist_id_kernels.cuh b/src/Grid/cuda/grid_dist_id_kernels.cuh new file mode 100644 index 00000000..9c4e1023 --- /dev/null +++ b/src/Grid/cuda/grid_dist_id_kernels.cuh @@ -0,0 +1,20 @@ +/* + * grid_dist_id_kernels.cuh + * + * Created on: Jun 25, 2019 + * Author: i-bird + */ + +#ifndef GRID_DIST_ID_KERNELS_CUH_ +#define GRID_DIST_ID_KERNELS_CUH_ + + +template +__global__ void grid_apply_functor(grid_type g, ite_gpu ite, func_t f, args_t ... args) +{ + f(g,ite,args...); +} + + + +#endif /* GRID_DIST_ID_KERNELS_CUH_ */ diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp index 0365a2e8..ded93aca 100644 --- a/src/Grid/grid_dist_id.hpp +++ b/src/Grid/grid_dist_id.hpp @@ -23,8 +23,11 @@ #include "grid_dist_id_comm.hpp" #include "HDF5_wr/HDF5_wr.hpp" #include "SparseGrid/SparseGrid.hpp" +#ifdef __NVCC__ +#include "SparseGridGpu/SparseGridGpu.hpp" +#include "cuda/grid_dist_id_kernels.cuh" +#endif -template struct Debug; //! Internal ghost box sent to construct external ghost box into the other processors template @@ -115,6 +118,9 @@ class grid_dist_id : public grid_dist_id_comm> cd_sm; + //! size of the insert pool on gpu + size_t gpu_insert_pool_size; + //! Communicator class Vcluster<> & v_cl; @@ -553,11 +559,11 @@ class grid_dist_id : public grid_dist_id_comm::meta_copy_(bv,loc_grid.get(i).getBackgroundValue());} } + /*! \brief set the background value + * + * You can use this function make sense in case of sparse in case of dense + * it does nothing + * + */ + template + void setBackgroundValue(const typename boost::mpl::at>::type & bv) + { + for (size_t i = 0 ; i < loc_grid.size() ; i++) + {loc_grid.get(i).template setBackgroundValue

(bv);} + } + /*! \brief Return the local total number of points inserted in the grid * * in case of dense grid it return the number of local points, in case of @@ -1992,6 +2011,16 @@ public: return loc_grid.get(v1.getSub()).remove_no_flush(v1.getKey()); } + + template + void flush(flush_type opt = flush_type::FLUSH_ON_HOST) + { + for (size_t i = 0 ; i < loc_grid.size() ; i++) + { + loc_grid.get(i).template flush(v_cl.getmgpuContext(),opt); + } + } + /*! \brief remove an element in the grid * * In case of dense grid this function print a warning, in case of sparse @@ -2654,6 +2683,73 @@ public: return this->ig_box; } +#ifdef __NVCC__ + + /*! \brief Set the size of the gpu insert buffer pool + * + * \param size of the insert pool + * + */ + void setInsertBuffer(size_t n_pool) + { + gpu_insert_pool_size = n_pool; + } + + template + void iterateGridGPU(it_t & it, args_t ... args) + { + while(it.isNextGrid()) + { + Box b = it.getGridBox(); + + size_t i = it.getGridId(); + + auto ite = loc_grid.get(i).getGridGPUIterator(b.getKP1(),b.getKP2()); + + loc_grid.get(i).setGPUInsertBuffer(ite.nblocks(),gpu_insert_pool_size); + loc_grid.get(i).initializeGPUInsertBuffer(); + + grid_apply_functor<<>>(loc_grid.get(i).toKernel(),ite,func_t(),args...); + + it.nextGrid(); + } + } + + template + void iterateGPU(args_t ... args) + { + for (int i = 0 ; i < loc_grid.size() ; i++) + { + auto & sp = loc_grid.get(i); + + // TODO Launch a kernel on every sparse grid GPU + } + } + + /*! \brief Move the memory from the device to host memory + * + */ + template void deviceToHost() + { + for (size_t i = 0 ; i < loc_grid.size() ; i++) + { + loc_grid.get(i).template deviceToHost(); + } + } + + /*! \brief Move the memory from the device to host memory + * + */ + template void hostToDevice() + { + for (size_t i = 0 ; i < loc_grid.size() ; i++) + { + loc_grid.get(i).template hostToDevice(); + } + } + +#endif + //! Define friend classes //\cond @@ -2665,4 +2761,9 @@ public: template using sgrid_dist_id = grid_dist_id,HeapMemory,sgrid_cpu>; +#ifdef __NVCC__ +template +using sgrid_dist_id_gpu = grid_dist_id,CudaMemory,SparseGridGpu>; +#endif + #endif diff --git a/src/Grid/tests/grid_dist_id_dlb_unit_test.cpp b/src/Grid/tests/grid_dist_id_dlb_unit_test.cpp new file mode 100644 index 00000000..bad7d6e1 --- /dev/null +++ b/src/Grid/tests/grid_dist_id_dlb_unit_test.cpp @@ -0,0 +1,242 @@ +/* + * grid_dist_id_dlb_unit_test.cpp + * + * Created on: May 4, 2018 + * Author: i-bird + */ + + +#define BOOST_TEST_DYN_LINK +#include + +#include "Point_test.hpp" +#include "Grid/grid_dist_id.hpp" +#include "data_type/aggregate.hpp" +#include "grid_dist_id_util_tests.hpp" +#include "Vector/vector_dist.hpp" + +BOOST_AUTO_TEST_SUITE( grid_dist_id_dlb_test ) + +template +void test_vector_grid_dlb() +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + Ghost<3,long int> g(1); + + size_t sz[3] = {37,37,37}; + + grid gdist(sz,domain3,g,DEC_GRAN(128)); + + + aggregate bck; + bck.template get<0>() = -57; + bck.template get<1>() = -90; + bck.template get<2>() = -123; + + gdist.setBackgroundValue(bck); + + vector vd(gdist.getDecomposition(),0); + + // we fill the grid with a gaussian + + auto it = gdist.getGridIterator(); + + while (it.isNext()) + { + auto p = it.get_dist(); + auto gkey = it.get(); + + gdist.template insert<0>(p) = gkey.get(0); + gdist.template insert<1>(p) = gkey.get(1); + gdist.template insert<2>(p) = gkey.get(2); + + ++it; + } + + // fill a sphere of particles + + auto it2 = vd.getGridIterator(sz); + + while (it2.isNext()) + { + auto gkey = it2.get(); + + vd.add(); + + float x = 0.2*gkey.get(0)*it2.getSpacing(0) + 0.05; + float y = 0.2*gkey.get(1)*it2.getSpacing(1) + 0.05; + float z = 0.2*gkey.get(2)*it2.getSpacing(2) + 0.05; + + vd.getLastPos()[0] = x; + vd.getLastPos()[1] = y; + vd.getLastPos()[2] = z; + + ++it2; + } + + size_t n_step = 50; + for (size_t i = 0; i < n_step ; i++) + { + vd.map(); + vd.addComputationCosts(); + vd.getDecomposition().decompose(); + vd.map(); + + gdist.getDecomposition() = vd.getDecomposition(); + gdist.map(); + + // Check + + bool check = true; + auto it = gdist.getDomainIterator(); + + + while (it.isNext()) + { + auto p = it.get(); + auto gkey = it.getGKey(p); + + check &= gdist.template get<0>(p) == gkey.get(0); + check &= gdist.template get<1>(p) == gkey.get(1); + check &= gdist.template get<2>(p) == gkey.get(2); + + ++it; + } + + BOOST_REQUIRE_EQUAL(check,true); + + // Calculate shift vector + + double t2 = 6.28*(double)(i+1)/n_step; + double t = 6.28*(double)i/n_step; + double v[3]; + v[0] = 0.7*fabs(sin(t2)) - 0.7*fabs(sin(t)); + v[1] = 0.7*fabs(sin(1.7*t2)) - 0.7*fabs(sin(1.7*t)); + v[2] = 0.7*fabs(sin(2.5*t2)) - 0.7*fabs(sin(2.5*t)); + + auto it2 = vd.getDomainIterator(); + + while (it2.isNext()) + { + auto p = it2.get(); + + vd.getPos(p)[0] += v[0]; + vd.getPos(p)[1] += v[1]; + vd.getPos(p)[2] += v[2]; + + ++it2; + } + vd.map(); + } +} + +struct GaussianDLB +{ + double t = 0.0; + + size_t resolution(const Point<3,float> & p) + { + float x0 = fabs(sin(t)); + float y0 = fabs(cos(t)); + float z0 = fabs(sin(3.0*t)); + + return 100 * exp( - ((p.get(0) - x0)*(p.get(0) - x0) + (p.get(1) - y0)*(p.get(1) - y0) + (p.get(2) - z0)*(p.get(2) - z0)) / 0.3); + } + + double distributionTol() + { + return 1.01; + } +}; + +template +void test_vector_grid_dlb_resolution() +{ + // Domain + Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0}); + + Ghost<3,long int> g(1); + + size_t sz[3] = {37,37,37}; + + grid gdist(sz,domain3,g,DEC_GRAN(128)); + + + aggregate bck; + bck.template get<0>() = -57; + bck.template get<1>() = -90; + bck.template get<2>() = -123; + + gdist.setBackgroundValue(bck); + + // we fill the grid with a gaussian + + auto it = gdist.getGridIterator(); + + while (it.isNext()) + { + auto p = it.get_dist(); + auto gkey = it.get(); + + gdist.template insert<0>(p) = gkey.get(0); + gdist.template insert<1>(p) = gkey.get(1); + gdist.template insert<2>(p) = gkey.get(2); + + ++it; + } + + GaussianDLB gdlb; + gdlb.t = 0.0; + + size_t n_step = 50; + for (size_t i = 0; i < n_step ; i++) + { + gdlb.t = (float)i/n_step; + gdist.addComputationCosts(gdlb); + gdist.getDecomposition().decompose(); + gdist.map(); + + gdist.write_frame("sgrid",i); + + // Check + + bool check = true; + auto it = gdist.getDomainIterator(); + + + while (it.isNext()) + { + auto p = it.get(); + auto gkey = it.getGKey(p); + + check &= gdist.template get<0>(p) == gkey.get(0); + check &= gdist.template get<1>(p) == gkey.get(1); + check &= gdist.template get<2>(p) == gkey.get(2); + + ++it; + } + + BOOST_REQUIRE_EQUAL(check,true); + } +} + +BOOST_AUTO_TEST_CASE( grid_dist_dlb_test ) +{ + typedef sgrid_dist_id<3,float,aggregate> grid_sparse; + typedef vector_dist<3,float,aggregate > particles; + + test_vector_grid_dlb(); +} + +BOOST_AUTO_TEST_CASE( grid_dist_dlb_test_resolution ) +{ + typedef sgrid_dist_id<3,float,aggregate> grid_sparse; + typedef vector_dist<3,float,aggregate > particles; + + test_vector_grid_dlb_resolution(); +} + +BOOST_AUTO_TEST_SUITE_END() + diff --git a/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu b/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu new file mode 100644 index 00000000..083830f6 --- /dev/null +++ b/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu @@ -0,0 +1,202 @@ +#define BOOST_TEST_DYN_LINK + +#include +#include "Grid/grid_dist_id.hpp" + +BOOST_AUTO_TEST_SUITE( sgrid_gpu_test_suite ) + +template +struct insert_kernel +{ + template + __device__ void operator()(SparseGridGpu_type & sg, ite_gpu & ite, float c) + { + sg.init(); + + const auto bDimX = blockDim.x; + const auto bDimY = blockDim.y; + const auto bDimZ = blockDim.z; + const auto bIdX = blockIdx.x; + const auto bIdY = blockIdx.y; + const auto bIdZ = blockIdx.z; + const auto tIdX = threadIdx.x; + const auto tIdY = threadIdx.y; + const auto tIdZ = threadIdx.z; + int x = bIdX * bDimX + tIdX; + int y = bIdY * bDimY + tIdY; + int z = bIdZ * bDimZ + tIdZ; + + if (x+ite.start.get(0) > ite.stop.get(0)) + {return;} + if (SparseGridGpu_type::d >= 2 && y+ite.start.get(1) > ite.stop.get(1)) + {return;} + if (SparseGridGpu_type::d >= 3 && z+ite.start.get(1) > ite.stop.get(2)) + {return;} + + grid_key_dx coord({x+ite.start.get(0), y+ite.start.get(1), z+ite.start.get(2)}); + + // size_t pos = sg.getLinId(coord); + // printf("insertValues: bDim=(%d,%d), bId=(%d,%d), tId=(%d,%d) : " + // "pos=%ld, coord={%d,%d}, value=%d\n", + // bDimX, bDimY, + // bIdX, bIdY, + // tIdX, tIdY, + // pos, + // x, y, + // x); //debug + + sg.template insert

(coord) = c; + + __syncthreads(); + + sg.flush_block_insert(); + + // Compiler avoid warning + y++; + z++; + } +}; + +template +struct stencil_kernel +{ + template + __device__ void operator()(SparseGridGpu_type & sg, ite_gpu & ite, float c) + { + // TODO + } +}; + +BOOST_AUTO_TEST_CASE( sgrid_gpu_test_base ) +{ + size_t sz[2] = {17,17}; + periodicity<2> bc = {PERIODIC,PERIODIC}; + + Ghost<2,long int> g(1); + + Box<2,float> domain({0.0,0.0},{1.0,1.0}); + + sgrid_dist_id_gpu<2,float,aggregate> gdist(sz,domain,g,bc); + + gdist.template setBackgroundValue<0>(666); + + /////// CPU insert + +/* auto it = gdist.getGridIterator(box.getKP1(),box.getKP2()); + + while (it.isNext()) + { + auto p = it.get_dist(); + + gdist.template insert<0>(p) = 1.0; + + ++it; + } + + gdist.template flush<>(); + + Box<2,size_t> box2({0,0},{15,15}); + auto it2 = gdist.getGridIterator(box2.getKP1(),box2.getKP2()); + + while (it2.isNext()) + { + auto p = it2.get_dist(); + + std::cout << gdist.template get<0>(p) << std::endl; + + ++it2; + }*/ + + /////// host to device + + /////// GPU insert + flush + + Box<2,size_t> box({1,1},{1,1}); + auto it = gdist.getGridIterator(box.getKP1(),box.getKP2()); + + /////// GPU Run kernel + + gdist.setInsertBuffer(128); + + float c = 5.0; + + gdist.template iterateGridGPU>(it,c); + gdist.template flush>(flush_type::FLUSH_ON_DEVICE); + + gdist.template deviceToHost<0>(); + + { + Box<2,size_t> box2({0,0},{15,15}); + + auto it = gdist.getGridIterator(box2.getKP1(),box2.getKP2()); + + while (it.isNext()) + { + auto p = it.get_dist(); + auto p2 = it.get(); + + if (p2.get(0) == box.getLow(0) && p2.get(1) == box.getLow(1)) + { + BOOST_REQUIRE_EQUAL(gdist.template get<0>(p), 5.0); + } + else + { + BOOST_REQUIRE_EQUAL(gdist.template get<0>(p), 666.0); + } + + ++it; + } + } + + // + + c = 3.0; + + Box<2,size_t> box3({3,3},{11,11}); + + auto it3 = gdist.getGridIterator(box3.getKP1(),box3.getKP2()); + gdist.setInsertBuffer(128); + + gdist.template iterateGridGPU>(it3,c); + gdist.template flush>(flush_type::FLUSH_ON_DEVICE); + gdist.template deviceToHost<0>(); + + { + Box<2,size_t> box2({0,0},{15,15}); + + auto it = gdist.getGridIterator(box2.getKP1(),box2.getKP2()); + + while (it.isNext()) + { + auto p = it.get_dist(); + auto p2 = it.get(); + + Point<2,size_t> p2_ = p2.toPoint(); + + if (box.isInside(p2_)) + { + BOOST_REQUIRE_EQUAL(gdist.template get<0>(p), 5.0); + } + else if (box3.isInside(p2_)) + { + BOOST_REQUIRE_EQUAL(gdist.template get<0>(p), 3.0); + } + else + { + BOOST_REQUIRE_EQUAL(gdist.template get<0>(p), 666.0); + } + + ++it; + } + } + + //////////////////////////////////// + + gdist.setInsertBuffer(128); + gdist.template iterateGPU>(); + + +} + + +BOOST_AUTO_TEST_SUITE_END() diff --git a/src/Grid/tests/sgrid_dist_id_unit_tests.cpp b/src/Grid/tests/sgrid_dist_id_unit_tests.cpp new file mode 100644 index 00000000..166c235e --- /dev/null +++ b/src/Grid/tests/sgrid_dist_id_unit_tests.cpp @@ -0,0 +1,250 @@ +/* + * sgrid_dist_id_unit_tests.cpp + * + * Created on: Nov 18, 2017 + * Author: i-bird + */ + + +#define BOOST_TEST_DYN_LINK +#include +#include "Grid/grid_dist_id.hpp" +#include "Point_test.hpp" + +////////////////////////////////////// THEESE TEST ARE BROKEN TO REMPOVE OR FIX //// + + +const int x = 0; +const int y = 1; +const int z = 2; + +BOOST_AUTO_TEST_SUITE( sgrid_dist_id_test ) + +BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test_2D) +{ + periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC}; + + // Domain + Box<2,double> domain({-0.3,-0.3},{1.0,1.0}); + + // grid size + size_t sz[2]; + sz[0] = 1024; + sz[1] = 1024; + + // Ghost + Ghost<2,double> g(0.01); + + sgrid_dist_id<2,double,Point_test> sg(sz,domain,g,bc); + + // create a grid iterator + + auto it = sg.getGridIterator(); + + while(it.isNext()) + { + auto gkey = it.get(); + auto key = it.get_dist(); + + + long int sx = gkey.get(0) - 512; + long int sy = gkey.get(1) - 512; + + if (sx*sx + sy*sy < 128*128) + { + sg.template insert<0>(key) = 1.0; + } + + ++it; + } + + sg.write("sg_test_write"); + + bool match = true; + auto it2 = sg.getGridIterator(); + + while(it2.isNext()) + { + auto gkey = it2.get(); + auto key = it2.get_dist(); + + long int sx = gkey.get(0) - 512; + long int sy = gkey.get(1) - 512; + + if (sx*sx + sy*sy < 128*128) + { + match &= (sg.template get<0>(key) == 1.0); + } + + ++it2; + } + + BOOST_REQUIRE_EQUAL(match,true); + + auto & gr = sg.getGridInfo(); + + auto it3 = sg.getDomainIterator(); + + while (it3.isNext()) + { + auto key = it3.get(); + auto gkey = it3.getGKey(key); + + sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1); + + ++it3; + } + + sg.ghost_get<0>(); + + // now we check the stencil + + bool good = true; + auto it4 = sg.getDomainIterator(); + + while (it4.isNext()) + { + auto key = it4.get(); + auto gkey = it4.getGKey(key); + + double lap; + + // Here we check that all point of the stencil are inside*/ + + long int sx = gkey.get(0) - 512; + long int sy = gkey.get(1) - 512; + + if (sx*sx + sy*sy < 126*126) + { + lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) + + sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) - + 4.0*sg.template get<0>(key); + + good &= (lap == 4.0); + + if (good == false) + { + int debug = 0; + + std::cout << sg.template get<0>(key.move(x,1)) << " " << sg.template get<0>(key.move(x,-1)) << " " << + sg.template get<0>(key.move(y,1))<< " " << sg.template get<0>(key.move(y,-1)) << " " << + 4.0*sg.template get<0>(key) << std::endl; + + debug++; + } + } + + ++it4; + } + + BOOST_REQUIRE_EQUAL(good,true); +} + +BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test) +{ + periodicity<3> bc = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC}; + + // Domain + Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0}); + + // grid size + size_t sz[3]; + sz[0] = 1024; + sz[1] = 1024; + sz[2] = 1024; + + // Ghost + Ghost<3,double> g(0.01); + + sgrid_dist_id<3,double,Point_test> sg(sz,domain,g,bc); + + // create a grid iterator over a bilion point + + auto it = sg.getGridIterator(); + + while(it.isNext()) + { + auto gkey = it.get(); + auto key = it.get_dist(); + + size_t sx = gkey.get(0) - 512; + size_t sy = gkey.get(1) - 512; + size_t sz = gkey.get(2) - 512; + + if (sx*sx + sy*sy + sz*sz < 128*128) + { + sg.template insert<0>(key) = 1.0; + } + + ++it; + } + + bool match = true; + auto it2 = sg.getGridIterator(); + + while(it2.isNext()) + { + auto gkey = it2.get(); + auto key = it2.get_dist(); + + size_t sx = gkey.get(0) - 512; + size_t sy = gkey.get(1) - 512; + size_t sz = gkey.get(2) - 512; + + if (sx*sx + sy*sy + sz*sz < 128*128) + { + match &= (sg.template get<0>(key) == 1.0); + } + + ++it2; + } + + auto & gr = sg.getGridInfo(); + + auto it3 = sg.getDomainIterator(); + + while (it3.isNext()) + { + auto key = it3.get(); + auto gkey = it3.getGKey(key); + + sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1) + gkey.get(2)*gkey.get(2); + + ++it3; + } + + sg.ghost_get<0>(); + // now we check the stencil + + bool good = true; + auto it4 = sg.getDomainIterator(); + + while (it4.isNext()) + { + auto key = it4.get(); + auto gkey = it4.getGKey(key); + + size_t sx = gkey.get(0) - 512; + size_t sy = gkey.get(1) - 512; + size_t sz = gkey.get(2) - 512; + + if (sx*sx + sy*sy + sz*sz < 125*125) + { + double lap; + + lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) + + sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) + + sg.template get<0>(key.move(z,1)) + sg.template get<0>(key.move(z,-1)) - + 6.0*sg.template get<0>(key); + + good &= (lap == 6.0); + } + + ++it4; + } + + BOOST_REQUIRE_EQUAL(match,true); +} + + +BOOST_AUTO_TEST_SUITE_END() diff --git a/src/SubdomainGraphNodes.hpp b/src/SubdomainGraphNodes.hpp index a3722bde..ec60446f 100755 --- a/src/SubdomainGraphNodes.hpp +++ b/src/SubdomainGraphNodes.hpp @@ -39,7 +39,7 @@ template struct nm_v { //! The node contain 3 unsigned long integer for communication computation memory and id - typedef boost::fusion::vector type; + typedef boost::fusion::vector type; //! type of the positional field typedef float s_type; -- GitLab