diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp
index 0e6a334477539f0632710a74a26b19a9236587e4..a2c9b7b250983ac08a397c71c3e9d62c88ff1323 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -9,6 +9,7 @@
* \subpage Vector_4_cp
* \subpage Vector_4_mp_cl
* \subpage Vector_5_md_vl_sym
+ * \subpage Vector_5_md_vl_sym_crs
* \subpage Vector_6_complex_usage
*
*/
diff --git a/example/Vector/5_molecular_dynamic_sym_crs/Makefile b/example/Vector/5_molecular_dynamic_sym_crs/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..cf0391fab74d27ace055d5904d8df2c180b2266d
--- /dev/null
+++ b/example/Vector/5_molecular_dynamic_sym_crs/Makefile
@@ -0,0 +1,24 @@
+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+
+OBJ_DORD = main.o
+
+all: md_sym
+
+%.o: %.cpp
+ $(CC) -O3 -g -c --std=c++11 $(OPT) -o $@ $< $(INCLUDE_PATH)
+
+md_sym: $(OBJ_DORD)
+ $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+run: md_sym
+ mpirun -np 3 ./md_sym
+
+.PHONY: clean all run
+
+clean:
+ rm -f *.o *~ core md_sym
+
diff --git a/example/Vector/5_molecular_dynamic_sym_crs/config.cfg b/example/Vector/5_molecular_dynamic_sym_crs/config.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..1eecbac3577c765edca7f90cf5f61cfb6b9f4880
--- /dev/null
+++ b/example/Vector/5_molecular_dynamic_sym_crs/config.cfg
@@ -0,0 +1,2 @@
+[pack]
+files = main.cpp Makefile
diff --git a/example/Vector/5_molecular_dynamic_sym_crs/main.cpp b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fe6bec6ba6215df5bab4fa418542dc06edd351a5
--- /dev/null
+++ b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp
@@ -0,0 +1,571 @@
+#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_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ *
+ * # Molecular dynamic with symmetric interactions (Crossing scheme)# {#md_e5_sym_crs}
+ *
+ * In a previous example we show how to build a parallel molecular dynamic simulation in the case
+ * of symmetric interaction.
+ *
+ * \see \ref Vector_5_md_vl_sym
+ *
+ * Symmetric interaction has one big advantage to cut by half the computation.
+ * Unfortunately has also one disadvantage. In the standard way it increase communication overhead.
+ * This effect can be seen in a simple way consider the non-symmetric case.
+ *
+ * \verbatim
+
+ Non-symmetric Symmetric
+
++-----------------+ +-----------------+ +-----------------+
+|XXXXXXXXXXXXXXXXX| |XXXXXXXXXXXXXXXXX| |XXXXXXXXXXXXXXXXX|
+|XX+-----------+XX| |XX+-----------+XX| |XX+-----------+XX|
+|XX| |XX| |XX| |XX| |XX| |XX|
+|XX| |XX| |XX| |XX| |XX| |XX|
+|XX| |XX| |XX| |XX| + |XX| |XX|
+|XX| |XX| |XX| |XX| |XX| |XX|
+|XX| |XX| |XX| |XX| |XX| |XX|
+|XX| |XX| |XX| |XX| |XX| |XX|
+|XX+-----------+XX| +-----------------+ +-----------------+
+|XXXXXXXXXXXXXXXXX|
++-----------------+ ghost_get ghost_put
+
+ ghost_get
+
+
+ * \endverbatim
+ *
+ * In the non-symmetric case we do a ghost-get communicating the area indicated
+ * with "X". After this we are able to calculate the forces inside the box.
+ * In the case of normal symmetric like explained in the previous example the first
+ * ghost_get communicate the area indicated by X and the ghost_put communicate
+ * an equivalent area. This mean that The communication area for the non-symmetric
+ * is smaller than the symmetric one. This mean that overall the standard symmetric
+ * requires more communication. In order to reduce the communication of the standard
+ * symmetric we can use the crossing scheme. The crossing scheme do the following
+ *
+ * \verbatim
+
+
+ 1 2 3 1 2 3
+ X 4 X 3 + X 4
+
+ \endverbatim
+ *
+ * If X is the Cell and the numbers are the neighborhood cells, the concept of the crossing scheme is to
+ * take the interaction X with "1" and shift by one to the right. This mean that for each cell "X"
+ * in our domain we have to do the interaction (X-1),(X-2),(X-3),(3,4). Note that in the case of standard
+ * symmetric, all the interactions are in the top,left and right cells of the domain cells. Because of this
+ * we can eliminate the bottom part in the ghost_get and in the subsequent ghost_put. In the case of
+ * the crossing scheme we have only interaction with the top and right cells, so we do not need the
+ * bottom and left part of the ghost.
+ *
+ * \verbatim
+
+ Symmetric crs
+
+ +--------------+ +--------------+
+ |XXXXXXXXXXXXXX| |XXXXXXXXXXXXXX|
+ +-----------+XX| +-----------+XX|
+ | |XX| | |XX|
+ | |XX| | |XX|
+ | |XX| + | |XX|
+ | |XX| | |XX|
+ | |XX| | |XX|
+ | |XX| | |XX|
+ +--------------+ +--------------+
+
+ ghost_get ghost_put
+
+
+ \endverbatim
+ *
+ * In this example we modify the molecular dynamic example to use
+ * this scheme.
+ *
+ */
+
+//! \cond [constants] \endcond
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+//! \cond [constants] \endcond
+
+/*!
+ *
+ * \page Vector_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ * ## Calculate forces ## {#md_e5_calc_force_crs}
+ *
+ * In this function we calculate the forces between particles. It require the vector of particles
+ * the Verlet-list and sigma factor for the Lennard-Jhones potential. The function is exactly the same
+ * as the normal symmetric
+ *
+ * \see \ref md_e5_sym_expl
+ *
+ * with the following changes
+ *
+ *
+ * * We use a different particle iterator compared to the getDomainIterator() this because the
+ * the particles to iterate are different from the simple domain particles.
+ *
+ *
+ * The calculation is the same as the normal symmetric case
+ *
+ *
+ * \snippet Vector/5_molecular_dynamic_sym_crs/main.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)
+{
+ // Reset force on the ghost
+
+ auto itg = vd.getDomainAndGhostIterator();
+
+ while (itg.isNext())
+ {
+ auto p = itg.get();
+
+ // Reset force
+ vd.getProp<force>(p)[0] = 0.0;
+ vd.getProp<force>(p)[1] = 0.0;
+ vd.getProp<force>(p)[2] = 0.0;
+
+ ++itg;
+ }
+
+ //! \cond [real and ghost] \endcond
+
+ // Get an iterator over particles
+ auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList());
+
+ //! \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.getNNIterator<NO_CHECK>(p);
+
+ // For each neighborhood particle ...
+ while (Np.isNext())
+ {
+ // ... q
+ auto q = Np.get();
+
+ // if (p == q) skip this particle
+ if (q == p) {++Np; continue;};
+
+ // Get the position of q
+ 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 [ghost_put] \endcond
+
+ // Sum the contribution to the real particles
+ vd.ghost_put<add_,force>();
+
+ //! \cond [ghost_put] \endcond
+}
+
+//! \cond [calc forces vl] \endcond
+
+
+/*!
+ *
+ * \page Vector_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ * ## Calculate energy ## {#md_e5_calc_one_crs}
+ *
+ * For the energy we use symmetric verlet-list in the same way as we did in the previous example.
+ * The changes are the following.
+ *
+ * * We use a different particle iterator compared to the getDomainIterator() this because the
+ * the particles to iterate are different from the simple domain particles.
+ *
+ * * Because we are iterating domain particles and ghost particles. For the Kinetic energy we have
+ * to filter-out ghost particles
+ *
+ * \snippet Vector/5_molecular_dynamic_sym/main.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 = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
+
+ // Get an iterator over particles
+ auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList());
+
+ // 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 particles of p
+ auto Np = NN.getNNIterator<NO_CHECK>(p);
+
+ 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) {++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 += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // To note that the crossing scheme go across the domain particles +
+ // some ghost particles. This mean that we have to filter out the ghost
+ // particles otherwise we count double energies
+ //
+ if (p < vd.size_local())
+ {
+ // 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_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ * ## Simulation ## {#md_e5_sym_sim_crs}
+ *
+ * The simulation is equal to the simulation explained in the example molecular dynamic
+ *
+ * \see \ref md_e5_sym
+ *
+ * The difference is that we create a symmetric Verlet-list for crossing scheme instead of a normal one
+ * \snippet Vector/5_molecular_dynamic_sym_crs/main.cpp sim verlet
+ *
+ * The rest of the code remain unchanged
+ *
+ * \snippet Vector/5_molecular_dynamic_sym_crs/main.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);
+ ghost.setLow(0,0.0);
+ ghost.setLow(1,0.0);
+ ghost.setLow(2,0.0);
+
+ vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost,BIND_DEC_TO_GHOST);
+
+ size_t k = 0;
+ size_t start = vd.accum();
+
+ 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;
+
+ k++;
+ ++it;
+ }
+
+ vd.map();
+ vd.ghost_get<>();
+
+ timer tsim;
+ tsim.start();
+
+ //! \cond [sim verlet] \endcond
+
+ // Get the Cell list structure
+ auto NN = vd.getVerletCrs(r_gskin);;
+
+ //! \cond [sim 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_CRS_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<>();
+
+
+ // We calculate the energy
+ double energy = calc_energy(vd,NN,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_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ * ## Finalize ## {#finalize_v_e5_md_sym_crs}
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/5_molecular_dynamic_sym_crs/main.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_5_md_vl_sym_crs Vector 5 molecular dynamic with symmetric Verlet list crossing scheme
+ *
+ * ## Full code ## {#full_code_v_e5_md_sym_crs}
+ *
+ * \include Vector/5_molecular_dynamic_sym_crs/main.cpp
+ *
+ */
+}
diff --git a/openfpm_data b/openfpm_data
index 155cc3537485efb4993a667330a01a7d2644449a..7e326c0134dedecc6cf06dfb9e3a466e66280300 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 155cc3537485efb4993a667330a01a7d2644449a
+Subproject commit 7e326c0134dedecc6cf06dfb9e3a466e66280300
diff --git a/openfpm_pdata.doc b/openfpm_pdata.doc
index 1375042ee7fb74340aa4414415ecafd645f56e6d..21a5e06b66dc77ca1b8be3e6157fc538b4f1de26 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.6.0
+PROJECT_NUMBER = 0.7.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/src/Decomposition/nn_processor.hpp b/src/Decomposition/nn_processor.hpp
index 54f3249ea23527ff52ad40b8636ea3ecd79c8e63..2bc8e30605dadebdad675366b6e2941184a6c245 100755
--- a/src/Decomposition/nn_processor.hpp
+++ b/src/Decomposition/nn_processor.hpp
@@ -623,7 +623,7 @@ public:
aBC=true;
- return add_box_periodic(domain,ghost,bc);
+ add_box_periodic(domain,ghost,bc);
}
/*! \brief Check if the nn_prcs contain the same information
diff --git a/src/Makefile.am b/src/Makefile.am
index 04449f1ce4088cf7e3dc42c46ecd8fe26fa5ffff..b2cad2a544a2d5c45102690be9e6589d9fe5cf2c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ pdata_CXXFLAGS = $(OPENMP_CFLAGS) $(AM_CXXFLAGS) $(LIBHILBERT_INCLUDE) $(PETSC_I
pdata_CFLAGS = $(CUDA_CFLAGS)
pdata_LDADD = $(LINKLIBS) -lparmetis -lmetis
nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartDecomposition_ext.hpp Decomposition/common.hpp Decomposition/Decomposition.hpp Decomposition/ie_ghost.hpp \
- Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \
+ Decomposition/Domain_NN_calculator_cart.hpp 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_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 \
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 97a433e09286ace21218022183dbf311786d765d..6080050f36d2c478ae34a5e7f2d07ada0e89f0cb 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -548,6 +548,43 @@ public:
return ver;
}
+ /*! \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> > getVerletCrs(St r_cut)
+ {
+ VerletList<dim,St,FAST,shift<dim,St>> ver;
+
+ // Processor bounding box
+ Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+ // Initialize the verlet list
+ ver.InitializeCrs(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);
+
+ // Get the internal cell list
+ auto & NN = ver.getInternalCellList();
+
+ // Shift
+ grid_key_dx<dim> cell_shift = NN.getShift();
+
+ // Shift
+ grid_key_dx<dim> shift = NN.getShift();
+
+ // Add padding
+ for (size_t i = 0 ; i < dim ; i++)
+ shift.set_d(i,shift.get(i) + NN.getPadding(i));
+
+ grid_sm<dim,void> gs = NN.getInternalGrid();
+
+ ver.createVerletCrs(r_cut,g_m,v_pos,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));
+
+ return ver;
+ }
+
/*! \brief for each particle get the verlet list
*
* \param r_cut cut-off radius
@@ -584,7 +621,27 @@ public:
*/
void updateVerlet(VerletList<dim,St,FAST,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
{
- ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
+ if (opt == VL_CRS_SYMMETRIC)
+ {
+ // Get the internal cell list
+ auto & NN = ver.getInternalCellList();
+
+ // Shift
+ grid_key_dx<dim> cell_shift = NN.getShift();
+
+ // Shift
+ grid_key_dx<dim> shift = NN.getShift();
+
+ // Add padding
+ for (size_t i = 0 ; i < dim ; i++)
+ shift.set_d(i,shift.get(i) + NN.getPadding(i));
+
+ grid_sm<dim,void> gs = NN.getInternalGrid();
+
+ ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));
+ }
+ else
+ ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
}
/*! \brief for each particle get the verlet list
diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index 4328e97caf882eec0d5a76354f59b1daff2e782c..8308909df42c16493101da543a2eb10c9a1016a5 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -564,6 +564,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
// ghost
Ghost<3,float> ghost(r_cut);
+ Ghost<3,float> ghost2(r_cut);
+ ghost2.setLow(0,0.0);
+ ghost2.setLow(1,0.0);
+ ghost2.setLow(2,0.0);
// Point and global id
struct point_and_gid
@@ -581,6 +585,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
// Distributed vector
vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
+ vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
size_t start = vd.init_size_accum(k);
auto it = vd.getIterator();
@@ -593,19 +598,32 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
vd.getPos(key)[1] = ud(eg);
vd.getPos(key)[2] = ud(eg);
+ vd2.getPos(key)[0] = vd.getPos(key)[0];
+ vd2.getPos(key)[1] = vd.getPos(key)[1];
+ vd2.getPos(key)[2] = vd.getPos(key)[2];
+
// Fill some properties randomly
vd.getProp<0>(key) = 0;
vd.getProp<1>(key) = 0;
vd.getProp<2>(key) = key.getKey() + start;
+ vd2.getProp<0>(key) = 0;
+ vd2.getProp<1>(key) = 0;
+ vd2.getProp<2>(key) = key.getKey() + start;
+
++it;
}
vd.map();
+ vd2.map();
// sync the ghost
vd.ghost_get<0,2>();
+ vd2.ghost_get<0,2>();
+
+ vd2.write("CRS_output");
+ vd2.getDecomposition().write("CRS_output_dec");
auto NN = vd.getCellList(r_cut);
auto p_it = vd.getDomainIterator();
@@ -653,20 +671,20 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
// We now try symmetric Cell-list
- auto NN2 = vd.getCellListSym(r_cut);
+ auto NN2 = vd2.getCellListSym(r_cut);
// In case of CRS we have to iterate particles within some cells
// here we define whichone
- auto p_it2 = vd.getParticleIteratorCRS(NN2);
+ auto p_it2 = vd2.getParticleIteratorCRS(NN2);
// For each particle
while (p_it2.isNext())
{
auto p = p_it2.get();
- Point<3,float> xp = vd.getPos(p);
+ Point<3,float> xp = vd2.getPos(p);
- auto Np = p_it2.getNNIteratorCSR(vd.getPosVector());
+ auto Np = p_it2.getNNIteratorCSR(vd2.getPosVector());
while (Np.isNext())
{
@@ -680,7 +698,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
// repulsive
- Point<3,float> xq = vd.getPos(q);
+ Point<3,float> xq = vd2.getPos(q);
Point<3,float> f = (xp - xq);
float distance = f.norm();
@@ -689,16 +707,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
if (distance < r_cut )
{
- vd.getProp<1>(p)++;
- vd.getProp<1>(q)++;
+ vd2.getProp<1>(p)++;
+ vd2.getProp<1>(q)++;
- vd.getProp<4>(p).add();
- vd.getProp<4>(q).add();
+ vd2.getProp<4>(p).add();
+ vd2.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);
+ vd2.getProp<4>(p).last().xq = xq;
+ vd2.getProp<4>(q).last().xq = xp;
+ vd2.getProp<4>(p).last().id = vd2.getProp<2>(q);
+ vd2.getProp<4>(q).last().id = vd2.getProp<2>(p);
}
++Np;
@@ -707,8 +725,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
++p_it2;
}
- vd.ghost_put<add_,1>();
- vd.ghost_put<merge_,4>();
+ vd2.ghost_put<add_,1>();
+ vd2.ghost_put<merge_,4>();
auto p_it3 = vd.getDomainIterator();
@@ -717,25 +735,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
{
auto p = p_it3.get();
- ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+ ret &= vd2.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();
-
-// if (v_cl.getProcessUnitID() == 0)
-// {
-// std::cerr << "Particle: " << p.getKey() << " Position: " << Point<3,float>(vd.getPos(p)).toString() << std::endl;
+ vd.getProp<3>(p).sort();
+ vd2.getProp<4>(p).sort();
- for (size_t i = 0 ; i < vd.getProp<4>(p).size() ; i++)
- {
-// std::cerr << "POSITION: " << vd.getProp<3>(p).get(i).xq.toString() << std::endl;
+ ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size();
-// std::cerr << "ID nn " << vd.getProp<3>(p).get(i).id << " " << vd.getProp<4>(p).get(i).id << std::endl;
- ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
- }
-// }
+ for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+ ret &= vd.getProp<3>(p).get(i).id == vd2.getProp<4>(p).get(i).id;
if (ret == false)
break;
@@ -1169,4 +1178,220 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
}
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
+{
+ Vcluster & v_cl = create_vcluster();
+
+ if (v_cl.getProcessingUnits() > 24)
+ return;
+
+ 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 / 4;
+ big_step = (big_step == 0)?1:big_step;
+
+ print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+ BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
+
+ Box<3,float> box({-L,-L,-L},{L,L,L});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ float r_cut = 100.0;
+
+ // ghost
+ Ghost<3,float> ghost(r_cut);
+ Ghost<3,float> ghost2(r_cut);
+ ghost2.setLow(0,0.0);
+ ghost2.setLow(1,0.0);
+ ghost2.setLow(2,0.0);
+
+ // Point and global id
+ struct point_and_gid
+ {
+ size_t id;
+ Point<3,float> xq;
+
+ bool operator<(const struct point_and_gid & pag) const
+ {
+ return (id < pag.id);
+ }
+ };
+
+ typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
+
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
+ vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
+ size_t start = vd.init_size_accum(k);
+
+ auto it = vd.getIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
+
+ vd2.getPos(key)[0] = vd.getPos(key)[0];
+ vd2.getPos(key)[1] = vd.getPos(key)[1];
+ vd2.getPos(key)[2] = vd.getPos(key)[2];
+
+ // Fill some properties randomly
+
+ vd.getProp<0>(key) = 0;
+ vd.getProp<1>(key) = 0;
+ vd.getProp<2>(key) = key.getKey() + start;
+
+ vd2.getProp<0>(key) = 0;
+ vd2.getProp<1>(key) = 0;
+ vd2.getProp<2>(key) = key.getKey() + start;
+
+ ++it;
+ }
+
+ vd.map();
+ vd2.map();
+
+ // sync the ghost
+ vd.ghost_get<0,2>();
+ vd2.ghost_get<0,2>();
+
+ auto NN = vd.getVerlet(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(p.getKey());
+
+ 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 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);
+ }
+
+ ++Np;
+ }
+
+ ++p_it;
+ }
+
+ // We now try symmetric Verlet-list Crs scheme
+
+ auto NN2 = vd2.getVerletCrs(r_cut);
+
+ // Because iterating across particles in the CSR scheme require a Cell-list
+ auto p_it2 = vd2.getParticleIteratorCRS(NN2.getInternalCellList());
+
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
+
+ Point<3,float> xp = vd2.getPos(p);
+
+ auto Np = NN2.getNNIterator<NO_CHECK>(p);
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ if (p == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ // repulsive
+
+ Point<3,float> xq = vd2.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ if (distance < r_cut )
+ {
+ vd2.getProp<1>(p)++;
+ vd2.getProp<1>(q)++;
+
+ vd2.getProp<4>(p).add();
+ vd2.getProp<4>(q).add();
+
+ vd2.getProp<4>(p).last().xq = xq;
+ vd2.getProp<4>(q).last().xq = xp;
+ vd2.getProp<4>(p).last().id = vd2.getProp<2>(q);
+ vd2.getProp<4>(q).last().id = vd2.getProp<2>(p);
+ }
+
+ ++Np;
+ }
+
+ ++p_it2;
+ }
+
+ vd2.ghost_put<add_,1>();
+ vd2.ghost_put<merge_,4>();
+
+ auto p_it3 = vd.getDomainIterator();
+
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
+
+ ret &= vd2.getProp<1>(p) == vd.getProp<0>(p);
+
+
+ vd.getProp<3>(p).sort();
+ vd2.getProp<4>(p).sort();
+
+ ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+ ret &= vd.getProp<3>(p).get(i).id == vd2.getProp<4>(p).get(i).id;
+
+ if (ret == false)
+ break;
+
+ ++p_it3;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+}
+
#endif /* SRC_VECTOR_VECTOR_DIST_CELL_LIST_TESTS_HPP_ */