#!/usr/bin/env python
# coding: utf-8

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

# 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)


# 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)


# 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)


# 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)


# In[7]:


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


# 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)


# 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)


# 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)


# 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]:


u_pad = TimeFunction(name='u_pad', grid=grid, space_order=2, padding=(0,2,2))
u_pad._data_allocated[:] = 0
u_pad.data_with_halo[:] = 1
u_pad.data[:] = 2
equation = Eq(u_pad.forward, u_pad + 2)
op = Operator(equation, opt='noop')
print(op)


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

# In[13]:


print(u_pad._data_allocated)


# Above, the __domain__ is filled with __2__, the __halo__ is filled with __1__, and the __padding__ is filled with __0__.