Iterative solution to the biharmonic equation in mixed form discretized by the Hybrid High-Order method

We consider the solution to the biharmonic equation in mixed form discretized by the Hybrid High-Order (HHO) methods. The two resulting second-order elliptic problems can be decoupled via the introduction of a new unknown, corresponding to the boundary value of the solution of the first Laplacian problem. This technique yields a global linear problem that can be solved iteratively via a Krylov-type method. More precisely, at each iteration of the scheme, two second-order elliptic problems have to be solved, and a normal derivative on the boundary has to be computed. In this work, we specialize this scheme for the HHO discretization. To this aim, an explicit technique to compute the discrete normal derivative of an HHO solution of a Laplacian problem is proposed. Moreover, we show that the resulting discrete scheme is well-posed. Finally, a new preconditioner is designed to speed up the convergence of the Krylov method. Numerical experiments assessing the performance of the proposed iterative algorithm on both two- and three-dimensional test cases are presented.

Equation (1) typically models the bending of a clamped plate supporting a load. By introducing the unknown ω := −∆ψ, (1a) can be rewritten into two second-order elliptic equations, yielding the mixed formulation This mixed form naturally arises in fluid dynamics, where ψ represents the stream and ω the vorticity. In the plate bending model, ψ represents the deflection and ω the bending moment or shear resultant force. There are several advantages in employing the mixed formulation (2) over the primal form (1). First, while the weak solution of the primal form (1) is to be found in H 2 (Ω), that of the mixed one lies in H 1 (Ω), for which approximation spaces are easier to be constructed. Second, the splitting (2) allows, after introducing an additional unknown, the use of fast, scalable solvers available for second-order elliptic equations. In particular, we will consider here the Hybrid High-Order (HHO) method [9,10], a non-conforming polyhedral discretization allowing arbitrary polynomial degrees of approximation,

HHO discretization of the Laplacian problem
In this section, we briefly recall (see, e.g., [8,Chap. 2] for extended details) the discrete HHO formulation of the following problem: find u : Ω → R such that where f : Ω → R and g D ∈ H 1/2 (∂Ω). A weak solution of (8) is obtained via the variational formulation: find where the bilinear form a is such that a(v, w) = (∇v, ∇w), for all v, w ∈ H 1 (Ω).

Mesh definition and notation
Let the couple (T h , F h ) define a mesh of the domain Ω ⊂ R d , d ∈ {1, 2, 3}: T h is a set of disjoint, open, polyhedral elements such that T ∈T h T = Ω; F h is the set of element faces; h := maxT ∈T h hT with hT denoting the diameter of T ∈ T h . The mesh is assumed to match the geometrical requirements of [8,Def. 1.4]. Let T B h the subset of T h collecting the elements located at the boundary of the domain. We also define the following subsets of F h : • F I h , collecting the interior faces; • F B h , collecting the boundary faces; • FT , collecting the faces of T , for all T ∈ T h . We denote by n ∂T the unit normal vector to ∂T pointing outward of T . For all T ∈ T h (resp. F ∈ F h ), we denote by (·, ·)T (resp. ⟨·, ·⟩F ) the standard inner product of L 2 (T ) (resp. L 2 (F )) or L 2 (T ) d . We define (·, ·)T h :

Local and broken polynomial spaces
The HHO method hinges on discrete unknowns representing polynomial functions local to elements and faces. So, for all m ∈ N0 and all T ∈ T h (resp. F ∈ F h ), we denote by P m (T ) (resp. P m (F )) the space spanned by the restriction to T (resp. F ) of d-variate polynomials of total degree ≤ m. From these local polynomial spaces, we can construct the following broken polynomial spaces supported by the mesh and its skeleton: respectively. The local space P m (FT ) is defined analogously for all T ∈ T h . For all cell or face X, we denote by π m X : L 2 (X) → P m (X) the local L 2 -orthogonal projector onto the space P m (X). By patching up those local projectors, we denote by π m T h : L 2 (Ω) → P m (T h ) the piecewise L 2 -orthogonal projector onto P m (T h ), and by π m

Discrete hybrid formulation
Given the polynomial degrees k ∈ N0 and l ∈ {k, k + 1}. The global and local spaces of hybrid variables are defined as respectively. For any v h ∈ U h , we denote by v T ∈ U T its restriction to T ∈ T h . Boundary data are strongly accounted for in the following subspaces: In particular, homogeneous Dirichlet conditions are strongly enforced in U h,0 . The global HHO bilinear form associated to the variational formulation of problem (8) is defined as a h : In this expression, the first term is responsible for consistency while the second is required to ensure stability of the scheme. The consistency term involves the local potential reconstruction operator p k+1 T reconstructs an approximation of v of degree k + 1. The stabilization bilinear form sT depends on its argument only through the difference operators δT : U T → P l (T ) and δT F : These operators capture the higher-order correction that the reconstruction p k+1 T adds to the element and face unknowns, respectively. A classical expression for sT : If l = k + 1, one can also use the simpler formula (used, e.g., in [4]) The global HHO problem reads: find u h ∈ U h,g D such that The final approximation is obtained through the post-processing step

Local problems and cell unknown recovery operator
It follows from this local construction that the cell unknowns are only locally coupled. For all T ∈ T h , we define the linear operator ΘT : We denote by ΘT the cell unknown recovery operator. In order to shorten the notation of the hybrid couple v T : Finally, the associated global operators ΘT h F h : P k (F h ) → P l (T h ) and Θ h : P k (F h ) → U h are defined locally such that for all T ∈ T h : Θ h verifies the following useful property: Proof. See Appendix A. (15) is solved through the equivalent, condensed formulation:

Discrete normal derivative
In this section, we derive an approximation, in the HHO context, of the normal derivative on ∂Ω. This formula will play a crucial role in the approximation scheme of the biharmonic problem. Let u be the solution of the boundary value problem (8) and let H : H 1/2 (∂Ω) → H 1 (Ω) be a linear operator such that for all v ∈ H 1/2 (∂Ω) defined on the boundary, Hv extends v in the interior of Ω. By Green's formula, it holds that ⟨∂nu, v⟩ = (∇u, ∇Hv) + (∆u, Hv) ∀v ∈ H 1/2 (∂Ω).
Using the bilinear form a(·, ·), and given that −∆u = f , the above equation becomes In the literature, the operator g D → ∂nu is called Poincaré-Steklov operator, or Dirichlet-to-Neumann map; see, e.g. [24]. Equation (23) is the well-known variational formula for the computation of the normal derivative. Now, given the solution u h of the discrete problem (15), we define its normal derivative on the boundary faces, denoted by , through a discrete counterpart of (23) in the HHO setting. Namely, where H h is an extension/lifting operator expressed in the hybrid setting, i.e.
Then, we set

Remark 2. (Condensed formula)
In order to facilitate the practical computation of ∂ n,h (u h ), equation (24) can be rewritten in terms of the condensed bilinear form a h (cf. (22)). First of all, it follows from the definition of Θ h that H h = Θ h HF h . Then, by using property (20) and the definition (22) of a h , we have

Equation (24) then becomes
Notice that compared to (24), which requires the knowledge of u h , formula (25) only involves uF h . Consequently, the following hold: (i) uT h does not have to be computed, (ii) the matrix used to compute the first term is smaller, as the matrix representation of a h is a Schur complement where the cell unknowns have been eliminated.
The computation of the discrete normal derivative as proposed here is numerically validated in Section 7.2. The experiments show a convergence O(h k+1 ) in L 2 -norm.

Discrete HHO problem
In this section, the global, discrete realization of problem (6) is formulated.

Bilinear form
In order to devise a stable problem, we will make use of the following bilinear form in the hybrid space Notice that (26) defines an inner product-like bilinear form based on the L 2 -inner product of the cell unknowns, to which a stabilizing term has been added. The stabilizing term is inspired from [4,Eq. (26)]. Notwithstanding the exponent of the hF factor, the latter reproduces, locally in T , the stabilization term (13) of the Laplacian bilinear form. The scaling factor hF is selected to ensure dimensional homogeneity with the consistency term it is added to. Note that (·, ·) ⋆ T h does not define an inner product in U h , due to the fact that only the boundary cells are involved in the stabilizing term. In spite of that, we allow ourselves the use of an inner product notation, since (·, ·) ⋆ T h shall undertake the role of the L 2 -inner product in the discrete problem and, thus, carries the same semantics.
As the stabilizing term of (26) is based on the same local bilinear form sT F used to stabilize the bilinear form a h , it naturally confers to (·, ·) ⋆ T h some useful properties. First, (·, ·) ⋆ T h remains symmetric and positive semi-definite. Second, if definiteness is not globally ensured, it is ensured "locally" in the cells where the stabilization term is applied. One might say that (·, ·) ⋆ T h enjoys a property of local stability in the boundary cells, in the sense that Let us now describe the discrete, variational form of the continuous operator (7).

Discrete scheme
The discrete formulation of problem (6) reads: We also define the linear operator L h : Since ℓ h is positive-definite, problem (38) is well-posed. Its solution exists and is unique. Once it is solved, we Finally, the discrete approximation of (ω, ψ) is given by T h is required within ℓ h for problem (38) to be well-posed, it has to be used neither for the enforcement of the source functions in (37) and (40), nor for the computation of ∂ n,h (ψ 0 h ) (via formula (24)) in the right-hand side of (38).

Preconditioner
In this section, we introduce the preconditioned iterative strategy for the solution of the algebraic realization of (38).

Algebraic setting
For a fixed basis for the discrete HHO spaces, let N B := dim P k (F B h ) . Introduce the algebraic counterpart L h : R N B → R N B of the linear operator L h defined by (39). We recall that L h can be viewed as a matrix in R N B ×N B , whose coefficients are not explicitly known, but whose application to a vector can be computed. The algebraic realization of problem (38) reads: find λ ∈ R N B such that where b ∈ R N B is an algebraic representation of ∂ n,h (ψ 0 h ) − g N in the chosen basis.

An approximate, sparse matrix
The implicit system (41) is solved using a preconditioned, flexible conjugate gradient (PFCG) algorithm. This section describes the construction of the preconditioner. The use of a flexible version of the Krylov method is justified by the preconditioner being non-symmetric [1,3]. Refer to Section 6.4 for a discussion about the symmetric property of the preconditioner. Introducing (ej) j=1,...,N B the canonical basis of R N B , the j th column of L h is given by the evaluation of L h (ej). Note that L h (ej) is a dense vector. Consequently, the explicit computation of L h would result in a dense matrix. The preconditioner proposed here consists in the explicit computation of a sparse matrix L h approximating L h . In a nutshell, while the computation of each column L h (ej) involves solving two Laplacian problems in the whole domain, our approximation L h (ej) consists in solving those problems in a restricted neighbourhood of the j th DoF.
Let j ∈ {1, . . . , N B } a fixed unknown index. Let F j ∈ F B h the face supporting the DoF associated to j, and T j ∈ T h the unique element owning F j . Define T j h ⊂ T h a neighbourhood of T j : roughly, a set of elements around T j , including T j . Let Ω j be the open, connected set such that Ω j := ∪ T ∈T j h T . Figure 1a illustrates (in grey) such a neighbourhood, where F j is represented by a thick line and T j by a darker triangle. Relatively to Ω j , one defines the sets of interior and boundary faces by Supposing that the j th basis function of P k (F B h ) corresponds to the k th basis function of P k ((F B h ) j ), we compute L j h ( e k ). As L j h ( e k ) yields the normal derivative only on ∂Ω j ∩ ∂Ω, we build the final, approximate column L h (ej) from L j h ( e k ) by setting zero coefficients where the domain boundary is not covered by the neighbourhood boundary. We claim that L j h (ej) yields a sparse approximation of L h (ej). To support our claim, in Figures 1b and 1c we plot, in iso-value curves, the solutions to the first and second Laplacian problems computed in Ω j . Moreover, Figures 1d and 1e plot the solutions of the Laplacian subproblems computed in the whole domain. We observe that, owing to the homogeneity of the first equation and the homogeneous Dirichlet condition everywhere except on one face, the solution culminates on that face and decreases towards zero as it goes further away. Roughly speaking, our strategy consists in approximating the solution locally in Ω j while imposing it to be zero in the rest of the domain. Then, this first truncated solution is used as the source function in the second problem, also solved in Ω j . Here, note that the quantity of interest is not the solution itself, but its normal derivative on the boundary. In particular, the neighbourhood of the considered face concentrates the most information, insofar as the solution flattens as it goes away from it. In that sense, L j h (ej) yields a reasonable approximation of L h (ej).

Computation in practice
Regarding the actual computation of L h , one can note that the HHO matrix blocks used to solve the Laplacian problems in Ω j and to compute the normal derivative can be extracted from the global ones. The sparsity of L h is controlled by the number of boundary faces in the chosen neighbourhood, while the accuracy of the approximation depends on the size of the neighbourhood. In particular, remark that choosing T j h := T h for all j yields the actual matrix L h . In practice, we construct T j h by adding successive layers of neighbours around T j . The number of layers is defined by the parameter α ∈ N0. T j h is defined as Note that this definition understands the neighbouring relationship as having at least one vertex in common. The neighbourhood represented in Figure 1a corresponds to α = 3.

Discussion on symmetry
Choosing a constant α for all j allows to preserve a symmetric sparsity pattern. Indeed, considering (i, j) ∈ {1, . . . , N B } 2 , ( L h )ij ̸ = 0 implies that F i ⊂ ∂Ω j . Consequently, ( L h )ji ̸ = 0 implies that, reciprocally, F j ⊂ ∂Ω i . However, symmetry itself is not preserved, inasmuch as ( L h )ij results from computations in the neighbourhood Ω j whereas ( L h )ji results from computations in Ω i ̸ = Ω j .
To enforce the symmetry of the matrix L h , one must ensure that for all couple (i, j) ∈ {1, . . . , N B } 2 , Ω i = Ω j . This requirement leads to a patch-based method: the set of boundary faces must be split into patches, i.e. connected subsets of boundary faces. For each patch, one single neighbourhood must be defined to process all DoFs included in that patch. This yields a block-diagonal matrix L h , one block corresponding to each patch. Algebraically, such a preconditioner is in fact an approximate block Jacobi preconditioner, insofar as each diagonal block approximates  the corresponding diagonal block of the exact matrix L h . This symmetric preconditioner was numerically tested. Even with the facts that a non-flexible Krylov method can be used and that the matrix L h is easier to factorize, the convergence of the method is significantly slower than with the unsymmetric preconditioner defined in Section 6.3. Therefore, only the latter has been used in the numerical experiments of Section 7.

Experimental setup
The implicit system (41) is solved iteratively by the PFCG method, using the preconditioning technique described in Section 6. Given the iterate λ ∈ R N B , the corresponding residual vector is defined as r := b − L h ( λ). The PFCG algorithm stops when ||r||2/||b||2 < ε, where ε > 0 is a fixed tolerance and || · ||2 denotes the Euclidean norm on R N B . All experiments are conducted choosing l = k, but we stress that the same qualitative results are obtained with l = k + 1. In 2D, the linear systems corresponding to the Laplacian subproblems are solved by Cholesky factorization; in 3D, we use a p-multigrid algorithm on top of the algebraic multigrid method designed in [11]. In that latter case, the same tolerance ε must be used to stop the Laplacian solvers and the global PFCG algorithm. To apply the preconditioner, we solve L h using the BiCGSTAB algorithm with tolerance ε. The computations are run on an 8-core processor (AMR M1 Pro) clocked at 3228 MHz.

Preliminary: convergence of the discrete normal derivative
We first assert the validity of the discrete normal derivative proposed in Section 4 by evaluating its order of convergence. Let Ω := (0, 1) 2 and u : (x, y) → sin(4πx) sin(4πy) the manufactured solution of the boundary value problem (8), where f and g D are defined accordingly. The problem is discretized on a sequence of successively refined meshes, and we consider ∂ n,h (u h ) computed by (25). The exact normal derivative ∂nu is known, and we assess the relative L 2 -error ∥∂nu − ∂ n,h (u h )∥ L 2 (∂Ω) /∥∂nu∥ L 2 (∂Ω) for k ∈ {0, 1, 2, 3}. Figures 2a and 2b present the experimental results on uniform Cartesian meshes and on unstructured, triangular Delaunay meshes, respectively. These experiments assess a convergence in O(h k+1 ) (also observed for l = k + 1). Note that in the case of a non-convex domain   with re-entrant corners, those convergence orders are expected to be reduced due to the corner singularities [23]. One way to mitigate this issue is the use of graded meshes (cf. [23]).

Convergence rate of the scheme
On the unit square, we consider ψ(x, y) = x sin(πy)e −xy and its corresponding ω = −∆ψ, the manufactured solution of (2), with f, g D and g N determined accordingly. The square is decomposed into a sequence of successively finer meshes, and the numerical scheme is applied for k ∈ {0, 1, 2, 3}. Figure 3 presents the evolution of the relative L 2 -error achieved by the approximation (ω h , ψ h ) with respect to (ω, ψ), on a sequence of Cartesian meshes. Similarly, Figure 4 shows the analogous results obtained with a sequence of polygonal meshes, each constructed from an initial Cartesian mesh by agglomeration and face collapsing. One such mesh is represented in Figure 6. We empirically interpret these results in light of Glowinski where C is a constant independent of h and ψ. In our HHO setting, the above error estimate would translate in the following one: if k ≥ 2, then The experimental results of Figure 3 validate the estimate (43), not only for k ≥ 2, but also for k ∈ {0, 1}. More precisely, while Figure 3b shows that the estimate is sharp for ψ (i.e. O(h k+2 )), Figure 3a shows, for ω, a faster convergence than indicated by (43). Namely, the convergence seems to be O(h k+1 ) for k = 0 and O(h k+1/2 ) for k ≥ 1.
The experiments of Figure 6 on polygonal meshes exhibit the same convergence orders, except for ω h with k = 0, where the error estimate reduces to O(h k+1/2 ), in accordance with the other values of k. Based on these experiments, we conjecture that the following error estimate is sharp for all k ≥ 0: Nonetheless, we keep in mind that k = 0 seems to be a special case for ω h , where superconvergence may be observed: besides the order 1 obtained on Cartesian meshes (Figure 3a), we observe the order 2 if we choose as an exact solution the polynomial function ψ(x, y) = x 4 (x − 1) 2 y 4 (y − 1) 2 , as illustrated in Figure 5. [17]) The HHO discretizations of [17] rely on the primal formulation (1) of the equation. They hinge on the local approximation space P k+2 (T ) × P k+1 (FT ) × P k (FT ) in 2D (resp. P k+2 (T ) × P k+2 (FT ) × P k (FT ) in 3D), where the associated set of DoFs aims at approximating (ψ |T , ψ |∂T , ∇ψ · n ∂T ). With this setting, the methods achieve a convergence order of k + 3 in L 2 -norm. In comparison, the present method requires the space P k+1 (T ) × P k+1 (FT ) to obtain the same convergence order, i.e. one less polynomial degree for the cell unknowns and no DoF dedicated to the approximation of ∇ψ · n ∂T . However, if our method is structurally lighter in DoFs than [17], one cannot conclude regarding the computational cost, since those DoFs have to be solved multiple times as part of an iterative process (vs. only once for [17]).

Preconditioner
We evaluate in this section the efficiency of the preconditioner designed in Section 6. We recall that the parameter α ∈ N0 is an indicator of the size of the neighbourhoods used for the approximation of the matrix.

Convergence vs. cost: the parameter α
The considered test case is the square domain partitioned into 256×256 Cartesian elements, with k = 1. The tolerance is set to ε = 10 −14 . Let us assess how the preconditioner improves the speed of convergence. Figure 7a plots the decay of the algebraic residual ||r||2/||b||2 without preconditioner (in red) and with preconditioner (in blue), considering increasing values of α. One can confirm that, as expected, the higher the value of α, the better the convergence rate. However, increasing α increases the setup cost of the preconditioner. Indeed, (i) the Laplacian problems are solved in larger subdomains, (ii) the approximate matrix L h is denser, which increases the cost of its factorization. Figure 7b presents the computational cost, measured in CPU time, to reach the algebraic solution. One can see that a trade-off between cost and convergence rate must be made to achieve optimal performance. If raising α is a good strategy up to α = 6, above that threshold, the better convergence rate does not make up for the cost of the setup phase, hence making the overall CPU time slightly increase. We stress that other choices of k yield the same qualitative results.

PFCG convergence rate: k-independence and h-dependency assessment
Problem (2) on the unit square with exact solution ψ(x, y) = x sin(πy)e −xy is solved on a sequence of 2D grids composed of N 2 Cartesian elements, N ∈ {32, 64, 128, 256, 512}. Table 1 presents, for k ∈ {0, 1, 2, 3}, the number of PFCG iterations executed to reach the convergence criterion with ε = 10 −8 . The parameter of the preconditioner is set to α = 8. For comparison, the number of iterations without preconditioning is reported in brackets. Firstly, one can remark that the use of the preconditioner allows a convergence rate independent of k, while a dependency is observed without preconditioner. Secondly, we approximate the asymptotic convergence rate with respect to the problem size by using the number of iterations measured on the two finer meshes. This yields a dependency in O(h −0.41 ). Additionally, Table 2 reports the CPU times of the setup and iteration phases. The setup includes the Cholesky factorization of the Laplacian matrix and the assembly of the preconditioning matrix. In brackets, the same  information without preconditioning is displayed. Based on these results, one can see that the compensation of setup cost by the gain in iteration time is more and more efficient as the problem grows larger. The equivalent experiment is performed in 3D (exact solution ψ(x, y) = xz sin(πy)e −xy ) with unstructured tetrahedral meshes. As the preconditioner is costlier in 3D, we choose α = 2, and set the tolerance to ε = 10 −5 . The number of PFCG iterations (independent of k) is presented in Table 3, which exhibits a dependency in O(h −0.37 ). The CPU times of setup and iteration phases are reported for k = 0.  Table 3: Test case: Cubic domain, unstructured tetrahedral mesh, ε = 10 −5 . Number of PFCG(α = 2) iterations, along with the setup and iteration CPU times (in seconds) for k = 0.