Stencil example

Stencil example and ghost

This example show how to move grid_key in order to create a Laplacian stencil, be careful, the function move is convenient, but not the fastest implementation. We also show how to do ghost communications

Define some convenient constants and types

 
constexpr size_t x = 0;
constexpr size_t y = 1;
constexpr size_t z = 2;
 
constexpr size_t A = 0;
constexpr size_t B = 1;
 

Initialization

Initialize the library and several objects

See also

Initialization

 
    openfpm_init(&argc,&argv);
 
    // domain
    Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // grid sizes
    size_t sz[3] = {100,100,100};
 
    // ghost extension
    Ghost<3,float> g(0.03);
 

Grid create

Create a distributed grid in 3D. With typedef we create an alias name for aggregate<float[3],float[3]>. In practice the type of grid_point == aggregate<float[3],float[3]>

See also

Grid instantiation

 
    // a convenient alias for aggregate<...>
    typedef aggregate<float,float> grid_point;
 
    grid_dist_id<3, float, grid_point> g_dist(sz,domain,g);
 

Loop over grid points

Get an iterator that go through the point of the domain (No ghost)

See also

Loop over grid points

 
    auto dom = g_dist.getDomainIterator();
 
    while (dom.isNext())
    {
 

 
        ++dom;
    }
 

Inside the cycle we get the local grid key

See also

Grid coordinates

 
        auto key = dom.get();
 

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

See also

Grid coordinates

 
        auto key_g = g_dist.getGKey(key);
 

we write on the grid point of position (i,j,k) the value i*i + j*j + k*k on the property A. Mathematically is equivalent to the function

\( f(x,y,z) = x^2 + y^2 + z^2 \)

 
        g_dist.template get<A>(key) = key_g.get(0)*key_g.get(0) + key_g.get(1)*key_g.get(1) + key_g.get(2)*key_g.get(2);
 

Ghost

Each sub-domain has an extended part, that is materially contained into another processor. In general is not synchronized ghost_get synchronize the property A in the ghost part

 
    g_dist.template ghost_get<A>();
    

Get again another iterator, iterate across all the domain points, calculating a Laplace stencil. Write the result on B

 
    auto dom2 = g_dist.getDomainIterator();
    
    while (dom2.isNext())
    {
        auto key = dom2.get();
 
        // Laplace stencil
        g_dist.template get<B>(key) = g_dist.template get<A>(key.move(x,1)) + g_dist.template get<A>(key.move(x,-1)) +
                                         g_dist.template get<A>(key.move(y,1)) + g_dist.template get<A>(key.move(y,-1)) +
                                         g_dist.template get<A>(key.move(z,1)) + g_dist.template get<A>(key.move(z,-1)) -
                                         6*g_dist.template get<A>(key);
 
        ++dom2;
    }
 

Finally we want a nice output to visualize the information stored by the distributed grid

See also

VTK and visualization

 
    g_dist.write("output");
 

Deinitialize the library

 
    openfpm_finalize();
 

Full code

#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
#include "Decomposition/CartDecomposition.hpp"
 
 
constexpr size_t x = 0;
constexpr size_t y = 1;
constexpr size_t z = 2;
 
constexpr size_t A = 0;
constexpr size_t B = 1;
 
 
int main(int argc, char* argv[])
{
 
    openfpm_init(&argc,&argv);
 
    // domain
    Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // grid sizes
    size_t sz[3] = {100,100,100};
 
    // ghost extension
    Ghost<3,float> g(0.03);
 
 
 
    // a convenient alias for aggregate<...>
    typedef aggregate<float,float> grid_point;
 
    grid_dist_id<3, float, grid_point> g_dist(sz,domain,g);
 
 
 
    auto dom = g_dist.getDomainIterator();
 
    while (dom.isNext())
    {
 
 
 
        auto key = dom.get();
 
 
 
        auto key_g = g_dist.getGKey(key);
 
 
 
        g_dist.template get<A>(key) = key_g.get(0)*key_g.get(0) + key_g.get(1)*key_g.get(1) + key_g.get(2)*key_g.get(2);
 
 
 
        ++dom;
    }
 
 
 
    g_dist.template ghost_get<A>();
    
 
 
    auto dom2 = g_dist.getDomainIterator();
    
    while (dom2.isNext())
    {
        auto key = dom2.get();
 
        // Laplace stencil
        g_dist.template get<B>(key) = g_dist.template get<A>(key.move(x,1)) + g_dist.template get<A>(key.move(x,-1)) +
                                         g_dist.template get<A>(key.move(y,1)) + g_dist.template get<A>(key.move(y,-1)) +
                                         g_dist.template get<A>(key.move(z,1)) + g_dist.template get<A>(key.move(z,-1)) -
                                         6*g_dist.template get<A>(key);
 
        ++dom2;
    }
 
 
 
    g_dist.write("output");
 
 
 
    openfpm_finalize();
 
 
}