main.cpp

#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"


void init(grid_dist_id<2,double,aggregate<double> > & g_dist, const size_t (& sz)[2])
{

    // Get the iterator
    auto it = g_dist.getDomainGhostIterator();

    // For each point in the grid
    while (it.isNext())
    {


        auto key = it.get();



        auto key_g = g_dist.getGKey(key);



        if (key_g.get(0) == 0 || key_g.get(0) == sz[0] ||
            key_g.get(1) == 0 || key_g.get(1) == sz[1])
        {
            // Boundary part
            g_dist.template get<0>(key) = 0.0;
        }
        else
        {
            // Internal part
            g_dist.template get<0>(key) = 0.0;
        }



        ++it;

    }

}


constexpr int x = 0;
constexpr int y = 1;

int main(int argc, char* argv[])
{

    openfpm_init(&argc,&argv);
    


    Box<2,double> domain({-1.0,-1.0},{1.0,1.0});
    size_t sz[2] = {64,64};
    
    periodicity<2> bc = {NON_PERIODIC,NON_PERIODIC};
    
    // Ghost in grid unit
    Ghost<2,long int> g(1);
    


    grid_dist_id<2, double, aggregate<double>> g_dist(sz,domain,g,bc);
    
    // spacing between two points
    double spacing[2] = {g_dist.spacing(0),g_dist.spacing(1)};



    init(g_dist,sz);



    // sync the ghost property 0
    g_dist.template ghost_get<0>();



    // flag that indicate if we are processing red or black
    // we start from red
    bool red_black = true;

    // 10000 iterations
    for (size_t i = 0 ; i < 10000 ; i++)
    {

        auto dom = g_dist.getDomainIterator();
        
        // Iterate over all the points
        while (dom.isNext())
        {

            // Get the local grid key
            auto key = dom.get();

            // Here we convert the local grid position, into global position, key_g contain 3 integers that identify the position
            // of the grid point in global coordinates
            auto key_g = g_dist.getGKey(key);


            //
            // If we are processing red and is odd jump to the next point
            // If we are processing black and is even jump to the next point
            //
            if (red_black == false && (key_g.get(0) + key_g.get(1)) % 2 == 0)
            {++dom; continue;}
            else if (red_black == true && (key_g.get(0) + key_g.get(1)) % 2 == 1)
            {++dom; continue;}

            //
            // Update the grid values
            //
            // P.S. The keyword template is removed, it is possible only if we are in a function
            // without template parameters (if you are unsure use the keyword template)
            //
            g_dist.get<0>(key) = (g_dist.get<0>(key.move(x,1)) + g_dist.template get<0>(key.move(x,-1)) +
                                 g_dist.get<0>(key.move(y,1)) + g_dist.template get<0>(key.move(y,-1)) +
                                 - 1.0)/4.0;

            //
            // next point (red/black)
            ++dom;
        }



        g_dist.template ghost_get<0>();

        // switch from red to black or black to red
        red_black = !red_black;

    }
    


    // It contain the error
    double error = 0.0;

    // Get the iterator
    auto dom = g_dist.getDomainIterator();

    // Iterate over all the points
    while (dom.isNext())
    {
        // same the the grid point and the global grid point
        auto key = dom.get();

        // Calculate the error on each point
        // The error is how much the solution does not respect the equation
        double error_tmp = abs((g_dist.get<0>(key.move(x,1)) + g_dist.get<0>(key.move(x,-1)) +
                           g_dist.get<0>(key.move(y,1)) + g_dist.get<0>(key.move(y,-1)) +
                           - 4.0*g_dist.get<0>(key)) - 1.0);

        // In the norm infinity the maximum error across all the point is
        // important
        if (error_tmp > error)
            error = error_tmp;

        // next point
        ++dom;

    }

    // Get the maximum across processor to calculate the norm infinity of the error
    // Norm infinity of the error is the maximum error across all the grid points
    Vcluster & v_cl = create_vcluster();
    v_cl.max(error);
    v_cl.execute();

    // Only master print the norm infinity
    if (v_cl.getProcessUnitID() == 0)
            std::cout << "Error norm infinity: " << error << std::endl;



    g_dist.write("output");

    


    openfpm_finalize();


}