OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
main.cpp
1 #include "Vector/vector_dist.hpp"
2 #include "Decomposition/CartDecomposition.hpp"
3 #include "data_type/aggregate.hpp"
4 #include "Plot/GoogleChart.hpp"
5 #include "Plot/util.hpp"
6 #include "timer.hpp"
7 
34 
36 constexpr int velocity = 0;
37 constexpr int force = 1;
38 
40 
77 
79 template<typename VerletList>
80 void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList & NN, double sigma12, double sigma6, double r_cut)
81 {
82  // Reset force on the ghost
83 
84  auto itg = vd.getDomainAndGhostIterator();
85 
86  while (itg.isNext())
87  {
88  auto p = itg.get();
89 
90  // Reset force
91  vd.getProp<force>(p)[0] = 0.0;
92  vd.getProp<force>(p)[1] = 0.0;
93  vd.getProp<force>(p)[2] = 0.0;
94 
95  ++itg;
96  }
97 
99 
100  // Get an iterator over particles
101  auto it2 = vd.getDomainIterator();
102 
104 
105  // For each particle p ...
106  while (it2.isNext())
107  {
108  // ... get the particle p
109  auto p = it2.get();
110 
111  // Get the position xp of the particle
112  Point<3,double> xp = vd.getPos(p);
113 
114  // Get an iterator over the neighborhood particles of p
115  // Note that in case of symmetric
116  auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
117 
118  // For each neighborhood particle ...
119  while (Np.isNext())
120  {
121  // ... q
122  auto q = Np.get();
123 
124  // if (p == q) skip this particle
125  if (q == p.getKey()) {++Np; continue;};
126 
127  // Get the position of q
128  Point<3,double> xq = vd.getPos(q);
129 
130  // Get the distance between p and q
131  Point<3,double> r = xp - xq;
132 
133  // take the norm of this vector
134  double rn = norm2(r);
135 
136  if (rn > r_cut * r_cut) {++Np;continue;}
137 
138  // Calculate the force, using pow is slower
139  Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
140 
141  // we sum the force produced by q on p
142  vd.template getProp<force>(p)[0] += f.get(0);
143  vd.template getProp<force>(p)[1] += f.get(1);
144  vd.template getProp<force>(p)[2] += f.get(2);
145 
147 
148  // we sum the force produced by p on q
149  vd.template getProp<force>(q)[0] -= f.get(0);
150  vd.template getProp<force>(q)[1] -= f.get(1);
151  vd.template getProp<force>(q)[2] -= f.get(2);
152 
154 
155  // Next neighborhood
156  ++Np;
157  }
158 
159  // Next particle
160  ++it2;
161  }
162 
164 
165  // Sum the contribution to the real particles
166  vd.ghost_put<add_,force>();
167 
169 }
170 
172 
173 
188 
190 template<typename VerletList>
191 double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList & NN, double sigma12, double sigma6, double r_cut)
192 {
193  double E = 0.0;
194 
195  double rc = r_cut*r_cut;
196  double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
197 
198  // Get the iterator
199  auto it2 = vd.getDomainIterator();
200 
201  // For each particle ...
202  while (it2.isNext())
203  {
204  // ... p
205  auto p = it2.get();
206 
207  // Get the position of the particle p
208  Point<3,double> xp = vd.getPos(p);
209 
210  // Get an iterator over the neighborhood of the particle p
211  auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
212 
213  double Ep = E;
214 
215  // For each neighborhood of the particle p
216  while (Np.isNext())
217  {
218  // Neighborhood particle q
219  auto q = Np.get();
220 
221  // if p == q skip this particle
222  if (q == p.getKey()) {++Np; continue;};
223 
224  // Get position of the particle q
225  Point<3,double> xq = vd.getPos(q);
226 
227  // take the normalized direction
228  double rn = norm2(xp - xq);
229 
230  if (rn >= r_cut*r_cut) {++Np;continue;}
231 
232  // potential energy (using pow is slower)
233  E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
234 
235  // Next neighborhood
236  ++Np;
237  }
238 
239  // Kinetic energy of the particle given by its actual speed
240  E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
241  vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
242  vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
243 
244  // Next Particle
245  ++it2;
246  }
247 
248  // Calculated energy
249  return E;
250 }
251 
253 
254 int main(int argc, char* argv[])
255 {
274 
276  double dt = 0.00025;
277  double sigma = 0.1;
278  double r_cut = 3.0*sigma;
279  double r_gskin = 1.3*r_cut;
280  double sigma12 = pow(sigma,12);
281  double sigma6 = pow(sigma,6);
282 
285 
286  openfpm_init(&argc,&argv);
287  Vcluster & v_cl = create_vcluster();
288 
289  // we will use it do place particles on a 10x10x10 Grid like
290  size_t sz[3] = {10,10,10};
291 
292  // domain
293  Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
294 
295  // Boundary conditions
296  size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
297 
298  // ghost, big enough to contain the interaction radius
299  Ghost<3,float> ghost(r_gskin);
300 
302 
303  size_t k = 0;
304  size_t start = vd.accum();
305 
306  auto it = vd.getGridIterator(sz);
307 
308  while (it.isNext())
309  {
310  vd.add();
311 
312  auto key = it.get();
313 
314  vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
315  vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
316  vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
317 
318  vd.template getLastProp<velocity>()[0] = 0.0;
319  vd.template getLastProp<velocity>()[1] = 0.0;
320  vd.template getLastProp<velocity>()[2] = 0.0;
321 
322  vd.template getLastProp<force>()[0] = 0.0;
323  vd.template getLastProp<force>()[1] = 0.0;
324  vd.template getLastProp<force>()[2] = 0.0;
325 
326  k++;
327  ++it;
328  }
329 
330  vd.map();
331  vd.ghost_get<>();
332 
333  timer tsim;
334  tsim.start();
335 
337 
338  // Get the Cell list structure
339  auto NN = vd.getVerletSym(r_gskin);
340 
342 
343  // calculate forces
344  calc_forces(vd,NN,sigma12,sigma6,r_cut);
345  unsigned long int f = 0;
346 
347  int cnt = 0;
348  double max_disp = 0.0;
349 
350  // MD time stepping
351  for (size_t i = 0; i < 10000 ; i++)
352  {
353  // Get the iterator
354  auto it3 = vd.getDomainIterator();
355 
356  double max_displ = 0.0;
357 
358  // integrate velicity and space based on the calculated forces (Step1)
359  while (it3.isNext())
360  {
361  auto p = it3.get();
362 
363  // here we calculate v(tn + 0.5)
364  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
365  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
366  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
367 
368  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});
369 
370  // here we calculate x(tn + 1)
371  vd.getPos(p)[0] += disp.get(0);
372  vd.getPos(p)[1] += disp.get(1);
373  vd.getPos(p)[2] += disp.get(2);
374 
375  if (disp.norm() > max_displ)
376  max_displ = disp.norm();
377 
378  ++it3;
379  }
380 
381  if (max_disp < max_displ)
382  max_disp = max_displ;
383 
384  // Because we moved the particles in space we have to map them and re-sync the ghost
385  if (cnt % 10 == 0)
386  {
387  vd.map();
388  vd.template ghost_get<>();
389  // Get the Cell list structure
390  vd.updateVerlet(NN,r_gskin,VL_SYMMETRIC);
391  }
392  else
393  {
394  vd.template ghost_get<>(SKIP_LABELLING);
395  }
396 
397  cnt++;
398 
399  // calculate forces or a(tn + 1) Step 2
400  calc_forces(vd,NN,sigma12,sigma6,r_cut);
401 
402  // Integrate the velocity Step 3
403  auto it4 = vd.getDomainIterator();
404 
405  while (it4.isNext())
406  {
407  auto p = it4.get();
408 
409  // here we calculate v(tn + 1)
410  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
411  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
412  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
413 
414  ++it4;
415  }
416 
417  // After every iteration collect some statistic about the confoguration
418  if (i % 100 == 0)
419  {
420  // We write the particle position for visualization (Without ghost)
421  vd.deleteGhost();
422  vd.write("particles_",f);
423 
424  // we resync the ghost
425  vd.ghost_get<>();
426 
427 
428  // We calculate the energy
429  double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut);
430  auto & vcl = create_vcluster();
431  vcl.sum(energy);
432  vcl.max(max_disp);
433  vcl.execute();
434 
435  // we save the energy calculated at time step i c contain the time-step y contain the energy
436  x.add(i);
437  y.add({energy});
438 
439  // We also print on terminal the value of the energy
440  // only one processor (master) write on terminal
441  if (vcl.getProcessUnitID() == 0)
442  std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
443 
444  max_disp = 0.0;
445 
446  f++;
447  }
448  }
449 
450  tsim.stop();
451  std::cout << "Time: " << tsim.getwct() << std::endl;
452 
454 
455  // Google charts options, it store the options to draw the X Y graph
456  GCoptions options;
457 
458  // Title of the graph
459  options.title = std::string("Energy with time");
460 
461  // Y axis name
462  options.yAxis = std::string("Energy");
463 
464  // X axis name
465  options.xAxis = std::string("iteration");
466 
467  // width of the line
468  options.lineWidth = 1.0;
469 
470  // Object that draw the X Y graph
471  GoogleChart cg;
472 
473  // Add the graph
474  // The graph that it produce is in svg format that can be opened on browser
475  cg.AddLinesGraph(x,y,options);
476 
477  // Write into html format
478  cg.write("gc_plot2_out.html");
479 
481 
493 
495  openfpm_finalize();
496 
498 
507 }
size_t get(size_t i, size_t j) const
Get the neighborhood element j for the particle i.
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
Class for Verlet list implementation.
std::string title
Title of the chart.
Definition: GoogleChart.hpp:27
This structure define the operation add to use with copy general.
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:22
double getwct()
Return the elapsed real time.
Definition: timer.hpp:108
Definition: Ghost.hpp:39
Implementation of VCluster class.
Definition: VCluster.hpp:36
std::string yAxis
Y axis name.
Definition: GoogleChart.hpp:29
Small class to produce graph with Google chart in HTML.
void start()
Start the timer.
Definition: timer.hpp:73
const T & get(size_t i) const
Get coordinate.
Definition: Point.hpp:142
size_t lineWidth
Width of the line.
Definition: GoogleChart.hpp:55
This class is a trick to indicate the compiler a specific specialization pattern. ...
Definition: memory_c.hpp:201
void write(std::string file)
It write the graphs on file in html format using Google charts.
Distributed vector.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:81
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:61
std::string xAxis
X axis name.
Definition: GoogleChart.hpp:31
Google chart options.
Definition: GoogleChart.hpp:24
Class for cpu time benchmarking.
Definition: timer.hpp:25
void stop()
Stop the timer.
Definition: timer.hpp:97