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

<span class=$81\times 81$ nodes" /><span class=$289 \times 289$ nodes" />

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. 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. Lastly, I apply the same dampened Jacobi iteration (same as step 1), with the solution $u_2$ as the initial guess. The final solution $u$ is obtained.

I think that this is a reasonable approach, and the same principle works for my test program (1D-FDM). Now, the issue that I am facing is that the solution after performing these steps is much lower than the corresponding reference solution.

The following two Figures illustrate the issue well, Direct solution Two-grid solution

The same tendencies are shown in both solution, but the two-grid solution is much smaller than the direct solution.

Comparing the $L^2$-errors of each solution, the following result is shown:

  • Two-grid solution: 0.0079995 (rel. error 0.24)
  • Direct solution: 0.000256 (rel. error 0.0077)

This is a huge difference. Adding more levels, or increasing the number of cycles (the principle above is only one cycle) makes the issue worse.

My question is, what am I doing wrong? This is a very odd behavior, but the implementation of the restriction and prolongation operators agree with what I have come up with using pen and paper.

My first thought is that there is some kind of scaling factor missing, since both solution have the same tendencies but have different sizes. I would assume that it is missing in the restriction operator, but adding any kind of factor there seems wrong, since the row-sum should be 1.

I am running out of ideas right now, and would appreciate any suggestion very much! Thanks in advance!

Edit: If I multiply the restriction operator by a scalar, say $4$, the $L^2$-error of the two-grid solution becomes $0.000612$ (rel.error 0.0184), which is significantly better. But I really cannot understand why this would be a better approach? Could the scalar be related to the number of levels, the ratio of the number of nodes between the levels?

Edit 2: These are the current results, for different attempts on the restriction operator:

  1. Restriction operator based only on injection:
    • $L^2$ error: 0.0327
  2. Same as (1) but with a factor 4 in front:
    • $L^2$ error: 0.0022
  3. Restriction operator as transpose of prolongation:
    • $L^2$ error: 0.0093
  4. Same as (3) but use injection on boundary nodes:
    • $L^2$ error: 0.0054
  • 2
    $\begingroup$ The usual choice is to make sure that the restriction operator is the transpose of the prolongation operator, or the other way around. That's not the case for you. $\endgroup$ Nov 2 '20 at 16:30
  • $\begingroup$ Thank you for your comment! Yes, I have seen that choice for the restriction operator in some literature. However, that approach has left me with 1) wrong boundaries, and 2) an odd behavior in the middle of the domain, a sudden decrease in the solution making it look like a "hole". However, it is of course possible that I had done a mistake in the implementation of that operator, but I could not find any. Therefore, I was recommended to attempt this easier restriction operator using only injection. $\endgroup$
    – BastiSand
    Nov 2 '20 at 16:34
  • $\begingroup$ I have also seen a recommendation of using a 7-point stencil for this kind of triangulation. But I cannot see how that would cause me less problems than this simple operator. Any thoughts? $\endgroup$
    – BastiSand
    Nov 2 '20 at 16:36
  • 1
    $\begingroup$ I'm most familiar with multigrid as a precondition, and in that case you need to make sure that the matrix $P{2\leftarrow 1}A_1R_{1\leftarrow 2}$ is symmetric. Here, $A_1$ is the matrix on the coarse level 1, $R_{1\leftarrow 2}$ is the matrix that restricts from level 2 (the "fine level") to the coarse level, and $P_{2\leftarrow 1}$ is the prolongation. You achieve this by choosing $R_{1\leftarrow 2} = P_{2\leftarrow 1}^T$ $\endgroup$ Nov 2 '20 at 22:45
  • 1
    $\begingroup$ The way you want to think about it is what happens if you already had the correct solution on the fine mesh: What would be the correction you get on the coarse mesh. Think, for example, if you wanted to solve $-\Delta u=0$ with $u=1$ on the boundary. The exact solution is $u=1$, i.e., the vector of nodal values $U$ has all ones. Move that vector through one multigrid cycle: You need to get the same $U$ out of it again. If you don't (on a piece of paper), then the operators you have defined are not useful. $\endgroup$ Nov 3 '20 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.