FAS-multigrid slower than linear defect correction?

I have implemented a V-Cycle multigrid solver using both a linear defect correction (LDC) and full approximation scheme (FAS).

My problem is the following: Using LDC the residual is reduced by a factor of ~0.03 per cycle. The FAS implementation does converge with a linear factor too, but the factor is only ~0.58. Thus FAS needs about 20 times the number of cycles.

Most of the code is shared, the only difference are the down/up calculations, LDC uses

down: $u_H:=0,\quad b_H:=I_h^H(b_h-L_hu_h)$

up: $u_h:=u_h+I_H^hu_H$

and FAS uses

down: $u_H:=I_h^Hu_h,\quad b_H:=I_h^Hb_h+L_HI_h^Hu_h-I_h^HL_hu_h$

up: $u_h:=u_h+I_H^h(u_H-I_h^Hu_h)$

My test setting is from Brigg's "A Multigrid Tutorial, Second Edition", p. 64, has the analytical solution

$u(x,y)=(x^2-x^4)(y^4-y^2) \quad$ with $x,y\in [0,1]^2$

and the equation is $Lv=\Delta u=:b$ using the typical linear 5-point stencil as the Laplace-Operator $L$. The initial guess is $v=0$.

Changing the test setting, e.g. to the trivial $u(x,y)=0$ using an initial guess of $v=1$ results in nearly the same convergence factors.

Since only the down/up code differs, the LDC results comply with the book and the FAS at least seems to work too, I have no Idea why it is so much slower in the same linear setting.

There is one odd behavior in both LDC and FAS that I cannot explain yet that only happens if the initial guess is bad (e.g. $=0$ but also in my full multigrid experiments where the interpolation to the new fine grid increases the residual from $10^{-15}$ to $10^{-1}$): If I increase to number of post correction relaxations to a very high number such that the solution is solved to machine precision on the coarse grid, it looses nearly all digits when going one step up to the next fine grid.

Since a picture says more than words:

// first cycle, levels 0-4
// DOWN
VCycle top 4, start               res_norm 3.676520e+02 // initial residual
VCycle top 4, cycle 0, current 4, res_norm 3.676520e+02
VCycle top 4, cycle 0, current 4, res_norm 1.520312e+02 // relaxed (2 iterations)
VCycle tau_norm 2.148001e+01 (DEBUG calculation)
VCycle top 4, cycle 0, current 3, res_norm 1.049619e+02 // restricted
VCycle top 4, cycle 0, current 3, res_norm 5.050392e+01 // relaxed (2 iterations)
VCycle top 4, cycle 0, current 2, res_norm 3.518764e+01 // restricted
VCycle top 4, cycle 0, current 2, res_norm 1.759372e+01 // relaxed (2 iterations)
VCycle top 4, cycle 0, current 1, res_norm 1.234398e+01 // restricted
VCycle top 4, cycle 0, current 1, res_norm 4.728777e+00 // relaxed (2 iterations)
VCycle top 4, cycle 0, current 0, res_norm 3.343750e+00 // restricted
// coarsest grid
VCycle top 4, cycle 0, current 0, res_norm 0.000000e+00 // solved
// UP
VCycle top 4, cycle 0, current 1, res_norm 3.738426e+00 // prolonged
VCycle top 4, cycle 0, current 1, res_norm 0.000000e+00 // relaxed (many iterations)
VCycle top 4, cycle 0, current 2, res_norm 1.509429e+01 // prolonged (loosing digits)
VCycle top 4, cycle 0, current 2, res_norm 2.512148e-15 // relaxed (many iterations)
VCycle top 4, cycle 0, current 3, res_norm 4.695979e+01 // prolonged (loosing digits)
VCycle top 4, cycle 0, current 3, res_norm 0.000000e+00 // relaxed (many iterations)
VCycle top 4, cycle 0, current 4, res_norm 1.469312e+02 // prolonged (loosing digits)
VCycle top 4, cycle 0, current 4, res_norm 9.172812e-24 // relaxed (many iterations)


I'm not sure if there can be only a few digits gained per cycle or if this indicates an error during the interpolation to the fine grid. If it is the latter case, how can the LDC achieve the by-the-book residual ratios of ~0.03 when using always 2 relaxations?

I don't directly know your answer as I mainly use FAS instead of correction since I do multigrid for nonlinear problems, but some thoughts you can look into:

• You're applying a linear correction scheme to a linear problem, so it's not shocking that it does very well.

• Consider your boundary conditions: make sure you're doing them correctly, and also note that complicated BCs could look completely different on the coarse grid, making corrections there not so useful.

• Double-check your treatment of the source term; I recall screwing up something in the prolongation stage related to that term when I wrote it for poisson.

• I've never seen the need to iterate to convergence on the coarsest grid. A solution there depends on the fine grid residual being correct, which it isn't. You're trying to push those errors out of the domain / smooth them out. If you're fully converged on the coarsest grid in an early iteration, your solution is naturally pretty far from the correct fine grid solution because your residuals there aren't up to date. This is almost certainly the reason for why you're seeing residuals jump up in the prolongation stage.

• Also, try a relaxation factor for both the restriction and prolongation operators, say 0.75.

If it helps, this is what my FAS residual history looked like for a poisson problem using single grid through full 7V cycles. I believe the relaxation factor was 0.75, and I was using a 3-stage RK scheme as a smoother with a single iteration at each grid level.

• Thanks for your reply, the linear case and simple BC (square border=0) are only the first step, testing real cases will be done after the "easy" settings work. I'm not sure if I understand what you mean with a relaxation factor for restriction and prolongation. I currently use bilinear interpolation for prolongation and half-weighting for restriction. Jul 10, 2013 at 19:20
• By relaxation I mean for your Prolongation stage alter it to: $u_h:=u_h+\alpha I_H^h(u_H-I_h^Hu_h)$ where $0< \alpha <1$ is a relaxation factor. Generally the more complicated the solution, the lower this factor needs to go. In solutions with a lot of discontinuities I sometimes have to push that down to around 0.6, but generally 0.75-0.85 works. Jul 10, 2013 at 21:52
• Good to know. In my setting this only slows down the convergence ratio by $(1-\alpha)$ but I'll keep it in mind when testing more complex data. Jul 12, 2013 at 18:40
• @Aurelius you mention that converging on the coarse grid is not necessary. I agree with your reasoning, but convergence proofs in the literature (for the linear case) do assume that the coarse grid solve is exact. I am not aware of any reference (for the linear or nonlinear case) in which it is stated that the coarse grid solve shouldn't be exact, and was wondering if you could cite a reference for this? I would be very interested in seeing this myself Nov 14, 2013 at 17:10
• @KeeranBrabazon I don't have a reference for this either, and I'm honestly not intimately familiar with the details of convergence proofs for multigrid. I'd suggest looking for any early literature that introduces that relaxation factor. That factor is common to all modern multigrid implementations I've seen, and it's intuitively true that it wouldn't be needed if exact course solutions were necessary/desirable. For an intuitive proof, I just imagine what boundary conditions look like for the coarsest grid vs the finest. It's easy to imagine them creating very different solutions. Nov 18, 2013 at 16:53

If you are using a vertex-centered discretization, then state restriction should be injection rather than the full-weighted residual restriction that it appears you use. That is, replace $I_h^H$ with $\hat I_h^H$ when restricting the state. Using full-weighted restriction for state produces aliasing of high-frequency components of the state which after applying $u^h \gets u^h (u^H - I_h^H u^h)$ results in new noise on the same scale as before the coarse correction (boundary conditions are especially likely culprits for this effect). Use injection, $\hat I_h^H u^h$, and this problem should go away.

Spectrally, state restriction only needs high secondary order (accurate preservation of low frequencies), but the primary order (aliasing of high frequencies) does not matter. Injection has primary order 0 and infinite secondary order. Meanwhile, residual restriction needs both primary and secondary order to be positive (at least). See section 4.3 of Achi Brandt's Multigrid Guide.

When designing MG methods, it's also good to look at error rather than residuals and to make sure you weight the norm appropriately.

• Good points, and I had failed to mention something along those lines. One important aspect of practical use of multigrid is the choice of smoother: you want one that dampens high frequency errors as fast as possible, which addresses the issue you describe. Jul 12, 2013 at 12:06
• @Aurelius From the log provided, you can see that the smoother is not the problem. Recall that Silpion is using the same smoother as for the defect correction MG, which converges properly. Jul 12, 2013 at 13:28
• Thanks for the link to Brandt's Multigrid Guide, I'll read it thoroughly after finishing Brigg's Multigrid Tutorial. I solved the problem now (see other answer) but I'm currently using full weighting for both state and residual restriction. Using injection doesn't seem to work in my setting, residual ratios change to $>0.8$ and both residual and error L2-norms stop decreasing very soon. Do you have any ideas why injection fails? Jul 12, 2013 at 19:02

I solved by problem now. I stored $u_{old}^H=I_h^Hu^h$ when going down during the V-cycle and reused it later in

$u^h\leftarrow u^h+I_H^h(u^H-I_h^Hu^h)=u^h+I_H^h(u^H-u_{old}^H)$.

The problem was before going down again from $H$ to $2H$, $u_{old}^H$ was relaxed in place. Storing a copy before the relaxation steps helped. Since $u_{old}^H$ was only required in FAS, it didn't show up in linear calculations.