OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
main_expr.cpp
1 #include "Vector/vector_dist.hpp"
2 #include "Plot/GoogleChart.hpp"
3 #include "Operators/Vector/vector_dist_operators.hpp"
4 #include "timer.hpp"
5 
6 constexpr int velocity = 0;
7 constexpr int force = 1;
8 
10 {
11  double sigma12,sigma6,r_cut2,shift;
12 
13  ln_potential(double sigma12_, double sigma6_, double r_cut2_, double shift_) {sigma12 = sigma12_; sigma6 = sigma6_; r_cut2 = r_cut2_;shift = shift_;}
14 
15  Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
16  {
17  double rn = norm2(xp - xq);
18  if (rn >= r_cut2) return 0.0;
19 
20  Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift,0.0});
21 
22  return E;
23  }
24 };
25 
26 struct ln_force
27 {
28  double sigma12,sigma6,r_cut2;
29 
30  ln_force(double sigma12_, double sigma6_, double r_cut2_) {sigma12 = sigma12_; sigma6 = sigma6_;r_cut2 = r_cut2_;}
31 
32  Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
33  {
34  Point<3,double> r = xp - xq;
35  double rn = norm2(r);
36 
37  if (rn > r_cut2) return 0.0;
38 
39  return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
40  }
41 };
42 
43 int main(int argc, char* argv[])
44 {
45  double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
46 
47  double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
48  double rc2 = r_cut * r_cut;
49  double shift = 2.0 * ( sigma12 / (rc2*rc2*rc2*rc2*rc2*rc2) - sigma6 / ( rc2*rc2*rc2) );
50 
53 
54 
55  openfpm_init(&argc,&argv);
56  Vcluster & vcl = create_vcluster();
57 
58  size_t sz[3] = {10,10,10};
59  Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
60  size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
61  Ghost<3,float> ghost(r_cut);
62 
63  ln_force lf(sigma12,sigma6,r_cut*r_cut);
64  ln_potential lp(sigma12,sigma6,r_cut*r_cut,shift);
65 
67 
68  auto v_force = getV<force>(vd);
69  auto v_velocity = getV<velocity>(vd);
70  auto v_pos = getV<PROP_POS>(vd);
71 
72  auto it = vd.getGridIterator(sz);
73 
74  while (it.isNext())
75  {
76  vd.add();
77 
78  auto key = it.get();
79 
80  vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
81  vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
82  vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
83 
84  ++it;
85  }
86 
87  v_force = 0;
88  v_velocity = 0;
89 
90  timer tsim;
91  tsim.start();
92 
93  auto NN = vd.getCellList(r_cut);
94 
95  vd.updateCellList(NN);
96  v_force = applyKernel_in_sim(vd,NN,lf);
97  unsigned long int f = 0;
98 
99  // MD time stepping
100  for (size_t i = 0; i < 10000 ; i++)
101  {
102  assign(v_velocity, v_velocity + 0.5*dt*v_force,
103  v_pos, v_pos + v_velocity*dt);
104 
105  vd.map();
106  vd.template ghost_get<>();
107 
108  // calculate forces or a(tn + 1) Step 2
109  vd.updateCellList(NN);
110  v_force = applyKernel_in_sim(vd,NN,lf);
111 
112  v_velocity = v_velocity + 0.5*dt*v_force;
113 
114  if (i % 100 == 0)
115  {
116  vd.deleteGhost();
117  vd.write("particles_",f);
118  vd.ghost_get<>();
119 
120  vd.updateCellList(NN);
121  Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
122 
123  vcl.sum(E.get(0));vcl.sum(E.get(1));
124  vcl.execute();
125 
126  // we save the energy calculated at time step i c contain the time-step y contain the energy
127  x.add(i);
128  y.add({E.get(0),E.get(1),E.get(0) - E.get(1)});
129 
130  if (vcl.getProcessUnitID() == 0)
131  std::cout << "Energy Total: " << E.get(0) << " Kinetic: " << E.get(1) << " Potential: " << E.get(0) - E.get(1) << std::endl;
132 
133  f++;
134  }
135  }
136 
137  tsim.stop();
138  std::cout << "Time: " << tsim.getwct() << std::endl;
139 
140  GCoptions options;
141  options.title = std::string("Energy with time");
142  options.yAxis = std::string("Energy");
143  options.xAxis = std::string("iteration");
144  options.lineWidth = 1.0;
145 
146  GoogleChart cg;
147  cg.AddLinesGraph(x,y,options);
148  cg.write("gc_plot2_out.html");
149 
150  openfpm_finalize();
151 }
void sum(T &num)
Sum the numbers across all processors and get the result.
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
size_t getProcessUnitID()
Get the process unit id.
void execute()
Execute all the requests.
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
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.
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