Example 5: Laplace equation

In this tutorial we will look at constructing the steady-state heat example using the Laplace equation. In contrast to the previous tutorials this example is entirely driven by the prescribed Dirichlet and Neumann boundary conditions, instead of an initial condition. We will also demonstrate how to use Devito to solve a steady-state problem without time derivatives and how to switch buffers explicitly without having to re-compile the kernel.

First, we again define our governing equation: $$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0$$

We are again discretizing second-order derivatives using a central difference scheme to construct a diffusion problem (see tutorial 3). This time we have no time-dependent term in our equation though, since there is no term $p_{i,j}^{n+1}$. This means that we are simply updating our field variable $p$ over and over again, until we have reached an equilibrium state. In a discretized form, after rearranging to update the central point $p_{i,j}^n$ we have $$p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}$$

And, as always, we first re-create the original implementation to see what we are aiming for. Here we initialize the field $p$ to $0$ and apply the following boundary conditions:

$p=0$ at $x=0$

$p=y$ at $x=2$

$\frac{\partial p}{\partial y}=0$ at $y=0, \ 1$

Developer note: The original tutorial stores the field data in the layout (ny, nx). Until now we have used (x, y) notation for creating our Devito examples, but for this one we will adopt the (y, x) layout for compatibility reasons.

OK, nice. Now, to re-create this example in Devito we need to look a little bit further under the hood. There are two things that make this different to the examples we covered so far:

• We have no time dependence in the p field, but we still need to advance the state of p in between buffers. So, instead of using TimeFunction objects that provide multiple data buffers for timestepping schemes, we will use Function objects that have no time dimension and only allocate a single buffer according to the space dimensions. However, since we are still implementing a pseudo-timestepping loop, we will need two objects, say p and pn, to act as alternating buffers.
• If we're using two different symbols to denote our buffers, any operator we create will only perform a single timestep. This is desired though, since we need to check a convergence criterion outside of the main stencil update to determine when we stop iterating. As a result we will need to call the operator repeatedly after instantiating it outside the convergence loop.

So, how do we make sure our operator doesn't accidentally overwrite values in the same buffer? Well, we can again let SymPy reorganise our Laplace equation based on pn to generate the stencil, but when we create the update expression, we set the LHS to our second buffer variable p.

Now we can add our boundary conditions. We have already seen how to prescribe constant Dirichlet BCs by simply setting values using the low-level notation. This time we will go a little further by setting a prescribed profile, which we create first as a custom 1D symbol and supply with the BC values. For this we need to create a Function object that has a different shape than our general grid, so instead of the grid we provide an explicit pair of dimension symbols and the according shape for the data.

Now we can create a set of expressions for the BCs again, where we wet prescribed values on the right and left of our grid. For the Neumann BCs along the top and bottom boundaries we simply copy the second row from the outside into the outermost row, just as the original tutorial did. Using these expressions and our stencil update we can now create an operator.

We can now use this single-step operator repeatedly in a Python loop, where we can arbitrarily execute other code in between invocations. This allows us to update our L1 norm and check for convergence. Using our pre-compiled operator now comes down to a single function call that supplies the relevant data symbols. One thing to note is that we now do exactly the same thing as the original NumPy loop, in that we deep-copy the data between each iteration of the loop, which we will look at after this.

One crucial detail about the code above is that the deep data copy between iterations will really hurt performance if we were to run this on a large grid. However, we have already seen how we can match data symbols to symbolic names when calling the pre-compiled operator, which we can now use to actually switch the roles of pn and p between iterations, eg. op(p=pn, pn=p). Thus, we can implement a simple buffer-switching scheme by simply testing for odd and even time-steps, without ever having to shuffle data around.