diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8c03e493a68d63a98059df51670da5ad03a9e82c..7c85736920c4e3e09a0ac5cbba09e11771300269 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.5.1] - Mid september
+## [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 54b4fccc64e8b3cb1f00dae622b9f67ac3d363d9..0e6a334477539f0632710a74a26b19a9236587e4 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -2,11 +2,14 @@
*
* \subpage Vector_0_simple
* \subpage Vector_1_celllist
+ * \subpage Vector_1_ghost_get
* \subpage Vector_2_expression
* \subpage Vector_3_md
* \subpage Vector_4_reo_root
* \subpage Vector_4_cp
* \subpage Vector_4_mp_cl
+ * \subpage Vector_5_md_vl_sym
+ * \subpage Vector_6_complex_usage
*
*/
@@ -37,7 +40,7 @@
* </div>
* \endhtmlonly
*
- * ## inclusion ## {#inclusion}
+ * ## inclusion ## {#e0_v_inclusion}
*
* In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
*
@@ -60,7 +63,7 @@ int main(int argc, char* argv[])
* Here we
* * Initialize the library
* * we create a Box that define our domain
- * * An array that define out boundary conditions
+ * * An array that define our boundary conditions
* * A Ghost object that will define the extension of the ghost part in physical units
*
*
@@ -345,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/1_celllist/main.cpp b/example/Vector/1_celllist/main.cpp
index 30cbb8945a25c9ac5b6c16add4242da78b4d1cc4..96aa95f7463ff4e30d7369746c8c17827e8130cb 100644
--- a/example/Vector/1_celllist/main.cpp
+++ b/example/Vector/1_celllist/main.cpp
@@ -173,6 +173,9 @@ int main(int argc, char* argv[])
* no properties. If we want to synchronize also properties we have to specify which one.
* For example with <0,1,2> here we synchronize the scalar property (0), the vector (1), and the rank 2 tensor (2)
*
+ * \warning Every ghost_get by default reset the status of the ghost so the information are lost
+ *
+ * \see \ref e1_gg_options
*
* \htmlonly
* <a href="#" onclick="hide_show('vector-video-5')" >Video 1</a>
diff --git a/example/Vector/3_molecular_dynamic/Makefile b/example/Vector/3_molecular_dynamic/Makefile
index 6d55806ec7a5643f4e95e499a145b277d5639604..d535a223f7b46c91e7f335459f98299debb699bb 100644
--- a/example/Vector/3_molecular_dynamic/Makefile
+++ b/example/Vector/3_molecular_dynamic/Makefile
@@ -7,9 +7,8 @@ LDIR =
OBJ = main.o
OBJ_EXPR = main_expr.o
OBJ_VL = main_vl.o
-OBJ_VL_SYM = main_vl_sym.o
-all: md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
+all: md_dyn md_dyn_expr md_dyn_vl
%.o: %.cpp
$(CC) -fext-numeric-literals -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
@@ -23,14 +22,11 @@ md_dyn_expr: $(OBJ_EXPR)
md_dyn_vl: $(OBJ_VL)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
-md_dyn_vl_sym: $(OBJ_VL_SYM)
- $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
-
run: all
- source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr; mpirun -np 3 ./md_dyn_vl; mpirun -np 3 ./md_dyn_vl_sym;
+ source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr; mpirun -np 3 ./md_dyn_vl;
.PHONY: clean all run
clean:
- rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
+ rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl
diff --git a/example/Vector/3_molecular_dynamic/main_vl_sym.cpp b/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
deleted file mode 100644
index 07136e19c202c479ad2559567ce58b1ea5fb419f..0000000000000000000000000000000000000000
--- a/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
+++ /dev/null
@@ -1,552 +0,0 @@
-#include "Vector/vector_dist.hpp"
-#include "Decomposition/CartDecomposition.hpp"
-#include "data_type/aggregate.hpp"
-#include "Plot/GoogleChart.hpp"
-#include "Plot/util.hpp"
-#include "timer.hpp"
-
-/*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
- *
- *
- * ## Variations for symmetric interaction ##
- *
- * In this example we show the variation in case of symmetric interations.
- *
- */
-
-//! \cond [constants] \endcond
-
-constexpr int velocity = 0;
-constexpr int force = 1;
-
-//! \cond [constants] \endcond
-
-/*!
- *
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
- *
- * ## Calculate forces ##
- *
- * In this function we calculate the forces between particles. It require the vector of particles
- * Cell list and sigma factor for the Lennard-Jhones potential. The function is exactly the same
- * as the original with the following changes
- *
- * \see \ref e3_md_vl_cf
- *
- * * We have to reset for force counter for each particles
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp reset
- *
- * * The iterator go through real AND ghost.
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp real and ghost
- *
- * * If we calculate the force for p-q we are also adding this force to q-p
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp add to q
- *
- * ** Explanation**
- *
- * The reason why symmetric Verlet exist is to reduce the computation. In a symmetric interaction
- * we have that the interaction p-q is is the same as q-p. This mean that we are calculating 2 times
- * the same quantity one when we are on particle p and one when we are on particle q. Symmetric verlet
- * avoid this. Every time we calculate a force p-q we also increment the force q-p.
- *
- * The reason instead why we have to go through ghost is the following
- *
- *
- \verbatim
-
- | |
- | Ghost|
- | |
- | |
- | |
- | |
- | |
- | +----+----+----+
- | | | | | |
- | | 1| | 2 | 3 |
- | | | | | |
- | +--------------+
- | | | * | 4 |
- | +--------------------
- | X | 5 | |
- | +---------+
- +---------------------------
-
- \endverbatim
- *
- * Where * is a particle and X is another particle on ghost
- *
- * The symmetric verlet list is constructed from the cell-List 1,2,3,4,5.
- * As you can see the interaction X-* is not present in the neighborhood of *
- * (cells 1,2,3,4,5). But is present if we draw the neighborhood of X. This
- * mean we have to iterate also across X if we want the interation X-*
- *
- *
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc forces vl
- *
- * The symmetric force calculation is faster than the normal one. but is not 2
- * time faster. The reason can be founded in 2 factors.
- *
- * * We have to reset the force
- *
- * * Every particle in the ghost
- * that interact with one in the domain is still calculated twice by the same
- * processor or some other.
- *
- * * The access pattern is worst we alternate read and write stressing more the Memory
- *
- * * We have to iterate across domain + ghost
- *
- * Overall the function calc_force is faster
- *
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc forces vl
- *
- */
-
-//! \cond [calc forces vl] \endcond
-
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
-{
- //! \cond [reset] \endcond
-
- // Reset the force for all particles
- auto it3 = vd.getDomainIterator();
-
- // For each particle p ...
- while (it3.isNext())
- {
- // ... get the particle p
- auto p = it3.get();
-
- // Reset the forice counter
- vd.template getProp<force>(p)[0] = 0.0;
- vd.template getProp<force>(p)[1] = 0.0;
- vd.template getProp<force>(p)[2] = 0.0;
-
- // Next particle p
- ++it3;
- }
-
- //! \cond [reset] \endcond
-
- //! \cond [real and ghost] \endcond
-
- // Get an iterator over particles
- auto it2 = vd.getDomainAndGhostIterator();
-
- //! \cond [real and ghost] \endcond
-
- // For each particle p ...
- while (it2.isNext())
- {
- // ... get the particle p
- auto p = it2.get();
-
- // Get the position xp of the particle
- Point<3,double> xp = vd.getPos(p);
-
- // Get an iterator over the neighborhood particles of p
- // Note that in case of symmetric
- auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
-
- // For each neighborhood particle ...
- while (Np.isNext())
- {
- // ... q
- auto q = Np.get();
-
- // if (p == q) skip this particle
- if (q == p.getKey() || (p.getKey() >= vd.size_local() && q >= vd.size_local())) {++Np; continue;};
-
- // Get the position of p
- Point<3,double> xq = vd.getPos(q);
-
- // Get the distance between p and q
- Point<3,double> r = xp - xq;
-
- // take the norm of this vector
- double rn = norm2(r);
-
- if (rn > r_cut * r_cut) {++Np;continue;}
-
- // Calculate the force, using pow is slower
- Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
-
- // we sum the force produced by q on p
- vd.template getProp<force>(p)[0] += f.get(0);
- vd.template getProp<force>(p)[1] += f.get(1);
- vd.template getProp<force>(p)[2] += f.get(2);
-
- //! \cond [add to q] \endcond
-
- // we sum the force produced by p on q
- vd.template getProp<force>(q)[0] -= f.get(0);
- vd.template getProp<force>(q)[1] -= f.get(1);
- vd.template getProp<force>(q)[2] -= f.get(2);
-
- //! \cond [add to q] \endcond
-
- // Next neighborhood
- ++Np;
- }
-
- // Next particle
- ++it2;
- }
-}
-
-//! \cond [calc forces vl] \endcond
-
-
-/*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
- *
- * ## Calculate energy ##
- *
- * No changes that been made to this code compared to the molecular dynamic example
- * with Verlet-list
- *
- * \see \ref e3_md_vl
- *
- *
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc energy vl
- *
- */
-
-//! \cond [calc energy vl] \endcond
-
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
-{
- double E = 0.0;
-
- double rc = r_cut*r_cut;
- double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
-
- // Get the iterator
- auto it2 = vd.getDomainIterator();
-
- // For each particle ...
- while (it2.isNext())
- {
- // ... p
- auto p = it2.get();
-
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
-
- // Get an iterator over the neighborhood of the particle p
- auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
-
- double Ep = E;
-
- // For each neighborhood of the particle p
- while (Np.isNext())
- {
- // Neighborhood particle q
- auto q = Np.get();
-
- // if p == q skip this particle
- if (q == p.getKey() && q < vd.size_local()) {++Np; continue;};
-
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
-
- // take the normalized direction
- double rn = norm2(xp - xq);
-
- if (rn >= r_cut*r_cut)
- {++Np;continue;}
-
- // potential energy (using pow is slower)
- E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
-
- // Next neighborhood
- ++Np;
- }
-
- // Kinetic energy of the particle given by its actual speed
- E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
- vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
- vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
-
- // Next Particle
- ++it2;
- }
-
- // Calculated energy
- return E;
-}
-
-//! \cond [calc energy vl] \endcond
-
-int main(int argc, char* argv[])
-{
- /*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
- *
- * ## Simulation ##
- *
- * The simulation is equal to the simulation explained in the example molecular dynamic
- *
- * \see \ref e3_md_vl
- *
- * The differences are that:
- *
- * * Create symmetric Verlet-list
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp sim verlet
- *
- * * For Energy calculation we have to create another verlet in this case a normal Verlet-list
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp norm verlet
- *
- * We need 2 verlet-list with one normal because we cannot use
- * symmetric Verlet to calculate total reductions like total energy
- *
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp simulation
- *
- */
-
- //! \cond [simulation] \endcond
-
- double dt = 0.00025;
- double sigma = 0.1;
- double r_cut = 3.0*sigma;
- double r_gskin = 1.3*r_cut;
- double sigma12 = pow(sigma,12);
- double sigma6 = pow(sigma,6);
-
- openfpm::vector<double> x;
- openfpm::vector<openfpm::vector<double>> y;
-
- openfpm_init(&argc,&argv);
- Vcluster & v_cl = create_vcluster();
-
- // we will use it do place particles on a 10x10x10 Grid like
- size_t sz[3] = {10,10,10};
-
- // domain
- Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
-
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
-
- // ghost, big enough to contain the interaction radius
- Ghost<3,float> ghost(r_gskin);
-
- vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
-
- auto it = vd.getGridIterator(sz);
-
- while (it.isNext())
- {
- vd.add();
-
- auto key = it.get();
-
- vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
- vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
- vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
-
- vd.template getLastProp<velocity>()[0] = 0.0;
- vd.template getLastProp<velocity>()[1] = 0.0;
- vd.template getLastProp<velocity>()[2] = 0.0;
-
- vd.template getLastProp<force>()[0] = 0.0;
- vd.template getLastProp<force>()[1] = 0.0;
- vd.template getLastProp<force>()[2] = 0.0;
-
- ++it;
- }
-
- timer tsim;
- tsim.start();
-
- //! \cond [sim verlet] \endcond
-
- // Get the Cell list structure
- auto NN = vd.getVerlet(r_gskin,VL_SYMMETRIC);
-
- //! \cond [sim verlet] \endcond
-
- //! \cond [norm verlet] \endcond
-
- auto NNE = vd.getVerlet(r_gskin);
-
- //! \cond [norm verlet] \endcond
-
- // calculate forces
- calc_forces(vd,NN,sigma12,sigma6,r_cut);
- unsigned long int f = 0;
-
- int cnt = 0;
- double max_disp = 0.0;
-
- // MD time stepping
- for (size_t i = 0; i < 10000 ; i++)
- {
- // Get the iterator
- auto it3 = vd.getDomainIterator();
-
- double max_displ = 0.0;
-
- // integrate velicity and space based on the calculated forces (Step1)
- while (it3.isNext())
- {
- auto p = it3.get();
-
- // here we calculate v(tn + 0.5)
- vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
- vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
- vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
-
- Point<3,double> disp({vd.template getProp<velocity>(p)[0]*dt,vd.template getProp<velocity>(p)[1]*dt,vd.template getProp<velocity>(p)[2]*dt});
-
- // here we calculate x(tn + 1)
- vd.getPos(p)[0] += disp.get(0);
- vd.getPos(p)[1] += disp.get(1);
- vd.getPos(p)[2] += disp.get(2);
-
- if (disp.norm() > max_displ)
- max_displ = disp.norm();
-
- ++it3;
- }
-
- if (max_disp < max_displ)
- max_disp = max_displ;
-
- // Because we moved the particles in space we have to map them and re-sync the ghost
- if (cnt % 10 == 0)
- {
- vd.map();
- vd.template ghost_get<>();
- // Get the Cell list structure
- vd.updateVerlet(NN,r_gskin,VL_SYMMETRIC);
- }
- else
- {
- vd.template ghost_get<>(SKIP_LABELLING);
- }
-
- cnt++;
-
- // calculate forces or a(tn + 1) Step 2
- calc_forces(vd,NN,sigma12,sigma6,r_cut);
-
- // Integrate the velocity Step 3
- auto it4 = vd.getDomainIterator();
-
- while (it4.isNext())
- {
- auto p = it4.get();
-
- // here we calculate v(tn + 1)
- vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
- vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
- vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
-
- ++it4;
- }
-
- // After every iteration collect some statistic about the confoguration
- if (i % 100 == 0)
- {
- // We write the particle position for visualization (Without ghost)
- vd.deleteGhost();
- vd.write("particles_",f);
-
- // we resync the ghost
- vd.ghost_get<>(SKIP_LABELLING);
-
- // update the verlet for energy calculation
- vd.updateVerlet(NNE,r_gskin);
-
- // We calculate the energy
- double energy = calc_energy(vd,NNE,sigma12,sigma6,r_cut);
- auto & vcl = create_vcluster();
- vcl.sum(energy);
- vcl.max(max_disp);
- vcl.execute();
-
- // we save the energy calculated at time step i c contain the time-step y contain the energy
- x.add(i);
- y.add({energy});
-
- // We also print on terminal the value of the energy
- // only one processor (master) write on terminal
- if (vcl.getProcessUnitID() == 0)
- std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
-
- max_disp = 0.0;
-
- f++;
- }
- }
-
- tsim.stop();
- std::cout << "Time: " << tsim.getwct() << std::endl;
-
- //! \cond [simulation] \endcond
-
- // Google charts options, it store the options to draw the X Y graph
- GCoptions options;
-
- // Title of the graph
- options.title = std::string("Energy with time");
-
- // Y axis name
- options.yAxis = std::string("Energy");
-
- // X axis name
- options.xAxis = std::string("iteration");
-
- // width of the line
- options.lineWidth = 1.0;
-
- // Object that draw the X Y graph
- GoogleChart cg;
-
- // Add the graph
- // The graph that it produce is in svg format that can be opened on browser
- cg.AddLinesGraph(x,y,options);
-
- // Write into html format
- cg.write("gc_plot2_out.html");
-
- //! \cond [google chart] \endcond
-
- /*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
- *
- * ## Finalize ## {#finalize_v_e3_md_sym}
- *
- * At the very end of the program we have always to de-initialize the library
- *
- * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp finalize
- *
- */
-
- //! \cond [finalize] \endcond
-
- openfpm_finalize();
-
- //! \cond [finalize] \endcond
-
- /*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic
- *
- * ## Full code ## {#code_v_e3_md_vl}
- *
- * \include Vector/3_molecular_dynamic/main_vl.cpp
- *
- */
-
- /*!
- * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
- *
- * ## Full code symmetric ## {#code_v_e3_md_sym}
- *
- * \include Vector/3_molecular_dynamic/main_vl_sym.cpp
- *
- */
-}
diff --git a/example/Vector/4_multiphase_celllist/main.cpp b/example/Vector/4_multiphase_celllist/main.cpp
deleted file mode 100644
index 60a2007c24246f1e0a369ad75bde9c0c980743fd..0000000000000000000000000000000000000000
--- a/example/Vector/4_multiphase_celllist/main.cpp
+++ /dev/null
@@ -1,240 +0,0 @@
-
-#include "Vector/vector_dist.hpp"
-#include "Decomposition/CartDecomposition.hpp"
-#include "data_type/aggregate.hpp"
-#include "NN/CellList/CellListM.hpp"
-
-/*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * [TOC]
- *
- *
- * # Vector Multi Phase cell-list # {#e4_ph_cl}
- *
- * This example show multi-phase cell lists for the distributed vector
- *
- * \warning BETA version
- *
- */
-
-int main(int argc, char* argv[])
-{
- /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * ## 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
- *
- * \snippet Vector/4_multiphase_celllist/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};
-
- // The domain
- Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
-
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
-
- float r_cut = 0.05;
-
- // 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>>(4096,box,bc,ghost) );
-
- // The other 3 phases
- phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
- phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
- phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
-
-
- //! \cond [Initialization and parameters] \endcond
-
-
- /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * ## Initialization ##
- *
- * We initialize all the phases with particle randomly positioned in the space
- *
- * \snippet Vector/4_multiphase_celllist/main.cpp rand dist
- *
- */
-
- //! \cond [rand dist] \endcond
-
- auto it = phases.get(0).getDomainIterator();
-
- // For all the particles
- while (it.isNext())
- {
- // for all phases
- for (size_t i = 0 ; i < phases.size() ; i++)
- {
- auto key = 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;
- }
- // next point
- ++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++)
- phases.get(i).map();
-
- //! \cond [rand dist] \endcond
-
- /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * ## Multi-phase cell-list construction ##
- *
- * 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.
- *
- * \snippet Vector/4_multiphase_celllist/main.cpp cl construction
- *
- */
-
- //! \cond [cl construction] \endcond
-
- // 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);
-
- CellListM<3,float,2> NN;
- 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).getDomainIterator();
- 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;
- }
- }
-
- //! \cond [cl construction] \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
- * 0 of the phase 0 the distance of the near particles from all the phases
- *
- * \snippet Vector/4_multiphase_celllist/main.cpp cl usage
- *
- */
-
- //! \cond [cl usage] \endcond
-
- vector_dist<3,float, aggregate<double,double> > & current_phase = phases.get(0);
-
- // Get the iterator of the particles of phase 0
- auto it2 = current_phase.getIterator();
-
- // For each particle ...
- while (it2.isNext())
- {
- // ... p
- auto p = it2.get();
-
- // Get the position of the particle p
- Point<3,float> xp = current_phase.getPos(p);
-
- // Get an iterator of all the particles neighborhood of p
- auto Np = NN.getNNIterator(NN.getCell(current_phase.getPos(p)));
-
- // For each particle near 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();
-
- Point<3,float> xq = phases.get(ph_q).getPos(q);
-
- // we accumulate all the distances
- current_phase.getProp<0>(p) = norm(xp - xq);
-
- ++Np;
- }
-
- // Next particle p
- ++it2;
- }
-
- //! \cond [cl usage] \endcond
-
- /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * ## Finalize ## {#finalize}
- *
- * At the very end of the program we have always de-initialize the library
- *
- * \snippet Vector/4_multiphase_celllist/main.cpp finalize
- *
- */
-
- //! \cond [finalize] \endcond
-
- openfpm_finalize();
-
- //! \cond [finalize] \endcond
-
- /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
- *
- * # Full code # {#code}
- *
- * \include Vector/4_multiphase_celllist/main.cpp
- *
- */
-}
-
-
-
-
diff --git a/example/Vector/4_multiphase_celllist/Makefile b/example/Vector/4_multiphase_celllist_verlet/Makefile
similarity index 69%
rename from example/Vector/4_multiphase_celllist/Makefile
rename to example/Vector/4_multiphase_celllist_verlet/Makefile
index 16168a458fb8b2d9c37bbe9e01d0b1ed7a3fbbaf..6043bb602cf6e5011ede7a7379252105c4d91db6 100644
--- a/example/Vector/4_multiphase_celllist/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/config.cfg b/example/Vector/4_multiphase_celllist_verlet/config.cfg
similarity index 100%
rename from example/Vector/4_multiphase_celllist/config.cfg
rename to example/Vector/4_multiphase_celllist_verlet/config.cfg
diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e426b889e26f9aec2d7d4a1a1e8c822a2e9c9809
--- /dev/null
+++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp
@@ -0,0 +1,575 @@
+
+#include "Vector/vector_dist.hpp"
+#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 and verlet
+ *
+ * [TOC]
+ *
+ *
+ * # Vector Multi Phase cell-list and Verlet # {#e4_ph_cl}
+ *
+ * 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.
+ *
+ *
+ */
+
+int main(int argc, char* argv[])
+{
+ /*!
+ * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+ *
+ * ## 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.
+ *
+ * \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 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});
+
+ // 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
+ phases.add( vector_dist<3,float, aggregate<double,double>>(4096,box,bc,ghost) );
+
+ // The other 3 phases
+ phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
+ phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
+ phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
+
+
+ //! \cond [Initialization and parameters] \endcond
+
+
+ /*!
+ * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+ *
+ * ## 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
+ *
+ * \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_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 of phase0
+ while (it.isNext())
+ {
+ // particle p
+ auto p = it.get();
+
+ // 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 particle
+ ++it;
+ }
+
+
+ // 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 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
+ 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)
+ *
+ * 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_verlet/main.cpp create multi-phase multi verlet
+ *
+ */
+
+ //! \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);
+
+ //! \cond [create multi-phase multi verlet] \endcond
+
+ /*!
+ * \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();
+
+ while (it.isNext())
+ {
+ // 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
+ CL_phase1 = phases.get(1).getCellListSym(r_cut);
+
+ // This function create a Verlet-list between phases 0 and 1
+ NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut);
+
+ // Get an iterator over the real and ghost particles
+ 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
+ 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
+ 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
+ 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())
+ {
+ // 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();
+
+ // 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;
+ }
+ }
+
+ // 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 ##
+ *
+ * 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_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
+ auto it2 = current_phase.getIterator();
+
+ // For each particle ...
+ while (it2.isNext())
+ {
+ // ... p
+ auto p = it2.get();
+
+ // 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 = CL_all.getNNIterator(CL_all.getCell(current_phase.getPos(p)));
+
+ // For each particle near 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();
+
+ Point<3,float> xq = phases.get(ph_q).getPos(q);
+
+ // we accumulate all the distances
+ current_phase.getProp<0>(p) = norm(xp - xq);
+
+ ++Np;
+ }
+
+ // Next particle p
+ ++it2;
+ }
+
+ //! \cond [cl usage] \endcond
+
+ /*!
+ * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+ *
+ * ## Finalize ## {#finalize}
+ *
+ * At the very end of the program we have always de-initialize the library
+ *
+ * \snippet Vector/4_multiphase_celllist_verlet/main.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+ *
+ * # Full code # {#code}
+ *
+ * \include Vector/4_multiphase_celllist_verlet/main.cpp
+ *
+ */
+}
+
+
+
+
diff --git a/example/Vector/6_complex_usage/Makefile b/example/Vector/6_complex_usage/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..98a3b9c708a9c4f1d4a3ecf84d83ca2caa43c686
--- /dev/null
+++ b/example/Vector/6_complex_usage/Makefile
@@ -0,0 +1,24 @@
+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+
+OBJ_DORD = main.o
+
+all: complex_use
+
+%.o: %.cpp
+ $(CC) -fext-numeric-literals -O3 -g -c --std=c++11 $(OPT) -o $@ $< $(INCLUDE_PATH)
+
+complex_use: $(OBJ_DORD)
+ $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+run: all_test
+ source $$HOME/openfpm_vars; mpirun -np 3 ./complex_use
+
+.PHONY: clean all run all_test on_test
+
+clean:
+ rm -f *.o *~ core complex_use
+
diff --git a/example/Vector/6_complex_usage/main.cpp b/example/Vector/6_complex_usage/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ecbd4b9a07f09ea12cab26888bdc1fdba1136255
--- /dev/null
+++ b/example/Vector/6_complex_usage/main.cpp
@@ -0,0 +1,409 @@
+#include "Vector/vector_dist.hpp"
+#include "data_type/aggregate.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Plot/util.hpp"
+#include "timer.hpp"
+
+/*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * [TOC]
+ *
+ * # Complex usage for validation and debugging # {#e6_cp_deb}
+ *
+ * In this example we show how the flexibility of the library can be used to perform complex
+ * tasks for validation and debugging. We will use a lot of feature that has been explained in
+ * the previous examples
+ *
+ * In a previous example we show in case of symmetric interaction between particles how to symmetric
+ * cell list could be used to speed up the calculation.
+ *
+ * \see not present
+ *
+ * In this example we validate that both computation match, in particular because the computation depend
+ * from the neighborhood, we will check and validate that the neighborhood of each
+ * particle is equivalent with symmetric cell list and normal cell-list.
+ *
+ */
+
+int main(int argc, char* argv[])
+{
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation
+ *
+ * ## Initialization ##
+ *
+ * The initialization is classical we define cutoff radius, domain, boundary conditions,
+ * ghost part and some useful constants. We also define a struct that will contain the debug information. This structure
+ * contain an id that is the global id of the particle and a point that is the the position of the
+ * neighborhood.
+ *
+ * \see \ref e0_s_init
+ *
+ * \snippet Vector/6_complex_usage/main.cpp initialization
+ *
+ */
+
+ //! \cond [initialization] \endcond
+
+ openfpm_init(&argc,&argv);
+
+ constexpr int gid = 0;
+ constexpr int nn_norm = 1;
+ constexpr int nn_sym = 2;
+
+ // used to define the domain
+ float L = 1000.0;
+
+ // Domain
+ Box<3,float> box({-L,-L,-L},{L,L,L});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // cut-off radius
+ float r_cut = 100.0;
+
+ // ghost
+ Ghost<3,float> ghost(r_cut);
+
+ // Point and global id
+ struct point_and_gid
+ {
+ // global id of the particle
+ size_t id;
+
+ // Position of the neighborhood particle
+ Point<3,float> xq;
+
+ // Used to reorder the neighborhood particles by id
+ bool operator<(const struct point_and_gid & pag)
+ {
+ return (id < pag.id);
+ }
+ };
+
+ //! \cond [initialization] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Particle ##
+ *
+ * For this example we will use a particle with 3 properties
+ *
+ * * The first property is a global index for the particle unique across processors
+ * * The second is a vector that contain for each neighborhood the global id of the particle
+ * and its position
+ * * The last one is again a vector equivalent to the previous but is produced with the symmetric
+ * cell-list
+ *
+ * \snippet Vector/6_complex_usage/main.cpp particle prop
+ *
+ */
+
+ //! \cond [particle prop] \endcond
+
+ // Particle properties list
+ typedef aggregate<size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
+
+ //! \cond [particle prop] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Distributed vector ##
+ *
+ * Here we create a distributed vector of 4096 particles. We initialize the particles
+ * randomly and we assign a global index to the first property of the particle
+ *
+ * The function accum return the following number
+ *
+ * \f$ \sum_{i < proc} size\_local(i) \f$
+ *
+ * where \f$ proc \f$ is the processor where we are computing the formula and \f$ size\_local(i) \f$
+ * is the number of particles the processor \f$ i \f$ has. This number is used to
+ * produce a global id of the particles
+ *
+ * \snippet Vector/6_complex_usage/main.cpp glob id part
+ *
+ */
+
+ //! \cond [glob id part] \endcond
+
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(4096,box,bc,ghost);
+
+ // used to calculate a global index for each particle
+ size_t start = vd.accum();
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ vd.getPos(key)[0] = 2.0*L*((float)rand()/RAND_MAX) - L;
+ vd.getPos(key)[1] = 2.0*L*((float)rand()/RAND_MAX) - L;
+ vd.getPos(key)[2] = 2.0*L*((float)rand()/RAND_MAX) - L;
+
+ vd.getProp<gid>(key) = key.getKey() + start;
+
+ ++it;
+ }
+
+ //! \cond [glob id part] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Redistribute ##
+ *
+ * Redistribute the particles and synchronize the ghosts. In this case we
+ * are only interested in synchronizing the global id property.
+ *
+ * \snippet Vector/6_complex_usage/main.cpp map and ghost
+ *
+ */
+
+ //! \cond [map and ghost] \endcond
+
+ vd.map();
+
+ // sync the ghost
+ vd.ghost_get<gid>();
+
+ //! \cond [map and ghost] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Calculate the neighborhood of each particle ##
+ *
+ * Here we calculate the neighborhood of each particles using a simple **cell list**.
+ * In case the particle q has a distance from the particle p smaller than r_cut.
+ * We add it in the first list of neighborhood
+ *
+ * \snippet Vector/6_complex_usage/main.cpp add nn particles
+ *
+ */
+
+ //! \cond [add nn particles] \endcond
+
+ auto NN = vd.getCellList(r_cut);
+ auto p_it = vd.getDomainIterator();
+
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ // skip self interaction
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // if the distance smalle than the cut-off radius add it to the neighborhood list
+ if (distance < r_cut )
+ {
+ vd.getProp<nn_norm>(p).add();
+ vd.getProp<nn_norm>(p).last().xq = xq;
+ vd.getProp<nn_norm>(p).last().id = vd.getProp<0>(q);
+ }
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // Next particle
+ ++p_it;
+ }
+
+ //! \cond [add nn particles] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Calculate the neighborhood of each particle with symmetric cell-list ##
+ *
+ * Here we calculate the neighborhood of each particle using instead a
+ * **symmetric cell list**. In case of symmetric cell-list if we find that a
+ * particle p is neighborhood of particle q we have to add p to the neighborhood
+ * of q and q to the neighborhood of p. Because q can be a ghost particle when
+ * we add the neighborhood p to q we have to transmit such information to the real
+ * owner of the particle. This can be done with the function **ghost_put**. In this
+ * case we use the operation **merge_** that add the already calculated neighborhood
+ * with the transmitted one. More in general it merge the information together.
+ *
+ * \snippet Vector/6_complex_usage/main.cpp add nn particles sym
+ *
+ */
+
+ //! \cond [add nn particles sym] \endcond
+
+ auto NN2 = vd.getCellListSym(r_cut);
+
+ auto p_it2 = vd.getDomainIterator();
+
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ // repulsive
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // Particle should be inside r_cut range
+
+ if (distance < r_cut )
+ {
+ vd.getProp<nn_sym>(p).add();
+ vd.getProp<nn_sym>(q).add();
+
+ vd.getProp<nn_sym>(p).last().xq = xq;
+ vd.getProp<nn_sym>(q).last().xq = xp;
+ vd.getProp<nn_sym>(p).last().id = vd.getProp<0>(q);
+ vd.getProp<nn_sym>(q).last().id = vd.getProp<0>(p);
+ }
+
+ ++Np;
+ }
+
+ ++p_it2;
+ }
+
+ vd.ghost_put<merge_,2>();
+
+ //! \cond [add nn particles sym] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Cheking for validation ##
+ *
+ * Here we check that the two calculated neighborhood match. In particular,
+ * because the order of the particles does not match, we have to first reorder
+ * by global-id, and than check the list. We cannot instead easily compare the
+ * position. The reason can be seen in this figure
+ *
+ \verbatim
+
+ +----------------+
+ | |
+ | 1 |
+ | * 2 |
+ | 3 * |
+ +<----- ghost * +<-------- real
+ |* <--- real |* <------- ghost
+ |4 |4
+ | |
+ | |
+ +----------------+
+
+
+
+ \endverbatim
+ *
+ * In the case we are calculating the neighborhood of the particle \f$ + \f$ in case of normal Cell-list.
+ * The real particle 1,2,3 are added, while the particle 4 is added with the ghost particle coordinates.
+ * In the case of symmetric cell-list the real 1,2,3 are added to \f$ + \f$. The particle 4
+ * instead is added by the ghost put. In particular 4 real add the particle to the \f$ + \f$ ghost particle
+ * and than ghost put merge the information. This mean that the symmetric add the particle 4 with the real coordinates
+ * of the particle 4
+ *
+ * \snippet Vector/6_complex_usage/main.cpp checking
+ *
+ */
+
+ //! \cond [checking] \endcond
+
+ auto p_it3 = vd.getDomainIterator();
+
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
+
+ vd.getProp<nn_norm>(p).sort();
+ vd.getProp<nn_sym>(p).sort();
+
+ ret &= vd.getProp<nn_norm>(p).size() == vd.getProp<nn_sym>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<1>(p).size() ; i++)
+ {
+ ret &= vd.getProp<nn_norm>(p).get(i).id == vd.getProp<nn_sym>(p).get(i).id;
+
+ if (box.isInside(vd.getProp<nn_norm>(p).get(i).xq) == true)
+ {
+ ret &= vd.getProp<nn_norm>(p).get(i).xq == vd.getProp<nn_sym>(p).get(i).xq;
+ }
+ }
+
+ ++p_it3;
+ }
+
+ if (ret != true)
+ {
+ std::cout << "ERROR" << std::endl;
+ exit(1);
+ }
+
+ //! \cond [checking] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * ## Finalize ##
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/6_complex_usage/main.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * # Full code # {#code}
+ *
+ * \include Vector/6_complex_usage/main.cpp
+ *
+ */
+}
diff --git a/openfpm_data b/openfpm_data
index fdf6629cc30bb52c73b703e2f66146439aebb7bc..ed732d7fe9e9644a9ad49f319923874b0c2c9e78 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit fdf6629cc30bb52c73b703e2f66146439aebb7bc
+Subproject commit ed732d7fe9e9644a9ad49f319923874b0c2c9e78
diff --git a/openfpm_devices b/openfpm_devices
index f46e31235688bd15c2a808ba7aa1bd413908bb8a..ecadd99fabe515ba06f00a5e7c656777c6bb9860 160000
--- a/openfpm_devices
+++ b/openfpm_devices
@@ -1 +1 @@
-Subproject commit f46e31235688bd15c2a808ba7aa1bd413908bb8a
+Subproject commit ecadd99fabe515ba06f00a5e7c656777c6bb9860
diff --git a/openfpm_io b/openfpm_io
index 62b8316404bf98c5cc5ad299563425656cfd3003..5cf769662d6e57a7a2fd12d1cfc835fd2793bc7b 160000
--- a/openfpm_io
+++ b/openfpm_io
@@ -1 +1 @@
-Subproject commit 62b8316404bf98c5cc5ad299563425656cfd3003
+Subproject commit 5cf769662d6e57a7a2fd12d1cfc835fd2793bc7b
diff --git a/openfpm_pdata.doc b/openfpm_pdata.doc
index 94101cdb7c8a63ebe6e1be42f6608f00bf1a879b..1375042ee7fb74340aa4414415ecafd645f56e6d 100644
--- a/openfpm_pdata.doc
+++ b/openfpm_pdata.doc
@@ -38,7 +38,7 @@ PROJECT_NAME = "OpenFPM_pdata"
# could be handy for archiving the generated documentation or if some version
# control system is used.
-PROJECT_NUMBER = 0.5.1
+PROJECT_NUMBER = 0.6.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
diff --git a/openfpm_vcluster b/openfpm_vcluster
index c65f473077ccf64b3e0848238fb235a38b443526..2d24d5c6dfab646fbeac429f53ef24a4e19631ac 160000
--- a/openfpm_vcluster
+++ b/openfpm_vcluster
@@ -1 +1 @@
-Subproject commit c65f473077ccf64b3e0848238fb235a38b443526
+Subproject commit 2d24d5c6dfab646fbeac429f53ef24a4e19631ac
diff --git a/script/install_PETSC.sh b/script/install_PETSC.sh
index e6537c72d9178fc24453c7c9d408363e9754a67b..ccb0c536fea18706c7cb2a2e1cefe2ba059feaba 100755
--- a/script/install_PETSC.sh
+++ b/script/install_PETSC.sh
@@ -181,6 +181,9 @@ if [ ! -d "$1/MUMPS" ]; then
$sed_command -i "/LIBBLAS\s=\s-lblas/c\LIBBLAS = -lopenblas" Makefile.inc
+ $sed_command -i "/INCPAR\s\+=\s\-I\/usr\/include/c\INCPAR =" Makefile.inc
+ $sed_command -i "/LIBPAR\s\+=\s\$(SCALAP)\s\-L\/usr\/lib\s\-lmpi/c\LIBPAR = \$(SCALAP)" Makefile.inc
+
make -j $2
if [ $? -eq 0 ]; then
diff --git a/src/Decomposition/ie_ghost.hpp b/src/Decomposition/ie_ghost.hpp
index 8b1bbe5fbe219033d624b7d9edc0293f26850bb5..e6d65ffc4479dbe391a4592c1260f01d380a701b 100755
--- a/src/Decomposition/ie_ghost.hpp
+++ b/src/Decomposition/ie_ghost.hpp
@@ -31,11 +31,10 @@ class ie_ghost
//! It store the same information of box_nn_processor_int organized by processor id
openfpm::vector< Box_dom<dim,T> > proc_int_box;
- // External ghost boxes for this processor, indicated with G8_0 G9_0 ...
+ //! External ghost boxes for this processor
openfpm::vector<p_box<dim,T> > vb_ext;
- // Internal ghost boxes for this processor domain, indicated with B8_0 B9_0 ..... in the figure
- // below as a linear vector
+ //! Internal ghost boxes for this processor domain
openfpm::vector<p_box<dim,T> > vb_int;
//! Cell-list that store the geometrical information of the internal ghost boxes
@@ -44,8 +43,10 @@ class ie_ghost
//! shift vectors
openfpm::vector<Point<dim,T>> shifts;
- // Temporal buffers to return information for ghost_processorID
+ //! Temporal buffers to return temporal information for ghost_processorID
openfpm::vector<std::pair<size_t,size_t>> ids_p;
+
+ //! Temporal buffers to return temporal information
openfpm::vector<size_t> ids;
@@ -188,6 +189,9 @@ protected:
* The geo cell list structure exist to speed up the labelling the points if they fall on some
* internal ghost
*
+ * \param domain where the cell list is defined
+ * \param div number of division of the cell list
+ *
*/
void Initialize_geo_cell(const Box<dim,T> & domain, const size_t (&div)[dim])
{
@@ -200,7 +204,11 @@ protected:
* For each sub-domain of the local processor it store the intersection between the enlarged
* sub-domain of the calling processor with the adjacent processors sub-domains (External ghost box)
*
+ * \param v_cl Virtual cluster
* \param ghost margins
+ * \param subdomains vector of local sundomains
+ * \param box_nn_processor it will store for each sub-domain the near processors
+ * \param nn_prcs contain the sub-domains of the near processors
*
* \note Are the G8_0 G9_0 G9_1 G5_0 boxes in calculateGhostBoxes
* \see calculateGhostBoxes
@@ -296,8 +304,8 @@ protected:
* \param v_cl Virtual cluster
* \param ghost margins
* \param sub_domains
- * \param box_nn_processors sub-domains of the adjacent processors
- * \param nn_prcs structure that store the adjacent processor information
+ * \param box_nn_processors sub-domains of the near processors
+ * \param nn_prcs structure that store the near processor sub-domains
* \param geo_cell Cell list that store the subdomain information
*
* \note Are the B8_0 B9_0 B9_1 B5_0 boxes in calculateGhostBoxes
@@ -508,6 +516,8 @@ public:
* This function return the set of shift vectors that determine such shift, for example
* in the example above the shift at position 5 will be (0,-1.0)
*
+ * \return the shift vectors
+ *
*/
const openfpm::vector<Point<dim,T>> & getShiftVectors()
{
@@ -667,6 +677,8 @@ public:
}
/*! \brief Given the internal ghost box id, it return the internal ghost box
+ *
+ * \param b_id internal ghost box id
*
* \return the internal ghost box
*
@@ -679,6 +691,8 @@ public:
/*! \brief Given the internal ghost box id, it return the near processor at witch belong
* or the near processor that produced this internal ghost box
*
+ * \param internal ghost box id
+ *
* \return the processor id of the ghost box
*
*/
@@ -698,6 +712,8 @@ public:
}
/*! \brief Given the external ghost box id, it return the external ghost box
+ *
+ * \param b_id external ghost box id
*
* \return the external ghost box
*
@@ -710,6 +726,8 @@ public:
/*! \brief Given the external ghost box id, it return the near processor at witch belong
* or the near processor that produced this external ghost box
*
+ * \param b_id external ghost box id
+ *
* \return the processor id of the external ghost box
*
*/
@@ -721,6 +739,7 @@ public:
/*! /brief Given a point it return the set of boxes in which the point fall
*
* \param p Point to check
+ *
* \return An iterator with the id's of the internal boxes in which the point fall
*
*/
@@ -757,6 +776,8 @@ public:
*
* \param return the processor ids (not the rank, the id in the near processor list)
*
+ * \return a vector of pairs containinf the requested infromation
+ *
*/
template <typename id1, typename id2> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(Point<dim,T> & p, const int opt = MULTIPLE)
{
@@ -806,6 +827,8 @@ public:
*
* \param return the processor ids
*
+ * \return a vector containing the requested information
+ *
*/
template <typename id> inline const openfpm::vector<size_t> ghost_processorID(Point<dim,T> & p, const int opt = MULTIPLE)
{
@@ -842,10 +865,12 @@ public:
/*! \brief Given a position it return if the position belong to any neighborhood processor ghost
* (Internal ghost)
*
- * \tparam id type of if to get box_id processor_id lc_processor_id
+ * \tparam id1 first index type to get box_id processor_id lc_processor_id
+ * \tparam id2 second index type to get box_id processor_id lc_processor_id
+ *
* \param p Particle position
*
- * \param return the processor ids
+ * \return a vector of pair containing the requested information
*
*/
template<typename id1, typename id2, typename Mem> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(const encapc<1,Point<dim,T>,Mem> & p, const int opt = MULTIPLE)
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index c16ae0f8fdf5843aabdeee22949c42995a31b2a1..208b654f9eb2edca1d9607b96ff49db1be20770d 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -200,7 +200,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
- ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
@@ -321,7 +320,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
- ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
diff --git a/src/Grid/grid_dist_id_unit_test.cpp b/src/Grid/grid_dist_id_unit_test.cpp
index ca7ad809a58fcfc07f360856d1e82988b31aba45..db7a81c8201d55f39a8c6199199278320a3c7312 100644
--- a/src/Grid/grid_dist_id_unit_test.cpp
+++ b/src/Grid/grid_dist_id_unit_test.cpp
@@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_domain_grid_unit_converter3D_test)
{
// Get the local hyper-cube
SpaceBox<3,float> sub = dec.getSubDomain(i);
- sub -= domain.getP1();
+// sub -= domain.getP1();
Box<3,size_t> g_box = g_dist.getCellDecomposer().convertDomainSpaceIntoGridUnits(sub,bc);
diff --git a/src/Grid/grid_dist_util.hpp b/src/Grid/grid_dist_util.hpp
index bba37eb0335b80bd7e99aeb03f9be68a185861d4..d16bf8bb1d073e85edd34b1a9d1778f61b823ea6 100644
--- a/src/Grid/grid_dist_util.hpp
+++ b/src/Grid/grid_dist_util.hpp
@@ -78,9 +78,7 @@ template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::ve
// Get the local sub-domain (Grid conversion must be done with the domain P1 equivalent to 0.0)
// consider that the sub-domain with point P1 equivalent to the domain P1 is a (0,0,0) in grid unit
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp = dec.getSubDomain(i);
- sp -= cd_sm.getOrig();
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp_g = dec.getSubDomainWithGhost(i);
- sp_g -= cd_sm.getOrig();
// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
SpaceBox<Decomposition::dims,long int> sp_t = cd_sm.convertDomainSpaceIntoGridUnits(sp,dec.periodicity());
diff --git a/src/Grid/staggered_dist_grid_util.hpp b/src/Grid/staggered_dist_grid_util.hpp
index c7d20ec6a1902d11264aa8ceeb7c67119ff368f8..938f1d9ab473645b22c4025043b49c17df9567fa 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 fff4127b0a3b696e5dc7eed7d653e86ddf9b38b1..3ba89adcd433d3aaf5eca6a0b02af994b68e499f 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.hpp b/src/Vector/vector_dist.hpp
index d2efb1f119d52055d5d0f880182f89bee7afb52d..ded398c37790ab18cb020533efe4ab45e7b36100 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -14,7 +14,6 @@
#include "Vector/vector_dist_iterator.hpp"
#include "Space/Shape/Box.hpp"
#include "Vector/vector_dist_key.hpp"
-#include "memory/PreAllocHeapMemory.hpp"
#include "memory/PtrMemory.hpp"
#include "NN/CellList/CellList.hpp"
#include "NN/CellList/CellListFast_hilb.hpp"
@@ -347,6 +346,36 @@ public:
return v_prp.template get<id>(g_m - 1);
}
+ /*! \brief Construct a cell list symmetric based on a cut of radius
+ *
+ * \tparam CellL CellList type to construct
+ *
+ * \param r_cut interation radius, or size of each cell
+ *
+ * \return the Cell list
+ *
+ */
+ template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellListSym(St r_cut)
+ {
+ // Cell list
+ CellL cell_list;
+
+ size_t pad = 0;
+ CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
+ cl_param_calculateSym(getDecomposition().getDomain(),cd_sm,getDecomposition().getGhost(),r_cut,pad);
+
+ // Processor bounding box
+ Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+ // Ghost padding extension
+ Ghost<dim,size_t> g_ext(0);
+ cell_list.Initialize(cd_sm,pbox,pad);
+
+ updateCellList(cell_list);
+
+ return cell_list;
+ }
+
/*! \brief Construct a cell list starting from the stored particles
*
* \tparam CellL CellList type to construct
@@ -392,21 +421,21 @@ public:
*/
template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > void updateCellList(CellL & cell_list)
{
- // Clear the cell list from the previous particles
- cell_list.clear();
-
- // for each particle add the particle to the cell list
-
- auto it = getIterator();
-
- while (it.isNext())
- {
- auto key = it.get();
+ populate_cell_list(v_pos,cell_list,g_m,CL_NON_SYMMETRIC);
- cell_list.add(this->getPos(key), key.getKey());
+ cell_list.set_gm(g_m);
+ }
- ++it;
- }
+ /*! \brief Update a cell list using the stored particles
+ *
+ * \tparam CellL CellList type to construct
+ *
+ * \param cell_list Cell list to update
+ *
+ */
+ template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > void updateCellListSym(CellL & cell_list)
+ {
+ populate_cell_list(v_pos,cell_list,g_m,CL_SYMMETRIC);
cell_list.set_gm(g_m);
}
@@ -481,15 +510,33 @@ public:
return cell_list;
}
+ /*! \brief for each particle get the symmetric verlet list
+ *
+ * \param r_cut cut-off radius
+ *
+ * \return the verlet list
+ *
+ */
+ VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
+ {
+ VerletList<dim,St,FAST,shift<dim,St>> ver;
+ // Processor bounding box
+ Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+ ver.InitializeSym(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);
+
+ return ver;
+ }
/*! \brief for each particle get the verlet list
*
* \param r_cut cut-off radius
- * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC
+ *
+ * \return a VerletList object
*
*/
- VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut, size_t opt = VL_NON_SYMMETRIC)
+ VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
{
VerletList<dim,St,FAST,shift<dim,St>> ver;
@@ -503,7 +550,7 @@ public:
// enlarge the box where the Verlet is defined
bt.enlarge(g);
- ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,opt);
+ ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);
return ver;
}
@@ -1040,6 +1087,9 @@ public:
pbox.enlarge(g);
cl_param_calculate(pbox, div,r_cut,enlarge);
+
+ // output the fixed domain
+ box = pbox;
}
/*! \brief It return the id of structure in the allocation list
@@ -1071,6 +1121,38 @@ public:
#endif
return v_cl;
}
+
+ /*! \brief return the position vector of all the particles
+ *
+ * \return the particle position vector
+ *
+ */
+ const openfpm::vector<Point<dim,St>> & getPosVector() const
+ {
+ return v_pos;
+ }
+
+ /*! \brief It return the sum of the particles in the previous processors
+ *
+ * \return the particles number
+ *
+ */
+ size_t accum()
+ {
+ openfpm::vector<size_t> accu;
+
+ size_t sz = size_local();
+
+ v_cl.allGather(sz,accu);
+ v_cl.execute();
+
+ sz = 0;
+
+ for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
+ sz += accu.get(i);
+
+ return sz;
+ }
};
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 0000000000000000000000000000000000000000..7d9eff81e52f5cbb5b5e2ad3b1951e69802ffbe0
--- /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_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index 0b747f83b14c439fb5afbadb3525a6e40c0ba1a4..bb878a10dff5af7066e4e83d726a5294064a4d91 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -332,339 +332,407 @@ BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
}
}
-
-BOOST_AUTO_TEST_CASE( vector_dist_sym_cell_list_test )
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
{
- long int k = 4096*create_vcluster().getProcessingUnits();
-
- long int big_step = k / 30;
- big_step = (big_step == 0)?1:big_step;
- long int small_step = 21;
+ Vcluster & v_cl = create_vcluster();
- if (create_vcluster().getProcessingUnits() > 48)
+ if (v_cl.getProcessingUnits() > 24)
return;
- print_test( "Testing vector symmetric cell-list k<=",k);
+ float L = 1000.0;
- // 3D test
- for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
- {
- double r_cut = 0.1;
+ // set the seed
+ // create the random generator engine
+ std::srand(0);
+ std::default_random_engine eg;
+ std::uniform_real_distribution<float> ud(-L,L);
- // domain
- Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ long int k = 4096 * v_cl.getProcessingUnits();
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ long int big_step = k / 4;
+ big_step = (big_step == 0)?1:big_step;
- // ghost, big enough to contain the interaction radius
- Ghost<3,float> ghost(r_cut);
+ print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+ BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
- vector_dist<3,double, aggregate<double> > vd(k,box,bc,ghost);
+ Box<3,float> box({-L,-L,-L},{L,L,L});
- {
- auto it = vd.getDomainIterator();
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
- while (it.isNext())
- {
- auto p = it.get();
+ float r_cut = 100.0;
- vd.getPos(p)[0] = (double)rand()/RAND_MAX;
- vd.getPos(p)[1] = (double)rand()/RAND_MAX;
- vd.getPos(p)[2] = (double)rand()/RAND_MAX;
+ // ghost
+ Ghost<3,float> ghost(r_cut);
- ++it;
- }
+ // Point and global id
+ struct point_and_gid
+ {
+ size_t id;
+ Point<3,float> xq;
+
+ bool operator<(const struct point_and_gid & pag)
+ {
+ return (id < pag.id);
}
+ };
- vd.map();
- vd.ghost_get<>();
+ typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
- // Get the Cell list structure
- auto NN = vd.getCellList(r_cut);
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+ size_t start = vd.init_size_accum(k);
- // Get an iterator over particles
- auto it2 = vd.getDomainAndGhostIterator();
+ auto it = vd.getIterator();
- openfpm::vector<openfpm::vector<size_t>> idx;
- idx.resize(vd.size_local_with_ghost());
+ while (it.isNext())
+ {
+ auto key = it.get();
- /////////// SYMMETRIC CASE CELL-LIST ////////
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
- // For each particle ...
- while (it2.isNext())
- {
- // ... p
- auto p = it2.get();
+ // Fill some properties randomly
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ vd.getProp<0>(key) = 0;
+ vd.getProp<1>(key) = 0;
+ vd.getProp<2>(key) = key.getKey() + start;
- // Get an iterator over the neighborhood of the particle p symmetric
- auto NpSym = NN.template getNNIteratorSym<NO_CHECK>(NN.getCell(vd.getPos(p)),p.getKey());
+ ++it;
+ }
- // For each neighborhood of the particle p
- while (NpSym.isNext())
- {
- // Neighborhood particle q
- auto q = NpSym.get();
+ vd.map();
- // if p == q skip this particle
- if (q == p.getKey() ) {++NpSym; continue;};
+ // sync the ghost
+ vd.ghost_get<0,2>();
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ auto NN = vd.getCellList(r_cut);
+ auto p_it = vd.getDomainIterator();
- // take the normalized direction
- double rn = norm2(xp - xq);
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
- // potential energy (using pow is slower)
- vd.getProp<0>(p) += rn;
- vd.getProp<0>(q) += rn;
+ Point<3,float> xp = vd.getPos(p);
- idx.get(p.getKey()).add(q);
- idx.get(q).add(p.getKey());
+ auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // Next neighborhood
- ++NpSym;
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
}
- // Next Particle
- ++it2;
- }
+ // repulsive
- /////////////// NON SYMMETRIC CASE ////////////////////////
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
- openfpm::vector<openfpm::vector<size_t>> idx2;
- idx2.resize(vd.size_local());
+ float distance = f.norm();
- auto it = vd.getDomainIterator();
+ // Particle should be inside 2 * r_cut range
- // For each particle ...
- while (it.isNext())
- {
- // ... p
- auto p = it.get();
+ if (distance < r_cut )
+ {
+ vd.getProp<0>(p)++;
+ vd.getProp<3>(p).add();
+ vd.getProp<3>(p).last().xq = xq;
+ vd.getProp<3>(p).last().id = vd.getProp<2>(q);
+ }
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ ++Np;
+ }
- // Get an iterator over the neighborhood of the particle p
- auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+ ++p_it;
+ }
- double Ep = 0.0;
+ // We now try symmetric Cell-list
- // For each neighborhood of the particle p
- while (Np.isNext())
- {
- // Neighborhood particle q
- auto q = Np.get();
+ auto NN2 = vd.getCellListSym(r_cut);
- // if p == q skip this particle
- if (q == p.getKey()) {++Np; continue;};
+ auto p_it2 = vd.getDomainIterator();
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
- // take the normalized direction
- double rn = norm2(xp - xq);
+ Point<3,float> xp = vd.getPos(p);
- idx2.get(p.getKey()).add(q);
+ auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
- // potential energy (using pow is slower)
- Ep += rn;
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // Next neighborhood
+ if (p.getKey() == q)
+ {
++Np;
+ continue;
}
- idx.get(p.getKey()).sort();
- idx2.get(p.getKey()).sort();
+ // repulsive
- bool ret = true;
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
- for (size_t i = 0 ; i < idx.get(p.getKey()).size() ; i++)
- ret &= idx.get(p.getKey()).get(i) == idx2.get(p.getKey()).get(i);
+ float distance = f.norm();
- BOOST_REQUIRE_EQUAL(ret,true);
+ // Particle should be inside r_cut range
- BOOST_REQUIRE_CLOSE(Ep,vd.getProp<0>(p),0.01);
+ if (distance < r_cut )
+ {
+ vd.getProp<1>(p)++;
+ vd.getProp<1>(q)++;
- // Next Particle
- ++it;
+ vd.getProp<4>(p).add();
+ vd.getProp<4>(q).add();
+
+ vd.getProp<4>(p).last().xq = xq;
+ vd.getProp<4>(q).last().xq = xp;
+ vd.getProp<4>(p).last().id = vd.getProp<2>(q);
+ vd.getProp<4>(q).last().id = vd.getProp<2>(p);
+ }
+
+ ++Np;
}
+
+ ++p_it2;
+ }
+
+ vd.ghost_put<add_,1>();
+ vd.ghost_put<merge_,4>();
+
+ auto p_it3 = vd.getDomainIterator();
+
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
+
+ ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+
+ vd.getProp<3>(p).sort();
+ vd.getProp<4>(p).sort();
+
+ ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+ ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
+
+ if (ret == false)
+ break;
+
+ ++p_it3;
}
+
+ BOOST_REQUIRE_EQUAL(ret,true);
}
-BOOST_AUTO_TEST_CASE( vector_dist_sym_verlet_list_test )
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
{
- if (create_vcluster().getProcessingUnits() > 48)
+ Vcluster & v_cl = create_vcluster();
+
+ if (v_cl.getProcessingUnits() > 24)
return;
- long int k = 4096*create_vcluster().getProcessingUnits();
+ float L = 1000.0;
+
+ // set the seed
+ // create the random generator engine
+ std::srand(0);
+ std::default_random_engine eg;
+ std::uniform_real_distribution<float> ud(-L,L);
+
+ long int k = 4096 * v_cl.getProcessingUnits();
- long int big_step = k / 30;
+ long int big_step = k / 4;
big_step = (big_step == 0)?1:big_step;
- long int small_step = 21;
- print_test( "Testing vector symmetric verlet-list k<=",k);
+ print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+ BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
- // 3D test
- for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
- {
- double r_cut = 0.1;
+ Box<3,float> box({-L,-L,-L},{L,L,L});
- // domain
- Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ float r_cut = 100.0;
- // ghost, big enough to contain the interaction radius
- Ghost<3,float> ghost(r_cut);
+ // ghost
+ Ghost<3,float> ghost(r_cut);
- vector_dist<3,double, aggregate<double> > vd(k,box,bc,ghost);
+ // Point and global id
+ struct point_and_gid
+ {
+ size_t id;
+ Point<3,float> xq;
+ bool operator<(const struct point_and_gid & pag)
{
- auto it = vd.getDomainIterator();
+ return (id < pag.id);
+ }
+ };
- while (it.isNext())
- {
- auto p = it.get();
+ typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
- vd.getPos(p)[0] = (double)rand()/RAND_MAX;
- vd.getPos(p)[1] = (double)rand()/RAND_MAX;
- vd.getPos(p)[2] = (double)rand()/RAND_MAX;
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+ size_t start = vd.init_size_accum(k);
- ++it;
- }
- }
+ auto it = vd.getIterator();
- vd.map();
- vd.ghost_get<>();
+ while (it.isNext())
+ {
+ auto key = it.get();
- // Get the Cell list structure
- auto NN = vd.getVerlet(r_cut,VL_SYMMETRIC);
- auto NNc = vd.getVerlet(r_cut);
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
- // Get an iterator over particles
- auto it2 = vd.getDomainAndGhostIterator();
+ // Fill some properties randomly
- openfpm::vector<openfpm::vector<size_t>> idx;
- idx.resize(vd.size_local_with_ghost());
+ vd.getProp<0>(key) = 0;
+ vd.getProp<1>(key) = 0;
+ vd.getProp<2>(key) = key.getKey() + start;
- /////////// SYMMETRIC CASE VERLET-LIST ////////
+ ++it;
+ }
- // For each particle ...
- while (it2.isNext())
- {
- // ... p
- auto p = it2.get();
+ vd.map();
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ // sync the ghost
+ vd.ghost_get<0,2>();
- // Get an iterator over the neighborhood of the particle p symmetric
- auto NpSym = NN.template getNNIterator<NO_CHECK>(p.getKey());
+ auto NN = vd.getVerlet(r_cut);
+ auto p_it = vd.getDomainIterator();
- // For each neighborhood of the particle p
- while (NpSym.isNext())
- {
- // Neighborhood particle q
- auto q = NpSym.get();
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
- // if p == q skip this particle
- if (q == p.getKey()) {++NpSym; continue;};
+ Point<3,float> xp = vd.getPos(p);
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ auto Np = NN.getNNIterator(p.getKey());
- // take the normalized direction
- double rn = norm2(xp - xq);
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // potential energy (using pow is slower)
- vd.getProp<0>(p) += rn;
- vd.getProp<0>(q) += rn;
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ // repulsive
- idx.get(p.getKey()).add(q);
- idx.get(q).add(p.getKey());
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+ float distance = f.norm();
- // Next neighborhood
- ++NpSym;
+ // Particle should be inside 2 * r_cut range
+
+ if (distance < r_cut )
+ {
+ vd.getProp<0>(p)++;
+ vd.getProp<3>(p).add();
+ vd.getProp<3>(p).last().xq = xq;
+ vd.getProp<3>(p).last().id = vd.getProp<2>(q);
}
- // Next Particle
- ++it2;
+ ++Np;
}
- /////////////// NON SYMMETRIC CASE VERLET-LIST ////////////////////////
+ ++p_it;
+ }
- openfpm::vector<openfpm::vector<size_t>> idx2;
- idx2.resize(vd.size_local());
+ // We now try symmetric Cell-list
- auto it = vd.getDomainIterator();
+ auto NN2 = vd.getVerletSym(r_cut);
- // For each particle ...
- while (it.isNext())
- {
- // ... p
- auto p = it.get();
+ auto p_it2 = vd.getDomainIterator();
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
- // Get an iterator over the neighborhood of the particle p
- auto Np = NNc.template getNNIterator<NO_CHECK>(p.getKey());
+ auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
- double Ep = 0.0;
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // For each neighborhood of the particle p
- while (Np.isNext())
+ if (p.getKey() == q)
{
- // Neighborhood particle q
- auto q = Np.get();
+ ++Np;
+ continue;
+ }
- // if p == q skip this particle
- if (q == p.getKey()) {++Np; continue;};
+ // repulsive
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
- // take the normalized direction
- double rn = norm2(xp - xq);
+ float distance = f.norm();
- idx2.get(p.getKey()).add(q);
+ // Particle should be inside r_cut range
- // potential energy (using pow is slower)
- Ep += rn;
+ if (distance < r_cut )
+ {
+ vd.getProp<1>(p)++;
+ vd.getProp<1>(q)++;
- // Next neighborhood
- ++Np;
+ vd.getProp<4>(p).add();
+ vd.getProp<4>(q).add();
+
+ vd.getProp<4>(p).last().xq = xq;
+ vd.getProp<4>(q).last().xq = xp;
+ vd.getProp<4>(p).last().id = vd.getProp<2>(q);
+ vd.getProp<4>(q).last().id = vd.getProp<2>(p);
}
- idx.get(p.getKey()).sort();
- idx2.get(p.getKey()).sort();
+ ++Np;
+ }
- bool ret = true;
+ ++p_it2;
+ }
- BOOST_REQUIRE_EQUAL(idx.get(p.getKey()).size(),idx2.get(p.getKey()).size());
+ vd.ghost_put<add_,1>();
+ vd.ghost_put<merge_,4>();
- for (size_t i = 0 ; i < idx.get(p.getKey()).size() ; i++)
- {
- ret &= idx.get(p.getKey()).get(i) == idx2.get(p.getKey()).get(i);
- }
+ auto p_it3 = vd.getDomainIterator();
- BOOST_REQUIRE_EQUAL(ret,true);
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
- BOOST_REQUIRE_CLOSE(Ep,vd.getProp<0>(p),0.01);
+ ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
- // Next Particle
- ++it;
- }
+ vd.getProp<3>(p).sort();
+ vd.getProp<4>(p).sort();
+
+ ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+ ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
+
+ if (ret == false)
+ break;
+
+ ++p_it3;
}
+
+ BOOST_REQUIRE_EQUAL(ret,true);
}
diff --git a/src/Vector/vector_dist_comm.hpp b/src/Vector/vector_dist_comm.hpp
index 00baad340d0631944fba8356b72f859aa317274e..6b43019bcacf22891e135488606dd3a732197f1b 100644
--- a/src/Vector/vector_dist_comm.hpp
+++ b/src/Vector/vector_dist_comm.hpp
@@ -11,6 +11,7 @@
#define V_SUB_UNIT_FACTOR 64
#define SKIP_LABELLING 512
+#define KEEP_PROPERTIES 512
#define NO_POSITION 1
#define WITH_POSITION 2
@@ -166,33 +167,6 @@ class vector_dist_comm
return end_id;
}
- /*! \brief It store for each processor the position and properties vector of the particles
- *
- * This structure is used in the map function
- *
- */
- struct pos_prop
- {
- //! position vector
- openfpm::vector<Point<dim, St>, PreAllocHeapMemory<2>, typename memory_traits_lin<Point<dim, St>>::type, memory_traits_lin, openfpm::grow_policy_identity> pos;
- //! properties vector
- openfpm::vector<prop, PreAllocHeapMemory<2>, typename memory_traits_lin<prop>::type, memory_traits_lin, openfpm::grow_policy_identity> prp;
- };
-
- /*! \brief for each processor store 2 vector containing the sending buffers
- *
- * This structure is used in the map_list function
- *
- */
- template <typename sel_prop>
- struct pos_prop_sel
- {
- //! position vector
- openfpm::vector<Point<dim, St>, PreAllocHeapMemory<2>, typename memory_traits_lin<Point<dim, St>>::type, memory_traits_lin, openfpm::grow_policy_identity> pos;
- //! properties vector
- openfpm::vector<sel_prop, PreAllocHeapMemory<2>, typename memory_traits_lin<sel_prop>::type, memory_traits_lin, openfpm::grow_policy_identity> prp;
- };
-
//! Flags that indicate that the function createShiftBox() has been called
bool is_shift_box_created = false;
@@ -275,8 +249,7 @@ class vector_dist_comm
// add this particle shifting its position
v_pos.add(p);
- v_prp.add();
- v_prp.last() = v_prp.get(key);
+ v_prp.get(lg_m+i) = v_prp.get(key);
}
}
@@ -393,7 +366,8 @@ class vector_dist_comm
// Create the shift boxes
createShiftBox();
- lg_m = v_prp.size();
+ if (!(opt & SKIP_LABELLING))
+ lg_m = v_prp.size();
if (box_f.size() == 0)
return;
@@ -749,7 +723,7 @@ public:
*
*/
vector_dist_comm(const vector_dist_comm<dim,St,prop,Decomposition,Memory> & v)
- :v_cl(create_vcluster()),dec(create_vcluster())
+ :v_cl(create_vcluster()),dec(create_vcluster()),lg_m(0)
{
this->operator=(v);
}
@@ -761,7 +735,7 @@ public:
*
*/
vector_dist_comm(const Decomposition & dec)
- :v_cl(create_vcluster()),dec(dec)
+ :v_cl(create_vcluster()),dec(dec),lg_m(0)
{
}
@@ -781,7 +755,7 @@ public:
*
*/
vector_dist_comm()
- :v_cl(create_vcluster()),dec(create_vcluster())
+ :v_cl(create_vcluster()),dec(create_vcluster()),lg_m(0)
{
}
@@ -845,10 +819,13 @@ public:
// send vector for each processor
typedef openfpm::vector<prp_object> send_vector;
- // reset the ghost part
- if (opt != NO_POSITION)
+ if (!(opt & NO_POSITION))
v_pos.resize(g_m);
- v_prp.resize(g_m);
+
+ // reset the ghost part
+
+ if (!(opt & SKIP_LABELLING))
+ v_prp.resize(g_m);
// Label all the particles
if ((opt & SKIP_LABELLING) == false)
@@ -859,14 +836,21 @@ public:
fill_send_ghost_prp_buf<send_vector, prp_object, prp...>(v_prp,g_send_prp);
// Create and fill the send buffer for the particle position
- if (opt != NO_POSITION)
+ if (!(opt & NO_POSITION))
fill_send_ghost_pos_buf(v_pos,g_pos_send);
prc_recv_get.clear();
recv_sz_get.clear();
- v_cl.SSendRecvP<send_vector,decltype(v_prp),prp...>(g_send_prp,v_prp,prc_g_opart,prc_recv_get,recv_sz_get);
- if (opt != NO_POSITION)
+ if (opt & SKIP_LABELLING)
+ {
+ op_ssend_gg_recv_merge opm(g_m);
+ v_cl.SSendRecvP_op<op_ssend_gg_recv_merge,send_vector,decltype(v_prp),prp...>(g_send_prp,v_prp,prc_g_opart,opm,prc_recv_get,recv_sz_get);
+ }
+ else
+ v_cl.SSendRecvP<send_vector,decltype(v_prp),prp...>(g_send_prp,v_prp,prc_g_opart,prc_recv_get,recv_sz_get);
+
+ if (!(opt & NO_POSITION))
{
prc_recv_get.clear();
recv_sz_get.clear();
@@ -1074,6 +1058,28 @@ public:
// Send and receive ghost particle information
op_ssend_recv_merge<op> opm(g_opart);
v_cl.SSendRecvP_op<op_ssend_recv_merge<op>,send_vector,decltype(v_prp),prp...>(g_send_prp,v_prp,prc_recv_get,opm,prc_recv_put,recv_sz_put);
+
+ // process also the local replicated particles
+
+ size_t i2 = 0;
+
+
+ if (v_prp.size() - lg_m != o_part_loc.size())
+ {
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Local ghost particles = " << v_prp.size() - lg_m << " != " << o_part_loc.size() << std::endl;
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Check that you did a ghost_get before a ghost_put" << std::endl;
+ }
+
+
+ for (size_t i = lg_m ; i < v_prp.size() ; i++)
+ {
+ auto dst = v_prp.get(o_part_loc.template get<0>(i2));
+ auto src = v_prp.get(i);
+ copy_cpu_encap_encap_op_prp<op,decltype(v_prp.get(0)),decltype(v_prp.get(0)),prp...> cp(src,dst);
+
+ boost::mpl::for_each_ref< boost::mpl::range_c<int,0,sizeof...(prp)> >(cp);
+ i2++;
+ }
}
};
diff --git a/src/Vector/vector_dist_multiphase_functions.hpp b/src/Vector/vector_dist_multiphase_functions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a3c6e32bd51a19a6d4fddb68b1f75d398401963
--- /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 88a7275a268faa0eebc47208390061bb33de8a76..07b360353afb58a8d2faced5010b340c070c84b6 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1239,8 +1239,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
vd.ghost_get<0>();
- vd.write("Debug_output");
-
// calculate the distance of the first, second and third neighborhood particle
// Consider that they are on a regular grid
@@ -1437,12 +1435,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
// ghost
Ghost<3,float> ghost(0.1);
- typedef aggregate<float> part_prop;
+ typedef aggregate<float,float,float> part_prop;
// Distributed vector
vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
@@ -1460,6 +1458,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
// Fill some properties randomly
vd.getProp<0>(key) = 0.0;
+ vd.getProp<1>(key) = vd.getPos(key)[0];
+ vd.getProp<2>(key) = vd.getPos(key)[0]*vd.getPos(key)[0];
++it;
}
@@ -1467,22 +1467,21 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
vd.map();
// sync the ghost
- vd.ghost_get<0>();
-
- openfpm::vector<size_t> list_idx;
- openfpm::vector<size_t> list_idx2;
+ vd.ghost_get<0,1,2>();
- auto it3 = vd.getGhostIterator();
- while (it3.isNext())
+ bool ret = true;
+ auto it2 = vd.getGhostIterator();
+ while (it2.isNext())
{
- auto key = it3.get();
+ auto key = it2.get();
- list_idx.add(key.getKey());
+ ret &= vd.getProp<1>(key) == vd.getPos(key)[0];
+ ret &= vd.getProp<2>(key) == vd.getPos(key)[0] * vd.getPos(key)[0];
- ++it3;
+ ++it2;
}
- list_idx.sort();
+ BOOST_REQUIRE_EQUAL(ret,true);
for (size_t i = 0 ; i < 10 ; i++)
{
@@ -1492,7 +1491,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
{
auto key = it.get();
- vd.getPos(key)[0] = ud(eg);
vd.getPos(key)[1] = ud(eg);
vd.getPos(key)[2] = ud(eg);
@@ -1505,7 +1503,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
vd.ghost_get<0>(SKIP_LABELLING);
- list_idx2.clear();
auto it2 = vd.getGhostIterator();
bool ret = true;
@@ -1513,40 +1510,36 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
{
auto key = it2.get();
- list_idx2.add(key.getKey());
ret &= vd.getProp<0>(key) == i;
+ ret &= vd.getProp<1>(key) == vd.getPos(key)[0];
+ ret &= vd.getProp<2>(key) == vd.getPos(key)[0] * vd.getPos(key)[0];
++it2;
}
- BOOST_REQUIRE_EQUAL(ret,true);
- BOOST_REQUIRE_EQUAL(list_idx.size(),list_idx2.size());
-
- list_idx2.sort();
-
- ret = true;
- for (size_t i = 0 ; i < list_idx.size() ; i++)
- ret &= list_idx.get(i) == list_idx2.get(i);
-
BOOST_REQUIRE_EQUAL(ret,true);
}
vd.map();
- vd.ghost_get<0>();
+ vd.ghost_get<0,1,2>();
- list_idx.clear();
+ // shift the particle position by 1.0
- auto it4 = vd.getGhostIterator();
- while (it4.isNext())
+ it = vd.getGhostIterator();
+ while (it.isNext())
{
- auto key = it4.get();
+ // Particle p
+ auto p = it.get();
- list_idx.add(key.getKey());
+ // we shift down he particles
+ vd.getPos(p)[0] = 10.0;
- ++it4;
- }
+ // we shift
+ vd.getPos(p)[1] = 17.0;
- list_idx.sort();
+ // next particle
+ ++it;
+ }
for (size_t i = 0 ; i < 10 ; i++)
{
@@ -1556,42 +1549,37 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
{
auto key = it.get();
- vd.getPos(key)[0] = ud(eg);
vd.getPos(key)[1] = ud(eg);
vd.getPos(key)[2] = ud(eg);
// Fill some properties randomly
vd.getProp<0>(key) = i;
+ vd.getProp<1>(key) = vd.getPos(key)[0];
+ vd.getProp<2>(key) = vd.getPos(key)[0]*vd.getPos(key)[0];
++it;
}
- vd.ghost_get<0>(SKIP_LABELLING);
+ vd.ghost_get<0>(SKIP_LABELLING | NO_POSITION);
- list_idx2.clear();
auto it2 = vd.getGhostIterator();
bool ret = true;
while (it2.isNext())
{
- auto key = it2.get();
+ // Particle p
+ auto p = it.get();
- list_idx2.add(key.getKey());
- ret &= vd.getProp<0>(key) == i;
+ ret &= vd.getPos(p)[0] == 10.0;
+ // we shift
+ ret &= vd.getPos(p)[1] == 17.0;
+
+ // next particle
++it2;
}
- BOOST_REQUIRE_EQUAL(ret,true);
- BOOST_REQUIRE_EQUAL(list_idx.size(),list_idx2.size());
-
- list_idx2.sort();
-
- ret = true;
- for (size_t i = 0 ; i < list_idx.size() ; i++)
- ret &= list_idx.get(i) == list_idx2.get(i);
-
BOOST_REQUIRE_EQUAL(ret,true);
}
}
@@ -1618,13 +1606,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
// 3D test
for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
{
+ float r_cut = 1.3 / k;
+ float r_g = 1.5 / k;
+
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
- Ghost<3,float> ghost(1.3/(k));
+ Ghost<3,float> ghost(r_g);
typedef aggregate<float> part_prop;
@@ -1656,7 +1647,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
vd.ghost_get<0>();
{
- auto NN = vd.getCellList(1.3/k);
+ auto NN = vd.getCellList(r_cut);
float a = 1.0f*k*k;
// run trough all the particles + ghost
@@ -1680,8 +1671,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
float dist = xp.distance(xq);
- if (dist < 1.0/k)
- vd.getProp<0>(q) += a*(-dist*dist+1.0/k/k);
+ if (dist < r_cut)
+ vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
++Np;
}
@@ -1712,8 +1703,18 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
BOOST_REQUIRE_EQUAL(ret,true);
}
+ auto itp = vd.getDomainAndGhostIterator();
+ while (itp.isNext())
+ {
+ auto key = itp.get();
+
+ vd.getProp<0>(key) = 0.0;
+
+ ++itp;
+ }
+
{
- auto NN = vd.getCellList(1.3/k);
+ auto NN = vd.getCellList(r_cut);
float a = 1.0f*k*k;
// run trough all the particles + ghost
@@ -1737,8 +1738,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
float dist = xp.distance(xq);
- if (dist < 1.0/k)
- vd.getProp<0>(q) += a*(-dist*dist+1.0/k/k);
+ if (dist < r_cut)
+ vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
++Np;
}
@@ -1817,6 +1818,7 @@ BOOST_AUTO_TEST_CASE( vector_of_vector_dist )
BOOST_REQUIRE_EQUAL(cnt,4*4096ul);
}
+
#include "vector_dist_cell_list_tests.hpp"
#include "vector_dist_NN_tests.hpp"
#include "vector_dist_complex_prp_unit_test.hpp"
diff --git a/src/debug.hpp b/src/debug.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d550941bfda0f57d0ba31de964b2c3671604e4bd
--- /dev/null
+++ b/src/debug.hpp
@@ -0,0 +1,36 @@
+/*
+ * debug.hpp
+ *
+ * Created on: Oct 7, 2016
+ * Author: i-bird
+ */
+
+#ifndef SRC_DEBUG_HPP_
+#define SRC_DEBUG_HPP_
+
+#include "Vector/map_vector.hpp"
+#include "Space/Shape/Point.hpp"
+
+Point<3,float> getPosPoint3f(openfpm::vector<Point<3,float>> & pos, size_t i)
+{
+ return Point<3,float>(pos.get(i));
+}
+
+Point<3,double> getPosPoint3d(openfpm::vector<Point<3,double>> & pos, size_t i)
+{
+ return Point<3,double>(pos.get(i));
+}
+
+Point<2,float> getPosPoint2f(openfpm::vector<Point<2,float>> & pos, size_t i)
+{
+ return Point<2,float>(pos.get(i));
+}
+
+Point<2,double> getPosPoint2d(openfpm::vector<Point<2,double>> & pos, size_t i)
+{
+ return Point<2,double>(pos.get(i));
+}
+
+
+
+#endif /* SRC_DEBUG_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 38a88d937269940f14d9282b2a8c655dda788382..c60b455f760aa5c64e8c240c5e273acc34046d6e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -22,6 +22,7 @@ int main(int argc, char* argv[])
return boost::unit_test::unit_test_main( &init_unit_test, argc, argv );
}
+#include "debug.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Point_test.hpp"
#include "Decomposition/CartDecomposition.hpp"
@@ -48,4 +49,5 @@ int main(int argc, char* argv[])
#include "Graph/DistGraphFactory.hpp"
#include "Decomposition/nn_processor_unit_test.hpp"
#include "Grid/staggered_grid_dist_unit_test.hpp"
+#include "Vector/vector_dist_MP_unit_tests.hpp"
//#include "antoniol_test_isolation.hpp"