# Multigrid solver

For the sake of simplicity, we will only consider one-dimensional boundary problems of the form

$\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

$A_h \widehat{u}_h = \widehat{q}_h,\quad h = \frac{1}{N+1}$

with

$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

$\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

$\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

$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

$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

## 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

$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

$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

$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 =$$(e0, e1, \dots, e{N+1}) = (u(x0) - u0, u(x1) - u1, \dots, u(x{N+1}) - u_{N+1})$. • Function evalResidual computes the residual as described above. • Function direct_solver is your code to fix. #include <cmath> #include <hpc/matvec/axpy.hpp> #include <hpc/matvec/copy.hpp> #include <hpc/matvec/densevector.hpp> #include <hpc/matvec/iterators.hpp> #include <hpc/matvec/print.hpp> #include <hpc/matvec/traits.hpp> //------------------------------------------------------------------------------ using namespace hpc; using namespace hpc::matvec; template <typename T> const T & assertEqual(const T &a, const T &b) { assert(a==b); return a; } //------------------------------------------------------------------------------ template < typename F, typename G, typename T, template<typename> class VectorF, template<typename> class VectorU, Require< Dense<VectorF<T>>, Dense<VectorU<T>> > = true > void discretizeProblem(const F &f, const G &g, VectorF<T> &f_h, VectorU<T> &u_h) { auto N = assertEqual(f_h.length(), u_h.length()) - 2; auto h = T(1)/(N+1); u_h(0) = g(0); u_h(N+1) = g(1); for (std::size_t i=1; i<N; ++i) { u_h(i) = 0; f_h(i) = f(i*h); } // actually we don't have to initialize f_h(0), f_h(N+1) so we store some // good and evil numbers there: f_h(0) = 42; f_h(N+1) = 666; } //------------------------------------------------------------------------------ template < typename T, template<typename> class Vector, Require< Dense<Vector<T>> > = true > T lInfNorm(const Vector<T> &x_h) { T res = 0; for (auto [i, xi] : x_h) { i=i; if (std::abs(xi) > res) { res = std::abs(xi); } } return res; } //------------------------------------------------------------------------------ template < typename U, typename T, template<typename> class VectorU, template<typename> class VectorE, Require< Dense<VectorU<T>>, Dense<VectorE<T>> > = true > void evalError(const U &u, const VectorU<T> &u_h, VectorE<T> &e_h) { auto N = assertEqual(u_h.length(), e_h.length()) - 2; auto h = T(1)/(N+1); for (auto [i, ei] : e_h) { ei = u(i*h) - u_h(i); } } //------------------------------------------------------------------------------ template < typename T, template<typename> class VectorU, template<typename> class VectorF, template<typename> class VectorR, Require< Dense<VectorF<T>>, Dense<VectorU<T>>, Dense<VectorR<T>> > = true > void evalResidual(const VectorU<T> &u_h, const VectorF<T> &f_h, VectorR<T> &r_h) { auto N = assertEqual(u_h.length(), r_h.length()) - 2; auto h = T(1)/(N+1); auto h2 = h*h; r_h(0) = r_h(N+1) = 0; for (std::size_t i=1; i<=N; ++i) { r_h(i) = f_h(i) - 1./h2 * (2*u_h(i) - u_h(i-1) - u_h(i+1)); } } //------------------------------------------------------------------------------ template < typename T, template<typename> class VectorF, template<typename> class VectorU, Require< Dense<VectorF<T>>, Dense<VectorU<T>> > = true > void direct_solver(const VectorF<T> &f_h, VectorU<T> &u_h) { auto N = assertEqual(f_h.length(), u_h.length()) - 2; auto h = T(1)/(N+1); auto h2 = h*h; // TODO: Your code } //------------------------------------------------------------------------------ int main() { auto f = [](double x) { return -2; }; auto g = [](double x) { return 1 + 3*x; }; auto u = [](double x) { return (x+1)*(x+1); }; std::size_t N = 10; DenseVector<double> f_h(N+2), u_h(N+2), r_h(N+2), e_h(N+2); discretizeProblem(f, g, f_h, u_h); fmt::printf("f_h =\n"); print(f_h); fmt::printf("u_h =\n"); print(u_h); fmt::printf("Calling solver: direct_solver(f_h, u_h);\n"); direct_solver(f_h, u_h); evalResidual(u_h, f_h, r_h); fmt::printf("u_h =\n"); print(u_h); fmt::printf("r_h =\n"); print(r_h); fmt::printf("Eval approximation error of u_h:\n"); evalError(u, u_h, e_h); fmt::printf("e_h =\n"); print(e_h); fmt::printf("\n\nIn short:\n"); fmt::printf("lInfNorm(r_h) = %.5e\n", lInfNorm(r_h)); fmt::printf("lInfNorm(e_h) = %.5e\n", lInfNorm(e_h)); }  ## Compiling the code theon$ g++ -Wall -I /home/numerik/pub/hpc/ws19/session31 -std=c++17 -o direct_solver direct_solver.cpp
direct_solver.cpp: In instantiation of 'void direct_solver(const VectorF<T>&, VectorU<T>&) [with T = double; VectorF = hpc::matvec::DenseVector; VectorU = hpc::matvec::DenseVector; typename std::common_type<typename std::enable_if<typename std::remove_reference<VectorX<T> >::type::is_DenseVector:: value, bool>::type, typename std::enable_if<typename std::remove_reference<MatrixB<T> >::type::is_DenseVector:: value, bool>::type>::type <anonymous> = 1]':
direct_solver.cpp:153:27:   required from here
direct_solver.cpp:129:10: warning: unused variable 'h2' [-Wunused-variable]
auto h2 = h*h;
^~
theon$./direct_solver f_h = 42.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 0.00 666.00 u_h = 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.00 Calling solver: direct_solver(f_h, u_h); u_h = 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.00 r_h = 0.00 119.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 -2.00 484.00 0.00 Eval approximation error of u_h: e_h = 0.00 1.19 1.40 1.62 1.86 2.12 2.39 2.68 2.98 3.31 3.64 0.00 In short: lInfNorm(r_h) = 4.84000e+02 lInfNorm(e_h) = 3.64463e+00 theon$ 

## 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.