Commit 8a479ba4 authored by Pietro Incardona's avatar Pietro Incardona

Example latest

parent b116804e
......@@ -6,6 +6,7 @@
* \subpage grid_0_simple
* \subpage Grid_1_stencil
* \subpage Grid_2_solve_eq
* \subpage Grid_3_gs
*
*/
......
......@@ -225,6 +225,8 @@ int main(int argc, char* argv[])
*
* Finally we want a nice output to visualize the information stored by the distributed grid
*
* \see \ref e0_s_VTK_vis
*
* \snippet Grid/1_stencil/main.cpp output
*
*/
......
This diff is collapsed.
This diff is collapsed.
......@@ -48,7 +48,7 @@ constexpr int force = 1;
//! \cond [calc forces] \endcond
void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
{
//! \cond [calc forces] \endcond
......@@ -125,12 +125,10 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
Point<3,double> r = xp - xq;
// take the norm of this vector
float rn = r.norm();
r /= rn;
rn *= L;
double rn = norm2(r);
// Calculate the force, using pow is slower
Point<3,double> f = 24.0*(2.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) - 1.0 / (rn*rn*rn*rn*rn*rn*rn)) * r;
Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
// we sum the force produced by q on p
vd.template getProp<force>(p)[0] += f.get(0);
......@@ -167,7 +165,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
//! \cond [calc energy] \endcond
double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
{
//! \cond [calc energy] \endcond
......@@ -240,28 +238,21 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
// Get position of the particle q
Point<3,double> xq = vd.getPos(q);
// get the distance vector r between p and q
Point<3,double> r = xp - xq;
// take the normalized direction
float rn = r.norm();
r /= rn;
// Multiply by L
rn *= L;
double rn = norm2(xp - xq);
// potential energy (using pow is slower)
E += 4.0 * ( 1.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) + 1.0 / ( rn*rn*rn*rn*rn*rn) );
// Kinetic energy of the particle given by its actual speed
E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
// Next neighborhood
++Np;
}
// Kinetic energy of the particle given by its actual speed
E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
// Next Particle
++it2;
}
......@@ -294,9 +285,11 @@ int main(int argc, char* argv[])
//! \cond [constants] \endcond
double dt = 0.005;
double L = 10.0;
double dt = 0.0005;
float r_cut = 0.3;
double sigma = 0.1;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
openfpm::vector<double> x;
openfpm::vector<openfpm::vector<double>> y;
......@@ -432,7 +425,7 @@ int main(int argc, char* argv[])
auto NN = vd.getCellList(r_cut);
// calculate forces
calc_forces(vd,NN,L);
calc_forces(vd,NN,sigma12,sigma6);
unsigned long int f = 0;
// MD time stepping
......@@ -464,7 +457,7 @@ int main(int argc, char* argv[])
vd.template ghost_get<>();
// calculate forces or a(tn + 1) Step 2
calc_forces(vd,NN,L);
calc_forces(vd,NN,sigma12,sigma6);
// Integrate the velocity Step 3
......@@ -489,15 +482,15 @@ int main(int argc, char* argv[])
vd.deleteGhost();
vd.write("particles_",f);
// we resync the ghost
vd.ghost_get<>();
// We calculate the energy
double energy = calc_energy(vd,NN,L);
double energy = calc_energy(vd,NN,sigma12,sigma6);
auto & vcl = create_vcluster();
vcl.sum(energy);
vcl.execute();
// we resync the ghost
vd.ghost_get<>();
// we save the energy calculated at time step i c contain the time-step y contain the energy
x.add(i);
y.add({energy});
......@@ -582,6 +575,19 @@ int main(int argc, char* argv[])
* \include Vector/3_molecular_dynamic/main.cpp
*
*/
/*!
* \page Vector_3_md Vector 3 molecular dynamic
*
* # Code with expression # {#code}
*
* Here we also show how we can simplify the example using expressions
*
* \see \ref Vector_2_expression
*
* \include Vector/3_molecular_dynamic/main_expr.cpp
*
*/
}
......
#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
constexpr int velocity = 0;
......@@ -10,21 +7,16 @@ constexpr int force = 1;
struct ln_potential
{
double L;
double sigma12,sigma6;
ln_potential(double L)
:L(L){}
ln_potential(double sigma12_, double sigma6_) {sigma12 = sigma12_, sigma6 = sigma6_;}
Point<3,double> value(Point<3,double> & xp, Point<3,double> xq, Point<3,double> & sp, Point<3,double> & sq)
Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
{
// Energy is a vector we calculate separately Potential and kinetic, and total
Point<3,double> E;
double rn = norm2(xp - xq);
float rn = norm(xp - xq) * L;
E.get(0)= 4.0 * ( 1.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) + 1.0 / ( rn*rn*rn*rn*rn*rn) );
E.get(1) = (sp * sp) / 2.0;
E.get(2) = E.get(0) + E.get(1);
Point<2,double> E({4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
0.0});
return E;
}
......@@ -32,76 +24,40 @@ struct ln_potential
struct ln_force
{
double L;
double sigma12,sigma6;
ln_force(double L)
:L(L){}
ln_force(double sigma12_, double sigma6_) {sigma12 = sigma12_; sigma6 = sigma6_;}
Point<3,double> value(Point<3,double> & xp, Point<3,double> xq)
Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
{
Point<3,double> r = xp - xq;
double rn = norm2(r);
// take the norm of this vector
float rn = r.norm();
r /= rn;
rn *= L;
return 24.0*(2.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) - 1.0 / (rn*rn*rn*rn*rn*rn*rn)) * r;
return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
}
};
void calc_forces(vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
{
vd.updateCellList(NN);
//! \cond [ucl] \endcond
auto v_force = getV<force>(vd);
ln_force lf(L);
v_force = applyKernel_in_sim(v_force,vd,NN,lf);
}
Point<3,double> calc_energy(vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
{
Point<3,double> E;
vd.updateCellList(NN);
auto v_velocity = getV<velocity>(vd);
ln_potential lf(L);
auto eE = applyKernel_reduce(v_velocity,vd,NN,lf);
eE.init();
E = eE.value(0);
return E;
}
int main(int argc, char* argv[])
{
double dt = 0.005;
double L = 10.0;
float r_cut = 0.3;
double dt = 0.0005, sigma = 0.1, r_cut = 0.3;
double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
openfpm::vector<double> x;
openfpm::vector<openfpm::vector<double>> y;
openfpm_init(&argc,&argv);
Vcluster & v_cl = create_vcluster();
Vcluster & vcl = create_vcluster();
size_t sz[3] = {10,10,10};
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
Ghost<3,float> ghost(r_cut);
ln_force lf(sigma12,sigma6);
ln_potential lp(sigma12,sigma6);
vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);
auto v_force = getV<force>(vd);
......@@ -128,24 +84,22 @@ int main(int argc, char* argv[])
auto NN = vd.getCellList(r_cut);
// calculate forces
calc_forces(vd,NN,L);
vd.updateCellList(NN);
v_force = applyKernel_in_sim(vd,NN,lf);
unsigned long int f = 0;
// MD time stepping
for (size_t i = 0; i < 10000 ; i++)
{
auto exp1 = v_velocity + 0.5*dt*v_force;
auto exp2 = v_pos + v_velocity*dt;
assign(v_velocity,exp1,
v_pos,exp2);
assign(v_velocity, v_velocity + 0.5*dt*v_force,
v_pos, v_pos + v_velocity*dt);
vd.map();
vd.template ghost_get<>();
// calculate forces or a(tn + 1) Step 2
calc_forces(vd,NN,L);
vd.updateCellList(NN);
v_force = applyKernel_in_sim(vd,NN,lf);
v_velocity = v_velocity + 0.5*dt*v_force;
......@@ -153,24 +107,19 @@ int main(int argc, char* argv[])
{
vd.deleteGhost();
vd.write("particles_",f);
vd.ghost_get<>();
// We calculate the energy
Point<3,double> energy = calc_energy(vd,NN,L);
auto & vcl = create_vcluster();
vcl.sum(energy.get(2));
vcl.execute();
Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
// we resync the ghost
vd.ghost_get<>();
vcl.sum(E.get(0));vcl.sum(E.get(1));
vcl.execute();
// we save the energy calculated at time step i c contain the time-step y contain the energy
x.add(i);
y.add({energy.get(2)});
y.add({E.get(0),E.get(1),E.get(0) - E.get(1)});
// We also print on terminal the value of the energy
// only one processor (master) write on terminal
if (vcl.getProcessUnitID() == 0)
std::cout << "Energy: " << energy.get(2) << std::endl;
std::cout << "Energy Total: " << E.get(0) << " Kinetic: " << E.get(1) << " Potential: " << E.get(0) - E.get(1) << std::endl;
f++;
}
......
......@@ -1627,13 +1627,27 @@ public:
*/
inline bool write(std::string out, size_t iteration, int opt = NO_GHOST)
{
// CSVWriter test
CSVWriter<openfpm::vector<Point<dim, St>>, openfpm::vector<prop> > csv_writer;
if ((opt & 0xFFFF0000) == CSV_WRITER)
{
// CSVWriter test
CSVWriter<openfpm::vector<Point<dim, St>>, openfpm::vector<prop> > csv_writer;
std::string output = std::to_string(out + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".csv"));
std::string output = std::to_string(out + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".csv"));
// Write the CSV
return csv_writer.write(output, v_pos, v_prp);
// Write the CSV
return csv_writer.write(output, v_pos, v_prp);
}
else
{
// VTKWriter for a set of points
VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
vtk_writer.add(v_pos,v_prp,g_m);
std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".vtk"));
// Write the VTK file
return vtk_writer.write(output);
}
}
/* \brief It return the id of structure in the allocation list
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment