### Example 3: Diffusion 

In this example we will look at the diffusion equation and how to handle second-order derivatives. For this we will introduce Devito's `.laplace` short-hand expression and demonstrate it using the examples from step 7 of the original tutorial.

So, the equation we are now trying to implement is

$$\frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2}$$

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 \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \\
&+ \frac{\nu \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.


```python
from examples.cfd import plot_field, init_hat
import numpy as np
%matplotlib inline

# Some variable declarations
nx = 31
ny = 31
nt = 17
nu = .05
dx = 2. / (nx - 1)
dy = 2. / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu
```

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


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

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


```python
#NBVAL_IGNORE_OUTPUT

# Initialise u with hat function
u = np.empty((nx, ny))
init_hat(field=u, dx=dx, dy=dy, value=2.)
plot_field(u, zmax=2.5)

diffuse(u, nt=10)
plot_field(u, zmax=2.5)
```


```python
#NBVAL_IGNORE_OUTPUT

# Initialise u with hat function
u = np.empty((nx, ny))
init_hat(field=u, dx=dx, dy=dy, value=2.)
 
diffuse(u, nt=50)
plot_field(u, zmax=2.5)
```

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 discretisation by setting `space_order=2`. Only then can we use the shorthand notation `u.dx2` and `u.dx2` to denote second derivatives.


```python
from devito import Grid, TimeFunction, Eq
from sympy import solve
from sympy.abc import a

# Initialize `u` for space order 2
grid = Grid(shape=(nx, ny), extent=(2., 2.))
u = TimeFunction(name='u', grid=grid, space_order=2)

# Create an equation with second-order derivatives
eq = Eq(u.dt, a * (u.dx2 + u.dy2))
stencil = solve(eq, u.forward)[0]
eq_stencil = Eq(u.forward, stencil)

print(eq_stencil)
```

    Eq(u(t + dt, x, y), (a*dt*h_x**2*(-2*u(t, x, y) + u(t, x, y - h_y) + u(t, x, y + h_y)) + a*dt*h_y**2*(-2*u(t, x, y) + u(t, x - h_x, y) + u(t, x + h_x, y)) + h_x**2*h_y**2*u(t, x, y))/(h_x**2*h_y**2))


Now, there is another trick here! Note how the above formulation explicitly uses `u.dx2` and `u.dy2` to denote the laplace operator, which makes this equation dependent on the spatial dimension. We can instead use the notation `u.laplace` to denote all second derivatives in space, allowing us to reuse this code for 2D and 3D examples.


```python
eq = Eq(u.dt, a * u.laplace)
stencil = solve(eq, u.forward)[0]
eq_stencil = Eq(u.forward, stencil)

print(eq_stencil)
```

    Eq(u(t + dt, x, y), (a*dt*h_x**2*(-2*u(t, x, y) + u(t, x, y - h_y) + u(t, x, y + h_y)) + a*dt*h_y**2*(-2*u(t, x, y) + u(t, x - h_x, y) + u(t, x + h_x, y)) + h_x**2*h_y**2*u(t, x, y))/(h_x**2*h_y**2))


Great. Now all that is left is to pu 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.


```python
#NBVAL_IGNORE_OUTPUT
from sympy import solve
from devito import Operator, Constant, Eq, INTERIOR

# Reset our data field and ICs
init_hat(field=u.data[0], dx=dx, dy=dy, value=2.)

# Create an operator with second-order derivatives
a = Constant(name='a')
eq = Eq(u.dt, a * u.laplace, region=INTERIOR)
stencil = solve(eq, u.forward)[0]
eq_stencil = Eq(u.forward, stencil)

# Create boundary condition expressions
x, y = grid.dimensions
t = grid.stepping_dim
bc = [Eq(u.indexed[t+1, 0, y], 1.)]  # left
bc += [Eq(u.indexed[t+1, nx-1, y], 1.)]  # right
bc += [Eq(u.indexed[t+1, x, ny-1], 1.)]  # top
bc += [Eq(u.indexed[t+1, x, 0], 1.)]  # bottom

op = Operator([eq_stencil] + bc)
op(time=10, dt=dt, a=nu)

plot_field(u.data[0], zmax=2.5)
```

And now let's use the same operator again to show the more diffused field. In fact, instead of resetting the previously used object `u`, we can also create a new `TimeFunction` object and tell our operator to use this.


```python
#NBVAL_IGNORE_OUTPUT
u2 = TimeFunction(name='u2', grid=grid, space_order=2)
init_hat(field=u2.data[0], dx=dx, dy=dy, value=2.)

op(u=u2, time=50, dt=dt, a=nu)

plot_field(u2.data[0], zmax=2.5)
```

21d7924b03928bafb32a6c78de4dc89295b271c0
ipynb
Jupyter Notebook
devito/examples/cfd/03_diffusion.ipynb
LukasMosser/stochastic_seismic_waveform_inversion
4976c3b9a39b8d246d3d220056f235df6fc7dbb3
MIT
2019-08-01T19:08:10.000Z
2022-03-20T05:04:06.000Z
devito/examples/cfd/03_diffusion.ipynb
esgomi/gan_seismic_waveform_inversion
4976c3b9a39b8d246d3d220056f235df6fc7dbb3
MIT
2019-12-05T13:00:14.000Z
2019-12-05T13:00:14.000Z
devito/examples/cfd/03_diffusion.ipynb
esgomi/gan_seismic_waveform_inversion
4976c3b9a39b8d246d3d220056f235df6fc7dbb3
MIT
2019-10-03T14:47:48.000Z
2021-09-22T14:06:47.000Z
Qwen/Qwen-72B
1. YES
2. YES