================ Multigrid solver [TOC] ================ For the sake of simplicity, we will only consider one-dimensional boundary problems of the form ---- LATEX --------------------------------------------------------------------- \begin{cases} -u''(x) = f(x), & x \in (0,1), \\ u(x) = g(x), & x \in \{0,1\}. \\ \end{cases} -------------------------------------------------------------------------------- This is called a one-dimensional Poisson equation with Dirichlet boundary conditions. For the discretization we will use the finite difference method (FDM) which leads to solving ---- LATEX --------------------------------------------------------------------- A_h \widehat{u}_h = \widehat{q}_h,\quad h = \frac{1}{N+1} -------------------------------------------------------------------------------- with ---- LATEX --------------------------------------------------------------------- A_h = \frac{1}{h^2} \left(\begin{array}{ccccc} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots& \ddots& \ddots& \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array}\right) \in \mathbb{R}^{N \times N} -------------------------------------------------------------------------------- and ---- LATEX --------------------------------------------------------------------- \widehat{q}_h = \left(\begin{array}{l} f_1 + \frac{1}{h^2} g(0) \\ f_2 \\ \vdots \\ f_{N-1} \\ f_N + \frac{1}{h^2} g(1) \\ \end{array}\right) \in \mathbb{R}^{N} -------------------------------------------------------------------------------- The numerical solution $\widehat{u}_h$ then will approximate the solution of the original problem as follows ---- LATEX --------------------------------------------------------------------- \widehat{u}_h \approx \left(\begin{array}{c} u(x_1) \\ \vdots \\ u(x_N) \\ \end{array}\right), \quad x_i = i \cdot h \;(1 \leq i \leq N) -------------------------------------------------------------------------------- Numerical representation of the grid functions ============================================== The so called *grid functions* $\widehat{u}_h$ and $\widehat{q}_h$ are $N$ dimensional vectors that give approximations of their continuous counterparts at the interior grid points $x_1, \dots, x_N$. However, for our numerical algorithms we will use $N+2$ dimensional vectors that also take the boundary node into account. We indicate this difference by skipping the _hat_ in the notation. For example, $\widehat{u}_h$ is represented as ---- LATEX --------------------------------------------------------------------- u_h \approx \left(\begin{array}{c} u_0 \\ u_1 \\ \vdots \\ u_N \\ u_{N+1} \\ \end{array}\right) -------------------------------------------------------------------------------- Because of the Dirichlet boundary conditions the values $u_0 = g(0)$ and $u_{N+1} = g(1)$ are known. We will also not make actual use of $\widehat{q}_h$ but instead use only the well-defined values (i.e. values not indicated as arbitrary by $*$) of ---- LATEX --------------------------------------------------------------------- f_h \approx \left(\begin{array}{c} * \\ f_1 \\ \vdots \\ f_N \\ f_{N+1} \\ \end{array}\right) = \left(\begin{array}{c} * \\ f(x_1) \\ \vdots \\ f(x_N) \\ * \\ \end{array}\right) -------------------------------------------------------------------------------- Ingredients needed for a multigrid solver ========================================= We will develop step by steps the ingredients needed for a multigrid solver: - Direct solver - Computing the residual - Smoother - Restriction - Prolongation We will start with computing the residual and the direct solver. Computing the residual ====================== For computing the residual we will use $u_h$ and $f_h$ and require: - Both vectors have length $N+1$ - In $u_h$ the $u_0$ and $u_{N+1}$ are initialized with $g(0)$ and $g(1)$ respectively. - In $f_h$ the components $f_1, \dots, f_N$ are set to $f(x_1), \dots, f(x_n)$. Then we can compute $\widehat{r}_h = \widehat{q}_h - A_h \widehat{u}_h$ as follows: - For $i = 1, \dots, N$ - $r_i = f_i - \frac{1}{h^2} (2 u_i - u_{i-1} - u_{i+1})$ In the code given below this algorithm is implemented in function evalResidual (which works on $r_h$ instead of $\widehat{r}_h)$. Direct solver ============= ---- LATEX --------------------------------------------------------------------- T_N = \left(\begin{array}{ccccc} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots& \ddots& \ddots& \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array}\right) -------------------------------------------------------------------------------- the Cholesky factorization is known to be $T_N = L_N D_N L_N^T$ with ---- LATEX --------------------------------------------------------------------- L_N = \left(\begin{array}{ccccc} 1 & & & & \\ \ell_1& 1 & & & \\ & \ddots& \ddots & & \\ & & \ell_{N-2}& 1 & \\ & & & \ell_{N-1}& 1 \\ \end{array}\right), \quad \ell_i = -\frac{i}{i+1} -------------------------------------------------------------------------------- and ---- LATEX --------------------------------------------------------------------- D_N = \text{diag}(d_1, \dots, d_N), \quad d_i = \frac{i+1}{i} -------------------------------------------------------------------------------- Code for exercise ================= Make yourself familiar with the ingredients: - Function discretizeProblem is used to initialize vectors $u_h$ and $f_h$ as follows: - Both vectos have length $N+1$. - For $f_h$ its components $f_1, \dots, f_N$ are initialized with function values $f(x_1), \dots, f(x_n)$. The components $f_0$ and $f_{N+1}$ are not initialized (and therefore should not be used). - For $u_h$ only $u_0$ and $u_{N+1}$ are initialized with $g(0)$ and $g(1)$ respectively. - Function lInfNorm computes the $\ell_\infty$ norm of a grid function. - Function evalError computes for an approximation $u_h$ and the exact, contiguous solution $u$ the vector $e_h =$(e_0, e_1, \dots, e_{N+1}) = (u(x_0) - u_0, u(x_1) - u_1, \dots, u(x_{N+1}) - u_{N+1})$. - Function evalResidual computes the residual as described above. - Function direct_solver is your code to fix. :import: session31/ex01/direct_solver.cpp Compiling the code ================== ---- SHELL(path=session31/ex01) ------------------------------------------------ g++ -Wall -I /home/numerik/pub/hpc/ws19/session31 -std=c++17 +++ -o direct_solver direct_solver.cpp ./direct_solver -------------------------------------------------------------------------------- Exercise ======== - Write down an algorithm that gets$f_h$and$u_h$and overwrites the unknown components of$u_h$: -$u_h$is a vector of length$N+2$and its values$u_0$and$u_{N+1}$are initialized with$g(0)$and$g(1)\$ respectively. - Complete the code with your implementation and test your algorithm. :navigate: up -> doc:index next -> doc:session31/page02