main_data_ord.cpp

#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.hpp"
#include "energy_force.hpp"
 
 
int main(int argc, char* argv[])
{
 
    double dt = 0.0001;
    double r_cut = 0.03;
    double sigma = r_cut/3.0;
    double sigma12 = pow(sigma,12);
    double sigma6 = pow(sigma,6);
 
    openfpm::vector<double> x;
    openfpm::vector<openfpm::vector<double>> y;
 
    openfpm_init(&argc,&argv);
    Vcluster<> & v_cl = create_vcluster();
 
    // we will use it do place particles on a 10x10x10 Grid like
    size_t sz[3] = {40,40,40};
 
    // domain
    Box<3,double> 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,double> ghost(r_cut);
 
    // Create vector
    vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
 
 
 
    auto it = vd.getGridIterator(sz);
 
    while (it.isNext())
    {
        vd.add();
 
        auto key = it.get();
 
        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);
 
        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;
    }
 
 
 
    // Get the Cell list structure
    auto NN = vd.getCellList(r_cut);
 
    // calculate forces
    calc_forces(vd,NN,sigma12,sigma6);
    unsigned long int f = 0;
 
    timer time2;
    time2.start();
 
    #ifndef TEST_RUN
    size_t Nstep = 30000;
    #else
    size_t Nstep = 300;
    #endif
 
    // MD time stepping
    for (size_t i = 0; i < Nstep ; i++)
    {
        // Get the iterator
        auto it3 = vd.getDomainIterator();
 
        // integrate velicity and space based on the calculated forces (Step1)
        while (it3.isNext())
        {
            auto p = it3.get();
 
            // 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;
 
            ++it3;
        }
 
        // Because we mooved the particles in space we have to map them and re-sync the ghost
        vd.map();
        vd.template ghost_get<>();
 
 
        if (i % 200 == 0)
            vd.reorder(5);
 
 
 
        timer time;
        if (i % 10 == 0)
            time.start();
 
 
        // calculate forces or a(tn + 1) Step 2
        calc_forces(vd,NN,sigma12,sigma6);
 
 
        if (i % 10 == 0)
        {
            time.stop();
            x.add(i);
            y.add({time.getwct()});
        }
 
 
        // Integrate the velocity Step 3
        auto it4 = vd.getDomainIterator();
 
        while (it4.isNext())
        {
            auto p = it4.get();
 
            // 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];
 
            ++it4;
        }
 
        // After every iteration collect some statistic about the confoguration
        if (i % 100 == 0)
        {
            // We calculate the energy
            double energy = calc_energy(vd,NN,sigma12,sigma6);
            auto & vcl = create_vcluster();
            vcl.sum(energy);
            vcl.execute();
 
            // 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++;
        }
    }
 
    time2.stop();
    std::cout << "Performance: " << time2.getwct() << std::endl;
 
 
 
    // Google charts options, it store the options to draw the X Y graph
    GCoptions options;
 
    // Title of the graph
    options.title = std::string("Force calculation time");
 
    // Y axis name
    options.yAxis = std::string("Time");
 
    // X axis name
    options.xAxis = std::string("iteration");
 
    // width of the line
    options.lineWidth = 1.0;
 
    // Object that draw the X Y graph
    GoogleChart cg;
 
    // Add the graph
    // The graph that it produce is in svg format that can be opened on browser
    cg.AddLinesGraph(x,y,options);
 
    // Write into html format
    cg.write("gc_plot2_out.html");
 
 
 
    openfpm_finalize();
 
 
}