main_expr.cpp

#include "Vector/vector_dist.hpp"
#include "Plot/GoogleChart.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "timer.hpp"

constexpr int velocity = 0;
constexpr int force = 1;

struct ln_potential
{
    double sigma12,sigma6,r_cut2,shift;

    ln_potential(double sigma12_, double sigma6_, double r_cut2_, double shift_) {sigma12 = sigma12_; sigma6 = sigma6_; r_cut2 = r_cut2_;shift = shift_;}

    Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
    {
        double rn = norm2(xp - xq);
        if (rn >= r_cut2)   return 0.0;

        Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift,0.0});

        return E;
    }
};

struct ln_force
{
    double sigma12,sigma6,r_cut2;

    ln_force(double sigma12_, double sigma6_, double r_cut2_) {sigma12 = sigma12_; sigma6 = sigma6_;r_cut2 = r_cut2_;}

    Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
    {
        Point<3,double> r = xp - xq;
        double rn = norm2(r);

        if (rn > r_cut2)    return 0.0;

        return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) -  sigma6 / (rn*rn*rn*rn)) * r;
    }
};

int main(int argc, char* argv[])
{
    double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;

    double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
    double rc2 = r_cut * r_cut;
    double shift = 2.0 * ( sigma12 / (rc2*rc2*rc2*rc2*rc2*rc2) - sigma6 / ( rc2*rc2*rc2) );

    openfpm::vector<double> x;
    openfpm::vector<openfpm::vector<double>> y;


    openfpm_init(&argc,&argv);
    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,r_cut*r_cut);
    ln_potential lp(sigma12,sigma6,r_cut*r_cut,shift);

    vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);

    auto v_force = getV<force>(vd);
    auto v_velocity = getV<velocity>(vd);
    auto v_pos = getV<PROP_POS>(vd);

    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);

        ++it;
    }

    v_force = 0;
    v_velocity = 0;

    timer tsim;
    tsim.start();

    auto NN = vd.getCellList(r_cut);

    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++)
    {
        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
        vd.updateCellList(NN);
        v_force = applyKernel_in_sim(vd,NN,lf);

        v_velocity = v_velocity + 0.5*dt*v_force;

        if (i % 100 == 0)
        {
            vd.deleteGhost();
            vd.write("particles_",f);
            vd.ghost_get<>();

            vd.updateCellList(NN);
            Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).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({E.get(0),E.get(1),E.get(0) - E.get(1)});

            if (vcl.getProcessUnitID() == 0)
                std::cout << "Energy Total: " << E.get(0) << "   Kinetic: " <<  E.get(1) << "   Potential: " << E.get(0) - E.get(1) << std::endl;

            f++;
        }
    }

    tsim.stop();
    std::cout << "Time: " << tsim.getwct() << std::endl;

    GCoptions options;
    options.title = std::string("Energy with time");
    options.yAxis = std::string("Energy");
    options.xAxis = std::string("iteration");
    options.lineWidth = 1.0;

    GoogleChart cg;
    cg.AddLinesGraph(x,y,options);
    cg.write("gc_plot2_out.html");

    openfpm_finalize();
}