A-posteriori-steered $p$-robust multigrid and domain decomposition methods with optimal step-sizes for mixed finite element discretizations of elliptic problems

TL;DR

This work tackles the challenge of solving saddle-point systems from mixed finite element discretizations of elliptic problems on highly graded meshes with arbitrary polynomial degree $por all $p\ge 0$. It introduces an a-posteriori-steered multigrid solver and an overlapping additive Schwarz domain decomposition method whose updates yield a computable a posteriori error estimator for the algebraic error, ensuring $p$-robust contraction at every iteration. A key theoretical contribution is a multilevel, $p$-robust stable decomposition for divergence-free Raviart–Thomas spaces in both 2D and 3D, enabling rigorous contraction proofs and estimator efficiency that do not degrade with $p$ or mesh grading. Numerical experiments corroborate the theory, showing $p$-robust iteration counts, effective localization of the algebraic error, and adaptive smoothing strategies that reduce computational effort on graded meshes. Overall, the methods offer scalable, provably robust solvers for mixed FEM on complex grids, with direct applicability to porous-media flow and related saddle-point problems.

Abstract

In this work, we develop algebraic solvers for linear systems arising from the discretization of second-order elliptic partial differential equations by saddle-point mixed finite element methods of arbitrary polynomial degree $p \ge 0$ on possibly highly graded simplicial meshes. We present a multigrid and a two-level domain decomposition approach in two and three space dimensions, steered by a posteriori estimators of the algebraic error. First, we extend [Miraçi, Papež, and Vohralík, SIAM J. Sci. Comput. 43 (2021), S117-S145] to the mixed finite element setting. Extending the multigrid procedure itself is rather natural. To obtain analogous theoretical results, however, a $p$-robust multilevel stable decomposition of the velocity space is needed. In two space dimensions, we can treat the velocity space as the curl of a stream-function Lagrange space, for which the previous results apply. In three space dimensions, we design a novel stable decomposition by combining a one-level high-order local stable decomposition of [Falk and Winther, Found. Comput. Math. (2025), DOI 10.1007/s10208- 025-09700-2] and a multilevel lowest-order stable decomposition of [Hiptmair, Wu, and Zheng, Numer. Math. Theory Methods Appl. 5 (2012), 297-332]. This allows us to prove that our multigrid solver contracts the algebraic error at each iteration $p$-robustly and, simultaneously, that the associated a posteriori estimator is $p$-robustly efficient. Next, we use this multilevel methodology to define a two-level domain decomposition method where the subdomains consist of overlapping patches of coarse-level elements sharing a common coarse-level vertex. We again establish a $p$-robust contraction of the solver and $p$-robust efficiency of the a posteriori estimator. Numerical results presented both for the multigrid approach and the domain decomposition method confirm the theoretical findings.

Figures (8)

• Figure 1: Illustration of the set $\mathcal{B}_j$; the refinement $\mathcal{T}_{j}$ (dotted lines) of the mesh $\mathcal{T}_{j-1}$ (full lines).
• Figure 2: Patch in the two-level overlapping additive Schwarz method: coarse grid $\mathcal{T}_{H}$ (left), fine grid $\mathcal{T}_{h}$ (right). The highlighted patch consists of four coarse elements of $\mathcal{T}_{H}$ which share a vertex and form a subdomain. The subdomains (coarse patches) are discretized with the fine grid $\mathcal{T}_{h}$.
• Figure 3: Left: mesh (zoomed in) used for the well wavefront test case \ref{['Wellwavefront']}, obtained after $J=10$ local adaptive refinements and consisting of 9814 elements. Center: variations of the coefficient $c(x, y)$ for the piecewise constant diffusion tensor ${{\bm K} = c(x,y)I}$ across the domain for the checkerboard test case \ref{['Checkerboard']}. Right: mesh used for the checkerboard test case \ref{['Checkerboard']}, obtained after $J=20$ local adaptive refinements and consisting of 1038 elements.
• Figure 4: Effectivity index $\eta^i_{\mathrm{alg}} / \|{\bm K}^{-1/2}(\bm{u}_{J} - \bm{u}_{J}^{i})\|$ of the a posteriori estimator $\eta^i_{\mathrm{alg}}$ for the algebraic error $\|{\bm K}^{-1/2}(\bm{u}_{J} - \bm{u}_{J}^{i})\|$. Multigrid solver of Algorithm \ref{['Definition_solver']} (left) and domain decomposition solver of Algorithm \ref{['AS_solver']} (right).
• Figure 5: Contraction factors of the solvers given by $\|{\bm K}^{-1/2}(\bm{u}_{J} - \bm{u}_{J}^{i+1})\| / \|{\bm K}^{-1/2}(\bm{u}_{J} - \bm{u}_{J}^{i})\|$. Multigrid solver of Algorithm \ref{['Definition_solver']} (left) and domain decomposition solver of Algorithm \ref{['AS_solver']} (right).

Theorems & Definitions (31)

• Remark 2.1: Other choices of discrete spaces
• Remark 4.1: Choice of patch subdomains
• Remark 4.3: Initialization and its cost
• Remark 4.4: Cost per iteration
• Remark 4.5: A-posteriori-steered solver
• Remark 4.6: Compact formulas
• Lemma 4.7: Optimal step-sizes
• proof
• Lemma 4.8: Norm of the levelwise corrections as sum of norms of the local corrections
• Remark 5.2: Initialization and its cost