Construction of Prolongation and Restriction Operator for Geometric Multigrid (2D-FEM): Resulting in a Decreasing Solution
Consider the following problem, $$ -\Delta u(x) = f(x), \qquad x \in \Omega \ u(x) = 0,\qquad x \in \partial \Omega$$ with $\Omega = [0,1]\times [0,1]$ being the domain and $\partial \Omega$ being the boundary. The source term $f(x) \in L^2(\Omega)$ is a smooth source term. The exact solution is known.
I have generated the two following triangulations of the domain, dimensioned $81\times 81$ and $289 \times 289$.
Using regular P1-FEM, I have constructed a linear system for each of these two meshes, $$ A_1 u_1 = f_1, \qquad A_2 u_2 = f_2.$$
The direct solutions of the two systems are working as intended, with the finer mesh having a relative error (of $L^2$-norm) of $0.0077$ compared to the known exact solution and the coarser mesh being a bit less accurate.
Now, rather than using a direct solver, I am attempting to construct a two-grid solution for this problem. I am using a V-cycle method, with the following structure:
1. Start with the fine grid, and perform 3 iterations of the dampened Jacobi method (weight $\omega = 2/3$), using zero as initial guess. The damped Jacobi method I have implemented has the following structure $$ u^{n+1} = u^n + \omega D^{-1} (f - Au^n),$$ where $D$ is the diagonal of $A$. This seem to work as intended, as the $L^2$-norm of the relaxed solution is lower than the $L^2$-norm for the initial guess. We can denote this relaxed solution by $u_2$.
2. Compute the residual of the relaxed solution by $r_2 = A_2 u_2 - f_2$.
3. Restrict the residual $r_2$ to the coarse mesh. Here, I am using an injection restriction operator, where I simply remove the finer nodes and keep the nodes which are part of both meshes. 935
The general structure of this operator is $$ R = \begin{bmatrix}1 & 0 & 0 & \cdots & & & 0\ 0 & 1 & 0 & \cdots & & & 0\ 0 & 0 & \ddots & \cdots & & & 0 \ 0 & 0 & 0 & 1 & 0 & \cdots & 0\end{bmatrix}, \qquad (81 \times 289)$$ The implementation of this operator seem to work as I intend it to. The dimensions become correct, and the "finer" nodes are removed from the residual. We can denote the restricted residual $r_2$ as $Rr_2$.
4. Solve the coarse system directly. Using Matlab notation and denoting the solution as $e$, $e = A_1 \setminus (Rr_2)$. This results in a $L^2$ error which is larger than the relaxed solution $u_2$, which seems odd.
5. Apply a prolongation operator to the solution $e$. The prolongation operator that I am using is a linear interpolation. The nodes which are shared between the coarse grid and the fine grid are left unchanged, while the finer nodes are computed as the average of the two closest coarse nodes (on the same edge). This results in the following structure of the prolongation operator, $$ I = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \ 0 & 1 & 0 & \cdots & 0 \ 0 & 0 & \ddots & \cdots & 0 \ 0 & 0 & \cdots & & 1 \ 0.5 & 0.5 & 0 & \cdots & 0 \ 0 & 0 & 0.5 & 0.5 & 0\ \vdots & \vdots & \vdots & \vdots & \vdots \ 0 & 0 & 0.5 & 0 & 0.5\end{bmatrix}, \qquad (289 \times 81).$$ Again, the dimension of this prolongation operator seem correct, and the vectors seem to update as intended. We can denote $u_1 = Ie$ as the prolonged solution $e$ onto the finer grid.
6. I update the original relaxed solution $u_2$ as $u_2 = u_2 - Ie$.
7. 490