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 
95 
97 constexpr int velocity = 0;
98 constexpr int force = 1;
99 
101 
128 
130 template<typename VerletList>
131 void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList & NN, double sigma12, double sigma6, double r_cut)
132 {
133  // Reset force on the ghost
134 
135  auto itg = vd.getDomainAndGhostIterator();
136 
137  while (itg.isNext())
138  {
139  auto p = itg.get();
140 
141  // Reset force
142  vd.getProp<force>(p)[0] = 0.0;
143  vd.getProp<force>(p)[1] = 0.0;
144  vd.getProp<force>(p)[2] = 0.0;
145 
146  ++itg;
147  }
148 
150 
151  // Get an iterator over particles
152  auto it2 = vd.getParticleIteratorCRS(NN);
153 
155 
156  // For each particle p ...
157  while (it2.isNext())
158  {
159  // ... get the particle p
160  auto p = it2.get();
161 
162  // Get the position xp of the particle
163  Point<3,double> xp = vd.getPos(p);
164 
165  // Get an iterator over the neighborhood particles of p
166  // Note that in case of symmetric
167  auto Np = NN.template getNNIterator<NO_CHECK>(p);
168 
169  // For each neighborhood particle ...
170  while (Np.isNext())
171  {
172  // ... q
173  auto q = Np.get();
174 
175  // if (p == q) skip this particle
176  if (q == p) {++Np; continue;};
177 
178  // Get the position of q
179  Point<3,double> xq = vd.getPos(q);
180 
181  // Get the distance between p and q
182  Point<3,double> r = xp - xq;
183 
184  // take the norm of this vector
185  double rn = norm2(r);
186 
187  if (rn > r_cut * r_cut) {++Np;continue;}
188 
189  // Calculate the force, using pow is slower
190  Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
191 
192  // we sum the force produced by q on p
193  vd.template getProp<force>(p)[0] += f.get(0);
194  vd.template getProp<force>(p)[1] += f.get(1);
195  vd.template getProp<force>(p)[2] += f.get(2);
196 
198 
199  // we sum the force produced by p on q
200  vd.template getProp<force>(q)[0] -= f.get(0);
201  vd.template getProp<force>(q)[1] -= f.get(1);
202  vd.template getProp<force>(q)[2] -= f.get(2);
203 
205 
206  // Next neighborhood
207  ++Np;
208  }
209 
210  // Next particle
211  ++it2;
212  }
213 
215 
216  // Sum the contribution to the real particles
217  vd.ghost_put<add_,force>();
218 
220 }
221 
223 
224 
244 
246 template<typename VerletList>
247 double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList & NN, double sigma12, double sigma6, double r_cut)
248 {
249  double E = 0.0;
250 
251  double rc = r_cut*r_cut;
252  double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
253 
254  // Get an iterator over particles
255  auto it2 = vd.getParticleIteratorCRS(NN);
256 
257  // For each particle ...
258  while (it2.isNext())
259  {
260  // ... p
261  auto p = it2.get();
262 
263  // Get the position of the particle p
264  Point<3,double> xp = vd.getPos(p);
265 
266  // Get an iterator over the neighborhood particles of p
267  auto Np = NN.template getNNIterator<NO_CHECK>(p);
268 
269  double Ep = E;
270 
271  // For each neighborhood of the particle p
272  while (Np.isNext())
273  {
274  // Neighborhood particle q
275  auto q = Np.get();
276 
277  // if p == q skip this particle
278  if (q == p) {++Np; continue;};
279 
280  // Get position of the particle q
281  Point<3,double> xq = vd.getPos(q);
282 
283  // take the normalized direction
284  double rn = norm2(xp - xq);
285 
286  if (rn >= r_cut*r_cut) {++Np;continue;}
287 
288  // potential energy (using pow is slower)
289  E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
290 
291  // Next neighborhood
292  ++Np;
293  }
294 
295  // To note that the crossing scheme go across the domain particles +
296  // some ghost particles. This mean that we have to filter out the ghost
297  // particles otherwise we count double energies
298  //
299  if (p < vd.size_local())
300  {
301  // Kinetic energy of the particle given by its actual speed
302  E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
303  vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
304  vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
305  }
306 
307  // Next Particle
308  ++it2;
309  }
310 
311  // Calculated energy
312  return E;
313 }
314 
316 
317 int main(int argc, char* argv[])
318 {
347 
349  double dt = 0.00025;
350  double sigma = 0.1;
351  double r_cut = 3.0*sigma;
352  double r_gskin = 1.3*r_cut;
353  double sigma12 = pow(sigma,12);
354  double sigma6 = pow(sigma,6);
355 
358 
359  openfpm_init(&argc,&argv);
360  Vcluster & v_cl = create_vcluster();
361 
362  // we will use it do place particles on a 10x10x10 Grid like
363  size_t sz[3] = {10,10,10};
364 
365  // domain
366  Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
367 
368  // Boundary conditions
369  size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
370 
372 
373  // ghost, big enough to contain the interaction radius
374  Ghost<3,float> ghost(r_gskin);
375  ghost.setLow(0,0.0);
376  ghost.setLow(1,0.0);
377  ghost.setLow(2,0.0);
378 
380 
382 
383  vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost,BIND_DEC_TO_GHOST);
384 
386 
387  size_t k = 0;
388  size_t start = vd.accum();
389 
390  auto it = vd.getGridIterator(sz);
391 
392  while (it.isNext())
393  {
394  vd.add();
395 
396  auto key = it.get();
397 
398  vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
399  vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
400  vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
401 
402  vd.template getLastProp<velocity>()[0] = 0.0;
403  vd.template getLastProp<velocity>()[1] = 0.0;
404  vd.template getLastProp<velocity>()[2] = 0.0;
405 
406  vd.template getLastProp<force>()[0] = 0.0;
407  vd.template getLastProp<force>()[1] = 0.0;
408  vd.template getLastProp<force>()[2] = 0.0;
409 
410  k++;
411  ++it;
412  }
413 
414  vd.map();
415  vd.ghost_get<>();
416 
417  timer tsim;
418  tsim.start();
419 
421 
422  // Get the Cell list structure
423  auto NN = vd.getVerletCrs(r_gskin);;
424 
426 
427  // calculate forces
428  calc_forces(vd,NN,sigma12,sigma6,r_cut);
429  unsigned long int f = 0;
430 
431  int cnt = 0;
432  double max_disp = 0.0;
433 
434  // MD time stepping
435  for (size_t i = 0; i < 10000 ; i++)
436  {
437  // Get the iterator
438  auto it3 = vd.getDomainIterator();
439 
440  double max_displ = 0.0;
441 
442  // integrate velicity and space based on the calculated forces (Step1)
443  while (it3.isNext())
444  {
445  auto p = it3.get();
446 
447  // here we calculate v(tn + 0.5)
448  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
449  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
450  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
451 
452  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});
453 
454  // here we calculate x(tn + 1)
455  vd.getPos(p)[0] += disp.get(0);
456  vd.getPos(p)[1] += disp.get(1);
457  vd.getPos(p)[2] += disp.get(2);
458 
459  if (disp.norm() > max_displ)
460  max_displ = disp.norm();
461 
462  ++it3;
463  }
464 
465  if (max_disp < max_displ)
466  max_disp = max_displ;
467 
468  // Because we moved the particles in space we have to map them and re-sync the ghost
469  if (cnt % 10 == 0)
470  {
471  vd.map();
472  vd.template ghost_get<>();
473  // Get the Cell list structure
474  vd.updateVerlet(NN,r_gskin,VL_CRS_SYMMETRIC);
475  }
476  else
477  {
478  vd.template ghost_get<>(SKIP_LABELLING);
479  }
480 
481  cnt++;
482 
483  // calculate forces or a(tn + 1) Step 2
484  calc_forces(vd,NN,sigma12,sigma6,r_cut);
485 
486  // Integrate the velocity Step 3
487  auto it4 = vd.getDomainIterator();
488 
489  while (it4.isNext())
490  {
491  auto p = it4.get();
492 
493  // here we calculate v(tn + 1)
494  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
495  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
496  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
497 
498  ++it4;
499  }
500 
501  // After every iteration collect some statistic about the confoguration
502  if (i % 100 == 0)
503  {
504  // We write the particle position for visualization (Without ghost)
505  vd.deleteGhost();
506  vd.write("particles_",f);
507 
508  // we resync the ghost
509  vd.ghost_get<>();
510 
511 
512  // We calculate the energy
513  double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut);
514  auto & vcl = create_vcluster();
515  vcl.sum(energy);
516  vcl.max(max_disp);
517  vcl.execute();
518 
519  // we save the energy calculated at time step i c contain the time-step y contain the energy
520  x.add(i);
521  y.add({energy});
522 
523  // We also print on terminal the value of the energy
524  // only one processor (master) write on terminal
525  if (vcl.getProcessUnitID() == 0)
526  std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
527 
528  max_disp = 0.0;
529 
530  f++;
531  }
532  }
533 
534  tsim.stop();
535  std::cout << "Time: " << tsim.getwct() << std::endl;
536 
538 
539  // Google charts options, it store the options to draw the X Y graph
540  GCoptions options;
541 
542  // Title of the graph
543  options.title = std::string("Energy with time");
544 
545  // Y axis name
546  options.yAxis = std::string("Energy");
547 
548  // X axis name
549  options.xAxis = std::string("iteration");
550 
551  // width of the line
552  options.lineWidth = 1.0;
553 
554  // Object that draw the X Y graph
555  GoogleChart cg;
556 
557  // Add the graph
558  // The graph that it produce is in svg format that can be opened on browser
559  cg.AddLinesGraph(x,y,options);
560 
561  // Write into html format
562  cg.write("gc_plot2_out.html");
563 
565 
577 
579  openfpm_finalize();
580 
582 
591 }
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