8 #include "Vector/vector_dist.hpp"
9 #include "Plot/GoogleChart.hpp"
10 #include "Plot/util.hpp"
35 constexpr
int velocity = 0;
36 constexpr
int force = 1;
59 template<
typename CellList>
void calc_forces(
vector_dist<3,
double,
aggregate<
double[3],
double[3]> > & vd,
CellList & NN,
double sigma12,
double sigma6,
double r_cut2)
79 vd.updateCellList(NN);
101 auto it2 = vd.getDomainIterator();
113 vd.template getProp<force>(p)[0] = 0.0;
114 vd.template getProp<force>(p)[1] = 0.0;
115 vd.template getProp<force>(p)[2] = 0.0;
118 auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
127 if (q == p.getKey()) {++Np;
continue;};
136 double rn = norm2(r);
142 Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
145 vd.template getProp<force>(p)[0] += f.
get(0);
146 vd.template getProp<force>(p)[1] += f.
get(1);
147 vd.template getProp<force>(p)[2] += f.
get(2);
179 template<
typename CellList>
double calc_energy(
vector_dist<3,
double,
aggregate<
double[3],
double[3]> > & vd,
CellList & NN,
double sigma12,
double sigma6,
double r_cut2)
182 double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
223 auto it2 = vd.getDomainIterator();
235 auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
244 if (q == p.getKey()) {++Np;
continue;};
250 double rn = norm2(xp - xq);
256 E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
263 E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
264 vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
265 vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
281 int main(
int argc,
char* argv[])
302 openfpm_init(&argc,&argv);
305 double r_cut = 3.0*sigma;
308 size_t sz[3] = {10,10,10};
314 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
320 double sigma12 = pow(sigma,12);
321 double sigma6 = pow(sigma,6);
363 auto it = vd.getGridIterator(sz);
375 vd.getLastPos()[0] =
key.get(0) * it.getSpacing(0);
376 vd.getLastPos()[1] =
key.get(1) * it.getSpacing(1);
377 vd.getLastPos()[2] =
key.get(2) * it.getSpacing(2);
380 vd.template getLastProp<velocity>()[0] = 0.0;
381 vd.template getLastProp<velocity>()[1] = 0.0;
382 vd.template getLastProp<velocity>()[2] = 0.0;
384 vd.template getLastProp<force>()[0] = 0.0;
385 vd.template getLastProp<force>()[1] = 0.0;
386 vd.template getLastProp<force>()[2] = 0.0;
438 auto NN = vd.getCellList<CELL_MEMBAL(3,
double)>(r_cut);
444 calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
445 unsigned long int f = 0;
448 for (
size_t i = 0; i < 10000 ; i++)
451 auto it3 = vd.getDomainIterator();
459 vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
460 vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
461 vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
464 vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
465 vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
466 vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
473 vd.template ghost_get<>();
476 calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
480 auto it4 = vd.getDomainIterator();
487 vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
488 vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
489 vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
499 vd.write(
"particles_",f);
505 double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
506 auto & vcl = create_vcluster();
516 if (vcl.getProcessUnitID() == 0)
517 std::cout <<
"Energy: " << energy << std::endl;
526 std::cout <<
"Time: " << tsim.
getwct() << std::endl;
549 options.
title = std::string(
"Energy with time");
552 options.
yAxis = std::string(
"Energy");
555 options.
xAxis = std::string(
"iteration");
561 options.
width = 1280;
567 options.
more = GC_ZOOM;
577 cg.
write(
"gc_plot2_out.html");
size_t heigh
height of the graph in pixels
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.
This class implement the point shape in an N-dimensional space.
double getwct()
Return the elapsed real time.
std::string yAxis
Y axis name.
Small class to produce graph with Google chart in HTML.
void start()
Start the timer.
const T & get(size_t i) const
Get coordinate.
size_t lineWidth
Width of the line.
This class is a trick to indicate the compiler a specific specialization pattern. ...
void write(std::string file)
It write the graphs on file in html format using Google charts.
size_t width
width of the graph in pixels
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Implementation of 1-D std::vector like structure.
Class for FAST cell list implementation.
std::string xAxis
X axis name.
Class for cpu time benchmarking.
void stop()
Stop the timer.