A posteriori error estimates and adaptivity for locally conservative methods. Inexpensive implementation and evaluation, polytopal meshes, iterative linearization and algebraic solvers, and applications to complex porous media flows

TL;DR

The paper develops a unified, guaranteed, and computable framework for a posteriori error estimation and adaptivity for locally conservative finite-volume-type methods on polytopal meshes. It leverages $H^1$-conforming potential reconstructions and $\bm{H}(\mathrm{div})$-conforming equilibrated flux reconstructions to bound the total error while identifying spatial, temporal, linearization, and algebraic components. The estimators are designed to be inexpensive to evaluate, enabling adaptive strategies that balance error components, enforce mass balance at every solver stage, and drive adaptive mesh/time refinement and solver stopping criteria. The approach applies across Poisson, steady linear and nonlinear diffusion (including Darcy flow), and unsteady multiphase compositional Darcy flows, with extensive numerical experiments in 2D and 3D demonstrating robustness, asymptotic exactness in favorable cases, and practical speedups for complex porous-media simulations.

Abstract

A posteriori estimates give bounds on the error between the unknown solution of a partial differential equation and its numerical approximation. We present here the methodology based on H1-conforming potential and H(div)-conforming equilibrated flux reconstructions, where the error bounds are guaranteed and fully computable. We consider any lowest-order locally conservative method of the finite volume type and treat general polytopal meshes. We start by a pure diffusion problem and first address the discretization error. We then progressively pass to more complicated model problems, up to complex multiphase multicomponent flow in porous media, and also take into account the errors arising in iterative linearization of nonlinear problems and in algebraic resolution of systems of linear algebraic equations. We focus on the ease of implementation and evaluation of the estimates. In particular, the evaluation of our estimates is explicit and inexpensive, since it merely consists in some local matrix-vector multiplications. Here, on each mesh element, the matrices are either directly inherited from the given numerical method, or easily constructed from the element geometry, while the vectors are the algebraic unknowns of the flux and potential approximations on the given element. Our mtehodology leads to an easy-to-implement and fast-to-run adaptive algorithm with guaranteed overall precision, adaptive stopping criteria for nonlinear and linear solvers, and adaptive space and time mesh refinements and derefinements. Progressively along the theoretical exposition, numerical experiments on academic benchmarks as well as on real-life problems in two and three space dimensions illustrate the performance of the derived methodology. The presentation is largely self-standing, developing all the details and recalling all necessary basic notions.

Figures (38)

• Figure 1: Example of a function belonging to the Sobolev space $H^1_0(\Omega)$ (left; the function actually lies in $C^0(\overline\Omega)$, but not in $C^1(\overline\Omega)$). Example of a function not belonging to the Sobolev space $H^1_0(\Omega)$ but belonging to the broken Sobolev space $H^1(\mathcal{T}_h)$ for a mesh $\mathcal{T}_h$ composed of two triangles (right; note that there is no trace-continuity on the common face of the two mesh elements)
• Figure 2: Example of a function belonging to the Sobolev space $\bm{H}(\mathrm{div},\Omega)$ (left; on the interface $x=0$ between the two subdomains $\Omega_1 = (-1,0)\times(0,1)$ and $\Omega_2 = (0,1)\times(0,1)$, $\bm{v}|_{\Omega_1} {\cdot} \bm{n} = \bm{v}|_{\Omega_2} {\cdot} \bm{n}$, with $\bm{n}=(1,0)^\mathrm{t}$; $\bm{v}$ is not continuous for each component (the $y$ component of $\bm{v}$ is discontinuous, as it passes from value $0$ in $\Omega_1$ to a nonzero value in $\Omega_2$) but $\bm{v}$ is normal-trace continuous). Example of a function not belonging to the Sobolev space $\bm{H}(\mathrm{div},\Omega)$ but belonging to the broken Sobolev space $\bm{H}(\mathrm{div},\mathcal{T}_h)$ (right, there is no normal trace-continuity across the mesh face at $x=0$).
• Figure 3: Simplicial mesh $\mathcal{T}_h$ for $d=2$ (left) and $d=3$ (right)
• Figure 4: Spaces ${\mathcal{P}}_k({K})$ and point values uniquely fixing a function $v_h \in {\mathcal{P}}_k({K})$ (the so-called Lagrange nodes) for a mesh element ${K} \in \mathcal{T}_h$. Polynomial degrees $k=1$ and $k=2$, space dimensions $d=2$ (left) and $d=3$ (right)
• Figure 5: Spaces $\bm{\mathcal{R\space T}}\space_0({K})$ and normal fluxes uniquely fixing a function $v_h \in \bm{\mathcal{R\space T}}\space_0({K})$ for a simplex ${K} \in \mathcal{T}_h$. Space dimensions $d=2$ (left) and $d=3$ (right)

Theorems & Definitions (71)

• Definition 3.1: Weak partial derivative
• Definition 3.2: Weak gradient
• Definition 3.3: The space $H^1(\Omega)$
• Definition 3.4: The space $H^1_0(\Omega)$
• Remark 3.5: Trace and trace continuity
• Definition 3.6: Weak divergence
• Definition 3.7: The space $\bm{H}(\mathrm{div},\Omega)$
• Remark 3.8: Normal trace and normal-trace continuity
• Theorem 3.9: Green theorem
• Lemma 3.10: Trace continuity over mesh faces