Multigrid Methods

Revant Kumar

► To cite this version:

Revant Kumar. Multigrid Methods. Engineering school. High Performance Computing for Engineering Problems, 11th Indo-German Winter Academy 2012, India. 2012, pp.41. [hal-04216283](https://hal.science/hal-04216283v1)

HAL Id: hal-04216283

https://hal.science/hal-04216283v1

Submitted on 24 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
MULTIGRID METHODS

REVANT KUMAR

MATHEMATICS AND COMPUTING
INDIAN INSTITUTE OF TECHNOLOGY, GUWAHATI

SUPERVISORS: PROF. ULRICH RÜDE, CHRISTIAN KUSCHEL
UNIVERSITY OF ERLANGEN - NUREMBERG (FAU), GERMANY
OUTLINE

□ Introduction

□ Basic Concepts

□ 1D Model Problem

□ Multigrid Algorithm

□ 2D Model Problem

□ Computational Complexity

□ Summary

□ References
# Introduction

## Multigrid Methods:

* Efficient and versatile set of algorithms for solving many types of computational problems.

* are a group of algorithms for solving differential and integral equations using a hierarchy of discretization.

* Applications:

* Mathematical Applications:
    * Differential Equations and Systems
    * Integral Equations
    * Elliptic Boundary Value Problems

* Physical Applications:
    * Optimization
    * Image Reconstruction
    * Computational Fluid Dynamics
# BASIC CONCEPTS

## Jacobi's player alignment problem

* $N + 1$ players are numbered from 0 to $N$.

* Players numbered 0 and $N$ stand at the left and right goal posts respectively, which are at a distance $L$ apart.

* Players numbered $i = 1, ..., N-1$ are to stand between the players numbered 0 and $N$ equally spaced.
# BASIC CONCEPTS

## Jacobi's player alignment problem : Global process

* Players numbered 0 and N stand at the left and right goal posts respectively.

* Player i, $i = 1, ..., N - 1$ could be ordered to move to the point on the goal line that is at a distance $iL/N$ from the left post.

* This is a global process which solves the problem directly.

<figure><img src="image_2.png" /></figure>
# BASIC CONCEPTS

## Jacobi's player alignment problem : Global process

* Players numbered 0 and N stand at the left and right goal posts respectively.

* Player i, $i = 1, ..., N - 1$ could be ordered to move to the point on the goal line that is at a distance $iL/N$ from the left post.

* This is a global process which solves the problem directly.

<figure><img src="image_3.png" /></figure>
# BASIC CONCEPTS

## Jacobi's player alignment problem : Global process

* Players numbered 0 and N stand at the left and right goal posts respectively.

* Player i, $i = 1, ..., N-1$ could be ordered to move to the point on the goal line that is at a distance $iL/N$ from the left post.

* This is a global process which solves the problem directly.

<figure><img src="image_4.png" /></figure>
# BASIC CONCEPTS

## Jacobi's player alignment problem : Global process

* Players numbered 0 and N stand at the left and right goal posts respectively.

* Player i, $i = 1, ..., N-1$ could be ordered to move to the point on the goal line that is at a distance $iL/N$ from the left post.

* This is a global process which solves the problem directly.

<figure><img src="image_5.png" /></figure>
# BASIC CONCEPTS

## Jacobi's player alignment problem : Local process

* To achieve equal spacing and ordering as in global process, a local approach can also be taken.

* Local approach taken is that each player is to stand between his two neighbours.

* At each iteration, after whistle is blown, player $i$, $i = 1, ..., N-1$, moves to the point at the centre between positions of player $i-1$ and player $i+1$ just before the whistle was blown.

* The error reduces in each iteration.
## Jacobi's player alignment problem

<figure><img src="image_6.png" /><figcaption>Figure 1. Jacobi's player-alignment iteration. Red disks show the current position, and blue disks show the previous position, before the last whistle blow: (a) initial position, (b) one whistle blow, (c) slow convergence, (d) fast convergence, (e) overshoot, and (f) damped iteration.</figcaption></figure>
# BASIC CONCEPTS

## Multiscale Concept

* The main idea behind multigrid methods is to use simple local process but to do so at all scales of the problem

<figure><img src="image_7.png" /><figcaption>Figure 2. Multiscale player alignment. Red disks show the current position, and blue disks show the previous position, before the last whistle blow: (a) large, (b) medium, and (c) small scales.</figcaption></figure>
# BASIC CONCEPTS

## Multiscale Concept

* In general, developing a multiscale solver for a given problem involves four main tasks:

* Choosing an appropriate local process

* Choosing appropriate coarse (large-scale) variables

* Choosing appropriate methods for transferring information across scales

* Developing appropriate equations (or processes) for the coarse variables

* Depending on the application, each of these tasks can be simple, moderately difficult, or extremely challenging.
# BASIC CONCEPTS

## Analysing Jacobi's Player-Alignment Problem

* Initial Error: $e^{(0)} = [e_1^{(0)}, ..., e_{N-1}^{(0)}]^T$, $e_0^{(0)} = e_N^{(0)} = 0$

* Error in the location of player $i$ after iteration $n$:

$$e_i^{(n)} = \frac{1}{2} (e_{i+1}^{(n-1)} + e_{i-1}^{(n-1)}), \quad i = 1, ..., N-1$$

$$e^{(n)} = Re^{(n-1)}$$

where $R$ is a $N-1 \times N-1$ tri-diagonal matrix

* Error $e$ is an eigenvector of $R$, satisfying $Re = \lambda e$ for some (scalar) eigenvalue $\lambda$

$$\lambda^{(k)} = \cos\left(\frac{k\pi}{N}\right), \quad \text{for } k = 1, ..., N-1$$
# BASIC CONCEPTS

## Analysing Jacobi's Player-Alignment Problem

* The $N - 1$ corresponding eigenvectors are given by:

$$v^{(k)} = \begin{bmatrix} \sin\left(\frac{k\pi}{N}\right) \\ \sin\left(\frac{2k\pi}{N}\right) \\ \vdots \\ \sin\left(\frac{(N-1)k\pi}{N}\right) \end{bmatrix}$$

<figure><img src="image_8.png" /><figcaption>Figure 3. Three eigenvectors. We associate scale with the wavelength, which is inversely proportional to $k$: (a) $v^{(1)}$, (b) $v^{(3)}$, and (c) $v^{(6)}$.</figcaption></figure>
# BASIC CONCEPTS

## Analysing Jacobi's Player-Alignment Problem

* The eigenvectors $v^{(k)}$ are linearly independent, spanning the space $\mathbf{R}^{N-1}$

$$ \therefore e^{(0)} = \sum_{k=1}^{N-1} c_k v^{(k)}, \quad \text{where } c_k \in \mathbb{R} $$

$$ e^{(n)} = \mathbf{R}^n e^{(0)} = \sum_{k=1}^{N-1} c_k \left[ \cos\left(\frac{k\pi}{N}\right) \right]^n v^{(k)} $$

$$ |\lambda^{(k)}| = |\cos(k\pi/N)| < 1, \quad \text{for } 1 \le k \le N-1 $$

$$ \therefore |\lambda^{(k)}|^n \to 0, \quad \text{as } n \to \infty $$
# BASIC CONCEPTS

## Damped Jacobi Iteration's Smoothing Factor

* To overcome the overshoot problem, we introduce a damping factor. If damping factor is 2/3, then

$$e_i^{(n)} = \frac{1}{3}e_i^{(n-1)} + \frac{2}{3}\left[\frac{1}{2}\left(e_{i+1}^{(n-1)} + e_{i-1}^{(n-1)}\right)\right], \quad i = 1, \dots, N-1$$

$$e^{(n)} = R_d e^{(n-1)} = \left(\frac{1}{3}I + \frac{2}{3}R\right) e^{(n-1)}$$

$R_d$ — Damped Jacobi Iteration matrix, $I$ — Identity Matrix

* The eigenvectors remain the same, but the eigenvalues of $R_d$ are given by

$$\lambda_d^{(k)} = \frac{1}{3} + \frac{2}{3}\lambda^{(k)} = \frac{1}{3} + \frac{2}{3}\cos\left(\frac{k\pi}{N}\right)$$

* These damped Jacobi's Iterations are efficient in reducing small-scale errors.
# BASIC CONCEPTS

## Damped Jacobi Iteration's Smoothing Factor

* Declare all eigenvectors with $1 \le k < N/2$ smooth (large scale)

* Declare all eigenvectors with $N/2 \le k < N$ rough (small scale)

$$ |\lambda_d^{(k)}| \le \frac{1}{3}, \quad \text{for } N/2 \le k < N $$

$$ \therefore \text{Smoothing Factor} = \frac{1}{3} $$

* Thus, damped Jacobi iteration (called relaxation) is efficient at smoothing the error, but is inefficient at reducing the error once it has been smoothed.

<figure><img src="image_9.png" /><figcaption>Figure 4. The smoothing property of the damped Jacobi relaxation. (a) Initial position, (b) after four relaxations, and (c) after eight relaxations.</figcaption></figure>
# 1D MODEL PROBLEM

## One-Dimensional (1D) Boundary Value Problem

$$Lu = \frac{d^2 u}{dx^2} = f(x), \quad x \in (0,1)$$

$$u(0) = u_l, u(1) = u_r$$

where $f$ is a given function, $u_l$ and $u_r$ are the boundary conditions, and $u$ is the solution we are seeking.

**GRID:** A uniform mesh is defined on the interval $[0,1]$, with mesh size $h = 1/N$ and $N + 1$ mesh points given by $x_i = ih$ for $i = 0, 1, ..., N$.

<figure><img src="image_10.png" /></figure>

$x_0 = 0$

$x_i = ih$

$x_N = 1$

□ Discrete approximation to $u$ and $f$:

$$u_i^h = u(x_i), f_i^h = f(x_i) \text{ for } i = 0, 1, ..., N$$

$$u_0^h = u_l, u_N^h = u_r$$
# 1D MODEL PROBLEM

## One-Dimensional (1D) Boundary Value Problem

* Using Taylor series differential equation can be written as

$$ (L^h u^h)_i \equiv \frac{u_{i-1}^h - 2u_i^h + u_{i+1}^h}{h^2} = f_i^h, \quad i = 1, \dots, N-1 $$

* Above system can also be written more concisely as

$$ L^h u^h = f^h $$

* Jacobi's algorithm can be easily adopted to solve the above equation.

* Let $u^{h(0)}$ denote an initial approximation to $u^h$ satisfying the boundary conditions.

* At the $n^{th}$ iteration, we set $u_i^{h(n)}, i = 1, \dots, N-1$, as follows

$$ u_i^{h(n)} = \frac{1}{2} (u_{i+1}^{h(n-1)} + u_{i-1}^{h(n-1)} - h^2 f_i^h), \quad i = 1, \dots, N-1 $$
# 1D MODEL PROBLEM

## One-Dimensional (1D) Boundary Value Problem

* The error in the approximation after $n^{th}$ iteration:

$$e^{h(n)} = u^h - u^{h(n)}$$

assuming $e_0^{h(n)} = e_N^{h(n)} = 0$,

since all the approximations satisfy the boundary conditions.

* Using the iteration equations in $u^h$, we get,

$$e_i^{h(n)} = \frac{1}{2} (e_{i+1}^{h(n-1)} + e_{i-1}^{h(n-1)}), \quad i = 1, \dots, N-1$$

* In the context of 1D boundary Value Problem, damped Jacobi relaxation, with a damping factor of 2/3, is defined by

$$u_i^{h(n)} = \frac{1}{3}(u_i^{h(n)}) + \frac{2}{3}\left[\frac{1}{2}(u_{i+1}^{h(n-1)} + u_{i-1}^{h(n-1)} - h^2 f_i^h)\right]$$

* Its error propagation properties are identical to those of the player-alignment problem.
# 1D MODEL PROBLEM

## Naïve Multiscale Algorithm

* The main point of the exercise is to demonstrate that to know what to do at the large scales, we must first propagate appropriate information from the small scales to the large scales.

* Assuming $N$ is a power of 2, perform Jacobi's iteration for the problem defined at the largest scale, which consists of the variables $[u_0^h, u_{N/2}^h, u_N^h]$.

$$u_{N/2}^h = \frac{1}{2}(u_0^h + u_N^h - 0.25 f_{N/2}^h)$$

* Next, based on the values of the largest scale, we perform a Jacobi iteration at the second largest scale, obtaining $u_{N/4}^h$ and $u_{3N/4}^h$.

* We continue similarly until all the variables have been determined. However, the algorithm is computationally inefficient.

* Unfortunately, the algorithm fails to give the correct answer (unless $f$ is trivial) because we are not solving the correct equations on the coarse grids.
# 1D MODEL PROBLEM

## Aliasing

* The coarse grid affords roughly half as many eigenvectors as the fine grid, so some pairs of fine grid eigenvectors must alias.

* Smooth eigenvectors are those that can be represented accurately on the coarse grid, whereas rough eigenvectors alias with smooth ones.

* Thus, constructing correct equations for coarse grids requires information from the fine grid.

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 5. Aliasing. Red disks correspond to the coarse grid: (a) $v^{(3)}$ and (b) $-v^{(13)}$.</figcaption></figure>
# 1D MODEL PROBLEM

The approach taken is:

* First apply damped Jacobi relaxation to smooth the algebraic error.

* Construct a linear system for this error vector and approximate it on a coarser grid.

* The above process continues recursively, such that, on each coarser grid, approximate equations for the next finer grid's algebraic error are solved.

Algebraic error is used as coarse-grid variable because:

* Limitation in computing an accurate approximation to $u^h$ on the coarse grid due to aliasing.

* Relaxation has smoothed the error, and thus it can be approximated accurately on the coarse grid.
# 1D MODEL PROBLEM

* Let $\tilde{u}^h$ denote our current approximation to $u^h$ after performing several damped Jacobi relaxations.

* Let $e^h = u^h - \tilde{u}^h$ denote the corresponding algebraic error, which is known to be smooth.

* The system of algebraic error:

$$ (L^h e^h)_i = f_i^h - (L^h \tilde{u}^h)_i \equiv r_i^h, \quad i = 1, \dots, N-1 $$

$$ e_0^h = 0, \quad e_N^h = 0 $$

where $r_i^h$ is the $i^{th}$ residual

* This system $(L^h e^h = r^h)$ is to be approximated on the coarser grid

* Once a good approximation to $e^h$ is computed, then it is added to $\tilde{u}^h$, and thus improve the approximation to $u^h$.
# MULTIGRID ALGORITHM

* To transfer information between grids, suitable operators are needed.

* Considering just two grids :

* *Fine Grid:* mesh size $h$ with $N-1$ mesh points excluding boundary points.

* *Coarse Grid:* mesh size $2h$ with $N/2-1$ mesh points.

* Later, we require two parameters, $v_1$ and $v_2$, that denote, respectively, the number of damped Jacobi relaxations performed before and after the coarse grid correction
# MULTIGRID ALGORITHM

□ **Restriction operator**, denoted by $I_h^{2h}$, for fine to coarse data transfer.

For a fine grid vector $g^h$,

By Injection:

$$ (I_h^{2h} g^h)_i = g_{2i}^h, \quad i = 1, \dots, \frac{N}{2} - 1 $$

<figure><img src="image_13.png" /></figure>
# MULTIGRID ALGORITHM

□ **Restriction operator**, denoted by $I_h^{2h}$, for fine to coarse data transfer.

For a fine grid vector $g^h$,

By Full Weighting (Apply Local Averaging):

$$ (I_h^{2h} g^h)_i = \frac{1}{4} (g_{2i-1}^h + 2g_{2i}^h + g_{2i+1}^h), \quad i = 1, \dots, \frac{N}{2} - 1 $$

<figure><img src="image_14.png" /></figure>
# MULTIGRID ALGORITHM

□ **Prolongation operator**, denoted by $I_{2h}^h$, for coarse to fine data transfer.

For a coarse grid vector $g^{2h}$,

Linear Interpolation:

$$ (I_{2h}^h g^{2h})_i = \begin{cases} g_{i/2}^{2h}, & \text{if } i \text{ is even} \\ \frac{1}{2}(g_{(i+1)/2}^{2h} + g_{(i-1)/2}^{2h}), & \text{if } i \text{ is odd} \end{cases} \quad i = 1, \dots, N-1 $$

<figure><img src="image_15.png" /></figure>
# MULTIGRID ALGORITHM

## Algorithm V-Cycle

$u^h = VCycle(u^h, h, N/2, f^h, v_1, v_2)$

1. If $N == 2$, solve $L^h u^h = f^h$ and return $u^h$;

2. $u^h = Relax(u^h, h, N, f^h, v_1)$; // relaxing $v_1$ times

3. $r^h = f^h - L^h u^h$; // Computing Residual Vector

4. $f^{2h} = I_h^{2h} r^h$; // Restrict Residual

5. $u^{2h} = 0$; // Setting $u^{2h}$ to 0, including boundary conditions

6. $u^{2h} = VCycle(u^{2h}, 2h, N/2, f^{2h}, v_1, v_2)$; //treat coarse-grid problem recursively

7. $u^h = u^h + I_{2h}^h u^{2h}$; // interpolate and add correction

8. $u^h = Relax(u^h, h, N, f^h, v_2)$; // relaxing $v_2$ times

9. Return $u^h$
MULTIGRID ALGORITHM

Coarse – Grid Correction

<figure><img src="image_16.png" /></figure>
# MULTIGRID ALGORITHM

## SCHEMATIC DESCRIPTION OF V-CYCLE

<figure><img src="image_17.png" /></figure>
# 2D MODEL PROBLEM

## Two-Dimensional (2D) Boundary Value Problem

$$Lu = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x, y), \quad (x, y) \in \Omega$$

$$u = b(x, y), \quad (x, y) \in \partial\Omega$$

where $f$ and $b$ are given functions and $u$ is the solution we're seeking.

* For simplicity, our domain is the unit square, $\Omega = (0,1)^2$, with $\partial\Omega$ denoting the boundary of $\Omega$.

* This is the well-known Poisson Problem with Dirichlet boundary conditions.

* Uniform grid: $h = 1/N$, $(x_i, y_j) = (ih, jh)$, $i, j = 0, 1, ..., N$

* Our vector of unknowns are $u^h, b^h$ and $f^h$
# 2D MODEL PROBLEM

## Two-Dimensional (2D) Boundary Value Problem

* By Taylor Series,

$$ (L^h u^h)_{i,j} = \frac{u_{i-1,j}^h + u_{i+1,j}^h + u_{i,j-1}^h + u_{i,j+1}^h - 4u_{i,j}^h}{h^2} = f_{i,j}^h, \quad i,j = 1, \dots, N-1 $$

$$ u_{i,j}^h = b_{i,j}^h, \quad i = 0 \text{ or } i = N \text{ or } j = 0 \text{ or } j = N $$

* The above equation can be written more concisely as $L^h u^h = f^h$

* Jacobi's Algorithm: At each iteration, we change our approximation to $u_{i,j}^h$, $i,j = 1, \dots, N-1$, to satisfy the $(i,j)^{th}$ equation, with neighbouring values taken from previous iteration.

* Damped Jacobi relaxation can be adapted in the case of 2D boundary value problem. It has been found that the best achievable smoothing factor is 0.6, obtained with a damping factor of 0.8.

* To apply the multigrid V-cycle to the 2D model problem, we coarsen naturally by eliminating all the mesh-points, with either *i* or *j* odd, such that only about ¼ of the points remain.
2D MODEL PROBLEM

For restriction, for a fine-grid vector, $g^h$,

*   **Simple Injection:** $(I_h^{2h} g^h)_{i,j} = g_{2i,2j}^h$, $i,j = 1, ..., N/2 - 1$

*   **Full weighting : Local Averaging applied**

$$ (I_h^{2h} g^h)_{i,j} = \frac{1}{16} \left( \begin{array}{l} g_{2i-1,2j-1}^h + g_{2i-1,2j+1}^h + g_{2i+1,2j-1}^h + g_{2i+1,2j+1}^h + \\ \qquad 2(g_{2i-1,2j}^h + g_{2i+1,2j}^h + g_{2i,2j-1}^h + g_{2i,2j+1}^h) + \\ \qquad 4g_{2i,2j}^h \end{array} \right), $$

where $i, j = 1, ..., N/2 - 1$

<figure><img src="image_18.png" /></figure>

We denote the operator by using a "give to" stencil, [ ]. Centered over a c-point, ●, it shows what fractions of the neighboring (●) f-points' value is contributed to the value at the c-point.

<figure><img src="image_19.png" /></figure>
2D MODEL PROBLEM

□ For prolongation, for a coarse-grid vector, $g^{2h}$,
Bilinear Interpolation:

$(I_{2h}^h g^{2h})_{i,j}$

$$
= \begin{cases}
g_{i/2, j/2}^{2h}, & \text{if i and j both are even} \\
\frac{1}{2} \left( g_{(i+1)/2, j/2}^{2h} + g_{(i-1)/2, j/2}^{2h} \right), & \text{if i is odd and j is even} \\
\frac{1}{2} \left( g_{i/2, (j+1)/2}^{2h} + g_{i/2, (j-1)/2}^{2h} \right), & \text{if i is even and j is odd} \\
\frac{1}{4} \left( g_{(i+1)/2, (j+1)/2}^{2h} + g_{(i-1)/2, (j+1)/2}^{2h} + g_{(i+1)/2, (j-1)/2}^{2h} + g_{(i-1)/2, (j-1)/2}^{2h} \right), & \text{if i is odd and j is odd}
\end{cases}
$$

<figure><img src="image_20.png" /></figure>

We denote the operator by using a "give to" stencil, []. Centered over a c-point, ●, it shows what fraction of the c-point's value is contributed to neighboring f-points, ●.

<figure><img src="image_21.png" /></figure>
COMPUTATIONAL COMPLEXITY

□ Number of Operations performed in the V-Cycle on each grid is proportional to the number of mesh-points.

□ Number of variables at each grid is approximately a constant fraction of that of the next finer grid (1/2 in 1D and 1/4 in 2D)

□ Hence the total number of operations per V-Cycle is greater than the number of operations performed on just the finest grid only by a constant factor.

□ In 1D, this factor is bounded by $1 + 1/2 + 1/4 + \cdots = 2$

□ In 2D, this factor is bounded by $1 + 1/4 + 1/16 + \cdots = 4/3$
# COMPUTATIONAL COMPLEXITY

□ **WORK UNIT (WU):** The amount of computation required for a single Jacobi relaxation on the finest grid.

□ V(2,1) cycle on the finest grid requires:

- 3 WUs for relaxation

- 1 WU for computing the residual

- The prolongation plus the restriction operations together require no more than 1 WU

□ The fine-grid work is roughly $3 + 1 + 1 = 5$ WUs per V-Cycle

□ Entire Work is approximately:

- $5 \times 2 = 10$ WUs per V-Cycle in 1D

- $5 \times 4/3 = 7$ WUs per V-Cycle in 2D

□ Finally, on each grid in the hierarchy, the operations (relaxation, prolongation, restriction) can be performed in parallel.
SUMMARY

☐ Multigrid has been proven to be efficient on a wide variety of problems, especially Elliptic PDEs and Integral Equations.

☐ With the right setup, multigrid is frequently an optimal solver.

☐ Multigrid algorithms are of great interest because they are one of the few set of scalable algorithms, and can be parallelized readily and efficiently.

☐ The overall complexity of the multigrid V cycles is $O(N)$.
# SUMMARY

## Drawbacks of Multigrid

* The fixed costs for setting up the problem, etc., are higher than for other methods. Therefore, Multigrid Method often pays off only for large problems.

* It's necessary to find coarse grids for the fine grid (often not straightforward).
REFERENCES

Yavneh L., "Why MultigridMethods are so Efficient". Computing in Science and Engineering, vol.8, no. 6, Dec'06.

Briggs W.L., Henson W.L., McCormick S.F., "A Multigrid Tutorial"
THANK YOU