In [1]:
using PyPlot
For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.
\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}We therefore want to solve the following PDE:
$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$For this third solution, we again use backward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:
$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^{n+1} $$Here, A is the finite differences matrix for the 2D Laplace operator. However, as $u^{n+1}$ is not known in advance, we have to solve this equation for it. The scheme becomes the following, which has to be solved for $u^{n+1}$:
$$ \left( I - D \cdot \Delta t \cdot A \right) u^{n+1} = u^{n} $$In this case, however, we don’t solve for $u^{n+1}$ directly (e.g. with $LU$-decomposition). Instead, we use a multigrid approach to solve the system iteratively. As the spatial approximation with finite differences has a simple, regular form, we can easily create the finite differences matrix for other, coarser grids as well. The matrices have the exact same form as for the finest grid, except that they have fewer entries. However, we impose that the number of intervals is a power of two to ensure that they can easily be mapped to a coarser grid.
In [2]:
immutable Grid2D
x_start::Real
x_end::Real
y_start::Real
y_end::Real
Nx::Int
Ny::Int
dx::Real
dy::Real
levels::Int
function Grid2D(x_start, x_end, y_start, y_end, Nx, Ny, levels)
if Nx/2^(levels-1) != int(Nx/2^(levels-1))
error("ERROR: The number of intervals in x-direction is not divisible by 2^i")
end
if Ny/2^(levels-1) != int(Ny/2^(levels-1))
error("ERROR: The number of intervals in y-direction is not divisible by 2^i")
end
new(x_start, x_end, y_start, y_end, Nx, Ny, (x_end-x_start)/Nx, (y_end-y_start)/Ny, levels)
end
end
In [3]:
immutable TimeSteps
t_start::Real
t_end::Real
Nt::Int
dt::Real
TimeSteps(t_start, t_end, Nt) = new(t_start, t_end, Nt, (t_end-t_start)/Nt)
end
As a very simple scenario, we use the following initial conditions:
$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$This leads to the following analytical solution:
$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$
In [4]:
function gen_exact(grid, C0, D)
x_size = grid.x_end - grid.x_start
y_size = grid.y_end - grid.y_start
exact(x,y,t) = C0 * sin((x-grid.x_start)/x_size * pi) * sin((y-grid.y_start)/y_size * pi) *
exp( - (1/x_size^2 + 1/y_size^2) * pi^2 * D * t)
end
Out[4]:
We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).
Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:
\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.
\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.
$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$As we have to create several different matrices for the multigrid matrix, we define a function that returns those matrices depending on the number and size of intervals.
In [5]:
function fdmatrix(grid::Grid2D, matrix_function::Function = A -> A)
# number of unknowns in x and y direction
ux = (grid.Nx-1)
uy = (grid.Ny-1)
# finite difference matrix for ∂²u/∂x²
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1))
FDx = blkdiag({FDx_local for i in 1:uy}...)
# finite difference matrix for ∂²u/∂y²
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux))
# return combined matrix for Δu
matrix_function(1/grid.dx^2 * FDx + 1/grid.dy^2 * FDy)
end
Out[5]:
In [6]:
# helper function to get coarser levels of a grid
function level(grid::Grid2D, lvl)
if lvl > grid.levels
error("ERROR: Grid doesn't support level $(level).")
end
Grid2D(grid.x_start, grid.x_end, grid.y_start, grid.y_end, int(grid.Nx/2^(lvl-1)), int(grid.Ny/2^(lvl-1)), grid.levels-(lvl-1))
end
Out[6]:
In [7]:
function weighted_jacobi(A, x, b)
weight = 2/3
(speye(A) - weight*A) * x + weight*b
end
Out[7]:
In [8]:
function restrictionmatrix(grid::Grid2D)
ux_old = grid.Nx - 1
uy_old = grid.Ny - 1
ux_new = div(grid.Nx,2) - 1
uy_new = div(grid.Ny,2) - 1
R = spzeros(ux_new*uy_new, ux_old*uy_old)
for j = 1:uy_new
for i = 1:ux_new
i_new = ux_new*(j-1) + i
R[i_new, ux_old*(2*j-1) + 2*i ] = 4
R[i_new, ux_old*(2*j-1) + 2*i - 1] = 2
R[i_new, ux_old*(2*j-1) + 2*i + 1] = 2
R[i_new, ux_old*(2*j-1 - 1) + 2*i ] = 2
R[i_new, ux_old*(2*j-1 - 1) + 2*i - 1] = 1
R[i_new, ux_old*(2*j-1 - 1) + 2*i + 1] = 1
R[i_new, ux_old*(2*j-1 + 1) + 2*i ] = 2
R[i_new, ux_old*(2*j-1 + 1) + 2*i - 1] = 1
R[i_new, ux_old*(2*j-1 + 1) + 2*i + 1] = 1
end
end
R /= 16
end
Out[8]:
In [9]:
function gen_mgsolver(matrix_generator, restriction_generator, smoother, grid::Grid2D, ν1::Int, ν2::Int, tol=1e-6, max_it = 100)
# create main matrices for all levels
A = Dict()
for i = 1:grid.levels
A[i] = matrix_generator(level(grid,i))
end
# create restriction matrices for all levels
R = Dict()
for i = 1:grid.levels-1
R[i] = restriction_generator(level(grid,i))
end
function mgsolver(b::Array, x0::Array = b, level::Int = 1)
if level == grid.levels
return A[level] \ b
end
x = x0
for it = 1:max_it
# smoothen problem on current level
for i = 1:ν1
x = smoother(A[level], x, b)
end
# restrict residual to next level and solve recursively
r = R[level] * (b - A[level] * x)
e = mgsolver(r, zeros(r), level+1)
# prolongate error to current level and add to solution
x += 4*R[level].' * e
# post-smoothening on current level
for i = 1:ν2
x = smoother(A[level], x, b)
end
if level > 1
return x
end
if norm(A[level] * x - b) < tol
#println("Converged with $(it) iterations")
return x
end
end
return x
end
return mgsolver
end
Out[9]:
In [10]:
function integrate(x::Array, integration_step::Function, time::TimeSteps)
for i in 1:time.Nt
x = integration_step(x)
end
x
end
Out[10]:
In [11]:
g = Grid2D(1,3,2,5,2^6,2^6,3)
Out[11]:
In [12]:
t = TimeSteps(0, 20, 100)
Out[12]:
In [13]:
C0 = 5
D = 0.001
exact = gen_exact(g, C0, D)
Out[13]:
In [14]:
C_exact = Float64[exact(x,y,t.t_end) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end]
pcolormesh(C_exact)
colorbar()
clim((0,5))
In [15]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,t.t_start) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end];
In [16]:
integration_step = gen_mgsolver(g -> speye((g.Nx-1)*(g.Ny-1)) - t.dt * D * fdmatrix(g),
restrictionmatrix, weighted_jacobi, g, 0, 1)
Out[16]:
In [17]:
x = vcat(C[2:end-1,2:end-1]...)
x = integrate(x, integration_step, t)
C[2:end-1,2:end-1] = reshape(x,g.Nx-1,g.Ny-1)
pcolormesh(C)
colorbar()
clim((0,C0))
In [18]:
vecnorm(C - C_exact)
Out[18]: