Schoeberl thesis online reading club - firedrakeproject/firedrake GitHub Wiki

Over the summer of 2020 we read through Joachim Schöberl's thesis on parameter robust multigrid methods.

In the analysis of the multigrid methods, the theory is that of Hackbusch (approximation and smoothing properties). These days, tighter bounds are available. Matt likes the paper of Brannick et al., which expands on theory originally developed in Falgout, Vassilevski, and Zikatanov (2005) and also covered in Xu and Zikatanov's algebraic multigrid review.

Schöberl considers the parameter-dependent problem requiring that A^ε is coercive. Patrick wonders if we can consider how the theory maps on to the case where we only have that the eigenvalues of A^ε are bounded away from zero.

A request that we carefully unpick, when going through proofs, all of the constants in the inequalities of the form a ≼ b (so that we see where all the bits came from).

We covered the basics of finite element approximation and showed how conforming discretisations can fail for parameter dependent problems. Ivan notes that his masters thesis contains an example of this for the case of the Poisson equation with varying coefficients (in section 3.3).

We got a little hung up on Hilbert space interpolation. Gonzalo suggests that this is going to be used in the multigrid convergence proof later (this theory is in Bramble). Brenner & Scott have a reasonably self-contained overview of Hilbert space interpolation and its relation to finite elements in Chapter 14 of their book.

We introduced more finite element theory. In particular Scott-Zhang interpolation, inverse inequalities, and partition of unity methods. We then started on theory for parameter-dependent problems and the migration to a mixed formulation. We showed that the mixed problem B^ε has parameter-dependent continuity constant in the parameter-dependent V×εQ-norm. We finished with discussion of the remark before Theorem 2.8 and the introduction of the dual ||·||_Q,0 norm, which is introduced such that ΛV has closed range with respect to this new norm. Gonzalo wrote up some more detailed notes explaining this point.

We covered the parameter-dependent case of the mixed problem and introduced some new norms, along with Fortin operators as a useful tool for moving from continuous to discrete inf-sup conditions.

We gathered our thoughts before embarking on examples, particularly, we reminded ourselves of the various different norms in play.

• ||u||_V = (u, u)_V^½
• ||u||_A^ε = A^ε(u, u)^½
  • A¹ is V-elliptic, so norm-equivalence ||u||_V ≃ ||u||_A¹
• ||p||_c = c(p, p)^½
• ||u, p||_V×εc = (||u||²_V + ε||p||_c²)^½
• ||p||_Q,0 ≃ sup_{v ∈ V} c(Λ v, p)/||v||_V
• ||p||_Q = (||p||_Q,0² + ε||p||_c²)^½
• ||u, p||_X = ||u||_V + ||p||_Q

Also have some norms on the discrete spaces.

• ||u_h||_Ah^ε = A_h^ε(u, u)^½
  • In general must check V_h ellipticity of A_h¹, if we have it, then ||u||_V ≃ ||u||_Ah¹
• ||p_h||_Qh,0 ≃ sup_{vh ∈ Vh} c(Λ_h v_h, p)/||v_h||_V
• ||p_h||_Qh = (||p_h||_Qh,0² + ε||p_h||_c²)^½
• ||u_h, p_h||_Xh = ||u_h||_V + ||p_h||_Qh

Existence of a Fortin operator is required so that we have equivalence of ||·||_Q and ||·||_Qh (and then ||·||_X and ||·||_Xh).

We then briefly summarised the abstract framework for showing that the discretisation will be parameter robust. Given a bilinear form in primal variables we must:

• Identify a, c, Λ, and ε
• Check A¹ for V-ellipticity
• Introduce dual variable p = ε^-1Λu
• Define ||·||_Q
• Show ||·||_Q ≃ ||·||_c

We must then show stability of the continuous mixed formulation. We then discretise and:

• Pick V_h ⊂ V and Q_h ⊂ Q
• Define Λ_h
• (Maybe) separately check V_h-ellipticity of A¹_h
• Pick c-stable space splitting Q_h = Q_0,h + Q_1,h
  • show ||p₀||_Q ≼ ε^½||p₀||_c for all p₀ ∈ Q_h,0
  • Construct a Fortin operator I^F_hwith
    1. c(ΛI^F_h u, q₁) = c(Λ u, q₁) for all q₁ ∈ Q_h,1, u ∈ V.
    2. ||I^F_h||_V ≼ 1

Some bibliographic comments. In the general case the continuity of the Fortin operator depends on both the polynomial degree of V_h and the shape regularity constants in the mesh. Many discretisations have continuity constants that degrade with (at least) the square root of the aspect ratio and/or discretisation degree. There is a brief summary of results where this is not the case in Apel, Kempf, Linke, and Merdon (2020). Note particularly that the Crouzeix–Raviart H¹-nonconforming element has a continuity constant of 1 on all simplex meshes (Apel, Nicaise, and Schöberl, 2001).

We went briefly recapped some of the abstract framework again and then applied it to the section on nearly incompressible materials, before starting on the more complicated results needed for the Reissner–Mindlin plate. The construction and demonstration of the regularity result in (2.87) was somewhat confusing, but it's actually just done by applying the Riesz representation theorem to find the Riesz representer of each linear functional.

Some more detailed notes are available.

We applied the abstract framework to the Reissner-Mindlin plate taking us to the end of chapter 2. Once again, Gonzalo provided some more detailed notes.

We went through the introduction of the abstract additive Schwarz framework. In particular the proof that the splitting norm is equal to the norm defined by the additive Schwarz preconditioner. This tool is used to determine spectral bounds for the operator. We looked at the upper bound and saw how to estimate it in terms of the number of overlapping subspaces. For a more detailed overview of these methods, Michael Holst has some nice notes.

We went through the spectral inequalities for the one-level domain decomposition method. We note a repeated technique in the proofs when using the splitting norm. We construct bounds on the energy norm by providing an explicit splitting that obeys one side of an inequality we want to show. To relate this to the preconditioner (Additive-Schwarz) norm, we go via the splitting norm and then apply the Additive-Schwarz Lemma (Theorem 3.1). Since this splitting norm is defined as inf-ing over all valid splittings, if we have a witness splitting then we obtain an inequality of the form |||u||| ≤ ∑ ||u_i||_Ah. This method is used, for example, in proving Theorem 3.4.

Noting that the one-level method does not provide optimal spectral constants, we then constructed a two-level system with a coarse grid. The crucial properties relating the coarse grid problem are collected in Lemma 3.5. While the two-level method provides optimal bounds, it is not scalable, since the coarse grid must grow with the overall problem size. We will therefore turn to multigrid methods to solve this.

We finished off chapter 3 by introducing the abstract framework for multigrid analysis. This uses the approach of Hackbusch (1982). In particular we introduced the norms involved in the approximation and smoothing properties. The theory is covered in a very general way that does not require full elliptic regularity, and hence is littered with fractional Sobolev norms. For second-order problems with full regularity the relevant norms simplify to H¹ for the energy norm, and a scaled L²-like norm for the local norm ||·||_l,0̅. The chapter finishes by collecting the convergence results of Hackbusch for abstract multigrid convergence. In particular that the contraction factor of a W-cycle with fixed numbers of smoothing iterations is bounded, and that the condition number of a variable V-cycle is also bounded (increasing numbers of iterations on coarse grids). For the former, one is referred to Hackbusch (1985), for the latter, Bramble (1993) and Bramble, Pasciak, and Xu (1991). The notation of Bramble (1993) is translated into that used here.

For a little more detail on the multigrid analysis, in addition to the papers referenced above, see these notes of Schöberl.

The next steps are to glue together the discretisation results of Chapter 2 and the multigrid theory of Chapter 3 to produce parameter robust preconditioners. This will entail producing space splittings and the subsequent approximation and smoothing properties such that the bounds have no dependence on the small parameter.

We kicked off with Chapter 4, which develops the theorems for parameter robust multigrid solvers. First we saw an example of how Jacobi preconditioning does not control the ε-dependence in the condition number of the Timoshenko beam. We then saw how an overlapping Schwarz smoother does control the ε-dependence. There is one wrinkle which will complicate later analysis is that the h-dependence gets worse: from h² for the Jacobi preconditioner to h⁴ for the block Jacobi version.

We then introduced and went through the proof of Theorem 4.1 which provides the conditions for a parameter-robust Schwarz smoother based on some space composition.

Essentially we need the three conditions (4.2–4.4) which describe stable splittings of functions in V_h, stable splittings of kernel functions, and a requirement on the bound of the c_h-norm in terms of the Q_h,0-norm.

There is an error in the statement of the first inequality in (4.5), it should read

(c₂(h) + c₁(h)c₃(h)²)^-1.

To prove the spectral inequalities we use the finite overlap lemma for the upper bound and the additive Schwarz lemma for the lower bound. This allows us to transition from the preconditioner norm to a sum over the space-decomposed pieces. Splitting the u_i into pieces, we then use the stable splitting assumptions. To get back into the V_h norm, we use the stability estimates coming from the Brezzi conditions (Theorem 2.8), along with the ellipticity of A_h to move from V_h to A_h-norms.

Next up, we'll see how to add a coarse grid to the preconditioner, obtaining optimal spectral bounds. This will require constructing a robust prolongation operator to satisfy the requirements of Lemma 3.5.

We went through Theorem 4.2 which sets out the requirements for an optimal, parameter-robust, two-grid preconditioner. The idea is to find a splitting of the space Q_h such that the prolongation of a coarse grid kernel function can be modified to live in the fine grid kernel by modifying only those dofs that live on "interior" fine grid entities. This way the problems decouple and can be solved efficiently. This is shown using an example of the Timoshenko beam in Figure 4.2. This figure is slightly misleading because it looks like both variables have the same scale. I wrote some code to show what is going (with the right scaling).

In constructing the space splittings, we need to ensure that the interpolation operator satisfies the conditions of Lemma 3.5, in particular that the splitting is stable for u_f := u_h - R^V_HI^F_Hu_h, constructed by round-tripping a fine grid function through the interpolation and prolongation.

This has to be done for each problem and discretisation in turn. This is somewhat involved, which motivates why it is often preferable to use the abstract framework of Hilbert complexes and exact sequences to construct the space decompositions. In these cases, when working on nested meshes, the trivial prolongation is continuous in the energy norm and we don't have to do this dance to satisfy the conditions of Lemma 3.5.

We finished up by going through, in somewhat less detail, 4.3--4.6 looking at multigrid convergence. Many of the details build on the results proving optimality of the two-level method. We must work harder to prove the approximation and smoothing properties than in the two-level scheme because the coarse grid is weaker. The approach is to transition to the mixed problem and restrict the analysis to a subspace in which the constraints are always satisfied. We prove equivalence of the primal and mixed algorithms when operating on this space. In contrast to parameter-independent multigrid the local norm one uses in the approximation property has ε-dependent terms that are used to control the continuity and energy norm and the size of the correction to the dual variable.

Two smoothers are analysed, one on the mixed form, which is shown to be equivalent to that of Braess and Sarazin (1997), the other in primal variables. For the mixed form one can show convergence of order O(m^-1), whereas for the primal form one only gets O(m^-1/2) as a function of the number of smoothing steps.

We skipped over a lot of the details of the proofs, but they rely heavily on the constraint space. A critical requirement of the space decomposition that is pointed out, and not enjoyed by most incarnations of Vanka relaxation, is that having picked a space decomposition for V_l and Q_l, we require that Λ_l V_l,i ⊂ Q_l,i. Just as Λ_l V_l ⊂ Q_l.

Finally we looked through the results which confirm that the scheme does indeed work numerically. Interestingly, despite the very weak theoretical convergence bounds, a W(2,2) cycle already produces mesh- and parameter-independent condition numbers, and a V(1,1) cycle is very close to mesh- and parameter-independent when used as preconditioner.