================ 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 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/ws18/session25/ -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 following code with your implementation and test your algorithm: