# A posteriori algebraic error estimates and nonoverlapping domain decomposition in mixed formulations: energy coarse grid balancing, local mass conservation on each step, and line search

Manuela Bastidas Olivares, Akram Beni Hamad, Martin Vohralík, Ivan Yotov

► To cite this version:

Manuela Bastidas Olivares, Akram Beni Hamad, Martin Vohralík, Ivan Yotov. A posteriori algebraic error estimates and nonoverlapping domain decomposition in mixed formulations: energy coarse grid balancing, local mass conservation on each step, and line search. 2025. hal-04843665v2

HAL Id: hal-04843665

https://inria.hal.science/hal-04843665v2

Preprint submitted on 13 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A posteriori algebraic error estimates and nonoverlapping domain decomposition in mixed formulations: energy coarse grid balancing, local mass conservation on each step, and line search

Manuela Bastidas Olivares†‡ Akram Beni Hamad†‡ Martin Vohralík†‡

Ivan Yotov §

June 13, 2025

Abstract

We consider iterative algebraic solvers for saddle-point mixed finite element discretizations of the model Darcy flow problem. We propose a posteriori error estimators of the algebraic error as well as a nonoverlapping domain decomposition algorithm. The estimators control the algebraic error from above and from below in a guaranteed and fully computable way. The distinctive feature of the domain decomposition algorithm is that it produces a locally mass conservative approximation on each iteration. Both the estimate and the algorithm rely on a coarse grid solver, a subdomain Neumann solver, and a subdomain Dirichlet solver. The algorithm also employs a line search to determine the optimal step size, leading to a Pythagoras formula for the algebraic error decrease in each iteration. We suppose that the fine mesh is a refinement of a coarse mesh, where both meshes need to be formed by simplices or rectangular parallelepipeds. Numerical experiments illustrate the theoretical developments and confirm the efficiency of the algebraic error estimates and of the domain decomposition algorithm.

**Keywords:** mixed finite elements, saddle-point system, iterative algebraic solver, algebraic error, a posteriori estimate, domain decomposition method, coarse grid solver, subdomain Neumann solver, subdomain Dirichlet solver, balancing, local mass conservation, line search.

# 1 Introduction

Saddle-point systems arise in many scientific and technical applications, such as fluid mechanics, structural mechanics, and electromagnetism. One typically seeks to find algebraic vectors **U** and **P** that satisfy the following system of linear algebraic equations:

$$ \begin{pmatrix} \mathbf{A} & \mathbf{B}^t \\ \mathbf{B} & 0 \end{pmatrix} \begin{pmatrix} \mathbf{U} \\ \mathbf{P} \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ \mathbf{F} \end{pmatrix}. \qquad (1.1) $$

Here, the matrix **A** is symmetric and positive definite and **B** has a full (row) rank, whereas **F** is a right-hand side vector. The overall system matrix is indefinite, of saddle-point type, posing difficulties for numerical solvers, cf. the overview paper Benzi, Golub, and Liesen [6].

Domain decomposition solvers allow for iterative algebraic approximation of the solution to (1.1) where smaller systems are solved on each iteration in parallel on smaller subdomains.

<footnote>†Inria, 48 rue Barrault, 75647 Paris, France, (manuela.bastidas-olivares@inria.fr, akram.beni-hamad@inria.fr, martin.vohralik@inria.fr).</footnote>

<footnote>‡CERMICS, Ecole nationale des ponts et chaussées, IP Paris, 77455 Marne la Vallée, France.</footnote>

<footnote>§Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA (yotov@math.pitt.edu). Partially supported by Inria Paris Visiting Professorship and NSF grants DMS 2111129 and DMS 2410686.</footnote>
In the context of saddle-point systems arising from mixed finite elements, the early works of Glowinski and Wheeler [19], Ewing and Wang [17], Mandel [22], Mathew [23], and Cowsar, Mandel, and Wheeler [11] laid important grounds. Since then, numerous related approaches appeared, see, e.g., Toselli [36], Arbogast et al. [4], Tu [37], Ganis and Yotov [18], Sousedík [35], Šístek, Březina, and Sousedík [33], Ciarlet, Jamelot, and Kpadonou [10], Oh et al. [27], Dobrev et al. [13], Solovský, Fučík, and Šístek [34], Boon et al. [8], and the references therein. However, to our knowledge, no available nonoverlapping domain decomposition solver preserves on each iteration the basic physical property behind (1.1): the local mass conservation expressed by $\mathbb{B}U = F$.

Whenever one considers (an iterative) approximate solution $(\mathbf{U}^j, \mathbf{P}^j)^t$ to the exact solution $(\mathbf{U}, \mathbf{P})^t$, there arises the misfit between them, the algebraic error. Algebraic error is a central notion in numerical linear algebra but, somehow surprisingly, computable and guaranteed a posteriori error estimates on it are seldom available, namely for bounding the error from above (bounds from below are typically easier); we refer to the discussion and references in Rey, Gosselet, and Rey [32] and Papež et al. [29, 28]. In the context of saddle-point systems and domain decomposition methods, this question was partly treated in Ali Hassan et al. [1], where an a posteriori error estimate on the overall error (distance between the piecewise polynomial $\mathbf{u}_h^j$ associated with $\mathbf{U}^j$ and the exact solution $\mathbf{u}$ of the underlying partial differential equation) has been derived, with a component associated with the algebraic error.

In our first main result, Theorem 7.1, we derive a guaranteed upper bound on the algebraic error

$$ (\mathbf{U} - \mathbf{U}^j)^t \mathbb{A}(\mathbf{U} - \mathbf{U}^j) = |||\mathbf{u}_h - \mathbf{u}_h^j|||^2 \leq (\eta^j)^2, \quad (1.2) $$

where $\eta^j$ is an a posteriori error estimate fully computable from $(\mathbf{U}^j, \mathbf{P}^j)^t$ (or, equivalently, the associated piecewise polynomials $(\mathbf{u}_h^j, p_h^j)$). This result is complemented in Theorem 8.3, where we also present a lower bound on the algebraic error of the form

$$ (\underline{\eta}^j)^2 \leq (\mathbf{U} - \mathbf{U}^j)^t \mathbb{A}(\mathbf{U} - \mathbf{U}^j) = |||\mathbf{u}_h - \mathbf{u}_h^j|||^2. \quad (1.3) $$

In our second main result, Algorithm 8.1, we design a nonoverlapping domain decomposition solver with the following properties: 1) on each iteration $j$, upon using Construction 5.1 in (8.2), we obtain an approximate solution satisfying the local mass conservation $\mathbb{B}U^j = F$, or, more precisely,

$$ \mathbf{u}_h^j \in V_{h,k}^{gN} \subset \mathbf{H}_{gN}(\mathrm{div}, \Omega), \quad \nabla \cdot \mathbf{u}_h^j = f, \quad (1.4) $$

on each iteration $j \geq 1$. 2) As per the Theorem 8.3, we have the guaranteed a posteriori algebraic error estimates (1.2)–(1.3) on each iteration $j \geq 1$. This in particular makes it suitable for use in mesh adaptive methods with inexact algebraic solvers, see, e.g., [15, 20] and the references therein. 3) The line search (8.3) gives the error decrease formula of Theorem 8.2 on each iteration $j \geq 1$. 4) Algorithm 8.1 essentially employs a coarse-grid solver in combination with subdomain Neumann and Dirichlet solvers in an additive Schwarz manner. All are formulated as quadratic constrained energy minimizations associated with the current residual, which is natural both from mathematical and physical perspective. 5) The built-in coarse-grid solver of step 2 of Construction 5.1 is a solution to the “balancing” problem for domain decomposition in mixed formulations [11]. It is closely related to the equilibration technique used in a posteriori analysis for numerical discretizations of partial differential equations [16, 28]. This also prevents the appearance of other issues such as “interior node compatibility” linked to Neumann conditions, nonconforming discretizations, or non-nested spaces. 6) Algorithm 8.1 uses functional writing by relying directly on the piecewise polynomial spaces and functions therein and is thus basis-independent, in contrast to a writing like (1.1), where the bases of the piecewise polynomial spaces need to be chosen first and potentially influence the solver behavior and properties. Congruently, we develop a basis-independent analysis where we exploit information and tools from function spaces, in contrast to being confined to linear algebraic information and tools.
This contribution is organized as follows. In Section 2, we introduce the model Darcy flow problem, the function spaces on the continuous level, and the weak formulation. Section 3 then settles the discrete level. Section 4 then presents the nonoverlapping domain decomposition configuration. Our key technical tool, the equilibrated flux of Construction 5.1 (yielding the local mass conservation), is presented in Section 5. The reconstruction of the scalar variable (potential) is then based on trace lifting of Construction 6.1 and elementwise postprocessing of Construction 6.2 of Section 6. Theorem 7.1 yielding the guaranteed upper bound on the algebraic error (1.2) is presented in Section 7, and our nonoverlapping domain decomposition solver of Algorithm 8.1 is detailed in Section 8. These theoretical results are numerically illustrated in Section 9. Finally, some technical details of elementwise postprocessing are collected in Appendix A.

# 2 Continuous setting

We introduce here the model Darcy flow problem, the function spaces on the continuous level, and the weak formulation.

## 2.1 Problem

Let $\Omega \subset \mathbb{R}^d$, $1 \le d \le 3$, be an open bounded connected domain (interval for $d=1$, polygon for $d=2$, polyhedron for $d=3$), with Lipschitz boundary $\partial\Omega$ such that $\partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N}$. We suppose that $\Gamma_D$ and $\Gamma_N$ are also Lipschitz, relatively open, and correspond to a set of mesh boundary faces as defined below. We denote by $\mathbf{n}_\Omega$ the outward unit normal vector to $\partial\Omega$ and consider the following diffusion problem: find $p: \Omega \to \mathbb{R}$ such that

$$ -\nabla \cdot (\mathbf{S}\nabla p) = f \quad \text{in } \Omega, \qquad (2.1a) $$

$$ p = 0 \quad \text{on } \Gamma_D, \qquad (2.1b) $$

$$ -\mathbf{S}\nabla p \cdot \mathbf{n}_{\Omega} = g_N \quad \text{on } \Gamma_N. \qquad (2.1c) $$

In (2.1), the primal unknown $p$ is called the potential, whereas the dual variable

$$ \mathbf{u} := -\mathbf{S}\nabla p \qquad (2.2) $$

is called the flux. This is a basic model in countless applications, namely thermo-diffusion processes and groundwater flow.

## 2.2 Assumptions on the data

For the sake of simplicity, with respect to the mesh $\mathcal{T}_h$ of Section 3.1 and with reference to the piecewise polynomial spaces from (3.3), (3.4) and (3.5a) below, we assume for the data in (2.1):

**Assumption 2.1 (Piecewise polynomial data).** The data $\mathbf{S}$, $f$, and $g_N$ are such that:

Ass1 The diffusion tensor $\mathbf{S}$ is bounded and symmetric positive definite. It is constant elementwise on simplices and constant elementwise and diagonal on rectangular parallelepipeds.

Ass2 The source term $f$ satisfies $f|_{K} \in W_{h,k}(K)$ for all $K \in \mathcal{T}_h$.

Ass3 The Neumann boundary condition $g_N$ satisfies $g_N \in (\mathbf{V}_{h,k}) \cdot \mathbf{n}_{\Omega}|_{\Gamma_N}$.
## 2.3 Domain and spaces

For a general open bounded connected domain $\mathcal{D} \subset \mathbb{R}^d$, we denote by $|\mathcal{D}|$ its Lebesgue measure and by $L^2(\mathcal{D})$ the space of scalar-valued and square-integrable functions. We let $(\cdot, \cdot)_\mathcal{D}$ represent the inner product on $L^2(\mathcal{D})$ and $\|\cdot\|_\mathcal{D}$ the corresponding norm; we set $L^2(\mathcal{D}) := [L^2(\mathcal{D})]^d$, where we use the same notation for the inner product and norm as in the scalar case. When $\mathcal{D} = \Omega$, we drop the subscripts. We let $\langle \cdot, \cdot \rangle_\gamma$ be the scalar product for $L^2(\gamma)$, where $\gamma = \partial\mathcal{D}$ or a subset of it.

Let $H^1(\mathcal{D})$ be the Sobolev space of $L^2(\mathcal{D})$ functions whose weak derivatives are square integrable, $H^1(\mathcal{D}) := \{v \in L^2(\mathcal{D}) : \nabla v \in L^2(\mathcal{D})\}$. Further, let $H(\text{div}, \mathcal{D})$ be the space of vector-valued functions in $L^2(\mathcal{D})$ whose weak divergences are square integrable, $H(\text{div}, \mathcal{D}) := \{v \in L^2(\mathcal{D}) : \nabla \cdot v \in L^2(\mathcal{D})\}$. We also use $H_{0,\mathcal{D}}^1(\Omega)$, the subspace of $H^1(\Omega)$ formed by the functions vanishing on $\Gamma_\mathcal{D}$ in the sense of traces. Similarly, define the set

$$ H_{g_N}(\text{div}, \Omega) := \{v \in H(\text{div}, \Omega) : v \cdot n_{\Omega} = g_N \text{ on } \Gamma_N\}, \quad (2.3) $$

where $v \cdot n_\Omega = g_N$ on $\Gamma_N$ means that $(v, \nabla\phi) + (\nabla \cdot v, \phi) = \langle g_N, \phi \rangle_{\partial\Omega} = \langle g_N, \phi \rangle_{\Gamma_N}$ for all functions $\phi \in H_{0,\mathcal{D}}^1(\Omega)$. Similarly, the space $H_0(\text{div}, \Omega)$ is defined as $H_{g_N}(\text{div}, \Omega)$ in (2.3) but with $g_N$ replaced by 0. We define the energy norm of a function $\psi \in L^2(\mathcal{D})$ by

$$ \|\psi\|_{\mathcal{D}} := \|S^{-\frac{1}{2}}\psi\|_{\mathcal{D}}. $$

We suppose

$$ (f, 1) = \langle g_N, 1 \rangle_{\partial\Omega} \quad (2.4) $$

when $\Gamma_N = \partial\Omega$. Congruently, $L_*^2(\Omega) := L^2(\Omega)$ when $\Gamma_N \neq \partial\Omega$; when $\Gamma_N = \partial\Omega$, $L_*^2(\Omega)$ stands for those $L^2(\Omega)$ functions whose mean value vanishes.

## 2.4 Weak formulation

In mixed form, problem (2.1) is equivalent to the following quadratic constrained minimization problem: seek for a function $u$ such that

$$ u = \arg\min_{\substack{v \in H_{g_N}(\text{div}, \Omega) \\ \nabla \cdot v = f}} \|v\|^2, \quad (2.5) $$

cf., e.g., [7]. The Euler-Lagrange optimality conditions for (2.5) read: find $u \in H_{g_N}(\text{div}, \Omega)$ with $\nabla \cdot u = f$ in $\Omega$ satisfying

$$ (S^{-1}u, v) = 0 \quad \forall v \in H_0(\text{div}, \Omega) \text{ with } \nabla \cdot v = 0 \text{ in } \Omega. $$

If we introduce a Lagrange multiplier associated with the constraint $\nabla \cdot u = f$ in $\Omega$, we can equivalently seek for the pair $(u, p) \in H_{g_N}(\text{div}, \Omega) \times L_*^2(\Omega)$ such that

$$ \begin{aligned} (S^{-1}u, v) - (p, \nabla \cdot v) &= 0 & \forall v \in H_0(\text{div}, \Omega), \\ (\nabla \cdot u, q) &= (f, q) & \forall q \in L_*^2(\Omega). \end{aligned} $$

When $\Gamma_N = \partial\Omega$, the Neumann compatibility condition (2.4) is important to guarantee the existence and uniqueness of the solutions.

# 3 Discrete setting

This section lays down the discrete setting: meshes, discrete (piecewise polynomial) spaces, and the mixed finite element discretization together with its equivalent hybridization.
## 3.1 Meshes

Let $\mathcal{T}_h$ be a partition of the domain $\Omega$ into a finite number of nonoverlapping elements $K$. We assume that the elements $K \in \mathcal{T}_h$ are either simplices or rectangular parallelepipeds. Let $|\mathcal{T}_h|$ be the number of elements in $\mathcal{T}_h$ and denote by $h_K$ the diameter of the element $K \in \mathcal{T}_h$. We assume that adjacent elements completely share their common vertex, edge, or face, so we avoid hanging nodes. Denote the set of all vertices (for $d=1$), edges (for $d=2$), or faces (for $d=3$) of $\mathcal{T}_h$ by $\mathcal{F}_h := \mathcal{F}_h^{\text{ext}} \cup \mathcal{F}_h^{\text{int}}$, where $\mathcal{F}_h^{\text{ext}}$ is the set of boundary vertices, edges, or faces lying on $\partial\Omega$ and $\mathcal{F}_h^{\text{int}}$ is the set of interior vertices, edges, or faces. Henceforth we only speak of *faces*. We assume that each boundary face is entirely contained either in $\bar{\Gamma}_D$ or in $\bar{\Gamma}_N$. Consequently, we denote $\mathcal{F}_h^D$ and $\mathcal{F}_h^N$ the sets of faces lying on $\bar{\Gamma}_D$ and $\bar{\Gamma}_N$, respectively. Moreover, for each $K \in \mathcal{T}_h$, we denote by $\mathcal{F}_K$ the set of its faces. For each $F \in \mathcal{F}_h$, $\mathbf{n}_F$ is a fixed face normal, coinciding with $\mathbf{n}_\Omega$ on $\partial\Omega$ and otherwise arbitrary but fixed. The shape-regularity parameter of the mesh $\mathcal{T}_h$ is the positive real number

$$ \kappa_{\mathcal{T}_h} := \max_{K \in \mathcal{T}_h} \frac{h_K}{\rho_K}, $$

where $\rho_K$ is the diameter of the largest ball contained in $K$.

For a sufficiently smooth function $\psi$ that is not necessarily univalued on the faces $F \in \mathcal{F}_h^{\text{int}}$ with $F \in \mathcal{F}_K \cap \mathcal{F}_{K'}$ and $\mathbf{n}_F$ pointing from $K$ to $K'$, with $K \in \mathcal{T}_h$ and $K' \in \mathcal{T}_h$, its jump and average are defined as follows:

$$ \llbracket \psi \rrbracket := \psi|_{K|_F} - \psi|_{K'|_F} \quad \text{and} \quad \llbracket \psi \rrbracket := \frac{1}{2} (\psi|_{K|_F} + \psi|_{K'|_F}). \qquad (3.1) $$

To achieve robustness with respect to diffusion inhomogeneities, we consider diffusivity-dependent weighted averages as in, e.g., [14]. For all $F \in \mathcal{F}_h^{\text{int}}$, we let

$$ \omega_F^- := \frac{\delta_K}{\delta_K + \delta_{K'}}, \quad \omega_F^+ := \frac{\delta_{K'}}{\delta_K + \delta_{K'}}, $$

where $\delta_K := \mathbf{n}_F \cdot \mathbf{S}|_K \mathbf{n}_F$, and define

$$ \llbracket \psi \rrbracket_w := \omega_F^- \psi|_{K|_F} + \omega_F^+ \psi|_{K'|_F}. \qquad (3.2) $$

## 3.2 Piecewise polynomial spaces

For $k \ge 0$, we define the piecewise polynomial pressure space

$$ W_{h,k} := \{ q_h \in L^2(\Omega) : q_h|_K \in W_k(K) \text{ for all } K \in \mathcal{T}_h \}, $$

$$ W_{h,k}^* := \{ q_h \in W_{h,k} : (q_h, 1) = 0 \text{ when } \Gamma_N = \partial\Omega \}, $$

where

$$ W_k(K) := \begin{cases} \mathcal{P}_k(K) & K \text{ simplex,} \\ \mathcal{Q}_{k,...,k}(K) & K \text{ rectangular parallelepiped,} \end{cases} \qquad (3.3) $$

with $\mathcal{P}_k(K)$ being the space of polynomials of total degree at most $k$ on $K$ and $\mathcal{Q}_{k_1,...,k_d}$ being the space of polynomials of degree at most $k_1$ in $x_1, \dots, k_d$ in $x_d$.

Similarly, for $k \ge 0$, the piecewise polynomial velocity space is

$$ \mathbf{V}_{h,k}^{\text{dc}} := \{ \mathbf{v}_h \in L^2(\Omega) : \mathbf{v}_h|_K \in \mathcal{RTN}_k(K) \text{ for all } K \in \mathcal{T}_h \}, $$
where the space $RTN_k(K)$ is the Raviart-Thomas-Nédélec mixed finite element space of order $k$ proposed in [31, 26] and defined on each element $K \in \mathcal{T}_h$ as follows:

$$ \mathcal{RTN}_k(K) := \left\{ \begin{array}{ll} \mathcal{P}_k(K) + x\mathcal{P}_k(K) & K \text{ simplex,} \\ Q_{k+1,k}(K) \times Q_{k,k+1}(K) & K \text{ rect. par., } d=2, \\ Q_{k+1,k,k}(K) \times Q_{k,k+1,k}(K) \times Q_{k,k,k+1}(K) & K \text{ rect. par., } d=3, \end{array} \right. \quad (3.4) $$

with $\mathcal{P}_k(K) := [\mathcal{P}_k(K)]^d$. This defines the so-called broken Raviart-Thomas-Nédélec spaces, without any normal continuity imposed.

The discrete subspaces/subsets of $H(\text{div}, \Omega)$, $H_0(\text{div}, \Omega)$, and $H_{gN}(\text{div}, \Omega)$ are then written as

$$ V_{h,k} := \{v_h \in H(\text{div}, \Omega) : v_h|_K \in \mathcal{RTN}_k(K) \text{ for all } K \in \mathcal{T}_h \}, \quad (3.5a) $$

$$ V_{h,k}^0 := \{v_h \in V_{h,k} : v_h \cdot n_{\Omega}|_{\Gamma_N} = 0\}, \quad (3.5b) $$

$$ V_{h,k}^{gN} := \{v_h \in V_{h,k} : v_h \cdot n_{\Omega}|_{\Gamma_N} = g_N\}. \quad (3.5c) $$

Finally, we define the pressure trace space

$$ \Psi_{h,k} := \prod_{F \in \mathcal{F}_h} \Psi_k(F), \quad (3.6) $$

where for all faces $F \in \mathcal{F}_h$, one has

$$ \Psi_k(F) := \left\{ \begin{array}{ll} \mathcal{P}_k(F) & K \text{ simplex,} \\ \mathcal{P}_k(F) & K \text{ rectangular parallelepiped, } d=2, \\ Q_{k,k}(F) & K \text{ rectangular parallelepiped, } d=3. \end{array} \right. \quad (3.7) $$

Similarly, we denote $\Psi_{h,k}^{dc}$ the broken pressure trace space

$$ \Psi_{h,k}^{dc} := \prod_{K \in \mathcal{T}_h} \prod_{F \in \mathcal{F}_K} \Psi_k(F). \quad (3.8) $$

The above spaces are constructed in a way that $\nabla \cdot \mathcal{RTN}_k(K) = W_k(K)$ together with $(\mathcal{RTN}_k(K) \cdot n_K)|_F = \Psi_k(F)$ for all $K \in \mathcal{T}_h$ and $F \in \mathcal{F}_K$. Later we will also use a postprocessing pressure space $M_{h,k} = \Pi_{K \in \mathcal{T}_h} M_{h,k}(K)$ as defined in [3, 5, 9, 38].

### 3.3 S-weighted $L^2(K)$-orthogonal projection

For every $K \in \mathcal{T}_h$ and $v \in L^2(K)$, the S-weighted $L^2(K)$-orthogonal projection onto the Raviart-Thomas-Nédélec space $\mathcal{RTN}_k(K)$ is a function $\Pi_k^{\mathcal{RTN}}v$ such that

$$ (S^{-1}(v - \Pi_k^{\mathcal{RTN}}v), w_h)_K = 0 \quad \forall w_h \in \mathcal{RTN}_k(K). \quad (3.9) $$

Below, we also apply it elementwise to $v \in L^2(\Omega)$.

### 3.4 Mixed finite elements

The discrete counterpart of the minimization problem (2.5) is, cf., [7],

$$ \begin{array}{l} u_h := \arg \min_{v_h \in V_{h,k}^{gN} \atop \nabla \cdot v_h = f} \|v_h\|^2. \quad (3.10) \\ v_h \in V_{h,k}^{gN} \end{array} $$
The Euler-Lagrange optimality conditions for (3.10) read: find $\mathbf{u}_h \in \mathbf{V}_{h,k}^{gN}$ with $\nabla \cdot \mathbf{u}_h = f$ in $\Omega$ and satisfying

$$ (\mathbf{S}^{-1}\mathbf{u}_h, \mathbf{v}_h) = 0 \quad \forall \mathbf{v}_h \in \mathbf{V}_{h,k}^0 \text{ with } \nabla \cdot \mathbf{v}_h = 0 \text{ in } \Omega. \qquad (3.11) $$

Problem (3.11) is equivalent to finding $(\mathbf{u}_h, p_h) \in \mathbf{V}_{h,k}^{gN} \times W_{h,k}^*$ such that

$$ (\mathbf{S}^{-1}\mathbf{u}_h, \mathbf{v}_h) - (p_h, \nabla \cdot \mathbf{v}_h) = 0 \quad \forall \mathbf{v}_h \in \mathbf{V}_{h,k}^0, \qquad (3.12a) $$

$$ (\nabla \cdot \mathbf{u}_h, q_h) = (f, q_h) \quad \forall q_h \in W_{h,k}^*. \qquad (3.12b) $$

The Lagrange multiplier $p_h$ approximates the potential $p$ on the mesh elements.

## 3.5 Hybridization

Imposing indirectly the normal trace continuity and Neumann boundary conditions, problem (3.12) is further equivalent to finding $(\mathbf{u}_h, p_h, \lambda_h) \in \mathbf{V}_{h,k}^{dc} \times W_{h,k}^* \times \Psi_{h,k}$ such that $\lambda_h = 0$ on $\Gamma_D$ and

$$ (\mathbf{S}^{-1}\mathbf{u}_h, \mathbf{v}_h) - (p_h, \nabla \cdot \mathbf{v}_h) + \sum_{K \in \mathcal{T}_h} \langle \lambda_h, \mathbf{v}_h \cdot \mathbf{n}_K \rangle_{\partial K} = 0 \quad \forall \mathbf{v}_h \in \mathbf{V}_{h,k}^{dc}, \qquad (3.13a) $$

$$ (\nabla \cdot \mathbf{u}_h, q_h) = (f, q_h) \quad \forall q_h \in W_{h,k}^*, \qquad (3.13b) $$

$$ \sum_{K \in \mathcal{T}_h} \langle \mathbf{u}_h \cdot \mathbf{n}_K, \mu_h \rangle_{\partial K} = \langle g_N, \mu_h \rangle_{\Gamma_N} \quad \forall \mu_h \in \Psi_{h,k} \text{ and } \mu_h = 0 \text{ on } \Gamma_D. \qquad (3.13c) $$

This is the so-called hybridized version of (3.12). The Lagrange multiplier $\lambda_h$ approximates the potential $p$ on the mesh faces.

# 4 Nonoverlapping domain decomposition configuration

Our results on a posteriori algebraic error estimation are generic but fit best the domain decomposition framework that we anyhow adopt for our iterative solver below.

## 4.1 Meshes

We henceforth suppose that the domain $\Omega$ is decomposed into $\mathcal{N}$ nonoverlapping subdomains $\Omega_i$, $i = 1, \dots, \mathcal{N}$. Here, the subdomains $\Omega_i$ are simplices or rectangular parallelepipeds forming a mesh of $\Omega$ in the sense of Section 3.1; we will call it a *coarse mesh*. We request that the subdomains $\Omega_i$ match with the mesh $\mathcal{T}_h$ and we denote the mesh of each subdomain $\Omega_i$ as $\mathcal{T}_{h_i} := \mathcal{T}_h|_{\Omega_i}$. For all $i = 1, \dots, \mathcal{N}$, let $\mathbf{n}_i$ be the unit normal pointing outward $\partial\Omega_i$, and denote the interface between two adjacent subdomains $\Omega_i$ and $\Omega_j$ by $\Gamma_{i,j}$. We also denote the set of all faces of $\mathcal{T}_{h_i}$ by $\mathcal{F}_{h_i}$, and in $\mathcal{F}_{h}^{\Gamma_{i,j}}$ we collect the faces that belong to the interface $\Gamma_{i,j}$. By $\mathcal{F}_H$, we denote the set of all faces of $\mathcal{T}_H$. All $\mathcal{T}_h$, $\mathcal{T}_{h_i}$, and $\mathcal{T}_H$ are matching meshes without hanging nodes, and we interpret $\mathcal{T}_h$ and $\mathcal{T}_{h_i}$ as a refinement of the coarse mesh $\mathcal{T}_H$, where $h, h_i$, and $H$ are the largest diameters of the elements in the meshes $\mathcal{T}_h, \mathcal{T}_{h_i}$, and $\mathcal{T}_H$ respectively.

## 4.2 Piecewise polynomial spaces

We define the local spaces

$$ W_{i,h,k} := \{ q_h \in L^2(\Omega_i) : q_h|_{K \in W_k(K)} \text{ for all } K \in \mathcal{T}_{h_i} \}, \quad W_{i,h,k}^* := W_{i,h,k} \cap L_*^2(\Omega_i), $$
where $L^2_*(\Omega_i)$ stands for those $L^2(\Omega_i)$ functions whose mean value vanishes over $\Omega_i$ when $|\partial\Omega_i \cap \Gamma_D| = 0$ ("interior" and "Neumann" subdomains) and $L^2(\Omega_i)$ otherwise. Moreover, on the coarse mesh $\mathcal{T}_H$, we will use the spaces

$$W_{H,0} := \{q_H \in L^2(\Omega) : q_H|_{\Omega_i} \in W_0(\Omega_i) \text{ for all } \Omega_i \in \mathcal{T}_H\}, \text{ and } W_{H,0}^* := W_{H,0} \cap L^2_*(\Omega).$$

Note that these are the lowest-order spaces in which $W_0(\Omega_i)$ denotes constant functions in $\Omega_i$. We let $\Pi^H$ denote the $L^2$-orthogonal projection onto $W_{H,0}^*$. We will also use

$$V_{i,h,k} := \{v_h \in \mathbf{H}(\mathrm{div}, \Omega_i) : v_h|_K \in \mathcal{RTN}_k(K) \text{ for all } K \in \mathcal{T}_{h_i}\},$$

$$V_{i,h,k}^0 := \{v_h \in \mathbf{H}(\mathrm{div}, \Omega_i) : v_h|_K \in \mathcal{RTN}_k(K) \text{ for all } K \in \mathcal{T}_{h_i} \text{ and } v_h \cdot \mathbf{n}_i = 0 \text{ on } \partial\Omega_i \setminus \Gamma_D\},$$

$$V_{H,0}^0 := \{v_H \in \mathbf{H}(\mathrm{div}, \Omega) : v_H|_{\Omega_i} \in \mathcal{RTN}_0(\Omega_i) \text{ for all } \Omega_i \in \mathcal{T}_H \text{ and } v_H \cdot \mathbf{n}_i = 0 \text{ on } \Gamma_N\}.$$

# 5 Equilibrated (balanced) flux reconstruction

Let $(\mathbf{u}_h^j, p_h^j) \in V_{h,k}^{\mathrm{dc}} \times W_{h,k}$ be an arbitrary approximation to the solution of (3.12). Note that it does not have to be included in the approximation space $V_{h,k}^{g_N} \times W_{h,k}^*$; in particular, $\mathbf{u}_h^j$ can be normal-trace discontinuous and $p_h^j$ does not have to satisfy the zero-mean value constraint $(p_h^j, 1) = 0$ when $\Gamma_N = \partial\Omega$. By the superscript $j \ge 0$, we precede that later on, in Section 8 below, these approximations will be obtained by an iterative domain decomposition algorithm with iteration index $j$. As the domain decomposition is nonoverlapping, normal-trace discontinuous iterates indeed naturally arise. In this section, we present our central tool, an equilibrated flux reconstruction and the corresponding potential reconstruction, formalizing the approach of [1, Section 5.3]. This yields $\mathcal{R}_F(\mathbf{u}_h^j, p_h^j) \in V_{h,k}^{g_N}$ with $\nabla \cdot (\mathcal{R}_F(\mathbf{u}_h^j, p_h^j)) = f$ and $\mathcal{R}_P(\mathbf{u}_h^j, p_h^j) \in W_{h,k}^*$, i.e., in particular a flux which is conforming (normal-trace continuous) and equilibrated (balanced, locally mass conservative).

## 5.1 Coarse–fine mesh construction of an equilibrated (balanced) flux and of a corresponding potential

Recall the definition of the average (3.1). The heart of our approach is the following construction:

**Construction 5.1** (Equilibrated (balanced) flux and the corresponding potential). Let $(\mathbf{u}_h^j, p_h^j) \in V_{h,k}^{\mathrm{dc}} \times W_{h,k}$ be arbitrary. Proceed in four steps:

1. Averaging on mesh faces.

From $\mathbf{u}_h^j \in V_{h,k}^{\mathrm{dc}}$, create an auxiliary function $\mathbf{u}_h^{j,1} \in V_{h,k}^{g_N}$ such that for each $K \in \mathcal{T}_h$, $\mathbf{u}_h^{j,1}$ satisfies

$$\mathbf{u}_h^{j,1} \cdot \mathbf{n}_F = \begin{cases} \{\{\mathbf{u}_h^j \cdot \mathbf{n}_F\}\} & F \in \mathcal{F}_h^{\mathrm{int}}, \\ \mathbf{u}_h^j \cdot \mathbf{n}_F & F \in \mathcal{F}_h^{\mathrm{D}}, \\ g_N & F \in \mathcal{F}_h^{\mathrm{N}} \end{cases} \quad (5.1a)$$

for all faces $F \in \mathcal{F}_K$ and, if the polynomial degree $k \ge 1$,

$$(S^{-1}\mathbf{u}_h^{j,1}, v_h)_K = (S^{-1}\mathbf{u}_h^j, v_h)_K \quad \forall v_h \in \mathcal{P}_{k-1}(K). \quad (5.1b)$$

The flux $\mathbf{u}_h^{j,1}$ is normal-component continuous, belongs to the set $V_{h,k}^{g_N}$, but a priori does not have the correct divergence, $\nabla \cdot \mathbf{u}_h^{j,1} \neq f$ in $\Omega$. If $\Gamma_N = \partial\Omega$, set $p_h^{j,1} := p_h^j - \frac{(p_h^j, 1)}{|\Omega|}$, so that $p_h^{j,1} \in W_{h,k}^*$. Otherwise, set $p_h^{j,1} := p_h^j$.
## 2. Coarse grid solver.

Solve the coarse grid residual problem: Seek for $(\delta_H^{j,2}, r_H^{j,2}) \in V_{H,0}^0 \times W_{H,0}^*$ such that

$$ (S^{-1}\delta_H^{j,2}, v_H) - (r_H^{j,2}, \nabla \cdot v_H) = (p_h^{j,1}, \nabla \cdot v_H) - (S^{-1}\mathbf{u}_h^{j,1}, v_H) \quad \forall v_H \in V_{H,0}^0, \\ (\nabla \cdot \delta_H^{j,2}, q_H) = (f - \nabla \cdot \mathbf{u}_h^{j,1}, q_H) \quad \forall q_H \in W_{H,0}^*. \qquad (5.2) $$

We update

$$ \mathbf{u}_h^{j,2} := \mathbf{u}_h^{j,1} + \delta_H^{j,2} \in V_{h,k}^{gN}, \\ p_h^{j,2} := p_h^{j,1} + r_H^{j,2} \in W_{h,k}^*. \qquad (5.3a) $$

This update now fulfils

$$ (\nabla \cdot \mathbf{u}_h^{j,2}, q_H) = (f, q_H) \quad \forall q_H \in W_{H,0}, \qquad (5.3b) $$

which is a weak divergence constraint on the coarse mesh. Moreover, using Lemma 5.3 below, if additionally $\nabla \cdot \mathbf{u}_h^{j,1} = f$, we get

$$ \|\mathbf{u}_h - \mathbf{u}_h^{j,2}\|^2 = \|\mathbf{u}_h - \mathbf{u}_h^{j,1}\|^2 - \|\delta_H^{j,2}\|^2. $$

## 3. Subdomain Neumann solver

On each subdomain $\Omega_i$, we seek for $(\delta_{h_i}^{j,3}, r_{h_i}^{j,3}) \in V_{i,h,k}^0 \times W_{i,h,k}^*$, such that

$$ (S^{-1}\delta_{h_i}^{j,3}, v_h)_{\Omega_i} - (r_{h_i}^{j,3}, \nabla \cdot v_h)_{\Omega_i} = (p_h^{j,2}, \nabla \cdot v_h)_{\Omega_i} - (S^{-1}\mathbf{u}_h^{j,2}, v_h)_{\Omega_i} \quad \forall v_h \in V_{i,h,k}^0, \\ (\nabla \cdot \delta_{h_i}^{j,3}, q_h)_{\Omega_i} = (f - \nabla \cdot \mathbf{u}_h^{j,2}, q_h)_{\Omega_i} \quad \forall q_h \in W_{i,h,k}^*. \qquad (5.4) $$

Setting $\delta_h^{j,3}|_{\Omega_i} := \delta_{h_i}^{j,3}$ and $r_h^{j,3}|_{\Omega_i} := r_{h_i}^{j,3}$, we update

$$ \mathbf{u}_h^{j,3} := \mathbf{u}_h^{j,2} + \delta_h^{j,3} \in V_{h,k}^{gN}, \\ p_h^{j,3} := p_h^{j,2} + r_h^{j,3} \in W_{h,k}^*. \qquad (5.5a) $$

and observe

$$ \nabla \cdot \mathbf{u}_h^{j,3} = f. \qquad (5.5b) $$

This flux update is thus normal-component continuous and satisfies the divergence constraint (is balanced). Moreover, using Lemma 5.4 below, if additionally $\nabla \cdot \mathbf{u}_h^{j,2} = f$, we get

$$ \|\mathbf{u}_h - \mathbf{u}_h^{j,3}\|^2 = \|\mathbf{u}_h - \mathbf{u}_h^{j,2}\|^2 - \|\delta_h^{j,3}\|^2. $$

## 4. Coarse grid correction

Compute $(\delta_H^{j,4}, r_H^{j,4}) \in V_{H,0}^0 \times W_{H,0}^*$ such that

$$ (S^{-1}\delta_H^{j,4}, v_H) - (r_H^{j,4}, \nabla \cdot v_H) = (p_h^{j,3}, \nabla \cdot v_H) - (S^{-1}\mathbf{u}_h^{j,3}, v_H) \quad \forall v_H \in V_{H,0}^0, \\ (\nabla \cdot \delta_H^{j,4}, q_H) = 0 \quad \forall q_H \in W_{H,0}^*. \qquad (5.6) $$

We finally update

$$ \mathcal{R}_F(\mathbf{u}_h^j, p_h^j) := \mathbf{u}_h^{j,3} + \delta_H^{j,4} \in V_{h,k}^{gN}, \\ \mathcal{R}_P(\mathbf{u}_h^j, p_h^j) := p_h^{j,3} + r_H^{j,4} \in W_{h,k}^*, \qquad (5.7a) $$

where we also denote $\mathcal{R}_{FP}(\mathbf{u}_h^j, p_h^j) := (\mathcal{R}_F(\mathbf{u}_h^j, p_h^j), \mathcal{R}_P(\mathbf{u}_h^j, p_h^j))$, and observe $\nabla \cdot \delta_H^{j,4} = 0$ and thus

$$ \nabla \cdot (\mathcal{R}_F(\mathbf{u}_h^j, p_h^j)) = f. \qquad (5.7b) $$

The flux reconstruction is thus normal-component continuous and satisfies the divergence constraint (is balanced). Moreover, using Lemma 5.3, we get

$$ \|\mathbf{u}_h - \mathcal{R}_F(\mathbf{u}_h^j, p_h^j)\|^2 = \|\mathbf{u}_h - \mathbf{u}_h^{j,3}\|^2 - \|\delta_H^{j,4}\|^2. $$
## 5.2 Correctness of Construction 5.1

Construction 5.1 is correct, as the following lemma establishes:

**Lemma 5.2 (Construction 5.1).** *Prescription (5.1) is well defined and problems (5.2), (5.4), and (5.6) are well-posed, satisfying (5.3), (5.5), and (5.7). In particular, the reconstructed flux $\mathcal{R}_F(\mathbf{u}_h^j, p_h^j)$ belongs to the set $V_{h,k}^{gN}$ with $\nabla \cdot \mathcal{R}_F(\mathbf{u}_h^j, p_h^j) = f$ in $\Omega$ and the reconstructed potential $\mathcal{R}_P(\mathbf{u}_h^j, p_h^j)$ belongs to the space $W_{h,k}^*$.*

*Proof.* 1) The formulas (5.1) are exactly the prescription of the degrees of freedom in the space $V_{h,k}$ from (3.5a), cf., e.g., [7], leading to $\mathbf{u}_h^{j,1} \in V_{h,k}^{gN}$.

2) When $\Gamma_D \neq \emptyset$ the problem (5.2) is well posed, there is no condition to verify. When $\Gamma_D = \emptyset$, the Neumann compatibility condition

$$ (f - \nabla \cdot \mathbf{u}_h^{j,1}, 1) = (f, 1) - \langle \mathbf{u}_h^{j,1} \cdot \mathbf{n}_{\Omega}, 1 \rangle_{\partial \Omega} \stackrel{(5.1a)}{=} (f, 1) - \langle g_{\mathrm{N}}, 1 \rangle_{\partial \Omega} \stackrel{(2.4)}{=} 0 \quad (5.8) $$

holds, which guarantees the existence and uniqueness of the corrector $\delta_H^{j,2}$ from (5.2). Moreover, (5.8) together with the second equation in (5.2) yield

$$ (\nabla \cdot \delta_H^{j,2}, q_H) = (f - \nabla \cdot \mathbf{u}_h^{j,1}, q_H) \quad \forall q_H \in W_{H,0} \quad (5.9) $$

in both cases $\Gamma_D \neq \emptyset$ and $\Gamma_D = \emptyset$.

3) On each subdomain $\Omega_i$ such that $|\partial\Omega \cap \Gamma_D| = 0$, problem (5.4) is a (homogeneous) Neumann problem and one needs to ensure that it satisfies the Neumann compatibility condition, i.e.,

$$ (f - \nabla \cdot \mathbf{u}_h^{j,2}, 1)_{\Omega_i} = 0. \quad (5.10) $$

This is immediate from (5.2)–(5.3) by construction: indeed, (5.9) ensures that, for each subdomain $\Omega_i$, we can take a test function $q_H$ which is one in $\Omega_i$ and zero elsewhere. Consequently, the solution $\delta_{h_i}^{j,3}$ of (5.4) exists and is unique. The reconstructed flux $\mathbf{u}_h^{j,3}$ given by (5.5) belongs to $\mathcal{RT}\mathcal{N}_k(K)$ for all $K \in \mathcal{T}_h$ and the normal fluxes are continuous for all $F \in \mathcal{F}_h^{\mathrm{int}}$ – here the zero normal fluxes of $\delta_{h_i}^{j,3}$ on $\partial\Omega_i$ are crucial. Therefore, $\mathbf{u}_h^{j,3} \in V_{h,k}^{gN}$. Moreover, the second equation in (5.4) together with (5.10) yield

$$ (\nabla \cdot \delta_{h_i}^{j,3}, q_h)_{\Omega_i} = (f - \nabla \cdot \mathbf{u}_h^{j,2}, q_h)_{\Omega_i} \quad \forall q_h \in W_{i,h,k}. \quad (5.11) $$

Thus, from (5.4)–(5.5a), there holds

$$ (\nabla \cdot \mathbf{u}_h^{j,3}, q_h)_{\Omega_i} = (f, q_h)_{\Omega_i} \quad \forall q_h \in W_{i,h,k} $$

for all subdomains $\Omega_i$, i.e., (5.5b).

4) When $\Gamma_D \neq \emptyset$ the problem (5.6) is well posed, there is no condition to verify. When $\Gamma_D = \emptyset$, the problem (5.6) is a (homogeneous) Neumann problem, and the solution $\delta_H^{j,4}$ exists and is unique since both the divergence constraint and the normal trace value are 0. As above, this allows to see that the second equation in (5.6) actually holds true for all $q_H \in W_{H,0}$ and thus $\delta_H^{j,4}$ is divergence-free.

5) Summarizing (5.1a), (5.3), (5.5), (5.7) leads to $\mathcal{R}_F(\mathbf{u}_h^j, p_h^j) \in V_{h,k}^{gN}$ with $\nabla \cdot \mathcal{R}_F(\mathbf{u}_h^j, p_h^j) = f$ and $\mathcal{R}_P(\mathbf{u}_h^j, p_h^j) \in W_{h,k}^*$, as stated in (5.7). $\square$

## 5.3 Error decrease formulas in Construction 5.1

We now show the boxed error decrease formulas in Construction 5.1. We use the equivalences recalled in Section 3.4
**Lemma 5.3** (Error decrease formula for the coarse solve). Let $\mathbf{u}_h \in \mathbf{V}_{h,k}^{gN}$ be the solution of problem (3.11), let $\mathbf{u}_h^{j,1} \in \mathbf{V}_{h,k}^{gN}$ such that $\nabla \cdot \mathbf{u}_h^{j,1} = f$ be arbitrary, and let $\delta_H^{j,2} \in \mathbf{V}_{H,0}^0$ be given by

$$ \delta_H^{j,2} := \arg \min_{\substack{\mathbf{v}_H \in \mathbf{V}_{H,0}^0 \\ \nabla \cdot \mathbf{v}_H = 0}} \|\|\mathbf{u}_h^j + \mathbf{v}_H\|\|^2, \qquad (5.12) $$

or, equivalently, $\delta_H^{j,2} \in \mathbf{V}_{H,0}^0$ such that $\nabla \cdot \delta_H^{j,2} = 0$ and

$$ (S^{-1}\delta_H^{j,2}, \mathbf{v}_H) = -(S^{-1}\mathbf{u}_h^j, \mathbf{v}_H) \quad \forall \mathbf{v}_H \in \mathbf{V}_{H,0}^0 \text{ such that } \nabla \cdot \mathbf{v}_H = 0. \qquad (5.13) $$

Then, we have

$$ \|\|\mathbf{u}_h - (\mathbf{u}_h^{j,1} + \delta_H^{j,2})\|\|^2 = \|\|\mathbf{u}_h - \mathbf{u}_h^{j,1}\|\|^2 - \|\|\delta_H^{j,2}\|\|^2. \qquad (5.14) $$

*Proof.* By choosing $\mathbf{v}_H = \delta_H^{j,2}$ in (5.13), we obtain

$$ \|\|\delta_H^{j,2}\|\|^2 = -(S^{-1}\mathbf{u}_h^{j,1}, \delta_H^{j,2}). \qquad (5.15) $$

Also, since $\delta_H^{j,2} \in \mathbf{V}_{H,0}^0 \subset \mathbf{V}_{h,k}^0$ and $\nabla \cdot \delta_H^{j,2} = 0$, one gets from (3.11)

$$ (S^{-1}\mathbf{u}_h, \delta_H^{j,2}) = 0. \qquad (5.16) $$

Therefore,

$$ \begin{align} \|\|\mathbf{u}_h - (\mathbf{u}_h^{j,1} + \delta_H^{j,2})\|\|^2 &= \|\|\mathbf{u}_h - \mathbf{u}_h^{j,1}\|\|^2 - 2(S^{-1}(\mathbf{u}_h - \mathbf{u}_h^{j,1}), \delta_H^{j,2}) + \|\|\delta_H^{j,2}\|\|^2 \nonumber \\ &\stackrel{(5.16)}{=} \|\|\mathbf{u}_h - \mathbf{u}_h^{j,1}\|\|^2 + 2(S^{-1}\mathbf{u}_h^{j,1}, \delta_H^{j,2}) + \|\|\delta_H^{j,2}\|\|^2 \tag{5.17} \\ &\stackrel{(5.15)}{=} \|\|\mathbf{u}_h - \mathbf{u}_h^{j,1}\|\|^2 - \|\|\delta_H^{j,2}\|\|^2. \nonumber \end{align} $$

□

**Lemma 5.4** (Error decrease formula for the subdomain Neumann solve). Let $\mathbf{u}_h \in \mathbf{V}_{h,k}^{gN}$ be the solution of problem (3.11), let $\mathbf{u}_h^{j,2} \in \mathbf{V}_{h,k}^{gN}$ be such that $\nabla \cdot \mathbf{u}_h^{j,2} = f$, and let $\delta_h^{j,3} \in \mathbf{V}_{h,k}^0$ and $\delta_{h_i}^{j,3} \in \mathbf{V}_{i,h,k}^0$ for each subdomain $\Omega_i$ be given by

$$ \delta_h^{j,3}|_{\Omega_i} := \delta_{h_i}^{j,3} := \arg \min_{\substack{\mathbf{v}_h \in \mathbf{V}_{i,h,k}^0 \\ \nabla \cdot \mathbf{v}_h = 0}} \|\|\mathbf{v}_h + \mathbf{u}_h^{j,2}\|\|^2_{\Omega_i} \qquad (5.18) $$

or, equivalently, $\delta_h^{j,3}|_{\Omega_i} = \delta_{h_i}^{j,3} \in \mathbf{V}_{i,h,k}^0$ such that $\nabla \cdot \delta_{h_i}^{j,3} = 0$ with

$$ (S^{-1}\delta_{h_i}^{j,3}, \mathbf{v}_h)_{\Omega_i} = -(S^{-1}\mathbf{u}_h^{j,2}, \mathbf{v}_h)_{\Omega_i} \quad \forall \mathbf{v}_h \in \mathbf{V}_{i,h,k}^0 \text{ such that } \nabla \cdot \mathbf{v}_h = 0. \qquad (5.19) $$

Then, we have

$$ \|\|\mathbf{u}_h - (\mathbf{u}_h^{j,2} + \delta_{h_i}^{j,3})\|\|^2_{\Omega_i} = \|\|\mathbf{u}_h - \mathbf{u}_h^{j,2}\|\|^2_{\Omega_i} - \|\|\delta_{h_i}^{j,3}\|\|^2_{\Omega_i}, \qquad (5.20a) $$

$$ \|\|\mathbf{u}_h - (\mathbf{u}_h^{j,2} + \delta_{h_i}^{j,3})\|\|^2 = \|\|\mathbf{u}_h - \mathbf{u}_h^{j,2}\|\|^2 - \|\|\delta_{h_i}^{j,3}\|\|^2. \qquad (5.20b) $$

*Proof.* On each subdomain $\Omega_i$, by taking $\mathbf{v}_h = \delta_{h_i}^{j,3}$ in (5.19), we get

$$ \|\|\delta_{h_i}^{j,3}\|\|^2_{\Omega_i} = -(S^{-1}\mathbf{u}_h^{j,2}, \delta_{h_i}^{j,3})_{\Omega_i}. \qquad (5.21) $$

In addition, since $\delta_{h_i}^{j,3} \in \mathbf{V}_{i,h,k}^0$ and $\nabla \cdot \delta_{h_i}^{j,3} = 0$, which we extend by zero outside of $\Omega_i$, one gets from (3.11)

$$ (S^{-1}\mathbf{u}_h, \delta_{h_i}^{j,3})_{\Omega_i} = 0. \qquad (5.22) $$
Therefore,

$$ \begin{aligned} |||\mathbf{u}_h - (\mathbf{u}_h^{j,2} + \delta_{h_i}^{j,3})|||_{\Omega_i}^2 &= |||\mathbf{u}_h - \mathbf{u}_h^{j,2}|||_{\Omega_i}^2 - 2(S^{-1}(\mathbf{u}_h - \mathbf{u}_h^{j,2}), \delta_{h_i}^{j,3})_{\Omega_i} + |||\delta_{h_i}^{j,3}|||_{\Omega_i}^2 \\ &\stackrel{(5.22)}{=} |||\mathbf{u}_h - \mathbf{u}_h^{j,2}|||_{\Omega_i}^2 + 2(S^{-1}\mathbf{u}_h^{j,2}, \delta_{h_i}^{j,3})_{\Omega_i} + |||\delta_{h_i}^{j,3}|||_{\Omega_i}^2 \quad (5.23) \\ &\stackrel{(5.21)}{=} |||\mathbf{u}_h - \mathbf{u}_h^{j,2}|||_{\Omega_i}^2 - |||\delta_{h_i}^{j,3}|||_{\Omega_i}^2. \end{aligned} $$

This is (5.20a); the result (5.20b) is obtained by summing over all subdomains $\Omega_i$. $\square$

# 6 Trace lifting and elementwise postprocessing

Let $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in V_{h,k}^{dc} \times W_{h,k}$ be arbitrary, typically an approximation to the solution of (3.12) obtained on step $j$ of an iterative algebraic domain decomposition solver as that in Section 8 below to which Construction 5.1 has been applied. Trace lifting and elementwise postprocessing will be based on the two following procedures.

## 6.1 Trace lifting

Recall the spaces $\Psi_k(F)$ from (3.7) and $\Psi_{h,k}^{dc}$ from (3.8). We first elementwise construct an approximation to the face pressures in the discontinuous trace space $\Psi_{h,k}^{dc}$.

**Construction 6.1 (Trace lifting).** Let $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in V_{h,k}^{dc} \times W_{h,k}$ be arbitrary. For each element $K \in \mathcal{T}_h$ and each face $F \in \mathcal{F}_K$ such that $F \notin \mathcal{F}_h^{D}$, define $(\lambda_h^j|_K)|_F \in \Psi_k(F)$, and consequently $\lambda_h^j \in \Psi_{h,k}^{dc}$ by

$$ \langle \lambda_h^j, \mathbf{v}_h \cdot \mathbf{n}_K \rangle_F := -(S^{-1}\mathbf{u}_h^j, \mathbf{v}_h)_K + (\mathbf{p}_h^j, \nabla \cdot \mathbf{v}_h)_K \quad (6.1) $$

for all $\mathbf{v}_h \in \mathcal{RT}\mathcal{N}_k(K)$ with $(\mathbf{v}_h \cdot \mathbf{n}_K)|_{F'} = 0$ for all $F' \in \mathcal{F}_K$, $F \neq F'$, and $(\mathbf{v}_h, \mathbf{r}_h)_K = 0$ for all $\mathbf{r}_h \in \mathcal{P}_{k-1}(K)$ if $k \ge 1$. On each $F \in \mathcal{F}_h^{D}$ the function $\lambda_h^j$ is defined as zero, i.e., $\lambda_h^j|_F := 0$.

Note that $\lambda_h^j$ are a priori double valued on all subdomain inter(faces) $F \in \mathcal{F}_h^{\Gamma_{i,j}}$. When, however, $(\mathbf{u}_h^j, \mathbf{p}_h^j)$ coincides with the solution of (3.12), then the definition of $\lambda_h^j$ in (6.1) is univalued on each face $F \in \mathcal{F}_h^{\text{int}}$ and coincides with $\lambda_h$ from (3.13).

## 6.2 Elementwise postprocessing

In this section, we give a recipe for an elementwise postprocessing $\tilde{\mathbf{p}}_h^j$ that is a higher-order piecewise polynomial on the mesh $\mathcal{T}_h$ and that is weakly continuous over the mesh faces. We more precisely construct $\tilde{\mathbf{p}}_h^j$ in some of the available spaces $\mathcal{M}_{h,k} = \Pi_{K \in \mathcal{T}_h} \mathcal{M}_{h,k}(K)$ from [3, 5, 9, 38], recalled in Appendix A below. We highlight that this elementwise postprocessing is independent of the algebraic method that produced $(\mathbf{u}_h^j, \mathbf{p}_h^j)$.

Recall the definition of the average (3.1). Then we define:

**Construction 6.2 (Elementwise postprocessing).** Let $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in V_{h,k}^{dc} \times W_{h,k}$ be arbitrary, and let $\lambda_h^j \in \Psi_{h,k}^{dc}$ be given by Construction 6.1. Then for each $K \in \mathcal{T}_h$, define $\tilde{\mathbf{p}}_h^j \in \mathcal{M}_{h,k}$ by prescribing $\tilde{\mathbf{p}}_h^j|_K \in \mathcal{M}_{h,k}(K)$ via

$$ (\tilde{\mathbf{p}}_h^j, q_h)_K = (\mathbf{p}_h^j, q_h)_K \qquad \forall q_h \in W_k(K), \quad (6.2a) $$

$$ \langle \tilde{\mathbf{p}}_h^j, \mu_h \rangle_F = \begin{cases} \langle \{\lambda_h^j\}, \mu_h \rangle_F & F \in \mathcal{F}_h^{\text{int}}, \\ 0 & F \in \mathcal{F}_h^{\text{D}}, \\ \langle \lambda_h^j, \mu_h \rangle_F & F \in \mathcal{F}_h^{\text{N}} \end{cases} \quad \forall \mu_h \in \Psi_k(F) \text{ and } F \in \mathcal{F}_K. \quad (6.2b) $$
The reconstructed potential $\tilde{p}_h^j$ in particular satisfies

$$ \langle [\tilde{p}_h^j], \mu_h \rangle_F = 0 \quad \text{for all } F \in \mathcal{F}_h^{\text{int}} \text{ and } \mu_h \in \Psi_k(F), \qquad (6.3a) $$

$$ \langle \tilde{p}_h^j, \mu_h \rangle_F = 0 \quad \text{for all } F \in \mathcal{F}_h^D \text{ and } \mu_h \in \Psi_k(F). \qquad (6.3b) $$

We define

$$ \tilde{\mathcal{R}}_P(\mathbf{u}_h^j, \mathbf{p}_h^j) := \tilde{p}_h^j. $$

# 7 A posteriori estimate of the algebraic error

Let $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in L^2(\Omega) \times L^2(\Omega)$ be arbitrary. Notice that here, the approximate solution $\mathbf{u}_h^j$ does not necessarily lie in the space $V_{h,k}^{gN}$ i.e. $\mathbf{u}_h^j$ can be normal-trace discontinuous, but it can also not be piecewise polynomial. Similarly, $\mathbf{p}_h^j$ does not need to satisfy the zero-mean value constraint and be piecewise polynomial. Later in Section 8, though, $(\mathbf{u}_h^j, \mathbf{p}_h^j)$ is still obtained on step j of a domain decomposition solver and belongs to $V_{h,k}^{dc} \times W_{h,k}$.

Recall the S-weighted $L^2$-orthogonal projection of Section 3.3. Our first main result states a guaranteed and fully computable a posteriori estimate of the algebraic error in $(\mathbf{u}_h^j, \mathbf{p}_h^j)$ with respect to $(\mathbf{u}_h, \mathbf{p}_h)$:

**Theorem 7.1 (A posteriori estimate of the algebraic error).** Let $(\mathbf{u}_h, \mathbf{p}_h) \in V_{h,k}^{gN} \times W_{h,k}^*$ be the solution of problem (3.12) and let $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in L^2(\Omega) \times L^2(\Omega)$ be arbitrary. Let $\sigma_h^j \in V_{h,k}^{gN}$ with $\nabla \cdot \sigma_h^j = f$ and $\tilde{p}_h^j \in \mathcal{M}_{h,k}$ satisfying (6.3) be arbitrary. Then, the following bound for the algebraic error holds

$$ \|\mathbf{u}_h - \mathbf{u}_h^j\|^2 \leq (\eta^j)^2 := (\eta_F^j)^2 + (\eta_P^j)^2, \qquad (7.1a) $$

where

$$ \eta_F^j := \|\mathbf{u}_h^j - \sigma_h^j\| \quad \text{and} \quad \eta_P^j := \|\mathbf{u}_h^j + \Pi_k^{\mathcal{RTN}}(S\nabla\tilde{p}_h^j)\|. \qquad (7.1b) $$

*Proof.* Since $\mathbf{u}_h^j$ does not necessarily belong to the space $V_{h,k}^{gN}$ and does not necessarily satisfy the divergence constraint, we construct, for the purpose of the proof, its orthogonal projection to $V_{h,k}^{gN}$, where divergence equal to $f$ is imposed. This is realized by the solution of the following quadratic constrained minimization problem

$$ \mathbf{w}_h := \arg \min_{\substack{\mathbf{v}_h \in V_{h,k}^{gN} \\ \nabla \cdot \mathbf{v}_h = f}} \|\mathbf{u}_h^j - \mathbf{v}_h\|^2. \qquad (7.2) $$

The minimization (7.2) is equivalent to finding $\mathbf{w}_h \in V_{h,k}^{gN}$ with $\nabla \cdot \mathbf{w}_h = f$ in $\Omega$ that satisfies

$$ (S^{-1}(\mathbf{w}_h - \mathbf{u}_h^j), \mathbf{v}_h) = 0 \quad \forall \mathbf{v}_h \in V_{h,k}^0 \text{ with } \nabla \cdot \mathbf{v}_h = 0 \text{ in } \Omega. \qquad (7.3) $$

The fact that $(\mathbf{w}_h - \mathbf{u}_h) \in V_{h,k}^0$ with $\nabla \cdot (\mathbf{w}_h - \mathbf{u}_h) = 0$ together with (7.3) implies that $(S^{-1}(\mathbf{w}_h - \mathbf{u}_h^j), \mathbf{w}_h - \mathbf{u}_h) = 0$. Consequently, using $\mathbf{w}_h$ from (7.2), we obtain the orthogonal decomposition

$$ \begin{aligned} \|\mathbf{u}_h^j - \mathbf{u}_h\|^2 &= \|\mathbf{u}_h^j - \mathbf{w}_h + \mathbf{w}_h - \mathbf{u}_h\|^2 \\ &= \|\mathbf{u}_h^j - \mathbf{w}_h\|^2 + \|\mathbf{w}_h - \mathbf{u}_h\|^2 + 2(S^{-1}(\mathbf{u}_h^j - \mathbf{w}_h), \mathbf{w}_h - \mathbf{u}_h) \quad (7.4) \\ &= \|\mathbf{u}_h^j - \mathbf{w}_h\|^2 + \|\mathbf{w}_h - \mathbf{u}_h\|^2. \end{aligned} $$

Now, we split the analysis into the two parts related to the two terms on the right-hand side of (7.4).
1. **Flux bound.** (First term in (7.4)) Since $w_h$ is the solution to (7.2), for an arbitrary $\sigma_h^j \in V_{h,k}^{gN}$ with $\nabla \cdot \sigma_h^j = f$ in $\Omega$, there holds

$$ \|\mathbf{u}_h^j - \mathbf{w}_h\| \le \|\mathbf{u}_h^j - \boldsymbol{\sigma}_h^j\| = \eta_F^j. $$

2. **Potential bound.** (Second term in (7.4)) Recall that $(w_h - u_h) \in V_{h,k}^0$ with $\nabla \cdot (w_h - u_h) = 0$. Then, the formulations (3.11) and (7.3) lead to

$$ \|w_h - u_h\|^2 = (S^{-1}(w_h - u_h), w_h - u_h) \stackrel{(3.11)}{=} (S^{-1}w_h, w_h - u_h) \stackrel{(7.3)}{=} (S^{-1}\mathbf{u}_h^j, w_h - u_h). \quad (7.5) $$

For any function $\tilde{\mathbf{p}}_h^j \in \mathcal{M}_{h,k}$, using Green's theorem separately on each mesh element $K \in \mathcal{T}_h$, one has, for all $\mathbf{v}_h \in V_{h,k}^0$ with $\nabla \cdot \mathbf{v}_h = 0$ in $\Omega$,

$$ \sum_{K \in \mathcal{T}_h} (\nabla \tilde{\mathbf{p}}_h^j, \mathbf{v}_h)_K = \sum_{K \in \mathcal{T}_h} (\langle \tilde{\mathbf{p}}_h^j, \mathbf{v}_h \cdot \mathbf{n}_K \rangle_{\partial K} - (\tilde{\mathbf{p}}_h^j, \nabla \cdot \mathbf{v}_h)_K) \\ = \sum_{F \in \mathcal{F}_h^{\text{int}}} \langle [\tilde{\mathbf{p}}_h^j], \mathbf{v}_h \cdot \mathbf{n}_F \rangle_F + \sum_{F \in \mathcal{F}_h^{\text{D}}} \langle \tilde{\mathbf{p}}_h^j, \mathbf{v}_h \cdot \mathbf{n}_F \rangle_F. $$

Moreover, if $\tilde{\mathbf{p}}_h^j$ satisfies (6.3) we obtain

$$ \sum_{K \in \mathcal{T}_h} (\nabla \tilde{\mathbf{p}}_h^j, \mathbf{v}_h)_K = 0 \quad (7.6) $$

for all $\mathbf{v}_h \in V_{h,k}^0$ with $\nabla \cdot \mathbf{v}_h = 0$ in $\Omega$. In particular, $(\nabla \tilde{\mathbf{p}}_h^j, w_h - u_h) = 0$. This, together with (7.5) and the definition of the projector $\Pi_k^{\mathcal{RTN}}$ from (3.9) leads to

$$ \|w_h - u_h\|^2 \stackrel{(7.5),(7.6)}{=} (S^{-1}\mathbf{u}_h^j + \nabla \tilde{\mathbf{p}}_h^j, w_h - u_h) = (S^{-1}(\mathbf{u}_h^j + \Pi_k^{\mathcal{RTN}}(S\nabla \tilde{\mathbf{p}}_h^j)), w_h - u_h). \quad (7.7) $$

Applying the Cauchy-Schwarz inequality finally yields

$$ \|w_h - u_h\| \le \|\mathbf{u}_h^j + \Pi_k^{\mathcal{RTN}}(S\nabla \tilde{\mathbf{p}}_h^j)\| = \eta_P^j. $$

□

**Remark 7.2** (Possible construction of $\boldsymbol{\sigma}_h^j$ and $\tilde{\mathbf{p}}_h^j$). A possible construction, when $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in V_{h,k}^{dc} \times W_{h,k}$, is $\boldsymbol{\sigma}_h^j := \mathcal{R}_F(\mathbf{u}_h^j, \mathbf{p}_h^j)$ from Construction 5.1 and then $\tilde{\mathbf{p}}_h^j = \tilde{\mathcal{R}}_P(\boldsymbol{\sigma}_h^j, \mathcal{R}_P(\mathbf{u}_h^j, \mathbf{p}_h^j))$ from Construction 6.2.

**Remark 7.3** (Consistency check). Note that if $\mathbf{u}_h^j = \mathbf{u}_h$ and $\mathbf{p}_h^j = p_h$, the estimator $\eta^j = 0$ for the choices $\boldsymbol{\sigma}_h^j = \mathbf{u}_h^j$ and $\tilde{\mathbf{p}}_h^j$ obtained from Construction 6.2. Indeed, on the one hand, if $\mathbf{u}_h^j \in V_{h,k}^{gN}$ with $\nabla \cdot \mathbf{u}_h^j = f$ in $\Omega$, then we can choose $\boldsymbol{\sigma}_h^j = \mathbf{u}_h^j$ to obtain $\eta_F^j = 0$. On the other hand, when $(\mathbf{u}_h^j, \mathbf{p}_h^j) \in V_{h,k}^{gN} \times W_{h,k}^*$ satisfies (3.12), the corresponding $\tilde{\mathbf{p}}_h^j \in \mathcal{M}_{h,k}$ is such that (6.3) holds. For this choice in (7.1), we conclude that $\eta_P^j = 0$.

## 8 Nonoverlapping domain decomposition algorithm with local mass conservation on each step

Recall the definition of the weighted average (3.2). We can now present our second main result, a domain decomposition algorithm producing a physically correct (locally mass conservative) approximation on each step:
**Algorithm 8.1** (Nonoverlapping domain decomposition algorithm with local mass conservation on each step). Let $\epsilon > 0$ be a user-specified tolerance. Let an arbitrary initial guess $(\mathbf{u}_h^0, p_h^0) \in V_{h,k}^{\text{dc}} \times W_{h,k}$ be given; recall that this means that the initial flux $\mathbf{u}_h^0$ can be normal-trace discontinuous with a wrong divergence and the initial potential $p_h^0$ can be non mean value-free if $\Gamma_N = \partial\Omega$. Proceed in initialization and four steps:

0. Initialization by flux equilibration.

Compute the reconstructed flux and potential using $(\mathbf{u}_h^1, p_h^1) := \mathcal{R}_{\text{FP}}(\mathbf{u}_h^0, p_h^0)$ from Construction 5.1. This gives a starting flux $\mathbf{u}_h^1$ which is normal-trace continuous and with the correct divergence and a starting potential which is of mean value zero if $\Gamma_N = \partial\Omega$,

$$ \begin{aligned} & \mathbf{u}_h^1 \in V_{h,k}^{gN}, \quad \nabla \cdot \mathbf{u}_h^1 = f, \\ & p_h^1 \in W_{h,k}^*. \end{aligned} \tag{8.1} $$

All the following iterations will retain these properties.

Set $j = 1$ and perform the following four steps:

1. Trace lifting.

From the solution $(\mathbf{u}_h^j, p_h^j)$, compute the associated intermediate Lagrange multiplier $\lambda_h^j \in \Psi_{h,k}^{\text{dc}}$ by Construction 6.1.

2. Subdomain Dirichlet solver.

Construct the intermediate flux $\delta_h^j \in V_{h,k}^{\text{dc}}$ and the intermediate potential $r_h^j \in W_{h,k}$ such that $\delta_h^j|_{\Omega_i} := \delta_{h_i}^j$ and $r_h^j|_{\Omega_i} := r_{h_i}^j$ where for all subdomains $\Omega_i$, $(\delta_{h_i}^j, r_{h_i}^j) \in V_{i,h,k} \times W_{i,h,k}$ is the solution of the subdomain Dirichlet problems

$$ \begin{aligned} (S^{-1}\delta_{h_i}^j, v_h)_{\Omega_i} - (r_{h_i}^j, \nabla \cdot v_h)_{\Omega_i} &= (p_h^j, \nabla \cdot v_h)_{\Omega_i} - (S^{-1}\mathbf{u}_h^j, v_h)_{\Omega_i} - \langle\langle\{\lambda_h^j\}\rangle\rangle_{\text{w}}, \mathbf{v}_h \cdot \mathbf{n}\rangle_{\partial\Omega_i} \partial\Omega \\ &\quad - \langle\lambda_h^j, \mathbf{v}_h \cdot \mathbf{n}\rangle_{\partial\Omega_i \cap \partial\Omega} \quad \forall \mathbf{v}_h \in V_{i,h,k}, \\ (\nabla \cdot \delta_{h_i}^j, q_{h_i})_{\Omega_i} &= 0 \quad \forall q_{h_i} \in W_{i,h,k}. \end{aligned} $$

3. Flux equilibration and line search.

From Steps 1 and 2, $\mathbf{u}_h^j + \delta_h^j$ would not have continuous normal component on the subdomain interfaces, which we now correct by Construction 5.1. We in particular first define the flux and potential as $(\hat{\mathbf{u}}_h^j, \hat{p}_h^j) := \mathcal{R}_{\text{FP}}(\hat{\mathbf{u}}_h^j + \delta_h^j, p_h^j + r_h^j)$ from Construction 5.1. Thus $\hat{\mathbf{u}}_h^j \in V_{h,k}^{gN}$ with $\nabla \cdot \hat{\mathbf{u}}_h^j = f$ and $\hat{p}_h^j \in W_{h,k}^*$ and $\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j \in V_{h,k}^{gN}$ with $\nabla \cdot (\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) = 0$ and $\hat{p}_h^j - p_h^j \in W_{h,k}^*$. We then finally update

$$ \begin{aligned} \mathbf{u}_h^{j+1} &:= \mathbf{u}_h^j + \alpha^j(\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j), \\ p_h^{j+1} &:= p_h^j + \alpha^j(\hat{p}_h^j - p_h^j), \end{aligned} \tag{8.2} $$

where

$$ \alpha^j := - \frac{(S^{-1}\mathbf{u}_h^j, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j)}{\|\|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|\|} \tag{8.3} $$

is the optimal step-length parameter obtained by line search. This in particular yields $\mathbf{u}_h^{j+1}$ which is normal-trace continuous and with the correct divergence (mass conservative) and $p_h^{j+1}$ which is of zero-mean value when $\Gamma_N = \partial\Omega$,

$$ \begin{aligned} & \mathbf{u}_h^{j+1} \in V_{h,k}^{gN}, \quad \nabla \cdot \mathbf{u}_h^{j+1} = f, \\ & p_h^{j+1} \in W_{h,k}^*. \end{aligned} \tag{8.4} $$
### 4. Algebraic error estimate and stopping criterion.

Estimate the algebraic error from below by

$$ (\underline{\eta}^j)^2 := \frac{(\mathbf{S}^{-1}\mathbf{u}_h^j, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j)^2}{\|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2} = (\alpha^j)^2 \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2. \quad (8.5) $$

If $\eta^j < \epsilon$, then stop the solver. Otherwise go to Step 1.

Thanks to the line search giving the optimal step size (8.3), we have the following link of the two consecutive errors:

**Theorem 8.2 (Error decrease formula).** There holds

$$ \|\mathbf{u}_h - \mathbf{u}_h^{j+1}\|^2 = \|\mathbf{u}_h - \mathbf{u}_h^j\|^2 - \frac{(\mathbf{S}^{-1}\mathbf{u}_h^j, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j)^2}{\|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2}. \quad (8.6) $$

*Proof.* This follows from the definition of $\alpha^j$ in (8.3). Indeed, developing from the first equation in (8.2),

$$ \begin{aligned} \|\mathbf{u}_h - \mathbf{u}_h^{j+1}\|^2 &= \|\mathbf{u}_h - (\mathbf{u}_h^j + \alpha^j(\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j))\|^2 \\ &= \|\mathbf{u}_h - \mathbf{u}_h^j\|^2 - 2\alpha^j (\mathbf{S}^{-1}(\mathbf{u}_h - \mathbf{u}_h^j), \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) + (\alpha^j)^2 \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2 \quad (8.7) \\ &= \|\mathbf{u}_h - \mathbf{u}_h^j\|^2 + 2\alpha^j (\mathbf{S}^{-1}\mathbf{u}_h^j, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) + (\alpha^j)^2 \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2, \end{aligned} $$

where the last step follows since $\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j \in V_{h,k}^0$ and $\nabla \cdot (\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) = 0$, so that (3.11) gives

$$ (\mathbf{S}^{-1}\mathbf{u}_h, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) = 0. \quad (8.8) $$

Therefore, inserting formula (8.3) into (8.7), we get the error decrease (8.6). Actually, (8.7) shows that $\alpha^j$ given by formula (8.3) gives the largest error decrease between all real parameters $\alpha$: (8.7) is a quadratic function in this parameter, so that the optimal step-size $\alpha^j$ is such that the derivative of this quadratic function vanishes,

$$ 2(\mathbf{S}^{-1}\mathbf{u}_h^j, \hat{\mathbf{u}}_h^j - \mathbf{u}_h^j) + 2\alpha^j \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2 = 0. \quad (8.9) $$

□

Combining Theorems 7.1 and 8.2, we also fully control the algebraic error by computable quantities:

**Theorem 8.3 (Guaranteed upper and lower bound on the algebraic error).** Let $\underline{\eta}^j$ be given by (8.5) and let $\eta^j$ be given by $\eta^j := \|\mathbf{u}_h^j + \mathbf{\Pi}_k^{\mathcal{RTN}}(\mathbf{S}\nabla\tilde{\mathbf{p}}_h^{j+1})\|$ with $\tilde{\mathbf{p}}_h^{j+1} := \tilde{\mathcal{R}}_P(\mathbf{u}_h^{j+1}, \mathbf{p}_h^{j+1})$ from Construction 6.2. Then, for $j \ge 1$, there holds

$$ \underline{\eta}^j \le \|\mathbf{u}_h - \mathbf{u}_h^j\| \le \eta^j. $$

Moreover, the a posteriori error estimators admit the following elementwise structure:

$$ (\underline{\eta}^j)^2 = (\alpha^j)^2 \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|^2 = \sum_{K \in \mathcal{T}_h} (\alpha^j)^2 \|\hat{\mathbf{u}}_h^j - \mathbf{u}_h^j\|_K^2, \quad (8.10) $$

$$ (\eta^j)^2 = \|\mathbf{u}_h^j + \mathbf{\Pi}_k^{\mathcal{RTN}}(\mathbf{S}\nabla\tilde{\mathbf{p}}_h^{j+1})\|^2 = \sum_{K \in \mathcal{T}_h} \|\mathbf{u}_h^j + \mathbf{\Pi}_k^{\mathcal{RTN}}(\mathbf{S}\nabla\tilde{\mathbf{p}}_h^{j+1})\|_K^2. \quad (8.11) $$
Proof. Using Theorem 7.1, where we choose $\sigma_h^j = u_h^j$ (this is possible since $u_h^j \in V_{h,k}^{\text{GN}}$ with $\nabla \cdot u_h^j = f$ for $j \ge 1$), we get the upper bound. Theorem 8.2, in turn, immediately gives the lower bound employing the definition $\underline{\eta}^j$ by (8.5). $\square$

Finally, similarly to [24, Theorem 5.3], we have:

**Theorem 8.4 (Equivalence of error contraction and efficiency of the lower a posteriori algebraic error estimate).** For $j \ge 1$, if $(u_h^{j+1}, p_h^{j+1})$ is constructed using one step of the iterative Algorithm 8.1, then there holds

$$ \|\|u_h - u_h^{j+1}\|\| \le C \|\|u_h - u_h^j\|\| \quad (8.12) $$

for some $0 \le C < 1$ if and only if

$$ \sqrt{(1-C^2)} \|\|u_h - u_h^j\|\| \le (\underline{\eta}^j). \quad (8.13) $$

*Proof.* By Theorem 8.2, we have

$$ \|\|u_h - u_h^{j+1}\|\|^2 = \|\|u_h - u_h^j\|\|^2 - (\underline{\eta}^j)^2. \quad (8.14) $$

Then

$$ \begin{align*} \|\|u_h - u_h^{j+1}\|\|^2 &\le C^2 \|\|u_h - u_h^j\|\|^2 \\ \stackrel{(8.14)}{\rightleftharpoons} &\|\|u_h - u_h^j\|\|^2 - (\underline{\eta}^j)^2 \le C^2 \|\|u_h - u_h^j\|\|^2 \\ \rightleftharpoons & (1-C^2) \|\|u_h - u_h^j\|\|^2 \le (\underline{\eta}^j)^2. \end{align*} $$

$\square$

**Remark 8.5 (Round-off errors).** Above, we have assumed an idealistic exact arithmetics without round-off errors. In practice, round-off errors, particularly when considering high ratio contrasts in the coefficients, can have an (important) influence on computational accuracy. A full analysis of this aspect is beyond the scope of the present study, and it is left as a consideration for future work.

**Remark 8.6 (Other models).** The addition of a zero-th order term or the introduction of time discretization in the Darcy flow problem (2.1) does not modify the structure of Algorithm 8.1. Incorporating a zero-order term modifies the variational formulation, affecting the calculation of a posteriori errors, the determination of the parameter $\alpha^j$ in the line search, and the calculation of the intermediate Lagrange multiplier $\lambda_h^j$ – the underlying scalar product needs to be changed. However, the basic properties of the method remain unchanged. Similarly, introducing time discretization allows the method to naturally extend to the time-dependent framework.

**Remark 8.7 (Pressure estimates).** Our theoretical results are limited to the velocities $u_h^j$, and we do not have corresponding theoretical results for the errors in the pressures $p_h^j$. These could be obtained but would not be constant-free (the inf-sup condition comes into play).

# 9 Numerical experiments

In this section, we numerically illustrate our a posteriori error estimates on the algebraic error of Theorem 8.3 and the nonoverlapping domain decomposition solver of Algorithm 8.1. We consider problem (2.5) on the domain $\Omega$ being the unit square $(0,1)^2$ in two situations: $S$ an identity diffusion tensor and $\mathcal{S}$ a piecewise constant diffusion tensor with a varying contrast. We consider $d=2$, triangular meshes, and the mixed finite element discretization (3.12) with polynomial degree $k=0$. For the spaces $M_{h,0}(K)$, we use (A.14) from Appendix A.2.3 below. The numerical experiments are performed with FreeFem++ https://freefem.org/, see Hecht [21].
<figure><img src="image_2.png" /><figcaption>Figure 1: [Test 1] Left: the coarse mesh $\mathcal{T}_H$ with 8 triangular subdomains $\Omega_i$. Right: the fine mesh $\mathcal{T}_h$ with 8192 triangles.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: [Test 1] Left: exact mixed finite element flux $u_h$ from (3.12). Right: exact mixed finite element potential $p_h$ from (3.12).</figcaption></figure>

## 9.1 Test 1: identity diffusion tensor

In this section, we suppose the diffusion tensor $S = I$ and problem (2.5) with the exact solution $u(x, y) := ((1 - 2x)(y^2 - y), (1 - 2y)(x^2 - x))$, corresponding to the source term $f(x, y) = -2(x^2 + y^2) + 2(x + y)$ and $\Gamma_D = \partial\Omega$; the exact potential writes as $p(x, y) = x(x - 1)y(y - 1)$.

To start with, we decompose $\Omega$ into 8 nonoverlapping triangular subdomains $\Omega_i$, $i = 1, \dots, 8$, to create the coarse mesh $\mathcal{T}_H$. We then refine each subdomain $\Omega_i$ individually. This gives the refined mesh $\mathcal{T}_h$ with $\mathcal{T}_{h_i} = \mathcal{T}_h \cap \Omega_i$. These coarse and fine meshes are illustrated in Figure 1; the number of unknowns in the linear system (3.12) is 20608. Later, we also refine $\mathcal{T}_H$. The exact solution of the mixed finite element discretization (3.12) is illustrated in Figure 2 and we consider Algorithm 8.1 with the initial guess $(u_h^0, p_h^0) = (0, 0)$ to approximate it.

### 9.1.1 The initialization step 0 of Algorithm 8.1

In the initialization step 0 of Algorithm 8.1, we apply Construction 5.1 to the potential $p_h^0$ and the flux $u_h^0$ in order to obtain a normal-trace continuous flux $u_h^1$ in $V_{h,k}$ with divergence equal to $f$ and the corresponding potential $p_h^1$ in $W_{h,k}$ as per (5.7). In Figures 3, 4, and 5 we illustrate the corrections made in Construction 5.1 by respectively the coarse-grid solver of step 2, the subdomain Neumann solver of step 3, and the coarse grid correction of step 4. We observe a
<figure><img src="image_4.png" /><figcaption>Figure 3: [Test 1] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $\delta_H^{0,2}$. Center: $\delta_h^{0,3}$. Right: $\delta_H^{0,4}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: [Test 1] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $u_h^{0,2} = \delta_H^{0,2}$. Center: $u_h^{0,3} = u_h^{0,2} + \delta_h^{0,3}$. Right: $\mathcal{R}_F(u_h^0, p_h^0) = u_h^1 = u_h^{0,3} + \delta_H^{0,4}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: [Test 1] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $p_h^{0,2} = r_H^{0,2}$. Center: $p_h^{0,3} = p_h^{0,2} + r_h^{0,3}$. Right: $\mathcal{R}_P(u_h^0, p_h^0) = p_h^1 = p_h^{0,3} + r_H^{0,4}$.</figcaption></figure>

nice progress from the initial zero guess to the target approximate solution of Figure 2 in all steps; the coarse grid correction of step 4 seems to primarily affect the potential.

### 9.1.2 Convergence

In Figure 6, we trace the algebraic error $\|u_h - u_h^j\|$ and the contraction factor defined as

$$ \frac{\|u_h - u_h^{j+1}\|}{\|u_h - u_h^j\|} \qquad (9.1) $$

as a function of the iteration index $j$. We see a substantial reduction of the algebraic error on each step of Algorithm 8.1 and a fast convergence, which validates the effectiveness of our domain decomposition algorithm.
<figure><img src="image_7.png" /><figcaption>Figure 6: [Test 1] Algebraic error $||u_h - u_h^j||$ (blue solid line, left scale). Contraction factor given by (9.1) (red dashed line, right scale).</figcaption></figure>

### 9.1.3 A posteriori estimates of the algebraic error

We now investigate the quality of the a posteriori error estimates of Theorem 8.3. In Figure 7, we report the lower bound error estimator $\eta^j$, the error $||u_h - u_h^j||$, and the upper bound error estimator $\eta^j$, again in function of the iteration index $j$. We observe a close match, especially for the lower bound. The quality of our a posteriori error estimators is then more closely examined in the subsequent Figure 8, where we report the lower bound and the upper bound effectivity indices given respectively by

$$ \underline{I}_{\text{eff}} := \frac{||u_h - u_h^j||}{\eta^j} \geq 1 \text{ and } \underline{I}_{\text{eff}} := \frac{\eta^j}{||u_h - u_h^j||} \geq 1. \qquad (9.2) $$

Especially $\underline{I}_{\text{eff}}$ is very close to the optimal value of 1 and confirms the sharpness of the estimates.

<figure><img src="image_8.png" /><figcaption>Figure 7: [Test 1] Algebraic error $||u_h - u_h^j||$ (blue solid line), upper bound $\eta^j$ (red dotted line), and lower bound $\eta^j$ (green dashed line).</figcaption></figure>

In Figure 9, we then plot the spatial distribution of the error $||u_h - u_h^1||$, the algebraic error lower estimators $\alpha^1||\hat{u}_h^1 - u_h^1||_K$, and the elementwise algebraic error upper estimators $||u_h^1 + \Pi_k^{RTN}(\nabla \tilde{p}_h^2)||_K$ on the 1st iteration, $j = 1$. We see a very close prediction of the distribution of the algebraic error, namely for the lower estimators. In Figure 10, the same results are plotted for iteration $j = 14$, again with a very close match. For illustration, we also plot the elementwise errors $||p_h - p_h^1||_K$ and $||p_h - p_h^{14}||_K$ in Figure 11.
<figure><img src="image_9.png" /><figcaption>Figure 8: [Test 1] Effectivity indices $I_{\text{eff}}$ (green solid line) and $I_{\text{eff}}$ (red dashed line) from (9.2).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: [Test 1] Left: elementwise 1st iteration flux errors $\|u_h - u_h^1\|_K$. Center: elementwise 1st iteration algebraic error lower estimators $\alpha^1 \|\hat{u}_h^1 - u_h^1\|_K$. Right: elementwise 1st iteration algebraic error upper estimators $\|u_h^1 + \Pi_k^{\mathcal{RTN}}(\nabla \tilde{p}_h^2)\|_K$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: [Test 1] Left: elementwise final flux errors $\|u_h - u_h^{14}\|_K$. Center: elementwise final algebraic error lower estimators $\alpha^{14} \|\hat{u}_h^{14} - u_h^{14}\|_K$. Right: elementwise final algebraic error upper estimators $\|u_h^{14} + \Pi_k^{\mathcal{RTN}}(\nabla \tilde{p}_h^{15})\|_K$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: [Test 1] Left: elementwise 1st iteration potential errors $||p_h - p_h^1||_K$. Right: elementwise final potential errors $||p_h - p_h^{14}||_K$.</figcaption></figure>

### 9.1.4 Behavior of Algorithm 8.1 for refinements of the coarse mesh $\mathcal{T}_H$ and of the fine mesh $\mathcal{T}_h$

We now perform a (laptop) scalability test of Algorithm 8.1: we fix the coarse mesh $\mathcal{T}_H$ and refine the fine mesh $\mathcal{T}_h$, fix the fine mesh $\mathcal{T}_h$ and refine the coarse mesh $\mathcal{T}_H$, and refine simultaneously the coarse mesh $\mathcal{T}_H$ and the fine mesh $\mathcal{T}_h$. In this section, we actually replace the calculation of $||u_h - u_h^j||$ by $||u - u_h^j||$, so that we avoid the “exact” solution of (3.12) for the evaluation of the algebraic error (the committed error is negligible).

#### Refining the fine mesh for a fixed coarse mesh

Here we fix the coarse mesh $\mathcal{T}_H$ to 512 triangular subdomains and consider the fine meshes $\mathcal{T}_h$ with respectively 32768, 294912, and 819200 triangles (82176, 738048, and 2049280 unknowns). We report the results in Figure 12. Similar convergence rates are achieved, though the contraction factor slightly increases with the fineness of the mesh $\mathcal{T}_h$.

<figure><img src="image_13.png" /><figcaption>Figure 12: [Test 1] Left: approximate algebraic error $||u - u_h^j||$. Right: the corresponding contraction factor given by (9.1). Coarse mesh $\mathcal{T}_H$ of 512 triangular subdomains fixed, fine meshes $\mathcal{T}_h$ with respectively 32768 (green dashed line), 294912 (blue solid line), and 819200 (red dotted line) triangles (respectively 82176, 738048, and 2049280 unknowns).</figcaption></figure>
## Refining the coarse mesh for a fixed fine mesh

Now we fix the fine mesh $\mathcal{T}_h$ to 819200 triangles (2049280 unknowns) and vary the coarse mesh $\mathcal{T}_H$ to respectively consist of 512, 2048, and 4608 triangular subdomains. In Figure 13, we observe that as we increase the number of subdomains, the convergence rates get slightly better.

<figure><img src="image_14.png" /><figcaption>Figure 13: [Test 1] Left: approximate algebraic error $\|u - u_h^j\|$. Right: the corresponding contraction factor given by (9.1). Fine mesh $\mathcal{T}_h$ of 819200 triangles (2049280 unknowns) fixed, coarse meshes $\mathcal{T}_H$ with respectively 512 (green dashed line), 2048 (blue solid line), and 4608 (red dotted line) triangular subdomains.</figcaption></figure>

## Refining simultaneously the fine and the coarse meshes

We suppose finally that we have three coarse meshes $\mathcal{T}_H$ with respectively 648, 5832, and 10368 triangular subdomains and associated three fine meshes $\mathcal{T}_h$ with respectively 209952, 472392, and 839808 triangles (respectively 525528, 1181952, and 2100816 unknowns). In Figure 14, the convergence rates (contraction factors) seem to be rather stable in this setting.

<figure><img src="image_15.png" /><figcaption>Figure 14: [Test 1] Left: approximate algebraic error $\|u - u_h^j\|$. Right: the corresponding contraction factor given by (9.1). Coarse meshes $\mathcal{T}_H$ and fine meshes $\mathcal{T}_h$ with respectively 648 triangular subdomains and 209952 triangles (green dashed line), 5832 triangular subdomains and 472392 triangles (blue solid line), and 10368 triangular subdomains and 839808 triangles (red dotted line) (respectively 525528, 1181952, and 2100816 unknowns).</figcaption></figure>
### 9.1.5 Behavior of Algorithm 8.1 for non-Cartesian coarse and fine meshes

In this test, we use non-Cartesian coarse and fine meshes. We still use the same data as in [Test 1]. The coarse and fine meshes are illustrated in Figure 15; the number of unknowns in the liner system (3.12) is 37650.

<figure><img src="image_16.png" /><figcaption>Figure 15: [Test 1, unstructured meshes] Left: the coarse mesh $\mathcal{T}_H$ with 6 triangular subdomains $\Omega_i$. Right: the fine mesh $\mathcal{T}_h$ with 15000 triangles.</figcaption></figure>

#### Convergence

In Figure 16, we trace the algebraic error $\|u_h - u_h^j\|$ and the contraction factor defined by (9.1). A comparison between Figure 16 and Figure 6 shows a similar behavior. Although the convergence is not as fast as in Figure 6, it remains uniform, confirming the effectiveness of our domain decomposition algorithm.

<figure><img src="image_17.png" /><figcaption>Figure 16: [Test 1, unstructured meshes] Algebraic error $\|u_h - u_h^j\|$ (blue solid line, left scale). Contraction factor given by (9.1) (red dashed line, right scale).</figcaption></figure>

#### A posteriori estimates of the algebraic error

Similarly to Figure 7, Figure 17 presents the lower bound error estimator $\eta^j$, the error $\|u_h - u_h^j\|$, and the upper bound error estimator $\eta^i$, all as functions of the iteration index $j$. We observe a close correspondence, particularly for the lower bound. The quality of our a posteriori error
estimators is illustrated more closely in Figure 18, where we show the effectivity indices of both the lower and upper bounds, as defined in (9.2). Here, we observe the influence of non-Cartesian subdomains on the error and a posteriori error estimators. Although a close correspondence is still observed, the Cartesian case, shown in Figures 7 and 8, yields better results.

<figure><img src="image_18.png" /><figcaption>Figure 17: [Test 1, unstructured meshes] Algebraic error $||u_h - u_h^j||$ (blue solid line), upper bound $\eta^j$ (red dotted line), and lower bound $\underline{\eta}^j$ (green dashed line).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: [Test 1, unstructured meshes] Effectivity indices $I_{\text{eff}}$ (green solid line) and $I_{\text{eff}}$ (red dashed line) from (9.2).</figcaption></figure>

## 9.2 Test 2: piecewise constant diffusion tensor aligned with the coarse mesh

We now consider a second test case with a piecewise constant diffusion tensor $S = c(x, y)I$, where the variations of $c(x, y)$ are illustrated in Figure 19. We first consider the values of $c(x, y)$ equal to 1 and $10^7$ but later we let this contrast vary in a robustness study. A similar test case is used in Anciaux-Sedrakian et al. [2, Section 5.3]. We set the source term $f = 1$ and the Dirichlet boundary condition $p(x, y) = 0$ on $\Gamma_D = \partial\Omega$. The coarse mesh $\mathcal{T}_H$ of 32 triangular subdomains of Figure 19 together with a fine mesh of 12800 triangles (32160 unknowns) will be used for all the numerical simulations below. The zero initial guess $(u_h^0, p_h^0) = (0, 0)$ will still be employed in Algorithm 8.1.
<figure><img src="image_20.png" /><figcaption>Figure 19: [Test 2] Variations of the coefficient $c(x, y)$ across the domain.</figcaption></figure>

### 9.2.1 Convergence for a fixed diffusion contrast

<figure><img src="image_21.png" /><figcaption>Figure 20: [Test 2] Left: exact mixed finite element flux $u_h$ from (3.12). Right: exact mixed finite element potential $p_h$ from (3.12).</figcaption></figure>

In Figure 20, we present the exact solution ($u_h, p_h$) of the mixed finite element discretization (3.12). We observe a formation of very narrow and strong Darcy flow fields in the subdomain tips caused by the huge contrast $1 : 10^7$.

In the initialization step 0 of Algorithm 8.1, the use of Construction 5.1 yields a normal-trace continuous flux $u_h^1$ lying in $V_{h,k}^{gN}$, with divergence equal to $f$, together with the corresponding potential $p_h^1$ in $W_{h,k}^*$. In Figures 21, 22, and 23, we illustrate the corrections made in Construction 5.1 by respectively the coarse-grid solver of step 2, the subdomain Neumann solver of step 3, and the coarse grid correction of step 4. As in Section 9.1.1, we observe a nice progress from the initial zero guess to the target approximate solution of Figure 20 in all steps.

In Figure 24, we trace the algebraic error $|||u_h - u_h^j|||$ and the corresponding contraction factor defined as in (9.1) (using the energy error). We see a substantial reduction of the algebraic error on each step of Algorithm 8.1 (except for the initialization) and a fast convergence, which validates the effectiveness of our algorithm.

In Figure 25, we report the lower bound error estimator $\eta^j$, the error $|||u_h - u_h^j|||$, and the upper bound error estimator $\eta^j$ of Theorem 8.3. In the consecutive Figure 26, we plot the lower bound effectivity index $I_{\text{eff}}$ as well as the upper bound effectivity index $I_{\text{eff}}$, defined as in (9.2). We observe very good results, even for the huge diffusion contrast $1 : 10^7$.

In Figure 27, we then plot the spatial distribution of the error $|||u_h - u_h^1|||$, the algebraic error lower estimators $\alpha^1|||\hat{u}_h^1 - u_h^1|||_K$, and the elementwise algebraic error upper estimators $|||u_h^1 + \Pi_k^{RTN}(S\nabla\tilde{p}_h^2)|||_K$ on the 1st iteration, $j = 1$. We again observe a very close prediction of the distribution of the algebraic error, even for the huge diffusion contrast $1 : 10^7$. In Figure 28,
<figure><img src="image_22.png" /><figcaption>Figure 21: [Test 2] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $\boldsymbol{\delta}_H^{0,2}$. Center: $\delta_h^{0,3}$. Right: $\delta_H^{0,4}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: [Test 2] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $\boldsymbol{u}_h^{0,2} = \boldsymbol{\delta}_H^{0,2}$. Center: $\boldsymbol{u}_h^{0,3} = \boldsymbol{u}_h^{0,2} + \boldsymbol{\delta}_h^{0,3}$. Right: $\mathcal{R}_F(\boldsymbol{u}_h^0, p_h^0) = \boldsymbol{u}_h^1 = \boldsymbol{u}_h^{0,3} + \boldsymbol{\delta}_H^{0,4}$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 23: [Test 2] Construction 5.1 in the initialization step 0 of Algorithm 8.1. Left: $p_h^{0,2} = r_H^{0,2}$. Center: $p_h^{0,3} = p_h^{0,2} + r_h^{0,3}$. Right: $\mathcal{R}_P(\boldsymbol{u}_h^0, p_h^0) = p_h^1 = p_h^{0,3} + r_H^{0,4}$.</figcaption></figure>

the same results are plotted for iteration $j = 14$, again with a very close match. For illustration, we also plot the elementwise errors $\|p_h - p_h^1\|_K$ and $\|p_h - p_h^{14}\|_K$ in Figure 29.

### 9.2.2 Robustness with respect to a varying diffusion contrast

We now let $c(x, y)$ be still piecewise constant as in Figure 19, but we let the diffusion contrast vary from $10^1$ to $10^8$. Table 1 shows the iteration numbers needed to reduce the initial algebraic error estimator $\eta^1$ by the factor $10^5$, i.e., five orders of magnitude. The number of iterations remains constant (or even decreases with increasing the diffusion contrast), clearly indicating that Algorithm 8.1 is in this case robust with respect to the diffusion contrast. This is an excellent property for a domain decomposition solver for applications in porous media.
<figure><img src="image_25.png" /><figcaption>Figure 24: [Test 2] Algebraic error $|||u_h - u_h^j|||$ (blue solid line, left scale). Contraction factor $|||u_h - u_h^{j+1}|||/|||u_h - u_h^j|||$ (red dashed line, right scale).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 25: [Test 2] Algebraic error $|||u_h - u_h^j|||$ (blue solid line), upper bound $\eta^j$ (red dotted line), and lower bound $\bar{\eta}^j$ (green dashed line).</figcaption></figure>

## 9.3 Test 3: piecewise constant diffusion tensor not aligned with the coarse mesh

In this test, we continue similarly as in the above Test 2, inspired by [2, 25]. The aim is to study two situations, the first where the discontinuities in $c(x, y)$ do not correspond to the coarse mesh, and the second where the interfaces of the coarse mesh cross the discontinuities. Figure 30 illustrates these two settings; we let $c(x, y)$ vary respectively in the range $1, \dots, 10^8$ and $10^5, \dots, 10^8$.

As in Figure 24, we plot in Figure 31 the algebraic error $|||u_h - u_h^j|||$ and the corresponding contraction factor for the two test cases. We again observe a substantial reduction in the algebraic error at each step of Algorithm 8.1, and a fast convergence, confirming the efficiency of our algorithm.

As in Figure 25, Figure 32 illustrates the lower-bound error estimator $\eta^j$, the error $|||u_h - u_h^j|||$, and the upper-bound error estimator $\bar{\eta}^j$ for the two test cases. Finally, similarly, as in Figure 26, Figure 33 represents the lower effectivity index $I_{\text{eff}}$ and the upper effectivity index $\bar{I}_{\text{eff}}$ for both test cases. Similar overall behavior is observed.
<figure><img src="image_27.png" /><figcaption>Figure 26: [Test 2] Effectivity indices $I_{\text{eff}} = |||u_h - \hat{u}_h^j|||$ / $\eta^j$ (green solid line) and $I_{\text{eff}} = \eta^j / |||u_h - u_h^j|||$ (red dashed line).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 27: [Test 2] Left: elementwise 1st iteration flux errors $|||u_h - \hat{u}_h^1|||$$_K$. Center: elementwise 1st iteration algebraic error lower estimators $\alpha^1 |||\hat{u}_h^1 - \hat{u}_h^1|||$$_K$. Right: elementwise 1st iteration algebraic error upper estimators $|||\hat{u}_h^1 + \Pi_k^{\text{RTN}}(S \nabla \tilde{p}_h^2)|||$$_K$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 28: [Test 2] Left: elementwise final flux errors $|||u_h - \hat{u}_h^{14}|||$$_K$. Center: elementwise final algebraic error lower estimators $\alpha^{14} |||\hat{u}_h^{14} - \hat{u}_h^{14}|||$$_K$. Right: elementwise final algebraic error upper estimators $|||\hat{u}_h^{14} + \Pi_k^{\text{RTN}}(S \nabla \tilde{p}_h^{15})|||$$_K$.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 29: [Test 2] Left: elementwise 1st iteration potential errors $||p_h - p_h^1||_K$. Right: elementwise final potential errors $||p_h - p_h^{14}||_K$.</figcaption></figure>

<table><caption>Table 1: [Test 2] Number of iterations needed to reduce the initial algebraic error estimator $\eta^1$ by the factor $10^5$, depending on the diffusion contrast.</caption><tr><td>Diffusion contrast</td><td>10<sup>1</sup></td><td>10<sup>2</sup></td><td>10<sup>3</sup></td><td>10<sup>4</sup></td><td>10<sup>5</sup></td><td>10<sup>6</sup></td><td>10<sup>7</sup></td><td>10<sup>8</sup></td></tr><tr><td>Number of iterations</td><td>19</td><td>16</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td></tr></table>

## 9.4 Test 4: piecewise constant diffusion tensor not aligned with the coarse mesh and with irregular interfaces

In this test, we still build upon the Test 2 above, using this time the coarse and fine meshes shown in Figure 1 while varying the parameter $c(x, y)$. The variations in $c(x, y)$ are depicted in Figure 34 – we consider values of $c(x, y)$ equal to 1 and $10^7$. The objective is to study a situation where the interfaces of the coarse mesh cross the discontinuities of $c(x, y)$ and where the interfaces of the discontinuities are irregularly shaped.

As in Figure 24 and Figure 31, we plot in Figure 35 the algebraic error $|||u_h - u_h^j|||$ and the corresponding contraction factor. We again observe a substantial reduction in the algebraic error at each step of Algorithm 8.1, and a fast convergence, though a less “smooth”.

As in Figures 25 and 32, Figure 36 illustrates on the left the lower-bound error estimator $\eta^j$ and the error $|||u_h - u_h^j|||$ and on the right the lower effectivity index $I_{\text{eff}}$. Similar overall behavior is observed. For the upper bound error estimator $\eta^j$, we were not able to obtain pertinent results with our code due to difficulties in computing the gradient across irregularly shaped discontinuity interfaces.

## 10 Conclusions and perspectives

In this paper, we have introduced a nonoverlapping domain decomposition algorithm for saddle-point mixed finite element discretizations of the Darcy flow problem. We have also designed a posteriori error estimators that control the algebraic error from above and from below in a fully computable way (assuming round-off errors are neglected), and this for any iterative algebraic solver. The contraction of the domain decomposition algorithm and the reliability and efficiency of a the posteriori error estimates have been demonstrated on different 2D numerical examples.

This work can be extended in several directions. Firstly, a proof of contraction of the domain decomposition algorithm and of the efficiency of the a posteriori error estimates is to be given. Next, it would be interesting to extend the numerical validation to 3D cases in order to confirm the theoretical results beyond the 2D examples presented in this work. Finally, it would be
<figure><img src="image_31.png" /><figcaption>Figure 30: [Test 3] Variations of the coefficient $c(x, y)$ (in $\log_{10}$ scale) across the domain for the Skyscraper test case: Left: the discontinuities do not correspond to the coarse mesh. Right: the interfaces of the coarse mesh cross the discontinuities.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 31: [Test 3] Algebraic error $|||u_h - u_h^j|||$ (blue solid line, left scale). Contraction factor $|||u_h - u_h^{j+1}|||/|||u_h - u_h^j|||$ (red dashed line, right scale): Left: the discontinuities do not correspond to the coarse mesh. Right: the interfaces cross the discontinuities.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 32: [Test 3] Algebraic error $|||u_h - u_h^j|||$ (blue solid line), upper bound $\eta^j$ (red dotted line), and lower bound $\bar{\eta}^j$ (green dashed line): Left: the discontinuities do not correspond to the coarse mesh. Right: the interfaces cross the discontinuities.</figcaption></figure>

interesting to investigate the use of a step in the proposed algorithm as a preconditioner for Krylov subspace methods [30], which may accelerate their convergence.
<figure><img src="image_34.png" /><figcaption>Figure 33: [Test 3] Effectivity indices $I_{\text{eff}} = |||u_h - u_h^j|||/\eta^j$ (green solid line) and $I_{\text{eff}} = \eta^j/|||u_h - u_h^j|||$ (red dashed line): Left: the discontinuities do not correspond to the coarse mesh. Right: the interfaces cross the discontinuities.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 34: [Test 4] Variations of the coefficient $c(x, y)$ (in log$_{10}$ scale) across the domain.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 35: [Test 4] Algebraic error $|||u_h - u_h^j|||$ (blue solid line, left scale). Contraction factor $|||u_h - u_h^{j+1}|||/|||u_h - u_h^j|||$ (red dashed line, right scale).</figcaption></figure>

## A The lowest-order space $M_{h,0}$ on simplices

In this appendix, we discuss several points related to the elementwise postprocessing space $M_{h,k}$ used Section 6.2; we focus on simplicial meshes. We consider

$$M_{h,k} := \Pi_{K \in \mathcal{T}_h} M_{h,k}(K), \qquad (A.1)$$
<figure><img src="image_37.png" /><figcaption>Figure 36: [Test 4] Left : Algebraic error |||u_h - u_h^j||| (blue solid line) and lower bound η^j (green dashed line). Right: Effectivity index I_{eff} = |||u_h - u_h^j|||/η^j.</figcaption></figure>

where $\mathcal{M}_{h,k}(K)$ is a polynomial space on each mesh element $K \in \mathcal{T}_h$ such that a function $\tilde{p}_h \in \mathcal{M}_{h,k}(K)$ is uniquely prescribed by

$$ (\tilde{p}_h, q_h)_K = (p_h, q_h)_K \quad \forall q_h \in W_k(K), \qquad (A.2a) $$

$$ \langle \tilde{p}_h, \mu_h \rangle_F = \langle \lambda_h, \mu_h \rangle_F \quad \forall \mu_h \in \Psi_k(F) \text{ and } F \in \mathcal{F}_K \qquad (A.2b) $$

for some $(p_h, \lambda_h) \in W_k(K) \times \prod_{F \in \mathcal{F}_K} \Psi_k(F)$. Suitable spaces $\mathcal{M}_{h,k}(K)$ are developed in [3, 5, 38]; we detail some prominent possibilities in the lowest-order case $k = 0$ in Section A.2. Before, in Section A.1, we state an equivalence result for a general class elementwise postprocessings we are interested in.

## A.1 Equivalent postprocessing properties

Recall the S-weighted $L^2(K)$-orthogonal projection of (3.9). The following equivalence holds true for an arbitrary function $\tilde{p}_h$ in the Sobolev space $H^1(K)$ (unrelated to any polynomial space $\mathcal{M}_{h,k}(K)$).

**Lemma A.1.** (Equivalent postprocessing properties) Let a mesh element $K \in \mathcal{T}_h$ be fixed and let $(u_h, p_h, \lambda_h) \in \mathcal{RTN}_k(K) \times W_k(K) \times \prod_{F \in \mathcal{F}_K} \Psi_k(F)$ be such that

$$ \langle \lambda_h, v_h \cdot n_K \rangle_{\partial K} = -(S^{-1}u_h, v_h)_K + (p_h, \nabla \cdot v_h)_K \qquad (A.3) $$

for all $v_h \in \mathcal{RTN}_k(K)$, as in (3.13a) or, when $k = 0$, in (6.1). Then, a function $\tilde{p}_h \in H^1(K)$ satisfies

$$ (\tilde{p}_h, q_h)_K = (p_h, q_h)_K \quad \forall q_h \in W_k(K), \qquad (A.4a) $$

$$ \langle \tilde{p}_h, \mu_h \rangle_F = \langle \lambda_h, \mu_h \rangle_F \quad \forall \mu_h \in \Psi_k(F) \text{ and } F \in \mathcal{F}_K \qquad (A.4b) $$

if and only if

$$ (\tilde{p}_h, q_h)_K = (p_h, q_h)_K \quad \forall q_h \in W_k(K), \qquad (A.5a) $$

$$ -(S^{-1}\Pi_k^{\mathcal{RTN}}(S\nabla\tilde{p}_h), v_h)_K = (S^{-1}u_h, v_h)_K \quad \forall v_h \in \mathcal{RTN}_k(K), \qquad (A.5b) $$

or, equivalently for (A.5b),

$$ -\Pi_k^{\mathcal{RTN}}(S\nabla\tilde{p}_h) = u_h. \qquad (A.6) $$
Proof. We follow [3, 5, 9, 38] and show two implications, using the definition of the projector $\Pi_k^{RTN}$ from (3.9) and Green's theorem, leading to

$$ -(S^{-1}\Pi_k^{RTN}(S\nabla\tilde{p}_h), v_h)_K \stackrel{(3.9)}{=} -(\nabla\tilde{p}_h, v_h)_K = (\tilde{p}_h, \nabla\cdot v_h)_K - \langle\tilde{p}_h, v_h \cdot n_K\rangle_{\partial K} \quad (A.7) $$

for all $v_h \in \mathcal{RTN}_k(K)$. Recall also that $\nabla \cdot \mathcal{RTN}_k(K) = W_k(K)$ and $(\mathcal{RTN}_k(K) \cdot n_K)|_F = \Psi_k(F)$ for all $F \in \mathcal{F}_K$ and that (A.4a) and (A.5a) coincide.

1. $\Rightarrow$. (A.4) implies (A.5) Let $v_h \in \mathcal{RTN}_k(K)$ be arbitrary. We have

$$ \begin{align*} -(S^{-1}\Pi_k^{RTN}(S\nabla\tilde{p}_h), v_h)_K &\stackrel{(A.7)}{=} (\tilde{p}_h, \nabla \cdot v_h)_K - \langle\tilde{p}_h, v_h \cdot n_K\rangle_{\partial K} \\ &\stackrel{(A.4)}{=} (p_h, \nabla \cdot v_h)_K - \langle\lambda_h, v_h \cdot n_K\rangle_{\partial K} \\ &\stackrel{(A.3)}{=} (S^{-1}u_h, v_h)_K, \end{align*} $$

which is (A.5b).

2. $\Leftarrow$. (A.5) implies (A.4) Let $F \in \mathcal{F}_K$ and take $v_h \in \mathcal{RTN}_k(K)$ such that $(v_h \cdot n_K)|_{F'} = 0$ for all $F' \in \mathcal{F}_K, F \neq F'$. We have

$$ \begin{align*} \langle \tilde{p}_h, v_h \cdot n_K \rangle_F &\stackrel{(A.7)}{=} (S^{-1}\Pi_k^{RTN}(S\nabla\tilde{p}_h), v_h)_K + (\tilde{p}_h, \nabla \cdot v_h)_K \\ &\stackrel{(A.5b)}{=} -(S^{-1}u_h, v_h)_K + (\tilde{p}_h, \nabla \cdot v_h)_K \\ &\stackrel{(A.5a)}{=} -(S^{-1}u_h, v_h)_K + (p_h, \nabla \cdot v_h)_K \\ &\stackrel{(A.3)}{=} \langle \lambda_h, v_h \cdot n_K \rangle_F \end{align*} $$

which is (A.4b).

□

## A.2 Choices for the space $M_{h,0}$ in the lowest-order case

Let $K \in \mathcal{T}_h$ be a simplex. From [3, 5, 9, 38], the following are prominent possibilities for the polynomial space $M_{h,0}(K)$:

### A.2.1 Nonconforming $P_2$-bubbles

From [3, 5] it is known that one choice for the lowest-order space $M_{h,0}(K)$ is

$$ M_{h,0}(K) := P_1(K) \oplus B_2(K). $$

Here, recall, $P_1(K)$ is the space of polynomials of total degree less than or equal to 1 defined in the simplex $K$. For the basis, we choose the Crouzeix-Raviart functions $\psi_F \in P_1(K)$ [12] satisfying

$$ \frac{\langle \psi_F, 1 \rangle_F}{|F|} = 1 \quad \text{and} \quad \frac{\langle \psi_F, 1 \rangle_F'}{|F'|} = 0 \quad \text{for all } \mathcal{F}' \in \mathcal{F}_K, F \neq F'. \quad (A.8) $$

The space $B_2(K)$ is then the span of the following nonconforming $P_2$-bubble function (see [3, Sections 4.1 and 10])

$$ b_{2,K}(\boldsymbol{x}) = \begin{cases} 2 - 3(\psi_1^2(\boldsymbol{x}) + \psi_2^2(\boldsymbol{x}) + \psi_3^2(\boldsymbol{x})) & \text{if } d = 2, \\ 1 - 2(\psi_1^2(\boldsymbol{x}) + \psi_2^2(\boldsymbol{x}) + \psi_3^2(\boldsymbol{x}) + \psi_4^2(\boldsymbol{x})) & \text{if } d = 3. \end{cases} $$
Here $\psi_i$ are the hat functions, i.e., the affine functions that take value one at the vertex $i$ and zero at the other vertices of $K$. The functions $b_{2,K}$ vanish at the Gauss points of each $F \in \mathcal{F}_K$ and are of mean value 1 in $K$, or equivalently

$$ \frac{\langle b_{2,K}, 1 \rangle_F}{|F|} = 0 \quad \text{and} \quad \frac{\langle b_{2,K}, 1 \rangle_K}{|K|} = 1 \quad \text{for all } \mathcal{F} \in \mathcal{F}_K. \qquad (A.9) $$

From the above, we want to construct a function $\tilde{p}_h$ by using the definition (A.2) for $k=0$. This function is such that

$$ \tilde{p}_h = \sum_{F \in \mathcal{F}_K} \alpha_F \psi_F + \beta_K b_{2,K} \qquad (A.10) $$

for some constants $\alpha_F$ and $\beta_K$. To determine the constants $\alpha_F$ and $\beta_K$, we use (A.8) and (A.9) and obtain

$$ \langle \tilde{p}_h, 1 \rangle_F = \alpha_F \langle \psi_F, 1 \rangle_F + \beta_K \langle b_{2,K}, 1 \rangle_F = \alpha_F |F| \quad \forall F \in \mathcal{F}_K, \\ (\tilde{p}_h, 1)_K = \left( \sum_{F \in \mathcal{F}_K} \alpha_F \psi_F, 1 \right)_K + \beta_K |K|. $$

This together with (A.2) yields

$$ \alpha_F \stackrel{(A.2b)}{=} \frac{\langle \lambda_h, 1 \rangle_F}{|F|} \quad \text{and} \quad \beta_K \stackrel{(A.2a)}{=} \frac{(p_h, 1)_K}{|K|} - \frac{\left( \sum_{F \in \mathcal{F}_K} \alpha_F \psi_F, 1 \right)_K}{|K|}. $$

If (A.3) holds, we know from Lemma A.1 that the construction of $\tilde{p}_h$ is equivalent to finding $\tilde{p}_h \in \mathcal{M}_{h,0}(K)$ such that

$$ (\tilde{p}_h, 1)_K = (p_h, 1)_K \quad \text{and} \quad -\Pi_k^{\mathcal{RTN}}(S \nabla \tilde{p}_h) = u_h. $$

### A.2.2 Potential space of $S^{-1}\mathcal{RTN}_0(K)$

From [38], we know that we can define the space $\mathcal{M}_{h,0}(K)$ such that

$$ S \nabla \mathcal{M}_{h,0}(K) = \mathcal{RTN}_0(K). \qquad (A.12) $$

This is indeed possible from Assumption 2.1, item Ass1, and yields

$$ \mathcal{M}_{h,0}(K) \subseteq \mathcal{P}_2(K). $$

If (A.3) holds, we know from Lemma A.1 that the construction of $\tilde{p}_h \in \mathcal{M}_{h,0}(K)$ using definition (A.2) is equivalent to finding $\tilde{p}_h \in \mathcal{M}_{h,0}(K)$ such that

$$ (\tilde{p}_h, 1)_K = (p_h, 1)_K \quad \text{and} \quad -S \nabla \tilde{p}_h = u_h, \qquad (A.13) $$

where we have crucially used (A.12) which gives $-\Pi_k^{\mathcal{RTN}}S \nabla \tilde{p}_h = -S \nabla \tilde{p}_h$.

In detail (for $d=2$), let $u_h = (\alpha_1 + \beta x, \alpha_2 + \beta y)$. In order to satisfy (A.13), the function $\tilde{p}_h = p_1 x^2 + p_2 y^2 + p_3 xy + p_4 x + p_5 y + p_6$ is such that

$$ p_1 = -\frac{\beta S_{22}}{2|S|}, \quad p_2 = -\frac{\beta S_{11}}{2|S|}, \quad p_3 = \frac{\beta S_{12}}{|S|}, \quad p_4 = \frac{\alpha_2 S_{12} - \alpha_1 S_{22}}{|S|}, \quad p_5 = \frac{\alpha_1 S_{12} - \alpha_2 S_{11}}{|S|}, $$

and $p_6$ is given so that $(\tilde{p}_h, 1)_K = (p_h, 1)_K$ is satisfied. Here $S_{ij}$ and $|S|$ denote the elements and the determinant of $S$, respectively. Similar algebraic manipulations give a unique function $\tilde{p}_h$ satisfying (A.13) also in the case $d=3$.
### A.2.3 $P_3$-bubbles

For $d=2$, another possible choice for the space $M_{h,0}(K)$ is, following [3, 5],

$$ M_{h,0}(K) := P_1(K) \oplus B_3(K). \qquad (A.14) $$

The space $P_1(K)$ is as defined in Section A.2.1 and $B_3(K)$ is the span of the following $P_3$-bubble function

$$ b_{3,K}(\boldsymbol{x}) = \psi_1(\boldsymbol{x})\psi_2(\boldsymbol{x})\psi_3(\boldsymbol{x}). $$

Note that the bubble function $b_{3,K}$ vanishes on each edge of the triangle $K$. We can thus use the same arguments as in Section A.2.1 to construct an appropriate $\tilde{p}_h \in M_{h,0}(K)$.

### A.2.4 Remarks

We finish by two remarks.

**Remark A.2 (Uniqueness).** As it follows from the above, the choice of the space $M_{h,0}(K)$ is not unique. However, for each fixed space $M_{h,0}(K)$, the construction of $\tilde{p}_h$ satisfying (A.2) is unique.

**Remark A.3 (Approaches of Sections A.2.1 and A.2.2).** In the case of $d=2$, $S = \alpha I$ with $\alpha \in \mathbb{R}$, and for equilateral triangles, the spaces $M_{h,0}(K)$ discussed in Sections A.2.1 and A.2.2 turn out to be equal, see [9, Lemma 4.2]. In general, though, approaches of Section A.2.1 and Section A.2.2 generate different spaces $M_{h,0}(K)$ and thus different functions $\tilde{p}_h$.

## References

[1] ALI HASSAN, S., JAPHET, C., KERN, M., AND VOHRALÍK, M. A posteriori stopping criteria for optimized Schwarz domain decomposition algorithms in mixed formulations. *Comput. Methods Appl. Math.* 18, 3 (2018), 495–519.

[2] ANCIAUX-SEDRAKIAN, A., GRIGORI, L., JORTI, Z., PAPEŽ, J., AND YOUSEF, S. Adaptive solution of linear systems of equations based on a posteriori error estimators. *Numer. Algorithms* 84, 1 (2020), 331–364.

[3] ARBOGAST, T., AND CHEN, Z. On the implementation of mixed methods as nonconforming methods for second-order elliptic problems. *Math. Comp.* 64, 211 (1995), 943–972.

[4] ARBOGAST, T., COWSAR, L. C., WHEELER, M. F., AND YOTOV, I. Mixed finite element methods on nonmatching multiblock grids. *SIAM J. Numer. Anal.* 37, 4 (2000), 1295–1315.

[5] ARNOLD, D. N., AND BREZZI, F. Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates. *RAIRO Modél. Math. Anal. Numér.* 19, 1 (1985), 7–32.

[6] BENZI, M., GOLUB, G. H., AND LIESEN, J. Numerical solution of saddle point problems. *Acta Numer.* 14 (2005), 1–137.

[7] BOFFI, D., BREZZI, F., AND FORTIN, M. *Mixed finite element methods and applications*, vol. 44 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013.

[8] BOON, W. M., GLÄSER, D., HELMIG, R., AND YOTOV, I. Flux-mortar mixed finite element methods with multipoint flux approximation. *Comput. Methods Appl. Mech. Engrg.* 405 (2023), Paper No. 115870, 28.
[9] CHEN, Z. Analysis of mixed methods using conforming and nonconforming finite element methods. *RAIRO Modél. Math. Anal. Numér.* **27**, 1 (1993), 9–34.

[10] CIARLET, JR., P., JAMELOT, E., AND KPADONOU, F. D. Domain decomposition methods for the diffusion equation with low-regularity solution. *Comput. Math. Appl.* **74**, 10 (2017), 2369–2384.

[11] COWSAR, L. C., MANDEL, J., AND WHEELER, M. F. Balancing domain decomposition for mixed finite elements. *Math. Comp.* **64**, 211 (1995), 989–1015.

[12] CROUZEIX, M., AND RAVIART, P.-A. Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I. *Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge* **7**, no. , no. R-3 (1973), 33–75.

[13] DOBREV, V., KOLEV, T., LEE, C. S., TOMOV, V., AND VASSILEVSKI, P. S. Algebraic hybridization and static condensation with application to scalable $H(\text{div})$ preconditioning. *SIAM J. Sci. Comput.* **41**, 3 (2019), B425–B447.

[14] ERN, A., STEPHANSEN, A. F., AND ZUNINO, P. A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity. *IMA J. Numer. Anal.* **29**, 2 (2009), 235–256.

[15] ERN, A., AND VOHRALÍK, M. Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs. *SIAM J. Sci. Comput.* **35**, 4 (2013), A1761–A1791.

[16] ERN, A., AND VOHRALÍK, M. Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations. *SIAM J. Numer. Anal.* **53**, 2 (2015), 1058–1081.

[17] EWING, R. E., AND WANG, J. Analysis of the Schwarz algorithm for mixed finite elements methods. *RAIRO Modél. Math. Anal. Numér.* **26**, 6 (1992), 739–756.

[18] GANIS, B., AND YOTOV, I. Implementation of a mortar mixed finite element method using a multiscale flux basis. *Comput. Methods Appl. Mech. Engrg.* **198**, 49-52 (2009), 3989–3998.

[19] GLOWINSKI, R., AND WHEELER, M. F. Domain decomposition and mixed finite element methods for elliptic problems. In *First International Symposium on Domain Decomposition Methods for Partial Differential Equations (Paris, 1987)*. SIAM, Philadelphia, 1988, pp. 144–172.

[20] HABERL, A., PRAETORIUS, D., SCHIMANKO, S., AND VOHRALÍK, M. Convergence and quasi-optimal cost of adaptive algorithms for nonlinear operators including iterative linearization and algebraic solver. *Numer. Math.* **147**, 3 (2021), 679–725.

[21] HECHT, F. New development in FreeFem++. *J. Numer. Math.* **20**, 3-4 (2012), 251–265.

[22] MANDEL, J. Balancing domain decomposition. *Comm. Numer. Methods Engrg.* **9**, 3 (1993), 233–241.

[23] MATHEW, T. P. Schwarz alternating and iterative refinement methods for mixed formula-tions of elliptic problems. I. Algorithms and numerical results. *Numer. Math.* **65**, 4 (1993), 445–468.

[24] MIRAÇI, A., PAPEŽ, J., AND VOHRALÍK, M. A multilevel algebraic error estimator and the corresponding iterative solver with p-robust behavior. *SIAM J. Numer. Anal.* **58**, 5 (2020), 2856–2884.
[25] MIRAÇI, A., PAPEŽ, J., AND VOHRALÍK, M. A-posteriori-steered *p*-robust multigrid with optimal step-sizes and adaptive number of smoothing steps. *SIAM J. Sci. Comput.* **43**, 5 (2021), S117-S145.

[26] NÉDÉLEC, J.-C. Mixed finite elements in $\mathbb{R}^3$. *Numer. Math.* **35**, 3 (1980), 315-341.

[27] OH, D.-S., WIDLUND, O. B., ZAMPINI, S., AND DOHRMANN, C. R. BDDC algorithms with deluxe scaling and adaptive selection of primal constraints for Raviart-Thomas vector fields. *Math. Comp.* **87**, 310 (2018), 659-692.

[28] PAPEŽ, J., RÜDE, U., VOHRALÍK, M., AND WOHLMUTH, B. Sharp algebraic and total a posteriori error bounds for *h* and *p* finite elements via a multilevel approach. Recovering mass balance in any situation. *Comput. Methods Appl. Mech. Engrg.* **371** (2020), 113243.

[29] PAPEŽ, J., STRAKOŠ, Z., AND VOHRALÍK, M. Estimating and localizing the algebraic and total numerical errors using flux reconstructions. *Numer. Math.* **138**, 3 (2018), 681-721.

[30] PEARSON, J. W., AND PESTANA, J. Preconditioners for Krylov subspace methods: an overview. *GAMM-Mitt.* **43**, 4 (2020), e202000015, 35.

[31] RAVIART, P.-A., AND THOMAS, J.-M. A mixed finite element method for 2nd order elliptic problems. In *Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975)*. Springer, Berlin, 1977, pp. 292-315. Lecture Notes in Math., Vol. 606.

[32] REY, V., GOSSELET, P., AND REY, C. Strict lower bounds with separation of sources of error in non-overlapping domain decomposition methods. *Internat. J. Numer. Methods Engrg.* **108**, 9 (2016), 1007-1029.

[33] ŠÍSTEK, J., BŘEZINA, J., AND SOUSEDÍK, B. BDDC for mixed-hybrid formulation of flow in porous media with combined mesh dimensions. *Numer. Linear Algebra Appl.* **22**, 6 (2015), 903-929.

[34] SOLOVSKÝ, J., FUČÍK, R., AND ŠÍSTEK, J. BDDC for MHFEM discretization of unsteady two-phase flow in porous media. *Comput. Phys. Commun.* **271** (2022), Paper No. 108199, 14.

[35] SOUSEDÍK, B. Nested BDDC for a saddle-point problem. *Numer. Math.* **125**, 4 (2013), 761-783.

[36] TOSELLI, A. Neumann-Neumann methods for vector field problems. *Electron. Trans. Numer. Anal.* **11** (2000), 1-24.

[37] TU, X. A BDDC algorithm for a mixed formulation of flow in porous media. *Electron. Trans. Numer. Anal.* **20** (2005), 164-179.

[38] VOHRALÍK, M. A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations. *SIAM J. Numer. Anal.* **45**, 4 (2007), 1570-1599.