Simple example of grid usage to solve an equation
This example show the usage of grid to solve the following equation
\(\frac{\partial^2 u}{\partial^2 x} + \frac{\partial^2 u}{\partial^2 y} = 1\)
\(u(x,y) = 0 \)
at the boundary. This lead to the solution shown in the picture
Field initialization
In order to initialize the field U, first we get an iterator that cover domain + Ghost to iterate all the grid points. Inside the cycle we initialize domain and border (Ghost part to 0.0)
- See Also
- Loop over grid points
Get the local grid key
- See Also
- Grid coordinates
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
- See Also
- Grid coordinates
auto key_g = g_dist.getGKey(
key);
Initialize to 0, domain + boundary
if (key_g.get(0) == 0 || key_g.get(0) == sz[0] ||
key_g.get(1) == 0 || key_g.get(1) == sz[1])
{
g_dist.template get<0>(
key) = 0.0;
}
else
{
g_dist.template get<0>(
key) = 0.0;
}
The full function look like this
{
auto it = g_dist.getDomainGhostIterator();
while (it.isNext())
{
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])
{
g_dist.template get<0>(
key) = 0.0;
}
else
{
g_dist.template get<0>(
key) = 0.0;
}
++it;
}
}
Initialization
Initialize the library
- See Also
- Initialization
openfpm_init(&argc,&argv);
Grid instantiation and initialization
Create
- A 2D box that define the domain
- an array of 2 unsigned integer that will define the size of the grid on each dimension
- A Ghost object that will define the extension of the ghost part for each sub-domain in grid point unit
Create a distributed grid in 2D (1° template parameter) space in with double precision (2° template parameter) each grid point contain a scalar (double),
Constructor parameters:
- sz: size of the grid on each dimension
- domain: where the grid is defined
- g: ghost extension
- bc: boundary conditions
- See Also
- Grid instantiation
double spacing[2] = {g_dist.
spacing(0),g_dist.spacing(1)};
Initialize U and fill the boundary conditions
- See Also
- e2_se_field_init
Ghost synchronization
Before the computation read the ghost point we have to guarantee that they have updated values.
g_dist.template ghost_get<0>();
Red-Black alghorithm
Do 10000 iteration of Red-Black Gauss-Siedel iterations
Get an iterator that go through the points of the grid (No ghost) To compute one iteration.
- See Also
- Loop over grid points
-
Grid coordinates
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key_g = g_dist.getGKey(
key);
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;}
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;
++dom;
}
Once an iteration is done we have to synchronize the ghosts to start a new iteration. Consider that we calculated the red points in the next iteration the red points in the ghost are used in reading. This mean that the ghost must have the updated values
g_dist.template ghost_get<0>();
red_black = !red_black;
The full Algorithm look like this
bool red_black = true;
for (size_t i = 0 ; i < 10000 ; i++)
{
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key_g = g_dist.getGKey(
key);
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;}
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;
++dom;
}
g_dist.template ghost_get<0>();
red_black = !red_black;
}
Solution statistic
Once we got the solution we want to check if it really satisfy the equation, and calculate the error (norm infinity)
double error = 0.0;
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
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);
if (error_tmp > error)
error = error_tmp;
++dom;
}
std::cout << "Error norm infinity: " << error << std::endl;
VTK Write and visualization
Finally we want a nice output to visualize the information stored by the distributed grid. The function write by default produce VTK files. One for each processor that can be visualized with the programs like paraview
- See Also
- VTK and visualization
Finalize
At the very end of the program we have always to de-initialize the library
Full code
#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
{
auto it = g_dist.getDomainGhostIterator();
while (it.isNext())
{
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])
{
g_dist.template get<0>(
key) = 0.0;
}
else
{
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);
size_t sz[2] = {64,64};
double spacing[2] = {g_dist.spacing(0),g_dist.spacing(1)};
init(g_dist,sz);
g_dist.template ghost_get<0>();
bool red_black = true;
for (size_t i = 0 ; i < 10000 ; i++)
{
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key_g = g_dist.getGKey(
key);
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;}
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;
++dom;
}
g_dist.template ghost_get<0>();
red_black = !red_black;
}
double error = 0.0;
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
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);
if (error_tmp > error)
error = error_tmp;
++dom;
}
std::cout << "Error norm infinity: " << error << std::endl;
g_dist.write("output");
openfpm_finalize();
}