An agglomeration-based massively parallel non-overlapping additive Schwarz preconditioner for high-order discontinuous Galerkin methods on polytopic grids
P. F. Antonietti P. Houston G. Pennesi E. Süli


Out-of-equilibrium operation of a quantum heat engine:
 The cost of thermal coupling control
    Joachim Ankerhold
    December 30, 2023
============================================================================================




In this article we design and analyze a class of two-level non-overlapping additive Schwarz preconditioners for the solution of the linear system of equations stemming from discontinuous Galerkin discretizations of second-order elliptic partial differential equations on polytopic meshes. The preconditioner is based on a coarse space and a non-overlapping partition of the computational domain where local solvers are applied in parallel. In particular, the coarse space can potentially be chosen to be non-embedded with respect to the finer space; indeed it can be obtained from the fine grid by employing agglomeration and edge coarsening techniques. We investigate the dependence of the condition number of the preconditioned system with respect to the diffusion coefficient and the discretization parameters, i.e., the mesh size and the polynomial degree of the fine and coarse spaces. Numerical examples are presented which confirm the theoretical bounds.










§ INTRODUCTION

The process of defining a computational grid characterized by standard triangular/tetrahedral or quadrilateral/hexahedral-shaped elements is one of the potential bottlenecks when traditional finite element methods are employed for the numerical approximation of problems characterized by strong complexity of the physical domain, such as, for example, in geophysical applications, fluid-structure interaction, or crack propagation problems. In order to overcome this issue, during the last decade a wide strand of literature has focused on the design of numerical methods that support the use of computational meshes composed of general po­ly­gonal and polyhedral elements. In the conforming setting we mention, for example, the Composite Finite Element Method <cit.>, the Mimetic Finite Difference
Method <cit.>, the Polygonal Finite Element Method <cit.>, the Extended Finite Element Method <cit.>, the Virtual Element Method 
<cit.>, and the Hybrid High-Order
Method <cit.>. A major issue in the design of conforming methods on such general polytopic meshes is the definition of a suitable space of continuous piecewise polynomial functions; in this context, this is far from being a trivial task, particularly for high-order approximations. An alternative strand of literature has focused on the non-conforming setting, where the ease of defining spaces of piecewise polynomial functions is naturally associated with the flexibility provided by polytopic meshes. Here, we mention, for example, Hybridizable Discontinuous Galerkin Methods <cit.>, non-conforming 
Virtual Element Methods <cit.>, Gradient Schemes <cit.>, and Discontinuous Galerkin (DG) Methods <cit.>. In particular, DG methods represent a class of powerful non-conforming numerical schemes in which the use of numerical grids characterized by general polytopic elements couples very well with the possibility to build the underlying discrete space in the physical frame, thereby avoiding the need to map polynomial spaces from a reference/canonical element.

However, as was shown in <cit.>, the condition number of the matrix in a system of linear equations stemming from DG methods may be prohibitively large; indeed, by writing h to denote the mesh-size and p the polynomial degree, the condition number of the DG approximation to Poisson's equation grows like [O](p^4/h^2) as h tends to zero and p tends to infinity. For this reason, in recent years, the development of fast solvers and preconditioners for systems of linear equations stemming from (high-order) DG discretizations has been an active area of research. A variety of two-level and multigrid/multilevel techniques have been proposed, both in the geometric and algebraic settings, for the solution of DG discretizations; see, for example, <cit.>. In particular, the availability of efficient geometric multilevel solvers is strongly related to the possibility of employing general-shaped polytopic grids; indeed, if polytopic grids can be employed, then the sequence of grids which are required within a multilevel iteration can be defined by agglomeration; see <cit.> for details. Besides multigrid, a recent strand in the literature has focused on Schwarz domain decomposition methods; see, for example, <cit.>, for a general abstract overview of these methods. In the DG setting where standard triangular/tetrahedral or quadrilateral/hexahedral grids are employed, one of the first contributions in terms of domain decomposition solvers was presented for the solution of elliptic problems in <cit.>, where bounds of order [O](H/δ) and [O](H/h) were obtained for the condition number of the preconditioned system in the framework of overlapping and non-overlapping Schwarz methods, respectively; here, H, h, and δ represent the size of the coarse grid, the fine grid, and the amount of overlap, respectively. 
Dryja and Sarkis proposed in <cit.> an additive Schwarz preconditioner for the solution of second-order elliptic problems with discontinuous coefficients. There, the authors showed that the condition number of the of the matrix of the preconditioned system is independent of the jumps of the coefficients across the substructure boundaries and outside a thin layer along the substructure boundaries. A further development of this algorithm, which is very well suited for parallel computation, can be found in <cit.>. 
Concerning the setting of high-order DG methods, we mention the work in <cit.>, where additive and multiplicative Schwarz preconditioners were introduced for efficiently solving systems of linear equations arising from the discretization of second-order symmetric elliptic boundary-value problems using hp-version DG methods; there, hp-spectral bounds of order [O](σ p^2 H/h) were derived for a class of domain decomposition preconditioners for DG discretizations, where σ is the coefficient of the interior penalty stabilization parameter, p is the polynomial approximation degree, and H, h is the size of the coarse and fine mesh, respectively. Recently, in <cit.>, this condition number estimate was improved to yield the optimal rate of [O](σ p^2 H/q h), where q denotes the polynomial approximation degree employed within the coarse grid solver; cf. <cit.> for related work. We also mention the recent work presented in <cit.>, where the influence of the penalty terms, as well as the choice of coarse mesh spaces, on the condition number of the matrix of the system of linear equations preconditioned with additive Schwarz methods were investigated.

The goal of this article is to design and analyze a class of two-level non-overlapping additive Schwarz preconditioners for hp-version DG discretizations of second-order elliptic problems on general polytopic grids. Given the DG discrete problem defined on a fine mesh of granularity h, the preconditioner is designed by introducing two additional partitions employed to define the local solver operators and the coarse space correction. On the one hand, the partition employed to build the local solvers is related to a suitable splitting of the DG space and hence it is assumed to be nested with respect to the fine polytopic mesh. On the other hand no conditions are imposed on the coarse partition, which can be non-nested with respect to the fine grid. In particular, we consider the massively parallel setting, whereby the partition employed for the local solvers are finer than the grid employed within the coarse solver. Here, we investigate the dependence of the condition number with respect to both the diffusion coefficient and the discretization parameters of the fine and coarse spaces, as well as the granularity of the partition for the local solver. We stress that our analysis is carried out in a very general setting, and, in particular, for nested meshes, it allows the computational domain where the model problem is posed to be non-convex; cf. <Ref> below.

The rest of the paper is organized as follows. In <Ref> we introduce the DG method on polytopic grids for the numerical approximation of second-order elliptic problems. In <Ref> we formulate the additive Schwarz preconditioner analyzed in this article. In <Ref> we then outline some key analytical results that are required for the analysis that follows. <Ref> is devoted to deriving some preliminary results required to obtain the desired bound on the condition number of the matrix of the preconditioned system stated in <Ref>. Numerical experiments are presented in <Ref> to confirm the theoretical bounds derived in this article.








§ DG METHOD ON POLYTOPIC GRIDS

In this article we consider the following second-order elliptic problem. Let 
⊂ℝ^d, d=2,3, be a polygonal/polyhedral domain with boundary ∂Ω and let f ∈ L^2() be a given function. We consider the following weak formulation: find u ∈ V = H^1_0()  such that

    𝒜(u,v) := ∫_Ωρ∇ u ·∇ v   = ∫_Ω f v    ∀  v ∈ V.

Here, ρ∈ L^∞() denotes the diffusion coefficient, which we assume to be such that 0< ρ_0 ≤ρ; here, we can assume that ρ_0 = 1, since (<ref>) can always be scaled by 1/ρ_0. Throughout this article, we use the notation x ≲ y to signify that there exists a positive constant C, independent of the diffusion coefficient ρ and the discretization parameters, such that x ≤ Cy. Similarly we write x ≳ y in lieu of x ≥ Cy, while x  y is used if both x ≲ y and x ≳ y hold.

Let [T][h] be a tessellation of  consisting of disjoint polytopic elements  of diameter h_ such that = ∪. Here, we denote by [F][h] the set of faces F, which are defined as the (d-1)-dimensional planar facets of the elements κ present in the mesh [T][h]. For d=3, we assume that each planar face of an element ∈[T][h] can be subdivided into a set of co-planar (d-1)-dimensional simplices and we refer to this set as the set of faces, cf. <cit.>. Moreover, we write [F][h]^B := { F ∈[F][h]: F ⊂∂} to denote the set of boundary faces and [F][h]^I := [F][h] ∖[F][h]^B the set of interior faces. We set h :=max_∈[T][h] h_ and, to ease the presentation, we assume that h  h_ for all ∈[T][h].

We adopt the hypothesis that the diffusion coefficient ρ is piecewise constant on each polytopic element ∈[T][] and write ρ_ = ρ|_ to denote its restriction to ; we refer to <cit.> for the more
general case when ρ violates this assumption.

We assume that [T][h] satisfies the following assumptions; cf. <cit.> for details.
 
For any ∈[T][h] there exists a set of non-overlapping d-di­mensional simplices _♭^F ⊂, for F ⊂∂, such that for any face F ⊂∂, we have that F = ∂∩∂_♭^F, ⋃_F ⊂∂_♭^F⊂, and the diameter h_ of  can be bounded by
h_≲d |_♭^F |/|F| for all F ⊂∂,
where |F| and |_♭^F| denote the Hausdorff measure of F and _♭^F, respectively.	 


We assume that there exists a covering [T][h]^#:= {𝒮_}_ of [T][h] consisting of shape-regular d-dimensional simplices 𝒮_, such that, for any ∈[T][h], there is an 𝒮_∈[T][h]^# satisfying ⊂𝒮_ and diam([S][]) ≲ h_. We also assume that

    max_∈[T][h]card{' ∈[T][h]: ' ∩𝒮_∅ , 𝒮_∈[T][h]^# such that ⊂𝒮_}≲ 1.



We write ⟨·⟩ to denote the harmonic average operator defined as follows: let η be a sufficiently smooth function; then, for any F ⊂∂, F ∈[F][h], we define

    ⟨η⟩|_F := 2 η_^+ η_^-/η_^+ + η_^-,     F∈[F][h]^I, F ⊂∂^+ ∩∂^-,
    η_,     F∈[F][h]^B, F ⊂∂∩∂Ω,

where η_, κ∈[T][h], denotes the trace of η on ∂.
Moreover, for sufficiently smooth vector- and scalar-valued functions τ and v, respectively, we define the following jump and weighted average operators across F ∈[F][]:
for F ⊂∂^+ ∩∂^-, F∈[F][h]^I, we write

    τ   := τ_^+·𝐧_^+ + τ_^-·𝐧_^-,    τ_ω   := ωτ_^+  + (1-ω) τ_^-,
    v   := v_^+𝐧_^+ + v_^-𝐧_^-,
       v_ω   := ω v_^+  + (1-ω) v_^-,

where 𝐧_^± denotes the unit outward normal vector to ^±, respectively. For  F ⊂∂, F∈[F][h]^B, we set τ_ω :=τ_, v := v_𝐧, cf. <cit.>. Here, ω∈ [0,1] represents the weight employed for the definition of ·_ω. Given [T][h]^s ⊆[T][h] and [F][h]^s ⊆[F][h], we write
	∫_[T][h]^sdx := ∑_∈[T][h]^s∫_dx,
	| v |_H^1([T][]^s)^2 := ∑_∈[T][]^s∇ v _L^2()^2,
	∫_[F][h]^sds := ∑_F∈[F][h]^s∫_F ds, and
	v _L^2([F][]^s)^2 := ∑_F ∈[F][]^s v _L^2(F)^2.
Finally, writing [P][p]() to denote the space of polynomials of total degree p ≥ 1 on , ∈[T][], the DG space is given by

    V_h :={v_h ∈ L^2():v_h|_∈[P][p]()  ∀ ∈[T][h]}.

With this notation, we introduce the Symmetric Interior Penalty DG (SIPDG) discretization of (<ref>), cf. <cit.>: find u_∈ V_ such that

    [h][u_h][v_h] = ∫_ f v_hdx  ∀  v_h∈ V_h,

where 𝒜_h: V_h × V_h →ℝ is the bilinear form of the DG method defined as


    [h][u_h][v_h]    := ∫_ρ[ ∇_h u_h ·∇_h v_h +  ∇_h u_h ·ℛ_ρ(v_h) + ∇_h v_h ·ℛ_ρ(u_h ) ] d
          + ∫_[F][h]σ_h,ρu_h·v_h ds,

and ∇_h denotes the piecewise gradient operator on [T][h]. Here, ℛ_ρ:[L^1([F][h])]^d → [V_h]^d denotes the lifting operator defined by

    ∫_[R][ρ](𝐪) ·η   := - ∫_[F][h]𝐪·η_ω ds   ∀ η∈ [V_h]^d,

where, we take ω|_F := ρ_^-/ρ_^+ + ρ_^- on each internal face F∈[F][h]^I shared by ^±.
In (<ref>), according to <cit.>, σ_h,ρ∈ L^∞([F][h]) denotes the interior penalty stabilization function, which is defined by

    σ_h,ρ|_F := C_σρ_ p^2 / h_  ∀  F ∈[F][h],

with C_σ>0 independent of ρ, p, |F|, || and h_. It is well known that the condition number of the operator [A][h] is potentially prohibitively large and depends on the size of the partition [T][h] and the polynomial degree p employed for the discretization; cf. <cit.> for standard triangular/tetrahedral/hexahedral grids and <cit.> for polytopic grids. Our goal is to introduce a massively parallel non-overlapping additive Schwarz preconditioner, which can be employed as a preconditioner to accelerate the convergence of iterative solvers such as the Conjugate Gradient method. 






§ NON-OVERLAPPING ADDITIVE SCHWARZ PRECONDITIONER

The definition of the additive Schwarz preconditioner requires the introduction of two additional partitions (besides [T][h]): a partition [T][] composed of disjoint polytopic subdomains where local solvers are applied in parallel and a non-overlapping partition [T][] employed for the coarse space correction. To this end, we introduce the following notation:


  * [T][] := {_1, …, _N_} of size := max_1 ≤ i ≤ N_{diam(_i)} such that each subdomain _i is the union of some elements ∈[T][]; we assume that diam(_i) for all i=1,…,N_. We also assume that a colouring property holds, i.e., there exists a positive integer N_𝕊 such that

    max_i=1,…,N_card{_j: ∂_i ∩∂_j ∅}≤ N_𝕊,

i.e., N_𝕊 represents the maximum number of neighbours that any subdomain _i ∈[T][] may possess.

  * [T][] := {[D]_1, …, [D]_N_} of size := max_1 ≤ j ≤ N_{diam([D]_j)} such that diam([D]_j) for all j=1,…,N_.

We remark that  the grids [T][] and [T][] are possibly non-nested, cf. <Ref>.



 
Given that [T][] is defined by agglomeration of fine-grid-elements ∈[T][], we write [T][] ⊆[T][], since, for all ∈[T][], there exists a [K] ∈[T][] such that ⊆[K]. However, we point out that no further assumptions are needed on the relationship between [T][] and [T][] for the definition of our method. Classical additive Schwarz methods have typically been defined based on the assumption that [T][] ⊆[T][] ⊆[T][]. In this article we take a different approach: firstly, we assume that the granularity of [T][] is finer than that of [T][]; indeed, we are particularly interested in the massively parallel case whereby [T][] = [T][], cf. <cit.>. Secondly, we also permit the use of non-nested coarse and fine partitions, i.e., when [T][] ⊈[T][].


The main ingredients of the additive Schwarz method are defined as follows.


Local Solvers. Consider the subdomain partition [T][] with cardinality N_. Then, for each subdomain _i ∈[T][] we define a local space V_i as the restriction of the DG finite element space V_ to _i, i.e., for i=1,…,N_,

    V_i := { v_i ∈ L^2(_i): v_i|_∈[P][p]() ∀ ∈[T][h], ⊆_i }≡ V_|__i.

The associated local bilinear form on V_i× V_i is defined by

    𝒜_i: V_i × V_i →ℝ,   𝒜_i(u_i,v_i) := [][R_i^⊤ u_i][R_i^⊤ v_i]   ∀  u_i,v_i ∈ V_i,

where R_i^⊤:V_i → V_ denotes the classical extension-by-zero operator from the local space V_i to the global space V_. The restriction operator R_i:V_→ V_i, i=1,…,N_, is defined as the transpose of R_i^⊤ with respect to the L^2(_i) inner product.


Coarse Solver. For 1 ≤ q ≤ p, the coarse solver is defined on the partition [T][]. To this end, let V_0 be the DG finite element space defined on [T][] given by

    V_0 ≡ V_ := {v_∈ L^2(): v_|_[D]_j∈[P][q]([D]_j), j=1,…,N_}.

Further, let R_0^⊤ be the L^2-prolongation operator from V_0 to V_, defined as:

    R_0^⊤: v_0 ∈ V_0 ⟼ R_0^⊤ v_0 ∈ V_:   ∫_ R_0^⊤ v_0 w_ :=  ∫_ v_0 w_  ∀  w_∈ V_.

In this way R_0^⊤ is well defined also when [T][] and [T][] are non-nested. Then, the bilinear form associated to V_0 is defined by

    𝒜_0: V_0 × V_0 →ℝ,   𝒜_0(u_0,v_0) := 𝒜_(R_0^⊤ u_0, R_0^⊤ v_0)   ∀  u_0,v_0 ∈ V_0.









When V_0 ≡ V_⊆ V_, i.e., when the coarse and fine grids [T][] and [T][], respectively, are nested, then the action of R_0^⊤ on a coarse function coincides with the action of the natural injection operator. Indeed, by contradiction, if there exists a v_0 ∈ V_0 such that R_0^⊤v_0 v_0, then, by employing the definition of R_0^⊤, we have

    0 <  R_0^⊤v_0 - v_0 _L^2() = min_w_h ∈ V_h w_h - v_0 _L^2()≤v_0 - v_0 _L^2() = 0,

which is a contradiction and hence R_0^⊤ v_0 = v_0 for all v_0 ∈ V_0 when V_0 ⊆ V_.

Introducing the projection operators P_i := R_i^⊤P_i: V_→ V_, i=0,1,…,N_, where

    P_i:V_→ V_i,  𝒜_i(P_i v_, w_i) := [h][v_][R_i^⊤ w_i]   ∀  w_i ∈ V_i,  i=1,…,N_,
       P_0:V_→ V_0,  𝒜_0(P_0 v_, w_0) := [h][v_][R_0^⊤ w_0]   ∀  w_0 ∈ V_0,

the additive Schwarz operator is defined by P_ad := ∑_i=0^N_ P_i. For an upper bound on the condition number of P_as, we refer to <Ref> below.


















§ ANALYTICAL BACKGROUND
 

Before we embark on developing the preliminary results needed to analyze the condition number of the additive Schwarz operator introduced in the previous section, we first state two key theorems, which are essential in the forthcoming analysis, and which may be of relevance in other more general settings. To this end, we write D⊂ℝ^d, d≥ 2, to denote a bounded, open, simply connected domain, with boundary ∂ D; in the proceeding analysis, D will be selected to be an element 𝒟_j, j=1,…,N_H, from the coarse mesh [T][]. Following <cit.>, cf. also <cit.>, we introduce the definition of a special Lipschitz domain, and the notion of a domain with a minimally smooth boundary.

(<cit.>,<cit.>)
	Let ϕ:ℝ^d-1→ℝ be a function that satisfies the Lipschitz condition
	
    |ϕ() - ϕ()| ≤ M |-|     ∀ ,∈ℝ^d-1.

	The smallest M for which (<ref>) holds is denoted by C_ϕ.
	Based on this function, we define the special Lipschitz domain it determines to be the set of points lying above the hypersurface y=ϕ() in ℝ^d, i.e.,
	
    ω = {∈ℝ^d: x_d > ϕ(x_1,x_2,…,x_d-1)}.

	The Lipschitz constant of the domain ω is defined by C_ω:=C_ϕ.

Equipped with this definition, we now introduce the concept of a minimally smooth boundary.
(<cit.>,<cit.>) 
The boundary ∂ D of D is said to be minimally smooth if there exists an ϵ>0, a positive integer N, an M>0, and a sequence U_1,U_2,… of open sets such that:

	
  * If ∈∂ D, then B(,ϵ) ⊂ U_i for some i, where B(,ϵ) denotes the ball with centre  and radius ϵ;
	
  * No point ∈ℝ^d is contained in more than N of the U_i's;
	
  * For each i there exists a special Lipschitz domain ω_i with C_ω_i≤ M such that
	
    U_i∩ D =U_i ∩ω_i.



Based on the previous definition, we now introduce the following classical extension operator.
(<cit.>, <cit.>) 
Let D be a domain with mi­ni­mally smooth boundary. Then, there exists a linear extension
operator E:H^s(D) → H^s(ℝ^d), s ∈ℕ_0, such that
Ev|_D=v and

    E v _H^s(ℝ^d)≤ C_E v _H^s(D),

where C_E is a positive constant depending only on s and the constants ϵ, N, and M defined in Definition <ref>, which characterize the boundary ∂ D.


We highlight that, crucially, the constant C_E appearing in Theorem <ref> is independent of the measure of the underlying domain D.


Secondly, we now study the regularity of the following Neumann boundary-value problem: find z such that

    - Δ z  = f    in  D,     ∇ z ·𝐧 = 0    on ∂ D,

with f ∈ L^2_0(D) := { v ∈ L^2(D) : ∫_D v  = 0 }.
 
	Let D be a bounded, open, convex (and therefore Lipschitz) domain. Then, there exists a unique solution z ∈ H^2(D) ∩ L^2_0(D) to the homogeneous Neumann problem (<ref>). Moreover, the following stability bound holds:
	
    z_H^2(D)≤√(6) ( f_L^2(D) + 1/π(D) ∇ z _L^2(D)).




	We defer the proof to Appendix <ref>.




An analogous bound to (<ref>), with constant equal to unity, has been derived in <cit.> for convex polygonal domains in ℝ^2; there the term
1/π(D) ∇ z_L^2(Ω) on the right-hand side of (<ref>) is replaced by z _H^1(D).










§ PRELIMINARY RESULTS

We first present some preliminary results, which will be employed within the analysis contained in <Ref>. For the sake of simplicity of the presentation, we assume that the grids [T][] and [T][] are nested; the extension of the theoretical analysis to the general case [T][] ⊈[T][] is deferred to <Ref>. 
Here, we introduce the following energy norm:

    [w][h,ρ][2]:=∫_[T][h]ρ |∇_h w|^2  + ∫_[F][h]σ_h,ρ|w|^2 ds.

The well-posedness of problem (<ref>) with respect to the norm (<ref>) is then established in the following lemma, cf. <cit.>.

Suppose that [T][h] satisfies <Ref>; then,

[h][u_h][v_h]   ≲[u_h][h,ρ][v_h][h,ρ]        ∀  u_h,v_h ∈V_h, 

[h][u_h][u_h]   ≳[u_h][h,ρ][2]        ∀  u_h ∈V_h. 

The second bound holds provided that C_σ appearing in (<ref>) is sufficiently large.

We also recall the following trace inequality on polytopic domains, introduced in <cit.>.

Suppose that [T][h] satisfies <Ref>; then, the following inequality holds:

    [v][2][∂][2]≲p^2/h_ [v][2][][2]   ∀  v ∈[P][p]() ∀ ∈[T][].




Next we also recall a result regarding the approximation operator presented in <cit.> and <cit.>, to which we refer for details. However, for the purposes of this article, we consider a slight generalization: given a connected subdomain D⊆Ω, we assume that D is formed from the union of a subset of elements κ∈[T][h]; we denote the collection of such elements by [T][h,D], i.e., D̅ := ∪_κ∈[T][h,D]κ̅. With this definition, the approximant relies on the properties of the extension operator 𝔈: H^s(D) → H^s(ℝ^d), s ∈ℕ_0, introduced in Theorem <ref>. Hence, following <cit.> and <cit.>, we deduce the following lemma.

 Suppose that <Ref> is satisfied, and let v ∈ L^2(D) be such that, for some k ≥ 0, v|_∈ H^k() and 𝔈v|_[S][]∈ H^k([S][]) for ∈[T][h,D], with [S][] ∈[T]^#_h as defined in <Ref>. Then, there exists a projection operator Π_h,: L^2(D)→[P][p]() such that

    [v-Π_h, v][l][]    ≲h^s-l/p^k-l[𝔈 v][k][[S][]],    0 ≤ l ≤ k.
     v-Π_h, v _L^2(∂)   ≲h^s - 1/2/p^k - 1/2[𝔈 v][k][[S][]],    k ≥ 1,

where s:=min{p+1,k} and p ≥ 1.



Given that <Ref> holds for all ∈[T][h,D] we can define the global approximation operator Π_h: L^2(D)→ V_h|_D such that Π_h|_ = Π_h,. If v ∈ H^k(D), k ≥ 0, then, by noting <Ref>, together with Theorem <ref>, the following bound holds:

    [v-Π_h v][l][[T][h,D]] ≲h^s-l/p^k-l[v][k][D],    0 ≤ l ≤ k.




A key ingredient in our analysis is the conforming approximant defined in <cit.>. In particular, to ensure that the preconditioner is scalable in the presence of jumps in the diffusion coefficient, here we define the conforming approximant in a slightly different manner, in order to obtain an approximation of discontinuous discrete functions v_h ∈ V_h on each local domain [D][j] ∈[T][], j=1,…,N_. To this end, we adopt the following assumption on the coarse mesh [T][].

We assume that [D][j] ∈[T][] is a convex polytope with Lipschitz boundary ∂[D][j], for any j=1,…,N_, and that |[D][j]| ^d, j=1,…,N_, where | [D][j] | represents the Hausdorff measure of [D][j].

















For the sake of the analysis, we define the following local grids generated from [T][h] and [T][]:







    [T][h,j] := {∈[T][]: ⊂[D][j], [D][j] ∈[T][]},
      [F][h,j]^I := { F ∈[F][h]: F ⊂[D][j], [D][j] ∈[T][]}, 
    
    [F][h,j]^B := {F∈[F][h]: F ⊂∂[D][j], [D][j] ∈[T][]}, [F][h,j] := [F][h,j]^I ∪[F][h,j]^B,

for j=1,…,N_.

Note that since the grids [T][] and [T][] are nested, i.e., [T][h,j] ⊆[T][], j=1,…,N_, [T][h,j] also satisfies Assumptions <ref> and <ref>, for all j=1,…,N_.

The local conforming approximant is then defined as follows.

Let [D][j] ∈[T][] satisfy <Ref>, j=1,…,N_, and let the discrete gradient operator of v_h ∈ V_h inside [D][j] be defined by the equality 𝒢_h,j(v_h) = (∇_h v_h + [R][1]([J][j](v_))), j=1,…,N_, where  is the characteristic function on [D][j], while [J][j]: V_→ [L^1([F][h])]^d is defined by

    [J][j](v_)|_F := v_|_F   if F∈[F][h,j]^I,     [J][j](v_)|_F := 0   otherwise,

where [F][h,j]^I is the set defined in (<ref>). Here, [R][1]: [L^1([F][h])]^d → [V_h]^d is the lifting operator with ρ=1 and ω = 1/2 in its definition given in (<ref>). Then, v_h,j is defined as the solution of the following problem: find v_h,j∈V_j = H^1([D][j]) ∩Ł^2_0([D][j]) such that

    ∫_[D][j]∇v_h,j·∇ w  = ∫_[D][j]𝒢_h,j(v_h) ·∇ w   ∀  w ∈V_j.





 Since v_h,j∈V_j and [D][j] satisfies <Ref>, we note that, cf. <cit.>, v_h,j_L^2([D][j])≤ C_p ∇v_h,j_L^2([D][j]), with C_p ≲ diam([D][j])^1+d/2 | [D][j] |^-1/2. Thereby, invoking <Ref> gives

    v_h,j_H^1([D][j])^2 ≤ (1 + C_p^2) ∇v_h,j_L^2([D][j])^2≲∇v_h,j_L^2([D][j])^2,

where we also made use of the fact that diam([D][j]) ≲ 1.


By proceeding as in <cit.> we prove the following approximation result.

Let [D][j] ∈[T][] satisfy <Ref>, j=1,…,N_. Given v_h ∈ V_h, we take v_h,j∈V_j, j=1,…,N_, to be the conforming approximant given in <Ref> and we define

    v_h,j := v_h,j + 1/|[D][j]|∫_[D][j] v_h .

Then, the following approximation and stability bounds hold:


    [v_h -  v_h,j ][2][[D][j]]    ≲h/p[σ_h,1^1/2v_h][2][[F][h,j]^I],
    
    | v_h,j |_H^1()^2    ≲∇_h v_h _L^2([D][j])^2 + [σ_h,1^1/2v_h][2][[F][h,j]^I]^2,


for j=1,…,N_.















Given [D][j] ∈[T][], j=1,…,N_H, let z ∈ H^2([D][j]) ∩ L^2_0([D][j]) be the solution of problem (<ref>) posed on the domain [D][j], i.e., D=[D][j], with f=v_h   - v_h,j. Then, by employing integration by parts we obtain 

    v_h   - v_h,j_L^2([D][j])^2     = - ∫_[D][j] (v_h - v_h,j ) Δ z 
        = ∫_[D][j] (∇_h v_h - ∇v_h,j ) ·∇ z  - ∫_[F][h,j]^I∇ z ·v_h  s
        = - ∫_[D][j][R][1]([J][j](v_)) ·∇ z    -∫_[F][h,j]^I∇ z ·v_h  s,

where we have also used the facts that ∇v_h,j = ∇v_h,j, ∇ z ·𝐧|_F = 0 if F ⊂∂[D][j], v_h,j|_F = 0 for all F∈[F][h,j]^I, since v_h,j∈ H^1([D][j]), and that [G][h,j](v_h) = (∇_h v_h + [R][1]([J][j](v_)) together with the definition of v_h,j, cf. (<ref>). Using the definition of [R][1] and [J][j], cf. (<ref>) and (<ref>), respectively, for any z_h ∈ V_h|_[D][j], we have

    v_h  - v_h,j_L^2([D][j])^2
        = - ∫_[F][h,j]^I∇ z ·v_h  s - ∫_[D][j][R][1]([J][j](v_)) ·∇ z 
        = - ∫_[F][h,j]^I∇ z ·v_h  s  + ∫_[F][h,j]^Iv_h·∇_h z_h_1/2  s 
           - ∫_[D][j][R][1]([J][j](v_)) · ( ∇ z - ∇_h z_h ) 
        = ∫_[F][h,j]^Iv_h·∇_h z_h - ∇ z_1/2  s 
           - ∫_[D][j][R][1]([J][j](v_)) · ( ∇ z - ∇_h z_h ) ;

here we have used that, since z ∈ H^2([D][j]), ∇ z_1/2|_F = ∇ z|_F, F∈[F][h,j]^I. Hence, we get


    v_h - v_h,j_L^2([D][j])^2    ≲[R][1]([J][j](v_)) _L^2([D][j])∇ z - ∇_h z_h _L^2([D][j]) 
          + σ_h,1^1/2v_h_L^2([F][h,j]^I)σ_h,1^-1/2∇ z - ∇_h z_h _1/2 _L^2([F][h,j]^I).


The first term on the right-hand side of (<ref>) can be written as follows:

    [R][1]([J][j](v_) ) _L^2([D][j])^2    ≤∫_[R][1]([J][j](v_)) ·[R][1]([J][j](v_)) 
        = - ∫_[F][h] ([J][j](v_) )·[R][1]([J][j](v_)) _1/2  s 
        = - ∫_[F][h,j]^Iv_h·[R][1]([J][j](v_)) _1/2  s,


where we have also used the definitions of [R][1] and [J][j], cf. (<ref>) and (<ref>), respectively. Then, from (<ref>) and the Cauchy–Schwarz inequality we obtain the following bound:

    [R][1]    ([J][j](v_)) _L^2([D][j])^2 ≤σ_h,1^1/2v_h_L^2([F][h,j]^I) σ_h,1^-1/2[R][1]([J][j](v_)) _1/2_L^2([F][h,j]^I).

The second term on the right-hand side of (<ref>) can be bounded by invoking Lemma <ref> as follows:

    σ_h,1^-1/2[R][1]([J][j](v_)) _1/2_L^2([F][h,j]^I)^2    ≤∑_∈[T][h,j]σ_h,1^-1/2[R][1]([J][j](v_)) _L^2(∂)^2 
        = C_σ^-1∑_∈[T][h,j]h_/ p^2[R][1]([J][j](v_)) _L^2(∂)^2 
       ≲∑_∈[T][h,j][R][1]([J][j](v_)) _L^2()^2,

where we have used that h_≤ 2 h_ for any ∈[T][]. By inserting (<ref>) into (<ref>) we obtain

    [R][1]([J][j](v_)) _L^2([D][j])≲σ_h,1^1/2v_h_L^2([F][h,j]^I).

Hence, by selecting z_h = Π_h z in (<ref>) and employing (<ref>), <Ref>, cf., also, <Ref>,  and Theorem <ref>, we obtain

    v_h - v_h,j_L^2([D][j])^2    ≲h/p σ_h,1^1/2v_h_L^2([F][h,j]^I)  z _H^2([D][j])
    
       ≲h/p σ_h,1^1/2v_h_L^2([F][h,j]^I) ( v_h - v_h,j_L^2([D][j]) + ∇ z _L^2([D][j])).

Here, we note that z also solves

    ∫_[D][j]∇ z ·∇ w  = ∫_[D][j] (v_ - v_h,j) w   ∀  w ∈V_j,

from which, by choosing w = z, upon application of the Cauchy–Schwarz inequality and Poincaré's inequality, cf. Remark <ref>, we get that ∇ z _L^2([D][j])≲ v_ - v_h,j_L^2([D][j]). By inserting this bound into (<ref>) we obtain (<ref>). In order to show (<ref>) we first select w = v_h,j∈V_j in (<ref>); then, using the Cauchy–Schwarz inequality we obtain:

    | v_h,j |_H^1([D][j])≲[G][h,j](v_h) _L^2([D][j]).

Then, from the definition of [G][h,j] given in <Ref> we have:

    [G][h,j](v_h) _L^2([D][j])^2  ≲ ∇_h v_h _L^2([T][h,j])^2 + [R][1]([J][j](v_)) _L^2([D][j])^2.

The bound (<ref>) is then obtained by inserting (<ref>) into (<ref>).



We are now ready to investigate the relationship between the spaces V_, V_, and V_ introduced above. The following result concerns the approximation of a function v_h ∈ V_ with a coarse function v_∈ V_; this represents an extension of the analogous result presented in <cit.>.


For any v_h ∈ V_h there exists a coarse function v_∈ V_ such that

    v_h - R_0^⊤ v__L^2([D][j])≲/q ( ∇_h v_h _L^2([T][h,j])^2 + σ_h,1^1/2v_h_L^2([F][h,j]^I)^2)^1/2, 
        | v_h - R_0^⊤ v_ |_H^1([T][h,j])≲ ( ∇_h v_h _L^2([T][h,j])^2 + σ_h,1^1/2v_h_L^2([F][h,j]^I)^2)^1/2,

for j=1,…,N_, where [T][h,j] and [F][h,j]^I are as defined in (<ref>).



Let v_h ∈ V_h and define v_ by v_|_[D][j] := (1/|[D][j]|∫_[D][j] v_h ) + (Π_ (v_h,j))|_[D][j], j=1,…,N_, where v_h,j is as defined in <Ref> and Π_ denotes the global variant of the hp-approximant introduced in <Ref>, cf. also <Ref>, defined on the coarse space V_. Then, by noting <Ref>, the application of the triangle inequality gives


    v_h - R_0^⊤ v__L^2([D][j])   =    v_h - v__L^2([D][j])
       ≲ v_h - v_h,j_L^2([D][j])+ v_h,j - v__L^2([D][j])
       ≲ v_h - v_h,j_L^2([D][j])+ v_h,j - Π_(v_h,j)_L^2([D][j]),

with v_h,j as defined in <Ref>.
Employing <Ref> together with <Ref>, cf. <Ref>, gives

    v_h -  R_0^⊤ v__L^2([D][j])≲ v_h - v_h,j_L^2([D][j])+ /q v_h,j_H^1([D][j]).

By applying Poincaré's inequality to v_h,j∈V_j, see also <Ref>, and noting the bounds given in <Ref>, we immediately deduce inequality (<ref>) by observing that ≤ and q ≤ p. In order to obtain (<ref>) we proceed as follows:

    | v_h - R_0^⊤ v_ |_H^1([T][h,j])   ≲ |v_h|_H^1([T][h,j]) + |R_0^⊤ v_|_H^1([T][h,j])
        =  |v_h|_H^1([T][h,j]) + |v_|_H^1([T][h,j]).

Thanks to the triangle inequality and observing that v_|_[D][j]∈[P][q]([D][j]) ⊂ H^1([D][j]) we have that

    | v_ |_H^1([T][,j])   ≲ | Π_(v_h,j) - v_h,j |_H^1([D][j]) + | v_h,j |_H^1([D][j])
       ≲ | Π_(v_h,j) - v_h,j |_H^1([T][]) + | v_h,j |_H^1()≲ | v_h,j |_H^1(),

where we have used the bound stated in <Ref> and Poincaré's inequality, cf. also <Ref>. Inserting (<ref>) into (<ref>) and noting <Ref> gives (<ref>).



By summing over all [D][j] ∈[T][], j=1,…,N_, the local bounds of <Ref> give rise to the following global estimates:

    v_h - R_0^⊤ v__L^2()≲/q [v_h][h,1],
          | v_h - R_0^⊤ v_ |_H^1([T][h])≲[v_h][h,1],

which are in agreement with the analogous results developed in <cit.>.



Before proceeding with the analysis of P_ad we also need the following result regarding the properties of the subdomain decomposition introduced in <Ref>.


Given v_h ∈ V_h, there exists a unique decomposition, i.e., v_h = ∑_i=1^N_ R_i^⊤ v_i, with v_i ∈ V_i i=1,…,N_, such that

    [A][h](v_h,v_h) = ∑_i=1^N_[A][i](v_i,v_i) + ∑_i,j=1, i j^N_[A][h](R_i^⊤ v_i, R_j^⊤ v_j),

and

    | ∑_i,j=1, i j^N_[A][h](R_i^⊤ v_i, R_j^⊤ v_j) |    ≲[√(ρ) ∇_h v_h][2][[T][h]]^2 + ∑_i=1^N_σ_h,ρ^1/2 v_h _L^2(∂_i)^2.




Given v_h ∈ V_h set v_i := R_i v_h, i=1,…,N_; then,
[A][h](R_i^⊤ v_i, R_j^⊤ v_j) = 0 if ∂_i ∩∂_j = ∅.
For i,j=1,…,N_, we have

    |  ∑_i,j=1, i j^N_[A][i](R_i^⊤ v_i,R_j^⊤ v_j) |  ≲∑_i,j=1, i j^N_ | [A][i](R_i^⊤ v_i, R_j^⊤ v_j) |   ∀  v_i ∈ V_i, v_j ∈ V_j.

Now let i  j be such that ∂_i ∩∂_j ∅ and write v̌_i = R_i^⊤ v_i and v̌_j = R_j^⊤ v_j; then,

    [h][v̌_i][v̌_j] = ∫_[ρ∇v̌_i ·ℛ_ρ(v̌_j ) + ρ∇v̌_j ·ℛ_ρ(v̌_i) ] + ∫_[F][h]σ_h,ρv̌_i·v̌_j s.

By recalling the definition of [R][ρ] given in (<ref>), the first term on the right-hand side of (<ref>) can be written as

    ∫_ρ∇_h v̌_i ·ℛ_ρ(v̌_j )     =  - ∫_[F][h]ρ∇_h v̌_i_ω·v̌_j s,

since ρ is piecewise constant and v̌_j ∈ R_i^⊤ V_i ⊆ V_. By observing that v̌_j|_F = 0 for all F∈[F][h] such that F ∩_j = ∅, and that ρ∇_h v̌_i_ω|_F = 0 when F ∩_i = ∅, we have

    ∫_ρ∇_h v̌_i  ·ℛ_ρ(v̌_j )     =  - ∫_∂_i ∩∂_jv̌_j·ρ∇_h v̌_i_ω s 
       ≤σ_h,ρ^1/2v̌_j_L^2(∂_i ∩∂_j)^2 + σ_h,ρ^-1/2ρ∇_h v̌_i_ω_L^2(∂_i ∩∂_j)^2 
       ≤σ_h,ρ^1/2v̌_j_L^2(∂_j)^2 + σ_h,ρ^-1/2ρ∇_h v̌_i_ω_L^2(∂_i)^2.

Here, the second term on the right-hand side of (<ref>) can be bounded by applying <Ref> and noting that ρ_≤ 2 ρ_ and h_≤ 2 h_, as follows:

    σ_h,ρ^-1/2ρ∇_h v̌_i_ω_L^2(∂_i)^2
       ≤∑_⊂_iσ_h,ρ^-1/2ρ_∇_h v̌_i _L^2(∂)^2 
        = C_σ^-1∑_⊂_ih_/ρ_ p^2ρ_^2 ∇_h v̌_i _L^2(∂)^2 
       ≲∑_⊂_i√(ρ_)∇_h v̌_i) _L^2()^2
     = √(ρ)∇_h v̌_i _L^2(_i)^2.

Inserting (<ref>) into (<ref>) gives

    ∫_ρ∇v̌_i ·ℛ_ρ(v̌_j ) ≲√(ρ) ∇_h v̌_i _L^2(_i)^2 + σ_h,ρ^1/2v̌_j_L^2(∂_j)^2.

Similarly, we have that

    ∫_ρ∇v̌_j ·ℛ_ρ(v̌_i ) ≲√(ρ) ∇_h v̌_j _L^2(_j)^2 +  σ_h,ρ^1/2v̌_i_L^2(∂_i)^2

and

    ∫_[F][h]σ_h,ρv̌_i·v̌_j s ≲σ_h,ρ^1/2v̌_i_L^2(∂_i)^2 + σ_h,ρ^1/2v̌_j_L^2(∂_j)^2.

Substituting (<ref>), (<ref>) and (<ref>) into (<ref>) we obtain

    [h][v̌_i][v̌_j]    ≲√(ρ) ∇_h v̌_i _L^2(_i)^2 + √(ρ) ∇_h v̌_j _L^2(_j)^2 
          + σ_h,ρ^1/2v̌_i_L^2(∂_i)^2 + σ_h,ρ^1/2v̌_j_L^2(∂_j)^2.

The result follows by summing over i,j=1,…,N_, i  j, and exploiting (<ref>).


For the forthcoming analysis we also require an extension of the trace-inverse inequality introduced by Feng and Karakashian in <cit.>; cf. also Smears <cit.>, to which we refer for the proof.

Let [T][h] and [T][] be a pair of nested polytopic grids. We assume that [T][] is obtained by agglomeration of elements of [T][] and that both [T][] and [T][] satisfy <Ref>. Moreover, we assume that for each _i ∈[T][], i=1,…,N_, there exists an 𝐱_0,i∈_i such that (𝐱 - 𝐱_0,i) ·𝐧_i ≳ for all 𝐱∈∂_i, where 𝐧_i is the unit outward normal vector to ∂_i. Then, for any v_h ∈ V_h, writing [F][h](_i) := { F ∈[F][h] such that F ⊂_i, F ⊄∂_i }, the following bound holds:

    v_h _L^2(∂_i)^2
       ≲∇_h v_h _L^2(_i) v_h _L^2(_i)
    + 1/ v_h _L^2(_i)^2 
          + σ_h,1^1/2v_h_L^2([F][h](_i)) v_h _L^2(_i).















































§ CONDITION NUMBER ESTIMATES

In this section we derive an upper bound on the condition number of P_ad by following the analysis presented in <cit.>; see also <cit.>. To this end, we show that the following three assumptions are satisfied.


[Local stability]
There exists an α∈ (0,2) such that

    [][R_i^⊤ v_i][R_i^⊤ v_i] ≤α[i][v_i][v_i]   ∀  v_i ∈ V_i, i=0,1,…,N_H.


We point out that in our setting <Ref> immediately follows with α = 1 from the definition of [i] given in <Ref>, cf. <cit.>.


[Strengthened Cauchy–Schwarz inequality]
There exist constants ϵ_ij∈ [0,1], for 1 ≤ i,j ≤ N_, such that

    | [h][R_i^⊤ v_i][R_j^⊤ v_j] | ≤ϵ_ij[h][R_i^⊤ v_i][R_i^⊤ v_i]^1/2[h][R_j^⊤ v_j][R_j^⊤ v_j]^1/2

for all v_i ∈ V_i, v_j ∈ V_j. Define Θ([𝐄]) to be the spectral radius of
([𝐄])_ij =  {ϵ_ij}_i,j=1,…,N_.

<Ref> immediately follows since each subdomain _i ∈[T][], i=1,…,N_, can possess only a finite number of neighbours, cf. (<ref>). In particular, by observing that if ∂_i ∩∂_j = ∅, then [h][R_i^⊤ v_i][R_j^⊤ v_j] = 0 for all v_i ∈ V_i,  v_j ∈ V_j, we deduce that
ϵ_ij = 0 if ∂_i ∩∂_j = ∅, ϵ_ij = 1, otherwise. Then Θ([𝐄]) is uniformly bounded by (N_𝕊+1), where N_𝕊 is the maximum number of neighbours that each subdomain may possess, cf. (<ref>). This result ensures that a stable (in the sense of the energy norm) decomposition can be found for the local spaces and the coarse one.


[Stable decomposition]
Each v_h ∈ V_h admits a decomposition of the form v_h = ∑_i=0^N_ R_i^⊤ v_i, v_i ∈ V_i, i=1,…,N_, and v_0 ∈ V_0, such that

    ∑_i=0^N_[i][u_i][u_i] ≤ C_♯^2 [h][u_h][u_h].



Following <cit.> the upper bound on the condition number of P_ad is stated in the following theorem.

Supposing that <Ref>–<Ref> hold, the condition number K(P_ad) of the additive Schwarz operator P_ad is bounded as follows:

    K(P_ad) ≲ C_♯^2 α (Θ([𝐄]) + 1),

where α, [𝐄], and C_♯ are as defined in Assumptions <ref>, <ref> and <ref>, respectively.


Next we prove that <Ref> holds.

Let v_h ∈ V_h, and assume that the grid [T][] satisfies Assumptions <ref> and <ref>. We also assume that [T][] is obtained by agglomeration of elements of [T][], [T][] is obtained by agglomeration of elements of [T][], and that both [T][] and [T][] satisfy Assumptions <ref> and <ref>. Then, <Ref> holds with

    C_♯^2 [ max_j=1,…,N_( ρ_j/ρ_j) ] (p^2/q/ + p^2/q^2^2/),

where ρ_j = min_∈[D][j] (ρ()) and ρ_j = max_∈[D][j] (ρ()).



Given v_h ∈ V_h, we select v_0 = v_, where v_∈ V_ is defined as in the proof of <Ref>. Then, by employing <Ref>, v_h - R_0^⊤ v_0 can be uniquely decomposed as
v_h - R_0^⊤ v_0 = ∑_i=1^N_ R_i^⊤ v_i,
where v_i = R_i (v_h - R_0^⊤ v_0), i=1,…,N_, and

    [h][v_h - R_0^⊤ v_0][ v_h - R_0^⊤ v_0] = ∑_i=1^N_[i][v_i][v_i] + ∑_i,j=1, i j^N_[h][R_i^⊤ v_i][R_j^⊤ v_j].

Adding [0][v_0][v_0] to both sides of (<ref>) we obtain the following inequality:





    | ∑_i=0^N_[i][v_i][v_i] |    ≤| [h][v_h - R_0^⊤ v_0][v_h - R_0^⊤ v_0] |   +    | [0][v_0][v_0]  | 
          + | ∑_i,j=1, i j^N_[h][R_i^⊤ v_i][R_j^⊤ v_j] | 
       ≡1 + 2 + 3.

From the definition of [0], cf. (<ref>), we have that





    2   ≤ |[h][R_0^⊤ v_0 - v_h][R_0^⊤ v_0]| + |[h][v_h][R_0^⊤ v_0]| 
       ≤ |[h][R_0^⊤ v_0 - v_h][R_0^⊤ v_0 - v_h]| + 2 |[h][R_0^⊤ v_0 - v_h][v_h]| + |[h][v_h][v_h]|.

Recalling the continuity of [A][h], cf. <Ref>, and applying the triangle inequality and Young's inequality gives

    |[h][R_0^⊤ v_0 - v_h][v_h]|    ≲[v_h - R_0^⊤ v_0][h,ρ] [v_h][h,ρ]
     ≲[v_h - R_0^⊤ v_0][h,ρ]^2 + [v_h][h,ρ]^2.

Then, by inserting the above bound into (<ref>) and using the continuity and coercivity of [A][h], cf. <Ref>, we deduce that




    1 + 2≲[v_h - R_0^⊤ v_0][h,ρ]^2 + [h][v_h][v_h].

In particular, we observe that from the definition of ·_h,ρ we have


    [v_h - R_0^⊤ v_0][h,ρ]^2  = √(ρ)∇_h (v_h - R_0^⊤v_0) _L^2([T][h])^2 + σ_h,ρ^1/2v_h - R_0^⊤ v_0_L^2([F][h])^2.


Writing [F][] to denote the set of faces of [T][], and observing that [F][] ⊆[F][] since [T][] ⊆[T][], the second term on the right-hand side
of (<ref>) can be bounded as follows:

    σ_h,ρ^1/2v_h  - R_0^⊤ v_0_L^2([F][h])^2
        = σ_h,ρ^1/2v_h - R_0^⊤ v_0_L^2([F][h]∖[F][])^2
     + σ_h,ρ^1/2v_h - R_0^⊤ v_0_L^2([F][])^2 
        =  σ_h,ρ^1/2v_h_L^2([F][h]∖[F][])^2
    + σ_h,ρ^1/2v_h - R_0^⊤ v_0_L^2([F][])^2
       ≤[v_h][h,ρ]^2 + ∑_j=1^N_σ_h,ρ^1/2 (v_h - R_0^⊤ v_0) _L^2(∂[D][j])^2 
       ≤[v_h][h,ρ]^2 + ∑_i=1^N_σ_h,ρ^1/2 (v_h - R_0^⊤ v_0) _L^2(∂_i)^2,


where we have used that R_0^⊤ v_0 = 0 on each face F ∈[F][h] ∖[F][] and, in the last step, the fact that [T][] ⊆[T][]; cf. <Ref>.
Hence, inserting (<ref>) into (<ref>) and employing <Ref>, inequality (<ref>) becomes





    1 + 2≲√(ρ) ∇_h(v_h - R_0^⊤ v_0)_L^2([T][h])^2
     + ∑_i=1^N_σ_h,ρ^1/2 (v_h - R_0^⊤ v_0) _L^2(∂_i)^2 + | [h][v_h][v_h] |.

From <Ref> we get




    3≲√(ρ) ∇_h(v_h - R_0^⊤ v_0)_L^2([T][h])^2 + ∑_i=1^N_σ_h,ρ^1/2 (v_h - R_0^⊤ v_0)_L^2(∂_i)^2.

Thereby, (<ref>) can be bounded as follows





    | ∑_i=0^N_[i][v_i][v_i] |    ≲ | [h][v_h][v_h] | + √(ρ) ∇_h(v_h - R_0^⊤ v_0) _L^2([T][h])^2 
          + ∑_i=1^N_σ_h,ρ^1/2 (v_h - R_0^⊤ v_0)_L^2(∂_i)^2 
       ≡4 + 5 + 6.

Thanks to <Ref> we have that







    5    = ∑_j=1^N_√(ρ)∇_h(v_h - R_0^⊤ v_0 )_L^2([D][j])^2 ≲∑_j=1^N_ρ_j ∇_h(v_h - R_0^⊤ v_0 )_L^2([D][j])^2 
       ≲∑_j=1^N_ρ_j  [∇_h v_h _L^2([T][h,j])^2 + σ_h,1^1/2v_h_L^2([F][h,j]^I)^2 ] 
    
    
    
       ≲∑_j=1^N_ρ_j/ρ_j[√(ρ) ∇_h v_h _L^2([T][h,j])^2 + σ_h,ρ^1/2v_h_L^2([F][h,j]^I)^2 ]
       ≲max_j=1,…,N_( ρ_j/ρ_j)  [v_h][h,ρ]^2 ≲max_j=1,…,N_( ρ_j/ρ_j) [h][v_h][v_h],

where we have used the coercivity bound from <Ref> in the last inequality. The bound on term 6 can be deduced by using the inverse trace inequality of <Ref>. To this end, we first observe that




    6≲∑_i=1^N_p^2 max_{⊂_i}ρ_/h v_h - R_0^⊤ v_0 _L^2(∂_i )^2,

where we have also employed the definition of σ_h,ρ and the fact that ρ_|_F≤ 2 ρ_^± for any F ⊂∂_i, F ⊂∂^±, for some ^±∈[T][], which implies that ρ_|_F ≤ 2 max_{⊂_i}ρ_ for all F ∈[F][] such that F ⊂∂_i. Then, by applying <Ref> to each _i ∈[T][], i=1,…,N_, from (<ref>) we obtain the following bound:






    6≲∑_i=1^N_   p^2 max_{⊂_i}ρ_/h[  ∇_h (v_h - R_0^⊤ v_0) _L^2(_i) v_h - R_0^⊤ v_0 _L^2(_i)
        + 1/  v_h - R_0^⊤ v_0 _L^2(_i)^2
        + ( ∑_F ∈[F][h](_i)σ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2(_i)].

Since [T][] ⊆[T][], we denote by [I][j] := { k: 1≤ k ≤ N_, _k ∈[T][] and _k ⊂[D][j]} the set of indices that correspond to the subdomains inside [D][j] ∈[T][], for all j=1,…,N_. Hence, [I][j] ∩[I][k] = ∅ for any j  k, 1≤ j,k ≤ N_, and ∪_j=1^N_[I][j] = {1,…,N_}. Then,


    6   ≲∑_j=1^N_∑_i ∈[I][j]p^2 max_{⊂_i}ρ_/h[  ∇_h (v_h - R_0^⊤ v_0) _L^2(_i) v_h - R_0^⊤ v_0 _L^2(_i)
            + 1/  v_h - R_0^⊤ v_0 _L^2(_i)^2
            + ( ∑_F ∈[F][h](_i)σ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2(_i)] 
    
    
    
       ≲∑_j=1^N_p^2 ρ_j/h[ ∑_i ∈[I][j]∇_h (v_h - R_0^⊤ v_0) _L^2(_i) v_h - R_0^⊤ v_0 _L^2(_i)
            + 1/ ∑_i ∈[I][j] v_h - R_0^⊤ v_0 _L^2(_i)^2 
            + ∑_i ∈[I][j]( ∑_F ∈[F][h](_i)σ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2(_i)].

We now proceed by bounding each term present in the bracket in (<ref>); to this end, using the Cauchy–Schwarz inequality for sums, we get



    ∑_i ∈[I][j]∇_h (v_h -     R_0^⊤ v_0) _L^2(_i) v_h - R_0^⊤ v_0 _L^2(_i)
       ≤( ∑_i ∈[I][j]∇_h (v_h - R_0^⊤ v_0) _L^2(_i)^2 )^1/2( ∑_i ∈[I][j] v_h - R_0^⊤ v_0 _L^2(_i)^2 )^1/2
    
        = | v_h - R_0^⊤ v_0 |_H^1([T][h,j]) v_h - R_0^⊤ v_0 _L^2([D][j]).

Similarly, by noting that [F][h](_i) is the set of faces F ∈[F][] strictly contained in _i, and therefore ∪_i ∈[I][j][F][h](_i) ⊂[F][h,j]^I, we deduce that


    ∑_i ∈[I][j](  ∑_F ∈[F][h](_i)   σ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2(_i)
       ≤( ∑_i ∈[I][j]σ_h,1^1/2v_h - R_0^⊤ v_0_L^2([F][h](_i))^2 )^1/2( ∑_i ∈[I][j] v_h - R_0^⊤ v_0 _L^2(_i)^2 )^1/2
       ≤( ∑_F ∈[F][h,j]^Iσ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2([D][j]).


Noting that
∑_i ∈[I][j] v_h - R_0^⊤ v_0 _L^2(_i)^2 =  v_h - R_0^⊤ v_0 _L^2([D][j])^2 gives

    6≲∑_j=1^N_   p^2 ρ_j/h[ |v_h - R_0^⊤ v_0|_H^1([T][h,j]) v_h - R_0^⊤ v_0 _L^2([D][j])
     + 1/  v_h - R_0^⊤ v_0 _L^2([D][j])^2
        + ( ∑_F ∈[F][h,j]^Iσ_h,1^1/2v_h - R_0^⊤ v_0_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2([D][j])].


The last term on the right-hand side of (<ref>) can rewritten as

    ( ∑_F ∈[F][h,j]^Iσ_h,1^1/2v_h - R_0^⊤ v_0   _L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2([D][j])
        = ( ∑_F ∈[F][h,j]^Iσ_h,1^1/2v_h_L^2(F)^2 )^1/2 v_h - R_0^⊤ v_0 _L^2([D][j])
        = σ_h,1^1/2v_h_L^2([F][h,j]^I) v_h - R_0^⊤ v_0 _L^2([D][j]);

 here we observe that R_0^⊤ v_0|_F = 0 on each F ∈[F][h,j]^I, since [T][] and [T][] are nested. 
Then, by employing the above estimate together with  <Ref>, we deduce that






    6   ≲∑_j=1^N_[ p^2 ρ_j/h( /q + 1/q^2^2/) (∇_h v_h _L^2([T][h,j])^2 + σ_h,1^1/2v_h_L^2([F][h,j]^I)^2 ) ] 
       ≲∑_j=1^N_ρ_j/ρ_j( p^2/q/h + p^2/q^2^2/h ) (√(ρ) ∇_h v_h _L^2([T][h,j])^2 + σ_h,ρ^1/2v_h_L^2([F][h,j]^I)^2) 
       ≲max_j=1,…,N_( ρ_j/ρ_j) (p^2/q/h + p^2/q^2^2/h ) [h][v_h][v_h],

where we have also made use of the coercivity bound of <Ref> in the last inequality. Inserting the estimates (<ref>) and (<ref>) into (<ref>) we obtain the desired result.












According to the statement of <Ref>, given that Assumptions <ref> and <ref> hold, using <Ref> we deduce that

    K(P_ad) ≲max_1≤ j ≤ N_( ρ_j/ρ_j) (p^2/q/h + p^2/q^2^2/h )(N_𝕊+1).

In particular, in the lowest order case, i.e., when p=q=1, we have
K_h(P_ad) ≲^2/h,
which is in agreement with the corresponding bound derived in <cit.>. On the other hand if the size of the coarse subdomain and fine meshes are fixed, we deduce that
K_p(P_ad) ≲p^2/q.
Moreover, we also observe that if the diffusion coefficient ρ is constant on each subdomain [D][j], j=1,…,N_, then the condition number is independent of the jump in ρ.




We remark that <Ref> is needed in order to obtain the local estimates of <Ref>, which allow to bound K(P_ad) with max_1≤ j ≤ N_ ( ρ_j/ρ_j), cf. (<ref>). However, if the diffusion coefficient ρ is constant on Ω, we point out that the analysis can be simplified by employing the global estimates of <Ref> without making <Ref>. We refer to <cit.> for further details.






























§ NUMERICAL RESULTS


In this section, we present a series of numerical experiments to demonstrate the sharpness of the condition number bounds stated in Remarks <ref> and <ref>. Throughout this section we solve (<ref>) by using the additive Schwarz Preconditioned Conjugate Gradient (ASPCG) method; here, we report the number of iterations needed to reduce the Euclidean norm of the relative residual vector below a tolerance of 10^-8, based on starting from the trivial initial guess. Furthermore, we estimate the condition number K(P_ad) by using the extreme-eigenvalue-estimate based on the ASPCG iterations. For the sake of simplicity of the presentation, we only consider the massively parallel case, i.e., when [T][] = [T][]. Furthermore, here we select the penalty parameter C_σ = 10.






































































































































































 §.§ Example 1








































































In this first example we investigate the dependence of the condition number of P_ad on the diffusion coefficient ρ. Based on <Ref>, we expect the condition number of the preconditioned system to be dependent on the choice of the coarse grid [T][]. More precisely, if [T][] is chosen to be aligned with the discontinuities of ρ, then K(P_ad) should be independent of the jump in the coefficient ρ; otherwise it will depend on the maximum ratio between the maximum and the minimum value of ρ present inside the subdomains [D][j] ∈[T][]. To verify this behavior, we consider two experiments based on fine/coarse grids [T][h]/[T][], respectively, where [T][] is a Voronoi polygonal grid on the L-shaped domain  depicted in Figure <ref> with 16 polygonal elements and [T][] is obtained by successive refinement of elements of [T][]. Here, we observe that the elements of [T][] are convex and satisfy <Ref>, cf. <Ref> (left). Moreover, we choose the polynomial degrees to be either p=q=1 or p=q=2. In the first experiment we fix ρ|_[D][j]  = ρ_o = 1 on the elements [D][j]∈[T][] with odd index j and set ρ|_[D][j]  = ρ_e∈{10^1,10^2,…,10^6}, in each test case, on the polygonal subdomains with even index j. The results shown in the first two lines of <Ref> confirm the independence with respect to the jumps of ρ when those jumps are aligned with the subdomains of [T][]. In the second experiment we proceed similarly, but here we take different values of ρ on odd and even polygonal elements ∈[T][h]: in this way [T][] is not aligned with the discontinuities of ρ, and hence the ratio between the maximum and the minimum value of ρ inside the polygonal subdomains [D][j] ∈[T][] is given by ρ_e. As expected from the theory, the results presented in the last two lines of <Ref> show that the condition number of P_ad grows linearly with ρ_e.
Finally, we repeat the same set of experiments by first selecting [T][] as a Voronoi polygonal grid consisting of 2000 elements and, subsequently, define the coarse grid [T][] by successive agglomerations of elements of [T][]. The agglomeration is undertaken based on employing Metis, cf. <cit.>. As shown in <Ref> (right), the elements of [T][] are clearly non-convex in this example. Although <Ref> is not satisfied in this case, the results of <Ref> illustrate analogous behavior to that observed in the previous setting when the coarse elements were convex.






















 §.§ Example 2




In this example, we investigate the performance of the proposed ASPCG algorithm on a set of
Voronoi polygonal fine grids [T][h], where =(0,1)^2 and ρ=1. For each grid [T][h] we construct a sequence of nested polygonal grids [T][] obtained by successive levels of agglomeration, cf. <cit.>. 
For each fine Voronoi grid of size h the agglomeration process has been performed in order to ensure that the size of the coarser partitions is approximately = 2h, 4h, …, cf. <Ref>, for example. In <Ref> we report the condition number and the iteration counts for the proposed ASPCG algorithm for p=q=1 and p=q=3, respectively. Here, we clearly observe that the condition number and the iteration counts grow quadratically and linearly, respectively, as  tends to zero for fixed . Moreover, if the ratio of  and  is kept fixed, then we observe that the condition number and iteration counts are approximately constant; cf. the diagonals and subdiagonals of <Ref>. This behavior is in agreement with the theoretical bound stated in <Ref>, where K(P_ad) = [O]( ^2/h^2 ).
















1.8












































































































 §.§ Example 3

We now consider the performance of the ASPCG algorithm on tetrahedral meshes in three dimensions. To this end, we set = (0,1)^3 and ρ=1; furthermore, the elements of the coarse mesh are general-shaped polyhedra obtained by successive agglomeration, cf. the previous example. The results for p=q=1 and p=q=3 are reported in <Ref>, respectively. Here, we have also added a line with the condition number of the operator A_:V_× V_→ V_ defined by (A_ u_, v_)_L^2() = [][u_][v_] for all u_,v_∈ V_, and, in parentheses, the iteration counts of the Conjugate Gradient method for solving (<ref>) without preconditioning. Analogous behaviour of the condition number and iteration counts to those presented in the previous example are observed. In particular, we observe that the condition number is roughly constant on the diagonals and subdiagonals of the two tables, while, along each row, i.e., when [T][] is fixed, the expected quadratic growth in K(P_ad) is observed. Similar considerations are also noted for the iteration counts.


1.5
 



































 §.§ Example 4




Given the definition of R_0^⊤, the proposed ASPCG algorithm naturally admits the use of non-nested coarse spaces, i.e., when V_⊈ V_h. In order to confirm the condition number bound stated in <Ref> when V_⊈ V_h, we consider a set of independently generated Voronoi polygonal tessellations of (0,1)^2 of size  and >h, respectively; in this way, [T][h] and [T][] are non-nested, cf. <Ref>.
The results (for ρ=1) shown in <Ref> for p=q=1 and p=q=3, respectively, illustrate analogous behavior to the results for the nested case presented in the previous examples; this is in agreement with the condition number bound stated in <Ref>.


1.5





































 §.§ Example 5

In this final example we investigate the dependence of the condition number on the polynomial degree p in both the nested and non-nested cases with ρ=1. For the nested case, we consider a total of four tests: two of them are characterized by quadrilateral fine grids with N_h = 256 and N_h = 1024 elements, while the two other tests are based on employing the polygonal fine grids depicted in <Ref>, where the fine meshes have N_h = 262 and N_h = 516 polygonal elements. For each test the coarse mesh [T][] is obtained by agglomeration of [T][] in order to guarantee h/4. Analogous fine meshes are also considered in the non-nested setting; however, here the coarse mesh [T][] is selected to be a Voronoi grid generated independently of [T][h], cf. <Ref>. In <Ref> we plot the condition number K(P_ad) on each set of grids as the polynomial degree p is increased. Here, we observe that, in the nested setting, i.e., when V_⊆ V_h, for a fixed mesh size K(P_ad) = [O](p) as p increases; however, when  V_⊈ V_h, then K(P_ad) = [O](p^2) as p increases. This behaviour is in agreement with the condition number bounds stated in <Ref>, respectively.








































































































































































































§ PROOF OF THEOREM <REF>
 
In this section we present the proof of Theorem <ref>. To this end, suppose that Ω⊂ℝ^n is a bounded, open, convex domain with boundary ∂Ω. Given f ∈ L_0^2(Ω), consider the homogeneous Neumann problem


    - Δ u  = f    in Ω,     ∇ u ·𝐧 = 0    on ∂Ω.


The weak formulation of (<ref>) is: find u ∈ H^1(Ω)/ℝ := H^1(Ω)∩ L^2_0(Ω), such that


    a(u,v):=∫_Ω∇ u ·∇ v   = ℓ(v) := ∫_Ω f v      ∀  v ∈ H^1(Ω)/ℝ,


with H^1(Ω)/ℝ equipped with the norm v_H^1(Ω)/ℝ := ∇ v_L^2(Ω). From <cit.> it follows that u∈ H^2(Ω); the proof of (<ref>) now proceeds with the following steps.

Step 1. [H^1(Ω) bound] The existence of a unique weak solution to (<ref>) follows by the Lax–Milgram lemma applied to the bilinear form a(·,·), which is bounded and coercive on H^1(Ω)/ℝ, and noting that the linear functional ℓ(·) is bounded on H^1(Ω)/ℝ. Indeed,


    ℓ(v) = ∫_Ω f v   = ∫_Ω f [v - v_Ω] ≤f_L^2(Ω) v - v_Ω_L^2(Ω),

where v_Ω := 1/|Ω|∫_Ω v.

By Poincaré's inequality,


    v - v_Ω_L^2(Ω)≤ C(Ω) ∇ v_L^2(Ω),


where C(Ω) is a positive constant. Recalling that Ω is a bounded, open, convex domain, it can be shown that C(Ω) ≤1/π (Ω), cf. <cit.>. Setting v=u in the weak formulation above, we then deduce that


    ∇ u^2_L^2(Ω) = ℓ(u) ≤ C(Ω) f_L^2(Ω)∇ u_L^2(Ω)≤1/π(Ω)f_L^2(Ω)∇ u_L^2(Ω).


Hence,

    ∇ u_L^2(Ω)≤1/π(Ω)f_L^2(Ω).




Step 2. Suppose that Ω⊂Ω_m, where Ω_m, m=1,2,…, is a sequence of bounded, open, convex C^2 domains, such that (Ω,Ω_m) ≤1/m.
The existence of such a sequence {Ω_m}_m=1^∞ follows from Eggleston's lemma, cf. <cit.>. Hence, in particular 0 ≤(Ω_m) - (Ω) → 0 as m →∞.



Step 3. With Ω_m, m=1,2,…, as in Step 2, we extend f by 0 from Ω to Ω_m for each m=1,2,…, and define


    f_m():= {[ f()   ,;   0     ].


Clearly, f_m ∈ L^2_0(Ω_m). We consider the following Neumann problem on Ω_m:


    - Δ u_m  = f_m    in Ω_m,     ∇ u_m ·𝐧_m = 0    on ∂Ω_m,

where 𝐧_m is the unit outward normal vector to ∂Ω_m. We have that the unique weak solution u_m ∈ H^1(Ω_m)/ℝ of the above Neumann problem satisfies the following elliptic regularity result: u_m ∈ H^2(Ω_m), cf. <cit.>, and


    u_m_H^2(Ω_m)≤√(6) ( f_m _L^2(Ω_m) +  u_m _L^2(Ω_m)),     m=1,2,…,


cf. <cit.>, with λ=1. Thereby, we deduce that


    u_m_H^2(Ω_m)≤√(6) ( f_L^2(Ω) +  u_m _L^2(Ω_m)),     m=1,2,….


Since Ω⊂Ω_m, upon application of Poincaré's inequality on the right-hand side of (<ref>),
followed by recalling the H^1(Ω) bound derived in Step 1, we get


    u_m_H^2(Ω)≤u_m_H^2(Ω_m)   ≤   √(6) ( f_L^2(Ω) + 1/π(Ω_m) ∇ u_m _L^2(Ω_m))
       ≤   √(6) ( f_L^2(Ω) + 1/π^2[(Ω_m)]^2 f_m_L^2(Ω_m))
       =   √(6) (1 + 1/π^2[(Ω_m)]^2) f_L^2(Ω),


m=1,2,…. Thanks to Step 2, lim_m →∞(Ω_m) = (Ω). As ((Ω_m))_m=1^∞ is a convergent sequence in ℝ it is automatically a bounded sequence, and therefore, because of (<ref>), there exists a positive constant C_0, independent of m, such that


    u_m_H^2(Ω)≤u_m_H^2(Ω_m)≤ C_0    .


Thus, (u_m)_m=1^∞ is a bounded sequence in H^2(Ω). Hence, there exists an element u_∞∈ H^2(Ω)
and a weakly convergent subsequence u_m_k⇀ u_∞ in H^2(Ω).

By weak lower semicontinuity of the norm function ·_H^2(Ω), we have that


    u_∞_H^2(Ω)≤lim inf_k →∞u_m_k_H^2(Ω).

Further, thanks to the compact Sobolev embedding H^2(Ω) ⋐ H^1(Ω) guaranteed by the Rellich–Kon­drashov theorem, by extracting a further subsequence (not indicated), we have that u_m_k→ u_∞ strongly
in H^1(Ω). Now,


    ∇ u_m_k_L^2(Ω_m_k)^2 = ∇ u_m_k^2_L^2(Ω) + ∇ u_m_k^2_L^2(Ω_m_k∖Ω).


Focusing on the second term on the right-hand side of (<ref>), by Hölder's inequality with conjugate exponents α=p/2 and α' = α/α-1 = p/(p-2), 2<p< 2n/(n-2) (where 2n/(n-2) is the critical Sobolev index), and (<ref>), we have that

    ∇ u_m_k_L^2(Ω_m_k∖Ω)
    
    
       ≤   ∇ u_m_k_L^p(Ω_m_k∖Ω) |Ω_m_k∖Ω|^p-2/2p
       ≤   ∇ u_m_k_L^p(Ω_m_k) |Ω_m_k∖Ω|^p-2/2p
       ≤    C(|Ω_m_k|,(Ω_m_k),n,p) u_m_k_H^2(Ω_m_k)  |Ω_m_k∖Ω|^p-2/2p
       ≤    C(n,p) u_m_k_H^2(Ω_m_k)  |Ω_m_k∖Ω|^p-2/2p
       ≤    C(n,p)  C_0   |Ω_m_k∖Ω|^p-2/2p→ 0     .


Passing to the limit k →∞ in (<ref>) we therefore have that


    lim_k →∞∇ u_m_k_L^2(Ω_m_k)^2 = lim_k →∞∇ u_m_k^2_L^2(Ω) = ∇ u_∞^2_L^2(Ω),


where the last equality follows from the strong convergence u_m_k→ u_∞ in H^1(Ω).

Recalling from (<ref>) that

    u_m_k_H^2(Ω)≤√(6) ( f_L^2(Ω) + 1/π(Ω_m_k) ∇ u_m_k_L^2(Ω_m_k)),


passage to the limit k →∞ in inequality (<ref>) using (<ref>), (<ref>), together with lim_m→∞(Ω_m) = (Ω), yields that

    u_∞_H^2(Ω)≤√(6) ( f_L^2(Ω) + 1/π(Ω) ∇ u_∞_L^2(Ω)).





Step 4. [Identification of u_∞] It remains to show that u_∞ = u, the weak solution of the original Neumann problem on Ω. To this end, we consider the weak formulation of the Neumann problem satisfied by u_m:


    ∫_Ω_m∇ u_m ·∇ v   = ∫_Ω_m f_m v      ∀  v ∈ H^1(Ω_m)/ℝ.


Thanks to the definition of f_m, this weak formulation is equivalent to


    ∫_Ω_m∇ u_m ·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω_m)/ℝ,


and therefore


    ∫_Ω∇ u_m ·∇ v  
    + ∫_Ω_m ∖Ω∇ u_m ·∇ v   =
    ∫_Ω f v      ∀  v ∈ H^1(Ω_m)/ℝ,


whereby


    ∫_Ω∇ u_m_k·∇ v  
    + ∫_Ω_m_k∖Ω∇ u_m_k·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω_m_k)/ℝ.


Equivalently, because ∫_Ω f()   = 0,


    ∫_Ω∇ u_m_k·∇ v  
    + ∫_Ω_m_k∖Ω∇ u_m_k·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω_m_k).


Consider a fixed bounded domain Ω_0 ⊂ℝ^n such that Ω_0 ⋑Ω_m ⊃Ω for all
m=1,2,…. Then,


    ∫_Ω∇ u_m_k·∇ v  
    + ∫_Ω_m_k∖Ω∇ u_m_k·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω_0).


By noting (<ref>), the strong convergence u_m_k→ u_∞ in H^1(Ω) as k →∞, we have that


    ∫_Ω∇ u_∞·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω_0),


hence also


    ∫_Ω∇ u_∞·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω),


since any element of v ∈ H^1(Ω_0) can be viewed as the extension of a v ∈ H^1(Ω) to the superset Ω_0.
Therefore, again since ∫_Ω f(x) dx = 0, also


    ∫_Ω∇ u_∞·∇ v   = ∫_Ω f v      ∀  v ∈ H^1(Ω)/ℝ.


Thus we have shown that u_∞ coincides with the unique weak solution u of the homogeneous Neumann problem posed on Ω.
Returning with this information to (<ref>), we have that the weak solution of the homogeneous Neumann problem on the
bounded, open, convex (and therefore Lipschitz) domain Ω satisfies


    u_H^2(Ω)   ≤√(6) ( f_L^2(Ω) + 1/π(Ω) ∇ u _L^2(Ω)),

as required.



§ CONDITION NUMBER ESTIMATES FOR NON-NESTED GRIDS

In this Appendix we provide a bound on the condition number of the additive Schwarz operator P_ad introduced in <Ref> when the fine and coarse grids [T][] and [T][], respectively, are non-nested. For the sake of simplicity, here we assume that ρ = 1 on  and consider the massively parallel case, i.e., when [T][] = [T][]. For the purposes of the proceeding analysis, we also assume that Ω is convex; moreover, we make the following additional assumption on [T][].
 (Coverability) For every polytopic element ∈[T][h], there exists a set of m_ overlapping shape-regular simplices [K]_i, i=1,…,m_, such that

    (,∂[K]_i) ≲ diam([K]_i)/p^2,  and   |[K]_i| ≳ |  |,

for all i=1,…,m_, see <cit.>.

Given that Assumption <ref> holds, we state the following inverse inequality, cf. <cit.>.


Suppose that v_∈ V_; then, the following bound holds:

    [∇ v_h][2][][2]≲  p^4 h_^-2[v_h][2][][2]     ∀ ∈[T][h].



We refer to <cit.> for the proof of this result.


We first provide a counterpart of <Ref>, which holds in the non-nested case and allows us to prove the validity of <Ref> also for non-nested spaces V_ and V_. The key aspect of our analysis is the construction of the conforming approximant introduced in <Ref>. In particular, we recall the following result.


Let 𝒢_h(v_h) := ∇_h v_h + [R][1](v_h) be the discrete gradient operator of v_h ∈ V_h defined as in <Ref> and let v_h∈ H^1_0() be such that

    ∫_∇v_h·∇ w  = ∫_𝒢_h(v_h) ·∇ w   ∀  w ∈ H^1_0().

Then, the following approximation and stability results hold:

    [v_h - v_h][2][]    ≲h/p[σ_h,1^1/2v_h][2][[F][h]],      | v_h |_H^1()≲[v_h][h,1].




<Ref> provides global bounds for v_∈ V_ in the L^2-norm. This result is a particular case of <Ref>, where local bounds on each coarse element [D][j] ∈[T][] are provided. We refer to <cit.> for the proof of <Ref>.


On the basis of the previous result, <Ref> can be generalized to non-nested spaces as follows.

For any v_h ∈ V_h there exists a coarse function v_∈ V_ such that

    v_h - R_0^⊤ v__L^2()≲/q[v_h][h,1], 
       [v_h - R_0^⊤ v_][h,1] ≲p^2/q/[v_h][h,1].




Let v_∈ V_ and let v_∈ V_ be defined as v_ = Π_v_h, with v_h as defined in <Ref> and where Π_ is the hp-approximant introduced in <Ref>. Then, by employing the triangle inequality we have

    v_h - R_0^⊤ v__L^2()≲ v_h - v_h_L^2()+ v_h - Π^0v_h_L^2()+ Π^0v_h - R_0^⊤( Π_v_h)_L^2(),

where Π^0:L^2() → V_ is the L^2-projection operator onto V_.






From the definition of R_0^⊤, we note that
Π^0 v_ = R_0^⊤ v_ for all v_∈ V_. Hence, exploiting <Ref> together with Assumption <ref>, cf. <Ref>, gives

    v_h -  R_0^⊤ v__L^2()   ≲ v_h - v_h_L^2() + v_h - Π^0v_h_L^2()  + Π^0( v_h - Π_v_h )_L^2()
       ≤ v_h - v_h_L^2() + v_h - Π_v_h_L^2()  + v_h - Π_v_h_L^2()
       ≲ v_h - v_h_L^2()+ h/p v_h_H^1() + /q v_h_H^1();

here we have used that Π^0v _L^2()≤ v _L^2() for all v ∈ L^2() and v - Π^0 v _L^2()≤ v - w _L^2() for all w ∈ L^2(). By applying Poincaré's inequality to v_h∈ H^1_0() and noting <Ref>, inequality (<ref>) immediately follows by observing that ≤ and q ≤ p. In order to obtain (<ref>) we proceed as follows:

    [v_h - R_0^⊤ v_][h,1]^2 = ∇_h (v_h - R_0^⊤ v_) _L^2([T][h])^2 + σ_h,1^1/2 v_h - R_0^⊤ v__L^2([F][])^2.

We bound the first term on the right-hand side of (<ref>) by means of <Ref> and (<ref>) as follows:

    ∇_h (v_h - R_0^⊤ v_) _L^2([T][h])^2      = ∑_∈[T][]∇_h (v_h - R_0^⊤ v_) _L^2()^2 
       ≲p^4/^2 v_h - R_0^⊤ v__L^2()^2 ≲p^4/q^2^2/^2[v_h][h,1]^2.

The second term on the right-hand side of (<ref>) can be bounded by recalling the definition of σ_h,1, <Ref> and (<ref>) as follows:

    σ_h,1^1/2 v_h - R_0^⊤ v__L^2([F][])^2    ≲p^2/h∑_∈[T][] v_h - R_0^⊤ v__L^2(∂)^2 
       ≲p^2/hp^2/h v_h - R_0^⊤ v__L^2()^2  ≲p^4/q^2^2/^2[v_h][h,1]^2.

Inserting (<ref>) and (<ref>) into (<ref>) we obtain (<ref>).


With <Ref> in hand, we can prove the following result.

 <Ref> holds with

    C_♯^2 ( p^4/q^2^2/^2).




Let v_∈ V_. Proceeding as in the proof of <Ref>, by selecting v_0 = v_ as in <Ref>, v_ can be decomposed as v_h = ∑_i=0^N_ R_i^⊤ v_i, with v_i = R_i (v_h - R_0^⊤ v_0) ∈ V_i, i=1,…,N_, so that

    | ∑_i=0^N_[i][v_i][v_i] |    ≲[v_h - R_0^⊤ v_0][h,1]^2 + [h][v_h][v_h],

where we have used (<ref>), (<ref>) and (<ref>) with the hypothesis [T][] = [T][]. The result then immediately follows by noting (<ref>) together with the coercivity of [A][h], cf. <Ref>.



Based on <Ref>, for non-nested coarse and fine spaces V_ and V_, respectively, the condition number K(P_ad) of the additive Schwarz operator P_ad can be bounded as follows:

    K(P_ad) ≲( p^4/q^2^2/^2) (N_𝕊+1).







abbrv

