```
from examples.cfd import plot_field, init_hat
import numpy as np
%matplotlib inline
# Some variable declarations
= 100
nx = 100
ny = 1000
nt
= 0.15 #the value of base viscosity
nu
= 1 # Used for field definition
offset
= np.full((nx, ny), nu) # Initialize viscosity
visc //4-offset:nx//4+offset, 1:-1] = 0.0001 # Adding a material with different viscosity
visc[nx1:-1,nx//4-offset:nx//4+offset ] = 0.0001
visc[3*nx//4-offset:3*nx//4+offset, 1:-1] = 0.0001
visc[
= visc[1:-1,1:-1]
visc_nb
= 2. / (nx - 1)
dx = 2. / (ny - 1)
dy = .25
sigma = sigma * dx * dy / nu
dt
# Initialize our field
# Initialise u with hat function
= np.empty((nx, ny))
u_init =u_init, dx=dx, dy=dy, value=1)
init_hat(field10:-10, 10:-10] = 1.5
u_init[
= 2.5 # zmax for plotting zmax
```

# Example 3 , part B: Diffusion for non uniform material properties

In this example we will look at the diffusion equation for non uniform material properties and how to handle second-order derivatives. For this, we will reuse Devito’s `.laplace`

short-hand expression outlined in the previous example and demonstrate it using the examples from step 7 of the original tutorial. This example is an enhancement of `03_diffusion`

in terms of having non-uniform viscosity as opposed to the constant \(\nu\). This example introduces the use of `Function`

in order to create this non-uniform space.

So, the equation we are now trying to implement is

\[\frac{\partial u}{\partial t} = \nu(x,y) \frac{\partial ^2 u}{\partial x^2} + \nu(x,y) \frac{\partial ^2 u}{\partial y^2}\]

In our case \(\nu\) is not uniform and \(\nu(x,y)\) represents spatially variable viscosity. To discretize this equation we will use central differences and reorganizing the terms yields

\[\begin{align} u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu(x,y) \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \\ &+ \frac{\nu(x,y) \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n) \end{align}\]

As usual, we establish our baseline experiment by re-creating some of the original example runs. So let’s start by defining some parameters.

We now set up the diffusion operator as a separate function, so that we can re-use if for several runs.

```
def diffuse(u, nt ,visc):
for n in range(nt + 1):
= u.copy()
un 1:-1, 1:-1] = (un[1:-1,1:-1] +
u[*dt / dy**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
visc*dt / dx**2 * (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
visc0, :] = 1
u[-1, :] = 1
u[0] = 1
u[:, -1] = 1 u[:,
```

Now let’s take this for a spin. In the next two cells we run the same diffusion operator for a varying number of timesteps to see our “hat function” dissipate to varying degrees.

```
#NBVAL_IGNORE_OUTPUT
# Plot material according to viscosity, uncomment to plot
import matplotlib.pyplot as plt
='Greys', interpolation='nearest')
plt.imshow(visc_nb, cmap
# Field initialization
= u_init
u
print ("Initial state")
=zmax)
plot_field(u, zmax
diffuse(u, nt , visc_nb )print ("After", nt, "timesteps")
=zmax)
plot_field(u, zmax
diffuse(u, nt, visc_nb)print ("After another", nt, "timesteps")
=zmax) plot_field(u, zmax
```

`Initial state`

`After 1000 timesteps`

`After another 1000 timesteps`

You can notice that the area with lower viscosity is not diffusing its heat as quickly as the area with higher viscosity.

```
#NBVAL_IGNORE_OUTPUT
# Field initialization
= u_init
u
diffuse(u, nt , visc_nb)print ("After", nt, "timesteps")
=zmax) plot_field(u, zmax
```

`After 1000 timesteps`

Excellent. Now for the Devito part, we need to note one important detail to our previous examples: we now have a second-order derivative. So, when creating our `TimeFunction`

object we need to tell it about our spatial discretization by setting `space_order=2`

. We also use the notation `u.laplace`

outlined previously to denote all second order derivatives in space, allowing us to reuse this code for 2D and 3D examples.

```
from devito import Grid, TimeFunction, Eq, solve, Function
from sympy.abc import a
# Initialize `u` for space order 2
= Grid(shape=(nx, ny), extent=(2., 2.))
grid
# Create an operator with second-order derivatives
= Function(name='a',grid = grid) # Define as Function
a = visc # Pass the viscosity in order to be used in the operator.
a.data[:]
= TimeFunction(name='u', grid=grid, space_order=2)
u
# Create an equation with second-order derivatives
= Eq(u.dt, a * u.laplace)
eq = solve(eq, u.forward)
stencil = Eq(u.forward, stencil)
eq_stencil
eq_stencil
```

`Eq(u(t + dt, x, y), dt*((Derivative(u(t, x, y), (x, 2)) + Derivative(u(t, x, y), (y, 2)))*a(x, y) + u(t, x, y)/dt))`

Great. Now all that is left is to put it all together to build the operator and use it on our examples. For illustration purposes we will do this in one cell, including update equation and boundary conditions.

```
#NBVAL_IGNORE_OUTPUT
from devito import Operator, Constant, Eq, solve, Function
# Reset our data field and ICs
=u.data[0], dx=dx, dy=dy, value=1.)
init_hat(field
# Field initialization
0] = u_init
u.data[
# Create an operator with second-order derivatives
= Function(name='a',grid = grid)
a = visc
a.data[:]
= Eq(u.dt, a * u.laplace, subdomain=grid.interior)
eq = solve(eq, u.forward)
stencil = Eq(u.forward, stencil)
eq_stencil
# Create boundary condition expressions
= grid.dimensions
x, y = grid.stepping_dim
t = [Eq(u[t+1, 0, y], 1.)] # left
bc += [Eq(u[t+1, nx-1, y], 1.)] # right
bc += [Eq(u[t+1, x, ny-1], 1.)] # top
bc += [Eq(u[t+1, x, 0], 1.)] # bottom
bc
= Operator([eq_stencil] + bc)
op =nt, dt=dt, a = a)
op(time
print ("After", nt, "timesteps")
0], zmax=zmax)
plot_field(u.data[
=nt, dt=dt, a = a)
op(timeprint ("After another", nt, "timesteps")
0], zmax=zmax) plot_field(u.data[
```

```
Operator `Kernel` ran in 0.01 s
After 1000 timesteps
```

```
Operator `Kernel` ran in 0.01 s
After another 1000 timesteps
```