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 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/traits.hpp>
#include <hpc/matvec/print.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;

}

//------------------------------------------------------------------------------

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/ws18/session25/ -std=c++17 -o direct_solver direct_solver.cpp direct_solver.cpp: In instantiation of 'void direct_solver(const VectorF&, VectorU&) [with T = double; VectorF = hpc::matvec::DenseVector; VectorU = hpc::matvec::DenseVector; typename std::common_type >::type::is_DenseVector:: value, bool>::type, typename std::enable_if >::type::is_DenseVector:: value, bool>::type>::type = 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.

• Complete the following code with your implementation and test your algorithm: