#!/usr/bin/env python # coding: utf-8 # ### Example 6: Poisson equation # # In the previous example we have seen how to solve the steady-state Laplace equation with a prescribed boundary condition and BCs that contain derivatives (Neumann BC). For this we used a pseudo-timestepping to illustrate the use of two explicit data symbols of type `Function` to derive at a steady state. In this example we are going to solve the Poisson equation (step 10 of the original tutorial series) with an initial condition to demonstrate how we can combine time-dependent functions of type `TimeFunction` with time-constant `Function` symbols. We will again implement a Python-driven outer loop for the diffusion kernel, before demonstrating how we can change the expression defining our operator to internalize the pseudo-timestepping and thus get the full performance advantage of compiled kernel loops. # # Our governing equation is the Poisson equation with a source term as the RHS: # $$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$ # # Discretizing and rearranging our stencil expression yields # $$p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}$$ # # The key difference to the Laplace example in the previous tutorial is that this example has simple Dirichlet boundary conditions $p = 0$ on all boundaries and is instead driven by the source term $b$ defined as # # $b_{i,j}=100$ at $i=nx/4,\ j=ny/4$ # # $b_{i,j}=-100$ at $i=3/4 nx,\ j=3/4 ny$ # # $b_{i,j}=0$ everywhere else. # # And of course, we will demonstrate the original example first using pure NumPy notation. # In[1]: from examples.cfd import plot_field, init_hat import numpy as np get_ipython().run_line_magic('matplotlib', 'inline') # Some variable declarations nx = 50 ny = 50 nt = 100 xmin = 0. xmax = 2. ymin = 0. ymax = 1. dx = (xmax - xmin) / (nx - 1) dy = (ymax - ymin) / (ny - 1) # In[2]: # Initialization p = np.zeros((nx, ny)) pd = np.zeros((nx, ny)) b = np.zeros((nx, ny)) # Source b[int(nx / 4), int(ny / 4)] = 100 b[int(3 * nx / 4), int(3 * ny / 4)] = -100 # In[3]: get_ipython().run_cell_magic('time', '', '#NBVAL_IGNORE_OUTPUT\nfor it in range(nt):\n pd = p.copy()\n p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +\n (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -\n b[1:-1, 1:-1] * dx**2 * dy**2) / \n (2 * (dx**2 + dy**2)))\n\n p[0, :] = 0\n p[nx-1, :] = 0\n p[:, 0] = 0\n p[:, ny-1] = 0\n \n') # In[4]: #NBVAL_IGNORE_OUTPUT plot_field(p, xmax=xmax, ymax=ymax, view=(30, 225)) # We can now pretty much use our previous implementation, although we will use `pd` instead of `pn` for consistency with the original. Our boundary conditions are even simpler since they are $0$ everywhere. Our source term is neatly wrapped in the symbol `b` and we can again use a Python-driven timestepping loop with switching buffers to repeatedly execute the operator we create. # In[5]: from devito import Grid, Function, TimeFunction, Operator, configuration, Eq, solve # Silence the runtime performance logging configuration['log-level'] = 'ERROR' grid = Grid(shape=(nx, ny), extent=(xmax, ymax)) p = Function(name='p', grid=grid, space_order=2) pd = Function(name='pd', grid=grid, space_order=2) p.data[:] = 0. pd.data[:] = 0. # Initialise the source term `b` b = Function(name='b', grid=grid) b.data[:] = 0. b.data[int(nx / 4), int(ny / 4)] = 100 b.data[int(3 * nx / 4), int(3 * ny / 4)] = -100 # Create Laplace equation base on `pd` eq = Eq(pd.laplace, b, subdomain=grid.interior) # Let SymPy solve for the central stencil point stencil = solve(eq, pd) # # Now we let our stencil populate our second buffer `p` eq_stencil = Eq(p, stencil) # Create boundary condition expressions x, y = grid.dimensions t = grid.stepping_dim bc = [Eq(p[x, 0], 0.)] bc += [Eq(p[x, ny-1], 0.)] bc += [Eq(p[0, y], 0.)] bc += [Eq(p[nx-1, y], 0.)] # Now we can build the operator that we need op = Operator([eq_stencil] + bc) # In[6]: get_ipython().run_cell_magic('time', '', '#NBVAL_IGNORE_OUTPUT\n\n# Run the outer loop explicitly in Python\nfor i in range(nt):\n # Determine buffer order\n if i % 2 == 0:\n _p = p\n _pd = pd\n else:\n _p = pd\n _pd = p\n\n # Apply operator\n op(p=_p, pd=_pd)\n') # In[7]: #NBVAL_IGNORE_OUTPUT # Plot result plot_field(p.data, xmax=xmax, ymax=ymax, view=(30, 225)) # Nice, we get the same spikes as the pure NumPy implementation. But we still drive the time loop though Python, which is slow. So let's try and utilise Devito's time-stepping mechanism that we saw earlier for this "pseudo-timestepping". # # The trick here is to again use `TimeFunction` symbols, despite the lack of a time-dependence. This will cause Devito to allocate two grid buffers instead of one, which we can directly address them via the terms `p` and `p.forward`. This will now create an internal time loop in the kernel that we can again control by supplying the `time` argument to the create operator during invocation. # # One caveat, however, is that a `TimeFunction` symbol has an implicit time dimension `t`, so the symbol is denoted as `p(t, x, y)`. This entails that, when creating our boundary condition expressions, we now need honour the leading time dimension in the "low-level" notation, and ensure that we update the index `t + 1` (the equivalent to `u.forward`). # In[8]: # Now with Devito we will turn `p` into `TimeFunction` object # to make all the buffer switching implicit p = TimeFunction(name='p', grid=grid, space_order=2) p.data[:] = 0. # Initialise the source term `b` b = Function(name='b', grid=grid) b.data[:] = 0. b.data[int(nx / 4), int(ny / 4)] = 100 b.data[int(3 * nx / 4), int(3 * ny / 4)] = -100 # Create Laplace equation base on `p` eq = Eq(p.laplace, b) # Let SymPy solve for the central stencil point stencil = solve(eq, p) # Let our stencil populate the buffer `p.forward` eq_stencil = Eq(p.forward, stencil) # Create boundary condition expressions # Note that we now add an explicit "t + 1" for the time dimension. bc = [Eq(p[t + 1, x, 0], 0.)] bc += [Eq(p[t + 1, x, ny-1], 0.)] bc += [Eq(p[t + 1, 0, y], 0.)] bc += [Eq(p[t + 1, nx-1, y], 0.)] # In[9]: #NBVAL_IGNORE_OUTPUT configuration['log-level'] = 'ERROR' # Create and execute the operator for a number of timesteps op = Operator([eq_stencil] + bc) get_ipython().run_line_magic('time', 'op(time=nt)') # In[10]: #NBVAL_IGNORE_OUTPUT plot_field(p.data[0], xmax=xmax, ymax=ymax, view=(30, 225))