Commit eaa40ac1 authored by Jochen Schulz's avatar Jochen Schulz
Browse files

added explanation for stencils (from code)

parent abe31bb1
%% Cell type:markdown id: tags:
# Poisson equation in 2d
## Poisson equation
Find $u \in C^2(\Omega)\cap C(\overline{\Omega})$ with
$$
\begin{cases}
-\Delta u = f & \text{in } \Omega \\
u = 0 & \mbox{on } \partial \Omega
\end{cases}
$$
for $\Omega = (0,1)^2$ and $f \in C(\Omega)$. For $f = 0$, the equation is known as *Laplace equation*.
## Laplace operator
$$
\Delta u := \sum_{i=1}^d \frac{\partial ^2 u}{\partial x_i^2}
$$
### Discretization
- Equidistant grid $\Omega_N = I_N \times I_N$, where
$$I_N := \left\{ x^N_j \colon j = 0, \ldots, N-1 \right\} \subset (0, 1), \qquad x^N_j := \left( j + \frac{1}{2} \right) h_N, \qquad h_N := \frac{1}{N}.$$
- Approximation of $\frac{\partial^2 u}{\partial x^2}$:
$$\frac{u(x-h, y) - 2 u(x, y) + u(x+h, y)}{h^2} = \frac{\partial^2 u}{\partial x^2} (x, y) + \mathcal{O}(h^2),$$
and similarly for $\frac{\partial^2 u}{\partial y^2}$.
- Approximation of Laplace operator:
$$\frac{u(x, y-h) + u(x-h, y) - 4 u(x, y) + u(x, y+h) + u(x+h, y)}{h^2} = \Delta u(x, y) + \mathcal{O}(h^2)$$
- Define
$$u^N_{jk} := u(x^N_j, x^N_k), \qquad f^N_{jk} := f(x^N_j, x^N_k),$$
for $j, k = 0, 1, \ldots, N-1$, and
$$(\Delta^N u^N)_{jk} := \frac{u^N_{j,k-1} + u^N_{j-1,k} - 4 u^N_{jk} + u^N_{j,k+1} + u^N_{j+1,k}}{h_N^2}.$$
Here, $u^N_{-1,k} = u^N_{N,k} = u^N_{j,-1} = u^N_{j,N} = 0$ (boundary condition).
- Solve
$$-\Delta^N u^N = f^N.$$
## Code
For performance reasons, the code is not written in the notebook. See files `01-hmatrix.hs`, `02-repa.hs` and `03-stencils.hs`.
%% Cell type:code id: tags:
``` haskell
negLaplace :: Array U DIM2 Double -> Array D DIM2 Double
negLaplace arr = fromFunction ext computeElem
where
!ext = extent arr
get !idx | inShape ext idx = unsafeIndex arr idx
| otherwise = 0
computeElem (Z:.(!jx):.(!jy)) = ( 4*get (Z:.jx:.jy)
- get (Z:.(jx-1):.jy) - get (Z:.(jx+1):.jy)
- get (Z:.jx:.(jy-1)) - get (Z:.jx:.(jy+1)) )
/ h**2
```
%% Cell type:markdown id: tags:
### Stencils
Operations like `laplace` in 02-repa.hs are a frequent pattern: each element
of the resulting array is a linear combination of the surrounding elements of
the input array.
Repa has special means for constructing such operations in 2-d using stencils. Stencils are constructed with
``` haskell
makeStencil2 :: Num a => Int -> Int -> (DIM2 -> Maybe a) -> Stencil DIM2 a
```
and applied to 2-d arrays with
``` haskell
mapStencil2 :: Source r a => Boundary a -> Stencil DIM2 a -> Array r DIM2 a -> Array PC5 DIM2 a
```
`Boundary a` is a type that handles the boundary conditions. It can be
- `BoundFixed x` returns a fixed value `x :: a` at the boundary,
- `BoundConst x` assumes pixels outside the array have the value `x :: a`,
- `BoundClamp` extends the boundary values to outside the boundary.
`PC5` is a special representation for stencil results; it is defined as (no
need to understand this)
``` haskell
type PC5 = P C (P (S D) (P (S D) (P (S D) (P (S D) X))))
```
`PC5` arrays can be converted to `D` arrays using the polymorphic function
``` haskell
delay :: (Shape sh, Source r e) => Array r sh e -> Array D sh e
```
%% Cell type:code id: tags:
``` haskell
negLaplaceStencil = makeStencil2 3 3 computeElem :: Stencil DIM2 Double
where
computeElem (Z:. -1 :. 0) = Just $ -c
computeElem (Z:. 1 :. 0) = Just $ -c
computeElem (Z:. 0 :. -1) = Just $ -c
computeElem (Z:. 0 :. 1) = Just $ -c
computeElem (Z:. 0 :. 0) = Just $ 4*c
computeElem _ = Nothing
c = 1/h**2
```
%% Cell type:markdown id: tags:
There is special syntax for stencils. The laplaceStencil above can be
produced more easily by
``` haskell
laplaceStencil = [stencil2| 0 -1 0
-1 4 -1
0 -1 0 |]
```
```haskell
laplace = R.map (/h**2) . mapStencil2 (BoundConst 0) laplaceStencil
```
This syntax requires the QuasiQuotes language extension. The stencil is
constructed at compile time and can therefore not depend on runtime values.
That's why we have to use R.map to do the division by h**2 here.
%% Cell type:code id: tags:
``` haskell
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment