Laplace’s Equation in Haskell: Using a DSL for Stencils

Introduction

Suppose we have a square thin plate of metal and we hold each of edges at a temperature which may vary along the edge but is fixed for all time. After some period depending on the conductivity of the metal, the temperature at every point on the plate will have stabilised. What is the temperature at any point?

We can calculate this using by solving Laplace’s equation in 2 dimensions. Apart from the preceeding motivation, a more compelling reason for doing so is that it is a moderately simple equation, in so far as partial differential equations are simple, that has been well studied for centuries.

In Haskell terms this gives us the opportunity to use the repa library and use hmatrix which is based on Lapack (as well as other libraries) albeit hmatrix only for illustratative purposes.

I had originally intended this blog to contain a comparison repa’s performance against an equivalent C program even though this has already been undertaken by the repa team in their various publications. And indeed it is still my intention to produce such a comparision. However, as I investigated further, it turned out a fair amount of comparison work has already been done by a team from Intel which suggests there is currently a performance gap but one which is not so large that it outweighs the other benefits of Haskell.

To be more specific, one way in which using repa stands out from the equivalent C implementation is that it gives a language in which we can specify the stencil being used to solve the equation. As an illustration we substitute the nine point method for the five point method merely by changing the stencil.

A Motivating Example: The Steady State Heat Equation

Fourier’s law states that the rate of heat transfer or the flux is proportional to the negative temperature gradient, as heat flows from hot to cold, and further that it flows in the direction of greatest temperature change. We can write this as

where is the temperature at any given point on the plate and is the conductivity of the metal.

Moreover, we know that for any region on the plate, the total amount of heat flowing in must be balanced by the amount of heat flowing out. We can write this as

Substituting the first equation into the second we obtain Laplace’s equation

For example, suppose we hold the temperature of the edges of the plate as follows

then after some time the temperature of the plate will be as shown in the heatmap below.

Notes:

1. Red is hot.

2. Blue is cold.

3. The heatmap is created by a finite difference method described below.

4. The -axis points down (not up) i.e. is at the bottom, reflecting the fact that we are using an array in the finite difference method and rows go down not up.

5. The corners are grey because in the five point finite difference method these play no part in determining temperatures in the interior of the plate.

Colophon

Since the book I am writing contains C code (for performance comparisons), I need a way of being able to compile and run this code and include it “as is” in the book. Up until now, all my blog posts have contained Haskell and so I have been able to use BlogLiteratelyD which allows me to include really nice diagrams. But clearly this tool wasn’t really designed to handle other languages (although I am sure it could be made to do so).

Using pandoc’s scripting capability with the small script provided

#!/usr/bin/env runhaskell
import Text.Pandoc.JSON

doInclude :: Block -> IO Block
doInclude cb@(CodeBlock ("verbatim", classes, namevals) contents) =
  case lookup "include" namevals of
       Just f     -> return . (\x -> Para [Math DisplayMath x]) =<< readFile f
       Nothing    -> return cb
doInclude cb@(CodeBlock (id, classes, namevals) contents) =
  case lookup "include" namevals of
       Just f     -> return . (CodeBlock (id, classes, namevals)) =<< readFile f
       Nothing    -> return cb
doInclude x = return x

main :: IO ()
main = toJSONFilter doInclude

I can then include C code blocks like this

~~~~ {.c include="Chap1a.c"}
~~~~

And format the whole document like this

pandoc -s Chap1.lhs --filter=./Include -t markdown+lhs > Chap1Expanded.lhs
BlogLiteratelyD Chap1Expanded.lhs > Chap1.html

Sadly, the C doesn’t get syntax highlighting but this will do for now.

PS Sadly, WordPress doesn’t seem to be able to handle \color{red} and \color{blue} in LaTeX so there are some references to blue and red which do not render.

Acknowledgements

A lot of the code for this post is taken from the repa package itself. Many thanks to the repa team for providing the package and the example code.

Haskell Preamble

> {-# OPTIONS_GHC -Wall                      #-}
> {-# OPTIONS_GHC -fno-warn-name-shadowing   #-}
> {-# OPTIONS_GHC -fno-warn-type-defaults    #-}
> {-# OPTIONS_GHC -fno-warn-unused-do-bind   #-}
> {-# OPTIONS_GHC -fno-warn-missing-methods  #-}
> {-# OPTIONS_GHC -fno-warn-orphans          #-}

> {-# LANGUAGE BangPatterns                  #-}
> {-# LANGUAGE TemplateHaskell               #-}
> {-# LANGUAGE QuasiQuotes                   #-}
> {-# LANGUAGE NoMonomorphismRestriction     #-}

> module Chap1 (
>     module Control.Applicative
>   , solveLaplaceStencil
>   , useBool
>   , boundMask
>   , boundValue
>   , bndFnEg1
>   , fivePoint
>   , ninePoint
>   , testStencil5
>   , testStencil9
>   , analyticValue
>   , slnHMat4
>   , slnHMat5
>   , testJacobi4
>   , testJacobi6
>   , bndFnEg3
>   , runSolver
>   , s5
>   , s9
>   ) where
>
> import Data.Array.Repa                   as R
> import Data.Array.Repa.Unsafe            as R
> import Data.Array.Repa.Stencil           as A
> import Data.Array.Repa.Stencil.Dim2      as A

> import Prelude                           as P

> import Data.Packed.Matrix
> import Numeric.LinearAlgebra.Algorithms

> import Chap1Aux

> import Control.Applicative

Laplace’s Equation: The Five Point Formula

We show how to apply finite difference methods to Laplace’s equation:

where

For a sufficiently smooth function (see (Iserles 2009, chap. 8)) we have

where the central difference operator is defined as

We are therefore led to consider the five point difference scheme.

We can re-write this explicitly as

Specifically for the grid point (2,1) in a grid we have

where blue indicates that the point is an interior point and red indicates that the point is a boundary point. For Dirichlet boundary conditions (which is all we consider in this post), the values at the boundary points are known.

We can write the entire set of equations for this grid as

A Very Simple Example

Let us take the boundary conditions to be

With our grid we can solve this exactly using the hmatrix package which has a binding to LAPACK.

First we create a matrix in hmatrix form

> simpleEgN :: Int
> simpleEgN = 4 - 1
>
> matHMat4 :: IO (Matrix Double)
> matHMat4 = do
>   matRepa <- computeP $ mkJacobiMat simpleEgN :: IO (Array U DIM2 Double)
>   return $ (simpleEgN - 1) >< (simpleEgN - 1) $ toList matRepa

ghci> matHMat4
  (2><2)
   [ -4.0, 1.0
   ,  1.0, 0.0 ]

Next we create the column vector as presribed by the boundary conditions

> bndFnEg1 :: Int -> Int -> (Int, Int) -> Double
> bndFnEg1 _ m (0, j) |           j > 0 && j < m = 1.0
> bndFnEg1 n m (i, j) | i == n && j > 0 && j < m = 2.0
> bndFnEg1 n _ (i, 0) |           i > 0 && i < n = 1.0
> bndFnEg1 n m (i, j) | j == m && i > 0 && i < n = 2.0
> bndFnEg1 _ _ _                                 = 0.0

> bnd1 :: Int -> [(Int, Int)] -> Double
> bnd1 n = negate .
>          sum .
>          P.map (bndFnEg1 n n)

> bndHMat4 :: Matrix Double
> bndHMat4 = ((simpleEgN - 1) * (simpleEgN - 1)) >< 1 $
>            mkJacobiBnd fromIntegral bnd1 3

ghci>  bndHMat4
  (4><1)
   [ -2.0
   , -3.0
   , -3.0
   , -4.0 ]

> slnHMat4 :: IO (Matrix Double)
> slnHMat4 = matHMat4 >>= return . flip linearSolve bndHMat4

ghci> slnHMat4
  (4><1)
   [               1.25
   ,                1.5
   , 1.4999999999999998
   , 1.7499999999999998 ]

The Jacobi Method

Inverting a matrix is expensive so instead we use the (possibly most) classical of all iterative methods, Jacobi iteration. Given and an estimated solution , we can generate an improved estimate . See (Iserles 2009, chap. 12) for the details on convergence and convergence rates.

The simple example above does not really give a clear picture of what happens in general during the update of the estimate. Here is a larger example

Sadly, WordPress does not seem to be able to render matrices written in LaTeX so you will have to look at the output from hmatrix in the larger example below. You can see that this matrix is sparse and has a very clear pattern.