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

``````

In [1]:

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:

``````

In [2]:

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.

``````

In [3]:

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.

``````

In [4]:

print(op)

``````
``````

#define _POSIX_C_SOURCE 200809L
#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

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

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)
{
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))
{
struct timeval start_section0, end_section0;
gettimeofday(&start_section0, NULL);
/* Begin 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;
}
}
/* End section0 */
gettimeofday(&end_section0, NULL);
timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
}
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.

``````

In [5]:

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.

``````

In [6]:

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.]]]

``````
``````

In [7]:

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.

``````

In [8]:

u_new = TimeFunction(name='u_new', grid=grid, space_order=(4, 3, 1))

``````
``````

In [9]:

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`.

``````

In [10]:

equation = Eq(u_new.forward, u_new + 2)
op = Operator(equation, opt='noop')
print(op)

``````
``````

#define _POSIX_C_SOURCE 200809L
#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

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

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_new_vec, const int time_M, const int time_m, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)
{
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))
{
struct timeval start_section0, end_section0;
gettimeofday(&start_section0, NULL);
/* Begin 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;
}
}
/* End section0 */
gettimeofday(&end_section0, NULL);
timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
}
return 0;
}

``````

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

``````

In [11]:

#NBVAL_IGNORE_OUTPUT
op.apply(time_M=2)
print(u_new.data_with_halo)

``````
``````

Operator `Kernel` run in 0.00 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:

``````

In [12]:

op = Operator(equation, opt='noop')
print(op)

``````
``````

#define _POSIX_C_SOURCE 200809L
#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

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

struct profiler
{
double section0;
} ;

int Kernel(struct dataobj *restrict u_pad_vec, const int time_M, const int time_m, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)
{
for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
struct timeval start_section0, end_section0;
gettimeofday(&start_section0, NULL);
/* Begin 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;
}
}
/* End section0 */
gettimeofday(&end_section0, NULL);
timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
}
return 0;
}

``````

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

``````

In [13]:

``````
``````

[[[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.