# Generalizing the Multiscale Hybrid-Mixed Method for Reactive-Advective-Diffusive Equations

Rodolfo Araya, Fabrice Jaillet, Diego Paredes, Frédéric Valentin

► To cite this version:

Rodolfo Araya, Fabrice Jaillet, Diego Paredes, Frédéric Valentin. Generalizing the Multiscale Hybrid-Mixed Method for Reactive-Advective-Diffusive Equations. Computer Methods in Applied Mechanics and Engineering, 2024, 428, 10.1016/j.cma.2024.117089. hal-04321639v2

HAL Id: hal-04321639

https://hal.science/hal-04321639v2

Submitted on 26 Jun 2024

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>
# Generalizing the Multiscale Hybrid-Mixed Method for Reactive-Advective-Diffusive Equations

Rodolfo Arayaª, Fabrice Jailletᵇ, Diego Paredesª,*, Frédéric Valentinᶜ

## Abstract

We propose a new family of multiscale hybrid mixed methods (MHM) for the reactive-advection-diffusive (RAD) equation in complex domains. It generalizes the MHM methods originally proposed in Harder, Paredes and Valentin (2013 and 2015) to polytopal meshes and covers all asymptotic regimes of the model within a single mathematical framework. As a result, the skeletal MHM method changes its structure automatically, from primal to mixed forms, depending on the asymptotic of local RAD solutions, which respond to multiscale basis functions at the element level. We establish the existence, uniqueness, and optimality of the MHM solution with respect to two-scale mesh parameters, relating it to the solution of a discrete primal hybrid version of the RAD model. Furthermore, we estimate the condition number of the matrices associated with the local problems responsible for upscaling, from which we establish upper limits for the condition number of the algebraic system associated with the MHM method. Numerical experiments validate theoretical results.

**Keywords:** Finite element method, multiscale problems, reaction-advection-diffusion model, polytopal meshes

## 1. Introduction

Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, be an open, bounded domain with polytopal Lipchitz boundary $\partial\Omega$. We consider the reaction-advection-diffusion problem of finding $u: \Omega \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \nabla \cdot (-\varepsilon \nabla u + \alpha u) + \sigma u = f & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega, \end{array} \right. \qquad (1) $$

where, $\varepsilon \in L^\infty(\Omega)^{d \times d}$, $\alpha \in W^{1,\infty}(\Omega)^d$ and $\sigma \in L^\infty(\Omega)$ are the diffusive, advective, and reactive coefficients, respectively, and $f \in L^2(\Omega)$ is the source term. Hereafter, we use the standard notation for Sobolev spaces. We assume that the diffusive tensor $\varepsilon := \{\varepsilon_{ij}\}$ is symmetric and uniformly elliptic, i.e., there exist positive constants $\varepsilon_{\min}$ and $\varepsilon_{\max}$, such that

$$ \varepsilon_{\min} |z|^2 \le \varepsilon_{ij}(x) z_i z_j \le \varepsilon_{\max} |z|^2 \quad \text{for all } z = \{z_i\} \in \mathbb{R}^d, \text{ a.e. } x \in \bar{\Omega}, \qquad (2) $$

where $| \cdot |$ is the Euclidean norm, and we use the Einstein summation convention, i.e., repeated indices indicate summation. The advective and reactive coefficients satisfy

$$ \frac{1}{2} \nabla \cdot \alpha(x) + \sigma(x) \ge 0, \quad \text{a.e. } x \text{ in } \bar{\Omega}. \qquad (3) $$

<footnote>ªDepartamento de Ingeniería Matemática and CI²MA, Universidad de Concepción, Concepción, Chile</footnote>

<footnote>ᵇUCBL, CNRS, INSA Lyon, LIRIS, UMR 5205, F-01000, Bourg-en-Bresse, France</footnote>

<footnote>ᶜDepartment of Mathematical and Computational Methods, National Laboratory for Scientific Computing, Petrópolis, RJ, Brazil</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: rodolfo.araya@udec.cl (Rodolfo Araya), fabrice.jaillet@liris.cnrs.fr (Fabrice Jaillet), dparedes@udec.cl (Diego Paredes), valentin@lncc.br (Frédéric Valentin)</footnote>
The exact solution of (1) presents boundary layers when the reaction or advection coefficients dominate the diffusion ones. Fast variations also characterize the solution of (1) in dominant diffusion regimes when the diffusion coefficient is multiscale. A typical example is the Poisson problem with a highly oscillatory coefficient. In those scenarios, the standard Galerkin method with piecewise polynomial interpolations needs fine meshes, which is costly, to approximate the exact solution accurately. Thus, alternative methods have been developed in the last decades to improve the accuracy of the finite element method on coarse meshes. Among them are the stabilized methods, which are based on adding local variational forms to the Galerkin method related to the residual of the Lagrange equation. They have been applied to dominant advection problems since the seminal works [12] and [20, 21], for example. Other possibilities are enriched methods, like Residual-Free Bubbles (RFB) and the (Petrov-)Galerkin enriched method ((P)GEM) [1, 19]. In those cases, the standard approximation spaces are enhanced with non-polynomial functions driven by the original operator at the element level (see [11, 10] for interconnection between stabilized and enriched methods).

On the other hand, there has been a vast literature on multiscale finite element methods since the seminal work [5]. Multiscale methods share the commonality of having global problems built on the solution of local problems that scale up the structures under the mesh. An attractive feature is that independent problems can locally compute the multiscale basis functions. Examples are the variational multiscale method (VMS) [27], the multiscale finite element method (MsFEM) and its generalization (GMsFEM) [17], the multiscale heterogeneous method (HMM) [16], multiscale mortar method [3], the local orthogonal decomposition (LOD) method [30], the hybrid localized spectral decomposition (LSD) method [29] and the higher order hybrid multiscale method (MsHHO) [15], to name a few. Enriched methods are also closely related to multiscale methods, as established, for instance, in [33].

The MHM method is a member of the multiscale finite element family. It was initially proposed in [24] for the Poisson equation, and a priori and a posteriori error estimates for the MHM method were established in [2]. It was extended to polygonal meshes in [7] and, recently, to general polytopal meshes in the context of the elasticity problem in [22]. The MHM method is a by-product of a hybrid formulation that starts at the continuous level defined on a coarse partition, which characterizes the exact solution in local and global contributions. When discretized, it dissociates local from global problems. The global formulation is responsible for the degrees of freedom on the mesh skeleton, and the local problems provide the multiscale basis functions. In this context, the MHM method was extended to deal with the reactive-adveptive-diffusive (RAD) model in [25] and numerically validated in two-dimensional problems. The multiscale basis functions satisfy local RAD problems with prescribed Robin boundary conditions, which were assumed to be well-posed. Indeed, the MHM method for the RAD model was formally built and validated numerically in singularly perturbed problems (e.g., advective and reactive dominate regimes), yet the numerical analysis of the method was absent in [25].

So, the first objective is to fill this theoretical gap and present a numerical analysis of the two-level method introduced in [25]. The two-level nomenclature means that the multiscale basis are locally approximated by a second level of discretization. For this, we establish the well-posedness and optimal convergence of a new family of two-level MHM methods for the RAD model in terms of the two-scale mesh parameters that include the method in [25] as a particular case. In addition to the numerical analysis of the method in [25], in this work we

* propose a novel MHM method (see (42)) that automatically changes its structure from a primal formulation to a mixed one so that the local RAD model responsible for the multiscale basis is well-posed in all asymptotic regimes. In the process, we establish sufficient conditions for the RAD model with Robin boundary condition to be well-posed, which is surprisingly absent in the literature to our knowledge. Thus, the new method generalizes the method in [25] and recovers the original method proposed in [24] in cases of null reaction and advection coefficients. Such a construction depends on a new space decomposition that differs from the one proposed in [25];

* establish a numerical analysis of the new method on general polytopal grids. Notably, we prove the existence and uniqueness under compatible conditions between the interpolation space used locally to approximate the solution of local problems and the interpolation space in the partition skeleton.
Furthermore, we establish optimal convergence results based on the mesh parameters (see Theorem 3), highlighting the interaction between the diameter of the mesh skeleton and the diameter of the sub-meshes used to approximate local problems. Particularly, we show that fixing the coarse mesh and refining the mesh skeleton convergence is achieved. Such analysis is the first for the MHM method applied to RAD equations, even for simplicial meshes. The proof is based on a pre-established equivalence between the two-level MHM method and the discretized primal hybrid formulation of the RAD equation. Such a technique differs from that used in [24] (see [6, 26] for similar approaches);

* calculate upper bounds for the condition numbers of the matrices underlying the local and global problems in the MHM algorithm regarding to the physical coefficients and mesh parameters. We show that, due to the spatial decomposition inherent to the proposed MHM method, we maintain a favorable matrix condition number when compared to that obtained by the standard Galerkin method in piecewise continuous polynomial spaces. This study is the first for MHM methods that prove to be useful for solving global problems via iterative solvers, aiming to reduce the computational cost of the MHM algorithm stage, which is not naturally parallelizable;

* validate the MHM method on polytope meshes through two and three-dimensional problems. The numerical evaluation extends the results presented in [25] and verifies convergence estimates and theoretical results related to matrix condition number.

The outline of this work is as follows: In Section 2, we present the continuous and discrete hybrid formulation defined on polytopal meshes and some preliminary results. Section 3 is dedicated to the introduction of the MHM method based on a new decomposition of spaces and its numerical analysis. We establish estimates for the condition numbers of matrices associated with local and global systems in Section 4, and present practical details of the MHM algorithm in Section 5. Also, Section 5 presents numerical validations of the method, and conclusions are addressed in Section 6.

## 2. Primal hybridization and preliminary results

### 2.1. Partitions

Let {$\mathcal{P}_H$}$_{H>0}$ be a family of conforming partitions of $\bar{\Omega}$, composed by closed, bounded, disjoint polytopes $K$ such that $\bar{\Omega} = \cup_{K \in \mathcal{P}_H} K$. The diameter of $K \in \mathcal{P}_H$ is $\mathcal{H}_K$, and without loss of generality, we shall use hereafter the terminology employed for three-dimensional domains. Each polytope $K$ has a boundary $\partial K$ consisting of facets $E$. We collect the boundaries associated with $\mathcal{P}_H$ in $\partial \mathcal{P}_H$ and its facets in $\mathcal{E}^H$, that is $\partial \mathcal{P}_H = \{\partial K : K \in \mathcal{P}_H\}$, and

$$ \mathcal{E}^H = \left\{ \begin{array}{l} E = \partial K \cap \partial K', \text{ or} \\ E = \partial K \cap \partial \Omega \end{array} \middle| \begin{array}{l} K, K' \in \mathcal{P}_H, K \neq K', \text{ and} \\ E \text{ is not an edge nor a point} \end{array} \right\}. $$

For each $E \in \mathcal{E}^H$ we associate a normal vector $n^E$ coinciding with $n$ if $E \subset \partial \Omega$, and we further denote by $n^K$ the outward normal vector on $\partial K$ for each $K \in \mathcal{P}_H$. We introduce {$\mathcal{E}_H^K$}$_{H>0}$ a shape regular family of meshes for $\mathcal{E}^H$, such as each $E \in \mathcal{E}^H$ is split into simplicial facets $F$ of diameter $H_F \le H := \max_{F \in \mathcal{E}_H} H_F$. By $\mathcal{E}^H(K)$ and $\mathcal{E}_H^K(K)$ we understand the restriction of the sets $\mathcal{E}^H$ and $\mathcal{E}_H^K$ to $K \in \mathcal{P}_H$. The following technical assumption on $\mathcal{E}_H^K(K)$ will require that neighbouring facets are not too dissimilar:

(A1): The mesh $\mathcal{E}_H^K(K)$ is such that for all $K \in \mathcal{P}_H$, a conforming and shape regular simplicial triangulation $\Xi_H(K)$ of $K$ can be constructed such that its restriction on $\partial K$ coincides with $\mathcal{E}_H^K(K)$.

The conforming triangulation $\Xi_H := \cup_{K \in \mathcal{P}_H} \Xi_H(K)$ will be useful in the proofs below but not used explicitly in the implementation of the method. We will also need the following definitions:

- for each $F \in \mathcal{E}_H^K$, we denote $\kappa_F^K$ the only element in $\Xi_H(K)$ such that $F = \partial \kappa_F^K \cap \partial K$ where $\partial \kappa_F^K$ is the boundary of $\kappa_F^K$. For simplicity, let us assume that, for two $F_1, F_2 \in \mathcal{E}_H$, $\kappa_{F_1}^K \neq \kappa_{F_2}^K$;
- for each $K \in \mathcal{P}_H$, we introduce a shape regular family of simplicial triangulations $\{\mathcal{T}_h^K\}_{h>0}$ built in the following way: For each $K$

(i) the triangulation $\Xi_H(K)$ undergoes a red refinement, i.e., a tetrahedron $K$ is refined into 8 tetrahedra by connecting the midpoints of the original edges (for details, refer to [28]). The resulting triangulation is called minimal triangulation;

(ii) the family $\{\mathcal{T}_h^K\}_{h>0}$ is formed by red refinements of minimal triangulation.

The diameter of $\kappa \in \mathcal{T}_h^K$ is denoted by $h_\kappa$, and $h := \max_{K \in \mathcal{P}_H} \max_{\kappa \in \mathcal{T}_h^K} h_\kappa$, and $\mathcal{T}_h := \cup_{K \in \mathcal{P}_H} \mathcal{T}_h^K$. It is important to highlight that, if $E = K_1 \cap K_2 \in \mathcal{E}^H$, where $K_1, K_2 \in \mathcal{P}_H$, then the traces of the two neighboring triangulations $\mathcal{T}_h^{K_1}$ and $\mathcal{T}_h^{K_2}$ do not necessarily coincide (see Figure 1 for a two-dimensional sketch).

<figure><img src="image_2.png" /><figcaption>Figure 1: The domain $\Omega$ is partitioned by a conforming polygonal mesh $\mathcal{P}_H$. Two sub-meshes (green and orange) discretize two different elements $K$ of $\mathcal{P}_H$ with different granularity, where $\kappa \in \mathcal{T}_h^K$. The red line represents an edge $E \in \mathcal{E}^H$, and the blue line an element $F \in \mathcal{E}_H^H$ of the skeleton mesh.</figcaption></figure>

## 2.2. Broken spaces and norms

Let $K \in \mathcal{P}_H$ and $m \ge 1$, we consider the local space $H^m(K)$ equipped with the semi-norm $\|\cdot\|_{m,K}$ and the norm $\|\cdot\|_{m,K}$, with their usual definitions. Thus, we define the broken space

$$ H^m(\mathcal{P}_H) := \{v \in L^2(\Omega) : v|_K \in H^m(K), \text{ for all } K \in \mathcal{P}_H\}, \quad (4) $$

equipped with the broken (semi-)norm,

$$ |v|_{m, \mathcal{P}_H} := \left( \sum_{K \in \mathcal{P}_H} |v|_{m,K}^2 \right)^{\frac{1}{2}} \quad \text{and} \quad \|v\|_{m, \mathcal{P}_H} := \left( \sum_{K \in \mathcal{P}_H} \|v\|_{m,K}^2 \right)^{\frac{1}{2}}. $$

Also, we set $V := H^1(\mathcal{P}_H)$ and equip it with the norm $\|v\|_V := \left(\sum_{K \in \mathcal{P}_H} \|v\|_{V(K)}^2\right)^{\frac{1}{2}}$, with $\|v\|_{V(K)} := (|v|_{1,K}^2 + \frac{1}{d_\Omega^2} \|v\|_{0,K}^2)^{\frac{1}{2}}$ for all $v \in H^1(K)$, where $d_\Omega := \text{diam}\,\Omega$. As usual, $(\cdot, \cdot)_D$ stands for the $L^2(D)$ inner product for a measurable set $D \subset \mathbb{R}^d$, and for all $u, v \in V$

$$ (u, v)_{\mathcal{P}_H} := \sum_{K \in \mathcal{P}_H} (u, v)_K \quad \text{and} \quad (u, v)_V := \frac{1}{d_\Omega^2} (u, v)_{\mathcal{P}_H} + (\nabla u, \nabla v)_{\mathcal{P}_H}. $$
We define the space of Lagrange multipliers as follows

$$ \Lambda := \{ \boldsymbol{q} \cdot \boldsymbol{n}^K |_{\partial K} , \text{ for all } K \in \mathcal{P}_H : \boldsymbol{q} \in H(\mathrm{div}; \Omega) \}, $$

and equip it with the norm

$$ \|\mu\|_{\Lambda} := \inf \{ \|\boldsymbol{q}\|_{\mathrm{div}} : \boldsymbol{q} \in H(\mathrm{div}; \Omega) \text{ and } \boldsymbol{q} \cdot \boldsymbol{n}^K|_{\partial K} = \mu|_{\partial K}, \text{ for all } K \in \mathcal{P}_H \}, $$

where

$$ \|\boldsymbol{q}\|_{\mathrm{div}} := ( \|\boldsymbol{q}\|_{0,\Omega}^2 + d_{\Omega}^2 \|\nabla \cdot \boldsymbol{q}\|_{0,\Omega}^2 )^{\frac{1}{2}} \quad \text{for all } \boldsymbol{q} \in H(\mathrm{div}; \Omega). $$

We adopt the following notation

$$ \langle \mu, v \rangle_{\partial \mathcal{P}_H} := \sum_{K \in \mathcal{P}_H} \langle \mu, v \rangle_{\partial K} \quad \text{for all } (\mu, v) \in \Lambda \times V, \qquad (5) $$

where $\langle \cdot, \cdot \rangle_{\partial K}$ stands for the duality product between $H^{-\frac{1}{2}}(\partial K)$ and $H^{\frac{1}{2}}(\partial K)$. Notice that, if $v \in H_0^1(\Omega)$ then

$$ \langle \mu, v \rangle_{\partial \mathcal{P}_H} = 0 \quad \text{for all } \mu \in \Lambda. \qquad (6) $$

Also, from [22], the following identity holds

$$ \|\mu\|_{\Lambda} = \sup_{v \in V} \frac{\langle \mu, v \rangle_{\partial \mathcal{P}_H}}{\|v\|_V} \quad \text{for all } \mu \in \Lambda, \qquad (7) $$

and then

$$ \langle \mu, v \rangle_{\partial \mathcal{P}_H} = \frac{\langle \mu, v \rangle_{\partial \mathcal{P}_H}}{\|v\|_V} \|v\|_V \le \sup_{w \in V} \frac{\langle \mu, w \rangle_{\partial \mathcal{P}_H}}{\|w\|_V} \|v\|_V \le \|\mu\|_{\Lambda} \|v\|_V. \qquad (8) $$

Above and hereafter, we lighten the notation and understand the supremum to be taken over sets excluding the zero function, even though this is not specifically indicated.

Let $k \ge 1$, and $K \in \mathcal{P}_H$. We introduce the following local finite element space

$$ V_h(K) := \{v_h \in C^0(K) : v_h|_{\kappa} \in \mathbb{P}_k(\kappa), \text{ for all } \kappa \in \mathcal{T}_h^K \}. \qquad (9) $$

The global counterpart of $V_h(K)$ on the partition $\mathcal{P}_H$ is denoted by $V_h := \prod_{K \in \mathcal{P}_H} V_h(K)$. Also, let $\ell \ge 0$ and consider the following finite element space on the skeleton mesh $\mathcal{E}_H^H$

$$ \Lambda_H := \{ \mu_H \in \Lambda : \mu_H|_F \in \mathbb{P}_\ell(F), \text{ for all } F \in \mathcal{E}_H^H \}. \qquad (10) $$

Above and hereafter, $\mathbb{P}_m(D)$ denotes the space of polynomials of degree up to $m \ge 0$ defined on $D$.

### 2.3. Space compatibility and interpolation mappings

This section describes the interpolation operators used throughout this work, for which the compatible conditions between $V_h$ and $\Lambda_H$ provided in [7, Section 3] and [22, Lemma 4.1] are required. For the sake of clarity, we remind them below:

**(A2):** Spaces $V_h$ and $\Lambda_H$ in (9) and (10) are compatible as follows:

- For $d=2$:

If $k=\ell$, we make one refinement on $\Xi_H(K)$ for $\ell \ge 2$ and two for $\ell=0, 1$;

If $k=\ell+1$, we make one refinement on $\Xi_H(K)$ for $\ell \ge 0$;

If $k \ge \ell+2$ then no further refinement on $\Xi_H(K)$ for $\ell \ge 0$.

- For $d=3$, we set $k \ge \ell+3$ with no extra refinements on $\Xi_H(K)$ for $\ell \ge 0$.
Consider that assumptions (A1) and (A2) hold. Using [7, Theorem 2] and [22, Lemma 4.1] there exists a mapping $\Pi_h : V \to V_h$ such that

$$ \int_F \Pi_h v \mu_H = \int_F v \mu_H \quad \text{for all } \mu_H \in \Lambda_H \quad \text{and} \quad F \in \mathcal{E}_H^H, \qquad (11) $$

$$ \|\Pi_h v\|_V \le c_F \|v\|_V, \qquad (12) $$

for all $v \in V$ and where $c_F > 0$ does not depend on mesh and physical parameters. Next, let $w \in H^{\ell+2}(\mathcal{P}_H) \cap H_0^1(\Omega)$ be such that $-\varepsilon \nabla w + \frac{1}{2} w \alpha \in H^{\ell+1}(\mathcal{P}_H)^d \cap H(\text{div}; \Omega)$, with $\ell \ge 0$, and $\mu \in \Lambda$ be defined by

$$ \mu|_E := \left( -\varepsilon \nabla w \cdot n^K + \frac{1}{2} w (\alpha \cdot n^K) \right) \Big|_E \quad \text{for all } E \in \mathcal{E}_H^H. $$

Following closely [7, Lemma 3], there exists $\mu_H \in \Lambda_H$ such that

$$ \|\mu - \mu_H\|_{\Lambda} \le C H^{\ell+1} \left| -\varepsilon \nabla w + \frac{1}{2} w \alpha \right|_{\ell+1, \mathcal{P}_H}, \qquad (13) $$

where C is a positive constant independent of the mesh and physical parameters.

## 2.4. Continuous hybrid formulation

The weak solution of problem (1) corresponds to: Find $u \in H_0^1(\Omega)$ such that

$$ (\varepsilon \nabla u, \nabla v)_{\Omega} + \frac{1}{2}(\alpha \cdot \nabla u, v)_{\Omega} - \frac{1}{2}(u, \alpha \cdot \nabla v)_{\Omega} + (\nu u, v)_{\Omega} = (f, v)_{\Omega}, \qquad (14) $$

for all $v \in H_0^1(\Omega)$, where $\nu \in L^{\infty}(\Omega)$ is a nonnegative function defined by

$$ \nu := \frac{1}{2} \nabla \cdot \alpha + \sigma. \qquad (15) $$

The well-posedness of the skew-symmetric formulation (14) follows from the Lax-Milgram's Lemma (see [18, Lemma 2.2]).

The primal hybrid formulation associated to (14) corresponds to: Find $(u, \lambda) \in V \times \Lambda$ such that

$$ \begin{aligned} a(u, v) + \langle \lambda, v \rangle_{\partial \mathcal{P}_H} &= (f, v)_{\mathcal{P}_H} && \text{for all } v \in V, \\ \langle \mu, u \rangle_{\partial \mathcal{P}_H} &= 0 && \text{for all } \mu \in \Lambda, \end{aligned} \qquad (16) $$

where the bilinear form $a: V \times V \to \mathbb{R}$ is

$$ a(u, v) := \sum_{K \in \mathcal{P}_H} a_K(u, v), \qquad (17) $$

and the local bilinear form $a_K: H^1(K) \times H^1(K) \to \mathbb{R}$, for all $K \in \mathcal{P}_H$, is

$$ a_K(u, v) := (\varepsilon \nabla u, \nabla v)_K + \frac{1}{2}(\alpha \cdot \nabla u, v)_K - \frac{1}{2}(u, \alpha \cdot \nabla v)_K + (\nu u, v)_K. \qquad (18) $$

Using (2), and (3), it holds

$$ \|a\| := \sup_{w,v \in V} \frac{a(w,v)}{\|w\|_V \|v\|_V} \le \max \{\varepsilon_{\max}, d_{\Omega} \|\alpha\|_{1,\infty,\Omega}, d_{\Omega}^2 \|\sigma\|_{0,\infty,\Omega}\}. \qquad (19) $$

The next lemma states that the hybrid formulation (16) is well-posed, and its solution coincides with the solution of (14).
**Lemma 1.** Function $u \in H_0^1(\Omega)$ solves (14) iff $(u, \lambda) \in V \times \Lambda$ solves (16) and, for all $K \in \mathcal{P}_H$,

$$ \lambda|_{\partial K} = \sigma \cdot n^K|_{\partial K} \quad \text{where} \quad \sigma := -\epsilon \nabla u + \frac{1}{2} u \alpha. \qquad (20) $$

*Proof.* Let $u \in H_0^1(\Omega)$ be the (unique) solution of (14) and consider the continuous linear functional on $V$

$$ L(v) := (f, v)_{\mathcal{P}_H} - a(u, v). $$

Then $L(\cdot)$ vanishes on $H_0^1(\Omega)$, and it follows from [32], that there exists a unique $\lambda \in \Lambda$ such that $L(v) = \langle \lambda, v \rangle_{\partial \mathcal{P}_H}$, for all $v \in V$. Hence, the pair $(u, \lambda) \in V \times \Lambda$ is the (unique) solution of (16). Now, since $f = \nabla \cdot (-\epsilon \nabla u + \alpha u) + \sigma u$, in a distributional sense, and using integration by parts, we get

$$ L(v) = \sum_{K \in \mathcal{P}_H} \left\langle -\epsilon \nabla u \cdot n^K + \frac{1}{2} u (\alpha \cdot n^K), v \right\rangle_{\partial K} = \langle \lambda, v \rangle_{\partial \mathcal{P}_H}, $$

which proves (20). Conversely, assume that $(u, \lambda) \in V \times \Lambda$ is the solution of (16). The characterization of space $H_0^1(\Omega)$ (c.f. [32]) with the second equation in (16) implies that $u \in H_0^1(\Omega)$. From the first equation in (16), and (6), we get

$$ (\epsilon \nabla u, \nabla v)_{\Omega} + \frac{1}{2}(\alpha \cdot \nabla u, v)_{\Omega} - \frac{1}{2}(u, \alpha \cdot \nabla v)_{\Omega} + (\nu u, v)_{\Omega} = (f, v)_{\Omega} $$

for all $v \in H_0^1(\Omega)$, thus $u$ solves (14). $\square$

## 2.5. Discrete hybrid formulation

Let $V_h$ and $\Lambda_H$ be given in (9) and (10), respectively. The discrete counterpart of (16) reads: Find $(u_h, \lambda_H) \in V_h \times \Lambda_H$

$$ \begin{aligned} a(u_h, v_h) + \langle \lambda_H, v_h \rangle_{\partial \mathcal{P}_H} &= (f, v_h)_{\mathcal{P}_H} && \text{for all } v_h \in V_h, \\ \langle \mu_H, u_h \rangle_{\partial \mathcal{P}_H} &= 0 && \text{for all } \mu_H \in \Lambda_H. \end{aligned} \qquad (21) $$

The next theorem establishes the well-posedness of the discrete hybrid method (21) under the compatible conditions given in (A2). Furthermore, it provides error estimates that take into account the two levels of space discretization.

**Theorem 2.** Assume that (A1) and (A2) hold. Then problem (21) is well-posed

$$ \|u_h\|_V \le \frac{c_B^2}{\epsilon_{min}} \|f\|_{0,\Omega} \quad \text{and} \quad \|\lambda_H\|_{\Lambda} \le c_F \left( 1 + \frac{c_B^2}{\epsilon_{min}} \|a\| \right) \|f\|_{0,\Omega}, \qquad (22) $$

where $c_B$ is a positive constant independent of mesh and physical parameters, but dependent on the shape regularity of $\mathcal{P}_H$, and $c_F$ is given in (12). Moreover, assume the exact solution $u \in H^{k+1}(\mathcal{P}_H) \cap H_0^1(\Omega)$ and $-\epsilon \nabla u + \frac{1}{2} u \alpha \in H^{\ell+1}(\mathcal{P}_H)^d \cap H(\text{div}; \Omega)$, with $k \ge 1$ and $\ell \ge 0$. Then, there are positive constants $C_1$ and $C_2$, independent of mesh parameters but dependent on the shape regularity of $\mathcal{P}_H$, such that

$$ \|u - u_h\|_V \le C_1 \left( H^{\ell+1} \left| -\epsilon \nabla u + \frac{1}{2} u \alpha \right|_{\ell+1, \mathcal{P}_H} + h^k |u|_{k+1, \mathcal{P}_H} \right), \qquad (23) $$

$$ \|\lambda - \lambda_H\|_{\Lambda} \le C_2 \left( H^{\ell+1} \left| -\epsilon \nabla u + \frac{1}{2} u \alpha \right|_{\ell+1, \mathcal{P}_H} + h^k |u|_{k+1, \mathcal{P}_H} \right). \qquad (24) $$

*Proof.* Let $\mathcal{N}_h$ be the following null space

$$ \mathcal{N}_h := \{v_h \in V_h : \langle \mu_H, v_h \rangle_{\partial \mathcal{P}_H} = 0, \text{ for all } \mu_H \in \Lambda_H\}. \qquad (25) $$
Let $v_h \in \mathcal{N}_h$, $E \in \mathcal{E}^H$, and $\mu_E \in \Lambda_H$ be defined by

$$ \mu_E|_{\partial K|_{E'}} = \begin{cases} \mathbf{n}^E \cdot \mathbf{n}^K|_{E'}, & E' = E \\ 0, & E' \neq E \end{cases}, \qquad (26) $$

on each $E' \in \mathcal{E}^H$, then it is clear that

$$ \int_E [v_h] = \langle \mu_E, v_h \rangle_{\partial \mathcal{P}_H} = 0 \quad \text{for all } E \in \mathcal{E}^H. \qquad (27) $$

Using the Poincaré–Friedrich inequality for piecewise $H^1$ functions (c.f. [8, equation (1.3)]) and (26), we have that there exists a positive constant $c_B$, that depends only on the shape of the polytope of $\mathcal{P}_H$ but not on $h$, $H$, or $\mathcal{H}$, such that

$$ \|v_h\|_V \le c_B |v_h|_{1, \mathcal{P}_H}, \quad \text{for all } v_h \in \mathcal{N}_h. \qquad (28) $$

Thus,

$$ a(v_h, v_h) \ge \varepsilon_{\min} |v_h|_{1, \mathcal{P}_H}^2 \ge \frac{\varepsilon_{\min}}{c_B^2} \|v_h\|_V^2 \quad \text{for all } v_h \in \mathcal{N}_h. \qquad (29) $$

Moreover, under the assumption (A2), we use (11)–(12) to prove

$$ c_F^{-1} \|\mu_H\|_{\Lambda} = c_F^{-1} \sup_{v \in V} \frac{\langle \mu_H, v \rangle_{\partial \mathcal{P}_H}}{\|v\|_V} \le \sup_{v_h \in V_h} \frac{\langle \mu_H, v_h \rangle_{\partial \mathcal{P}_H}}{\|v_h\|_V}. \qquad (30) $$

Hence, using [18, Theorem 2.34] we obtain (22) and the well-posedness of problem (21). Next from [18, Lemma 2.44], we obtain the best approximation results

$$ \begin{aligned} \|u - u_h\|_V &\le c_1 \inf_{v_h \in V_h} \|u - v_h\|_V + c_2 \inf_{\rho_H \in \Lambda_H} \|\lambda - \rho_H\|_{\Lambda}, \\ \|\lambda - \lambda_H\|_{\Lambda} &\le c_3 \inf_{v_h \in V_h} \|u - v_h\|_V + c_4 \inf_{\rho_H \in \Lambda_H} \|\lambda - \rho_H\|_{\Lambda}, \end{aligned} \qquad (31) $$

where, $c_1 := (1 + \frac{c_B^2}{\varepsilon_{\min}^2}) \|a\|$) (1 + c_F)$, $c_2 := \frac{c_B^2}{\varepsilon_{\min}^2}$, $c_3 := c_1 \|a\| c_F$, and $c_4 := 1 + c_F + c_2 \|a\| c_F$. The error estimates (24) and (23) arise using (13), (31), and classical Lagrange interpolation estimates. $\square$

**Remark 1** (Convergence with fixed $\mathcal{H}$). Under conditions (A1)–(A2), the error estimates (23) and (24) indicate that convergence can be achieved by fixing the coarse mesh (i.e. $\mathcal{H}$ is fixed) and refining the facet meshes (e.g. $H \to 0$). In this case, no assumption about the shape regularity of $\mathcal{P}_H$ is necessary. On the other hand, if we look at the standard convergence with respect to $\mathcal{H}$, then we assume that $\{\mathcal{P}_H\}_{\mathcal{H}>0}$ satisfies the shape regularity conditions given in [8, Corollary 6.3 and 7.3] in order to preserve the optimality of the error estimates (23)-(24).

**Remark 2** (Second level numerical pollution). Due to the relationship between the polynomial degrees $\ell$ and $k$ given in (A2), we note that the main error is of the order $O(H^{\ell+1})$ unless $\ell = k$ is used where special attention must be given to refining the second level so that the error associated with the second level of discretization (i.e. $O(h^k)$) does not pollute the first one (i.e. $O(H^{\ell+1})$).

### 3. The Multiscale Hybrid-Mixed Method

The formal construction of the two-level MHM method starts from the discrete hybrid formulation given in (21) using a decomposition of the space $V_h$. The objective is to locate calculations through element-by-element RAD problems that are responsible for calculating the multiscale basis. The occurrence of such spatial decomposition depends on the advection and reaction coefficients to ensure that these local problems are well-posed. We describe the strategy in this section, propose and analyze the two-level MHM methods.
### 3.1. Space decomposition

Let $V_0$ be the piecewise constant polynomial space

$$V_0 := \{v_0 \in V_h : v_0|_K \in \mathbb{P}_0(K) \text{ for all } K \in \mathcal{P}_H\}, \quad (32)$$

and consider the subspace $W_0$ of $V_0$

$$W_0 := \{v_0 \in V_0 : v_0|_K = 0 \text{ for all } K \in \mathcal{P}_H^0\}, \quad (33)$$

where

$$\mathcal{P}_H^0 := \{K \in \mathcal{P}_H : \text{ there exists } \xi_h \in V_h(K) \text{ s.t. } a_K(1_K, \xi_h) > 0\}. \quad (34)$$

We decompose the discrete space $V_h$ given in (9) by

$$V_h = W_0 \oplus \tilde{W}_h, \quad (35)$$

with respect to the inner product $(\cdot, \cdot)_V$, and then $\tilde{W}_h$ corresponds to

$$\tilde{W}_h := \{\tilde{w}_h \in V_h : \tilde{w}_h|_K \in L_0^2(K), \text{ for all } K \in \mathcal{P}_H \setminus \mathcal{P}_H^0\}, \quad (36)$$

where $L_0^2(K)$ is the space of functions in $L^2(K)$ with zero average in $K$.

We observe then that the dimension of (33) depends on the physical coefficients. For instance, $W_0 = V_0$ in the pure diffusion case, and $W_0$ is trivial when the reaction coefficient is present.

**Remark 3 (Characterisation of $W_0$).** Notice that $W_0$ coincides with the null space

$$\{v_0 \in V_h : a(v_0, v_h) = 0 \text{ for all } v_h \in V_h\}, \quad (37)$$

and also is equal to

$$\{v_0 \in V_h : a(v_h, v_0) = 0 \text{ for all } v_h \in V_h\}. \quad (38)$$

This follows by first observing that $W_0$ is a subspace of (37). For the opposite inclusion, take $v_0$ in (37) and note that $a(v_0, v_0) \ge \varepsilon_{min}|v_0|_1^2_{\mathcal{P}_H}$. Therefore, $v_0|_K \in \mathbb{P}_0(K)$ for all $K \in \mathcal{P}_H$ which implies that (37) is a subspace of $W_0$. Writing $v_0|_K = 1_K c_K$, where $c_K \in \mathbb{R}$, we have $c_K a_K(1_K, \xi_h) = 0$ for all $\xi_h \in V_h(K)$, and therefore, $c_K = 0 \Rightarrow v_0|_K = 0$ if $K \in \mathcal{P}_H^0$ which implies $v_0 \in W_0$. Then (37) is also a subspace of $W_0$. The equivalence (37) and (38) follows by noting that if there exists $\xi_h \in V_h(K)$ such that $a_K(1_K, \xi_h) > 0$, then we can propose $\xi_h^* \in V_h(K)$ using $\xi_h$ such that $a_K(\xi_h^*, 1_K) > 0$.

### 3.2. The method

Let $T_h \in \mathcal{L}(\Lambda, \tilde{W}_h)$ and $\hat{T}_h \in \mathcal{L}(L^2(\Omega), \tilde{W}_h)$ be two mappings defined such that, given $\rho \in \Lambda$ and $q \in L^2(\Omega)$, the functions $T_h \rho$ and $\hat{T}_h q$ satisfy

$$a(T_h \rho, \tilde{v}_h) = - \langle \rho, \tilde{v}_h \rangle_{\partial \mathcal{P}_H} \quad \text{for all } \tilde{v}_h \in \tilde{W}_h, \quad (39)$$

$$a(\hat{T}_h q, \tilde{v}_h) = (q, \tilde{v}_h)_{\mathcal{P}_H} \quad \text{for all } \tilde{v}_h \in \tilde{W}_h. \quad (40)$$

The operators $T_h$ and $\hat{T}_h$ are well-defined element-by-element as the bilinear form $a(\cdot, \cdot)$ given in (17)–(18) induces a bijective operator on $\tilde{W}_h$. It follows from [18, Proposition 2.21] that there is a positive constant $\gamma$ such that

$$\inf_{w_h \in \tilde{W}_h} \sup_{v_h \in \tilde{W}_h} \frac{a(w_h, v_h)}{\|w_h\|_V \|v_h\|_V} \ge \gamma, \quad (41)$$

and then, from (39)–(41) the mappings are bounded with $\|T_h\| \le \frac{1}{\gamma}$ and $\|\hat{T}_h\| \le \frac{d_\Omega}{\gamma}$.
The two-level MHM method writes: Find $(\lambda_H, u_0^H) \in \Lambda_H \times W_0$ such that

$$ -\langle \mu_H, T_h \lambda_H \rangle_{\partial \mathcal{P}_H} - \langle \mu_H, u_0^H \rangle_{\partial \mathcal{P}_H} = \langle \mu_H, \hat{T}_h f \rangle_{\partial \mathcal{P}_H} \quad \text{for all } \mu_H \in \Lambda_H, \\ -\langle \lambda_H, v_0 \rangle_{\partial \mathcal{P}_H} = -(f, v_0)_{\mathcal{P}_H} \quad \text{for all } v_0 \in W_0. \qquad (42) $$

The local problems $T_h\lambda_H$ and $\hat{T}_h f$ intend to approximate the fine and under-mesh structure of $u$ using the coarse and global MHM's solution $(\lambda_H, u_0^H)$. As such, the exact solution $u$ is approximated by

$$ u \sim u_{H,h} := u_0^H + T_h \lambda_H + \hat{T}_h f. \qquad (43) $$

As a result, we have the characterization

$$ u_0^H|_K = \frac{1}{|K|} \int_K u_{H,h} d\boldsymbol{x} \sim \frac{1}{|K|} \int_K u d\boldsymbol{x} \quad \text{for all } K \in \mathcal{P}_H \setminus \mathcal{P}_H^0, $$

and $\lambda_H$ represents the approximate flux on the elements boundary which is in local equilibrium with external force, i.e.,

$$ \int_{\partial K} \lambda_H ds = \int_K f d\boldsymbol{x} \quad \text{for all } K \in \mathcal{P}_H \setminus \mathcal{P}_H^0. $$

**Remark 4 (Variants of the method (42)).** The two-level MHM method (42) recovers previous MHM methods for the RAD model as a function of the structure of the coarse $W_0$ space. Notably:

* $W_0 = \{0\}$ case: We find the method originally proposed in [25], i.e., find $\lambda_H \in \Lambda_H$

$$ -\langle \mu_H, T_h \lambda_H \rangle_{\partial \mathcal{P}_H} = \langle \mu_H, \hat{T}_h f \rangle_{\partial \mathcal{P}_H} \quad \text{for all } \mu_H \in \Lambda_H. \qquad (44) $$

This case occurs, for instance, when $\nu \ge \nu_0 > 0$ in $\Omega$, or when $\nu = 0$ and there exists $\xi_K \in V_h(K)$ such that $(1_K, \boldsymbol{\alpha} \cdot \nabla \xi_K)_K > 0$ for all $K \in \mathcal{P}_H$;

* $W_0 = V_0$ case: We find the method originally proposed in [2, 7], that is, the method (42) with $W_0$ replaced by $\bar{V}_0$. This scenario occurs in the pure diffusive case, for example. Furthermore, it occurs when $\nu = 0$ and $\boldsymbol{\alpha}$ is such that $(\boldsymbol{\alpha}, \nabla v_h)_K = 0$ for all $v_h \in V_h$ and $K \in \mathcal{P}_H$. A typical example is when the advection field satisfies $\nabla \cdot \boldsymbol{\alpha}|_K = 0$ and $\boldsymbol{\alpha} \cdot \boldsymbol{n}^K|_{\partial K} = 0$ for all $K \in \mathcal{P}_H$.

### 3.3. Well-posedness and convergence

We follow the strategy proposed in [26, 6] to analyze the well-posedness and convergence of the MHM method (42). Specifically, and unlike the original MHM work, we directly estimate the error in (43), without first estimating the error between the flow variable $\lambda$ and $\lambda_H$ (see [13] for a similar strategy).

**Theorem 3.** Assume that (A1) and (A2) hold. Then the method (42) has a unique solution $(\lambda_H, u_0^H) \in \Lambda_H \times W_0$. Moreover, if $u_h \in V_h$ is obtained by solving the discrete hybrid problem (21), then

$$ u_h = u_{H,h} := u_0^H + T_h \lambda_H + \hat{T}_h f, $$

where $T_h$ and $\hat{T}_h$ are given in (39)–(40). Also the error estimates (23)–(24), with $u_h = u_{H,h}$ given in (43), hold.

*Proof.* Let $(\rho_H, w_0)$ be an element of $\Lambda_H \times W_0$ such that

$$ -\langle \mu_H, w_0 \rangle_{\partial \mathcal{P}_H} - \langle \mu_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} = 0 \quad \text{for all } \mu_H \in \Lambda_H, \qquad (45) $$

$$ -\langle \rho_H, v_0 \rangle_{\partial \mathcal{P}_H} = 0 \quad \text{for all } v_0 \in W_0. \qquad (46) $$
Testing (45) with $(\mu_H, v_0) = (\rho_H, w_0)$ and subtracting (46) from (45), we get

$$0 = - \langle \rho_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} = a(T_h \rho_H, T_h \rho_H) \geq \varepsilon_{\min} |T_h \rho_H|_{1, \mathcal{P}_H}^2 \Rightarrow T_h \rho_H \in V_0.$$

On the other hand, (45) implies that $w_0 + T_h \rho_H \in \mathcal{N}_h$ and using the definition of $\tilde{W}_h$ given in (36), and (28) we get

$$\|w_0\|_V^2 + \|T_h \rho_H\|_V^2 = \|w_0 + T_h \rho_H\|_V^2 \leq c_B^2 |T_h \rho_H|_{1, \mathcal{P}_H}^2 = 0 \Rightarrow w_0 = T_h \rho_H = 0.$$

We conclude $\rho_H = 0$ from the injectivity of $T_h$ using (A2), thus the well-posedness of (42) follows. Next, notice that the (unique) solution $u_h$ of (21) belongs to the null space $\mathcal{N}_h$ and that the discrete hybrid formulation (21) is equivalent to find $u_h \in \mathcal{N}_h$

$$a(u_h, v_h) = (f, v_h)_{\mathcal{P}_H} \quad \text{for all } v_h \in \mathcal{N}_h.$$

Also, observe that $u_{H,h}$ given in (43) belongs to $\mathcal{N}_h$ by taking $v_0 = 0$ in (42). Let $v_h$ be an element of $\mathcal{N}_h$ decomposed as $v_h = v_0 + \tilde{v}_h$, with $v_0 \in W_0$ and $\tilde{v}_h \in \tilde{W}_h$. Using the definitions of the space $W_0$ in (37), mappings $T_h$ and $\tilde{T}_h$ in (39)–(40), and the second equation in (42), we get

$$
\begin{align*}
a(u_{H,h}, v_h) &= a(u_0^H + T_h \lambda_H + \tilde{T}_h f, v_0) + a(u_0^H + T_h \lambda_H + \tilde{T}_h f, \tilde{v}_h) \\
&= a(T_h \lambda_H + \tilde{T}_h f, \tilde{v}_h) \\
&= -\langle \lambda_H, \tilde{v}_h \rangle_{\partial \mathcal{P}_H} + (f, \tilde{v}_h)_{\mathcal{P}_H} \\
&= -\langle \lambda_H, \tilde{v}_h \rangle_{\partial \mathcal{P}_H} - \langle \lambda_H, v_0 \rangle_{\partial \mathcal{P}_H} + (f, v_0)_{\mathcal{P}_H} + (f, \tilde{v}_h)_{\mathcal{P}_H} \\
&= -\langle \lambda_H, v_h \rangle_{\partial \mathcal{P}_H} + (f, v_h)_{\mathcal{P}_H} = (f, v_h)_{\mathcal{P}_H}.
\end{align*}
$$

Hence, $u_h = u_{H,h}$ by the uniqueness of the solution of (21). Finally, the error estimates (23) and (24) follow from Theorem 2. $\square$

**Remark 5** (Strong Dirichlet boundary conditions). The numerical analysis of the MHM method with weakly imposed Dirichlet boundary conditions can be easily extended to the case where the Dirichlet condition is strongly imposed (c.f. [25]). To this end, we redefine some geometric and functional settings. First, we distinguish between internal and external elements as follows

$$\mathcal{P}_H^{\text{ext}} := \{K \in \mathcal{P}_H : |\partial K \cap \partial \Omega| > 0\} \quad \text{and} \quad \mathcal{P}_H^{\text{int}} := \mathcal{P}_H \setminus \mathcal{P}_H^{\text{ext}},$$

and denote the internal boundaries of $K \in \mathcal{P}_H$ by $\delta K := \partial K \setminus \partial \Omega$, note that $\delta K = \partial K$ for all $K \in \mathcal{P}_H^{\text{int}}$. Next, we define the space of Lagrange multipliers as

$$\bar{\Lambda} := \{ \boldsymbol{q} \cdot \boldsymbol{n}^K |_{\delta K}, \text{ for all } K \in \mathcal{P}_H : \boldsymbol{q} \in H(\text{div}; \Omega) \},$$

equipped it with the norm

$$\|\mu\|_{\bar{\Lambda}} := \inf \{ \|\boldsymbol{\sigma}\|_{\text{div}} : \boldsymbol{\sigma} \in H(\text{div}; \Omega) \text{ and } \boldsymbol{\sigma} \cdot \boldsymbol{n}^K|_{\delta K} = \mu, \forall K \in \mathcal{P}_H \}.$$

We also define $\bar{V} := \{v \in V : v|_{\partial \Omega} = 0\}$, the duality product $\langle \cdot, \cdot \rangle_{\partial \mathcal{P}_H} : \bar{\Lambda} \times \bar{V} \to \mathbb{R}$, by

$$\langle \mu, v \rangle_{\partial \mathcal{P}_H} := \sum_{K \in \mathcal{P}_H^{\text{int}}} \langle \mu, v \rangle_{H^{-\frac{1}{2}}(\partial K) \times H^{\frac{1}{2}}(\partial K)} + \sum_{K \in \mathcal{P}_H^{\text{ext}}} \langle \mu, v \rangle_{H^{-\frac{1}{2}}(\delta K) \times H_{00}^{\frac{1}{2}}(\delta K)},$$

for all $(\mu, v) \in \bar{\Lambda} \times \bar{V}$. The Sobolev spaces $H^{-\frac{1}{2}}(D)$, $H^{\frac{1}{2}}(D)$ and $H_{00}^{\frac{1}{2}}(D)$ have their usual meaning (c.f [23] and [34, pp. 121]). Observe that

$$\langle \bar{\mu}, v \rangle_{\partial \mathcal{P}_H} = \langle \mu, v \rangle_{\partial \mathcal{P}_H} \quad \forall \mu \in \Lambda, \forall \bar{\mu} \in \bar{\Lambda}, \forall v \in \bar{V}. \qquad (47)$$
The identity (47) implies that the use of the dual product $\langle \cdot, \cdot \rangle_{\partial \mathcal{P}_H}$ can be mimicked by $\langle \cdot, \cdot \rangle_{\overline{\partial \mathcal{P}_H}}$ making trivial adaptations. Notably, following [31, Lemma 1], we can write

$$H_0^1(\Omega) := \{v \in \bar{V} : \langle \mu, v \rangle_{\overline{\partial \mathcal{P}_H}} = 0, \text{ for all } \mu \in \bar{\Lambda}\},$$

and consider the following hybrid variational problem as the starting point of the construction of the MHM method: Find $(u, \lambda) \in \bar{V} \times \bar{\Lambda}$ such that

$$\left\{
\begin{array}{ll}
a(u, v) + \langle \lambda, v \rangle_{\overline{\partial \mathcal{P}_H}} = (f, v)_{\mathcal{P}_H} & \text{for all } v \in \bar{V} \\
\langle \mu, u \rangle_{\overline{\partial \mathcal{P}_H}} = 0 & \text{for all } \mu \in \bar{\Lambda}.
\end{array}
\right.
\qquad (48)$$

The proof of the well-posedness of the problem (48) as well as its equivalence with the problem (14) follow closely those given in Lemma 1 and Theorem 2, respectively. Then, the other theoretical results remain valid. Strong conditions are used in the computational simulations in [25, Fig. 6] and in Section 5.

## 4. Condition number Analysis

Solving an algebraic system of equations accurately in a computationally efficient way involves theoretical knowledge of the condition number of the matrices. This is more dramatic when iterative solvers are employed, which are attractive for solving the global MHM system that is coupled, unlike the local problems that are already solved in an “embarrassingly parallel” way. It turns out that the local and global linear systems that arise in the two-level MHM method change from mixed to positive defined forms as a function of physical and mesh parameters, with an immediate impact on the condition numbers of the matrices. How physical and mesh parameters impact the condition number of global and local linear systems is the topic of this section.

For completeness, we recall some fundamental concepts related to condition numbers used in the rest of this section. Let $\mathcal{M} \in \mathbb{R}^{N \times N}$ be a non-singular matrix. The 2-norm of a vector $\mathbf{x} \in \mathbb{R}^N$ is defined as $\|\mathbf{x}\|_2 := (x_1^2 + \cdots + x_N^2)^{\frac{1}{2}}$. The 2-norm of a matrix $\mathcal{M}$ is given by $\|\mathcal{M}\|_2 := \sup_{\mathbf{x} \in \mathbb{R}^N} \frac{\|\mathcal{M}\mathbf{x}\|_2}{\|\mathbf{x}\|_2}$, and the condition number of $\mathcal{M}$ is defined by

$$\mathrm{cond}(\mathcal{M}) := \|\mathcal{M}\|_2 \|\mathcal{M}^{-1}\|_2. \qquad (49)$$

Consider a finite-dimensional linear space $Z$ equipped with the norm $\| \cdot \|_Z$. We set $\{z_1, \dots, z_N\}$ as a basis for $Z$. Then, given $x := x_1 z_1 + \cdots + x_N z_N \in Z$, we denote $\mathbf{x} = [x_1 \cdots x_N]^T \in \mathbb{R}^{N \times 1}$ and the following equivalence of norms holds

$$c_1 \|\mathbf{x}\|_2 \le \|\mathbf{x}\|_Z \le c_2 \|\mathbf{x}\|_2, \qquad (50)$$

where $c_1, c_2$ are positive constants dependent only on $N$ and the choice of the basis $\{z_1, \dots, z_N\}$. Also we consider a bounded bilinear form $b: Z \times Z \to \mathbb{R}$ and its induced norm $\|b\|$, and assume that the following inf-sup condition holds

$$\inf_{x \in Z} \sup_{y \in Z} \frac{b(x, y)}{\|x\|_Z \|y\|_Z} \ge \gamma_b > 0.$$

Then, given a matrix $\mathcal{M}$ with entries $\mathcal{M}_{ij} := b(z_j, z_i)$, for $i, j = 1, \dots, N$, the condition number of $\mathcal{M}$ in (49) is bounded as follows

$$\mathrm{cond}(\mathcal{M}) \le \left(\frac{c_2}{c_1}\right)^2 \frac{\|b\|}{\gamma_b}. \qquad (51)$$

In what follows, we use (51) to estimate the condition number of the local and global matrices from the two-level MHM methods.
## 4.1. Settings

To avoid unnecessary technicalities and simplify the presentation, we assume the following constraint on the advection field:

(A3): There exists $\xi_h \in V_h$, with $\max_{\kappa \in \mathcal{T}_h} \|\xi_h\|_{1,\infty,\kappa} \le 1$, and a positive constant $\omega$, independent of mesh parameters, such that

$$ -\frac{1}{2}\alpha \cdot \nabla \xi_h + \nu \xi_h \ge \omega, \quad \text{in } \Omega_0, \qquad (52) $$

where $\Omega_0$ is the interior of the union of polytopes $\bigcup_{K \in \mathcal{P}_H^0} K$.

**Remark 6.** Note that the (A3) assumption is not very restrictive. In fact, if $\nu|_{\Omega_0} \ge \nu_0 > 0$, the condition (52) is valid taking $\xi_h = 1$ in $\Omega_0$. The zero reaction situation is more involved. In this case, we can assume that a continuous counterpart of (A3) is valid to prove (52). Specifically, let us assume that

$$ \|\xi\|_{1,\infty,\Omega_0} \le 1 \quad \text{and} \quad -\frac{1}{2}\alpha \cdot \nabla \xi + \nu \xi \ge \bar{\omega} \quad \text{in } \Omega_0, \qquad (53) $$

for a given $\xi \in W^{m,\infty}(\Omega_0)$, $m \ge 2$, and positive constant $\bar{\omega} > 0$ independent of mesh parameters. Observe that condition (53) is fulfilled, for example, if the advective field has no closed curves or stationary points in $\Omega_0$ (see [4] for details). Now, from (53), we prove that (A3) is valid. Indeed, let $\mathcal{J}_h^k: W^{k+1,\infty}(\Omega_0) \to V_h$ be the Lagrange interpolation operator, with $1 \le k \le m-1$. Using [18, Corollary 1.109], we have

$$ \|v - \mathcal{J}_h^k v\|_{0,\infty,\Omega_0} + h \max_{\kappa \in \mathcal{J}_h} |v - \mathcal{J}_h^k v|_{1,\infty,\kappa} \le C_k h^{k+1} |v|_{k+1,\infty,\Omega_0}. \qquad (54) $$

Let us define $\omega_1 := 1 + C_k h_0 (h_0 + 1) |\xi|_{k+1,\infty,\Omega_0}$, with $h_0 := \left(\frac{\bar{\omega}}{2C_k |\xi|_{k+1,\infty,\Omega_0}}\right)^{\frac{1}{k}}$. Then, condition (A3) holds with $\omega := \frac{\bar{\omega}}{2\omega_1}$ and $\xi_h := \frac{1}{\omega_1} \mathcal{J}_h^k(\xi)$ for $h \in (0, h_0]$. Note that (A3) implies the existence of $\xi_K^* \in V_h(K)$ such that

$$ \frac{a_K(1_K, \xi_K^*)}{|K|^{\frac{1}{2}} \|\xi_K^*\|_{V(K)}} \ge \omega \quad \text{for all } K \in \mathcal{P}_H^0. \qquad (55) $$

Additionally, in this section, we assume the following mesh constraints:

(i) The families of coarse partitions {$\mathcal{P}_H$}$_{H>0}$, facet partitions {$\mathcal{E}_H^H$}$_{H>0}$ and the local meshes {$\mathcal{T}_h$}$_{h>0}$ are quasi-uniform. Notably, for the family of partitions {$\mathcal{P}_H$}$_{H>0}$ we assume the existence of $\bar{\omega} > 0$ such that $\mathcal{H}_K \ge \bar{\omega} \mathcal{H}$ for all $K \in \mathcal{P}_H$ and for all $\mathcal{H} > 0$;

(ii) The elements $K \in \mathcal{P}_H$ are such that the local Poincaré-Wirtinger inequality holds

$$ \|v_h\|_{V(K)}^2 \le (1 + c_K^2) \|\nabla v_h\|_{0,K}^2 \le c_W^2 \|\nabla v_h\|_{0,K}^2, \qquad (56) $$

for all $v_h \in \tilde{V}_h(K)$ and where $c_W > 0$ does not depend on $K$ or $\mathcal{H}$. For instance, this is the case when $K$ is convex;

(iii) There exists $\vartheta > 0$ such that for all $\mathcal{H} > 0$ the polytopes $K \in \mathcal{P}_H$ are star-shaped with respect to a ball of radius $\zeta_K \ge \vartheta \mathcal{H}_K$ and center $x_K$.

**Remark 7.** Constraints above implies two uniform bounds with respect to $K \in \mathcal{P}_H$ and for all $\mathcal{H} > 0$. First, the quasi-uniformity of the family of meshes {$\mathcal{T}_h$}$_{h>0}$ in (i) gives the inverse inequality in [18, Corollary 1.141] that reads as follows

$$ \|\nabla v_h\|_{0,K} \le c_I h^{-1} \|v_h\|_{0,K}, \quad \text{for all } v_h \in V_h(K), \qquad (57) $$

with $c_I > 0$ not depending on $\mathcal{H}$ or $h$. On the other hand, from assumption (iii) and as a straightforward consequence of [9, eq. 2.18] we arrive at the existence of $c_t > 0$ depending only on $\vartheta$ such that

$$ \|v\|_{0,\partial K} \le c_t \mathcal{H}^{-\frac{1}{2}} \|v\|_{V(K)}, \quad \text{for all } v \in H^1(K). \qquad (58) $$
On the other hand, as mentioned at the beginning of this section, constants $c_1$ and $c_2$ in bound (51) arise from the norm equivalence (50) and depend on the choice of the basis $\{z_1, \dots, z_N\}$. Therefore, it is necessary to set the basis for the spaces $V_h$, $V_0$, and $\Lambda_H$. This is done as follows:

* Let $K \in \mathcal{P}_H$ be an arbitrary element. The chosen basis for $V_h(K)$ is the set $\{\phi_1, \dots, \phi_{\dim V_h(K)}\}$ composed of Lagrange polynomials of degree up to $k$, defined on the triangulation $\mathcal{T}_h^K$.

* For $V_0$, we set the basis $\{1_{K_1}, \dots, 1_{K_{\dim V_0}}\}$, where $1_{K_i}$ is the characteristic function for $K_i \in \mathcal{P}_H$, with $1 \le i \le \#\mathcal{P}_H$.

* To set the basis $\{\psi_1, \dots, \psi_{\dim \Lambda_H}\}$ for $\Lambda_H$, we consider an arbitrary $K \in \mathcal{P}_H$ and collect the basis functions with support on $\partial K$ into the set $\{\psi_1^K, \dots, \psi_{M_K}^K\}$. We define their elements as $\psi_i^K|_E = (n^E \cdot n^K) \varphi_i|_E$ for $E \in \mathcal{E}^H$ such that $E \subset \partial K$, where $\{\varphi_1^K, \dots, \varphi_{M_K}^K\}$ are the Lagrange polynomials of degree up to $\ell$, defined on the skeleton mesh $\mathcal{E}_h^H(K)$.

Hereafter, the basis for the spaces $V_h$, $V_0$, and $\Lambda_H$ are as defined above.

## 4.2. Equivalence of Norms

To use (51) estimation in practice, it is necessary to estimate the constants $c_1$ and $c_2$ in (50) in relation to the norms induced by the local and global weak forms of the MHM method. This is covered in this section and summarized in the following lemma.

**Lemma 4.** Consider fixed basis for the spaces $V_h$, $V_0$, and $\Lambda_H$. The following equivalence of norms holds:

(1) Let $K \in \mathcal{P}_H$, under a quasi-uniformity assumption on $\{\mathcal{T}_h\}_{h>0}$ there exist positive constants $c_1$ and $c_2$, independent of $h$, such that

$$c_1 h^{\frac{d}{2}} \|v_h\|_2 \le \|v_h\|_{V(K)} \le c_2 h^{\frac{d}{2}-1} \|v_h\|_2 \quad \text{for all } v_h \in V_h(K). \quad (59)$$

(2) Under a quasi-uniformity assumption on $\{\mathcal{P}_H\}_{H>0}$ there exist positive constants $c_1$ and $c_2$, independent of $H$, such that

$$c_1 \mathcal{H}^{\frac{d}{2}} \|v_0\|_2 \le \|v_0\|_V \le c_2 \mathcal{H}^{\frac{d}{2}} \|v_0\|_2 \quad \text{for all } v_0 \in V_0. \quad (60)$$

(3) Under a quasi-uniformity assumption on $\{\mathcal{E}_H^H\}_{H>0}$ there exist positive constants $c_1$ and $c_2$, independent on $H$ and $H$, such that

$$c_1 H^{\frac{d-1}{2}} \|\mu_H\|_2 \le \left( \sum_{E \in \mathcal{E}_H^H} \|\mu_H\|_{0,E}^2 \right)^{\frac{1}{2}} \le c_2 H^{\frac{d-1}{2}} \|\mu_H\|_2 \quad \text{for all } \mu_H \in \Lambda_H. \quad (61)$$

(4) Under a quasi-uniformity assumption on $\{\mathcal{E}_H^H\}_{H>0}$ and assuming that $\mathcal{P}_H$ can be generated from one reference polytopal through an affine transformation, there exist positive constants $c_1$ and $c_2$, independent of $H$ and $H$, such that

$$c_1 H^{\frac{d}{2}} \|\mu_H\|_2 \le \|\mu_H\|_{\Lambda} \le c_2 \frac{H^{\frac{d-1}{2}}}{\mathcal{H}^{\frac{d}{2}}} \|\mu_H\|_2 \quad \text{for all } \mu_H \in \Lambda_H. \quad (62)$$

*Proof.* The item (1) is a straightforward application of [18, Lemma 9.7] joint to the inverse inequality (57). Indeed, [18, Lemma 9.7] gives the existence of positive constants $C_1$ and $C_2$ not depending on $K$ or $h$ and such that

$$C_1^2 h^d \|v_h\|_2 \le \|v_h\|_{0,K}^2 \le C_2^2 h^d \|v_h\|_2,$$
then, applying (57) in the upper bound we obtain

$$C_1^2 \frac{h^d}{d_{\Omega}^2} \|\mathbf{v}_h\|_2^2 \le \|v_h\|_{V(K)}^2 \le C_2^2 \left(c_I^2 + \frac{h^2}{d_{\Omega}^2}\right) h^{d-2} \|\mathbf{v}_h\|_2^2.$$

For item (2), we consider $v_0 \in V_0$ and characterize it as $v_0 = \sum_{K \in \mathcal{P}_H} v_K 1_K$, with $v_K \in \mathbb{R}$ for all $K \in \mathcal{P}_H$, thus $\|v_0\|_2^2 := \sum_{K \in \mathcal{P}_H} v_K^2$. We immediately get the desired result due to the quasi-uniformity of the family $\{\mathcal{P}_H\}_{H>0}$.

The proof of item (3) starts by denoting the basis for $\Lambda_H$ as $\mathbf{B} := \{\psi_1, \dots, \psi_{\dim \Lambda_H}\}$, and define $\mathbf{B}_E := \{\psi^E \in \mathbf{B} : \text{supp}(\psi^E) \subseteq E\}$ for $E \in \mathcal{E}^H$. Then we observe that $\mu_H \in \Lambda_H$ can be expressed as

$$\mu_H = \sum_{E \in \mathcal{E}^H} \sum_{i=1}^{\#\mathbf{B}_E} c_i^E \psi_i^E,$$

with $c_1^E, \dots, c_{\#\mathbf{B}_E}^E \in \mathbb{R}$, and, $\psi_1^E, \dots, \psi_{\#\mathbf{B}_E}^E \in \mathbf{B}_E$. Following the proof of [18, Lemma 9.7] there exist constants $c_1 > 0$ and $c_2 > 0$, independent of $H$ or $|E|$, such that

$$c_1 H^{d-1} \sum_{i=1}^{\#\mathbf{B}_E} (c_i^E)^2 \le \|\mu_H\|_{0,E}^2 \le c_2 H^{d-1} \sum_{i=1}^{\#\mathbf{B}_E} (c_i^E)^2,$$

and summing over $\mathcal{E}^H$ we arrive to (61).

We begin the proof of item (4) by considering a simplex of reference $\hat{F} \subset \mathbb{R}^{d-1}$ and define an affine isomorphism $S_F : \hat{F} \to F$ for each $F \in \mathcal{E}_H^H$. In addition we consider a bubble function $\hat{p} = l_1 \cdots l_d \in \mathbb{P}_d(\hat{F})$, where $l_i : \hat{F} \to \mathbb{R}$ is the $i$-th linear Lagrange polynomial with $i \in \{1, \dots, d\}$. As $\hat{\mu} \mapsto \|\hat{p}^{\frac{1}{2}} \hat{\mu}\|_{0,\hat{F}}$, for all $\hat{\mu} \in \mathbb{P}_\ell(\hat{F})$, is a norm on $\mathbb{P}_\ell(\hat{F})$, thus we can establish the existence of positive constants $\hat{c}_1$ and $\hat{c}_2$, depending only on $\ell$, $\hat{p}$, and $\hat{F}$, such that

$$\hat{c}_1 \|\hat{\mu}\|_{0,\hat{F}} \le \|\hat{p}^{\frac{1}{2}} \hat{\mu}\|_{0,\hat{F}} \le \hat{c}_2 \|\hat{\mu}\|_{0,\hat{F}}, \quad \text{for all } \hat{\mu} \in \mathbb{P}_\ell(\hat{F}). \qquad (63)$$

Then, for any $F \in \mathcal{E}_H^H$, the polynomial $\hat{p} \circ S_F^{-1} : F \to \mathbb{R}$ is also a bubble function and, moreover, by setting $\hat{\mu} = \mu \circ S_F \in \mathbb{P}_\ell(\hat{F})$ for all $\mu \in \mathbb{P}_\ell(F)$, and using standard scaling arguments, we arrive at

$$\hat{c}_1 \|\mu\|_{0,F} \le \|(\hat{p} \circ S_F^{-1})^{\frac{1}{2}} \mu\|_{0,F} \le \hat{c}_2 \|\mu\|_{0,F}, \quad \text{for all } \mu \in \mathbb{P}_\ell(F), \qquad (64)$$

where $\hat{c}_1$ and $\hat{c}_2$ are the constants in (63). Let $K \in \mathcal{P}_H$ and consider the triangulation $\Xi_H(K)$ defined in Assumption (A1). Also, let $\hat{\kappa} \subset \mathbb{R}^d$ be a simplex of reference and $R_\kappa : \hat{\kappa} \to \kappa$ an affine isomorphism, for all $\kappa \in \Xi_H(K)$. Then, by setting $\hat{v} := v \circ R_\kappa \in H^1(\hat{\kappa})$ for all $v \in H^1(\kappa)$, and using standard scaling arguments, we prove the existence of positive constants $\hat{c}_3$ and $\hat{c}_4$, independent of $H$, such that

$$\hat{c}_3 \|\hat{v}\|_{1,\hat{\kappa}} \le \left( |v|_{1,\kappa}^2 + \frac{1}{d_{\Omega}^2} \|v\|_{0,\kappa}^2 \right)^{\frac{1}{2}} \le \hat{c}_4 H^{\frac{d}{2}-1} \|\hat{v}\|_{1,\hat{\kappa}}, \qquad (65)$$

for all $v \in H^1(\kappa)$. Now, let $\mu_H \in \Lambda_H$ and consider the finite element space

$$V_H(K) := \{v_H \in \mathcal{C}^0(K) : v_H|_{\kappa} \in \mathbb{P}_{\ell+d}(\kappa) \text{ for all } \kappa \in \Xi_H(K)\}.$$

Observe that the product function $g_F := (\hat{p} \circ S_F^{-1}) \mu_H|_{\partial K}|_F$ belongs to $\mathbb{P}_{\ell+d}(F)$ and $g_F|_{\partial F} = 0$, for all $F \in \mathcal{E}_H^H$ such that $F \subset \partial K$, and let $w_F \in V_H(K)$ be the extension by zero of $g_F$ to $K$. Note that $w_F$ has support in a single $\kappa_F \in \Xi_H(K)$ from the definition of the mesh $\Xi_H(K)$. This implies that the functions
in {$w_F \in V_H(K) : F \in \mathcal{E}_H^H$, and, $F \subset \partial K$} are $H^1(K)$-orthogonal. In addition, from (65), [35, Proposition 3.37], and the quasi-uniformity of {$\mathcal{E}_H^H$}$_{H>0}$, we arrive at

$$
\begin{align*}
\|w_F\|_{V(K)} &= \left( |w_F|^2_{1,\kappa_F} + \frac{1}{d_{\Omega}^2} \|w_F\|_{0,\kappa_F}^2 \right)^{\frac{1}{2}} \le \hat{c} H^{\frac{d}{2}-1} \|\hat{w}_F\|_{1,\hat{\kappa}} \le \hat{c} H^{\frac{d}{2}-1} \|\hat{\mu}_H\|_{0,\hat{F}} \\
&= \hat{c} H^{\frac{d}{2}-1} \frac{|\hat{F}|^{\frac{1}{2}}}{|F|^{\frac{1}{2}}} \|\mu_H\|_{0,F} \le \frac{\hat{c}}{H^{\frac{1}{2}}} \|\mu_H\|_{0,F}.
\end{align*}
$$

Thus, by defining $w_K := \sum_{F \in \mathcal{E}_H^H(K)} w_F$ we obtain

$$
\|w_K\|_{V(K)}^2 = \sum_{F \in \mathcal{E}_H^H(K)} \|w_F\|_{V(K)}^2 \le \sum_{F \in \mathcal{E}_H^H(K)} \frac{\hat{c}^2}{H} \|\mu_H\|_{0,F}^2 \le \frac{\hat{c}^2}{H} \|\mu_H\|_{0,\partial K}^2. \quad (66)
$$

For $F \in \mathcal{E}_H^H(K)$, we use standard scaling arguments and (64) to get

$$
\int_F \mu_H w_K = \int_F (\hat{p} \circ S_F^{-1}) \mu_H^2 \ge \hat{c} \|\mu_H\|_{0,F}^2. \quad (67)
$$

Let define now $w \in V$ such that $w|_K = w_K$ for all $K \in \mathcal{P}_H$. Then, from (67) and (66) we arrive at

$$
\begin{align}
\sup_{v \in V} \frac{\langle \mu_H, v \rangle_{\partial \mathcal{P}_H}}{\|v\|_V} &\ge \frac{\sum_{K \in \mathcal{P}_H} \sum_{F \in \mathcal{E}_H^H(K)} \int_F \mu_H w_K}{\left( \sum_{K \in \mathcal{P}_H} \|w_K\|_{V(K)}^2 \right)^{\frac{1}{2}}} \ge \hat{c} H^{\frac{1}{2}} \frac{\sum_{K \in \mathcal{P}_H} \sum_{F \in \mathcal{E}_H^H(K)} \|\mu_H\|_{0,F}^2}{\left( \sum_{K \in \mathcal{P}_H} \|\mu_H\|_{0,\partial K}^2 \right)^{\frac{1}{2}}} \nonumber \\
&\ge \hat{c} H^{\frac{1}{2}} \left( \sum_{E \in \mathcal{E}_H^H} \|\mu_H\|_{0,E}^2 \right)^{\frac{1}{2}}. \tag{68}
\end{align}
$$

On the other hand, from the Cauchy-Schwarz inequality and (58), it holds

$$
\begin{align}
\sup_{v \in V} \frac{\langle \mu_H, v \rangle_{\partial \mathcal{P}_H}}{\|v\|_V} &= \sup_{v \in V} \frac{\sum_{K \in \mathcal{P}_H} \int_{\partial K} \mu_H v}{\|v\|_V} \le \sup_{v \in V} \frac{\sum_{K \in \mathcal{P}_H} \|\mu_H\|_{0,\partial K} \|v\|_{0,\partial K}}{\|v\|_V} \nonumber \\
&\le 2 \left( \sum_{E \in \mathcal{E}_H^H} \|\mu_H\|_{0,E}^2 \right)^{\frac{1}{2}} \sup_{v \in V} \frac{\left( \sum_{K \in \mathcal{P}_H} \|v\|_{0,\partial K}^2 \right)^{\frac{1}{2}}}{\|v\|_V} \nonumber \\
&\le \frac{c_t}{H^{\frac{1}{2}}} \left( \sum_{E \in \mathcal{E}_H^H} \|\mu_H\|_{0,E}^2 \right)^{\frac{1}{2}}, \tag{69}
\end{align}
$$

where $c_t$ is a positive constant independent on $\mathcal{H}$ or $H$. Finally, we use (68), (69) and (61) to obtain (62). $\square$

4.3. Conditioning of local problems

Let $\mathbf{a}^K \in \mathbb{R}^{\dim V_h(K) \times \dim V_h(K)}$ and $\tilde{\mathbf{a}}^K \in \mathbb{R}^{(\dim V_h(K)-1) \times (\dim V_h(K)-1)}$, with $K \in \mathcal{P}_H^0$ and $K \in \mathcal{P}_H \setminus \mathcal{P}_H^0$, respectively, the local matrices induced by left-hand side of problems (39)–(40) acting on the basis functions defined in Section 4.1. A detailed definition of such matrices is deferred until Section 5.1. This section aims to propose upper bounds for these local matrices. The proof is postponed to the end of this section.
**Theorem 5.** Assume that (A1)-(A3) hold. Then, there exist positive constants $c_1$ and $c_2$, independent of meshes and physical parameters, such that

$$ \mathrm{cond}(\mathbf{a}^K) \le c_1 \left( \frac{\omega + \|a\|}{\omega} \right)^2 \frac{\|a\|}{\varepsilon_{\mathrm{min}}} \frac{1}{h^2} \quad \text{for all } K \in \mathcal{P}_{\mathcal{H}}^0, \qquad (70) $$

$$ \mathrm{cond}(\tilde{\mathbf{a}}^K) \le c_2 \frac{\|a\|}{\varepsilon_{\mathrm{min}}} \frac{1}{h^2} \quad \text{for all } K \in \mathcal{P}_{\mathcal{H}} \setminus \mathcal{P}_{\mathcal{H}}^0. \qquad (71) $$

In addition, if $K \in \mathcal{P}_{\mathcal{H}}^0$, with $\nu|_K \ge \nu_0 > 0$, there exists a positive constant $c_3$, independent of meshes and physical parameters, such that

$$ \mathrm{cond}(\mathbf{a}^K) \le c_3 \frac{\|a\|}{h^2} \left\{ \begin{array}{ll} \nu_0 d_{\Omega}^2 & \text{if } \varepsilon_{\mathrm{min}} \ge \nu_0 d_{\Omega}^2, \\ \frac{c_I^2 \varepsilon_{\mathrm{min}} + \nu_0 h^2}{c_I^2 + \frac{h^2}{d_{\Omega}^2}} & \text{if } \varepsilon_{\mathrm{min}} < \nu_0 d_{\Omega}^2, \end{array} \right. \qquad (72) $$

where $c_I > 0$ is given in (57).

Before proving such a result, and since conditioning is closely related to the (inverse of) inf-sup constant $\gamma$, we first estimate a lower bound for it.

**Lemma 6.** Under assumptions (A1) and (A3), the constant $\gamma$ given in (41) satisfies the following lower bound

$$ \gamma \ge \left( \frac{\omega}{\omega + \|a\|} \right)^2 \frac{\varepsilon_{\mathrm{min}}}{c_W^2}, \qquad (73) $$

with $c_W$ given in (56) and $\omega$ in (A3).

*Proof.* Let $0 \ne w_h$ be an arbitrary function in $\tilde{W}_h$ and consider the decomposition $w_h = w_0 + \tilde{w}_h \in V_0 \oplus \tilde{V}_h$. Note that $w_0|_K = 0$ for all $K \in \mathcal{P}_{\mathcal{H}} \setminus \mathcal{P}_{\mathcal{H}}^0$ from (36). Using the existence of $\omega$ and $\xi_K^*$ in (55), we define $\bar{v}_h \in \tilde{W}_h$ as

$$ \bar{v}_h|_K := \left\{ \begin{array}{ll} 0, & K \in \mathcal{P}_{\mathcal{H}} \setminus \mathcal{P}_{\mathcal{H}}^0, \\ \frac{\|w_0\|_{V(K)}}{\|\xi_K^*\|_{V(K)}} \xi_K^*, & K \in \mathcal{P}_{\mathcal{H}}^0, \end{array} \right. \qquad (74) $$

and then $\|\bar{v}_h\|_V = \|w_0\|_V$ and $a(w_0, \bar{v}_h) \ge \omega \|w_0\|_V^2$. Next, we consider $v_h^* \in \tilde{W}_h$

$$ v_h^*|_K := \left\{ \begin{array}{ll} \left(\omega + \frac{\|a\|^2}{\omega}\right) \frac{\omega c_W^2}{\varepsilon_{\mathrm{min}}} w_h, & K \in \mathcal{P}_{\mathcal{H}} \setminus \mathcal{P}_{\mathcal{H}}^0, \\ \left(\omega + \frac{\|a\|^2}{\omega}\right) \frac{c_W^2}{\varepsilon_{\mathrm{min}}} w_h + 2 \bar{v}_h, & K \in \mathcal{P}_{\mathcal{H}}^0. \end{array} \right. \qquad (75) $$

For $K \in \mathcal{P}_{\mathcal{H}}^0$, and using $\|w_h\|_{V(K)}^2 = \|\tilde{w}_h\|_{V(K)}^2 + \|w_0\|_{V(K)}^2$, it holds

$$ \begin{aligned} a_K(w_h, v_h^*) &= \left(\omega + \frac{\|a\|^2}{\omega}\right) \frac{c_W^2}{\varepsilon_{\mathrm{min}}} a_K(w_h, w_h) + 2a_K(w_0, \bar{v}_h) + 2a_K(\tilde{w}_h, \bar{v}_h) \\ &\ge \left(\omega + \frac{\|a\|^2}{\omega}\right) \|\tilde{w}_h\|_{V(K)}^2 + 2\omega \|w_0\|_{V(K)}^2 - 2\|a\| \|\tilde{w}_h\|_{V(K)} \|w_0\|_{V(K)} \\ &= \omega \|w_h\|_{V(K)}^2 + \left(\frac{\|a\|}{\omega^{\frac{1}{2}}} \|\tilde{w}_h\|_{V(K)} - \omega^{\frac{1}{2}} \|w_0\|_{V(K)}\right)^2 \ge \omega \|w_h\|_{V(K)}^2, \end{aligned} $$

and then $a(w_h, v_h^*) \ge \omega \|w_h\|_V^2$. On the other hand, from the definition of $v_h^*$ in (75) we get

$$ \|v_h^*\|_V \le \frac{c_W^2}{\varepsilon_{\mathrm{min}}} \left( \frac{\omega^2 + 2 \frac{\varepsilon_{\mathrm{min}}}{c_W^2} \omega + \|a\|^2}{\omega} \right) \|w_h\|_V \le \frac{c_W^2}{\varepsilon_{\mathrm{min}}} \frac{(\omega + \|a\|)^2}{\omega} \|w_h\|_V, $$
then it follows

$$ \frac{a(w_h, v_h^*)}{\|v_h^*\|_V} \ge \left( \frac{\omega}{\omega + \|a\|} \right)^2 \frac{\varepsilon_{\min}}{c_W^2} \|w_h\|_V, \qquad (76) $$

and the result follows.

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

Furthermore, we can particularize (73) for the three different asymptotic regimes present in the local problems defined in $K \in \mathcal{P}_H$, to obtain more precise estimates. Specifically:

* Assume $K \in \mathcal{P}_H^0$ for an arbitrary $\nu$. In this case (76) implies that for all $v_h \in V_h(K)$ there exists $v_h^*$ such that

$$ a_K(v_h, v_h^*) \ge \left( \frac{\omega}{\omega + \|a\|} \right)^2 \frac{\varepsilon_{\min}}{c_W^2} \|v_h\|_{V(K)} \|v_h^*\|_{V(K)}. \qquad (77) $$

* Assume $K \in \mathcal{P}_H \setminus \mathcal{P}_H^0$. Then,

$$ a_K(\tilde{v}_h, \tilde{v}_h) \ge \frac{\varepsilon_{\min}}{c_W^2} \|\tilde{v}_h\|_{V(K)}^2 \quad \text{for all } \tilde{v}_h \in V_h(K) \cap L_0^2(K), \qquad (78) $$

where $c_W$ is given in (56).

* Assume $K \in \mathcal{P}_H^0$ with $\nu$ such that $\nu|_K \ge \nu_0 > 0$. In this case, we have

$$ a_K(v_h, v_h) \ge \varepsilon_{\min} \|\nabla v_h\|_{0,K}^2 + \nu_0 \|v_h\|_{0,K}^2 \quad \text{for all } v_h \in V_h(K). $$

In the asymptotic regime $\varepsilon_{\min} \ge \nu_0 d_{\Omega}^2$ we arrive at

$$ a_K(v_h, v_h) \ge \nu_0 d_{\Omega}^2 \|v_h\|_{V(K)}^2 \quad \text{for all } v_h \in V_h(K). \qquad (79) $$

Conversely, by assuming $\varepsilon_{\min} < \nu_0 d_{\Omega}^2$ and using (57), we get

$$ a_K(v_h, v_h) \ge \frac{c_I^2 \varepsilon_{\min} + \nu_0 h^2}{c_I^2 + \frac{h^2}{d_{\Omega}^2}} \|v_h\|_{V(K)}^2, \quad \text{for all } v_h \in V_h(K), \qquad (80) $$

where $c_I > 0$ is independent of $h$ or any data parameter.

Now, we are ready to prove Theorem 5.

*Proof of Theorem 5.* Applying (51) to the matrices induced by the left side of (39) and (40), and using the norm equivalence in (59) and the estimates in (77)–(80), the result follows. □

### 4.4. Conditioning of the global problem

This section is dedicated to prove upper bound for the condition number of matrices associated with the left-hand side of the MHM method (42). For that, we define the mapping $\mathcal{A} : \Lambda_H \to \Lambda'_H$, $\mathcal{B} : \Lambda_H \to W'_0$ and $\mathcal{B}^T : W_0 \to \Lambda'_H$ given by

$$ \langle \mathcal{A} \rho_H, \mu_H \rangle_{\Lambda'_H, \Lambda_H} := - \langle \mu_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} \qquad (81) $$

$$ \langle \mathcal{B} \mu_H, w_0 \rangle_{W'_0, W_0} = \langle \mathcal{B}^T w_0, \mu_H \rangle_{\Lambda'_H, \Lambda_H} := - \langle \mu_H, v_0 \rangle_{\partial \mathcal{P}_H}. \qquad (81) $$

We also define $A \in \mathbb{R}^{\dim(\Lambda_H) \times \dim(\Lambda_H)}$ and $B \in \mathbb{R}^{\dim(W_0) \times \dim(\Lambda_H)}$, as the matrices associated to $\mathcal{A}$ and $\mathcal{B}$, respectively. Specifically, we obtain such matrices by applying the operators to the canonical basis of $W_0$ and the Lagrange polynomials defined on $\mathcal{E}_H^H$ as the basis for $\Lambda_H$, for more details, see Section 5.1.
**Theorem 7.** Assume that (A1)-(A3) hold. Then there exist positive constants $C_i$, $i = 1, \dots, 4$, independent of physical parameters, and mesh parameters $\mathcal{H}$ and $\mathcal{H}$, but which may depend on $h$, such that

$$ \mathrm{cond}(A) \leq \frac{C_1}{\gamma \theta} \frac{1}{\mathcal{H} H} \quad \text{if } W_0 = \{0\}, \qquad (82) $$

$$ \mathrm{cond}\left(\begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix}\right) \leq \frac{C_2}{\gamma \varrho} \max\left\{C_3 \frac{1}{\mathcal{H} H}, C_4 \left(\frac{\mathcal{H}}{H}\right)^d\right\} \quad \text{if } W_0 \neq \{0\}, \qquad (83) $$

where $\gamma$ is bounded in Lemma 6, and

$$ \theta := \frac{\varepsilon_{\min}}{\|a\| (\theta_1 \|a\| + \theta_2 \varepsilon_{\min})} \quad \text{and} \quad \varrho \geq \frac{\beta \theta \gamma^2}{d_{\Omega} ((1 + \gamma + \theta)^2 \gamma^2 \beta^2 + (1 + \beta)^2 \theta^2)^{\frac{1}{2}}}, $$

where $\theta_1$, $\theta_2$ and $\beta$ are positive constant independent of mesh and physical parameters.

The proof of Theorem 7 is postponed until the end of this section. Before, we first estimate lower bounds for the inf-sup constants associated with the matrices $A$ and $B$.

**Lemma 8.** Assume (A1) and (A2) hold. There exist positive constants $\theta_1$ and $\theta_2$, independent of mesh or physical parameters, such that if we define $\theta := \frac{\varepsilon_{\min}}{\|a\|(\theta_1 \|a\| + \theta_2 \varepsilon_{\min})}$, then

$$ \inf_{\rho_H \in \mathcal{N}_H} \sup_{\mu_H \in \mathcal{N}_H} \frac{\langle \mathcal{A} \rho_H, \mu_H \rangle_{\Lambda'_H, \Lambda_H}}{\|\rho_H\|_{\Lambda} \|\mu_H\|_{\Lambda}} \geq \theta, \qquad (84) $$

where

$$ \mathcal{N}_H := \{ \mu_H \in \Lambda_H : \langle \mathcal{B} \mu_H, w_0 \rangle_{W'_0, W_0} = 0 \text{ for all } w_0 \in W_0 \}. \qquad (85) $$

Moreover, if $W_0 \neq \{0\}$, there exists a positive constant $\beta$, independent of any mesh or physical parameters, such that

$$ \inf_{w_0 \in W_0} \sup_{\mu_H \in \Lambda_H} \frac{\langle \mathcal{B} \mu_H, w_0 \rangle_{W'_0, W_0}}{\|w_0\|_V \|\mu_H\|_{\Lambda}} \geq \beta. \qquad (86) $$

*Proof.* Let $\rho_H \in \mathcal{N}_H$, and consider $T_h \rho_H \in \tilde{W}_h$ with its orthogonal decomposition $T_h \rho_H = v_0 + \tilde{v}_h \in V_0 \oplus \tilde{V}_h$. Consider the global lowest order Raviart-Thomas space defined on the virtual mesh $\Xi_H$, i.e.,

$$ X_H := \{ \tau_H \in H(\mathrm{div}; \Omega) : \tau_H|_{\kappa} \in RT_0(\kappa) \text{ for all } \kappa \in \Xi_H \}. \qquad (87) $$

Given $v_0 \in W_0$, there exists $\tilde{\tau}_H \in X_H$ such that $\nabla \cdot \tilde{\tau}_H = -v_0$ in $\Omega$ and a positive constant $\beta$ independent of mesh parameters

$$ \beta \|\tilde{\tau}_H\|_{\mathrm{div}} \leq d_{\Omega} \|v_0\|_{0,\Omega} \leq d_{\Omega}^2 \|v_0\|_V. \qquad (88) $$

Thus,

$$ \begin{aligned} (\tilde{\tau}_H, \nabla T_h \rho_H)_{\mathcal{P}_H} &\leq \|\tilde{\tau}_H\|_{0,\Omega} |T_h \rho_H|_{1,\mathcal{P}_H} \leq d_{\Omega}^2 \frac{1}{\beta} \|v_0\|_V |T_h \rho_H|_{1,\mathcal{P}_H} \\ &= d_{\Omega}^2 \frac{1}{\beta} \|v_0\|_V |\tilde{v}|_{1,\mathcal{P}_H} \leq d_{\Omega}^2 \frac{1}{\beta} \|v_0\|_V \|\tilde{v}\|_V. \end{aligned} \qquad (89) $$

Then, take $\tilde{\mu}_H \in \Lambda_H$ be such that $\tilde{\mu}_H := \tilde{\tau}_H \cdot n^K|_F$ on every $F \in \mathcal{E}_H^H$. Using integration by parts, the definition of $\tilde{\tau}_H$, and (89), we arrive at

$$ \begin{aligned} \langle \tilde{\mu}_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} &= (\tilde{\tau}_H, \nabla T_h \rho_H)_{\mathcal{P}_H} + (\nabla \cdot \tilde{\tau}_H, T_h \rho_H)_{\mathcal{P}_H} \\ &\leq d_{\Omega}^2 \frac{1}{\beta} \|v_0\|_V \|\tilde{v}\|_V - d_{\Omega}^2 \|v_0\|_V^2. \end{aligned} $$
Let $w_0 \in W_0$, and observe that $w_0|_K = 0$ on every $K \in \mathcal{P}_H^0$ (see (33)). On the other hand, as a consequence of $(v_0 + \tilde{v}_h)|_K = T_h \rho_H|_K \in V_h(K) \cap L_0^2(K)$ we have that $v_0|_K = 0$ on every $K \in \mathcal{P}_H \setminus \mathcal{P}_H^0$, then $(w_0, v_0)_{\mathcal{P}_H} = 0$. As a result,

$$ \langle \tilde{\mu}_H, w_0 \rangle_{\partial \mathcal{P}_H} = \sum_{K \in \mathcal{P}_H} \langle \tilde{\tau}_H \cdot n^K, w_0 \rangle_{\partial K} = (\nabla \cdot \tilde{\tau}_H, w_0)_{\mathcal{P}_H} = -(w_0, v_0)_{\mathcal{P}_H} = 0, $$

and, thus, $\tilde{\mu}_H \in \mathcal{N}_H$. We define now $\rho_H^* := \frac{c_W^2}{\varepsilon_{\min}} \left(1 + \frac{1}{\beta^2}\right) \rho_H + 2 \frac{1}{d_\Omega^2} \tilde{\mu}_H \in \mathcal{N}_H$, then

$$ \begin{aligned} - \langle \rho_H^*, T_h \rho_H \rangle_{\partial \mathcal{P}_H} &= - \frac{c_W^2}{\varepsilon_{\min}} \left(1 + \frac{1}{\beta^2}\right) \langle \rho_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} - 2 \frac{1}{d_\Omega^2} \langle \tilde{\mu}_H, T_h \rho_H \rangle_{\partial \mathcal{P}_H} \\ &\geq \left(1 + \frac{1}{\beta^2}\right) \|\tilde{v}_h\|_V^2 - \frac{2}{\beta} \|v_0\|_V \|\tilde{v}_h\|_V + 2 \|v_0\|_V^2 \\ &= \|T_h \rho_H\|_V^2 + \left(\frac{1}{\beta} \|\tilde{v}_h\|_V - \|v_0\|_V\right)^2 \geq \|T_h \rho_H\|_V^2. \end{aligned} \quad (90) $$

Next, using (30) and (46) we arrive at

$$ \begin{aligned} c_F^{-1} \|\rho_H\|_{\Lambda} &\le \sup_{v_h \in V_h} \frac{\langle \rho_H, v_h \rangle_{\partial \mathcal{P}_H}}{\|v_h\|_V} \\ &\le \sup_{\tilde{v}_h \in \tilde{W}_h} \frac{\langle \rho_H, \tilde{v}_h \rangle_{\partial \mathcal{P}_H}}{\|\tilde{v}_h\|_V} \\ &= \sup_{\tilde{v}_h \in \tilde{W}_h} \frac{a(T_h \rho_H, \tilde{v}_h)}{\|\tilde{v}_h\|_V} \\ &\le \|a\| \|T_h \rho_H\|_V, \end{aligned} $$

thus $\|\rho_H\|_{\Lambda} \le c_F \|a\| \|T_h \rho_H\|_V$. Using (88) and the definition of $\rho_H^*$, we obtain

$$ \begin{aligned} \|\rho_H^*\|_{\Lambda} &\le \frac{c_W^2}{\varepsilon_{\min}} \left(1 + \frac{1}{\beta^2}\right) \|\rho_H\|_{\Lambda} + 2 \frac{1}{d_{\Omega}^2} \|\tilde{\mu}_H\|_{\Lambda} \\ &\le \frac{c_W^2}{\varepsilon_{\min}} \left(1 + \frac{1}{\beta^2}\right) c_F \|a\| \|T_h \rho_H\|_V + 2 \frac{1}{\beta} \|v_0\|_V \\ &\le \left(\frac{c_W^2}{\varepsilon_{\min}} \left(1 + \frac{1}{\beta^2}\right) c_F \|a\| + 2 \frac{1}{\beta}\right) \|T_h \rho_H\|_V, \end{aligned} $$

and then,

$$ \|\rho_H\|_{\Lambda} \|\rho_H^*\|_{\Lambda} \le \left( \theta_1 \frac{\|a\|}{\varepsilon_{\min}} + \theta_2 \right) \|a\| \|T_h \rho_H\|_V^2, \quad (91) $$

where the constants $\theta_1$ and $\theta_2$ are independent of mesh and data parameters. The result (84) follows by replacing (90) into (91). Using $W_0$ as a subspace of $V_0$, the inf-sup condition (86) follows from [7, eq. (28)]. □

The next result follows the analysis proposed in [36, Theorem 3] by combining (86) and (84).

**Corollary 9.** Assume (A1) and (A2) hold and $W_0 \ne \{0\}$. Then,

$$ \inf_{(\rho_H, w_0) \in \Lambda_H \times W_0} \sup_{(\mu_H, v_0) \in \Lambda_H \times W_0} \frac{\langle A \rho_H, \mu_H \rangle_{\Lambda'_H, \Lambda_H} + \langle B^T w_0, \mu_H \rangle_{\Lambda'_H, \Lambda_H} + \langle B \mu_H, w_0 \rangle_{W'_0, W_0}}{\left( \|\rho_H\|_{\Lambda}^2 + \|w_0\|_V^2 \right)^{\frac{1}{2}} \left( \|\mu_H\|_{\Lambda}^2 + \|v_0\|_V^2 \right)^{\frac{1}{2}}} \ge \varrho $$

where $\varrho$ satisfies

$$ \varrho \ge \frac{\beta \theta \gamma^2}{d_{\Omega} ((1 + \gamma + \theta)^2 \gamma^2 \beta^2 + (1 + \beta)^2 \theta^2)^{\frac{1}{2}}}, \quad (92) $$

with $\gamma$ bounded in Lemma 6, and $\theta$ and $\beta$ in Lemma 8.
Finally, we can prove Theorem 7.

**Proof of Theorem 7.** Regarding the condition number of the matrix that defines the global problem (81), it is first necessary to have an estimate for the norm of the bilinear form that defines (42). Therefore, using (81) and $||T_h|| \le \frac{1}{\gamma}$ we note that

$$ \sup_{\substack{\mu_H \in \Lambda_H \\ \rho_H \in \Lambda_H}} \frac{\langle \mathcal{A} \rho_H, \mu_H \rangle_{\Lambda'_H, \Lambda_H}}{\|\rho_H\|_{\Lambda} \|\mu_H\|_{\Lambda}} \le \sup_{\rho_H \in \Lambda_H} \frac{\|T_h \rho_H\|_V}{\|\rho_H\|_{\Lambda}} \le \frac{1}{\gamma}. \quad (93) $$

For the case $W_0 \ne \{0\}$, we get

$$ \sup_{\substack{(\mu_H, v_0) \in \Lambda_H \times W_0 \\ (\rho_H, w_0) \in \Lambda_H \times W_0}} \frac{\langle \mathcal{A} \rho_H, \mu_H \rangle_{\Lambda'_H, \Lambda_H} + \langle \mathcal{B}^T w_0, \mu_H \rangle_{\Lambda'_H, \Lambda_H} + \langle \mathcal{B} \mu_H, w_0 \rangle_{W'_0, W_0}}{(\|\rho_H\|_{\Lambda}^2 + \|w_0\|_V^2)^{\frac{1}{2}} (\|\mu_H\|_{\Lambda}^2 + \|v_0\|_V^2)^{\frac{1}{2}}} \le 1 + \frac{1}{\gamma}. \quad (94) $$

Thus, the bilinear form in the left-hand side of (42) is bounded by $1 + \frac{1}{\gamma}$.

Therefore, we prove (82) using (51) for the matrix A associated to A in (81), along with the equivalence given in (62), Lemma 8 and estimate (93). Following the same idea, we prove (83) using (51) for each of the block diagonal matrices related to (81) along with the equivalence (60) and Corollary 9, and estimate (94). □

**Remark 8.** As was highlighted in Remark 4, if $\nu \ge \nu_0 > 0$ in $\Omega$, then (42) simplifies to (44). In addition from (79) and (80) we arrive at the coercivity $a(v_h, v_h) \ge \gamma ||v_h||_V^2$ for all $v_h \in V_h$, with

$$ \gamma := \begin{cases} \nu_0 d_{\Omega}^2 & \text{if } \varepsilon_{\min} \ge \nu_0 d_{\Omega}^2 \\ \frac{c_I^2 \varepsilon_{\min} + \nu_0 h^2}{c_I^2 + \frac{h^2}{d_{\Omega}^2}} & \text{if } \varepsilon_{\min} < \nu_0 d_{\Omega}^2 \end{cases}. \quad (95) $$

Then, we also obtain the coercivity

$$ - \langle \mu_H, T_h \mu_H \rangle_{\partial \mathcal{P}_H} = a(T_h \mu_H, T_h \mu_H) \ge \frac{\gamma}{c_F \|a\|} \|\mu_H\|_{\Lambda}^2, \text{ for all } \mu_H \in \Lambda_H. $$

Therefore, $\rho = \theta \ge \frac{\gamma}{c_F \|a\|}$. This provides a lower bound for $\theta$ (and $\rho$) which is robust with respect to small values of the physical parameter $\varepsilon_{\min}$ and it does not depend on $\omega$.

# 5. Computational Aspects

## 5.1. The two-level MHM Algorithm

We present the matrix versions of the two-level MHM method (39)–(40) and (42) and the parallel algorithm underlying it. Let K be an element of the partition $\mathcal{P}_H$, consider the basis for $V_h(K)$ defined in Section 4.1 and set $N_K := \dim V_h(K)$. We define the local matrix $\mathbf{a}^K \in \mathbb{R}^{N_K \times N_K}$ with entries

$$ \mathbf{a}_{ij}^K := a_K(\phi_j, \phi_i), \text{ for all } i, j \in \{1, \dots, N_K\}. \quad (96) $$

Also, let $\tilde{\mathbf{a}}^K \in \mathbb{R}^{(N_K-1) \times (N_K-1)}$ be the reduced matrix, where $\tilde{\mathbf{a}}_{ij}^K := \mathbf{a}_{ij}^K$ for all $i, j \in \{1, \dots, N_K-1\}$, and $\mathbf{p}^K \in \mathbb{R}^{N_K \times (N_K-1)}$ a rectangular matrix with entries

$$ \mathbf{p}_{ij}^K := \delta_{ij} - \frac{1}{|K|} \int_K \phi_j \, d\mathbf{x}, \quad \text{for } (i, j) \in \{1, \dots, N_K\} \times \{1, \dots, N_K-1\}, $$
where $\delta_{ij}$ represents the Kronecker delta.

Let $\{\psi_1, \dots, \psi_{\dim \Lambda_H}\}$ be the basis for $\Lambda_H$ defined in Section 4.1. We observe that the restrictions $\{\psi_1^K, \dots, \psi_{M_K}^K\}$, for each $K \in \mathcal{P}_H$, naturally defines a local-global mapping

$$i_K: \{1, \dots, M_K\} \mapsto \{1, \dots, \dim \Lambda_H\}, \qquad (97)$$

that identify the global degrees of freedom with the local numbering. As such, we have $\psi_j^K = \psi_{i_K(j)}$ for $j \in \{1, \dots, M_K\}$. Then, we express the right-hand side of local problems (39)-(40) by defining $\mathfrak{b}^K \in \mathbb{R}^{N_K \times (M_K+1)}$ through

$$\begin{aligned} \mathfrak{b}_{ij}^K &:= -\langle \psi_j^K, \phi_i \rangle_{\partial K}, && \text{for all } (i,j) \in \{1, \dots, N_K\} \times \{1, \dots, M_K\}, \\ \mathfrak{b}_{i,M_K+1}^K &:= (f, \phi_i)_K, && \text{for all } i \in \{1, \dots, N_K\}. \end{aligned}$$

We also consider a reduced matrix $\tilde{\mathfrak{b}}^K \in \mathbb{R}^{(N_K-1) \times (M_K+1)}$ which is obtained by excluding the last row of the matrix $\mathfrak{b}^K$. Owing to the previous constructions, we proceed to compute the basis functions matrix $\mathfrak{c}^K$ as follows:

$$\left\{ \begin{array}{ll} \text{Solve: } \mathfrak{a}^K \mathfrak{c}^K = \mathfrak{b}^K & \text{for } K \in \mathcal{P}_H^0 \\ \text{Solve: } \tilde{\mathfrak{a}}^K \tilde{\mathfrak{c}}^K = \tilde{\mathfrak{b}}^K & \text{and } \mathfrak{c}^K = \mathfrak{p}^K \tilde{\mathfrak{c}}^K \text{ for } K \in \mathcal{P}_H \setminus \mathcal{P}_H^0. \end{array} \right. \qquad (98)$$

Thus, solution of, (39) for $j = 1, \dots, M_K$, and (40) are locally computed as

$$T_h \psi_j^K|_K = \sum_{i=1}^{N_K} \mathfrak{c}_{i,j}^K \phi_i, \quad \text{and} \quad \hat{T}_h f|_K = \sum_{i=1}^{N_K} \mathfrak{c}_{M_K+1,i}^K \phi_i. \qquad (99)$$

It remains to write the linear system associated with (42). First, we construct a matrix $A \in \mathbb{R}^{\dim \Lambda_H \times \dim \Lambda_H}$ along with a vector $\mathbf{f} \in \mathbb{R}^{\dim \Lambda_H}$, with entries $A_{i,j} := -\langle \psi_j, T_h \psi_i \rangle_{\partial \mathcal{P}_H}$ and $\mathbf{f}_j := \langle \psi_j, \hat{T}_h f \rangle_{\partial \mathcal{P}_H}$ for $i, j \in \{1, \dots, \dim \Lambda_H\}$. The procedure to treat local problems detailed in (98) allows to construct A and f locally, through the product of matrices $\mathbf{q} := (\mathfrak{c}^K)^T \mathfrak{b}^K$ on each $K \in \mathcal{P}_H$. Then, by using the mapping (97), we assemble A and f as follows

$$\mathfrak{q}_{i,j}^K \xrightarrow{\text{contributes to}} A_{i_K(j),i_K(i)}, \quad \text{and} \quad -\mathfrak{q}_{M_K+1,j}^K \xrightarrow{\text{contributes to}} \mathbf{f}_{i_K(j)},$$

for $i, j \in \{1, \dots, N_K\}$. If $W_0$ is non-trivial let $\{w_1, \dots, w_{\dim W_0}\}$ be its canonical basis. In this case, we need to construct the matrix $B \in \mathbb{R}^{\dim W_0 \times \dim \Lambda_H}$ and the vector $\mathbf{g} \in \mathbb{R}^{\dim W_0}$ with entries $B_{i,j} := -\langle w_j, T_h \psi_i \rangle_{\partial \mathcal{P}_H}$ and $\mathbf{g}_j := -(f, w_j)_\Omega$, for $i \in \{1, \dots, \dim \Lambda_H\}$ and $j \in \{1, \dots, \dim W_0\}$, respectively. To localize those computations, we consider $K \in \mathcal{P}_H \setminus \mathcal{P}_H^0$ and set $n_K$ as the unique index in $\{1, \dots, \dim W_0\}$ such that $w_{n_K} = 1_K$. Furthermore, we set $\mathbf{n}_K = [1, \dots, 1] \in \mathbb{R}^{N_K}$ and $\mathbf{r}^K = \mathbf{n}_K \mathfrak{b}^K$, thus for $j \in \{1, \dots, N_K\}$ we have

$$\mathbf{r}_j^K \xrightarrow{\text{contributes to}} B_{n_K, i_K(j)}, \quad \text{and,} \quad -\mathbf{r}_{M_K+1}^K \xrightarrow{\text{contributes to}} \mathbf{g}_{n_K}.$$

Finally, the matrix form (42) writes

$$\begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} \lambda_H \\ \mathbf{u}_0 \end{bmatrix} = \begin{bmatrix} \mathbf{f} \\ \mathbf{g} \end{bmatrix}. \qquad (100)$$

Therefore, $u_0^H \in W_0$ is recovered by $u_0^H|_K = (\mathbf{u}_0)_{n_K}$ for all $K \in \mathcal{P}_H \setminus \mathcal{P}_H^0$, and the full discrete MHM's solution $u_{H,h}$ restricted to $K \in \mathcal{P}_H$ writes

$$u_{H,h}|_K = u_0^H|_K + \sum_{j=1}^{M_K} (\lambda_H)_{i_K(j)} T_h \psi_j^K|_K + \hat{T}_h f|_K. \qquad (101)$$

The process of building the linear system in (100) serves as the foundation for implementing the MHM method. A crucial point concerns the practical implementation of the $\mathcal{P}_H^0$ partition criterion. An element $K$ is categorized into $\mathcal{P}_H \setminus \mathcal{P}_H^0$ if $\|\mathbf{a}^K \mathbf{n}_K^T\|_2$ is below a predefined numerical tolerance, otherwise it belongs to $\mathcal{P}_H^0$. The entire process is described in Algorithm 1.
<caption>**Algorithm 1** Two-Level MHM Algorithm</caption>

for $K \in \mathcal{P}_H$ do (in parallel)
  $n_K = [1, \dots, 1] \in \mathbb{R}^{N_K}$
  if $(||a^K n_K^T||_2 > \text{tol})$ then
    solve $a^K c^K = b^K$
    set $u_0^H|_K = 0$
  else
    solve $\tilde{a}^K \tilde{c}^K = \tilde{b}^K$
    set $c^K = p^K \tilde{c}^K$
    $n_K b^K \xrightarrow{\text{contributes to}} B$ and $g$
  end if
  $(c^K)^T b^K \xrightarrow{\text{contributes to}} A$ and $f$
end for

solve $\begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} \lambda_H \\ u_0 \end{bmatrix} = \begin{bmatrix} f \\ g \end{bmatrix}$

for $K \in \mathcal{P}_H$ do (in parallel)
  process $u_{H,h}|_K$ from (99) and (101)
end for

## 5.2. Numerical experiments

We present two-dimensional and three-dimensional numerical experiments verifying the theoretical results for the two-level MHM method (42) using solutions of the RAD equation (1) with highly oscillatory behavior or boundary layers. They complement the extensive numerical validation presented in [25]. We distinguish two types of convergence in numerical tests based on Theorem 3:

* $\mathcal{H} \to 0$. This is the usual mesh convergence, called *mesh-based*;

* $H \to 0$ with fixed $\mathcal{H}$, called *space-based*.

We adopt second-level simplicial meshes to approximate the solution of local problems and select the diameter $h$ in order to avoid second-level pollution in the convergence according to the error estimates (23). Furthermore, we verify the theoretical upper bounds for local and global matrix condition numbers in terms of physical and mesh parameters.

### 5.2.1. An analytical oscillatory case

Consider the diffusion-advection problem defined in a unit square, with diffusive coefficient $\varepsilon = 10^{-1} \mathcal{I}$ and advective field $\alpha = (1, 0)$ if $d = 2$, and $\alpha = (1, 0, 0)$ if $d = 3$. We define the right-hand side such that the exact solution to the problem is given by

$$ 
\begin{align*}
u(x, y) &= \sin(2n\pi x) \sin(2m\pi y), \quad \text{with } n, m \in \mathbb{N}, \text{ if } d = 2, \\
u(x, y, z) &= \sin(2n\pi x) \sin(2m\pi y) \sin(2l\pi z), \quad \text{with } n, m, l \in \mathbb{N}, \text{ if } d = 3.
\end{align*}
$$

We assess the error estimates in the two-dimensional case using meshes composed by the following polygons: triangle, square, rhombus, hexagon, and L-shaped. The second-level local meshes correspond to the minimum triangulations allowed by the stability results in (A2) for the case $k = \ell + 2$. The convergence profiles for $m = 7$ and $n = 3$ are presented in Figure 2, where we can observe that the errors tend to zero as predicted by theory for the mesh-based case (see Theorem 3). For the space-based case, we found that the numerical solution is superconvergent with an extra $O(H^{1/2})$ rate. Such behavior was recently demonstrated in [14] for the Poisson model. We illustrate the isolines of the solution in some of the polygonal meshes in Figure 3.

For the 3D tests, we set $\ell = 1$ and $k = 4$, and the mesh-based strategy in the MHM method. In this case, we adopt $l = 1, m = 2, n = 3$ and tetrahedral meshes to compare our results with those obtained
<figure><img src="image_4.png" /><figcaption>Figure 2: Convergence rates for mesh-based (left) and space-based (right) approach.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Isolines of the solution using $\mathcal{H} = H$, $\ell = 1$ and $k = 3$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Isovalues of the solution at z = 0.25, using non-regular tetrahedral meshes. Galerkin approximation $u_H$ computed with quadratic Lagrange polynomials (left); MHM approximations with coarse (center) and fine (right) second-level refinements.</figcaption></figure>

by the classical Galerkin method using quadratic interpolations. In Figure 4, we present the isolines of the solutions. We notice that the MHM method achieves a high approximation accuracy with a coarser mesh at the first-level and a similar amount of degrees of freedom for the global problem compared to the Galerkin method. Refining the second-level, the recovery of the multiscale features is improved; however, the effect on the error in norm $||·||_V$ is negligible (see Figure 4c) reproducing (23).

### 5.2.2. An analytical boundary layer case

The domain $\Omega = (0, 1)^2$, and the physical coefficients are $\varepsilon = \varepsilon \mathcal{I}$, $\varepsilon \in \mathbb{R}^+$, $\alpha = (a, 0)$, $a \in \mathbb{R}$, and $\sigma = 0$. The source term $f = 1$, and the Dirichlet boundary conditions at $x = 0$ and $x = 1$, and the Neumann boundary conditions at $y = 0$ and $y = 1$ are chosen such that the exact solution is given by

$$u(x, y) := \frac{1}{a} \left( x - \frac{\sinh\left(\frac{a}{2\varepsilon}x\right)}{\sinh\left(\frac{a}{2\varepsilon}\right)} e^{\frac{a}{2\varepsilon}(x-1)} \right).$$

In Figure 5, we note the high accuracy of the MHM method for approximating boundary layers without spurious oscillations in polygonal meshes (hexagons). In Figure 6 we present convergence curves with respect

<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of elevation between the Galerkin method with $k = 2$ (left), and MHM with $\ell = 1, k = 3$ on hexagonal elements (center). Profiles of Galerkin (blue) and MHM (red) solutions and the exact solution (black dot) at $y = 0.4375$ (right). Here $a = 1$ and $\varepsilon = 10^{-2}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison between Galerkin and MHM method. We denote by $u_H$ the approximate solution given by the Galerkin method with quadratic Lagrange polynomials. For MHM, we set $\ell = 1$, $k = 3$ using mesh- and space-based strategies with hexagonal elements. Here $a = 1$ and $\epsilon = 10^{-2}$.</figcaption></figure>

to mesh parameters for the mesh- and space-based strategies of the MHM method, and compare them with those obtained from the Galerkin method. We also compare the MHM and Galerkin methods in terms of accuracy for different numbers of degrees of freedom (DoF). We observe that errors decrease drastically when using the MHM method instead of the Galerkin method, if we pay the price of computing multiscale basis functions in a completely parallel way.

### 5.3. Checking theoretical conditioning of matrices

This section is dedicated to validate the estimates presented in Theorem 5 and Theorem 7. Consider problem (1) with $\Omega = (0, 1)^2$, $f = 1$, $\epsilon = \epsilon I$, $\epsilon \in \mathbb{R}^+$, and $\sigma \in L^\infty(\Omega)$ defined by

$$ \sigma(x, y) = \begin{cases} \nu_1, & x \in (0.0, 0.5] \\ \nu_0, & x \in (0.5, 1.0) \end{cases}, $$

where $\nu_0, \nu_1 \in \mathbb{R}^+ \cup \{0\}$. To define $\alpha$, we partition $\Omega$ into 16 equal square regions $K_1, \dots, K_{16}$ and consider a vector function $\beta: \Omega \to \mathbb{R}^2$ such that $\beta|_{K_i} \in \mathbb{Q}_2(K_i)^2$ for all $i = 1, \dots, 16$, satisfying $\nabla \cdot \beta|_{K_i} = 0$ and $\beta \cdot n^{K_i} = 0$. Specifically, let $R_i: \hat{K} \to K_i$ be an affine transformation with $\hat{K} = [0, 1]^2$, define $\hat{\beta} \in \mathbb{Q}_2(\hat{K})$ as $\hat{\beta}(s, t) := ((s^2 - s)(1 - 2t), (1 - 2s)(t - t^2))$ for all $(s, t) \in \hat{K}$. Then, we set $\beta|_{K_i} := s_i \hat{\beta} \circ R_i^{-1}$, where $s_i$ is defined as 1 or -1 to preserve the continuity of $\beta$ on $\Omega$. Thus we set $\alpha(x, y) := \chi\beta(x, y) + (\delta(x, y), 0)$ for all $(x, y) \in \Omega$, with $\chi \in \mathbb{R}$ and $\delta: \Omega \to \mathbb{R}$ given by

$$ \delta(x, y) := \begin{cases} \omega, & (x, y) \in (0, 1) \times (0.0, 0.5] \\ (3 - 4y)(\omega - \omega_0) + \omega_0, & (x, y) \in (0, 1) \times (0.5, 0.75) \\ \omega_0, & (x, y) \in (0, 1) \times [0.75, 1.0) \end{cases}, $$

where $\omega, \omega_0 \in \mathbb{R}$. Then, we approximate the solution of (1) using the method (42) in a partition $\mathcal{P}_H$ composed of the squares $K_1, \dots, K_{16}$ as described in Figure 7. To cover a representative set of scenarios, we consider five cases as described in Table 1, for which we check the estimates provided by Theorem 5 and Theorem 7 in relation to $\epsilon, \omega$ and $H$. For each case of Table 1, we highlight the structure of the MHM method, that is, whether it has a mixed form such as (42) or primal given in (44). Next, the analysis of local problems is limited to the problems defined in the square highlighted in Figure 7. Our initial study, summarized in Figure 8, focuses on analyzing the dependence of condition numbers on $\epsilon$. We observed that cases 1, 3, and 5 corroborate the estimates (70) regarding local problems. As predicted, the condition number increases as $\epsilon^{-1}$ for small values of $\epsilon$, and for large values of $\epsilon$ the condition number increases as $\epsilon^2$ (see Figure 7 on left). According to (73) and (83), we anticipated that the global problem condition number for cases 1 and 2 would follow a rate $\epsilon^{-3}$ for small values of $\epsilon$. However, both numerical data sets
<figure><img src="image_9.png" /><figcaption>Figure 7: Illustration of the partition $\mathcal{P}_H$ with an example of $\alpha$ taking the values $\omega = 1$ and $\omega_0 = 0$. The red element is used to calculate the condition number of matrices associated with local problems.</figcaption></figure>

<table><caption>Table 1: Summary of physical cases with $\epsilon$ and $\omega$ to be set.</caption><thead><tr><th></th><th>χ</th><th>ω<sub>0</sub></th><th>ν<sub>0</sub></th><th>ν<sub>1</sub></th><th>Method's Type</th></tr></thead><tbody><tr><td>Case 1</td><td>7.5 × 10<sup>1</sup></td><td>0</td><td>0</td><td>0</td><td>Mixed</td></tr><tr><td>Case 2</td><td>1.0 × 10<sup>0</sup></td><td>0</td><td>0</td><td>1.0 × 10<sup>0</sup></td><td>Mixed</td></tr><tr><td>Case 3</td><td>1.0 × 10<sup>0</sup></td><td>5.0 × 10<sup>-1</sup></td><td>0</td><td>0</td><td>Primal</td></tr><tr><td>Case 4</td><td>1.0 × 10<sup>0</sup></td><td>0</td><td>1.0 × 10<sup>0</sup></td><td>1.0 × 10<sup>1</sup></td><td>Primal</td></tr><tr><td>Case 5</td><td>2.5 × 10<sup>1</sup></td><td>1.0 × 10<sup>0</sup></td><td>1.0 × 10<sup>0</sup></td><td>0</td><td>Primal</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure 8: Conditioning of matrices related to local (left) and global (right) problems with respect to $\epsilon$. Here, $\ell = 0$, $k = 2$, $h = H = \mathcal{H}$, and $\omega = 1$.</figcaption></figure>

consistently exhibit a $\epsilon^{-2}$ rate. As for larger values of $\epsilon$, both cases show an anticipated rate of $\epsilon^2$, which corroborates the theoretical results (see Figure 7 on right).

As proposed by Theorem 7, in the case $W_0 = \{0\}$, the condition number of the global problem (42) depends on the inverse of the product $\theta\gamma$. We observe that $\theta = O(\epsilon)$ and $\gamma = O(\epsilon)$ for small values of $\epsilon$. Then, we expect a rate of $\epsilon^{-2}$, that is accurately fulfilled by case 3. We notice that case 5 shows a better order $(\epsilon^{-1})$ than predicted by the theory. On the other hand, in case 4, using Remark 8, we can see that the condition number does not depend on $\epsilon$, which is corroborated by the numerical results (see Figure 8). Similarly, for $W_0 \neq \{0\}$ the condition number of the global problem (42) depends on the inverse of the product $\varrho\gamma$. We observe that $\varrho = O(\epsilon^2)$ for small values of $\epsilon$ thus, we expect that the condition number of the global problem has a rate of $\epsilon^{-3}$. The numerical results in Figure 8, for cases 1 and 2, show a better order of dependence of the condition number with respect to $\epsilon$, than the predicted by the theory.

The theoretical behavior of the condition number for local problems relative to small values of $\omega$ follows from the estimate (70). In cases 1, 3, and 5, a rate of $\omega^{-2}$ is predicted, while cases 2 and 4 are expected to exhibit stable behavior as predicted by (72). Notably, these theoretical estimates closely match the numerical results shown in Figure 9 (left). For large values of $\omega$, both theoretical estimates and numerical results show a linear increase with respect to $\omega$.

Figure 9 right shows that the condition number associated with the global problem increases as $\omega^2$ for large values of $\omega$ in all five cases. This behavior was anticipated by (82) for cases 3 and 5, but for cases 1 and 2, the theoretical prediction of (83) suggested a rate of $\omega^6$. This may indicate that the (83) estimate is pessimistic and can be improved.

Similarly, using (82) joint to (73) we can see that, cases 3 and 5 define global problems with a condition number following a rate of $\omega^{-2}$, for $\omega$ small. That behavior is corroborated by the numerical results in Figure 9 (right). On the other side, (83) with the definition of $\varrho$ in Theorem 7, predicts that the condition number of the global problem increases with a rate of $\omega^{-6}$, for small values of $\omega$. However, the numerical results in Figure 9 show a rate of $\omega^{-2}$ for cases 1 and 2. Notably, in case 4, the numerical experiments show that the condition number is robust with respect to small values of $\omega$, as predicted by (82) using Remark 8.

Estimates for the condition number of the global system, given in (82) and (83), may depend on the small length-scale parameter $h$ through the constant $C_1$. In Figure 10 (left), we numerically study the dependence of $\text{cond}(\mathcal{A})$ in terms of $h$, and one finds that it remains stable over a large range of $h$, indicating that the constant $C_1$ may depend only slightly on $h$. This is particularly attractive when using the MHM method in
<figure><img src="image_11.png" /><figcaption>Figure 9: Conditioning of matrices related to local (left) and global (right) problems with respect to $\omega$. Here, $\ell = 0, k = 2, h = H = \mathcal{H}$ and $\epsilon = 1$.</figcaption></figure>

highly multiscale problems where the second-level mesh must be refined.

Also in Figure 10 (right), we validate (82) and (83) in terms of $H$ for a fixed $\mathcal{H}$. We observe that the numerical experiments for the cases 3, 4 and 5 exactly recover the theoretical estimate given in (82). Interestingly, we find that the estimate behaves like $H^{-1}$ which is better than the theoretical estimate of order $H^{-2}$ given in (83). Additional studies are needed to verify whether (83) can be improved.

Finally, we illustrate the advantages of the proposed method in comparison to the classical Galerkin method in Figure 11. Specifically, we compare the performance of the Galerkin method using $\mathbb{P}_1$ elements with the MHM method (with $\ell = 0$) for the problem in the Section 5.2.2 with $\epsilon = a = 1$, and considering different choices for the second level approximation. The idea is to measure how the condition number of the matrices underlying the Galerkin and MHM methods behaves in relation to the error you want to achieve. In all cases (see Figure 11), we observe a clear advantage of the MHM method. Note that for a given target error, this advantage arises when we refine $H$, which impacts the condition number with a rate $H^{-1}$ in the case of the MHM method, in contrast to a rate of $h^{-2}$ observed in the Galerkin method.

## 6. Conclusions

The present work proposed a general methodology to construct and analyze a new family of MHM methods for the RAD model for polytopal elements, where the methods in [25] and [7] are particular cases. In this sense, the present work is a companion article to [25], which proves the well-posedness and optimal convergence of the method proposed in [25], while also extending it to polytopal meshes. The functioning of the upscaling procedure in the novel MHM methods depends on the balance among the diffusion, reaction, and advection coefficients at the elementary level, impacting the mathematical structure of the method. Notably, the method recovers the original mixed form of the MHM method for the Poisson model in [2] since the multiscale basis functions satisfy local Neumann problems with a zero mean value constraint. On the other hand, when the reaction coefficient is present (but the advective term is absent), the local problems are still of the Newman type, and the MHM method induces a coercive operator in the $H^1$ space. The presence of an advection field in the model is more involved. In this case, the MHM method relies on local RAD problems with Robin boundary conditions. Since the mathematical structure of the MHM method depends on the well-posedness of local problems, whenever a local RAD problem is identified as ill-posed, the MHM method becomes of mixed type, which recovers the well-posedness of the local problem.
<figure><img src="image_12.png" /><figcaption>Figure 10: Study of condition number related to the global system in terms of $h$ for Case 2 (left) and in terms of $H$ with $\ell = 0$, $k = 2$, and $h = H$ (right). Here $\epsilon = \omega = 1$.</figcaption></figure>

Such an adaptive mechanism naturally impacts the matrix conditioning of linear algebraic systems (local and global) underlying the MHM method. We prove theoretical upper bounds of the condition number for the local and global matrices associated with the method and study their asymptotic behavior from the point of view of physical coefficients and different spatial scales. The main conclusions are that the associated global matrix is better conditioned (i.e. $O(H^{-1})$ instead of $O(H^{-2})$ for the classical Galerkin method), when we use the method in its space-based form ($H$ fixed, $H \to 0$) and is also insensitive to local refinement of sub-meshes ($h \to 0$). This last feature is particularly attractive when dealing with multi-scale problems where the local mesh must be refined.

Furthermore, we prove that the new MHM methods are well-posed under the same compatibility conditions between the interpolation spaces presented in [7] and [22] and achieve optimal convergence in terms of the skeletal mesh diameter assuming local regularity. Interestingly, such proofs did not follow the original strategy proposed in [2] but are the result of the equivalence between the two-level MHM methods and a discrete version of the primal hybrid formulation following the idea proposed in [6, 26] for the Poisson equation.

Extensive numerical tests validated the theoretical results on meshes containing L-type and hexagonal elements, for example, applied to problems with solutions containing boundary layers and oscillatory behaviors in two-dimensional and three-dimensional domains, showing the robustness of the methods in terms of physical parameters. On the other hand, the theoretical demonstration of such robustness with regard to the impact of physical parameters on error estimates remains an open problem and deserves further investigation.

## Acknowledgments

The first author was partially supported by ANID-Chile through grant AFB170001 of the PIA Program: Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal and grant FONDECYT-1211649. The third author was partially supported by Project EOLIS (MATH-AMSUD 21-MATH-04) and ANID-Chile through grant FONDECYT-1181572. The fourth author was partially supported by CNPq/Brazil No. 309173/2020-5 and FAPERJ E-26/201.182.2021, Project EOLIS (MATH-AMSUD 21-MATH-04) and Inria/France under the Inria International Chair.
<figure><img src="image_13.png" /><figcaption>Figure 11: Error versus condition number of the matrix associated with the Galerkin method and the global problem matrix associated with the MHM method.</figcaption></figure>

## References

[1] Araya, R., Barrenechea, G.R., Franca, L.P., Valentin, F., 2009. Stabilization arising from PGEM: A review and further developments. Appl. Numer. Math. 59, 2065-2081. doi:10.1016/j.apnum.2008.12.004.

[2] Araya, R., Harder, C., Paredes, D., Valentin, F., 2013. Multiscale Hybrid-Mixed Method. SIAM J. Numer. Anal. 51, 3505-3531. doi:10.1137/120888223.

[3] Arbogast, T., Pencheva, G., Wheeler, M.F., Yotov, I., 2007. A multiscale mortar mixed finite element method. Multiscale Model. Simul. 6, 319-346. doi:10.1137/060662587.

[4] Ayuso, B., Marini, L.D., 2009. Discontinuous Galerkin methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. 47, 1391-1420. doi:10.1137/080719583.

[5] Babuška, I., Osborn, J.E., 1983. Generalized finite element methods: Their performance and their relation to mixed methods. SIAM J. Numer. Anal. 20, 510-536. doi:10.1137/0720034.

[6] Barrenechea, G.R., A. Gomes, A.T., Paredes, D., 2022. A Multiscale Hybrid Methods. URL: https://hal.archives-ouvertes.fr/hal-03907060. preprint.

[7] Barrenechea, G.R., Jaillet, F., Paredes, D., Valentin, F., 2020. The multiscale hybrid mixed method in general polygonal meshes. Numer. Math. 145, 197-237. doi:10.1007/s00211-020-01103-5.

[8] Brenner, S.C., 2003. Poincaré-Friedrichs inequalities for piecewise $H^1$ functions. SIAM J. Numer. Anal. 41, 306-324. doi:10.1137/S0036142902401311.

[9] Brenner, S.C., Sung, L.Y., 2018. Virtual element methods on meshes with small edges or faces. Math. Models Methods Appl. Sci. 28, 1291-1336. doi:10.1142/S0218202518500355.

[10] Brezzi, F., Franca, L.P., Russo, A., 1998. Further considerations on Residual-Free Bubbles for advective-diffusive equations. Comput. Methods Appl. Mech. Engrg. 166, 25-33. doi:10.1016/S0045-7825(98)00080-2.

[11] Brezzi, F., Russo, A., 1994. Choosing bubbles for advection-diffusion problems. Math. Models Methods Appl. Sci. 4, 571-587. doi:10.1142/S0218202594000327.

[12] Brooks, A.N., Hughes, T.J.R., 1982. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg. 32, 199-259. doi:10.1016/0045.

[13] Chaumont-Frelet, T., Ern, A., Lemaire, S., Valentin, F., 2022a. Bridging the multiscale hybrid-mixed and multiscale hybrid high-order methods. ESAIM: M2AN 56, 261-285. doi:10.1051/m2an/2021082.

[14] Chaumont-Frelet, T., Paredes, D., Valentin, F., 2022b. Flux approximation on unfitted meshes and application to multi-scale hybrid-mixed methods. URL: https://hal.inria.fr/hal-03834748. preprint.

[15] Cicuttin, M., Ern, A., Lemaire, S., 2018. A hybrid high-order method for highly oscillatory elliptic problems. Comput. Methods Appl. Math. 19, 723-748. doi:10.1515/cmam-2018-0013.

[16] E, W., Engquist, B., 2003. The heterogeneous multiscale methods. Comm. Math. Sci. 1, 87-133. doi:10.4310/CMS.2003.v1.n1.a8.
[17] Efendiev, Y., Hou, T., Wu, X., 2000. Convergence of a nonconforming multiscale finite element method. SIAM J. Numer. Anal. 37, 888-910. doi:10.1137/S0036142997330329.

[18] Ern, A., Guermond, J.L., 2004. Theory and Practice of Finite Elements. Number 159 in Applied Mathematical Sciences, Springer, New York, NY. doi:10.1007/978-1-4757-4355-5.

[19] Fernando, H., Harder, C., Paredes, D., Valentin, F., 2012. Numerical multiscale methods for a reaction-dominated model. Comput. Methods Appl. Mech. Engrg. 201-204, 228-244. doi:10.1016/j.cma.2011.09.007.

[20] Franca, L., Frey, S., Hughes, T.J., 1992. Stabilized finite element methods. I. Application to the advective-diffusive model. Comput. Methods Appl. Mech. Engrg. 95, 253-276. doi:10.1016/0045-7825(92)0143-8.

[21] Franca, L.P., Valentin, F., 2000. On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation. Comput. Methods Appl. Mech. Engrg. 190, 1785-1800. doi:10.1016/S0045-7825(00)00190-0.

[22] Gomes, A.T.A., Pereira, W.S., Valentin, F., 2022. The MHM method for linear elasticity on polytopal meshes. IMA J. Numer. Anal. doi:10.1093/imanum/drac041.

[23] Grisvard, P., 2011. Elliptic Problems in Nonsmooth Domains. Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia. doi:10.1137/1.9781611972030.

[24] Harder, C., Paredes, D., Valentin, F., 2013. A family of multiscale hybrid-mixed finite element methods for the Darcy equation with rough coefficients. J. Comput. Phys. 245, 107-130. doi:10.1016/j.jcp.2013.03.019.

[25] Harder, C., Paredes, D., Valentin, F., 2015. On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogenous coefficients. Multiscale Model. Simul. 13, 491-518. doi:10.1137/130938499.

[26] Harder, C., Valentin, F., 2016. Foundations of the MHM method, in: Barrenechea, G.R., Brezzi, F., Cangiani, A., Georgoulis, E.H. (Eds.), Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Springer, Switzerland. pp. 401-433. doi:10.1007/978-3-319-41640-3.

[27] Hughes, T.J.R., Feijoo, G.R., Mazzei, L., Quincy, J., 1998. The variational multiscale method - a paradigm for computational mechanics. Comput. Methods Appl. Mech. Engrg. 166, 3-24. doi:10.1016/S0045-7825(98)00079-6.

[28] Korotov, S., Krížek, M., 2013. On simplicial red refinement in three and higher dimensions. URL: https://dml.cz/handle/10338.dmlcz/702939.

[29] Madureira, A.L., Sarkis, M., 2021. Hybrid localized spectral decomposition for multiscale problems. SIAM J. Numer. Anal. 59, 829-863. doi:10.1137/20M1314896.

[30] Mälqvist, A., Peterseim, D., 2014. Localization of elliptic multiscale problems. Math. Comp. 83, 2583-2603. doi:10.1090/S0025-5718-2014-02868-8.

[31] Raviart, P.A., Thomas, J.M., 1977a. A mixed finite element method for 2nd order elliptic problems, in: Mathematical Aspects of Finite Element Methods. Springer Berlin Heidelberg, Berlin. volume 606 of Lecture Notes in Mathematics, pp. 292-315. doi:10.1007/BFb0064470.

[32] Raviart, P.A., Thomas, J.M., 1977b. Primal hybrid finite element methods for 2nd order elliptic equations. Math. Comp. 31, 391-413. doi:10.1090/S0025-5718-1977-0431752-8.

[33] Sangalli, G., 2003. Capturing small scales in elliptic problems using a Residual-Free Bubbles finite element method. SIAM Multiscale Model. Simul. 1, 485-503. doi:10.1137/S1540345902411402.

[34] Sayas, F.J., Brown, T.S., Hassell, M.E., 2019. Variational techniques for elliptic partial differential equations. Theoretical tools and advanced applications. CRC Press, Boca Raton, FL. doi:10.1201/9780429507069.

[35] Verfürth, R., 2013. A posteriori error estimation techniques for finite element methods. Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford. doi:10.1093/acprof:oso/9780199679423.001.0001.

[36] Xu, J., Zikatanov, L., 2003. Some observations on Babuška and Brezzi theories. Numer. Math. 94, 195-202. doi:10.1007/s002110100308.