Skip to content

GitLab

Commit 509a43a2 authored by incardon's avatar incardon
Browse files
Options

Adding missing files

parent 88c20d93
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1531 additions and 716 deletions
example/Vector/3_molecular_dynamic_gpu_opt/Makefile
View file @ 509a43a2
include ../../example.mk
CC_SCOREP=scorep --nocompiler --cuda --mpp=mpi nvcc
CC=${CC_SCOREP}
CC=nvcc #${CC_SCOREP}
CC_MPI=mpic++
LDIR =
OBJ = main.o
OBJ_GPU = main_gpu.o
OBJ_CPU = main_cpu.o
OBJ_CPU_BEST = main_cpu_best.o
all: md_dyn
all: md_dyn_gpu md_dyn_cpu md_dyn_cpu_best
%.o: %.cu
$(CC) -O3 -g -c -isystem=/home/i-bird/MPI/include --std=c++11 -o $@ $< $(INCLUDE_PATH_NVCC)
md_dyn: $(OBJ)
%.o: %.cpp
$(CC_MPI) -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
md_dyn_gpu: $(OBJ_GPU)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH_NVCC) $(LIBS) -L/home/i-bird/MPI/lib -L/usr/local/cuda/lib64 -lcudart -lmpi -L/usr/local/cuda/extras/CUPTI/lib64 -lhdf5
md_dyn_cpu: $(OBJ_CPU)
$(CC_MPI) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
md_dyn_cpu_best: $(OBJ_CPU_BEST)
$(CC_MPI) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
run: all
mpirun -np 3 ./md_dyn && mpirun -np 3 ./md_dyn_expr && mpirun -np 3 ./md_dyn_vl;
mpirun -np 3 ./md_dyn_gpu && mpirun -np 3 ./md_dyn_cpu && mpirun -np 3 ./md_dyn_cpu_best;
.PHONY: clean all run
clean:
rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl
rm -f *.o *~ core md_dyn_gpu md_dyn_cpu md_dyn_cpu_best
example/Vector/3_molecular_dynamic_gpu_opt/config.cfg
View file @ 509a43a2
[pack]
files = main.cu Makefile
files = main.cu main_cpu.cpp Makefile
example/Vector/3_molecular_dynamic_gpu_opt/main.cu deleted 100644 → 0
View file @ 88c20d93
/*! \page Vector_3_md Vector 3 molecular dynamic
*
* \subpage Vector_3_md_dyn
* \subpage Vector_3_md_vl
*
*/
#include "Vector/vector_dist.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.hpp"
typedef float real_number;
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list on GPU
*
* [TOC]
*
* # Molecular Dynamic with Lennard-Jones potential # {#e3_md}
*
* This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime.
* Particle feel each other by the potential. In this example we are going to use GPU instead of CPU
*
* \f$ V(x_p,x_q) = 4( (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6 ) \f$
*
* ## Constants ##
*
* Here we define some useful constants
*
* \snippet Vector/3_molecular_dynamic/main.cpp constants
*
*/
//! \cond [constants] \endcond
constexpr int velocity = 0;
constexpr int force = 1;
constexpr int energy = 2;
//! \cond [constants] \endcond
/*!
*
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Calculate forces ## {#e3_md_cf}
*
* In this function we calculate the forces between particles. It require the vector of particles
* Cell list and scaling factor for the Lennard-Jhones potential.
*
* \snippet Vector/3_molecular_dynamic/main.cpp calc forces
* \snippet Vector/3_molecular_dynamic/main.cpp calc forces2
*
*
* In the following we are going into the detail of this function
*
*/
//! \cond [calc forces] \endcond
template<typename vector_dist_type,typename NN_type>
__global__ void calc_force_gpu(vector_dist_type vd, NN_type NN, real_number sigma12, real_number sigma6, real_number r_cut2)
{
auto p = GET_PARTICLE(vd);
// Get the position xp of the particle
Point<3,real_number> xp = vd.getPos(p);
// Reset the force counter
vd.template getProp<force>(p)[0] = 0.0;
vd.template getProp<force>(p)[1] = 0.0;
vd.template getProp<force>(p)[2] = 0.0;
// Get an iterator over the neighborhood particles of p
auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
// For each neighborhood particle ...
while (Np.isNext())
{
// ... q
auto q = Np.get_sort();
// if (p == q) skip this particle
if (q == p) {++Np; continue;};
// Get the position of p
Point<3,real_number> xq = vd.getPos(q);
// Get the distance between p and q
Point<3,real_number> r = xp - xq;
// take the norm of this vector
real_number rn = norm2(r);
if (rn > r_cut2)
{++Np; continue;};
// Calculate the force, using pow is slower
Point<3,real_number> 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);
vd.template getProp<force>(p)[1] += f.get(1);
vd.template getProp<force>(p)[2] += f.get(2);
// Next neighborhood
++Np;
}
}
template<typename vector_dist_type>
__global__ void update_velocity_position(vector_dist_type vd, real_number dt)
{
auto p = GET_PARTICLE(vd);
// here we calculate v(tn + 0.5)
vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
// here we calculate x(tn + 1)
vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
}
template<typename vector_dist_type>
__global__ void update_velocity(vector_dist_type vd, real_number dt)
{
auto p = GET_PARTICLE(vd);
// here we calculate v(tn + 1)
vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
}
template<typename vector_dist_type,typename NN_type>
__global__ void particle_energy(vector_dist_type vd, NN_type NN, real_number sigma12, real_number sigma6, real_number shift, real_number r_cut2)
{
auto p = GET_PARTICLE(vd);
// Get the position of the particle p
Point<3,real_number> xp = vd.getPos(p);
// Get an iterator over the neighborhood of the particle p
auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
real_number E = 0;
// For each neighborhood of the particle p
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
// if p == q skip this particle
if (q == p) {++Np; continue;};
// Get position of the particle q
Point<3,real_number> xq = vd.getPos(q);
// take the normalized direction
real_number rn = norm2(xp - xq);
if (rn > r_cut2)
{++Np;continue;}
// potential energy (using pow is slower)
E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
// Next neighborhood
++Np;
}
// Kinetic energy of the particle given by its actual speed
vd.template getProp<energy>(p) = 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;
}
template<typename CellList> void calc_forces(vector_dist_gpu<3,real_number, aggregate<real_number[3],real_number[3],real_number> > & vd, CellList & NN, real_number sigma12, real_number sigma6, real_number r_cut2)
{
//! \cond [calc forces] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* This function in called several time and require neighborhood of each particle. In order to speed-up the
* Cell-list construction we can use updateCellList function to reuse the memory of the previous cell-list.
* updateCellList can be faster than createCellList
*
* \see \ref e1_part_celllist
*
* \snippet Vector/3_molecular_dynamic/main.cpp ucl
*
*/
//! \cond [ucl] \endcond
vd.updateCellList(NN);
//! \cond [ucl] \endcond
/*!
*
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
* and calculate the force based on the Lennard-Jhones potential given by
*
* \f$ F(x_p,x_q) = 24( \frac{2 \sigma^{12}}{r^{13}} - \frac{\sigma^6}{r^{7}}) \hat{r} \f$
*
* \see \ref e0_s_assign_pos
*
* \snippet Vector/3_molecular_dynamic/main.cpp force calc
*
*/
//! \cond [force calc] \endcond
// Get an iterator over particles
auto it2 = vd.getDomainIteratorGPU();
calc_force_gpu<<<it2.wthr,it2.thr>>>(vd.toKernel_sorted(),NN.toKernel(),sigma12,sigma6,r_cut2);
//! \cond [force calc] \endcond
//! \cond [calc forces2] \endcond
}
//! \cond [calc forces2] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Calculate energy ## {#e3_md_ce}
*
* We also need a function to calculate energy, this function require the same parameter as calculate forces
*
* \snippet Vector/3_molecular_dynamic/main.cpp calc energy
* \snippet Vector/3_molecular_dynamic/main.cpp calc energy2
*
*/
//! \cond [calc energy] \endcond
template<typename CellList> real_number calc_energy(vector_dist_gpu<3,real_number, aggregate<real_number[3],real_number[3],real_number> > & vd, CellList & NN, real_number sigma12, real_number sigma6, real_number r_cut2)
{
real_number rc = r_cut2;
real_number shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
//! \cond [calc energy] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Reset the counter for the energy counter and
* update the cell list from the actual particle configuration
*
* \snippet Vector/3_molecular_dynamic/main.cpp up cell ene
*
* In the following we are going into the detail of this function
*
*/
//! \cond [up cell ene] \endcond
vd.updateCellList(NN);
//! \cond [up cell ene] \endcond
/*!
*
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* First we get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
* and calculate the energy based on the Lennard-Jhones potential given by
*
* \f$ V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) \f$
*
* \see \ref e0_s_assign_pos
*
* \snippet Vector/3_molecular_dynamic/main.cpp energy calc comp
*
*/
//! \cond [energy calc comp] \endcond
auto it2 = vd.getDomainIteratorGPU();
particle_energy<<<it2.wthr,it2.thr>>>(vd.toKernel(),NN.toKernel(),sigma12,sigma6,shift,r_cut2);
//! \cond [energy calc comp] \endcond
//! \cond [calc energy2] \endcond
// Calculated energy
return reduce<energy>(vd);
}
//! \cond [calc energy2] \endcond
int main(int argc, char* argv[])
{
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Initialization ## {#e3_md_init}
*
* After we defined the two main function calc forces and calc energy, we Initialize
* the library, we create a Box that define our domain, boundary conditions and ghost.
* Than we define important parameters of the simulation, time step integration,
* size of the box, and cut-off radius of the interaction. We also define 2 vectors
* x and y (they are like std::vector) used for statistic
*
* \see \ref e0_s_init
*
* \snippet Vector/3_molecular_dynamic/main.cpp constants run
*
*/
//! \cond [constants run] \endcond
openfpm_init(&argc,&argv);
real_number sigma = 0.01;
real_number r_cut =3.0*sigma;
// we will use it do place particles on a 10x10x10 Grid like
size_t sz[3] = {100,100,100};
// domain
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost, big enough to contain the interaction radius
Ghost<3,float> ghost(r_cut);
real_number dt = 0.00005;
real_number sigma12 = pow(sigma,12);
real_number sigma6 = pow(sigma,6);
openfpm::vector<real_number> x;
openfpm::vector<openfpm::vector<real_number>> y;
//! \cond [constants run] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* Than we define a distributed vector in 3D, containing 2 vectorial properties the
* first is the actual velocity of the particle the other is the force
*
* \see \ref e0_s_vector_inst
*
* \snippet Vector/3_molecular_dynamic/main.cpp vect create
*
*/
//! \cond [vect create] \endcond
vector_dist_gpu<3,real_number, aggregate<real_number[3],real_number[3],real_number> > vd(0,box,bc,ghost);
//! \cond [vect create] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Particles on a grid like position ## {#e3_md_gl}
*
* We define a grid iterator, to create particles on a grid like way. In the same cycle we also reset
* force and velocity
*
* \see \ref e1_cl_gr_it
*
* \snippet Vector/3_molecular_dynamic/main.cpp vect grid
*
*/
//! \cond [vect grid] \endcond
// We create the grid iterator
auto it = vd.getGridIterator(sz);
while (it.isNext())
{
// Create a new particle
vd.add();
// key contain (i,j,k) index of the grid
auto key = it.get();
// The index of the grid can be accessed with key.get(0) == i, key.get(1) == j ...
// We use getLastPos to set the position of the last particle added
vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
// We use getLastProp to set the property value of the last particle we added
vd.template getLastProp<velocity>()[0] = 0.0;
vd.template getLastProp<velocity>()[1] = 0.0;
vd.template getLastProp<velocity>()[2] = 0.0;
vd.template getLastProp<force>()[0] = 0.0;
vd.template getLastProp<force>()[1] = 0.0;
vd.template getLastProp<force>()[2] = 0.0;
++it;
}
vd.hostToDevicePos();
vd.hostToDeviceProp<velocity,force>();
vd.map(RUN_ON_DEVICE);
vd.ghost_get<>(RUN_ON_DEVICE);
//! \cond [vect grid] \endcond
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Molecular dynamic steps ## {#e3_md_vi}
*
* Here we do 10000 MD steps using verlet integrator
*
* The verlet integration stepping look like this
*
* \f[ \vec{v}(t_{n}+1/2) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f]
* \f[ \vec{x}(t_{n}+1) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f]
*
* calculate the forces from \f$ \vec{a} (t_{n}) \f$ finally
*
* \f[ \vec{v}(t_{n+1}) = \vec{v}_p(t_n+1/2) + \frac{1}{2} \delta t \vec{a}(t_n+1) \f]
*
* The cell-list structure is required to calculate forces. As demonstration
* purpose instead of using the standard Cell-list with (getCellList) we use the CELL_MEMBAL
* type. The impact in performance of using the CELL_MEMBAL instead of CELL_MEMFAST is less
* than 1% on the other hand CELL_MEMFAST use 16 Megabyte of memory
*
* \see \ref e1_cls_types
*
* Inside this cycle we are using several features that has been explained before in particular
*
* \see \ref e0_s_assign_pos
*
* \see \ref e0_s_map
*
* \see \ref e0_s_reduce
*
* \see \ref e1_part_ghost
*
* \see \ref e0_s_vis_vtk
*
* \snippet Vector/3_molecular_dynamic/main.cpp md steps
*
*/
timer tsim;
tsim.start();
//! \cond [md steps] \endcond
// Get the Cell list structure
auto NN = vd.getCellListGPU(r_cut);
// The standard
// auto NN = vd.getCellList(r_cut);
// calculate forces
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
unsigned long int f = 0;
// MD time stepping
for (size_t i = 0; i < 1000 ; i++)
{
// Get the iterator
auto it3 = vd.getDomainIteratorGPU();
update_velocity_position<<<it3.wthr,it3.thr>>>(vd.toKernel(),dt);
// Because we moved the particles in space we have to map them and re-sync the ghost
vd.map(RUN_ON_DEVICE);
vd.template ghost_get<>(RUN_ON_DEVICE);
// calculate forces or a(tn + 1) Step 2
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
std::cout << "STEP " << std::endl;
// Integrate the velocity Step 3
auto it4 = vd.getDomainIteratorGPU();
update_velocity<<<it4.wthr,it4.thr>>>(vd.toKernel(),dt);
// After every iteration collect some statistic about the configuration
if (i % 100 == 0)
{
vd.deviceToHostPos();
vd.deviceToHostProp<0,1,2>();
// We write the particle position for visualization (Without ghost)
vd.deleteGhost();
vd.write_frame("particles_",f);
// we resync the ghost
vd.ghost_get<>(RUN_ON_DEVICE);
// We calculate the energy
real_number energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
auto & vcl = create_vcluster();
vcl.sum(energy);
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});
// 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 << std::endl;
f++;
}
}
//! \cond [md steps] \endcond
tsim.stop();
std::cout << "Time: " << tsim.getwct() << std::endl;
/*!
* \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
*
* ## Plotting graphs ## {#e3_md_pg}
*
* After we terminate the MD steps our vector x contains at which iteration we calculated the energy
* while y contains the energy value at that time-step. We can produce a graph X Y
*
* \note The graph produced is an svg graph that can be view with a browser. From the browser we can
* also easily save the graph into pure svg format
*