```
In [1]:
```from devito import *

The main objective of this tutorial is to demonstrate how Devito and its SymPy-powered symbolic API can be used to solve partial differential equations using the finite difference method with highly optimized stencils in a few lines of Python. We demonstrate how computational stencils can be derived directly from the equation in an automated fashion and how Devito can be used to generate and execute, at runtime, the desired numerical scheme in the form of optimized C code.

Before we can begin creating finite-difference (FD) stencils we will need to give Devito a few details regarding the computational domain within which we wish to solve our problem. For this purpose we create a `Grid`

object that stores the physical `extent`

(the size) of our domain and knows how many points we want to use in each dimension to discretise our data.

```
In [2]:
```grid = Grid(shape=(5, 6), extent=(1., 1.))
grid

```
Out[2]:
```

To express our equation in symbolic form and discretise it using finite differences, Devito provides a set of `Function`

types. A `Function`

object:

- Behaves like a
`sympy.Function`

symbol - Manages data associated with the symbol

To get more information on how to create and use a `Function`

object, or any type provided by Devito, we can take a look at the documentation.

```
In [3]:
```print(Function.__doc__)

```
```

`name`

or `grid`

when creating `Function`

objects.

```
In [4]:
```f = Function(name='f', grid=grid)
f

```
Out[4]:
```

```
In [5]:
```f.data

```
Out[5]:
```

`Function`

objects use the spatial dimensions `(x, y)`

for 2D grids and `(x, y, z)`

for 3D grids. To solve a PDE over several timesteps a time dimension is also required by our symbolic function. For this Devito provides an additional function type, the `TimeFunction`

, which incorporates the correct dimension along with some other intricacies needed to create a time stepping scheme.

```
In [6]:
```g = TimeFunction(name='g', grid=grid)
g

```
Out[6]:
```

`TimeFunction`

is `1`

, the shape of `f`

is `(2, 5, 6)`

, i.e. Devito has allocated two buffers to represent `g(t, x, y)`

and `g(t + dt, x, y)`

:

```
In [7]:
```g.shape

```
Out[7]:
```

The functions we have created so far all act as `sympy.Function`

objects, which means that we can form symbolic derivative expressions from them. Devito provides a set of shorthand expressions (implemented as Python properties) that allow us to generate finite differences in symbolic form. For example, the property `f.dx`

denotes $\frac{\partial}{\partial x} f(x, y)$ - only that Devito has already discretised it with a finite difference expression. There are also a set of shorthand expressions for left (backward) and right (forward) derivatives:

Derivative | Shorthand | Discretised | Stencil |
---|---|---|---|

$\frac{\partial}{\partial x}f(x, y)$ (right) | `f.dxr` |
$\frac{f(x+h_x,y)}{h_x} - \frac{f(x,y)}{h_x}$ | |

$\frac{\partial}{\partial x}f(x, y)$ (left) | `f.dxl` |
$\frac{f(x,y)}{h_x} - \frac{f(x-h_x,y)}{h_x}$ |

A similar set of expressions exist for each spatial dimension defined on our grid, for example `f.dy`

and `f.dyl`

. Obviously, one can also take derivatives in time of `TimeFunction`

objects. For example, to take the first derivative in time of `g`

you can simply write:

```
In [8]:
```g.dt

```
Out[8]:
```

```
In [9]:
```g.dt.evaluate

```
Out[9]:
```

`g(t+dt, x, y)`

and `g(t-dt, x, y)`

.

```
In [10]:
```g.forward

```
Out[10]:
```

```
In [11]:
```g.backward

```
Out[11]:
```

And of course, there's nothing to stop us taking derivatives on these objects:

```
In [12]:
```g.forward.dt

```
Out[12]:
```

```
In [13]:
```g.forward.dy

```
Out[13]:
```

**Note:** The following example is derived from step 5 in the excellent tutorial series CFD Python: 12 steps to Navier-Stokes.

In this simple example we will show how to derive a very simple convection operator from a high-level description of the governing equation. We will go through the process of deriving a discretised finite difference formulation of the state update for the field variable $u$, before creating a callable `Operator`

object. Luckily, the automation provided by SymPy makes the derivation very nice and easy.

The governing equation we want to implement is the linear convection equation: $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\partial y} = 0.$$

Before we begin, we must define some parameters including the grid, the number of timesteps and the timestep size. We will also initialize our velocity `u`

with a smooth field:

```
In [14]:
```from examples.cfd import init_smooth, plot_field
nt = 100 # Number of timesteps
dt = 0.2 * 2. / 80 # Timestep size (sigma=0.2)
c = 1 # Value for c
# Then we create a grid and our function
grid = Grid(shape=(81, 81), extent=(2., 2.))
u = TimeFunction(name='u', grid=grid)
# We can now set the initial condition and plot it
init_smooth(field=u.data[0], dx=grid.spacing[0], dy=grid.spacing[1])
init_smooth(field=u.data[1], dx=grid.spacing[0], dy=grid.spacing[1])
plot_field(u.data[0])

```
```

Next, we wish to discretise our governing equation so that a functional `Operator`

can be created from it. We begin by simply writing out the equation as a symbolic expression, while using shorthand expressions for the derivatives provided by the `Function`

object. This will create a symbolic object of the dicretised equation.

Using the Devito shorthand notation, we can express the governing equations as:

```
In [15]:
```eq = Eq(u.dt + c * u.dxl + c * u.dyl)
eq

```
Out[15]:
```

`solve`

, built on top of SymPy's `solve`

, to rearrange our equation so that it represents a valid state update for $u$. Here, we use `solve`

to create a valid stencil for our update to `u(t+dt, x, y)`

:

```
In [16]:
```stencil = solve(eq, u.forward)
update = Eq(u.forward, stencil)
update

```
Out[16]:
```

The right-hand side of the 'update' equation should be a stencil of the shape

Once we have created this 'update' expression, we can create a Devito `Operator`

. This `Operator`

will basically behave like a Python function that we can call to apply the created stencil over our associated data, as long as we provide all necessary unknowns. In this case we need to provide the number of timesteps to compute via the keyword `time`

and the timestep size via `dt`

(both have been defined above):

```
In [17]:
```op = Operator(update)
op(time=nt+1, dt=dt)
plot_field(u.data[0])

```
```

`Operator`

, it will automatically generate and compile the optimized C code. We can look at this code (noting that this is not a requirement of executing it) via:

```
In [18]:
```print(op.ccode)

```
```

In the above example only a combination of first derivatives was present in the governing equation. However, second (or higher) order derivatives are often present in scientific problems of interest, notably any PDE modeling diffusion. To generate second order derivatives we must give the `devito.Function`

object another piece of information: the desired discretisation of the stencil(s).

First, lets define a simple second derivative in `x`

, for which we need to give $u$ a `space_order`

of (at least) `2`

. The shorthand for this second derivative is `u.dx2`

.

```
In [19]:
```u = TimeFunction(name='u', grid=grid, space_order=2)
u.dx2

```
Out[19]:
```

```
In [20]:
```u.dx2.evaluate

```
Out[20]:
```

We can increase the discretisation arbitrarily if we wish to specify higher order FD stencils:

```
In [21]:
```u = TimeFunction(name='u', grid=grid, space_order=4)
u.dx2

```
Out[21]:
```

```
In [22]:
```u.dx2.evaluate

```
Out[22]:
```

`Grid`

object and use this to re-define our `Function`

's:

```
In [23]:
```grid_3d = Grid(shape=(5, 6, 7), extent=(1., 1., 1.))
u = TimeFunction(name='u', grid=grid_3d, space_order=2)
u

```
Out[23]:
```

`u`

with a different `space_order`

argument to change the discretisation order of the stencil expression created. For example, we can derive an expression of the 12th-order Laplacian $\nabla^2 u$:

```
In [24]:
```u = TimeFunction(name='u', grid=grid_3d, space_order=12)
u.laplace

```
Out[24]:
```

The same expression could also have been generated explicitly via:

```
In [25]:
```u.dx2 + u.dy2 + u.dz2

```
Out[25]:
```

```
In [26]:
```u = TimeFunction(name='u', grid=grid, space_order=2)
v = TimeFunction(name='v', grid=grid, space_order=2, time_order=2)

```
In [27]:
```v.dt2 + u.laplace

```
Out[27]:
```

```
In [28]:
```(v.dt2 + u.laplace).dx2

```
Out[28]:
```

Which can, depending on the chosen discretisation, lead to fairly complex stencils:

```
In [29]:
```(v.dt2 + u.laplace).dx2.evaluate

```
Out[29]:
```