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