OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
main_expr_paper.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
7void Init_grid(size_t (& sz)[3],
9{
10 auto it = particles.getGridIterator(sz);
11
12 while (it.isNext())
13 {
14 particles.add();
15
16 auto key = it.get();
17
18 particles.getLastPos()[0] = key.get(0) * it.getSpacing(0);
19 particles.getLastPos()[1] = key.get(1) * it.getSpacing(1);
20 particles.getLastPos()[2] = key.get(2) * it.getSpacing(2);
21
22 ++it;
23 }
24}
25
26double sigma12;
27double sigma6;
28double r_cut2;
29
30constexpr int velocity_prop = 0;
31constexpr int force_prop = 1;
32
34DEFINE_INTERACTION_3D(ln_force)
35 Point<3,double> r = xp - xq;
36 double rn = norm2(r);
37
38 if (rn > r_cut2) return 0.0;
39
40 return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
41END_INTERACTION
42
43int main(int argc, char* argv[])
44{
46 openfpm_init(&argc,&argv);
47
49 double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
50 sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
51 double rc2 = r_cut * r_cut;
52
55 size_t sz[3] = {10,10,10};
56 Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
57 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
58 Ghost<3,double> ghost(r_cut);
59
61 ln_force lennard_jones;
62
67 Init_grid(sz,particles);
68
69 // Take aliases for the particle force, velocity and position
70 // we use them to write math expressions
71 auto force = getV<force_prop>(particles);
72 auto velocity = getV<velocity_prop>(particles);
73 auto position = getV<PROP_POS>(particles);
74
75 // set the velocity to zero (a math expression)
76 velocity = 0;
77
80 auto NN = particles.getCellListSym(r_cut);
81 force = applyKernel_in_sim(particles,NN,lennard_jones);
82
83 // Molecula Dynamic time stepping
84 for (size_t i = 0; i < 10000 ; i++)
85 {
86 // 1-step Verlet velocity
87 // v(t + 1/2*dt) = v(t) + 1/2*force(t)*dt
88 // x(t + dt) = x(t) + v(t + 1/2*dt)
89 velocity = velocity + 0.5*dt*force;
90 position = position + velocity*dt;
91
92 // redistribute particles and fill ghost
93 particles.map();
95
96 // Calculate the force at t + dt
98 force = applyKernel_in_sim(particles,NN,lennard_jones);
99
100 // 2-step Verlet velocity
101 // v(t+dt) = v(t + 1/2*dt) + 1/2*force(t+dt)*dt
102 velocity = velocity + 0.5*dt*force;
103 }
104
105 openfpm_finalize();
106}
107
This class represent an N-dimensional box.
Definition Box.hpp:61
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Distributed vector.
grid_dist_id_iterator_dec< Decomposition > getGridIterator(const size_t(&sz)[dim])
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor.
CellL getCellListSym(St r_cut)
Construct a cell list symmetric based on a cut of radius.
auto getLastPos() -> decltype(v_pos.template get< 0 >(0))
Get the position of the last element.
void add()
Add local particle.
void updateCellListSym(CellL &cell_list)
Update a cell list using the stored particles.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...