Solving the 2D Heat Equation using DevitoΒΆ
Consider the 2D Heat Equation defined by
This equation can be solved numerically using Finite Difference approximations.
First of all, we need to allocate the grid and set its initial condition:
from devito import TimeData
u = TimeData(name = 'u', shape = 'nx, ny', time_order = 1, space_order = 2)
u.data[0, :] = ui[:]
TimeData is a Devito data object used to store and manage time-varying data.
We initialise our grid to be of size nx * ny
for some nx
and ny
. time_order
and space_order
represent the discretization
order for time and space respectively. The initial configuration is given as a
numpy.array
in u.data
.
The next step is to generate the stencil to be solved by a
devito.operator.Operator
:
from devito import Operator
from sympy import Eq, solve, symbols
a, h, s = symbols("a h s")
eqn = Eq(u.dt, a * (u.dx2 + u.dy2))
stencil = solve(eqn, u.forward)[0]
The stencil is generated according to Devito conventions. It uses a sympy
equation to represent the 2D Heat equation and store it in eqn
.
Devito makes easy to represent the equation by providing properties dt
,
dx2
, and dx2
that represent the derivatives.
We then generate the stencil by solving eqn
for u.forward
, a
symbol for the time-forward state of the function.
We plug the stencil in an Operator, as shown, and define the values of the
thermal conductivity a
, the spacing between cells h
and the
temporal spacing s
.:
op = Operator(Eq(u.forward, stencil),
subs={h: spacing, s: dt, a: tc})
To execute the generated Operator, we simply call op.apply(u=u,
t=timesteps)
. The results will then be found in u.data[1, :]
.
For a complete example of this code, please see examples/diffusion/example_diffusion.py. A more comprehensive set of CFD tutorials based on the excellent 12 steps to Navier-Stokes tutorial is currently under construction and will be published here soon.