{'fraction_non_alphanumeric': 0.04190890459659542, 'fraction_numerical': 0.01505709147183069, 'mean_word_length': 4.2851458885941645, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Petrov-Galerkin multilayer discretization to second order elliptic boundary value problems

We study in this paper a multilayer discretization of second order elliptic problems, aimed at providing reliable multilayer discretizations of shallow fluid flow problems with diffusive effects. This discretization is based upon the formulation by transposition of the equations. It is a Petrov-Galerkin discretization in which the trial functions are piecewise constant per horizontal layers, while the trial functions are continuous piecewise linear, on a vertically shifted grid. We prove the well posedness and optimal error order estimates for this discretization in natural norms, based upon specific inf-sup conditions. We present some numerical tests with parallel computing of the solution based upon the multilayer structure of the discretization, for academic problems with smooth solutions, with results in full agreement with the theory developed.


Introduction and motivation
This paper deals with the numerical approximation of the Poisson and related problems by means of layer-wise discontinuous solutions. It is motivated by the construction of multilayer discretizations of fluid flow equations on shallow domains. This produces a reduction of the dimensionality of the approximated problem, yielding in practice a domain decomposition discretization that is solved by parallel procedures. This is the case for instance of the works of Fernndez-Nieto and co-workers [1,2] and of Sainte-Marie and co-workers [3] to solve the Navier-Stokes equations. In these papers the velocity is discretized by layer-wise smooth functions, discontinuous across the layer boundaries. The conservation of momentum needs the continuity of diffusive fluxes across layer boundaries. This is ensured in these papers by specific techniques based upon finite-difference discretizations in the vertical direction. Both papers perform stability analysis for the discretizations introduced, however no convergence proofs are reported. This mainly occurs because the discretizations are purely finite-difference approximations, that are not linked to variational formulation of the targeted equations.
In this paper we look for stable multilayer discretizations for which convergence proofs are reachable. More concretely, we look for stable discretizations that are related to variational formulation of the equations. A possible discretization meeting these criteria could be provided by the Discontinuous Galerkin (DG) method. This method was introduced in the 90's by Cockburn and Shu for hyperbolic conservation laws (cf. [4]), and later adapted to elliptic problems, by re-writing as two coupled first-order conservation laws for the unknown and its gradient (cf. [5,6] and references therein). DG methods for elliptic problems are based upon the flux formulation of the elliptic problem, that includes both the original unknown so as its gradient as unknowns of the discretization. This allows to approximate both the unknown and its gradient by discontinuous finite elements. The stability of the formulation strongly depends on the choice of the numerical fluxes of both unknowns.
An alternative discretization could be provided by the Discontinuous Petrov Galerkin methodology, introduced in a series of papers by Sacco and co-workers (cf. e.g. [7,8]) and systematically studied by Demkowicz and Gopalakrishnan (see [9,10,11,12] and references therein). The DPG methodology approximates the unknown by either continuous or discontinuous trial functions, while the tests functions necessarily need to be discontinuous. This formulation admits three equivalent interpretations: a Petrov-Galerkin method with test functions that achieve the supremum in the inf-sup condition, a minimum-residual method with residual measured in a suitable dual norm, and a mixed formulation where one solves simultaneously for the Riesz representation of the residual. The stability of the formulation is a direct consequence of the first of these three interpretations.
Here we propose a discretization specifically adapted to multilayer discretizations, for elliptic problems related to the Poisson equations. The main idea is to start from the solution by transposition of the equations, much as the ultraweak formulation considered in the works [7,8] by Sacco and co-workers. The solution by transposition naturally belongs to L 2 spaces, and thus admits piecewise discontinuous approximations. Based upon this procedure, we propose a Petrov-Galerkin discretization for cylindrical domains, in which the solution is a layer-wise constant function, while the test functions are continuous piecewise affine polynomials in the vertical direction, with knots in a shifted grid. We derive a single "recepy" to build this kind of discretizations: Approximate the vertical derivative of the unknown by the vertical derivative of its interpolate on the test functions space.
We prove the well-posedness of this discretization, given by an inf-sup condition satisfied by the bilinear form appearing in the discrete problem. Also, we prove optimal order error estimates for smooth solutions. We further extend the multilayer discretization technique introduced to domains which are vertical deformations of cylindrical domains, which appear in shallow water problems when a flat surface is deformed. We also prove well-posedness and optimal order error estimates for this extension. We further extend the method to Neumann boundary conditions, by slightly changing the test space.
We finally present some 3D numerical tests by parallel solution of the resulting linear system, taking advantage of the multilayer structure of the discretization, what allows to only solve 2D linear systems. These tests present speeds-ups rates ranging from 20 to 50, while presenting optimal convergence orders for smooth functions, in agreement with the theoretical expectations.
The paper is organized as follows. Section 2 introduces the motivation and the basic multilayer discretization that we consider, for the Poisson problem in a cylindric domain. This discretization is studied in Section 3, where stability, convergence and obtention of error estimates are analyzed. In Section 4 the discretization is extended to Poisson problems in domains with non-flat upper boundary, also analyzing stability, convergence and error estimates. Section 5 is devoted to the extension of the discretization to Neumann boundary conditions. Finally Section 6 present several 3D numerical tests for each of the cases considered, for smooth solutions.

Multilayer approximation
Let us consider a cylindrical domain Ω = ω × (0, L) where ω ⊂ R d is bounded domain, and d ≥ 1 is an integer number. Let us consider the homogeneous Dirichlet Poisson problem as a model problem: As is standard, this problem can be written in variational form: Given f ∈ with a (v, ϕ) = Ω ∇v ∇ϕ dΩ, and L (ϕ) = < f, ϕ > Ω .
We assume ∂Ω = Γ b ∪ Γ s ∪ Γ l , where Γ b and Γ s are the boundaries of the bottom and the top, respectively, and Γ l is the vertical boundary of the domain. We shall denote by z b = z 1 2 = 0 and z s = z N + 1 2 = L the equations of the bottom and the top interfaces Γ b and Γ s , respectively. Notice that ∂Ω α = Then the vertical boundary is split as Γ l = N α=1 Γ l,α . We shall assume from now on that ω is polygonal. Let us now introduce discretization space for the unknown v, where κ α is the characteristic function of the interval (z α−1/2 , z α+1/2 ) and V 0,k (ω) is a finite element sub-space of H 1 0 (ω). The sub-index k denotes the horizontal grid size. We assume that V 0,k (ω) = V k (ω) ∩ H 1 0 (ω) where V k (ω) is a Lagrange finite element space, constructed when d = 2 with either triangular or quadrilateral elements, i. e., there exists a grid M k such that for some integer l ≥ 1; where P l (K) is the space of polynomials on K of global degree less than or equal to l, and Q l (K) is the space of polynomials on K of degree less than or equal to l in each variable. When d = 1, the triangulation is formed by segments K, and R l (K) = P l (K).
We approximate the solution v of problem (4) by some v h ∈ X h . Observe that for all ϕ ∈ H 2 (Ω) ∩ H 1 0 (Ω), As v α h | ∂ω = 0 , then it holds where ∂ z denotes the vertical derivative. Thus, where a h is a the bilinear form defined by Next, we introduce the discrete test functions space, where σ 1 , σ 2 , · · · , σ N are piecewise affine 1D functions on the intervals [z b , z 1 ], , for all i, j = 1, 2, ..., N , as we see in Figure 2. Then Y h is a subspace of H 1 0 (Ω) formed by piecewise P 1 polynomials in the variable z with coefficients in H 1 0 (ω). Note that the intervals [z b , z 1 ], and [z N , z s ] have length h/2. Introducing them We now are in a position to state the multilayer discretization of (2) that we study in this paper, Notice that although the space Y h is not a sub-set of the form a h can be written on Y h as where we denote ∇ = (∇ H , ∂ z ) with ∇ H = (∂ 1 , · · · , ∂ d−1 ) the horizontal gradient.

Stability and convergence analysis
In this section we study the well-posedness of problem (8) in the sense of Hadamard: We prove that it admits a unique solution that depends continuously on the data. This will be based upon the Banach-Neĉas-Babuska theorem (Cf. [14]).
We consider the space X h endowed with the following H 1 0 (Ω) discrete norm: and the space Y h endowed with the H 1 0 (Ω) norm. It holds the Lemma 1. The mapping T h is an isomorphism between normed spaces. Moreover, it holds Proof: By construction T h is linear and bijective. To prove the equivalence of Also, as it follows To obtain the first inequality in (11), (13) and Young's inequality yield Analogously, to obtain the first inequality in (11), again (13) and Young's in- Observe that this results justifies the choice of the discrete norm (10) if X h is considered as an external approximation of H 1 0 (Ω). The stability of the multilayer problem (8) is stated next. Theorem 1. The form a h (·, ·) satisfies the following properties: 2. Proof: Consequently, we can express a h (v h , ϕ h ), as To obtain (14), By Young's inequality, Then, from (18), it follows 2. Let v h ∈ X h and ϕ h ∈ Y h , starting from (17) and using (11) and Holder's inequality, we have As a conclusion we deduce Proof: It is a direct consequence of estimate (14).
This implies the well posedness of problem (8): Proof: The form a h is stable by Theorem 1, and satisfies the inf-sup condition then necessarily ϕ h = 0. Consequently the hypotheses of the Banach-Necas-Babuska Theorem hold (Cf. [14]). This ensures that problem (8) admits a unique solution that depends continuously on the data f . The constant 2 in estimate (21) also follows from (14).

Convergence analysis
In this section we prove a convergence result for general solutions of problem (2), so as optimal order error estimates for smooth solutions. For that, let us construct an interpolation operator on X h for functions defined on Ω. Let P h be the prismatic grid of Ω constructed by vertical displacements of the grid T k , located at the nodes z b , z 1 , z 2 · · · , z N −1 , z N , z s . The geometric elements of P h are constructed as If {x i } i∈I is the set of Lagrange interpolation nodes for V 0,k (ω), for some set of indices I, then the set of Lagrange interpolation nodes for Y h is Following Bernardi et al. [15], we introduce a nodal Lagrange interpolation operator Π k : where the α i,j are the nodal Lagrange interpolation basis functions on V 0,k (ω) corresponding interpolation nodes r i,j , and thew i,j are averaged value of w on an element to which the node r i,j belongs.
The operator Π k satisfies the standard finite element interpolation error estimates, We shall assume that the operator Π k also is defined componentwise for vector functions, without change of notations. We next define an interpolation operator We at first obtain the optimal order error estimate Theorem 4. Assume that the weak solution of (2) satisfies v ∈ H l+1 (Ω) ∩ for some constant C > 0 a constant.
Proof: Define the consistency error ε h ∈ Y h by It is sufficient show that, there exits a constant C > 0 such that Indeed, by the inf-sup condition, we have and we obtain the error estimate with C = 2 C .
Observation 2. The overall order of the discretization (8) is then one. Thus, to minimize the number of degrees of freedom the best choice is l = 1 and k = h.
We now can deduce the general convergence result,

Corollary 5.
Let v ∈ H 1 0 (Ω) the solution of problem (2). Then it holds where v h is the solution of the multilayer discretization (8).
Proof: As in the proof of Theorem 4, it holds Operators Π h and T h • Π h are stable, in the sense that Indeed, is straightforward to prove that ∇ H (Π h w) 0,Ω ≤ Π k w 1,Ω and then the first estimate follows as it because Π k is stable in H 1 (Ω) norm. The second one is a standard property of the Clément interpolation operator (cf. [16]).
Then a standard argument, using the density of D(Ω) in H 1 0 (Ω), proves that the r. h. s. of (31) tends to zero as h → 0.

Extension to domains with non-flat upper boundary
In this section we extend the multilayer discretization of the Poisson problem (1) to domains with non-flat upper boundary, that arise in geophysical flows.
Concretely, we consider domainsΩ that can be obtained as vertical transformations of the domain Ω = ω × (0, 1) (that we call in this section "reference" domain) by a change of variables of the form: with x ∈ ω ⊂ R 2 , z ∈ (0, 1) and where η is an smooth and strictly positive function defined in ω (the surface equations are thus z = η(x, y)). We assume that η ∈ W 1,∞ (ω) is a function such that for some constant η 0 . Note that he restriction in (33) on the gradient of η is consistent with the need of having η > 0 on ω.
We intend to solve the Poisson problem inΩ: To do it, we use the change of variables (32), to transform problem (34) in the following elliptic problem in the reference domain Ω in variational form: f (x,ẑ) and K is a symmetric and positive definite matrix given by K 3H = K T H3 and K 33 = 1 η 2 + |K H3 | 2 Here our objective is to deduce a multilayer approach of problem (34).

Multilayer discretization
To obtain our new multilayer system, we consider the same vertical decomposition of domain Ω and the same definitions and notations as those used in Section 2.
Similarly to the definition of the form a h by (9), we propose the following multilayer discretization of problem (34): where, the bilinear form a K,h (·, ·) is given by

Analysis of the multilayer discretization
We now study some properties of the multilayer problem (36). The stability of the multilayer discretization is achieved, under some restrictions on the gradient of function η, as follows: Theorem 6. We assume that function η satisfies (33). Then, the bilinear form a K,h (·, ·) satisfies the following properties: 1. There exists a constant C 1,K > 0 that depends on η, such that 2. There exists a constant C 2,K > 0 that depends on η, such that Proof: Now, using that η > 0, Young's inequality a b ≥ − 1 2 ( a 2 + 1 b 2 ) for the last integral and the estimates within the proof of (11) in Lemma 1 for the Then, if we take ε = 4 (1 + β) with β > 0 and we use that ∇ H η ∞,Ω < 1, we have where for the last inequality we use (11).

2.
Let v h ∈ X h and ϕ h ∈ Y h , then as η > 0, (11) and Holders inequality, we Thus, we obtain (39) with As consequence of this theorem we deduce the following well posedness result for the multilayer problem (36): Corollary 7. We assume that function η satisfies (33), then the multilayer problem (36) admits a unique solution v h ∈ X h that satisfies the estimate where C 3,K > 0 is a constant that depends on η.
The convergence of the multilayer discretization (36) is stated as follows: Theorem 8. We assume that η verifies (33) and that the weak solution of (35) satisfies v ∈ H l+1 (Ω) ∩ H 1 0 (Ω), with l ≥ 1. Then, the solution v h ∈ X h of the discrete problem (36) verifies the estimate ,Ω with C 4,K a constant that depends on η and Π h the operator defined by (24).
The proof of this result is an extension of that of Theorem 4 and Corollary 5 by standard techniques, we omit it for brevity.

Neumann boundary conditions
In this section we extend the multilayer discretization to the Poisson problem with Dirichlet boundary conditions (1) to the following Poisson problem with Dirichlet-Neumann boundary conditions: The variational form of this problem is: Given f ∈ H −1 (Ω) and g ∈ H −1/2 (Γ s ), with where Here we consider the same space of semi-discrete solutions in the vertical direction as in Section 2, that we endow with the following H 1 bl (Ω) discrete norm: This new discrete norm is motivated by the following new semi-discrete test functions space, as we see in Figure 3. Note thatŶ h is a subspace of H 1 bl (Ω) such that dim(Ŷ h ) = dim(X h ) and the interval [z N −1 , z s ] has length 3h/2.
We look for an approximation of the solution of (42) of the form v v h + g h with g h = g k (z − z s ) κ N and v h is the solution of the following multilayer discretization of (42): whereâ h is a the bilinear form defined bŷ We assume in particular g ∈ H 1 0 (ω) in order to have this last term well defined.
Consequently, the mappingT h is too an isomorphism between the normed spaces X h andŶ h . This justifies the choice of the new norm on space X h .

Well-posedness and convergence analysis
In this section we first study, analogously as we see in Section 3, the wellposedness of the multilayer problem (45): Theorem 9. We assume that g ∈ H 1 0 (ω), then the multilayer problem (45) admits a unique solution v h ∈ X h that satisfies the estimate Proof: (Sketch) Reasoning as in the proof of Theorem 1 for the multilayer problem (8), we deduce that the formâ h is stable and satisfies the inf-sup This new inf-sup condition is proved similarly to (20), as a consequence of estimates (46) and taking into account the differences between the forms a h andâ h mentioned in Observation 3.
Then problem (45) admits a unique solution that depends continuously on the data f and g.
for some constant C > 0 a constant.
Moreover, if the solution v of problem (42) is only in H 1 bl (Ω), then it holds

Numerical results
In this section we present some numerical 3D experiments, to analyze the computing time reduction obtained by the parallel computation of each multilayer discretization, so as to test the error estimates. We have used the FreeFem++ software (cf. [17]).
Concretely, for the Poisson problem (1), we solve the multilayer problem (8) by an iterative procedure. Taking the test functions ϕ h = ϕ α ⊗ σ α , 1 ≤ α ≤ N in (8), this problem is equivalent to the following linear system, with unknowns v α h , α = 1, · · · , N ∈ V 0,k (ω):  Figure 4 shows for each of these two tests two ratios, on the one hand, the ratio between the CPU time to solve the sequential 3D problem by the Conjugate Gradient Solver (CG) versus the mean CPU time used by processors to build the matrix and the second member of the 2D system that we solve in each layer/processor and solve it by the affine parallel GMRES algorithm (lines with legend "CPU TIME init+resol"). On the other hand, the ratio between the CPU time to solve the sequential problem by the CG Solver versus the mean CPU time used by processors only to solve the 2D system in each layer/processor by the affine parallel GMRES algorithm (lines with legend "CPU TIME resol") We notice that we obtain similar results for both tests. As the number of layers increases, there is a large gain of computing time, up to an optimal rate around 30/40 layers, that deteriorates as N Z increases from this value.
Furthermore, we test the error estimates. For that, Table 1 shows the obtained results for the relative errors in the L 2 (Ω) and the H 1 (Ω) discrete norms introduced in (10) We obtain first order convergence in the H 1 (Ω) discrete norm, in full agreement with our theoretical expectations. We also obtain second order convergence in the L 2 (Ω) discrete norm, as arises for Galerkin finite element solutions of regular elliptic problems.
We obtain quite similar convergence orders if we use unstructured horizontal meshes, that we not display for brevity.
that is the exact solution for the Laplace problem (34) inΩ, and we consider the previous function η with = 0.15. Figure 5 shows the global 3D Galerkin solution of (34) with h = 1/36.  Figure 6 shows the same two CPU ratios for each of the two tests shown in Figure 4 We observe a similar behavior as in Test 1, there is an increasing computational gain as the number of processors increases up to 40/50 layers (lines with legend "CPU TIME init+resol"), that further progressively deteriorates.The gain is somewhat higher than in the previous test, possibly due to the larger computational cost to obtain the 3D sequential global solution in the computational domainΩ.
Also, Table 2 display the norms of the relative errors, estimated orders of convergence, as well as the number of iterations used by GMRES, for N Z = 10, 20, 40, 80. We observe a quite similar behavior as in Test 1, we recover the theoretical first order in the H 1 (Ω) discrete norm, and second order in the L 2 (Ω) discrete norm. Here we do the tests for structured horizontal meshes, but we again obtain quite similar results if we use unstructured horizontal meshes.

Test 3: Dirichlet-Neumann boundary conditions
In this third test, we consider an exact solution of (41) in Ω with Neumann boundary conditions on Γ s and homogeneous Dirichlet boundary conditions on  Γ b ∪ Γ l : v(x, y, z) = 64 sin πz 2 x 2 (1 − x) y (1 − y) 2 . Figure 7 shows the global 3D Galerkin solution of (41) with h = 1/36. Figure 8 shows shows the same two CPU ratios for each of the two tests as in Figure 4 and Figure 6 We observe a similar behavior, there is an increasing computational gain as the number of processors increases to around 30/40 layers, that further progressively deteriorates. The gain is somewhat smaller as for the Test 1 due to  additional calculations required in the upper layers.

Conclusions
We have studied in this paper a multilayer discretization of second order elliptic problems, aimed at providing reliable multilayer discretizations of shallow fluid flow problems with diffusive effects. This is a Petrov-Galerkin discretization in which the trial functions are piecewise constant per horizontal layers, while the trial functions are continuous piecewise linear, on a vertically shifted grid.
We have introduced the discretization for the Poisson problem with Dirichlet boundary conditions in cylindric domains, and extended it to Neumann boundary conditions in cylindric domains, so as to domains with variable surface, that appear in free-surface fluid flow problems. We have proved the well posedness and optimal error order estimates for these three discretizations in natural norms, based upon specific inf-sup conditions.
We have performed numerical tests with parallel computing of the solution by a block-Jacobi algorithm, based upon the multilayer structure of the discretization, for academic problems with smooth solutions. We recover the theoretical optimal error order convergence, and observe a high increase of the speed of computations for a moderate number of processors. This confirms the interest of applying the technique introduced to multilayer discretizations of fluid flow problems.
In further steps we shall improve the parallelization procedure to obtain a good scaling of the CPU computing time for large number of processors. Also, we will apply the discretization introduced to convection-diffusion problems, and then to multilayer discretizations of Navier-Stokes and related equations.
We will also extend the discretization to higher-order approximations of the unknowns. These works are now in progress.