# Domain, Halo and Padding regions

In this tutorial we will learn about data regions and how these impact the `Operator` construction. We will use a simple time marching example.

``````from devito import Eq, Grid, TimeFunction, Operator

grid = Grid(shape=(3, 3))
u = TimeFunction(name='u', grid=grid)
u.data[:] = 1``````

At this point, we have a time-varying 3x3 grid filled with 1’s. Below, we can see the domain data values:

``print(u.data)``
``````[[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]

[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]]``````

We now create a time-marching `Operator` that, at each timestep, increments by `2` all points in the computational domain.

``````from devito import configuration
configuration['language'] = 'C'

eq = Eq(u.forward, u+2)
op = Operator(eq, opt='noop')``````

We can print `op` to get the generated code.

``print(op)``
``````#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

struct dataobj
{
void *restrict data;
unsigned long * size;
unsigned long * npsize;
unsigned long * dsize;
int * hsize;
int * hofs;
int * oofs;
void * dmap;
} ;

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)
{
float (*restrict u)[u_vec->size[1]][u_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]][u_vec->size[2]]) u_vec->data;

for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 2;
}
}
STOP(section0,timers)
}

return 0;
}
``````

When we take a look at the constructed expression, `u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 2`, we see several `+1` were added to the `u`’s spatial indices.

This is because the domain region is actually surrounded by ‘ghost’ points, which can be accessed via a stencil when iterating in proximity of the domain boundary. The ghost points define the halo region. The halo region can be accessed through the `data_with_halo` data accessor. As we see below, the halo points correspond to the zeros surrounding the domain region.

``print(u.data_with_halo)``
``````[[[0. 0. 0. 0. 0.]
[0. 1. 1. 1. 0.]
[0. 1. 1. 1. 0.]
[0. 1. 1. 1. 0.]
[0. 0. 0. 0. 0.]]

[[0. 0. 0. 0. 0.]
[0. 1. 1. 1. 0.]
[0. 1. 1. 1. 0.]
[0. 1. 1. 1. 0.]
[0. 0. 0. 0. 0.]]]``````

By adding the `+1` offsets, the Devito compiler ensures the array accesses are logically aligned to the equation’s physical domain. For instance, the `TimeFunction` `u(t, x, y)` used in the example above has one point on each side of the `x` and `y` halo regions; if the user writes an expression including `u(t, x, y)` and `u(t, x + 2, y + 2)`, the compiler will ultimately generate `u[t, x + 1, y + 1]` and `u[t, x + 3, y + 3]`. When `x = y = 0`, therefore, the values `u[t, 1, 1]` and `u[t, 3, 3]` are fetched, representing the first and third points in the physical domain.

By default, the halo region has `space_order` points on each side of the space dimensions. Sometimes, these points may be unnecessary, or, depending on the partial differential equation being approximated, extra points may be needed.

``````u0 = TimeFunction(name='u0', grid=grid, space_order=0)
u0.data[:] = 1
print(u0.data_with_halo)``````
``````[[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]

[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]]``````
``````u2 = TimeFunction(name='u2', grid=grid, space_order=2)
u2.data[:] = 1
print(u2.data_with_halo)``````
``````[[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]]

[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 1. 1. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]]]``````

One can also pass a 3-tuple `(o, lp, rp)` instead of a single integer representing the discretization order. Here, `o` is the discretization order, while `lp` and `rp` indicate how many points are expected on left and right sides of a point of interest, respectivelly.

``u_new = TimeFunction(name='u_new', grid=grid, space_order=(4, 3, 1))``
``````u_new.data[:] = 1
print(u_new.data_with_halo)``````
``````[[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 0. 0. 0. 0.]]

[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 1. 1. 1. 0.]
[0. 0. 0. 0. 0. 0. 0.]]]``````

Let’s have a look at the generated code when using `u_new`.

``````equation = Eq(u_new.forward, u_new + 2)
op = Operator(equation, opt='noop')
print(op)``````
``````#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

struct dataobj
{
void *restrict data;
unsigned long * size;
unsigned long * npsize;
unsigned long * dsize;
int * hsize;
int * hofs;
int * oofs;
void * dmap;
} ;

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_new_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)
{
float (*restrict u_new)[u_new_vec->size[1]][u_new_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_new_vec->size[1]][u_new_vec->size[2]]) u_new_vec->data;

for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
u_new[t1][x + 3][y + 3] = u_new[t0][x + 3][y + 3] + 2;
}
}
STOP(section0,timers)
}

return 0;
}
``````

And finally, let’s run it, to convince ourselves that only the domain region values will be incremented at each timestep.

``````#NBVAL_IGNORE_OUTPUT
op.apply(time_M=2)
print(u_new.data_with_halo)``````
``Operator `Kernel` ran in 0.01 s``
``````[[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 5. 5. 5. 0.]
[0. 0. 0. 5. 5. 5. 0.]
[0. 0. 0. 5. 5. 5. 0.]
[0. 0. 0. 0. 0. 0. 0.]]

[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 7. 7. 7. 0.]
[0. 0. 0. 7. 7. 7. 0.]
[0. 0. 0. 7. 7. 7. 0.]
[0. 0. 0. 0. 0. 0. 0.]]]``````

The halo region, in turn, is surrounded by the padding region, which can be used for data alignment. By default, there is no padding. This can be changed by passing a suitable value to `padding`, as shown below:

``````u_pad = TimeFunction(name='u_pad', grid=grid, space_order=2, padding=(0,2,2))
op = Operator(equation, opt='noop')
print(op)``````
``````#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

struct dataobj
{
void *restrict data;
unsigned long * size;
unsigned long * npsize;
unsigned long * dsize;
int * hsize;
int * hofs;
int * oofs;
void * dmap;
} ;

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_pad_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)
{

for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
u_pad[t1][x + 2][y + 2] = u_pad[t0][x + 2][y + 2] + 2;
}
}
STOP(section0,timers)
}

return 0;
}
``````

Although in practice not very useful, with the (private) `_data_allocated` accessor one can see the entire domain + halo + padding region.

``print(u_pad._data_allocated)``
``````[[[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]]

[[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 2. 2. 2. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]]]``````

Above, the domain is filled with 2, the halo is filled with 1, and the padding is filled with 0.