OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
6constexpr int velocity = 0;
7constexpr 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
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
43int 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,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
60 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
61 Ghost<3,double> 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_frame("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).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}
This class represent an N-dimensional box.
Definition Box.hpp:61
Small class to produce graph with Google chart in HTML.
void write(std::string file)
It write the graphs on file in html format using Google charts.
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
void execute()
Execute all the requests.
void sum(T &num)
Sum the numbers across all processors and get the result.
size_t getProcessUnitID()
Get the process unit id.
Implementation of VCluster class.
Definition VCluster.hpp:59
Implementation of 1-D std::vector like structure.
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
Distributed vector.
Google chart options.
std::string xAxis
X axis name.
size_t lineWidth
Width of the line.
std::string title
Title of the chart.
std::string yAxis
Y axis name.