OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
main_data_ord.cpp
1 #include "Vector/vector_dist.hpp"
2 #include "data_type/aggregate.hpp"
3 #include "Plot/GoogleChart.hpp"
4 #include "Plot/util.hpp"
5 #include "timer.hpp"
6 #include "energy_force.hpp"
7 
22 
24 int main(int argc, char* argv[])
25 {
41 
43  double dt = 0.0001;
44  float r_cut = 0.03;
45  double sigma = r_cut/3.0;
46  double sigma12 = pow(sigma,12);
47  double sigma6 = pow(sigma,6);
48 
51 
52  openfpm_init(&argc,&argv);
53  Vcluster & v_cl = create_vcluster();
54 
55  // we will use it do place particles on a 10x10x10 Grid like
56  size_t sz[3] = {40,40,40};
57 
58  // domain
59  Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
60 
61  // Boundary conditions
62  size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
63 
64  // ghost, big enough to contain the interaction radius
65  Ghost<3,float> ghost(r_cut);
66 
67  // Create vector
69 
71 
85 
87  auto it = vd.getGridIterator(sz);
88 
89  while (it.isNext())
90  {
91  vd.add();
92 
93  auto key = it.get();
94 
95  vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
96  vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
97  vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
98 
99  vd.template getLastProp<velocity>()[0] = 0.0;
100  vd.template getLastProp<velocity>()[1] = 0.0;
101  vd.template getLastProp<velocity>()[2] = 0.0;
102 
103  vd.template getLastProp<force>()[0] = 0.0;
104  vd.template getLastProp<force>()[1] = 0.0;
105  vd.template getLastProp<force>()[2] = 0.0;
106 
107  ++it;
108  }
109 
111 
268 
270  // Get the Cell list structure
271  auto NN = vd.getCellList(r_cut);
272 
273  // calculate forces
274  calc_forces(vd,NN,sigma12,sigma6);
275  unsigned long int f = 0;
276 
277  timer time2;
278  time2.start();
279 
280  #ifndef TEST_RUN
281  size_t Nstep = 30000;
282  #else
283  size_t Nstep = 300;
284  #endif
285 
286  // MD time stepping
287  for (size_t i = 0; i < Nstep ; i++)
288  {
289  // Get the iterator
290  auto it3 = vd.getDomainIterator();
291 
292  // integrate velicity and space based on the calculated forces (Step1)
293  while (it3.isNext())
294  {
295  auto p = it3.get();
296 
297  // here we calculate v(tn + 0.5)
298  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
299  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
300  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
301 
302  // here we calculate x(tn + 1)
303  vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
304  vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
305  vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
306 
307  ++it3;
308  }
309 
310  // Because we mooved the particles in space we have to map them and re-sync the ghost
311  vd.map();
312  vd.template ghost_get<>();
313 
315 
316  if (i % 200 == 0)
317  vd.reorder(5);
318 
320 
322 
323  timer time;
324  if (i % 10 == 0)
325  time.start();
326 
328 
329  // calculate forces or a(tn + 1) Step 2
330  calc_forces(vd,NN,sigma12,sigma6);
331 
333 
334  if (i % 10 == 0)
335  {
336  time.stop();
337  x.add(i);
338  y.add({time.getwct()});
339  }
340 
342 
343  // Integrate the velocity Step 3
344  auto it4 = vd.getDomainIterator();
345 
346  while (it4.isNext())
347  {
348  auto p = it4.get();
349 
350  // here we calculate v(tn + 1)
351  vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
352  vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
353  vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
354 
355  ++it4;
356  }
357 
358  // After every iteration collect some statistic about the confoguration
359  if (i % 100 == 0)
360  {
361  // We calculate the energy
362  double energy = calc_energy(vd,NN,sigma12,sigma6);
363  auto & vcl = create_vcluster();
364  vcl.sum(energy);
365  vcl.execute();
366 
367  // We also print on terminal the value of the energy
368  // only one processor (master) write on terminal
369  if (vcl.getProcessUnitID() == 0)
370  std::cout << "Energy: " << energy << std::endl;
371 
372  f++;
373  }
374  }
375 
376  time2.stop();
377  std::cout << "Performance: " << time2.getwct() << std::endl;
378 
380 
395 
397  // Google charts options, it store the options to draw the X Y graph
398  GCoptions options;
399 
400  // Title of the graph
401  options.title = std::string("Force calculation time");
402 
403  // Y axis name
404  options.yAxis = std::string("Time");
405 
406  // X axis name
407  options.xAxis = std::string("iteration");
408 
409  // width of the line
410  options.lineWidth = 1.0;
411 
412  // Object that draw the X Y graph
413  GoogleChart cg;
414 
415  // Add the graph
416  // The graph that it produce is in svg format that can be opened on browser
417  cg.AddLinesGraph(x,y,options);
418 
419  // Write into html format
420  cg.write("gc_plot2_out.html");
421 
423 
435 
437  openfpm_finalize();
438 
440 
449 }
450 
451 
452 
453 
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
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
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.
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