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

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.

💡 Research Summary

This paper addresses the solution of linear systems that arise from mixed finite element (MFE) discretizations of second‑order elliptic partial differential equations in two and three spatial dimensions. The discretizations use Raviart–Thomas (RT) or Brezzi–Douglas–Marini (BDM) spaces of arbitrary polynomial degree p ≥ 0 on possibly highly graded simplicial meshes. The authors develop two algebraic solvers— a multigrid (MG) method and a two‑level overlapping domain decomposition (DD) method—both steered by a posteriori estimators of the algebraic error.

The central technical contribution is a polynomial‑degree‑robust (p‑robust) multilevel stable decomposition of the divergence‑free velocity space. In two dimensions the velocity space is identified as the curl of a stream‑function Lagrange space, allowing the authors to reuse the p‑robust decomposition from their earlier work on conforming elements. In three dimensions they combine a recent one‑level high‑order local stable decomposition (Falk & Winther, 2025) with a multilevel lowest‑order decomposition (Hiptmair, Wu & Zheng, 2012) to obtain a novel p‑robust splitting. This decomposition guarantees that the constants in the stability estimates are independent of p, which is essential for the subsequent convergence analysis.

The multigrid algorithm consists of a V‑cycle with zero pre‑smoothing and a single post‑smoothing step performed by an additive Schwarz (block‑Jacobi) smoother. At each level j a line‑search is carried out to compute an optimal step size λᵢⱼ. The resulting Pythagorean identity

‖K⁻¹ᐟ²(uᴶ−uⁱ⁺¹ᴶ)‖² = ‖K⁻¹ᐟ²(uᴶ−uⁱᴶ)‖² − Σⱼ λᵢⱼ‖K⁻¹ᐟ²ρᵢⱼ‖²

implies a built‑in a posteriori error estimator

ηᵢₐₗg² = Σⱼ λᵢⱼ‖K⁻¹ᐟ²ρᵢⱼ‖².

The estimator provides a guaranteed lower bound for the algebraic error and, together with the Pythagorean relation, yields a contraction estimate

‖errorⁱ⁺¹‖ ≤ α‖errorⁱ‖, 0<α<1,

with α independent of p. Moreover, the estimator is p‑robustly efficient: β‖errorⁱ‖ ≤ ηᵢₐₗg ≤ ‖errorⁱ‖ for a constant β>0 independent of p.

The two‑level DD method uses overlapping patches of coarse‑level elements that share a common coarse vertex as subdomains. Each subdomain solves a fine‑level problem with the same block‑Jacobi smoother, while a coarse‑grid solve enforces global consistency. The same a posteriori estimator is employed, and the authors prove a p‑robust contraction analogous to the multigrid case. No additional smoothing steps or damping parameters are required, making the method parameter‑free.

The theoretical results are organized as follows: (1) construction of the p‑robust multilevel stable decomposition; (2) reliability and efficiency of the a posteriori estimator; (3) convergence of the multigrid solver; (4) convergence of the domain decomposition solver; (5) extension of all results to highly graded meshes generated by newest‑vertex bisection, provided the meshes remain shape‑regular.

Numerical experiments cover polynomial degrees p = 0,…,5, a variety of mesh refinements (including strongly non‑uniform meshes), heterogeneous diffusion tensors, and mixed boundary conditions. The experiments confirm that (i) the estimator ηᵢₐₗg closely matches the true algebraic error, (ii) both solvers converge in a number of iterations that does not grow with p, and (iii) adaptive selection of the number of smoothing steps further reduces computational cost.

In summary, the paper delivers a unified framework for p‑robust, a posteriori‑driven multigrid and domain decomposition solvers for mixed finite element discretizations of elliptic problems. The approach is mathematically rigorous, works on arbitrary simplicial meshes (including highly graded ones), and requires no user‑tuned parameters. The authors suggest that the methodology can be extended to nonlinear, time‑dependent, and multiphysics problems, and that future work will explore non‑nested meshes, adaptive p‑refinement, and large‑scale parallel scalability.