From 5fdb12a668508333a2dbea2ebb159082c7c022d3 Mon Sep 17 00:00:00 2001 From: Pietro Incardona <incardon@mpi-cbg.de> Date: Thu, 20 Oct 2016 13:45:01 +0200 Subject: [PATCH] Adding missing files --- CHANGELOG.md | 27 +- example/Vector/0_simple/main.cpp | 2 +- .../4_multiphase_celllist_verlet/Makefile | 8 +- .../4_multiphase_celllist_verlet/main.cpp | 448 +++++++++++++++--- src/Grid/staggered_dist_grid_util.hpp | 22 + src/Makefile.am | 2 +- src/Vector/vector_dist_MP_unit_tests.hpp | 377 +++++++++++++++ .../vector_dist_multiphase_functions.hpp | 126 +++++ src/Vector/vector_dist_unit_test.hpp | 36 +- 9 files changed, 978 insertions(+), 70 deletions(-) create mode 100644 src/Vector/vector_dist_MP_unit_tests.hpp create mode 100644 src/Vector/vector_dist_multiphase_functions.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index bfc246e8c..7c8573692 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,12 +1,30 @@ # Change Log All notable changes to this project will be documented in this file. +## [0.6.0] - End October 2016 + +### Added +- Symmetric cell-list/verlet list +- Multi-phase cell-list and Multi-phase cell-list +- Added ghost_get that keep properties +- Examples: 1_ghost_get_put it show how to use ghost_get and put with the new options + 4_multiphase_celllist_verlet completely rewritten for new Cell-list and multiphase verlet + 5_molecular_dynamic use case of symmetric cell-list and verlet list with ghost put + 6_complex_usage It show how the flexibility of openfpm can be used to debug your program + +### Fixed +- Option NO_POSITION was untested + +### Changes + + ## [0.5.1] - 27 September 2016 ### Added - ghost_put support for particles - Full-Support for complex property on vector_dist (Serialization) - Added examples for serialization of complex properties 4_Vector +- improved speed of the iterators ### Fixed - Installation PETSC installation fail in case of preinstalled MPI @@ -131,17 +149,14 @@ All notable changes to this project will be documented in this file. - Algebraic Multigrid solver - Parallel VTK, improved visualization +- Asynchronous communication -## [0.6.0] - Middle of October +## [0.7.0] - December of October ### Added +- Asynchronous communication -- Symmetric Cell list and Verlet (15 days) -- Semantic communication (??) -- Improved Finite difference interface (15 days) -## [0.6.0] - Beginning of September -- Complex properties and serialization interface (15 days) diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp index aeba11f0a..0e6a33447 100644 --- a/example/Vector/0_simple/main.cpp +++ b/example/Vector/0_simple/main.cpp @@ -348,7 +348,7 @@ int main(int argc, char* argv[]) * * ## Finalize ## {#finalize_e0_sim} * - * At the very end of the program we have always to de-initialize the library + * At the very end of the program we have always de-initialize the library * * \snippet Vector/0_simple/main.cpp finalize * diff --git a/example/Vector/4_multiphase_celllist_verlet/Makefile b/example/Vector/4_multiphase_celllist_verlet/Makefile index 16168a458..6043bb602 100644 --- a/example/Vector/4_multiphase_celllist_verlet/Makefile +++ b/example/Vector/4_multiphase_celllist_verlet/Makefile @@ -9,16 +9,16 @@ OBJ = main.o %.o: %.cpp $(CC) -fext-numeric-literals -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH) -cell: $(OBJ) +multip: $(OBJ) $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS) -all: cell +all: multip run: all - source $$HOME/openfpm_vars; mpirun -np 2 ./cell + source $$HOME/openfpm_vars; mpirun -np 2 ./multip .PHONY: clean all run clean: - rm -f *.o *~ core cell + rm -f *.o *~ core multip diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp index c353f0d8d..bf6223ca4 100644 --- a/example/Vector/4_multiphase_celllist_verlet/main.cpp +++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp @@ -3,16 +3,18 @@ #include "Decomposition/CartDecomposition.hpp" #include "data_type/aggregate.hpp" #include "NN/CellList/CellListM.hpp" +#include "Vector/vector_dist_multiphase_functions.hpp" /*! - * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list and verlet * * [TOC] * * - * # Vector Multi Phase cell-list # {#e4_ph_cl} + * # Vector Multi Phase cell-list and Verlet # {#e4_ph_cl} * - * This example show multi-phase cell lists for the distributed vector + * This example show how to use multi-phase cell lists and Verlet-list.More in general + * it show how to construct Verlet and Cell-list between multiple vector_dist. * * */ @@ -25,19 +27,23 @@ int main(int argc, char* argv[]) * ## Initialization ## * * Here we Initialize the library, and we create a set of distributed vectors all forced to have the same - * decomposition. Each vector identify one phase + * decomposition. Each vector identify one phase. * - * \snippet Vector/4_multiphase_celllist/main.cpp Initialization and parameters + * \note Be carefull on how you initialize the other phases. All the other phases + * must be forced to use the same decomposition. In order to do this we have + * to use the special constructor where we pass the decomposition from the first + * phase. The second parameter is just the the number of particles + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp Initialization and parameters * */ //! \cond [Initialization and parameters] \endcond openfpm_init(&argc,&argv); - Vcluster & v_cl = create_vcluster(); - // we will place the particles on a grid like way with 128 particles on each direction - size_t sz[3] = {128,128,128}; + // Vcluster for general usage + Vcluster & v_cl = create_vcluster(); // The domain Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0}); @@ -45,11 +51,13 @@ int main(int argc, char* argv[]) // Boundary conditions size_t bc[3]={PERIODIC,PERIODIC,PERIODIC}; + // cut-off radius float r_cut = 0.05; // ghost, big enough to contain the interaction radius Ghost<3,float> ghost(r_cut); + // The set of phases openfpm::vector< vector_dist<3,float, aggregate<double,double>> > phases; // first phase @@ -67,108 +75,431 @@ int main(int argc, char* argv[]) /*! * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list * - * ## Initialization ## + * ## Create phases ## + * + * We initialize all the phases with particle randomly positioned in the space. Completed this + * iteration we redistribute the particles using the classical map, and we synchronize the ghost * - * We initialize all the phases with particle randomly positioned in the space + * \warning we cannot use the same iterator for all the phases even if the number of particles for + * each phase is the same. Consider that the number of particles (per-processor) can be + * different. The case of the initialization ("before map") is the only case where we are + * sure that the number of particles per processor is the same for each phase * - * \snippet Vector/4_multiphase_celllist/main.cpp rand dist + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp rand dist * */ //! \cond [rand dist] \endcond + // An iterator over the particles of the phase0 auto it = phases.get(0).getDomainIterator(); - // For all the particles + // For all the particles of phase0 while (it.isNext()) { - // for all phases - for (size_t i = 0 ; i < phases.size() ; i++) - { - auto key = it.get(); + // particle p + auto p = it.get(); - phases.get(i).getPos(key)[0] = (float)rand() / RAND_MAX; - phases.get(i).getPos(key)[1] = (float)rand() / RAND_MAX; - phases.get(i).getPos(key)[2] = (float)rand() / RAND_MAX; + // Assign the position of the particles to each phase + for (size_t i = 0 ; i < phases.size(); i++) + { + // Assign random position + phases.get(i).getPos(p)[0] = (float)rand() / RAND_MAX; + phases.get(i).getPos(p)[1] = (float)rand() / RAND_MAX; + phases.get(i).getPos(p)[2] = (float)rand() / RAND_MAX; } - // next point + + // Next particle ++it; } - // Redistribute all the phases in the mean while also take the iterators of all the phases - typedef decltype(phases.get(0).getIterator()) iterator; - openfpm::vector<iterator> phase_it; - - for (size_t i = 0 ; i < phases.size() ; i++) + // Redistribute and sync all the phases + for (size_t i = 0 ; i < 4 ; i++) + { phases.get(i).map(); + phases.get(i).ghost_get<>(); + } //! \cond [rand dist] \endcond /*! * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list * - * ## Multi-phase cell-list construction ## + * ## Multi phase Verlet-list ## + * + * In general verlet-list can be constructed from the vector itself using **getVerlerList()** + * + * \see \ref Vector_3_md_vl + * + * In the multi-phase case if we use such function on phase0 such function + * produce a Verlet list where for each particle of the phase0 we get the neighborhood + * particles within the phase0. Most of time we also want to construct Verlet-list across + * phases, for example between phase0 and phase1. Let suppose now that we want to construct + * a verlet-list that given the particles of the phase0 it return the neighborhood particles + * in phase1. In order to do this we can create a Cell-list from phase1. Once we have the + * Cell-list for phase 1 we give to **createVerlet()** + * the particle set of phase0, the Cell-list of phase1 and the cut-off radius. + * The function return a VerletList between phase0 and 1 that can be used to do computation. + * + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create multi-phase verlet + * + */ + + //! \cond [create multi-phase verlet] \endcond + + // Get the cell list of the phase1 + auto CL_phase1 = phases.get(1).getCellList(r_cut); + + // This function create a Verlet-list between phases 0 and 1 + auto NN_ver01 = createVerlet(phases.get(0),CL_phase1,r_cut); + + //! \cond [create multi-phase verlet] \endcond + + /*! + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * + * Once we have a Verlet-list, we can do easily computation with that. We can use + * the function **getNNIterator()** to get an iterator of the neighborhood particles + * for a specified particle. In this case for each particle of the phase0 we count + * the neighborhood particles in phase1 + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp count part from phase0 to 1 + * + */ + + //! \cond [count part from phase0 to 1] \endcond + + // Get an iterator of the particles of the phase0 + auto it = phases.get(0).getDomainIterator(); + + // For each particle of the phase0 + while (it.isNext()) + { + // Get the particle p + auto p = it.get(); + + // reset the counter + phases.get(0).getProp<0>(p) = 0; + + // Get an iterator over the neighborhood particles for the particle p + auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey()); + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Neighborhood particle q + auto q = Np.get(); + + // Count the number of particles + // Here in general we can do our computation + phases.get(0).getProp<0>(p)++; + + // Next particle + ++Np; + } + + // Next particle + ++it; + } + + //! \cond [count part from phase0 to 1] \endcond + + /*! + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * + * ## From one phase to all ## + * + * It is also possible to construct a Verlet-list from one phase to all the other phases. + * To do this we use the function **createCellListM<2>()** to first create a + * multi-phase cell-list. The template parameter is required to indicate how many bit + * to reserve for the phase information. + * + * \note In case of + * * 2 bit mean that you can store up to 4 phases (2^62 number of particles for each phase) + * * 3 bit mean that you can store up to 8 phases (2^61 number of particles for each phase) * - * In this part we construct the Multi-phase cell list. The multiphase cell list has 3 parameters - * * one is the dimensionality (3) - * * The precision of the coordinates (float), - * * How many bit to use for the phase. - * Multi-phase cell-list try to pack into a 64bit number information about the particle id and the - * the phase id. + * Once created the multiphase-cell list from the phases. We can use the function + * **createVerletM()** to create a verlet-list from one phase to all the others. + * The function requires the particle for which we are constructing the verlet list, in this case + * the phase0, the Cell-list containing all the other phases and the cut-off radius. * - * \snippet Vector/4_multiphase_celllist/main.cpp cl construction + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create multi-phase multi verlet * */ - //! \cond [cl construction] \endcond + //! \cond [create multi-phase multi verlet] \endcond + + // This function create an "Empty" Multiphase Cell List from all the phases + auto CL_all = createCellListM<2>(phases,r_cut); + + // This create a Verlet-list between phase0 and all the other phases + auto NNver0_all = createVerletM<2>(phases.get(0),CL_all,r_cut); - // Construct one single Multi-phase cell list to use in the computation - // in 3d, precision float, 2 bit dedicated to the phase for a maximum of 2^2 = 4 (Maximum number of phase) - // - // - - size_t div[3]; - Box<3,float> box_cl; - phases.get(0).getCellListParams(r_cut,div,box_cl); + //! \cond [create multi-phase multi verlet] \endcond - CellListM<3,float,2> NN; - NN.Initialize(box_cl,div); + /*! + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * + * Compute on a multiphase verlet-list is very similar to a non multi-phase. + * The only difference is that the function **get()** is substituted by two + * other functions **getP()** and **getV()** The first return the phase from where + * it come the particle the second it return the particle-id + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp compute multi-phase multi verlet + * + */ + //! \cond [compute multi-phase multi verlet] \endcond + // Get an iterator for the phase0 particles + it = phases.get(0).getDomainIterator(); - // for all the phases i - for (size_t i = 0; i < phases.size() ; i++) + while (it.isNext()) { - // iterate across all the particle of the phase i - auto it = phases.get(i).getDomainIterator(); + // Get the particle p + auto p = it.get(); + + // Get an interator over the neighborhood of the particle p + auto Np = NNver0_all.template getNNIterator<NO_CHECK>(p.getKey()); + + // reset the counter + phases.get(0).getProp<0>(p) = 0; + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Get the particle q near to p + auto q = Np.getP(); + + // Get from which phase it come from + auto ph_q = Np.getV(); + + // count + phases.get(0).getProp<0>(p)++; + + // Next particle + ++Np; + } + + // Next particle + ++it; + } + + //! \cond [compute multi-phase multi verlet] \endcond + + /*! + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * + * ## Symmetric interaction case ## + * + * The same functions exist also in the case we want to construct + * symmetric-verlet list. + * + * \see \ref Vector_5_md_vl_sym For more details on how to use symmetric + * verlet-list in a real case. + * + * In general the main differences can be summarized in + * * we have to reset the **forces** or **interaction** variable(s) + * * calculate the interaction p-q and apply the result to p and q at the same time + * * merge with ghost_put the result is stored in the ghost part with the + * real particles + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp compute sym multi-phase two phase + * + */ + + //! \cond [compute sym multi-phase two phase] \endcond + + // Get the cell list of the phase1 + auto CL_phase1 = phases.get(1).getCellListSym(r_cut); + + // This function create a Verlet-list between phases 0 and 1 + auto NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut); + + // Get an iterator over the real and ghost particles + auto it = phases.get(0).getDomainAndGhostIterator(); + + // For each particles + while (it.isNext()) + { + // Get the particle p + auto p = it.get(); + + // Reset the counter + phases.get(0).getProp<0>(p) = 0; + + // Next particle + ++it; + } + + // Compute interaction from phase0 to phase1 + + // Get an iterator over the real particles of phase0 + auto it = phases.get(0).getDomainIterator(); + + // For each particle of the phase0 + while (it.isNext()) + { + // get particle p + auto p = it.get(); + + // Get the neighborhood of particle p + auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey()); + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Neighborhood particle q of p + auto q = Np.get(); + + // increment the counter for the phase 0 and 1 + phases.get(0).getProp<0>(p)++; + phases.get(1).getProp<0>(q)++; + + // Next particle + ++Np; + } + + // Next particle + ++it; + } + + // Merge the information of the ghost with the real particles + phases.get(0).ghost_put<add_,0>(); + phases.get(1).ghost_put<add_,0>(); + + //! \cond [compute sym multi-phase two phase] \endcond + + /*! + * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list + * + * For the case of a Verlet-list between the phase0 and all the other phases. + * The code remain very similar to the previous one the only differences is that + * we have to create a Multi-phase cell list from the vector of the phases. + * When we compute with the verlet-list list now the simple function **get** + * is substituted by the function **getP** and **getV**. In this example we are + * creating 4 multiphase symmetric verlet-list + * + * * 0 to all + * * 1 to all + * * 2 to all + * * 3 to all + * + * The computation is an all to all + * + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create sym multi-phase multi verlet + * + */ + + //! \cond [create sym multi-phase multi verlet] \endcond + + // Get an iterator over the phase0 + auto it = phases.get(0).getDomainAndGhostIterator(); + + // For each particle of the phase 0 + while (it.isNext()) + { + // get the particle p + auto p = it.get(); + + // Reset the counter for the domain and ghost particles + phases.get(0).getProp<0>(p) = 0; + phases.get(1).getProp<0>(p) = 0; + phases.get(2).getProp<0>(p) = 0; + phases.get(3).getProp<0>(p) = 0; + + // next particle + ++it; + } + + // This function create an "Empty" Multiphase Cell List + auto CL_all = createCellListSymM<2>(phases,r_cut); + + // Type of the multiphase Verlet-list + typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type; + + // for each phase we create one Verlet-list that contain the neighborhood + // from all the phases + verlet_type NNver_all[4]; + + // Here we create a Verlet-list between each phases + + // 0 to all + NNver_all[0] = createVerletSymM<2>(phases.get(0),CL_all,r_cut); + + // 1 to all + NNver_all[1] = createVerletSymM<2>(phases.get(1),CL_all,r_cut); + + // 2 to all + NNver_all[2] = createVerletSymM<2>(phases.get(2),CL_all,r_cut); + + // 3 to all + NNver_all[3] = createVerletSymM<2>(phases.get(3),CL_all,r_cut); + + // For each phase + for (size_t i = 0 ; i < phases.size() ; i++) + { + // Get an iterator over the particles of that phase + it = phases.get(i).getDomainIterator(); + + // for each particle of the phase while (it.isNext()) { - auto key = it.get(); + // Get the particle p + auto p = it.get(); + + // Get an iterator for neighborhood of the particle p + auto Np = NNver_all[i].template getNNIterator<NO_CHECK>(p.getKey()); + + // For each neighborhood particle + while (Np.isNext()) + { + // Get the particle q near to p + auto q = Np.getP(); - // Add the particle of the phase i to the cell list - NN.add(phases.get(i).getPos(key), key.getKey(), i); + // Get from which phase it come from + auto ph_q = Np.getV(); + // increment the counter on both p and q + phases.get(i).getProp<0>(p)++; + phases.get(ph_q).getProp<0>(q)++; + + // Next particle + ++Np; + } + + // Next particle ++it; } } - //! \cond [cl construction] \endcond + // Merge the ghost part with the real particles + phases.get(0).ghost_put<add_,0>(); + phases.get(1).ghost_put<add_,0>(); + phases.get(2).ghost_put<add_,0>(); + phases.get(3).ghost_put<add_,0>(); + + //! \cond [create sym multi-phase multi verlet] \endcond /*! * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list * * ## Multi-phase cell-list usage ## * - * After construction we show how to use the Cell-list. In this case we accumulate on the property + * The multiphase cell-list that we constructed before to create a Verlet-list can be also used + * directly. In this case we accumulate on the property * 0 of the phase 0 the distance of the near particles from all the phases * - * \snippet Vector/4_multiphase_celllist/main.cpp cl usage + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp cl usage * */ //! \cond [cl usage] \endcond + // we create a reference to phase0 particle for convenience vector_dist<3,float, aggregate<double,double> > & current_phase = phases.get(0); // Get the iterator of the particles of phase 0 @@ -183,8 +514,11 @@ int main(int argc, char* argv[]) // Get the position of the particle p Point<3,float> xp = current_phase.getPos(p); + // Reset to zero the propety 0 + current_phase.getProp<0>(p) = 0.0; + // Get an iterator of all the particles neighborhood of p - auto Np = NN.getNNIterator(NN.getCell(current_phase.getPos(p))); + auto Np = CL_all.getNNIterator(NN.getCell(current_phase.getPos(p))); // For each particle near p while (Np.isNext()) @@ -216,7 +550,7 @@ int main(int argc, char* argv[]) * * At the very end of the program we have always de-initialize the library * - * \snippet Vector/4_multiphase_celllist/main.cpp finalize + * \snippet Vector/4_multiphase_celllist_verlet/main.cpp finalize * */ @@ -231,7 +565,7 @@ int main(int argc, char* argv[]) * * # Full code # {#code} * - * \include Vector/4_multiphase_celllist/main.cpp + * \include Vector/4_multiphase_celllist_verlet/main.cpp * */ } diff --git a/src/Grid/staggered_dist_grid_util.hpp b/src/Grid/staggered_dist_grid_util.hpp index c7d20ec6a..938f1d9ab 100644 --- a/src/Grid/staggered_dist_grid_util.hpp +++ b/src/Grid/staggered_dist_grid_util.hpp @@ -66,11 +66,13 @@ struct vtk_write<ele,vtk,false> template<typename T> struct extends { + //! number of elements static inline size_t mul() { return 1; } + //! number of indexes static inline size_t dim() { return 0; @@ -81,11 +83,13 @@ struct extends template<typename T,size_t N1> struct extends<T[N1]> { + //! number of elements static inline size_t mul() { return N1; } + //! number of indexes static inline size_t dim() { return 1; @@ -96,11 +100,13 @@ struct extends<T[N1]> template<typename T,size_t N1,size_t N2> struct extends<T[N1][N2]> { + //! number of elements static inline size_t mul() { return N1 * N2; } + //! number of indexes static inline size_t dim() { return 2; @@ -111,11 +117,13 @@ struct extends<T[N1][N2]> template<typename T,size_t N1,size_t N2,size_t N3> struct extends<T[N1][N2][N3]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3; } + //! number of indexes static inline size_t dim() { return 3; @@ -126,11 +134,13 @@ struct extends<T[N1][N2][N3]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4> struct extends<T[N1][N2][N3][N4]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4; } + //! number of indexes static inline size_t dim() { return 4; @@ -141,11 +151,13 @@ struct extends<T[N1][N2][N3][N4]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5> struct extends<T[N1][N2][N3][N4][N5]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5; } + //! number of indexes static inline size_t dim() { return 5; @@ -156,11 +168,13 @@ struct extends<T[N1][N2][N3][N4][N5]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6> struct extends<T[N1][N2][N3][N4][N5][N6]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5 * N6; } + //! number of indexes static inline size_t dim() { return 6; @@ -171,11 +185,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7> struct extends<T[N1][N2][N3][N4][N5][N6][N7]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5 * N6 * N7; } + //! number of indexes static inline size_t dim() { return 7; @@ -186,11 +202,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8> struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8; } + //! number of indexes static inline size_t dim() { return 8; @@ -201,11 +219,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9> struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9; } + //! number of indexes static inline size_t dim() { return 9; @@ -216,11 +236,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]> template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9, size_t N10> struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9][N10]> { + //! number of elements static inline size_t mul() { return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9 * N10; } + //! number of indexes static inline size_t dim() { return 10; diff --git a/src/Makefile.am b/src/Makefile.am index fff4127b0..3ba89adcd 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -9,7 +9,7 @@ nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartD Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \ Graph/CartesianGraphFactory.hpp \ Grid/grid_dist_id.hpp Grid/grid_dist_id_iterator_dec.hpp Grid/grid_dist_util.hpp Grid/grid_dist_id_iterator_sub.hpp Grid/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \ - Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \ + Vector/vector_dist_multiphase_functions.hpp Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \ config/config.h \ example.mk \ Decomposition/Distribution/metis_util.hpp Decomposition/Distribution/parmetis_dist_util.hpp Decomposition/Distribution/parmetis_util.hpp Decomposition/Distribution/MetisDistribution.hpp Decomposition/Distribution/ParMetisDistribution.hpp Decomposition/Distribution/DistParMetisDistribution.hpp dec_optimizer.hpp SubdomainGraphNodes.hpp \ diff --git a/src/Vector/vector_dist_MP_unit_tests.hpp b/src/Vector/vector_dist_MP_unit_tests.hpp new file mode 100644 index 000000000..7d9eff81e --- /dev/null +++ b/src/Vector/vector_dist_MP_unit_tests.hpp @@ -0,0 +1,377 @@ +/* + * vector_dist_MP_unit_tests.hpp + * + * Created on: Oct 14, 2016 + * Author: i-bird + */ + +#ifndef SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_ +#define SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_ + +#include "Vector/vector_dist_multiphase_functions.hpp" + +BOOST_AUTO_TEST_SUITE( vector_dist_multiphase_test ) + +BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_test ) +{ + if (create_vcluster().getProcessingUnits() > 24) + return; + + size_t sz[3] = {60,60,40}; + + // The domain + Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0}); + + // Boundary conditions + size_t bc[3]={PERIODIC,PERIODIC,PERIODIC}; + + float r_cut = 51.0; + + // ghost, big enough to contain the interaction radius + Ghost<3,float> ghost(r_cut); + + openfpm::vector< vector_dist<3,float, aggregate<double,double>> > phases; + + // first phase + phases.add( vector_dist<3,float, aggregate<double,double>>(0,box,bc,ghost) ); + + // The other 3 phases + phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) ); + phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) ); + phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) ); + + // Fill the phases with particles + + auto g_it = phases.get(0).getGridIterator(sz); + + while (g_it.isNext()) + { + auto key = g_it.get(); + + // Add a particle to all the phases + phases.get(0).add(); + phases.get(1).add(); + phases.get(2).add(); + phases.get(3).add(); + + phases.get(0).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(0).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(0).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(1).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(1).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(1).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(2).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(2).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(2).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(3).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(3).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(3).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + ++g_it; + } + + // Sync all phases + for (size_t i = 0 ; i < 4 ; i++) + { + phases.get(i).map(); + phases.get(i).ghost_get<>(); + } + + // Get the cell list of the phase 0 and 1 + auto CL_phase0 = phases.get(0).getCellList(r_cut); + auto CL_phase1 = phases.get(1).getCellList(r_cut); + + // This function create a Verlet-list between phases 0 and 1 + auto NN_ver01 = createVerlet(phases.get(0),CL_phase1,r_cut); + + // Check NNver0_1 + + bool ret = true; + auto it = phases.get(0).getDomainIterator(); + + while (it.isNext()) + { + auto p = it.get(); + auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey()); + + size_t nn_count = 0; + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Neighborhood particle q + auto q = Np.get(); + + // Count the number of particles + nn_count++; + + ++Np; + } + + ret &= nn_count == 7ul; + + ++it; + } + + BOOST_REQUIRE_EQUAL(ret,true); + + // Sync all phases + for (size_t i = 0 ; i < 4 ; i++) + { + phases.get(i).map(); + phases.get(i).ghost_get<>(); + } + + // NN_ver0_all + + // This function create an "Empty" Multiphase Cell List + auto CL_all = createCellListM<2>(phases,r_cut); + + // This create a Verlet-list between phase 0 and all the other phases + auto NNver0_all = createVerletM<2>(phases.get(0),CL_all,r_cut); + + it = phases.get(0).getDomainIterator(); + + while (it.isNext()) + { + auto p = it.get(); + auto Np = NNver0_all.template getNNIterator<NO_CHECK>(p.getKey()); + + size_t nn_cout[4] = {0,0,0,0}; + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Get the particle q near to p + auto q = Np.getP(); + + // Get from which phase it come from + auto ph_q = Np.getV(); + + nn_cout[ph_q]++; + + ++Np; + } + + ret &= nn_cout[0] == 7; + ret &= nn_cout[1] == 7; + ret &= nn_cout[2] == 7; + ret &= nn_cout[3] == 7; + + ++it; + } + + BOOST_REQUIRE_EQUAL(ret,true); +} + + +BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test ) +{ + if (create_vcluster().getProcessingUnits() > 24) + return; + + size_t sz[3] = {60,60,40}; + + // The domain + Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0}); + + // Boundary conditions + size_t bc[3]={PERIODIC,PERIODIC,PERIODIC}; + + float r_cut = 51.0; + + // ghost, big enough to contain the interaction radius + Ghost<3,float> ghost(r_cut); + + openfpm::vector< vector_dist<3,float, aggregate<size_t>> > phases; + + // first phase + phases.add( vector_dist<3,float, aggregate<size_t>>(0,box,bc,ghost) ); + + // The other 3 phases + phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) ); + phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) ); + phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) ); + + // Fill the phases with particles + + auto g_it = phases.get(0).getGridIterator(sz); + + while (g_it.isNext()) + { + auto key = g_it.get(); + + // Add a particle to all the phases + phases.get(0).add(); + phases.get(1).add(); + phases.get(2).add(); + phases.get(3).add(); + + phases.get(0).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(0).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(0).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(1).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(1).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(1).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(2).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(2).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(2).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + phases.get(3).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0); + phases.get(3).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1); + phases.get(3).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2); + + ++g_it; + } + + // Sync all phases + for (size_t i = 0 ; i < 4 ; i++) + { + phases.get(i).map(); + phases.get(i).ghost_get<>(); + } + + // Get the cell list of the phase 0 and 1 + auto CL_phase0 = phases.get(0).getCellListSym(r_cut); + auto CL_phase1 = phases.get(1).getCellListSym(r_cut); + + // This function create a Verlet-list between phases 0 and 1 + auto NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut); + + // Check NNver0_1 + + bool ret = true; + auto it = phases.get(0).getDomainIterator(); + + while (it.isNext()) + { + auto p = it.get(); + auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey()); + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Neighborhood particle q + auto q = Np.get(); + + phases.get(0).getProp<0>(p)++; + phases.get(1).getProp<0>(q)++; + + ++Np; + } + + ++it; + } + + phases.get(0).ghost_put<add_,0>(); + phases.get(1).ghost_put<add_,0>(); + + it = phases.get(0).getDomainIterator(); + while (it.isNext()) + { + auto p = it.get(); + + ret &= phases.get(0).getProp<0>(p) == 7; + ret &= phases.get(1).getProp<0>(p) == 7; + + ++it; + } + + BOOST_REQUIRE_EQUAL(ret,true); + + // Sync all phases + for (size_t i = 0 ; i < 4 ; i++) + { + phases.get(i).map(); + phases.get(i).ghost_get<>(); + } + + // Reset counter on all phases + + for (size_t i = 0 ; i < phases.size() ; i++) + { + it = phases.get(i).getDomainAndGhostIterator(); + while (it.isNext()) + { + auto p = it.get(); + + phases.get(i).getProp<0>(p) = 0; + + ++it; + } + } + + // NN_ver0_all + + // This function create an "Empty" Multiphase Cell List + auto CL_all = createCellListSymM<2>(phases,r_cut); + + typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type; + + verlet_type NNver_all[4]; + + // This create a Verlet-list between phase all phases to all the other phases + NNver_all[0] = createVerletSymM<2>(phases.get(0),CL_all,r_cut); + NNver_all[1] = createVerletSymM<2>(phases.get(1),CL_all,r_cut); + NNver_all[2] = createVerletSymM<2>(phases.get(2),CL_all,r_cut); + NNver_all[3] = createVerletSymM<2>(phases.get(3),CL_all,r_cut); + + // all phases to all phases + + for (size_t i = 0 ; i < phases.size() ; i++) + { + it = phases.get(i).getDomainIterator(); + + while (it.isNext()) + { + auto p = it.get(); + auto Np = NNver_all[i].template getNNIterator<NO_CHECK>(p.getKey()); + + // For each neighborhood of the particle p + while (Np.isNext()) + { + // Get the particle q near to p + auto q = Np.getP(); + + // Get from which phase it come from + auto ph_q = Np.getV(); + + phases.get(i).getProp<0>(p)++; + phases.get(ph_q).getProp<0>(q)++; + + ++Np; + } + + ++it; + } + } + + phases.get(0).ghost_put<add_,0>(); + phases.get(1).ghost_put<add_,0>(); + phases.get(2).ghost_put<add_,0>(); + phases.get(3).ghost_put<add_,0>(); + + it = phases.get(0).getDomainIterator(); + while (it.isNext()) + { + auto p = it.get(); + + ret &= phases.get(0).getProp<0>(p) == 32; + ret &= phases.get(1).getProp<0>(p) == 32; + ret &= phases.get(2).getProp<0>(p) == 32; + ret &= phases.get(3).getProp<0>(p) == 32; + + ++it; + } + + BOOST_REQUIRE_EQUAL(ret,true); +} + +BOOST_AUTO_TEST_SUITE_END() + +#endif /* SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_ */ diff --git a/src/Vector/vector_dist_multiphase_functions.hpp b/src/Vector/vector_dist_multiphase_functions.hpp new file mode 100644 index 000000000..4a3c6e32b --- /dev/null +++ b/src/Vector/vector_dist_multiphase_functions.hpp @@ -0,0 +1,126 @@ +/* + * vector_dist_multiphase_functions.hpp + * + * Created on: Oct 14, 2016 + * Author: i-bird + */ + +#ifndef SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_ +#define SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_ + +#include "NN/CellList/CellListM.hpp" +#include "NN/VerletList/VerletListM.hpp" + +template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerlet(Vector & v, CL & cl, T r_cut) +{ + VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver; + + ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local()); + + return ver; +} + +template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletM(Vector & v, CL & cl, T r_cut) +{ + VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver; + + ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local()); + + return ver; +} + +template<unsigned int nbit, typename Vector, typename T> CellListM<Vector::dims,typename Vector::stype,nbit> createCellListM(openfpm::vector<Vector> & phases, T r_cut) +{ + size_t div[3]; + Box<Vector::dims,typename Vector::stype> box_cl; + + CellListM<Vector::dims,typename Vector::stype,nbit> NN; + if (phases.size() == 0) + return NN; + + box_cl = phases.get(0).getDecomposition().getProcessorBounds(); + phases.get(0).getCellListParams(r_cut,div,box_cl); + + NN.Initialize(box_cl,div); + + // for all the phases i + for (size_t i = 0; i < phases.size() ; i++) + { + // iterate across all the particle of the phase i + auto it = phases.get(i).getDomainAndGhostIterator(); + while (it.isNext()) + { + auto key = it.get(); + + // Add the particle of the phase i to the cell list + NN.add(phases.get(i).getPos(key), key.getKey(), i); + + ++it; + } + } + + return NN; +} + + +/////// Symmetric version + +template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerletSym(Vector & v, CL & cl, T r_cut) +{ + VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver; + + ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local()); + + return ver; +} + +template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletSymM(Vector & v, CL & cl, T r_cut) +{ + VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver; + + ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local(),VL_SYMMETRIC); + + return ver; +} + +template<unsigned int nbit, typename Vector, typename T> CellListM<Vector::dims,typename Vector::stype,nbit> createCellListSymM(openfpm::vector<Vector> & phases, T r_cut) +{ + Box<Vector::dims,typename Vector::stype> box_cl; + + CellListM<Vector::dims,typename Vector::stype,nbit> NN; + if (phases.size() == 0) + return NN; + + // Calculate the Cell list division for this CellList + CellDecomposer_sm<Vector::dims,typename Vector::stype,shift<Vector::dims,typename Vector::stype>> cd_sm; + + size_t pad = 0; + cl_param_calculateSym(phases.get(0).getDecomposition().getDomain(),cd_sm,phases.get(0).getDecomposition().getGhost(),r_cut,pad); + + // Processor bounding box + Box<Vector::dims, typename Vector::stype> pbox = phases.get(0).getDecomposition().getProcessorBounds(); + + // Ghost padding extension + Ghost<Vector::dims,size_t> g_ext(0); + NN.Initialize(cd_sm,pbox,pad); + + // for all the phases i + for (size_t i = 0; i < phases.size() ; i++) + { + // iterate across all the particle of the phase i + auto it = phases.get(i).getDomainAndGhostIterator(); + while (it.isNext()) + { + auto key = it.get(); + + // Add the particle of the phase i to the cell list + NN.add(phases.get(i).getPos(key), key.getKey(), i); + + ++it; + } + } + + return NN; +} + +#endif /* SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_ */ diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp index 8c3c800e3..866868b36 100644 --- a/src/Vector/vector_dist_unit_test.hpp +++ b/src/Vector/vector_dist_unit_test.hpp @@ -1523,6 +1523,24 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering ) vd.map(); vd.ghost_get<0,1,2>(); + // shift the particle position by 1.0 + + it = vd.getGhostIterator(); + while (it.isNext()) + { + // Particle p + auto p = it.get(); + + // we shift down he particles + vd.getPos(p)[0] = 10.0; + + // we shift + vd.getPos(p)[1] = 17.0; + + // next particle + ++it; + } + for (size_t i = 0 ; i < 10 ; i++) { auto it = vd.getDomainIterator(); @@ -1543,7 +1561,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering ) ++it; } - vd.ghost_get<0>(SKIP_LABELLING); + vd.ghost_get<0>(SKIP_LABELLING | NO_POSITION); auto it2 = vd.getGhostIterator(); bool ret = true; @@ -1560,6 +1578,22 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering ) ++it2; } + it2 = vd.getGhostIterator(); + while (it2.isNext()) + { + // Particle p + auto p = it.get(); + + // we shift down he particles + ret &= vd.getPos(p)[0] == 10.0; + + // we shift + ret &= vd.getPos(p)[1] == 17.0; + + // next particle + ++it2; + } + BOOST_REQUIRE_EQUAL(ret,true); } } -- GitLab