III. Soving the diffusion problem with implicit Euler integration with multigrid & finite differences


In [254]:
using PyPlot

Problem description

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 fourth solution, we still use a second-order finite differences method in space. In time, however, we use an explicit singly diagonal implicit Runge-Kutta method. The details of this integration step are mentioned in the corresponding section below.

For the spatial derivatives, 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.

Types for spatial and temporal domain


In [255]:
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 [256]:
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

Generator for exact solution

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 [257]:
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[257]:
gen_exact (generic function with 1 method)

Matrix for finite differences

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 [258]:
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[258]:
fdmatrix (generic function with 2 methods)

Multigrid for spatial derivatives


In [259]:
# 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[259]:
level (generic function with 1 method)

In [260]:
function weighted_jacobi(A, x, b)
    weight = 2/3
    (speye(A) - weight*A) * x + weight*b
end


Out[260]:
weighted_jacobi (generic function with 1 method)

In [261]:
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[261]:
restrictionmatrix (generic function with 1 method)

In [262]:
function gen_mgsolver(matrix_generator, restriction_generator, smoother, grid::Grid2D, ν1::Int, ν2::Int, tol=1e-12, 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[262]:
gen_mgsolver (generic function with 3 methods)

Define ESDIRK method

We define an ESDIRK method with the following coefficients. As the implicit coefficient σ is the same for all stages, we always have the same matrix in the system $Ax=b$ that we solve – the only thing that changes is the right-hand side. To simplify the formulation, we use the following notation:

$$ FD = \text{finite differences matrix} \\ A = \Delta t \cdot D \cdot FD $$

This way, $\Delta t \cdot f(u) = A \cdot u$ and we get the following formulation for the intermediate stages $k_i = \Delta t \cdot f(\xi_i)$:

$$ k_1 = A \cdot u^n \\ k_2 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{21} \cdot k_1 \right) \\ k_3 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{31} \cdot k_1 + a_{32} \cdot k_2 \right) \\ k_4 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{41} \cdot k_1 + a_{42} \cdot k_2 + a_{43} \cdot k_3 \right) $$

Finally, we can compute the next time step of $u$:

$$ u^{n+1} = u^n + b_1 \cdot k_1 + b_2 \cdot k_2 + b_3 \cdot k_3 + b_4 \cdot k_4 $$

In [263]:
function gen_esdirk3(D::Real, grid::Grid2D, time::TimeSteps, matrix_generator::Function, solver::Function)
    
    # generate matrix for RHS: f(u) = Au
    A = time.dt * D * matrix_generator(grid)
    
    # coefficients for ESDIRK3
    c  = [0 1767732205903/2027836641118 3/5 1] # not used
    σ = 1767732205903/4055673282236
    a1 = [0 0 0 0]
    a2 = [σ σ 0 0]
    a3 = [2746238789719/10658868560708 -640167445237/6845629431997 σ 0]
    a4 = [1471266399579/7840856788654 -4482444167858/7529755066697 11266239266428/11593286722821 σ]
    b  = a4
    
    #=
    order_conditions = [
        sum(b) - 1
        b * c.' - 1/2
        b * (c.^2).' / 2 - 1/6
        b[1] * sum(a1 .* c) + b[2] * sum(a2 .* c) + b[3] * sum(a3 .* c) + b[4] * sum(a4 .* c) - 1/6
    ]
    println(order_conditions)
    =#
    
    function esdirk3(u)
        k1 = A * u
        k2 = A * solver(u + σ * k1)
        k3 = A * solver(u + a3[1] * k1 + a3[2] * k2)
        k4 = A * solver(u + a4[1] * k1 + a4[2] * k2 + a4[3] * k3)
        u += b[1] * k1 + b[2] * k2 + b[3] * k3 + b[4] * k4
    end
    
    return esdirk3
end


Out[263]:
gen_esdirk3 (generic function with 1 method)

Do integration


In [264]:
function integrate(x::Array, integration_step::Function, time::TimeSteps)
    for i in 1:time.Nt
        x = integration_step(x)
    end
    x
end


Out[264]:
integrate (generic function with 1 method)

Define and solve the problem


In [265]:
g = Grid2D(1,3,2,5,2^5,2^5,3)


Out[265]:
Grid2D(1,3,2,5,32,32,0.0625,0.09375,3)

In [266]:
t = TimeSteps(0, 1, 200)


Out[266]:
TimeSteps(0,1,200,0.005)

In [267]:
C0 = 5
D  = 0.1
exact = gen_exact(g, C0, D)


Out[267]:
exact (generic function with 1 method)

Exact solution


In [268]:
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))


Approximate solution


In [269]:
# 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 [270]:
σ = 1767732205903/4055673282236
solver = gen_mgsolver(g -> speye((g.Nx-1)*(g.Ny-1)) - t.dt * σ * D * fdmatrix(g),
                      restrictionmatrix, weighted_jacobi, g, 0, 1)
integration_step = gen_esdirk3(D, g, t, fdmatrix, solver)


Out[270]:
esdirk3 (generic function with 1 method)

In [271]:
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))


Compare to exact solution


In [272]:
norm(C - C_exact)


Out[272]:
0.016032013696206993