<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The 2D diffusion equation on GPUs, in minutes
This notebook implements, for a given initial density profile, a solver for the 2D diffusion equation using an explicit finite difference scheme with 'do-nothing' conditions on the boundaries (and hence will not provide a reasonable solution once the profile has diffused to a boundary).
Step1:  Solver implementation
We start by creating a Cartesian Grid representing the computational domain
Step2:  To represent the density, we use a TimeFunction -- a scalar, discrete function encapsulating space- and time-varying data. We also use a Constant for the diffusion coefficient.
Step3:  The 2D diffusion equation is expressed as
Step4:  From this diffusion equation we derive our time-marching method -- at each timestep, we compute u at timestep t+1, which in the Devito language is represented by u.forward. Hence
Step5:  OK, it's time to let Devito generate code for our solver!
Step6:  Before executing the Operator we must first specify the initial density profile. Here, we place a "ring" with a constant density value in the center of the domain.
Step7:  We're now ready to execute the Operator. We run it with a diffusion coefficient of 0.5 and for a carefully chosen dt. Unless specified otherwise, the simulation runs for 199 timesteps as specified in the definition of u (i.e. the function was defined with save=200 the initial data + 199 new timesteps).
Step8:  Initial conditions and snapshots every 40 timesteps
Step9:  GPU-parallel solver
Let us now generate a GPU implementation of the same solver. It's actually straightforward!
Step10:  That's it! We can now run it exactly as before
Step11:  We should see a big performance difference between the two runs. We can also inspect op to see what Devito has generated to run on the GPU
<ASSISTANT_TASK:>
Python Code:
# Some imports we will need below
import numpy as np
from devito import *
import matplotlib.pyplot as plt 
%matplotlib inline
Explanation: The 2D diffusion equation on GPUs, in minutes
This notebook implements, for a given initial density profile, a solver for the 2D diffusion equation using an explicit finite difference scheme with 'do-nothing' conditions on the boundaries (and hence will not provide a reasonable solution once the profile has diffused to a boundary).
End of explanation
nx, ny = 100, 100
grid = Grid(shape=(nx, ny))
Explanation: Solver implementation
We start by creating a Cartesian Grid representing the computational domain:
End of explanation
u = TimeFunction(name='u', grid=grid, space_order=2, save=200)
c = Constant(name='c')
Explanation: To represent the density, we use a TimeFunction -- a scalar, discrete function encapsulating space- and time-varying data. We also use a Constant for the diffusion coefficient.
End of explanation
eqn = Eq(u.dt, c * u.laplace)
Explanation: The 2D diffusion equation is expressed as:
End of explanation
step = Eq(u.forward, solve(eqn, u.forward))
Explanation: From this diffusion equation we derive our time-marching method -- at each timestep, we compute u at timestep t+1, which in the Devito language is represented by u.forward. Hence:
End of explanation
op = Operator([step])
Explanation: OK, it's time to let Devito generate code for our solver!
End of explanation
xx, yy = np.meshgrid(np.linspace(0., 1., nx, dtype=np.float32),
                     np.linspace(0., 1., ny, dtype=np.float32))
r = (xx - .5)**2. + (yy - .5)**2.
# Inserting the ring
u.data[0, np.logical_and(.05 <= r, r <= .1)] = 1.
Explanation: Before executing the Operator we must first specify the initial density profile. Here, we place a "ring" with a constant density value in the center of the domain.
End of explanation
stats = op.apply(dt=5e-05, c=0.5)
Explanation: We're now ready to execute the Operator. We run it with a diffusion coefficient of 0.5 and for a carefully chosen dt. Unless specified otherwise, the simulation runs for 199 timesteps as specified in the definition of u (i.e. the function was defined with save=200 the initial data + 199 new timesteps).
End of explanation
plt.rcParams['figure.figsize'] = (20, 20)
for i in range(1, 6):
    plt.subplot(1, 6, i)
    plt.imshow(u.data[(i-1)*40])
plt.show()
Explanation: Initial conditions and snapshots every 40 timesteps
End of explanation
# Instead of `platform=nvidiaX`, you may run your Python code with
# the environment variable `DEVITO_PLATFORM=nvidiaX`
# We also need the `gpu-fit` option to tell Devito that `u` will definitely
# fit in the GPU memory. This is necessary every time a TimeFunction with
# `save != None` is used. Otherwise, Devito could generate code such that
# `u` gets streamed between the CPU and the GPU, but for this advanced 
# feature you will need `devitopro`.
op = Operator([step], platform='nvidiaX', opt=('advanced', {'gpu-fit': u}))
Explanation: GPU-parallel solver
Let us now generate a GPU implementation of the same solver. It's actually straightforward!
End of explanation
# Uncomment and run only if Devito was installed with GPU support.
# stats = op.apply(dt=5e-05, c=0.5)
Explanation: That's it! We can now run it exactly as before
End of explanation
print(op)
Explanation: We should see a big performance difference between the two runs. We can also inspect op to see what Devito has generated to run on the GPU
End of explanation
<END_TASK>