OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
main_comp_ord.cpp
1 
2 #include "Vector/vector_dist.hpp"
3 #include "Decomposition/CartDecomposition.hpp"
4 #include "data_type/aggregate.hpp"
5 #include "Plot/GoogleChart.hpp"
6 #include "Plot/util.hpp"
7 #include "timer.hpp"
8 
41 
43 constexpr int velocity = 0;
44 constexpr int force = 1;
45 
47 
48 template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
49 {
50 
52 
53 
54  // Uodate the cell-list
55  vd.updateCellList(NN);
56 
58 
59  // Get an iterator over particles
60  auto it2 = NN.getIterator();
61 
63 
64  // For each particle p ...
65  while (it2.isNext())
66  {
67  // ... get the particle p
68  auto p = it2.get();
69 
70  // Get the position xp of the particle
71  Point<3,double> xp = vd.getPos(p);
72 
73  // Reset the force counter
74  vd.template getProp<force>(p)[0] = 0.0;
75  vd.template getProp<force>(p)[1] = 0.0;
76  vd.template getProp<force>(p)[2] = 0.0;
77 
78  // Get an iterator over the neighborhood particles of p
79  auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
80 
81  // For each neighborhood particle ...
82  while (Np.isNext())
83  {
84  // ... q
85  auto q = Np.get();
86 
87  // if (p == q) skip this particle
88  if (q == p) {++Np; continue;};
89 
90  // Get the position of p
91  Point<3,double> xq = vd.getPos(q);
92 
93  // Get the distance between p and q
94  Point<3,double> r = xp - xq;
95 
96  // take the norm of this vector
97  double rn = norm2(r);
98 
99  // Calculate the force, using pow is slower
100  Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
101 
102  // we sum the force produced by q on p
103  vd.template getProp<force>(p)[0] += f.get(0);
104  vd.template getProp<force>(p)[1] += f.get(1);
105  vd.template getProp<force>(p)[2] += f.get(2);
106 
107  // Next neighborhood
108  ++Np;
109  }
110 
111  // Next particle
112  ++it2;
113  }
114 }
115 
117 
141 
143 template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
144 {
145  double E = 0.0;
146 
147  // update cell-list
148  vd.updateCellList(NN);
149 
150  // Get the iterator
151  auto it2 = NN.getIterator();
152 
153  // For each particle ...
154  while (it2.isNext())
155  {
156  // ... p
157  auto p = it2.get();
158 
159  // Get the position of the particle p
160  Point<3,double> xp = vd.getPos(p);
161 
162  // Reset the force
163  vd.template getProp<force>(p)[0] = 0.0;
164  vd.template getProp<force>(p)[1] = 0.0;
165  vd.template getProp<force>(p)[2] = 0.0;
166 
167  // Get an iterator over the neighborhood of the particle p
168  auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
169 
170  // For each neighborhood of the particle p
171  while (Np.isNext())
172  {
173  // Neighborhood particle q
174  auto q = Np.get();
175 
176  // if p == q skip this particle
177  if (q == p) {++Np; continue;};
178 
179  // Get position of the particle q
180  Point<3,double> xq = vd.getPos(q);
181 
182  // take the normalized direction
183  double rn = norm2(xp - xq);
184 
185  // potential energy (using pow is slower)
186  E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
187 
188  // Next neighborhood
189  ++Np;
190  }
191 
192  // Kinetic energy of the particle given by its actual speed
193  E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
194  vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
195  vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
196 
197  // Next Particle
198  ++it2;
199  }
200 
201  // Calculated energy
202  return E;
203 }
204 
206 
207 
208 int main(int argc, char* argv[])
209 {
226 
228  double dt = 0.0001;
229  float r_cut = 0.03;
230  double sigma = r_cut/3.0;
231  double sigma12 = pow(sigma,12);
232  double sigma6 = pow(sigma,6);
233 
236 
237  openfpm_init(&argc,&argv);
238  Vcluster & v_cl = create_vcluster();
239 
240  // we will use it do place particles on a 40x40x40 Grid like
241  size_t sz[3] = {40,40,40};
242 
243  // domain
244  Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
245 
246  // Boundary conditions
247  size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
248 
249  // ghost, big enough to contain the interaction radius
250  Ghost<3,float> ghost(r_cut);
251 
253 
255 
269 
271  auto it = vd.getGridIterator(sz);
272 
273  while (it.isNext())
274  {
275  vd.add();
276 
277  auto key = it.get();
278 
279  vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
280  vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
281  vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
282 
283  vd.template getLastProp<velocity>()[0] = 0.0;
284  vd.template getLastProp<velocity>()[1] = 0.0;
285  vd.template getLastProp<velocity>()[2] = 0.0;
286 
287  vd.template getLastProp<force>()[0] = 0.0;
288  vd.template getLastProp<force>()[1] = 0.0;
289  vd.template getLastProp<force>()[2] = 0.0;
290 
291  ++it;
292  }
293 
295 
412 
414  // Get the Cell list structure
415  auto NN = vd.getCellList_hilb(r_cut);
416 
417  // calculate forces
418  calc_forces(vd,NN,sigma12,sigma6);
419  unsigned long int f = 0;
420 
421  timer time2;
422  time2.start();
423 
424  #ifndef TEST_RUN
425  size_t Nstep = 30000;
426  #else
427  size_t Nstep = 300;
428  #endif
429 
430  // MD time stepping
431  for (size_t i = 0; i < Nstep ; i++)
432  {
433  // Get the iterator
434  auto it3 = vd.getDomainIterator();
435 
436  // integrate velicity and space based on the calculated forces (Step1)
437  while (it3.isNext())
438  {
439  auto p = it3.get();
440 
441  // here we calculate v(tn + 0.5)
442  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
443  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
444  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
445 
446  // here we calculate x(tn + 1)
447  vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
448  vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
449  vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
450 
451  ++it3;
452  }
453 
454  // Because we mooved the particles in space we have to map them and re-sync the ghost
455  vd.map();
456  vd.template ghost_get<>();
457 
458  timer time;
459  if (i % 10 == 0)
460  time.start();
461 
462  // calculate forces or a(tn + 1) Step 2
463  calc_forces(vd,NN,sigma12,sigma6);
464 
465  if (i % 10 == 0)
466  {
467  time.stop();
468  x.add(i);
469  y.add({time.getwct()});
470  }
471 
472  // Integrate the velocity Step 3
473  auto it4 = vd.getDomainIterator();
474 
475  while (it4.isNext())
476  {
477  auto p = it4.get();
478 
479  // here we calculate v(tn + 1)
480  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
481  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
482  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
483 
484  ++it4;
485  }
486 
487  // After every iteration collect some statistic about the confoguration
488  if (i % 100 == 0)
489  {
490  // We write the particle position for visualization (Without ghost)
491  vd.deleteGhost();
492  vd.write("particles_",f);
493 
494  // we resync the ghost
495  vd.ghost_get<>();
496 
497  // We calculate the energy
498  double energy = calc_energy(vd,NN,sigma12,sigma6);
499  auto & vcl = create_vcluster();
500  vcl.sum(energy);
501  vcl.execute();
502 
503  // We also print on terminal the value of the energy
504  // only one processor (master) write on terminal
505  if (vcl.getProcessUnitID() == 0)
506  std::cout << std::endl << "Energy: " << energy << std::endl;
507 
508  f++;
509  }
510  }
511 
512  time2.stop();
513  std::cout << "Performance: " << time2.getwct() << std::endl;
514 
516 
532 
534  // Google charts options, it store the options to draw the X Y graph
535  GCoptions options;
536 
537  // Title of the graph
538  options.title = std::string("Force calculation time");
539 
540  // Y axis name
541  options.yAxis = std::string("Time");
542 
543  // X axis name
544  options.xAxis = std::string("iteration");
545 
546  // width of the line
547  options.lineWidth = 1.0;
548 
549  // Object that draw the X Y graph
550  GoogleChart cg;
551 
552  // Add the graph
553  // The graph that it produce is in svg format that can be opened on browser
554  cg.AddLinesGraph(x,y,options);
555 
556  // Write into html format
557  cg.write("gc_plot2_out.html");
558 
560 
573 
575  openfpm_finalize();
576 
578 
588 }
589 
590 
591 
592 
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
std::string title
Title of the chart.
Definition: GoogleChart.hpp:27
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
auto get(size_t cell, size_t ele) -> decltype(this->Mem_type::get(cell, ele))
Get an element in the cell.
Definition: CellList.hpp:794
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
Class for FAST cell list implementation.
Definition: CellList.hpp:269
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