#
On preconditioning saddle point systems with trace constraints
coupling 3 and 1 domains – applications to matching and nonmatching FEM
discretizations
^{1}

^{1}

## Abstract

Multiscale or multiphysics problems often involve the coupling of partial differential equations posed on domains of different dimensionality. In this work we consider a simplified model problem of a 3-1 coupling and the main motivation is to construct algorithms that may utilize standard multilevel algorithms for the 3 domain, which has the dominating computational complexity. Preconditioning for a system of two elliptic problems posed, respectively, in a three dimensional domain and an embedded one dimensional curve and coupled by the trace constraint is discussed. Investigating numerically the properties of the well-defined discrete trace operator, it is found that negative fractional Sobolev norms are suitable preconditioners for the Schur complement. The norms are employed to construct a robust block diagonal preconditioner for the coupled problem.

SISCxxxxxxxx–x

reconditioning, saddle-point problem, Lagrange multipliers, trace

65F08

## 1 Introduction

Let be a bounded domain in 3, while represents a 1 structure inside , and consider the following coupled problem

(1a) | |||||

(1b) | |||||

(1c) |

Here the term is to be understood as a Dirac measure such that for a continuous function . We remark that from a mathematical point of view the trace of required in (1c) is in the continuous case not well-defined unless the functions are sufficiently regular. For simplicity, the system shall be considered with homogeneous Neumann boundary conditions.

The system (1) is relevant in numerous biological applications where the embedded (three dimensional) structure is such that order reduction techniques can be used to capture its response by a one dimensional model. Equation (1a) then models processes in the bulk, while (1c) is the coupling between the domains. A typical example of such a system is a vascular network surrounded by a tissue and the order reduction is due to assumption of radii of the arteries being negligible in comparison to their lengths. To list a few concrete applications, the 3-1 models have been used, e.g., in [17, 21, 16, 31] to study blood and oxygen transport in the brain or in [11] to describe fluid exchange between microcirculation and tissue interstitium. Efficiency of cancer therapies delivered through microcirculation was studied in [10], and hyperthermia as a cancer treatment in [27]. We note that the employed models are more involved than (1), but that the system still qualifies as a relevant model problem.

Due to the measure term and the three-to-one dimensional trace operator, the problem (1) is not standard and establishing its well-posedness is a delicate issue. In fact, considering (1a) with a known and homogeneous Dirichlet boundary conditions, the equation is not solvable in , as may be unbounded in the neighborhood of . A similar problem was studied in [14], where two elliptic problems were coupled via a measure source term, and a unique weak solution was found using weighted Sobolev spaces. In particular, the weighted spaces ensured that the trace could be defined as a bounded operator. A corresponding finite element method (FEM) for the problem was discussed in [13], where optimal convergence in the weighted Sobolev norm was shown using graded meshes. For an elliptic problem with measure data, it was shown in [19] that FEM with regular meshes yields optimal convergence in the norm outside of the fixed neighborhood of the singularity. We note that the more application oriented works [11, 10, 27], that build on the analysis in [14, 13], relied on incomplete LU preconditioning.

In the current paper we shall assume that (1) is well posed and the focus is then on the construction of optimal preconditioners for the linear system due to (1) and FEM. Because the computational complexity of the 3 problem dominates the 1 problem, we put focus on preconditioners that are composed of standard multilevel algorithms for the 3 problem. This means that the weighted Sobolev spaces, or extra regularity in the equation in the 3 domain (1a), is disregarded, and that we rather add an extra requirement to (1c). We note that shall be approximated within conforming finite element spaces and as such the approximation has a well defined trace.

The current paper is an extension of [20], where a system similar to (1a)–(1c) was analyzed for the case a bounded domain in 2 and a structure of codimension one. Therein, robust preconditioners were established, based on the operator preconditioning framework [26], in which preconditioners are constructed as approximate Riesz mappings in properly chosen Hilbert spaces. The framework often allows for construction of order-optimal preconditioners, with convergence independent of material and discretization parameters, directly from the analysis of the continuous system of equations. In particular, in [20] it was shown that the proper preconditioning relied on a nonstandard fractional inner product. Crucial for the analysis was the fact that the trace operator is a well-defined mapping between and , when is of codimension one with respect to . Furthermore, for the finite element approximation in [20], it was assumed that the discrete meshes representing and matched in the sense that the cells of the mesh of were edges in the mesh of . Finally, only continuous linear Lagrange elements were used.

This paper utilizes ideas presented in [20] for the construction of the preconditioner, but here we go beyond what was theoretically established. In particular, we consider the case where is of codimension two with respect to . As the trace operator is not well-defined in the continuous case, we do not attempt to provide mathematically rigorous proofs, as this would require additional regularity of the solution. Instead, we study for which the inner product provides numerically stable behaviour. In addition, we consider the case where the discretizations of and do not match. Finally, an approximation by discontinuous elements is discussed.

Our work is structured as follows. In §2 the theoretical background is presented. Section 3 discusses numerical experiments using spectral and finite element discretizations that identify suitable norms for the discrete 3-1 trace operator. In §4 the identified norms are employed to construct optimal preconditioners for coupled model 3-1 problems discretized with FEM and matched discretizations of and . In §5 this restriction is lifted, the corresponding inf-sup condition is discussed, and we present numerical experiments that suggest the identified norms lead to good preconditioners. Finally, conclusions are drawn in §6.

## 2 Notation and preliminaries

Let be a Hilbert space of functions defined on a domain , . The norm of the space is denoted by , while is the duality pairing between and its dual space . We let denote the inner product of , while, to simplify the notation, is the inner product. The Sobolev space of functions with square integrable derivatives is . Finally, denotes the closure of the space of smooth functions with compact support in in the norm.

We use normal capital font to denote operators over infinite dimensional spaces, e.g. . For a discrete subspace , , the subscript is used to distinguish the finite dimensional operator due to the Galerkin method, e.g., defined by

For a given basis, of , the matrix representation of the operator is denoted by sans serif font. Thus is represented by with entries

Finally, the function is represented in the basis by a coefficient vector , where (summation convention invoked).

### 2.1 Properties of the trace operator

We consider an open connected domain with Lipschitz boundary and a Lipschitz submanifold of codimension one or two in . The trace operator is defined by for .

In case the codimension of is one, the properties of the trace operator are well known. In particular, is bounded and surjective, see, e.g., [1, ch. 7], where is a fractional Sobolev space equipped with the norm

Moreover, is bounded, but not surjective. To define the trace over as a surjective operator, the range is given as ,

and is some suitable extension of , e.g., , in which case . We refer to [20] for these results.

The integral norms of and can be expensive to compute. For construction of efficient numerical algorithms, it is therefore more suitable to relate the spaces to interpolation spaces, see [22, 5] or [20]. For the sake of completeness, we review here the presentation from [20]. Let . For fixed is in and by the Riesz-Fréchet theorem there is a unique such that for any . The operator is injective and compact and thus the eigenvalue problem (no summation implied) is well-defined. In addition, is self-adjoint and positive-definite such that the eigenvalues form a nonincreasing sequence and . By definition, the eigenvectors satisfy

or equivalently

(2) |

Further, the set of eigenvectors forms a basis of , which is orthogonal in the inner product of and orthonormal in the inner product. Finally, for we define the -norm of as

(3) |

The space is defined as the closure of the in the -norm, while is then defined analogically to with in the construction. We remark that , , and the norms of the spaces are equal. Moreover and with the equivalence of norms.

Following the approach in [20], a weak formulation of the homogeneous Dirichlet problem for (1)–(1c) with , of codimension one, using the method of Lagrange multipliers, reads: Find such that

(4) | ||||||

Letting , the well-posedness of (4) is guaranteed as the Brezzi conditions are satisfied with , see [20] for the proof in 2-1 setting, which immediately generalizes to 3-2. Crucial for the well-posedness is the fact that is an isomorphism. Consequently, [26] is invoked to yield a block diagonal preconditioner for the discretized problem where individual blocks are conceived as approximations of the corresponding Riesz mappings.

Let now . Considering (1) on the finite dimensional spaces, we obtain a variational problem: Find such that

(5) | ||||||

Here the discrete trace operator is well-defined as the functions in are continuous. The continuous problem, on the other hand, is not well defined since does not to permit a bounded trace in , see [14]. In turn we cannot directly follow the steps of [20] and employ operator preconditioning [26] to construct an optimal preconditioner. Instead, we shall reason about the properties of the discrete system.

From a linear algebra point of view, the problem (5) is a saddle-point system

Block diagonal preconditioners for such problems can be constructed as an approximate inverse of the matrix , where should be spectrally equivalent with and should be spectrally equivalent with the Schur complement , see, e.g., [32, 33]. Considering (5), the key question is thus whether it is possible (in an efficient and systematic manner) to construct an operator that is spectrally equivalent with the Schur complement. Motivated by the 2-1, the operator shall be based on the fractional -norm (3).

Following [20], the discrete approximation of the -norm shall be constructed by mirroring the continuous eigenvalue problem (2). More specifically, let and matrices , be the representations of , ; the Galerkin approximations of operators , from (2). Then there exists an invertible matrix and diagonal, positive-definite matrix satisfying . Moreover, the product is an identity such that the columns of form an orthogonal and orthonormal basis of . In order to define the discrete norm, we let be a symmetric, positive-definite matrix

(6) |

The matrices are defined analogically to (6), using the eigenvalue problem for the Laplace operator with homogeneous Dirichlet boundary conditions. For represented in the basis of the space by a coefficient vector , let be the representation of in the basis of eigenvectors, that is, . We then set

(7) |

The generalized eigenvalue problem required for evaluating the discrete -norm (7) becomes trivial if the approximation space is such that , i.e. the basis is formed by the eigenvectors of the continuous problem (2). Such a discretization is practically limited to Cartesian domains, however, it will prove useful in studying the trace operator when the codimension of in is two. The technique is introduced in Example 2.1.

###### Example 2.1 (Spectral method for 2-1 coupled problem)

Let , and consider the task of finding which minimizes subject to and on the boundary (in the sense of traces). Introducing , the problem is formulated as a saddle point system for , satisfying

(8) | |||||

Well-posedness of (8) is readily established by verifying the Brezzi conditions [9]. In particular, the inf-sup condition can be shown to hold, see, e.g., Appendix A. By operator preconditioning [26] the canonical preconditioner for (8) is the Riesz map with respect to the inner product inducing the norm of .

The Galerkin approximation of (8) is defined with spaces , such that and . Recall that functions are the eigenfunctions of (2) satisfying on . For greater readability, let us introduce . In the basis of eigenfunctions, the discrete trace operator is represented by a trace matrix with entries

Here is a column index . With matrix and vectors , defined in a natural way, the preconditioned linear system from (8) is

(9) |

where the first matrix is the discretization of the canonical preconditioner. We note that matrices are diagonal. Consequently .

For stability of the solution obtained by solving (9), it is required that the discrete inf-sup condition holds. Note that for the trace matrix does not have a full row rank and thus is necessary. We shall set and show that for this choice the inf-sup condition is satisfied.

By, e.g., [24] or [7], the constant of the discrete inf-sup condition is the smallest eigenvalue of the generalized eigenvalue problem for the negative Schur complement of the system matrix, i.e.,

where is the sought eigenpair. The simple structure of the involved matrices allows us to compute all the eigenvalues of the problem analytically. In fact,

Then and the lower bound can be evaluated. Using Mathematica [36], we have obtained as the discrete inf-sup constant. For the largest eigenvalue, the following estimate can be established

and in turn the theoretical spectral condition number for the preconditioned Schur complement is . We remark that the remaining Brezzi constants are both equal to one and the condition number of (9) can be determined from spectral bounds presented in [32].

The theoretical findings about the condition number of the preconditioned Schur complement are confirmed by numerical experiments, with results shown in Figure 1 and Table 1. The figure shows that there is a range of exponents for which the condition numbers are stable. Interestingly, in this range gives the largest condition number while for a slightly smaller value is observed, cf. Table 1.

1.9848 | 1.9959 | 1.9987 | 1.9997 | 1.9999 | 2.0000 | 2.0000 | |

1.7189 | 1.7190 | 1.7190 | 1.7190 | 1.7190 | 1.7190 | 1.7190 |

## 3 Norms for the discrete 3-1 trace

In Example 2.1 a priori knowledge of the trace space lead to an optimal preconditioner for the model problem (8). In particular, the norm of the trace space was used to construct a spectrally equivalent operator to the Schur complement of (9). For and a one dimensional curve, the trace space is not a priori known and we shall therefore attempt to characterize it numerically. To this end, we shall at first use the spectral discretization and search for the -norm (7) for which the condition number of the preconditioned Schur complement is bounded in the discretization parameter. We note that the condition is motivated by the fact that convergence of the preconditioned conjugate gradient method is estimated in terms of the condition number, see, e.g., [35]. For suitable the linear system with the Schur complement could thus be solved efficiently. We also note that the condition is weaker than spectral equivalence. In fact, if such exists, the matrix is spectrally equivalent with the Schur complement if and only if one of the extremal eigenvalues is bounded by a constant.

### 3.1 Trace operator with spectral discretization

Let , and , cf. Example 2.1. We consider the problem of minimizing , , subject to on the boundary and the constraint on , where the trace operator restricts either to or respectively. The weak formulation of the problem reads

(10) | |||||

and is equivalent with a linear system

(11) |

In (11) the trace matrix for curve is sparse with entries

Here was introduced for readability. Note that for the matrix does not have a full row rank and the system is singular. We therefore set . For the trace matrix is sparse with a more involved sparsity pattern and at most four nonzero entries per row

Finally, we consider the generalized eigenvalue problem for the Schur complement of (11) and matrices where such exponents are of interest, for which the spectral condition number is bounded in the discretization parameter. Note that with the Schur complement is a diagonal matrix ,

(12) |

For the matrix is dense and shall be computed from assembled terms. As such a smaller is explored in this configuration.

The results of the numerical experiments with are summarized in Figure 2. We observe that values yield bounded condition numbers for . The condition numbers are not quite converged for the other configuration, however, it is possible to identify unstable exponents . Moreover, the values close to appear to be stable also in this configuration. This fact is easier to appreciate in Table 2, which shows , and as functions of the discretization parameter for . With the condition number is evidently constant, while for the number appears to be bounded. Note that with both configurations the smallest and largest eigenvalues are not bounded and thus is not spectrally equivalent with the Schur complement with the bounds independent of . However, any of , (or their linear combinations) define a mesh-depenedent scale that yields spectral equivalence. Such scale, however, is not easily computable in general.

0.6218 | 2.0696 | 3.3285 | 0.8476 | 1.9916 | 2.3496 | ||
---|---|---|---|---|---|---|---|

0.9167 | 3.0511 | 3.3285 | 1.0298 | 2.4283 | 2.3581 | ||

1.3514 | 4.4982 | 3.3285 | 1.2513 | 2.9491 | 2.3569 | ||

1.9923 | 6.6315 | 3.3285 | 1.5201 | 3.5804 | 2.3553 | ||

0.0648 | 1.2167 | 18.7767 | 0.1939 | 1.2180 | 6.2807 | ||

0.0648 | 1.3270 | 20.4792 | 0.1938 | 1.4080 | 7.2655 | ||

0.0648 | 1.4373 | 22.1818 | 0.1938 | 1.5985 | 8.2487 | ||

0.0648 | 1.5476 | 23.8843 | 0.1938 | 1.7893 | 9.2312 |

In the numerical experiment the range of exponents was limited to and the upper bound yielded condition numbers independent of the discretization parameter, cf. Figure 2. The observation raises a question about the suitablity of , i.e. considering the multiplier space with the norm. It is shown in Remark 3.1 that the choice leads to a condition number with logarithmic growth.

###### Remark 3.1

We consider (11) with . Since is (due to the employed discretization) an identity, the values in (12) are the eigenvalues of the preconditioned Schur complement, where is the preconditioner. We have and observe that the lower bound sums terms that are at most in magnitude. Thus is bounded from below by a constant. On the other hand the upper bound grows as .

Note that for the 2-1 trace and we have, cf. Example 2.1,

while the lower bound as a sum of terms with magnitude decays as . Thus leads to a linearly growing condition number.

The estimates for are confirmed by numerical experiment summarized in Table 2. In particular, the constant lower bound and the upper bound growing as , are visible for both configurations.

Experiments with the spectral discretization suggest that there exists an exponent , independent of , such that the discrete trace operator defined over can be controlled by the -norm (7). However, the space considered thus far consisted of infinitely smooth functions. We proceed to show that a similar conjecture holds if the discrete spaces are obtained by FEM. In particular, the space shall be constructed using the conforming continuous linear Lagrage elements.

### 3.2 Trace operator with FEM discretizaton

Let . Further, let and be, respectively, the basis and degrees of freedom/dual basis nodal with respect to of the finite element space over . The trace mapping shall be defined by interpolation so that is represented in the basis by vector ,

(13) |

Equivalently we have where and the matrix representing the trace operator has entries

where are the basis functions of .

[Discrete trace operator by projection] Let be given and be the projection

Further let be defined via (13). Then is necessary and sufficient for . {proof} To verify the assertion let be the Riesz representation of , i.e. , , and arbitrary. Then by definition and

by the property of the Riesz basis , nodality of the basis and definition of . It follows that . Note that was required to apply the Riesz theorem. {definition}[-matching spaces] Let be a manifold in and , the finite element spaces over the respected domains. The spaces are called -matching if (i) and are constructed from the same elements and (ii) meshes of and are matched.

###### Remark 3.2 (Equivalence of interpolation and projection trace)

The condition from Lemma 3.2 is satisfied with if and are -matching.

Finally, note that the interpolation trace is in general cheaper to construct than the trace due to projection. We shall employ (13) throughout the rest of the paper.

Let now , be a pair of -matching spaces constructed from
continuous linear Lagrange elements. Further, the discretization of the geometry
shall be such that the mesh of is finer at/near than in
the rest of the domain, cf. Table 10 in Appendix B
and Figure 4.
This way the dimensionality of is increased. Finally, we consider the Schur complement^{6}

Figure 3 explores the condition numbers for . It is evident, cf. the zoom-out plot, that for , is not a good preconditioner for the Schur complement. For both configurations there are exponents in that lead to bounded condition numbers. For several values of in this interval, the condition numbers observed on a sequence of uniformly refined meshes are reported in Table 3. Therein can be observed to lead to bounded . Exponent , i.e. the norm, leads to a slight growth in with both and .

L\ | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|

-0.16 | -0.14 | -0.12 | -0.1 | 0 | -0.16 | -0.14 | -0.12 | -0.1 | 0 | |

1 | 4.568 | 4.932 | 5.517 | 6.531 | 19.530 | 5.760 | 6.316 | 7.064 | 8.129 | 24.484 |

2 | 3.883 | 4.282 | 4.804 | 5.545 | 17.525 | 5.743 | 6.300 | 7.085 | 8.175 | 25.253 |

3 | 4.023 | 4.400 | 4.927 | 5.710 | 19.713 | 5.192 | 5.744 | 6.488 | 7.525 | 25.386 |

4 | 4.062 | 4.477 | 5.045 | 5.781 | 21.561 | 5.381 | 5.926 | 6.698 | 7.798 | 28.731 |

We note that in both configurations the behaviour of the eigenvalues is similar to the spectral case. In particular, and grow for , whereas for only grows while is bounded by a constant, see Table 4. Since the extremal eigenvalues are in general not bounded by a constant, is not a discretization of an operator spectrally equivalent to the Schur complement with constants independent of the discretization parameter. However, the relation observed in the experiments

(14) |

suggests existence of a mesh dependent scale in which spectral equivalence can be achieved. In particular, rescaling the -norm matrix as leads to constant bounds, cf. observed constant spectral condition number. We remark that is bounded away from zero for all , in fact the eigenvalue increases with , and in this sense the discrete inf-sup constant never approaches zero.

Based on the mesh-dependent -norm a block-diagonal preconditioner

could be analysed and
shown to be optimal using the results of [32, 33] (see
also the review paper [29]). However, obtaining the scale is
computationally expensive. We shall therefore proceed with (7)
only. In particular, the exponents identified previously shall
be used to construct preconditioners for several 3-1 constrained problems. We note
that the bounds (14) enter estimates for convergence of
iterative solvers, see, e.g., [33], and since the bounds here are
not constant, the proposed preconditioners are theoretically suboptimal. Nevertheless,
the number of iterations in the studied examples will be bounded. We remark
that the smallest and largest eigenvalues are never far from unity in our examples.

L | ||||
---|---|---|---|---|

1 | (0.290, 1.433) | (0.051, 1.000) | (0.207, 1.310) | (0.041, 1.000) |

2 | (0.420, 1.799) | (0.059, 1.040) | (0.256, 1.610) | (0.041, 1.026) |

3 | (0.502, 2.208) | (0.059, 1.161) | (0.342, 1.965) | (0.045, 1.145) |

4 | (0.603, 2.701) | (0.059, 1.276) | (0.401, 2.379) | (0.044, 1.265) |

## 4 Trace coupled problems

The previous experiments revealed a range of exponents for which matrices behaved similarly to the Schur complement, in terms of stability of the condition number, of the related generalized eigenvalue problem. To simplify the discussion, we shall in the following employ . The exponent shall be used to construct preconditioners for two model 3-1 coupled problems.

### 4.1 Babuška’s problem

Let , be a pair of -matching spaces constructed by continuous linear Lagrange elements and consider the problem: Find , such that

(15) | |||||

The system (15) is a Lagrange multiplier formulation of the minimization problem for , and the constraint on . We note that the problem is considered with homogeneous Neumann boundary conditions. A similar problem with and was first studied in [6] to introduce Lagrange multipliers as means of prescribing boundary data.

Similar to the Schur complement study in Section 3.2, the problem shall be considered with two different curves . Moreover, for each configuration we consider three different sequences of uniformly refined meshes, to investigate numerically whether the construction of the preconditioner relies on a quasi-uniform mesh, or if shape-regular elements are sufficient. In a uniform discretization the characteristic mesh size of and are identical and the tessellation of is structured. In finer and coarser discretizations the mesh is unstructured and is either finer or coarser near than in the rest of the domain. The example meshes are pictured in Figure 4. Information about the parameters of the discretizations and sizes of the corresponding finite element spaces are then summarized in Table 10.

Since (15) is considered with Neumann boundary conditions, the block diagonal preconditioner for the system shall have the multiplier block based on (not ). We propose the following preconditioned linear system

(16) |

where and are, respectively, the mass matrices of and . We remark that the proposed preconditioner is not theoretically optimal because of the estimate (14).

In our implementation the leading block of the preconditioner is inverted by
algebraic multigrid from the Hypre^{7}

The recorded iterations counts are reported in Table 5. It can be seen that the proposed preconditioner results in a bounded number of iterations for all the considered geometrical configurations and their discretizations. In the table we also report iteration counts for the preconditioner that employs for the multiplier block. Recall that with and spectral discretization, the spectral condition number of the preconditioned Schur complement showed a logarithmic growth, cf. Table 2. Using FEM, the growth was less evident (see Table 3), however, the condition number was significantly larger than for . The iteration counts agreee with this observation; the norm leads to at least 20 more iterations. We remark that the norms in which the convergence criterion is measured differ between the two cases.

L | ||||||
---|---|---|---|---|---|---|

uniform | finer | coarser | uniform | finer | coarser | |

2 | (28, 59) | (53, 81) | (44, 46) | (29, 57) | (73, 107) | (62, 71) |

3 | (27, 68) | (52, 82) | (49, 58) | (27, 59) | (69, 103) | (64, 81) |

4 | (25, 70) | (52, 83) | (47, 62) | (25, 61) | (69, 105) | (67, 88) |

5 | (23, 70) | (53, 83) | (51, 71) | (25, 62) | (70, 105) | (67, 91) |