main.cu

#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
#include "timer.hpp"
 
#ifdef __NVCC__
 
constexpr int U = 0;
constexpr int V = 1;
constexpr int U_next = 2;
constexpr int V_next = 3;
 
typedef sgrid_dist_id_gpu<3,double,aggregate<double,double,double,double> > sgrid_type;
 
void init(sgrid_type & grid, Box<3,double> & domain)
{
    auto it = grid.getGridIterator();
    Point<3,double> p1;
    Point<3,double> p2;
    Point<3,double> p3;
 
    // Shere1
    for (int i = 0 ; i < 3 ; i++)
    {p1.get(i) = 2.0;}
 
        // Shere2
        for (int i = 0 ; i < 3 ; i++)
        {p2.get(i) = 1.0;}
 
        // Shere3
        for (int i = 0 ; i < 3 ; i++)
        {p3.get(i) = 0.5;}
 
    Sphere<3,double> sph1(p1,0.3);
    Sphere<3,double> sph2(p2,0.3);
    Sphere<3,double> sph3(p3,0.3);
 
    Point<3,double> u({1.0,1.0,1.0});
    Box<3,double> channel_box(p3,p1);
 
    double spacing_x = grid.spacing(0);
    double spacing_y = grid.spacing(1);
    double spacing_z = grid.spacing(2);
 
    typedef typename GetAddBlockType<sgrid_type>::type InsertBlockT;
 
    grid.addPoints([spacing_x,spacing_y,spacing_z,u,sph1,sph2,sph3,channel_box] __device__ (int i, int j, int k)
                    {
                        Point<3,double> pc({i*spacing_x,j*spacing_y,k*spacing_z});
                        Point<3,double> vp;
 
                        // calculate the distance from the diagonal
                        vp.get(0) = pc.get(1)*u.get(2) - pc.get(2)*u.get(1);
                        vp.get(1) = pc.get(2)*u.get(0) - pc.get(0)*u.get(2);
                        vp.get(2) = pc.get(0)*u.get(1) - pc.get(1)*u.get(0);
 
                        double distance = vp.norm() / sqrt(3.0f);
 
                        // Check if the point is in the domain
                        if (sph1.isInside(pc) || sph2.isInside(pc) || sph3.isInside(pc) || (distance < 0.1 && channel_box.isInside(pc)) )
                        {return true;}
 
                        return false;
                    },
                    [] __device__ (InsertBlockT & data, int i, int j, int k)
                    {
                        data.template get<U>() = 1.0;
                        data.template get<V>() = 0.0;
                    }
                    );
 
 
    grid.template flush<smax_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
 
    long int x_start = grid.size(0)*1.95f/domain.getHigh(0);
    long int y_start = grid.size(1)*1.95f/domain.getHigh(1);
    long int z_start = grid.size(1)*1.95f/domain.getHigh(2);
 
    long int x_stop = grid.size(0)*2.05f/domain.getHigh(0);
    long int y_stop = grid.size(1)*2.05f/domain.getHigh(1);
    long int z_stop = grid.size(1)*2.05f/domain.getHigh(2);
 
    grid_key_dx<3> start({x_start,y_start,z_start});
    grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
 
        grid.addPoints(start,stop,[] __device__ (int i, int j, int k)
                                {
                                                return true;
                                },
                                [] __device__ (InsertBlockT & data, int i, int j, int k)
                                {
                                        data.template get<U>() = 0.5;
                                        data.template get<V>() = 0.24;
                                }
                                );
 
    grid.template flush<smin_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
}
 
 
int main(int argc, char* argv[])
{
    openfpm_init(&argc,&argv);
 
    // domain
    Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
    
    // grid size
        size_t sz[3] = {256,256,256};
 
    // Define periodicity of the grid
    periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
    
    // Ghost in grid unit
    Ghost<3,long int> g(1);
    
    // deltaT
    double deltaT = 0.1;
 
    // Diffusion constant for specie U
    double du = 2*1e-5;
 
    // Diffusion constant for specie V
    double dv = 1*1e-5;
 
#ifdef TEST_RUN
        // Number of timesteps
        size_t timeSteps = 300;
#else
    // Number of timesteps
        size_t timeSteps = 150000;
#endif
 
    // K and F (Physical constant in the equation)
        double K = 0.053;
        double F = 0.014;
 
    sgrid_type grid(sz,domain,g,bc);
 
    grid.template setBackgroundValue<0>(-0.5);
    grid.template setBackgroundValue<1>(-0.5);
    grid.template setBackgroundValue<2>(-0.5);
    grid.template setBackgroundValue<3>(-0.5);
    
    // spacing of the grid on x and y
    double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
 
    init(grid,domain);
 
    // sync the ghost
    grid.template ghost_get<U,V>(RUN_ON_DEVICE);
 
    // because we assume that spacing[x] == spacing[y] we use formula 2
    // and we calculate the prefactor of Eq 2
    double uFactor = deltaT * du/(spacing[0]*spacing[0]);
    double vFactor = deltaT * dv/(spacing[0]*spacing[0]);
 
    grid.template deviceToHost<U,V>();
    grid.write("Init_condition");
 
    timer tot_sim;
    tot_sim.start();
 
    for (size_t i = 0; i < timeSteps; ++i)
    {
 
                typedef typename GetCpBlockType<decltype(grid),0,1>::type CpBlockType;
 
                auto func = [uFactor,vFactor,deltaT,F,K] __device__ (double & u_out, double & v_out,
                                                           CpBlockType & u, CpBlockType & v,
                                                           int i, int j, int k){
 
                        double uc = u(i,j,k);
                        double vc = v(i,j,k);
 
                        double u_px = u(i+1,j,k);
                        double u_mx = u(i-1,j,k);
 
                        double u_py = u(i,j+1,k);
                        double u_my = u(i,j-1,k);
 
                        double u_pz = u(i,j,k+1);
                        double u_mz = u(i,j,k-1);
 
                        double v_px = v(i+1,j,k);
                        double v_mx = v(i-1,j,k);
 
                        double v_py = v(i,j+1,k);
                        double v_my = v(i,j-1,k);
 
                        double v_pz = v(i,j,k+1);
                        double v_mz = v(i,j,k-1);
 
                        // U fix
 
                        if (u_mx < -0.1 && u_px < -0.1)
                        {
                            u_mx = uc;
                            u_px = uc;
                        }
 
                        if (u_mx < -0.1)
                        {u_mx = u_px;}
 
                        if (u_px < -0.1)
                        {u_px = u_mx;}
 
                        if (u_my < -0.1 && u_py < -0.1)
                        {
                            u_my = uc;
                            u_py = uc;
                        }
 
                        if (u_my < -0.1)
                        {u_my = u_py;}
 
                        if (u_py < -0.1)
                        {u_py = u_my;}
 
                        if (u_mz < -0.1 && u_pz < -0.1)
                        {
                            u_mz = uc;
                            u_pz = uc;
                        }
 
                        if (u_mz < -0.1)
                        {u_mz = u_pz;}
 
                        if (u_pz < -0.1)
                        {u_pz = u_mz;}
 
                        // V fix
 
                        if (v_mx < -0.1 && v_px < -0.1)
                        {
                            v_mx = uc;
                            v_px = uc;
                        }
 
                        if (v_mx < -0.1)
                        {v_mx = v_px;}
 
                        if (v_px < -0.1)
                        {v_px = v_mx;}
 
                        if (v_my < -0.1 && v_py < -0.1)
                        {
                            v_my = uc;
                            v_py = uc;
                        }
 
                        if (v_my < -0.1)
                        {v_my = v_py;}
 
                        if (v_py < -0.1)
                        {v_py = v_my;}
 
                        if (v_mz < -0.1 && v_pz < -0.1)
                        {
                            v_mz = uc;
                            v_pz = uc;
                        }
 
                        if (v_mz < -0.1)
                        {v_mz = v_pz;}
 
                        if (v_pz < -0.1)
                        {v_pz = v_mz;}
 
                        u_out = uc + uFactor *(u_mx + u_px +
                                                               u_my + u_py +
                                                               u_mz + u_pz - 6.0*uc) - deltaT * uc*vc*vc
                                                               - deltaT * F * (uc - 1.0);
 
 
                        v_out = vc + vFactor *(v_mx + v_px +
                                                               v_py + v_my +
                                                               v_mz + v_pz - 6.0*vc) + deltaT * uc*vc*vc
                                           - deltaT * (F+K) * vc;
 
                        };
 
 
                if (i % 2 == 0)
                {
                    grid.conv2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
 
                    // After copy we synchronize again the ghost part U and V
 
                    grid.ghost_get<U_next,V_next>(RUN_ON_DEVICE | SKIP_LABELLING);
                }
                else
                {
                    grid.conv2<U_next,V_next,U,V,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
 
                    // After copy we synchronize again the ghost part U and V
                    grid.ghost_get<U,V>(RUN_ON_DEVICE | SKIP_LABELLING);
                }
 
 
        // After copy we synchronize again the ghost part U and V
 
        // Every 500 time step we output the configuration for
        // visualization
/*      if (i % 500 == 0)
        {
            grid.save("output_" + std::to_string(count));
            count++;
        }*/
 
                std::cout << "STEP: " << i  << std::endl;
/*                if (i % 300 == 0)
                {
                    grid.template deviceToHost<U,V>();
                        grid.write_frame("out",i);
                }*/
    }
    
    tot_sim.stop();
    std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;
 
    grid.template deviceToHost<U,V>();
    grid.write("Final");
 
 
 
    openfpm_finalize();
 
 
}
 
#else
 
int main(int argc, char* argv[])
{
        return 0;
}
 
#endif