6 Iterative Solvers

As model problem, we consider the homogeneous Dirichlet boundary value problem ( u(x) = 0    $\forall$x $\in$ $\partial$$\Omega$) for the Poisson equation in the unit square $\Omega$ : = (0, 1)² in its weak formulation:

Find u $\in$ H¹₀($\Omega$) such that

$$\int_{{\Omega}}^{} \nabla^{T}_{} \nabla \int_{{\Omega}}^{} \forall \bf v \in \Omega$$ (1)

We use linear finite elements for the discretization and achieve the linear system of equations

$$\underline{{u}} \underline{{f}}$$ (2)

6.1 $\omega$-Jacobi solver

Let us denote the diagonal of matrix K by D = diag(K). Now, we can formulate the $\omega$-Jacobi iteration

$$\underline{{u}}^{{k+1}}_{} \underline{{u}}^{{k}}_{} \omega \underline{{f}} \underline{{u}}^{k}_{}$$ (3)

You will find a sequential version of the Jacobi solver in the directory Solution/jacseqc with the following functions in addition to the functions from P 5.

GetMatrix(nx, ny, xl, xr, yb, yt, sk, id, ik, f)

in which matrix K and right hand side $\underline{{f}}$ are calculated using the function FunctF(x,y) for describing $\bf f$(x). Note, that only coordinates and element connectivities are related to the rectangular domain - all other parts in this routine are written for general 2d-domains.
The Dirichlet boundary conditions are set in SetU(nx, ny, u). Alternatively, one could use FunctU(x,y). These b.c. are applied in

ApplyDirichletBC(nx, ny, neigh, u, sk, id, ik, f)

via penalty method.
The solver itself is implemented in

JacobiSolve(nx, ny, sk, id, ik, f, u )

and uses

GetDiag(nx, ny, sk, id, d)

to get the diagonal from matrix K. Matrix-times-vector is realized in

CrsMult(iza, ize, w, u, id, ik, sk, alfa) .

A vector $\underline{{u}}$ can be saved in file named name by calling

SaveVector(name, u, nx, ny, xl, xr, yb, yt, ierr) ,

such that

         gnuplot jac.dem

will give an impression of that vector.

Implement a parallel version of the sequential code !

6.2 Multigrid solver

Here, we will solve (2) by means of a multigrid iteration. You will find a sequential version of the multigrid solver using a $\omega$-Jacobi iterations a smoother and, for simplicity, also as a coarse grid solver in the directory Solution/mgseqc . The following functions are added to the functions from previous sections.

AllocLevels(nlevels, nx, ny, pp, ierr)

allocates the memory for all levels and stores the pointers in variable pp which is a structure PointerStruct. For details see mg.h . The matrices and right hand sides will be calculated in

IniLevels(xl, xr, yb, yt, nlevels, neigh, pp, ierr) .

FreeLevels(nlevels, pp, ierr)

frees the allocated memory at the end. Note, that we need for multigrid a different procedure to apply the Dirichlet b.c.

ApplyDirichletBC1(nx, ny, neigh, u, sk, id, ik, f)

The multigrid solver

MGMSolver(pp, crtl, level, neigh)

uses crtl as predefined structure ControlStruct to control the multigrid cycle (mg.h). The same holds for each mg-iteration

MGM(pp, crtl, level, neigh) .

The mg components smoothing, coarse grid solver, interpolation and restriction are defined in

JacobiSmooth(nx, ny, sk, id, ik, f, u, dd, aux, maxiter )

JacobiSolve2(nx, ny, sk, id, ik, f, u, dd, aux )

Interpolate(nx, ny, w, nxc, nyc, wc, neigh)

Restrict(nx, ny, w, nxc, nyc, wc, neigh) .

A vector $\underline{{u}}$ on level level can be saved in file named name by calling

SaveVector(name, pp->u[level], pp->nx[level], pp->ny[level],

u, nx, ny, xl, xr, yb, yt, ierr) ,

such that

         gnuplot mgm.dem

will give an impression of that vector.

Implement a parallel version of the sequential code !