We'll be taking GitLab down for maintenance around 22 in the evening on the 15th of September, so this Sunday. Let us know (tt.mpi-cbg.de) if you experience any issues with it after the maintenance period.

Commit c3556b53 authored by Pietro Incardona's avatar Pietro Incardona

Ghost put working

parent d9335a87
......@@ -29,7 +29,6 @@ int main(int argc, char* argv[])
//
// Get the vcluster object and the number of processor
//
Vcluster & v_cl = create_vcluster();
long int N_prc = v_cl.getProcessingUnits();
......@@ -37,7 +36,7 @@ int main(int argc, char* argv[])
// ### WIKI 3 ###
//
// We find the maximum of the processors rank, that should be the Number of
// processora minus one, only processor 0 print on terminal
// processor minus one, only processor 0 print on terminal
//
long int id = v_cl.getProcessUnitID();
......
......@@ -344,7 +344,7 @@ int main(int argc, char* argv[])
/*!
* \page Vector_0_simple Vector 0 simple
*
* ## Finalize ## {#finalize}
* ## Finalize ## {#finalize_e0_sim}
*
* At the very end of the program we have always to de-initialize the library
*
......@@ -361,7 +361,7 @@ int main(int argc, char* argv[])
/*!
* \page Vector_0_simple Vector 0 simple
*
* # Full code # {#code}
* ## Full code ## {#code_e0_sim}
*
* \include Vector/0_simple/main.cpp
*
......
......@@ -12,7 +12,7 @@ OBJ_VL_SYM = main_vl_sym.o
all: md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
%.o: %.cpp
$(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
$(CC) -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
md_dyn: $(OBJ)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
......@@ -27,7 +27,7 @@ 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; mpirun -np 3 ./md_dyn_vl_sym;
.PHONY: clean all run
......
/*! \page Vector_3_md Vector 3 molecular dynamic
*
* \subpage Vector_3_md_dyn
* \subpage Vector_3_md_vl
*
*/
#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
......@@ -6,7 +13,7 @@
#include "timer.hpp"
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* [TOC]
*
......@@ -31,7 +38,7 @@ constexpr int force = 1;
/*!
*
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Calculate forces ## {#e3_md_cf}
*
......@@ -54,7 +61,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
//! \cond [calc forces] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* This function in called several time and require neighborhood of each particle. In order to speed-up the
* Cell-list construction we can use updateCellList function to reuse the memory of the previous cell-list.
......@@ -74,7 +81,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
/*!
*
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md Vector 3 molecular dynamic with cell-list
*
* Get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
* and calculate the force based on the Lennard-Jhones potential given by
......@@ -155,7 +162,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
//! \cond [calc forces2] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Calculate energy ## {#e3_md_ce}
*
......@@ -176,7 +183,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
//! \cond [calc energy] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Reset the counter for the energy counter and
* update the cell list from the actual particle configuration
......@@ -196,7 +203,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
/*!
*
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* First we get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
* and calculate the energy based on the Lennard-Jhones potential given by
......@@ -273,7 +280,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
int main(int argc, char* argv[])
{
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Initialization ## {#e3_md_init}
*
......@@ -282,11 +289,11 @@ int main(int argc, char* argv[])
* size of the box, and cut-off radius of the interaction. We also define 2 vectors
* x and y (they are like std::vector) used for statistic
*
* \snippet Vector/3_molecular_dynamic/main.cpp constants
* \snippet Vector/3_molecular_dynamic/main.cpp constants run
*
*/
//! \cond [constants] \endcond
//! \cond [constants run] \endcond
double dt = 0.0005;
double sigma = 0.1;
......@@ -297,10 +304,10 @@ int main(int argc, char* argv[])
openfpm::vector<double> x;
openfpm::vector<openfpm::vector<double>> y;
//! \cond [constants] \endcond
//! \cond [constants run] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Here we Initialize the library, we create a Box that define our domain, boundary conditions and ghost
*
......@@ -330,7 +337,7 @@ int main(int argc, char* argv[])
//! \cond [init] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Than we define a distributed vector in 3D, containing 2 vectorial properties the
* first is the actual velocity of the particle the other is the force
......@@ -348,7 +355,7 @@ int main(int argc, char* argv[])
//! \cond [vect create] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Particles on a grid like position ## {#e3_md_gl}
*
......@@ -389,7 +396,7 @@ int main(int argc, char* argv[])
//! \cond [vect grid] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Molecular dynamic steps ## {#e3_md_vi}
*
......@@ -397,16 +404,16 @@ int main(int argc, char* argv[])
*
* The verlet integration stepping look like this
*
* $$ \vec{v}(t_{n+1/2}) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) $$ // Step 1
* $$ \vec{x}(t_{n}) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) $$ // Step 1
* \f[ \vec{v}(t_{n+1/2}) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f]
* \f[ \vec{x}(t_{n}) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f]
*
* calculate the forces $$ \vec{a} (t_{n}) $$ from $$ \vec{x} (t_{n}) $$ // Step 2
* calculate the forces from \f$ \vec{a} (t_{n}) \f$ finally
*
* $$ \vec{v}(t_{n+1}) = \vec{v}_p(t_n+1/2) + \frac{1}{2} \delta t \vec{a}(t_n+1) $$ // Step 3
* \f[ \vec{v}(t_{n+1}) = \vec{v}_p(t_n+1/2) + \frac{1}{2} \delta t \vec{a}(t_n+1) \f]
*
* The cell-list structure is required to calculate forces
*
* Inside this cycle we are using several features that has been explained before in particuilar
* Inside this cycle we are using several features that has been explained before in particular
*
* \see \ref e0_s_assign_pos
*
......@@ -516,7 +523,7 @@ int main(int argc, char* argv[])
std::cout << "Time: " << tsim.getwct() << std::endl;
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Plotting graphs ## {#e3_md_pg}
*
......@@ -560,13 +567,13 @@ int main(int argc, char* argv[])
//! \cond [google chart] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Finalize ## {#finalize}
* ## Finalize ## {#finalize_v_e3_md}
*
* At the very end of the program we have always to de-initialize the library
*
* \snippet Vector/1_celllist/main.cpp finalize
* \snippet Vector/3_molecular_dynamic/main.cpp finalize
*
*/
......@@ -577,18 +584,18 @@ int main(int argc, char* argv[])
//! \cond [finalize] \endcond
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* # Full code # {#code}
* ## Full code ## {#code_v_e3_md}
*
* \include Vector/3_molecular_dynamic/main.cpp
*
*/
/*!
* \page Vector_3_md Vector 3 molecular dynamic
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* # Code with expression # {#code}
* ## Code with expression ## {#code_v_e3_md_expr}
*
* Here we also show how we can simplify the example using expressions
*
......
openfpm_data @ 36003173
Subproject commit 9301fe459ca2ef81115c642ee7a54c983ddfd6f1
Subproject commit 360031739584423b20497b1364d2c3bcf26c5143
......@@ -844,6 +844,18 @@ public:
this->template ghost_get_<prp...>(v_pos,v_prp,g_m,opt);
}
/*! \brief It synchronize the properties and position of the ghost particles
*
* \tparam op which kind of operation to apply
* \tparam prp list of properties to get synchronize
*
*
*/
template<template<typename,typename> class op, int ... prp> inline void ghost_put()
{
this->template ghost_put_<op,prp...>(v_pos,v_prp,g_m);
}
/*! \brief Remove a set of elements from the distributed vector
*
* \warning keys must be sorted
......
......@@ -407,7 +407,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_sym_cell_list_test )
auto q = NpSym.get();
// if p == q skip this particle
if (q == p.getKey()) {++NpSym; continue;};
if (q == p.getKey() ) {++NpSym; continue;};
// Get position of the particle q
Point<3,double> xq = vd.getPos(q);
......@@ -664,262 +664,5 @@ BOOST_AUTO_TEST_CASE( vector_dist_sym_verlet_list_test )
}
}
struct couple_contrib
{
size_t i;
size_t j;
size_t ir;
size_t jr;
double E;
bool operator<(couple_contrib & ct)
{
return E < ct.E;
}
};
BOOST_AUTO_TEST_CASE( vector_dist_sym_tot_red_computation )
{
if (create_vcluster().getProcessingUnits() > 12)
return;
long int k = 4096*create_vcluster().getProcessingUnits();
size_t base = create_vcluster().getProcessUnitID() * 4096;
long int big_step = k / 30;
big_step = (big_step == 0)?1:big_step;
long int small_step = 21;
print_test( "Testing vector symmetric verlet-list k<=",k);
openfpm::vector<couple_contrib> cct;
openfpm::vector<couple_contrib> cct2;
// 3D test
/* for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
{*/
double r_cut = 0.05;
// 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};
// ghost, big enough to contain the interaction radius
Ghost<3,float> ghost(r_cut);
vector_dist<3,double, aggregate<size_t> > vd(k,box,bc,ghost);
{
auto it = vd.getDomainIterator();
while (it.isNext())
{
auto p = it.get();
vd.getPos(p)[0] = (double)rand()/RAND_MAX;
vd.getPos(p)[1] = (double)rand()/RAND_MAX;
vd.getPos(p)[2] = (double)rand()/RAND_MAX;
vd.getProp<0>(p) = p.getKey() + base;
++it;
}
}
vd.map();
vd.ghost_get<0>();
// Get the Cell list structure
auto NN = vd.getVerlet(r_cut,VL_SYMMETRIC_RED);
auto NNc = vd.getVerlet(r_cut);
// Get an iterator over particles
auto it2 = vd.getDomainAndGhostIterator();
/////////// SYMMETRIC CASE VERLET-LIST REDUCTION ////////
{
auto it2 = vd.getDomainAndGhostIterator();
// DEBUG
// For each particle ...
while (it2.isNext())
{
// ... p
auto p = it2.get();
if( vd.getProp<0>(p) == 2524 || vd.getProp<0>(p) == 4135 )
{
Vcluster & v_cl = create_vcluster();
Point<3,double> xp = vd.getPos(p);
std::cout << "Particle " << p.getKey() << " " << xp.toString() << " " << v_cl.getProcessUnitID() << std::endl;
}
++it2;
}
}
openfpm::vector<couple_contrib> ctt;
double tot_sym = 0.0;
size_t tot_n_sym = 0;
// 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 symmetric
auto NpSym = NN.template getNNIterator<NO_CHECK>(p.getKey());
// For each neighborhood of the particle p
while (NpSym.isNext())
{
// Neighborhood particle q
auto q = NpSym.get();
// if p == q skip this particle
if (q == p.getKey()) {++NpSym; continue;};
// Get position of the particle q
Point<3,double> xq = vd.getPos(q);
// take the normalized direction
double rn = norm2(xp - xq);
////// debug ////////////
couple_contrib cc;
cc.i = vd.getProp<0>(p);
cc.j = vd.getProp<0>(q);
cc.ir = p.getKey();
cc.jr = q;
cc.E = rn;
cct.add(cc);
/////////////////////////
// potential energy (using pow is slower)
tot_sym += rn;
tot_n_sym++;
// Next neighborhood
++NpSym;
}
// Next Particle
++it2;
}
/////////////// NON SYMMETRIC CASE VERLET-LIST ////////////////////////
double tot = 0.0;
size_t tot_n = 0;
auto it = vd.getDomainIterator();
// For each particle ...
while (it.isNext())
{
// ... p
auto p = it.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 = NNc.template getNNIterator<NO_CHECK>(p.getKey());
// 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()) {++Np; continue;};
// Get position of the particle q
Point<3,double> xq = vd.getPos(q);
// take the normalized direction
double rn = norm2(xp - xq);
// potential energy (using pow is slower)
tot += rn/2.0;
tot_n++;
////// debug ////////////
couple_contrib cc;
cc.i = vd.getProp<0>(p);
cc.j = vd.getProp<0>(q);
cc.E = rn;
cct2.add(cc);
/////////////////////////
// Next neighborhood
++Np;
}
// Next Particle
++it;
}
Vcluster & v_cl = create_vcluster();
v_cl.sum(tot_n);
v_cl.sum(tot_n_sym);
v_cl.sum(tot);
v_cl.sum(tot_sym);
v_cl.execute();
openfpm::vector<couple_contrib> cc_collect;
openfpm::vector<couple_contrib> cc_collect2;
v_cl.SGather(cct,cc_collect,0);
v_cl.SGather(cct2,cc_collect2,0);
if (v_cl.getProcessUnitID() == 0)
{
cc_collect.sort();
cc_collect2.sort();
for (size_t i = 0 ; i < cc_collect2.size() ; i++)
{
if (i % 10000 == 0)
std::cout << i << " " << cc_collect2.size() << std::endl;
for (size_t k = i+1 ; k < i+2 && k < cc_collect2.size() ; k++)
{
if (cc_collect2.get(i).i == cc_collect2.get(k).j && cc_collect2.get(i).j == cc_collect2.get(k).i)
{
cc_collect2.remove(k);
}
}
}
for (size_t i = 0 ; i < cc_collect.size() && i < cc_collect2.size() ; i++)
std::cout << "E1: " << cc_collect.get(i).E << " " << cc_collect.get(i).i << "(" << cc_collect.get(i).ir << ")" << " " << cc_collect.get(i).j << "(" << cc_collect.get(i).jr << ")" << " E2: " << cc_collect2.get(i).E << " " << cc_collect2.get(i).i << " " << cc_collect2.get(i).j << std::endl;
}
// BOOST_REQUIRE_EQUAL(2*tot_n_sym,tot_n);
// BOOST_REQUIRE_CLOSE(tot,tot_sym,0.00001);
/* }*/
}
#endif /* SRC_VECTOR_VECTOR_DIST_CELL_LIST_TESTS_HPP_ */
This diff is collapsed.
......@@ -1573,6 +1573,124 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
}
}
BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
{
Vcluster & v_cl = create_vcluster();
long int k = 25*25*25*create_vcluster().getProcessingUnits();
k = std::pow(k, 1/3.);
if (v_cl.getProcessingUnits() > 48)
return;
print_test("Testing 3D periodic vector with ghost buffering k=",k);
BOOST_TEST_CHECKPOINT( "Testing 3D periodic with ghost buffering k=" << k );
long int big_step = k / 30;
big_step = (big_step == 0)?1:big_step;
long int small_step = 21;
// 3D test
for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
{
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));
typedef aggregate<float> part_prop;
// Distributed vector
vector_dist<3,float, part_prop > vd(0,box,bc,ghost);
auto it = vd.getGridIterator({k,k,k});
while (it.isNext())
{
auto key = it.get();
vd.add();
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);
// Fill some properties randomly
vd.getLastProp<0>() = 0.0;
++it;
}
vd.map();
// sync the ghost
vd.ghost_get<0>();
auto NN = vd.getCellList(1.3/k);
float a = 1.0f*k*k;
// run trough all the particles + ghost
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
// particle p
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
// Get an iterator over the neighborhood particles of p
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
// For each neighborhood particle ...
while (Np.isNext())
{
auto q = Np.get();
Point<3,float> xq = vd.getPos(q);
float dist = xp.distance(xq);
if (dist < 0.05)
vd.getProp<0>(q) += a*(-dist*dist+1.0/k/k);
++Np;
}
++it2;
}
vd.ghost_put<add_,0>();
bool ret = true;
auto it3 = vd.getDomainIterator();
float constant = vd.getProp<0>(it3.get());
float eps = 0.001;
while (it3.isNext())
{
float constant2 = vd.getProp<0>(it3.get());
if (fabs((constant - constant2)/constant > eps))
{
std::cout << Point<3,float>(vd.getPos(it3.get())).toString() << " " << constant2 << " " << v_cl.getProcessUnitID() << std::endl;
ret = false;
break;
}
++it3;
}
BOOST_REQUIRE_EQUAL(ret,true);
}
}
#include "vector_dist_cell_list_tests.hpp"
#include "vector_dist_NN_tests.hpp"
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment