Preconditioning trace coupled 3d-1d systems using fractional Laplacian
aa r X i v : . [ m a t h . NA ] A p r PRECONDITIONING TRACE COUPLED 3 D -1 D SYSTEMS USINGFRACTIONAL LAPLACIAN ∗ MIROSLAV KUCHTA † , KENT-ANDRE MARDAL †‡ , AND
MIKAEL MORTENSEN † Abstract.
Multiscale or multiphysics problems often involve coupling of partial differentialequations posed on domains of different dimensionality. In this work we consider a simplified modelproblem of a 3 d -1 d coupling and the main objective is to construct algorithms that may utilize stan-dard multilevel algorithms for the 3 d domain, which has the dominating computational complexity.Preconditioning for a system of two elliptic problems posed, respectively, in a three dimensionaldomain 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 thatnegative fractional Sobolev norms are suitable preconditioners for the Schur complement of the sys-tem. The norms are employed to construct a robust block diagonal preconditioner for the coupledproblem. Key words. preconditioning, saddle-point problem, Lagrange multipliers, trace
AMS subject classifications.
1. Introduction.
Let Ω be a bounded domain in 3 d , while Γ represents a 1 d structure inside Ω, and consider the following coupled problem − ∆ u + u + pδ Γ = f in Ω , (1.1a) − ∆ v + v − p = g on Γ , (1.1b) T u − v = h on Γ . (1.1c)Here the term pδ Γ is to be understood as a Dirac measure such that R Ω p ( x ) δ Γ w ( x ) d x = R Γ p ( t ) w ( t ) d t for a continuous function w . We remark that from a mathematical pointof view the trace T of u required in (1.1c) is in the continuous case not well-definedunless the functions are sufficiently regular. For simplicity of implementation thesystem shall be considered with homogeneous Neumann boundary conditions.The system (1.1) is relevant in numerous biological applications where the em-bedded (three dimensional) structure is such that order reduction techniques can beused to capture its response by a one dimensional model. Equation (1.1a) then modelsprocesses in the bulk, while (1.1c) is the coupling between the domains. A typicalexample of such a system is a vascular network surrounded by a tissue and the orderreduction is due to the employed assumption of radii of the arteries being negligiblein comparison to their lengths. To list a few concrete applications, the 3 d -1 d modelshave been used, e.g., in [18, 25, 17, 33] to study blood and oxygen transport in thebrain or in [11] to describe fluid exchange between microcirculation and tissue inter-stitium. Efficiency of cancer therapies delivered through microcirculation was studiedin [10], and hyperthermia as a cancer treatment in [31]. We note that the employedmodels are more involved than (1.1), but that the system still qualifies as a relevantmodel problem.Due to the Dirac measure term and the three-to-one dimensional trace operator,the problem (1.1) is not standard and establishing its well-posedness is a delicateissue. In fact, considering (1.1a) with a known p and homogeneous Dirichlet boundary † Department of Mathematics, Division of Mechanics, University of Oslo { mirok, kent-and,mikaem } @math.uio.no ‡ Center for Biomedical Computing, Simula Research Laboratory, ∗ Submitted to Numerical Methods for Partial Differential Equations1
Preconditioning for 3 d -1 d coupled problems conditions, the equation is not solvable in H (Ω), as ∇ u may be unbounded in theneighborhood of Γ. A similar problem was studied in [14], where two elliptic problemswere coupled via a Dirac measure source term, and a unique weak solution was foundusing weighted Sobolev spaces. In particular, the weighted spaces that include adistance function ensured that the trace could be defined as a bounded operator. Acorresponding finite element method (FEM) for the problem was discussed in [13],where optimal convergence in the weighted Sobolev norm was shown using gradedmeshes. Optimal convergence of FEM with regular meshes is proved in [22] and [21]for the elliptic problems with point singular data and line singular data and the L norm outside of the fixed neighborhood of the singularity. Therein, the existence ofthe weak solution relies on spaces W ,p (Ω), 1 ≤ p < W ,p (Ω) in [22, 21] were introduced in the analysis of the continuous problems, how-ever, standard finite elements were used in the implementation.While standard finite element methods provide accurate discretization in the alter-native norms of the above mentioned nonstandard Sobolev spaces, this does not implythat standard preconditioning algorithms will be efficient. In fact, as described e.g. in[29], the construction of preconditioners is deeply connected with the mapping proper-ties of the underlying continuous differential operators and to the authors knowledgethe efficiency in the Sobolev spaces with distance functions have not been analyzed.Hence, the use of the weighted Sobolev spaces has prevented the construction of ef-ficient solution algorithms and the more application oriented works [10, 11, 31], thatbuild on the analysis in [13, 14], relied on incomplete LU preconditioning. To resolvethis problem, we have in this paper taken an alternative approach where standardmultilevel algorithms for elliptic problems are reused for the 3 d problem. This ap-proach does however require that novel algorithms are developed for the 1 d problem.The special construction of algorithms for the 1 d problem is justified by the fact that,in general, the computational complexity of a 1 d problem is low compared to a 3 d problem. We shall illustrate this fact by several numerical experiments.The current paper is an extension of [23], where a system similar to (1.1a)–(1.1c)was analyzed for the case Ω a bounded domain in 2 d and Γ a structure of codimensionone. Therein, robust preconditioners were established, based on the operator precon-ditioning framework [29], in which preconditioners are constructed as approximateRiesz mappings in properly chosen Hilbert spaces. The framework often allows forconstruction of order-optimal preconditioners, with convergence independent of mate-rial and discretization parameters, directly from the analysis of the continuous systemof equations. In particular, in [23] it was shown that the proper preconditioning reliedon a nonstandard fractional H − inner product. Crucial for the analysis was the factthat the trace operator T is a well-defined mapping between H (Ω) and H (Γ), whenΓ is of codimension one with respect to Ω.The case when the trace operator T maps functions defined on Ω to Γ and Γis of codimension two is challenging as the properties of the trace operator are notestablished from a theoretical point of view. If we assume some additional regularitysuch that u ∈ H ǫ then T u ∈ H ǫ for ǫ >
0, see e.g. [15]. However, the resultis known to break down in the limit when ǫ = 0. We therefore propose to weakenthe requirements on T and instead consider T as a mapping between H (Ω) and H s (Γ) for some s < s we perform a comprehensive numerical study with various discretizations; that is,finite element methods with conforming/non-conforming elements and matching/non-matching meshes, and by considering the Galerkin method with eigenbasis of Laplaceoperator. We demonstrate that all of these different methods point to the constructionof the same preconditioning operator, namely ( − ∆) s , where s ∈ ( − . , − .
1) andthe range seems to be independent of the discretization method. We demonstratenumerically that this choice defines a good preconditioner for problems with complex1 d geometries and for 3 d meshes that are both highly refined or rather coarse close tothe 1 d mesh as long as the mesh is shape-regular and the discrete problem is invertible.Our work is structured as follows. In § d -1 d trace operator. In § d -1 d problems discretized with FEM and matched discretizations of Ω and Γ. In § §
2. Notation and preliminaries.
Let X be a Hilbert space of functions definedon a domain D ⊂ R d , d = 1 , ,
3. The norm of the space is denoted by k·k X , while h· , ·i X ′ ,X is the duality pairing between X and its dual space X ′ . We let ( · , · ) X denote the inner product of X , while, to simplify the notation, ( · , · ) D is the L innerproduct. The Sobolev space of functions with m square integrable derivatives is H m ( D ). Finally, H m ( D ) denotes the closure of the space of smooth functions withcompact support in D in the H m ( D ) norm. We will also employ Sobolev spaces withfractional derivatives, which are more precisely defined later.We use normal capital font to denote operators over infinite dimensional spaces,e.g. A : X → X ′ . If A : X → Y is a bounded operator we let A ′ : Y ′ → X ′ denote theadjoint operator h y ′ , Ax i Y ′ ,Y = h A ′ y, x i X ′ ,X , y ∈ Y ′ , x ∈ X . For a discrete subspace X h ⊂ X , dim X h = n , the subscript h is used to distinguish the finite dimensionaloperator due to the Galerkin method, e.g., A h : X h → X ′ h defined by h A h u h , v h i X ′ ,X = h Au, v h i X ′ ,X u h , v h ∈ X h and u ∈ X. For a given basis, { φ i } ni =1 of X h , the matrix representation of the operator is denotedby sans serif font. Thus A h is represented by A ∈ R n × n with entries A i,j = h A h φ j , φ i i X ′ ,X . The function u h ∈ X h is represented in the basis by a coefficient vector u ∈ R n , where u h = u i φ i (summation convention invoked). Finally, for the inner product of vectors u , v in R n shall be denoted as u ⊤ v . We consider Ω ⊂ R d an open con-nected domain with Lipschitz boundary ∂ Ω and Γ a Lipschitz submanifold of codi-mension one or two in Ω. The trace operator T is defined by T u = u | Γ for u ∈ C (Ω).In case the codimension of Γ is one, the properties of the trace operator are wellknown. In particular, T : H s (Ω) → H s − (Γ) is bounded and surjective, see, e.g.,[1, ch. 7] for s > . As a direct consequence we then have that the trace to Γ ofcodimension two is well behaved as mapping from H ǫ (Ω) to H ǫ (Γ) for any ǫ >
0, cf.[34] and [15] for the case of unbounded and bounded domains respectively. However,for ǫ = 0, it is known that the trace operator is unbounded as a mapping between H (Ω) and L (Γ). We therefore conjecture that the trace is well-behaved between H (Ω) and H s (Γ), s < Preconditioning for 3 d -1 d coupled problems The fractional Sobolev space H s (Γ) shall be defined by interpolation [26, 5]. Forthe sake of completeness, we review here the presentation from [23]. Let u, v ∈ X = H (Γ). For u fixed v ( u, v ) Γ is in X ′ and by the Riesz-Fr´echet theorem there is aunique w ∈ X such that ( w, v ) X = ( u, v ) Γ for any v ∈ X . The operator S : u → w is injective and compact and thus the eigenvalue problem Sφ i = λ i φ i (no summationimplied) is well-defined. In addition, S is self-adjoint and positive-definite such thatthe eigenvalues form a nonincreasing sequence 0 < λ k +1 ≤ λ k and λ k →
0. Bydefinition, the eigenvectors satisfy( φ i , v ) X = λ − i ( φ i , v ) Γ v ∈ X, or equivalently Aφ i = λ − i M φ i with h Au, v i X ′ ,X = ( u, v ) X and h M u, v i X ′ ,X = ( u, v ) Γ . (2.1)Further, the set of eigenvectors { φ k } ∞ k =1 forms a basis of X , which is orthogonal inthe inner product of X and orthonormal in the L (Γ) inner product. Finally, for s ∈ [ − ,
1] we define the s -norm of u = c k φ k ∈ span { φ k } ∞ k =1 as k u k H s (Γ) = q c k λ − sk . (2.2)The space H s (Γ) is finally defined as the closure of the span { φ k } ∞ k =1 in the s -norm,while H s (Γ) is then defined analogously to H s (Γ) with X = H (Γ) in the construction.Following the approach in [23], a weak formulation of the homogeneous Dirichletproblem for (1.1a)–(1.1c) with Ω ∈ R , Γ ⊂ Ω of codimension one or two , using themethod of Lagrange multipliers, reads: Find ( u, v, p ) ∈ H (Ω) × H (Γ) × Q such that( ∇ u, ∇ φ ) Ω + ( u, φ ) Ω + ( p, T φ ) Γ = ( f, φ ) Ω ∀ φ ∈ H (Ω) , ( ∇ v, ∇ ψ ) Γ + ( v, ψ ) Γ − ( p, ψ ) Γ = ( g, ψ ) Γ ∀ ψ ∈ H (Γ) , ( χ, T u − v ) Γ = ( h, χ ) Γ ∀ χ ∈ Q. (2.3)For Γ of codimension one, the problem is well-posed with Q = H − owing to the factthat T : H (Ω) → H (Γ) is an isomorphism. Similarly, for Γ of codimension two, thewell-posedness hinges on whether T : H (Ω) → Q ′ is an isomorphism for some space Q . Hovewer, to the best of the authors knowledge this result is not known. In thispaper, we therefore conjecture that the space Q is closely related to H s (Γ) for somesuitable s < A defined by (2.3) A uvp = I Ω − ∆ Ω T ′ I Γ − ∆ Γ − I Γ T − I Γ uvp = fgh . (2.4)is an isomorphism mapping H (Ω) × H (Γ) × Q to its dual space and a proper pre-conditioner can be formed as B = ( I Ω − ∆ Ω ) − I Γ − ∆ Γ ) −
00 0 R Q , (2.5)where R Q is the Riesz mapping between the dual of Q and Q , cf. [29]. As the Rieszmapping is not easily obtained from the analysis of the continuous problem, we shallin the following resort to investigating the mapping properties of the trace operatorof codimension two by a series of numerical experiments with different spaces for Q .Let now V h ⊂ H (Ω). Considering (1.1) on the finite dimensional spaces, weobtain a variational problem: Find ( u h , v h , p h ) ∈ V h × W h × Q h such that( ∇ u h , ∇ φ h ) Ω + ( u h , φ h ) Ω + ( p h , T h φ h ) Γ = ( f, φ h ) Ω φ h ∈ V h , ( ∇ v h , ∇ ψ h ) Γ + ( v h , ψ h ) Γ − ( p h , ψ h ) Γ = ( g, ψ h ) Γ ψ h ∈ W h , h χ h , T h u h − v h i Γ = ( h, χ h ) Γ χ h ∈ Q h . (2.6)Here the discrete trace operator T h is well-defined as the functions in V h are con-tinuous. In the absence of existence result for the continuous problem the discretepreconditioner cannot be constructed within the framework of operator precodition-ing, i.e. as a discretization of a suitable Riesz mapping. We therefore adapt a differentframework, namely the matrix Schur complement preconditining [8, 30]. That is, weattempt to construct the preconditioner for (2.6) by reasoning directly about theproperties of the discrete systems.From a linear algebra point of view, the problem (2.6) is a saddle-point system (cid:20) A B ⊤ B (cid:21) (cid:20) xy (cid:21) = (cid:20) bc (cid:21) , with A a symmetric positive definite matrix. In case B has a full row rank, the discreteproblem is uniquely solvable and block diagonal preconditioner can be constructedas an approximate inverse of the matrix diag( K , L ), where K should be spectrallyequivalent with A and L should be spectrally equivalent with the Schur complement BA − B ⊤ , see, e.g., [35, 36]. Considering (2.6), the key question is thus whether itis possible (in an efficient and systematic manner) to construct an operator that isspectrally equivalent with the Schur complement. Motivated by the 2 d -1 d problemand our conjectured mapping properties of the trace the operator shall be based onthe norm of the H s (Γ) space (2.2).Following [23], the discrete approximation of the s -norm shall be constructed bymirroring the continuous eigenvalue problem (2.1). More specifically, let X h ⊂ X and matrices A , M be the representations of A h , M h ; the Galerkin approximations ofoperators A , M from (2.1). Then there exists an invertible matrix U and diagonal,positive-definite matrix Λ satisfying AU = MUΛ . Moreover, the product U ⊤ MU is anidentity such that the columns of U form an A orthogonal and M orthonormal basis of R n . In order to define the discrete norm, we let H s be a symmetric, positive-definitematrix H s = ( MU ) ⊤ Λ s ( MU ) . (2.7)The matrices H s, are defined analogously to (2.7), using the eigenvalue problem forthe Laplace operator with homogeneous Dirichlet boundary conditions. For u h ∈ X h represented in the basis of the space by a coefficient vector u , let c be the representationof u in the basis of eigenvectors, that is, u = Uc . We then set k u h k H s (Γ) = p u ⊤ H s u = √ c ⊤ Λ s c . (2.8)
3. Norms for the discrete 3 d -1 d trace. The matrices H s shall be employed toconstruct a preconditioner for the Schur complement of the system (2.6). Considering(2.3), the matrix is a sum of two parts which correspond respectively to operators Preconditioning for 3 d -1 d coupled problems T ( − ∆ Ω + I Ω ) − T ′ and I Γ ( − ∆ Γ + I Γ ) − I ′ Γ . As the matrix stemming from the latterterm is by definition spectrally equivalent with H − we shall next focus only on theformer trace term. We note that if our conjucture on the mapping properties ofthe trace operator holds, that is T : H (Ω) → H s (Γ) is bounded and sujective forsome s <
0, then T ( − ∆ Ω + I Ω ) − T ′ : H s (Γ) ′ → H s (Γ) is an isomorphism and thepreconditioner could be realized by the fractional norm matrix.To investiagate the conjectured spectral equivalence of the trace term, let V, Q bethe spaces of continuous functions over Ω and Γ respectively and consider the problemof minimizing v ( ∇ v, ∇ v ) Ω − f, v ) Ω , v ∈ V , subject to v = 0 on the boundaryand the constraint T v = g on Γ. The minimization problem leads to the variationalproblem for u ∈ V , p ∈ Q satisfying( ∇ u, ∇ v ) Ω + ( p, T v ) Γ = ( f, v ) Ω ∀ v ∈ V, ( q, T u ) Γ = ( q, g ) Γ ∀ q ∈ Q. (3.1)The Schur complement of (3.1) is thus closely related to the critical trace term in theSchur complement of (2.3).Using finite dimensional subspaces of V h ⊂ V and Q h ⊂ Q the problem (3.1) isequivalent to the linear system (cid:20) A T ⊤ T (cid:21) (cid:20) up (cid:21) = (cid:20) fg (cid:21) . (3.2)and we wish to find computational evidence for the following claim. Conjecture 3.1.
There exist s < and constants < λ ∗ ≤ λ ∗ such that forany h > λ ∗ ≤ x ⊤ ( TA − T ⊤ ) xx ⊤ H s, x ≤ λ ∗ . (3.3)In addition to spectral equivalence condition (3.3) we shall also consider a weakerrequirement, where we wish to find s for which the condition number of the precon-ditioned Schur complement is bounded in h, H for some s <
0. More precisely, let0 < λ min ( s, h ) ≤ λ max ( s, h ) be the smallest and largest eigenvalues of the generalizedeigenvalue problem ( TA − T ⊤ ) p = λ H s p . (3.4) Conjecture 3.2.
There exist s < such that the condition number κ ( s, h ) = λ max ( s, h ) λ min ( s, h ) ≤ C ∀ h > , (3.5) for some constant C . We note that the condition (3.5) is motivated by the fact that convergence ofthe preconditioned conjugate gradient method is estimated in terms of the conditionnumber, see, e.g., [38]. For suitable s the linear system with the Schur complementcould thus be solved efficiently. We also note that the condition is weaker than spectralequivalence (3.3). We use the same subscript to signify that the function spaces cannot be arbitrary and insteadmust satisfy inf-sup compatibility condition.
To investigate conjectures 3.1, 3.2 we let Ω = [0 , and choose Γ as simplestraight lines; Γ = { ( t, , ); t ∈ [0 , } and Γ = { ( t, t, t ); t ∈ [0 , } . For dis-cretization of (3.1) the discrete subspaces shall be first constructed using the basis ofeigenfunctions of the Laplacian. Let { φ k } k ≥ be the setof eigenvectors of the Laplace operator on unit interval with homogeneous Dirichletboundary conditions and set Q m = span { φ k } mk =1 while the n dimensional space V n of functions on Ω shall be defined as a tensor product.Considering (3.1) with spaces V n , Q m the matrix A in (3.2) diagonal. The tracematrix T ∈ R m × n for curve Γ is sparse with entries T j, ( i,k,l ) = ( k or l even( − k +1 ( − l +1 δ ij otherwise . Note that for m > n the matrix does not have a full row rank and the system issingular. We therefore set m = n . For Γ the trace matrix is sparse with a moreinvolved sparsity pattern and at most four nonzero entries per row T j, ( i,k,l ) = 4 √ Z sin jπt sin iπt sin kπt sin lπt d t. Having defined the terms in (3.2) we consider the generalized eigenvalue problem(3.4) for different values of s and the discretization parameter n . Observe that in caseof Γ the Schur complement can be obtained in a closed form. Indeed, the matrix isdiagonal S j δ ij (no summation implied) with entries S j = 4 π n X l,m odd 1 j + l + m . (3.6)For Γ the matrix is dense and shall be computed from the definition TA − T ⊤ . Assuch a smaller n is explored in this configuration.The results of the numerical experiments with s ∈ [ − . , − .
1] are summarizedin Figure 3.1. We observe that values s ∈ [ − . , − .
1] yield bounded conditionnumbers for Γ . The condition numbers are not quite converged for the other config-uration, however, it is possible to identify unstable exponents s < − .
18. Moreover,the values close to s = − .
14 appear to be stable also in this configuration. This factis easier to appreciate in Table 3.1, which shows λ min , λ max and κ as functions of thediscretization parameter for s = − .
14. For Γ the condition number is evidently con-stant, while for Γ the number appears to be bounded. The observation are thereforesupportive of conjecture 3.2.For neither of the configurations and any of the considered vales s the smallestand largest eigenvalues are bounded and thus, contrary to conjecture 3.1, the matrices H s, are not spectrally equivalent with the Schur complement. However, taking e.g. s = − .
14, either of λ min ( n ), λ max ( n ) defines a mesh-depenedent scale τ ( n ) H − . , that yields spectral equivalence. Such scale, however, is not easily computable ingeneral as it involves the inverse of the 3 d problem.In the numerical experiment the range of exponents was limited to s ∈ [ − . , − . Preconditioning for 3 d -1 d coupled problems − . − . − . − . − . − . s κ n = 2 n = 2 n = 2 n = 2 − . − . − . − . − . − . . s κ n = 2 n = 2 n = 2 n = 2 Fig. 3.1: Spectral condition numbers (3.5) computed from the generalized eigenvalueproblem for Schur complement of (3.2) and matrices H s, , see (2.8). (Left) Theconstraint is considered on Γ = { ( t, , ); t ∈ [0 , } . (Right) Γ = { ( t, t, t ); t ∈ [0 , } is considered.Table 3.1: Smallest and largest eigenvalues λ min , λ max and the spectral conditionnumbers κ of (3.4). (Top) The preconditioner is H − . , . While the eigenvalues areunbounded the condition number is bounded in n . (Bottom) Matrix H , (identitymatrix) is used as the preconditioner. In agreement with the analysis in Remark 3.1,constant λ min and λ max with a logarithmic growth are observed. Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } log n λ min λ max κ log n λ min λ max κ
10 0.6218 2.0696 3.3285 6 0.8476 1.9916 2.349612 0.9167 3.0511 3.3285 7 1.0298 2.4283 2.358114 1.3514 4.4982 3.3285 8 1.2513 2.9491 2.356916 1.9923 6.6315 3.3285 9 1.5201 3.5804 2.355311 0.0648 1.2167 18.7767 6 0.1939 1.2180 6.280712 0.0648 1.3270 20.4792 7 0.1938 1.4080 7.265513 0.0648 1.4373 22.1818 8 0.1938 1.5985 8.248714 0.0648 1.5476 23.8843 9 0.1938 1.7893 9.2312 s = 0, i.e. considering the multiplier space Q m with the L norm. It is shown inRemark 3.1 that the choice leads to a condition number with logarithmic growth. Remark 3.1.
We consider (3.2) with Γ . Since H , is (due to the employeddiscretization) an identity, the values S j in (3.6) are the eigenvalues of the precon-ditioned Schur complement, where H , is the preconditioner. We have S j ≥ S n andobserve that the lower bound sums O ( n ) terms that are at most n − in magnitude.Thus S n is bounded from below by a constant. On the other hand the upper bound S j ≤ S grows as log n .The estimates for Ω ⊂ R are confirmed by numerical experiments summarizedin Table 3.1. In particular, the constant lower bound and the upper bound growingproportianaly to log n , are visible for both configurations. Experiments with the spectral discretization suggest that there exists a rangeof negative exponents s , independent of Γ, such that the discrete trace operator T h defined over V h can be controlled by the s -norm (2.8) in the sense of (3.5) andconjecture 3.2. However, the space V h considered thus far consisted of infinitelysmooth functions. We proceed to show that the statement holds if the discrete spacesare obtained by FEM. In particular, the space V h shall be constructed using the H conforming continuous linear Lagrange elements. Let V h ⊂ H (Ω). Further,let { ψ k } mk =1 and { L j } mj =1 be, respectively, the basis and degrees of freedom/dual basisnodal with respect to { ψ k } mk =1 of the finite element space Q h over Γ. The tracemapping T h : V h → Q h shall be defined by interpolation so that p h = T h u h isrepresented in the basis by vector p ∈ R m , p j = h L j , u h | Γ i . (3.7)Equivalently we have p = Tu where u ∈ R n and the matrix representing the traceoperator has entries T i,j = h L i , φ j | Γ i , (3.8)where { φ j } nj =1 are the basis functions of V h . Lemma 3.1 (
Discrete trace operator by projection ). Let u h ∈ V h be given and ˜ p h ∈ Q h be the L projection (˜ p h , q ) Γ = ( u h | Γ , q ) Γ , q ∈ Q h . Further let p h ∈ Q h be defined via (3.7) . Then V h | Γ ⊆ Q h is necessary and sufficientfor p h = ˜ p h .Proof . To verify the assertion let q k ∈ Q h be the Riesz representation of L k , i.e.( q k , v ) Γ = h L k , v i , v ∈ Q h , and u h ∈ V h arbitrary. Then by definition ( p h , q k ) Γ = h L i , u h | Γ i ( ψ i , q k ) Γ and h L i , u h | Γ i ( ψ i , q k ) Γ = ( q i , u h | Γ ) Γ h L k , ψ i i = ( q k , u h | Γ ) Γ = ( q k , ˜ p h ) Γ by the property of the Riesz basis { q k } mk =1 , nodality of the basis { ψ i } mi =1 and definitionof ˜ p h . It follows that ( p h − ˜ p h , q k ) Γ = 0. Note that u h | Γ ∈ Q h was required to applythe Riesz theorem. The above result ensures that T ⊤ has full column rank, andconsequently the matrix TAT ⊤− is non-singular. Definition 3.2 (
Γ-matching spaces ). Let Γ be a manifold in Ω and Q h , V h thefinite element spaces over the respective domains. The spaces are called Γ -matchingif (i) V h and Q h are constructed from the same elements and (ii) meshes of Ω and Γ are matched. Remark 3.2 (
Equivalence of interpolation and projection trace ). The conditionfrom Lemma 3.1 is satisfied with V h | Γ = Q h if V h and Q h are Γ -matching. Finally,note that the interpolation trace is in general cheaper to construct than the trace dueto projection. We shall employ (3.7) throughout the rest of the paper. Consequentlythe trace matrix T in (3.2) is a product of the mass matrix of the space Q h and (3.8).Let now V h , Q h be a pair of Γ-matching spaces constructed from continuous linearLagrange elements. Further, the discretization of the geometry shall be such that themesh of Ω is finer at/near Γ than in the rest of the domain, cf. Table A.1 in AppendixA and Figure 4.1. This way the dimensionality of Q h is increased. Finally, we considerthe Schur complement of (3.1) preconditioned by different matrices H s, . Recall that The Schur comeplement is computed from its definition, where the components T , A areassembled using FEniCS [27, 2] and PETSc [7] libraries. The Laplacian matrix is then inverted byconjugate gradient method with algebraic multigrid (AMG) preconditioner from Hypre library [16].Relative tolerance 10 − was set as a convergence criterion. Preconditioning for 3 d -1 d coupled problems − . − . − . − .
05 005101520 s κ n = 161 n = 321 n = 641 n = 1281 − . − .
25 00100200300400500 − . − . − . − .
05 00102030 s κ n = 187 n = 373 n = 742 n = 1481 − . − .
25 0050100150200250300
Fig. 3.2: Condition numbers (3.5) of (3.4) with finite element discretization, n =dim Q h , and different preconditioners H s, . (Left) the curve is Γ . (Right) the curveis Γ . The zoomed out plot shows that s < − .
25 yields unbounded κ . For bothconfigurations exponents from the interval around s = − . s . The finite ele-ment discretization is considered on a sequence of uniformly refined meshes, see TableA.1. For each discretization the mesh is finer near the curve than in the rest of thedomain. Exponent s = − .
14 observed in the spectral discretization, cf. Table 3.1,yields bounded κ also with discrization by FEM. Note that similar to the spectraldiscretization there is a slight growth of κ for s = 0. L \ s Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } -0.16 -0.14 -0.12 -0.1 0 -0.16 -0.14 -0.12 -0.1 01 4.568 4.932 5.517 6.531 19.530 5.760 6.316 7.064 8.129 24.4842 3.883 4.282 4.804 5.545 17.525 5.743 6.300 7.085 8.175 25.2533 4.023 4.400 4.927 5.710 19.713 5.192 5.744 6.488 7.525 25.3864 4.062 4.477 5.045 5.781 21.561 5.381 5.926 6.698 7.798 28.731 previously global trigonometric polynomial basis functions were used with (3.1) and − . < s ≤ − . s ∈ [ − . , s < − . H s, is not a good preconditioner for the Schurcomplement. For both configurations there are exponents in ( − . , − .
1) that leadto bounded condition numbers. For several values of s in this interval, the conditionnumbers observed on a sequence of uniformly refined meshes are reported in Table3.2. Therein s ≤ − . κ . Exponent s = 0, i.e.the L norm, leads to a slight growth in κ with both Γ and Γ .We note that in both configurations the behaviour of the eigenvalues is similarto the spectral case. In particular, λ max and λ min grow for s ≤ − .
1, whereas for s = 0 only λ max grows while λ min is bounded by a constant, see Table 3.3. Sincethe extremal eigenvalues are in general unbounded H s, is not a discretization of an1Table 3.3: Smallest and largest eigenvalues of the H s, preconditioned Schur comple-ment considered in Table 3.2. Similar to spectral discretization both the extremaleigenvalues grow for s = − .
14 while the lower bound is constant and the upper onegrows for s = 0. L Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } s = − . s = 0 s = − . s = 01 (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) operator spectrally equivalent to the Schur complement and, similar to § < λ min ( s, h ) ≤ x ⊤ TA − T ⊤ xx ⊤ H s, x ≤ λ max ( s, h ) ∀ x ∈ R m (3.9)suggests existence of a mesh dependent scale in which spectral equivalence can beachieved, cf. also results of § s -norm matrix as λ min ( s, h ) H s, leads to constant bounds, cf. observed constant spectral conditionnumber. We remark that λ s, min is bounded away from zero for all h and s observed,in fact the eigenvalue increases with h − , and in this sense the discrete inf-sup constantnever approaches zero.Computational results with the spectral basis and FEM both suggest to the con-struction of the Schur complement preconditioner based on the mesh dependent s -norm λ min ( s, h ) H s, . However, as noted before, obtaining the scaling factor is compu-tationally expensive and we shall therefore proceed with (2.8) only and not include thescale. In particular, the exponents s identified previously shall be used to constructpreconditioners for several 3 d -1 d constrained problems. We note that the bounds(3.9) enter estimates for convergence of iterative solvers, see, e.g., [36], and since thebounds here are not constant, the proposed preconditioners are theoretically subopti-mal. 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 ourexamples.
4. Trace coupled problems.
The previous experiments revealed a range ofnegative exponents s for which matrices H s behaved similarly to the Schur com-plement, in terms of stability of the condition number, of the related generalizedeigenvalue problem. To simplify the discussion, we pick s = − .
14 and employ theexponent to construct preconditioners for two model 3 d -1 d coupled problems. We notethat this choice is somewhat arbitrary and based on § s = − . Let V h , Q h be a pair of Γ-matching spaces con-structed by continuous linear Lagrange elements and consider the problem: Find u ∈ V h ⊂ H (Ω), p ∈ Q h ⊂ H (Γ) such that( ∇ u, ∇ v ) Ω + ( u, v ) Ω + ( p, T v ) Γ = ( f, v ) Ω v ∈ V h , ( q, T u ) Γ = ( q, g ) Γ q ∈ Q h . (4.1)The system (4.1) is a Lagrange multiplier formulation of the minimization problemfor v
7→ k v k H (Ω) − f, v ) Ω , with the constraint T v − g = 0 on Γ. The problem2 Preconditioning for 3 d -1 d coupled problems Fig. 4.1: Domains used in experiments with matching discretization. The one di-mensional curve Γ is drawn in blue with element boundaries signified by red dots.(Left) The curve is, respectively, a horizontal or diagonal segment. The triangulationof Ω is either refined or coarsened at Γ. (Right) The curve contains branches andbifurcations, thus capturing some of the features of complex vascular systems.is considered with homogeneous Neumann boundary conditions. A similar problemwith Ω ⊂ R and Γ ⊂ ∂ Ω was first studied in [6] to introduce Lagrange multipliers asmeans of prescribing boundary data.Similar to the Schur complement study in Section 3.2, the problem shall be consid-ered with two different curves Γ. Moreover, for each configuration we consider threedifferent sequences of uniformly refined meshes, to investigate numerically whetherthe construction of the preconditioner relies on a quasi-uniform mesh, or if shape-regular elements are sufficient. In a uniform discretization the characteristic meshsize 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.1. In-formation about the parameters of the discretizations and sizes of the correspondingfinite element spaces are then summarized in Table A.1.Since (4.1) is considered with Neumann boundary conditions, the block diagonalpreconditioner for the system shall have the multiplier block based on H s (not H s, ).We propose the following preconditioned linear system (cid:20) A + M H − . (cid:21) − (cid:20) A + M ( M Γ T ) ⊤ ( M Γ T ) (cid:21) (cid:20) up (cid:21) = (cid:20) A + M H − . (cid:21) − (cid:20) fg (cid:21) , (4.2)where M and M Γ are, respectively, the mass matrices of V h and Q h . We remark thatthe proposed preconditioner is not theoretically optimal because of the estimate (3.9).In our implementation the leading block of the preconditioner is realized by a3Table 4.1: Iteration counts for preconditioned Babuˇska’s problem (4.1) with precon-ditioners based on (2.7) and s = − .
14 or s = 0 (discrete L norm). Two geometricconfigurations and their different discretizations ( L denotes the refinement level) areconsidered cf. Figure 4.1 and Table A.1. Both preconditioners yield bounded numberof iterations. The L norm leads to a less efficient preconditioner. L Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } uniform finer coarser uniform finer coarser2 (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) single V cycle of algebraic multigrid from the Hypre library [16]. The system isthen solved iteratively with the minimal residual method (MINRES) implemented incbc.block [28] and requiring a preconditioned residual norm smaller than 10 − forconvergence. The initial vectors were random.The recorded iterations counts are reported in Table 4.1. It can be seen thatthe proposed preconditioner results in a bounded number of iterations for all theconsidered geometrical configurations and their discretizations. In the table we alsoreport iteration counts for the preconditioner that employs H = M Γ for the multiplierblock. Recall that with s = 0 and spectral discretization, the spectral conditionnumber of the preconditioned Schur complement showed a logarithmic growth, cf.Table 3.1. Using FEM, the growth was less evident (see Table 3.2), however, thecondition number was significantly larger than for s = − .
14. The iteration countsagreee with this observation; the L norm leads to at least 20 more iterations. Weremark that the norms in which the convergence criterion is measured differ betweenthe two cases. Building upon the Babuˇska problem wenext consider a model multiphysics problem (1.1). A similar problem with Ω ⊂ R and Γ a manifold of codimension one was previously studied by the authors in [23].Therein it was found that the problem is well posed with the Lagrange multiplier inthe intersection space H − (Γ) ∩ H − (Γ). The structure of the space was mirrored bythe preconditioner, which used ( H − . + H − ) − in the corresponding block.We note that the exponent − was dictated by the properties of the continuoustrace operator. In the 3 d -1 d case, which is of interest here, we shall instead base theexponent/preconditioner on the previous numerical experiments. More specifically,the linear system obtained by considering (2.6) on finite dimensional finite elementsubspaces A Ω + M Ω ( M Γ T ) ⊤ A Γ + M Γ − M Γ ( M Γ T ) − M Γ uwp = fgh (4.3)shall be considered with the preconditioner A Ω + M Ω A Γ + M Γ H − . + H − − . (4.4) We have used default values of all the parameters. Preconditioning for 3 d -1 d coupled problems Table 4.2: Iteration counts for the model problem (4.3) with preconditioner (4.4).Spatial configurations and disretizations from Table 4.1 are considered. In all thecases the number of iterations is bounded.
L Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } uniform finer coarser uniform finer coarser2 51 45 42 44 62 623 49 45 48 43 59 624 47 43 47 43 59 645 46 43 49 42 59 66 Note that in (4.4) the structure of the trailing block mimics the related 2 d -1 d problem.We remark that in the implementation, the remaining two blocks are realized by AMG.Moreover the discrete spaces are such that W h = Q h and V h , Q h are Γ-matching. As inthe previous example, continuous linear Lagrange elements are used. To demonstratethe performance of the preconditioner, (2.6) is considered on the same geometricalconfigurations and their discretizations as (4.1). The preconditioned system is thensolved by MINRES, starting from a random initial vector and terminating if thepreconditioned residue is less than 10 − in magnitude. As can be seen in Table 4.2,the preconditioner yields bounded iteration counts. Interestingly, the convergence isfaster on the finer discretization than on the coarser one. We note that the systemson the latter discretization are in general of smaller size and have more than a factor10 fewer degrees of freedom in Q h . However, dim Q h ≪ dim V h is a desirable featureof the model order reduction which was applied to obtain the problem on Γ.In the examples presented thus far, Γ was always a straight segment. To showthat the preconditioner (4.4) (or the general idea of H s based preconditioners for 3 d -1 d problems) is not limited to such simple curves, we shall in the final example consider(2.6) with Γ having a more complicated stucture. The considered domain, pictured inthe right pane of Figure 4.1, is inspired by biomechnical applications and is intendedto mimic some of the features of the vasculature. In particular, the domain consistsof numerous branches and contains multiple bifurcations.Repeating the setup of the previous experiment, Table 4.3 reports the iterationcounts for the (4.4) preconditioned linear system (4.3), obtained by considering (2.6)on the complex Γ. The number of iterations is clearly bounded.The good performance of the proposed preconditioner in all the considered ex-amples brings in the question of practicability of its construction. Here, the questionshall be addressed by considering the setup costs of the preconditioner for the domainwith complex Γ. The choice is motivated by the fact that (i) the domain is potentiallyrelevant for practical applications and (ii) the large (relative to dim V h ) number ofdegrees of freedom of Q h puts the emphasis on the construction of (2.7). We notethat the costs are expected to be determined by the multigrid setup and the solutiontime of the generalized eigenvalue problem (2.7). As in [23] the eigenvalue problem issolved by the DSYGVD routine from LAPACK [3].The timings obtained on a Linux machine with a single Intel Xeon E5-2680 CPUwith 2.5GHz and 32GB of RAM are reported in Table 4.3. The observed costs of theeigenvalue solve are 3-4 times smaller than that of the multigrid setup, and thus thespectral construction does not present a bottleneck. Morover, both AMG and GEVPare expected to scale roughly as dim Q h . However, due to the cubic scaling, thesystem/preconditioner is unlikely to be assembled/setup in serial. For such a case,5Table 4.3: Iteration counts and setup costs (in seconds) for system (4.3) and precon-ditioner (4.4). Both operators are assembled for the complex Γ pictured in Figure4.1. The number of iterations is bounded in the discretization parameter. In the con-sidered example, the eigenvalue (GEVP) based construction (2.7) does not present abottleneck as it is 3-4 times cheaper than setting up the algebraic multigrid (AMG). dim V h dim Q h s ] GEVP [ s ]18K 817 86 0.2 0.1100K 1605 81 1.9 0.6634K 3193 76 15.0 4.24.8M 6381 68 141.6 36.4 a scalable parallel implementation, for the construction of (2.7), remains an issue,and approaches that provide the approximate action of H s matrices may offer betterperformance. Examples of such approaches are the [5, 4] and [20] where polynomialand rational function approximations are constructed, fast Fourier transforms [32] ormethods [19, 9] based on integral definitions of fractional Laplacian [24].
5. Nonmatching discrete trace.
The numerical examples presented thus farhave always employed Γ-matching finite element spaces. We note that in [23] thisconstruction is shown to imply that the discrete inf-sup condition holds for problems(4.1) and (2.6) considered with Ω ⊂ R and Γ a one dimensional curve. However,the assumption of matched discretizations of Ω and Γ can be too limiting, e.g, if fineresolution is requested on the curve. In this section we present numerical examplesusing the Babuˇska problem (4.1), which demonstrate that the matching discretiza-tion assumption is not necessary and to the extent given by the new inf-sup conditionthe discretizations can be independent. Using such stable discretizations and pre-conditioners based on characterization of the trace the observed number of Kryloviterations will remain bounded. Consider (4.1) with Ω ⊂ R . For Γ ⊂ ∂ Ω, the finiteelement discretization of the problem requires that the spaces V h , Q H (we use differentsubscripts to indicate the difference in underlying triangulations) are such that h ≤ cH for some c <
1. Here h is understood as a mesh size of V h on Γ. The inequality ensuresthat the discrete inf-sup condition is satisfied [37, 12].Let now Γ be a curve, contained in Ω, where the domains are discretized suchthat the condition from the previous paragraph is met. Further, the space V h shallbe discretized by continuous linear Lagrange elements, while, for the construction of Q H , either the same elements or piecewise constant Lagrange elements are employed.We note that with the latter choice the eigenvalue problem for the discrete s -normsimplifies, since the mass matrix is diagonal in this case.Table 5.1 reports the number of MINRES iterations on the system (4.1), usingdiag(AMG( A + M ) , H − . − ) as the preconditioner. The iterations are started froma random vector using 10 − as the stopping tolerance for the magnitude of thepreconditioned residual. With both considered finite element discretizations of themultiplier space the number of iterations is bounded indicating (i) that the inf-supcondition is satisfied and (ii) the optimality of the preconditioner. We note that for h > H , the inf-sup condition is violated and in turn the the iterations are unbounded(not reported here). An example of a pair of inf-sup stable and unstable discretizationsis shown in Figure 5.1.6 Preconditioning for 3 d -1 d coupled problems Fig. 5.1: Domains used in experiments with nonmatching discretization. (Left) Thespaces V h and Q H are inf-sup stable for (4.1) if h ≤ cH , c <
1. The condition issatisfied/violated in the top/bottom configurations. (Right) The 3 d -1 d experimentsuse two curves Γ. The mesh of Ω is obtained by first subdividing the domain intoodd number of cubes in each direction. Thus degrees of freedom of V h , Q H are notassociated with identical spatial points. Moreover h ≪ H is ensured in the refinement.Table 5.1: Iteration counts and error convergence for (4.1) and Ω a unit square andΓ a circle. The spaces V h and Q H are formed either by continuous linear Lagrangeelements or Q H uses discontinuous piecewise constant Lagrange elements. Note thatΓ is closed and thus Q H has the same dimension with either of the elements. Theinequality h ≤ cH , c < H (Ω) norm of the error u − u h . We note that the exactsolution is smooth. The error of the Lagrange multiplier measured in the s = − norm (2.8) (computed on Q H ) norm decays with order 1.5. dim V h dim Q H Q H continuous Q H discontinuous k u − u h k V k p − p h k Q k u − u h k V k p − p h k Q
22K 136 52 9.54E-02 5.28E-03 47 9.54E-02 3.68E-0387K 272 52 4.78E-02 1.71E-03 48 4.78E-02 1.15E-03348K 544 51 2.39E-02 5.77E-04 49 2.39E-02 4.18E-041.4M 1088 51 1.19E-02 1.87E-04 50 1.19E-02 1.49E-04
Due to the difficulties with the trace operator for Γ amanifold of codimension two, cf. §
2, the functional setting of (4.1) is not clear andtherefore corresponding discrete inf-sup conditions for the problem is not available.However, we shall assume that the inequality h ≤ cH , c <
1, which was cruacial forthe 2 d -1 d problems, plays a role also in the 3 d -1 d case and discretize the domainsaccordingly.The problem (4.1) is considered with two carefully constructed curves Γ, seeFigure 5.1, and Ω a unit cube discretized such that the inequality is ensured. As be-fore, the spaces Q H are constructed from continuous piecewise linear or discontinuouspiecewise constant Lagrange elements. We note that dim Q h ≪ dim V h . Further, the7MINRES iterations use the same initial and convergence conditions as in § A + M ) , H − . − ) is used as the preconditioner. In Table 5.2 we observethat the discretization and the preconditioner lead to bounded iteration counts. Wenote that if the discretization of Γ violates the inequality h < cH , the number ofiterations cannot be bounded anymore.Table 5.2: Iteration counts for (4.2) posed on Ω ⊂ R and the two curves picturedin Figure 5.1. For each domain, Q H from continuous linear (first column) or discon-tinuous constant (second column) Lagrange elements is considered. The domains arediscretized such that h ≤ cH , c <
1. In all the cases, the number of iterations isbounded. dim V h Square Spiraldim Q H Q H Q H Q H
6. Conclusions.
We have discussed preconditioning of a model multiphysicsproblem (1.1), where two elliptic subproblems were coupled by a trace constraint,bridging the dimensionality gap of size two. In order to facilitate the re-use of standardmultilelevel preconditioners for the 3 d domain we considered the trace as a mappingfrom H (Ω) to H s (Γ) for some s < − ∆) s . Using a simplerproblem (3.1) the spectral equivalence was investigated by a series of numerical exper-iments revealing for s ∈ ( − . , − .
1) existence of a mesh-dependent scale τ ( s, h ) suchthat τ ( s, h )( − ∆ h ) s is a robust preconditioner for the Schur complement. As the scaleis, in general, impractical to compute only the fractional Laplacian was further used inpreconditioning the coupled problem (1.1). Robustness of the proposed preconditionerwas demonstrated by numerical experiments with curves of different complexity andvarious shape-regular meshes using, at first, the assumption V h | Γ = Q h and finallywith spaces V h , Q H satisfying the compatibility condition h ≤ cH , c < d -1 d problems [37, 12]. Appendix A. Geometrical configurations and their discretization.
Nu-merical experiments with the Schur complement in § § = { ( t, , ); t ∈ [0 , } or Γ = { ( t, t, t ); t ∈ [0 , } . For eachcase the domains are discretized in three ways: ( uniform ) the meshes for Ω, Γ havethe same characteristic size, ( finer ) the mesh of Ω is finer at Γ than in the rest ofthe domain, ( coarser ) the mesh of Ω is coarser at Γ than in the rest of the domain.Parameters of the meshes for each refinement level are summarized in Table A.1. REFERENCES[1]
R. A. Adams and J. F. Fournier , Sobolev spaces , vol. 140, Academic press, 2003.[2]
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson,J. Ring, M. Rognes, and G. Wells , The FEniCS project version 1.5 , Archive of Numer-ical Software, 3 (2015). Preconditioning for 3 d -1 d coupled problems Table A.1: Sizes of FEM spaces and mesh parameters for different levels of refine-ments ( L ). The length of the largest cell in the mesh of Γ i is denoted by H . Forreadability the reported value is H × . Lengths of smallest/largest edges of cells ofthe mesh for Ω \ Γ i are respectively h min and h max . (Top) In uniform discretizationthe characteristic mesh size of Ω and Γ i triangulations are identical. (Middle) Finer discretization uses finer mesh near Γ i . (Bottom) In the coarser cases the mesh of Ωis coarser near the curve. L Γ = { ( t, , ); t ∈ [0 , } Γ = { ( t, t, t ); t ∈ [0 , } dim V h dim Q H h min
H h max H H dim V h dim Q H h min
H h max H H E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz,A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen , LAPACK Users’Guide , Society for Industrial and Applied Mathematics, Philadelphia, PA, third ed., 1999.[4]
M. Arioli, D. Kourounis, and D. Loghin , Discrete fractional Sobolev norms for domaindecomposition preconditioning , IMA Journal of Numerical Analysis, (2012), p. drr024.[5]
M. Arioli and D. Loghin , Discrete interpolation norms with applications , SIAM Journal onNumerical Analysis, 47 (2009), pp. 2924–2951.[6]
I. Babuˇska , The finite element method with Lagrangian multipliers , Numerische Mathematik,20 (1973), pp. 179–192.[7]
S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G.Knepley, L. C. McInnes, B. F. Smith, and H. Zhang , PETSc users manual , Tech.Report ANL-95/11 - Revision 3.4, Argonne National Laboratory, 2013.[8]
Michele Benzi, Gene H Golub, and J¨org Liesen , Numerical solution of saddle point prob-lems , Acta numerica, 14 (2005), pp. 1–137.[9]
Andrea Bonito and Joseph E. Pasciak , Numerical approximation of fractional powers ofelliptic operators , Math. Comp., 84 (2015), pp. 2083–2110.[10]
L. Cattaneo and P. Zunino , A computational model of drug delivery through microcirculationto compare different tumor treatments , International Journal for Numerical Methods inBiomedical Engineering, 30 (2014), pp. 1347–1371.[11] ,
Computational models for fluid exchange between microcirculation and tissue intersti-tium , Networks and Heterogeneous Media, 9 (2014), pp. 135–159.[12]
W. Dahmen and A. Kunoth , Appending boundary conditions by Lagrange multipliers: Anal-ysis of the lbb condition , Numerische Mathematik, 88 (2001), pp. 9–42.[13]
C. D’Angelo , Finite element approximation of elliptic problems with Dirac measure termsin weighted spaces: applications to one-and three-dimensional coupled problems , SIAMJournal on Numerical Analysis, 50 (2012), pp. 194–215.[14]
C. D’Angelo and A. Quarteroni , On the coupling of 1D and 3D diffusion-reaction equations:Application to tissue perfusion problems , Mathematical Models and Methods in AppliedSciences, 18 (2008), pp. 1481–1504.[15]
Zhonghai Ding , A proof of the trace theorem of Sobolev spaces on Lipschitz domains , Pro-ceedings of the American Mathematical Society, 124 (1996), pp. 591–600. [16] R. D. Falgout and U. Meier Yang , hypre: A library of high performance preconditioners ,in Computational Science ICCS 2002, P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan,and J. J. Dongarra, eds., vol. 2331 of Lecture Notes in Computer Science, Springer BerlinHeidelberg, 2002, pp. 632–641.[17] Q. Fang, S. Sakadˇzi´c, L. Ruvinskaya, A. Devor, A. M. Dale, and D. A. Boas , Oxygen ad-vection and diffusion in a three-dimensional vascular anatomical network , Optics express,16 (2008), pp. 17530–17541.[18]
L. Grinberg, E. Cheever, T. Anor, J. R. Madsen, and G. E. Karniadakis , Modeling bloodflow circulation in intracranial arterial networks: a comparative 3D/1D simulation study ,Annals of biomedical engineering, 39 (2011), pp. 297–309.[19]
N. Hale, N. J. Higham, and L. N. Trefethen , Computing a α , log A , and related matrixfunctions by contour integrals , SIAM Journal on Numerical Analysis, 46 (2008), pp. 2505–2523.[20] Stanislav Harizanov, Raytcho Lazarov, Pencho Marinov, Svetozar Margenov, andYavor Vutov , Optimal solvers for linear systems with fractional powers of sparse spdmatrices , arXiv preprint arXiv:1612.04846, (2016).[21]
Tobias K¨oppl, Ettore Vidotto, and Barbara Wohlmuth , A local error estimate for thePoisson equation with a line source term , in Numerical Mathematics and Advanced Ap-plications ENUMATH 2015, Springer, 2016, pp. 421–429.[22]
T. Koppl and B. Wohlmuth , Optimal a priori error estimates for an elliptic problem withDirac right-hand side , SIAM Journal on Numerical Analysis, 52 (2014), pp. 1753–1769.[23]
M. Kuchta, M. Nordaas, J. C. G. Verschaeve, M. Mortensen, and K.-A. Mardal , Pre-conditioners for saddle point systems with trace constraints coupling 2d and 1d domains ,SIAM Journal on Scientific Computing, 38 (2016), pp. B962–B987.[24]
Mateusz Kwa´snicki , Ten equivalent definitions of the fractional Laplace operator , FractionalCalculus and Applied Analysis, 20 (2017), pp. 7–51.[25]
A. A. Linninger, I. G. Gould, T. Marinnan, C.-Y. Hsu, M. Chojecki, and A. Alaraj , Cerebral microcirculation and oxygen tension in the human secondary cortex , Annals ofbiomedical engineering, 41 (2013), pp. 2264–2284.[26]
J. L. Lions and E. Magenes , Non-homogeneous boundary value problems and applications ,vol. 1, Springer Science & Business Media, 2012.[27]
A. Logg, K.-A. Mardal, and G. Wells , Automated solution of differential equations by thefinite element method: The FEniCS book , vol. 84, Springer Science & Business Media,2012.[28]
K.-A. Mardal and J. B. Haga , Block preconditioning of systems of PDEs , in AutomatedSolution of Differential Equations by the Finite Element Method, G. N. Wells et al. A. Logg,K.-A. Mardal, ed., Springer, 2012.[29]
K.-A. Mardal and R. Winther , Preconditioning discretizations of systems of partial differ-ential equations , Numerical Linear Algebra with Applications, 18 (2011), pp. 1–40.[30]
M. F. Murphy, G. H. Golub, and A. J. Wathen , A note on preconditioning for indefinitelinear systems , SIAM J. Sci. Comput., 21 (1999), pp. 1969–1972.[31]
M. Nabil and P. Zunino , A computational study of cancer hyperthermia based on vascularmagnetic nanoconstructs , Open Science, 3 (2016).[32]
P. Peisker , On the numerical solution of the first biharmonic equation , ESAIM: MathematicalModelling and Numerical Analysis - Modlisation Mathmatique et Analyse Numrique, 22(1988), pp. 655–676.[33]
J. Reichold, M. Stampanoni, A. L. Keller, A. Buck, P. Jenny, and B. Weber , Vasculargraph model to simulate the cerebral blood flow in realistic vascular networks , Journal ofCerebral Blood Flow & Metabolism, 29 (2009), pp. 1429–1443.[34]
M. Renardy and R.C. Rogers , An Introduction to Partial Differential Equations , Texts inApplied Mathematics, Springer New York, 2006.[35]
T. Rusten and R. Winther , A preconditioned iterative method for saddlepoint problems ,SIAM J. Matrix Anal. Appl., 13 (1992), pp. 887–904.[36]
D. Silvester and A. Wathen , Fast iterative solution of stabilised Stokes systems part ii: usinggeneral block preconditioners , SIAM Journal on Numerical Analysis, 31 (1994), pp. 1352–1367.[37]
O. Steinbach , Numerical Approximation Methods for Elliptic Boundary Value Problems: Fi-nite and Boundary Elements , Texts in applied mathematics, Springer New York, 2007.[38]