Multigrid algorithm

Introduction. Solving linear systems

Large, sparse linear systems of equations Ax = b arise in many applications e.g. Many approaches can be used to find x on GPU A few Jacoby iterations are enough to get good looking visual effect in 2D fluid simulations but long wave relaxation is very slow on fine grid. Therefore we need multigrid algorithm for accurate simulation of heart tissue contraction. Yang Chenglin sent me kindly link on his FFT-based Droplet Simulation. Conjugate gradients solver will be very interesting too.

1D weighted Jacobi relaxation algorithm

Consider one-dimensional boundary value problem
    u" = f(x),     0 < x < 1,     u(0) = u(1) = 0,
with the finite difference we get
    ui - 1 - 2ui + ui + 1 = h2fi ,     u0 = un = 0,     h = 1/n,     xi = i h,     i = 0, 1, ... , n.
We can solve these equations with respect to ui for Jacobi relaxation [1]
    uinew = (ui-1old + ui+1old - h2fi )/2.
Jacobi iterations converge to solution of the problem. If we use weighted value uiold + ω(uinew - uiold) we get weighted Jacobi relaxation (we put fi = 0 further)
    uinew = (1 - ω)ui + ω(ui-1 + ui+1 )/2,     (*)
Press the Step button below to see relaxation process for fi = 0,   n = 64 and initial condition
    ui = [c1 sin(πi/n) + c3 sin(3πi/n) + c5 sin(5πi/n)]/3.
it c1 c3 c5 ω Ymin Ymax

Error log( max |ui | ) is plotted against iteration number (the red curve).

Note that at first fast oscillations decay quickly, then the main harmonic decays slowly at constant rate.
It is not difficult to check by direct substitution that
    ui(m) = sin πmi/n = sin km i
is eigenvector of (*) with the eigenvalue
    λm = (1 - ω) + ω cos km = 1 - 2ω sin2πm/2n.
λm determines relaxation rate for weighted Jacobi iterations. For ω = 1/3, 1/2, 2/3, 1 they are plotted below (-1 ≤ λm ≤ 1 and 0 ≤ m ≤ n, see also [1]). For ω = 2/3 (the red curve) and n/2 < m < n we get m| ≤ 1/3 therefore high harmonics decay at least 3 times after every iteration.

For small m/n we get approximately
    λm = 1 - ω(πm/n)2/2,
so relaxation is very slow for m/n < 1/100 values. As since convergence is faster for small n we need caurse grids to accelerate relaxation. We consider 2D Poisson equation next.

[1] William L. Briggs, Van Emden Henson, Steve F. McCormick   A Multigrid Tutorial (look at "Tutorial Slides")


Simulations on GPU     updated 6 Aug 2012