Robust preconditioning for coupled Stokes-Darcy problems with the Darcy problem in primal form

The coupled Darcy-Stokes problem is widely used for modeling fluid transport in physical systems consisting of a porous part and a free part. In this work we consider preconditioners for monolitic solution algorithms of the coupled Darcy-Stokes problem, where the Darcy problem is in primal form. We employ the operator preconditioning framework and utilize a fractional solver at the interface between the problems to obtain order optimal schemes that are robust with respect to the material parameters, i.e. the permeability, viscosity and Beavers-Joseph-Saffman parameter. Our approach is similar to our earlier work, but since the Darcy problem is in primal form, the mass conservation at the interface introduces some challenges. These challenges will be specifically addressed in this paper. Numerical experiments illustrating the performance are provided. The preconditioner is posed in non-standard Sobolev spaces which may be perceived as an obstacle for its use in applications. However, we detail the implementational aspects and show that the preconditioner is quite feasible to realize in practice.



There are no comments yet.


page 12


Robust monolithic solvers for the Stokes-Darcy problem with the Darcy equation in primal form

We construct mesh-independent and parameter-robust monolithic solvers fo...

Robust preconditioning of monolithically coupled multiphysics problems

In many applications, one wants to model physical systems consisting of ...

Analysis of the diffuse interface method for the Stokes-Darcy coupled problem

We consider the interaction between a free flowing fluid and a porous me...

Parameter-robust methods for the Biot-Stokes interfacial coupling without Lagrange multipliers

In this paper we advance the analysis of discretizations for a fluid-str...

Mass-Conserving Implicit-Explicit Methods for Coupled Compressible Navier-Stokes Equations

Earth system models are composed of coupled components that separately m...

Dual-Primal Isogeometric Tearing and Interconnecting methods for the Stokes problem

We are interested in a fast solver for linear systems obtained by discre...

The art of coarse Stokes: Richardson extrapolation improves the accuracy and efficiency of the method of regularized stokeslets

The method of regularised stokeslets is widely used in microscale biolog...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Let , where is the domain of the viscous flow, is the domain of the porous media and their common interface. Further let the domain boundaries be decomposed as and , where subscripts signify respectively that Dirichlet and Neumann boundary conditions are prescribed on the part of the boundary. The boundary of , i.e., the intersection of and is denoted by . An illustration is given in Figure 1.
The Stokes problem reads:

(1) (2)
while the Darcy problem in primal form reads:
Figure 1: Schematic domain of Darcy-Stokes problem. Dirichlet conditions shown in dashed line, and interface in red. Here, are the unknown velocity and pressure for the Stokes problem (1)-(2) in , is the unknown pressure of the Darcy problem (3) in . The material parameters are the fluid viscosity and the permeability . Here we shall consider the problem with the Dirichlet boundary conditions

and Neumann conditions

where , are the outer unit normals of the respective subdomains. In particular we assume that and for . Moreover, the coupled problem must be equipped with interface conditions expressing the continuity of stress as well as mass balance. We postpone their description until we describe the weak formulation of the problem.

The discretization of the coupled Darcy-Stokes problem with the Darcy problem in a mixed form is challenging since the Darcy and Stokes problems, respectively, call for different schemes. For example, typical finite element methods for the Darcy problem, like the Raviart-Thomas or Brezzi-Douglas-Marini elements, are not stable for Stokes problem as the discretization of the flux specifically targets the properties of rather than which is natural for Stokes discretizations. For this reason, a wide range of methods have been proposed over the last decade that address this particular challenge. For example, new elements robust for both the Darcy and Stokes problem have been proposed in Arbogast and Brunson (2007); Johnny and Michael (2012); Karper et al. (2009); Mardal et al. (2002); Zhang et al. (2009). Alternatively, stabilization or modifications of standard methods may be used as in Burman and Hansbo (2007); Feng et al. (2010); Xie et al. (2008). In this work we will consider the coupled problem with the Darcy equation in a primal form. Standard elements in both the Darcy and the Stokes domain will be used together with a Lagrange multiplier to couple the unknowns appropriately at the interface.

The well-posedness of the Darcy-Stokes problem coupled together through the use of a Lagrange multiplier is well-known when the Darcy problem is in mixed form Layton et al. (2002); Galvis and Sarkis (2007), where both the continuous setting and various discretizations were proposed. Other solution and discretization algorithms for the coupled problem are presented in e.g. Rivière and Yotov (2005); Gatica et al. (2008), see Discacciati and Quarteroni (2009); Layton et al. (2002) for an overview. For the mixed formulation we have, in our previous work Holter et al. , developed monolithic solvers that are robust with respect to all material parameters by utilizing fractional solvers on the interface. Here, we continue with the same type of approach, but address the difficulty of the Darcy problem in primal form. We remark that the problem to be studied further is symmetric and includes an explicit variable, the Lagrange multiplier, on . In this respect it differs from the more common primal formulation, which leads to a non-symmetric system to be solved for , and . Well-posedness of the latter problem was established in Discacciati and Quarteroni (2009) with efficient solvers proposed and analyzed e.g. in Discacciati and Quarteroni (2003); Cai et al. (2009); Discacciati et al. (2007).

An outline of the paper is as follows: Section 2 describes the notation, introduces the symmetric primal Darcy-Stokes problem and illustrates the difficulties in its preconditioning. The main challenge for the solver construction, i.e. the proper posing of the coupling operator, is addressed in Section 3. Parameter robust preconditioners are then established in Section 4.

2 Preliminaries

Let be a bounded Lipschitz domain in , =2 or 3, and denote its boundary by . We denote by the Lebesgue space of square integrable functions, with the norm , and by the Sobolev space of functions with first derivative in with norm . Note that the spaces are both Hilbert spaces, with the standard inner products. These spaces are defined in the same way when

is a vector field, in which case we will write

in boldface. We also define the subspace to be the completion in of , the space of smooth functions on whose restriction to is zero.

For a Lipschitz domain with , we can define a trace operator by for smooth . This can be extended to a bounded, surjective and right-invertible operator (cf. e.g.  Ding (1996)), where the space will be defined later. Given a subset of , we let , or for readability just , be the subspace of for which the restriction to is zero, where the restriction is defined in terms of the trace operator. Typically, will be the subset of on which Dirichlet conditions are prescribed. We also define the semi-norm on to be the norm of the tangential component of at . In 2D, this is just where is a tangent unit vector, while in 3D it is more conveniently written as .

For any inner product space , we let denote its inner product. When , we will omit the subscript if there is no cause for confusion. We write the space of continuous linear operators from to as , or just as if . For any two Sobolev spaces both contained in a common ambient space, we define the intersection and sum spaces and in terms of the norms

For any , we define the scaled space to be just as a set, but with the inner product . Its norm is trivially equivalent to , but because the equivalence constant depends on , the distinction between the two norms becomes important when we need to establish the independence of bounds with respect to problem parameters.

We define the fractional space following Kuchta et al. (2016). Let be the operator such that for all . We can then find a basis of

of orthonormal eigenfunctions


with eigenvalues

. Writing in this basis, we define the norm for any . Further, let the space be the completion of with respect to . We also define the space in the same manner, except that we then apply Dirichlet boundary conditions by choosing in . Furthermore, is the completion of rather than .

For the sake of completeness we review here the construction of a matrix realization of fractional operators given in Kuchta et al. (2016). To this end let , be a finite dimensional finite element subspace with basis functions , and , be the symmetric positive definite (stiffness and mass) matrices such that

In case and piecewise constant (P0) discretization is used we let

where is a set of all the facets of the finite element mesh. Further the (facet) average and jump operators are defined as , with and the two cells sharing facet . When is an exterior facet, we define , where is the unique cell with as facet.

It follows that the generalized eigenvalue problem

has only positive eigenvalues and a complete set of eigenvectors that form the basis of

so that the powers of are well defined. For we then set . Letting

be the vector of degrees of freedom of

, i.e. , we finally have

When is a vector function, we define the normal trace using the trace operator component-wise. As such is a continuous map . Moreover, we let be the tangential trace operator. We remark that in 2D and 3D the operator maps to scalar, respectively vector fields. The normal derivative, , is more challenging to define properly in this context. Let us therefore briefly sketch an approach, which at least in the authors’ opinion at first glance seems like a natural starting point. However, as we will show, the approach does not yield robust preconditioners in our context. First, notice that if we impose additional regularity on and require that then is well defined. In detail, let and be a (harmonic) extension operator. Then clearly lies in because

This extra regularity assumption is, however, hard to express in the operator preconditioning framework. In particular, to the author’s knowledge, there are no standard finite elements that would enable us to exploit the extra regularity. A possible approach could be NURBS Hughes et al. (2005) or discretizations developed for fourth order problems. However, the latter often show poor performance for second order problems Nilssen et al. (2001).

Alternatively, we may attempt to define as a composition of the first order derivative operator, , with the 1/2 order normal trace operator . The composition could then be expected to be a 3/2 operator . From an operator preconditioning point of view, this would be feasible to realize, as we will see below. However, as we will demonstrate, robustness will not be obtained if we realize as a 3/2 operator. In fact, robustness is only obtained if is a first order operator, . We remark here that while the operator in a continuous setting is , in the discrete setting we will include a scaling parameter, i.e. the mesh size, because we use the finite element method. To see that this is reasonable, notice that for finite elements, the mass matrix, as representation of the identity, is differently scaled in different dimensions. In Example 3.2 we detail the scaling in a simplified example.

In order to demonstrate why posing the operator properly is required, let us now formulate the coupled Darcy-Stokes problem, where the Darcy problem is in primal form. As a starting point, let the Lagrangian of the coupled problem be,

Note that the sign of has been changed from (1). Here, the Lagrange multiplier in is used to ensure mass conservation, while the extra term , where , corresponds to the Beavers-Joseph-Saffman condition Mikelic and Jäger (2000).

The corresponding weak formulation is obtained by the first order optimality conditions of the Lagrangian, that is; , , , and . A variational formulation hence reads: Find such that


where the bilinear forms , are defined as


We shall refer to (4) as the (primal) Darcy-Stokes problem. Note that the resulting formulation is symmetric.

While appropriate function spaces are readily available for and their corresponding test functions, it is less clear what the appropriate requirements are for and . This will be addressed below.

Example 2.1.

Preconditioner for coupled Darcy-Stokes problem assuming . Let us assume that is a 3/2 operator so that for . Next, observe that since then . Per assumption the coupling term is so that the dual variable . In turn, we consider the following weak formulation: Find such that


The coefficient matrix associated with (4) reads


Assuming that the proposed spaces indeed lead to well-posed operator , the operator preconditioning framework Mardal and Winther (2011) yields as a preconditioner the Riesz mapping


In order to test the preconditioner, we solve problem (6) on , where and and the Dirichlet boundary domains are and , cf. Figure 1. The mesh is a uniform triangular mesh, consisting of equally sized isosceles triangles. To discretize (4), we use lowest order (P2-P1) Taylor-Hood elements for the Stokes velocity and pressure, while piecewise quadratic elements (P2) were used for the Darcy pressure and piecewise constant elements (P0) for the Lagrange multiplier. Discretization is carried out in the FEniCS library Logg et al. (2012), with coupling maps between the interface and domains and the fractional Laplacians being implemented by the extension FEniCSii Kuchta (2019).

Approximation of the preconditioner (8) is constructed by using single sweep of -cycle of algebraic multigrid BoomerAMG from the Hypre library Falgout and Yang (2002) for all the blocks except for the interface block, which is inverted exactly. Starting from a random initial vector, we count the number of iterations required to solve the preconditioned linear system using the MINRES solver from the PETSc library Balay et al. (2019) with convergence criterion based on relative tolerance of and absolute tolerance of . Additionally, the condition numbers of are computed using an iterative solver from the SLEPc library Hernandez et al. (2005). In the condition number computations the operator is computed exactly, that is, all the blocks are inverted by LU. We remark that the solver setup should be used also in the subsequent examples.

The results of the experiment are plotted in Figure 2. By the failure of the iteration counts to stabilize, we see that using as multiplier space does not lead to a robust preconditioner over the whole parameter range. Note, however, that in the regime where is significantly smaller than (i.e. the lower left region of the plots in Figure 2), iteration counts and condition numbers appear to be stable as the mesh is refined. In this regime, the norm of the multiplier space is dominated by the part from , which is determined by posing of the trace operator. This suggests that the choice of , i.e. wrong posing of the operator, is responsible for the lack of boundedness.

Figure 2: Mesh refinement vs. iteration counts (left) and condition numbers (right) for Example 2.1. All subplots share - and -axes. For fixed , the x-axis range in the iterations subplot extends from to . In the condition number plots the range is from to . In all cases, .

3 Approximating the trace normal gradient operator

A crucial step in the analysis of the Darcy-Stokes problem will be the mapping properties of the operator . As a computationally practical choice of space for the Darcy pressure is , we immediately run into the problem discussed in the preliminaries because cannot be defined on all of . This necessitates either an assumption of extra regularity or an alternative approach.

Motivated by the observation in Holter et al. (2017), that in a discrete finite element setting the trace operator is stable as a map , we propose an alternative approach to construct the preconditioners. We start off by outlining the construction of an operator which will be an approximation to . Suppose is a sufficiently regular subset of , and that is of co-dimension 1 in . The -thick envelope is a higher-dimensional approximation of . For any ,


where is a test function in .

Note that although the integral over is not well-defined for a general , the integral over is. Provided is sufficiently regular and sufficiently small, we assume that there exists a vector field on which approximates the normal vector of at . Using , we further assume that we can define a bounded extension along for which for any . Provided and can be defined, then for any we can define by

for any , thus defining the required map approximating . We assume that the resulting operator is both surjective and bounded, with , and that has a bounded right inverse.

We emphasize that is just an analytical tool constructed for the analysis in the continuous setting and that is not related to the mesh size . In fact, we can choose far smaller than the mesh size and for any practical purposes in computations we assume that will be practically identical to . We summarize the assumption as follows:

Assumption 1.

Given a sufficiently regular , is a bounded surjection which approximates on the subspace of on which can be defined. Further, has a bounded right inverse.

Although characterizing the conditions under which creftypecap 1 holds is beyond the scope of this paper, we motivate the existence of the required constructions in a few simple examples below.

Example 3.1.

Let be the axis, and be the positive half-plane. The construction of is then given by and for we let . This continuously extends to all of . Clearly as for . Given any , define by . Then the map continuously extends to a right inverse of , as by linearity .

Next, suppose is the unit disk, and its boundary. By parametrizing with e.g. polar coordinates, this case can be effectively translated to the above. is now the unit radial vector , and for any , . Again, this definition of extends to all of . Because and as , we again have as for . Analogously to the previous case, a right inverse can be defined by sending any to .

Before considering the Darcy-Stokes problem, we justify Assumption 1. First we consider a simplified example in order to illustrate how the scaling of mass matrices in different dimensions affect preconditioners constructed via the application of trace operators. Then, in Example 3.3 we construct preconditioners for a Poisson problem with a -constraint which is to be enforced by a Lagrange multiplier, cf. the Babuška problem Babuška (1973) involving the trace operator.

Example 3.2.

Trace constrained projection. Let be a bounded domain with and . We then consider the problem


Letting denote the Lagrange multiplier associated with the boundary constraint, the extrema , of the Lagrangian of (10) satisfy the variational problem: Find and such that


The operator of the preconditioned continuous problem then reads


where is to be constructed such that the condition number is bounded in the discretization parameter . Here we shall consider three constructions. We remark that when using the finite element method, the identity or the mass matrix has eigenvalues such that both the smallest and the largest eigenvalues scale as on uniform mesh. First we consider , with eigenvalues . Then, following Holter et al. (2017), we let , i.e., a matrix with eigenvalues . Finally, the choice of is included to show that the relevant trace space in (12) is not (by viewing the trace as an order 1/2 operator) so that dual variable would reside in .

We remark that the first two operators are in practical computations assembled as weighted mass matrices where the weights for the respective operators are 1 and inverse cell volume. Recalling Preliminaries §2 the matrix representation of the fractional operator is .

To compare the three preconditioners, we let be a unit square, . Further, the domain shall be discretized uniformly into isosceles triangles with size , see Figure 3. Considering finite element discretization by P2-P1 elements Table 3.1 lists spectral condition numbers of (12). It can be seen that only the preconditioner leads to results independent of .

The growth of the condition number in Table 3.1 due to the preconditioner with power indeed confirms that is not appropriate in our setting. An attempt to establish the trace space could be based on viewing the trace as an 1/2 operator. Starting from a formal calculation then leads to the space and as the multipler space. While we do not include here the results for we remark that the condition number behaves practically identically to .

8.72 24.08 4.63 12.11 47.84 4.63 16.91 95.15 4.63 23.70 189.9 4.63 33.31 379.2 4.63 46.90 758.0 4.63 66.12 1515 4.63 Table 3.1: Condition numbers of (12) with different preconditioners and discretization by P2-P1 elements on (us) mesh from Figure 3. Boundedness is obtained with the Schur complement preconditioner . P2-P1 P2-P0 (us) (uu) (nu) (us) (uu) (nu) 1 4.63 4.63 4.10 4.63 4.63 3.98 2 4.63 4.06 4.32 4.63 4.07 4.33 3 4.63 4.20 4.28 4.63 4.20 4.31 4 4.63 4.29 4.31 4.63 4.32 4.34 5 4.63 4.45 4.50 4.63 4.43 4.45 6 4.63 4.25 4.28 4.63 4.32 4.37 7 4.63 4.25 4.36 4.63 4.28 4.39 Table 3.2: Condition numbers of (12) with preconditioner using . Boundedness with different types of triangulations, cf. Figure 3, and discretizations can be observed.

In order to verify that the properties of preconditioner are not due to the highly structured mesh, we consider two additional discretizations of shown in Figure 3. In particular, the triangulations are obtained as refinements of the unstructured meshes where in one case the mesh size is uniform while in the other one the mesh is finer close to the multiplier domain . Moreover, using these triangulations, problem (12) shall be discretizated by P2-P1 elements as well P2-P0 elements to provide more evidence for the preconditioner construction. Indeed, Table 3.2 shows that the condition numbers of (12) are bounded irrespective of the underlying mesh and the finite element discretization considered.

Example 3.3.

Babuška problem with Neumann boundary conditions. Let be a bounded domain with the boundary partitioned into non-overlapping subdomains such that and . We will consider both the case that and later the case that . Let and consider the problem


With the Lagrange multiplier associated with -constraint (13) leads to a variational problem: Find and such that


The preconditioned continuous problem then reads


Following the preliminaries where was regarded as a 3/2 operator we let . Alternatively, is set following the Assumption 1. Finally is considered. Matrix realization of the operators shall be identical to Example 3.2. We shall also use the tessellations described in Example 3.2 as well as identical eigenvalue solvers.

To compare the three preconditioners we let be a unit square and and we consider first the (Neumann) case where , i.e. where the multiplier domain intersects the part of boundary with Neumann boundary conditions. Using the uniform meshes (marked as (us) Figure in 3) and P2-P1 elements, Table 3.3 shows the spectral condition numbers of (15). As in Example 3.2 only preconditioner (based on creftypecap 1) leads to results independent of .

11.99 6.70 4.88
14.55 9.27 4.88
18.47 12.89 4.88
24.44 18.01 4.88
33.25 25.26 4.88
45.96 35.52 4.88
64.10 50.02 4.88
Table 3.3: Condition numbers of (15) discretized by P2-P1 elements on uniform refinements of (us) mesh in Figure 3. Boundednes in discretization is obtained only with .
Figure 3: Parent meshes for uniform refinement. From left to right: uniform structured(us), uniform unstructured(uu), non-uniform unstructured(nu). Non-uniform mesh has finer (by factor 3) mesh size close to .

Table 3.4 shows that the performance of in (15) remains robust if different tessellations and finite element discretizations are used.

P2-P1 P2-P0
(us) (uu) (nu) (us) (uu) (nu)
1 4.88 4.77 6.64 3.49 3.49 3.06
2 4.88 5.98 6.56 3.49 3.04 3.37
3 4.88 5.78 5.67 3.49 3.24 3.36
4 4.88 6.31 6.67 3.49 3.40 3.40
5 4.88 5.25 5.68 3.49 3.44 3.48
6 4.88 5.71 5.89 3.49 3.41 3.44
7 4.88 6.14 6.61 3.49 3.35 3.47
P2-P1 P2-P0 (us) (uu) (nu) (us) (uu) (nu) 1 5.34 5.25 6.67 3.48 3.45 3.04 2 5.34 6.25 6.67 3.49 2.99 3.37 3 5.34 5.94 5.84 3.49 3.24 3.36 4 5.34 6.52 6.93 3.49 3.40 3.40 5 5.34 5.56 6.07 3.49 3.44 3.48 6 5.34 5.91 6.17 3.49 3.41 3.44 7 5.34 6.40 6.85 3.49 3.35 3.47
Table 3.4: Condition numbers of (15) using preconditioner discretized on uniform refinements of parent meshes in Figure 3 using two element types. Refinement level is indicated by . (Left) intersects . (Right) intersects .

In the context of multiscale problems, compatibility of boundary conditions of the multiplier space and the boundary conditions prescribed on the domain intesecting is known to present an issue, cf. e.g. Galvis and Sarkis (2007). Here, we address this problem by considering (15) with , i.e. we let intersect only the Dirichlet boundary. We remark that until this point only intersection with Neumann boundary was considered.

In Table 3.4 the Dirichlet problem is considered with an unmodified preconditioner. In particular, with P2-P1 discretization we impose no boundary conditions on the multiplier space. Using this construction the condition numbers can be seen to remain bounded on all the meshes and with both finite element discretizations.

We remark that the preconditioner is equally unaffected by the Dirichlet boundary conditions on in the trace-constrained projection problem (10) with , cf. Example 3.2, in contrast to the problems considered in Holter et al. , where the appropriate preconditioner was or depending on whether the interface intersected the Dirichlet boundary or not. We remark that in the continuous setting boundary values have measure zero and this may then be perceived as the space being the correct one in our discrete setting. Of course, the counterargument in the continuous setting is that then the trace cannot be defined. However, in the discrete setting, this can be done.

Without including the simulation results we comment here that the condition numbers of the Dirichlet problem are practically identical to those presented in Tables 3.1 and 3.2. In addition, with the two preconditioners and on the unstructured meshes a growth of condition numbers with is observed similar to Table 3.3.

We remark that the stability of the preconditioner in Example 3.3 provides numerical evidence for well-posedness of (14), i.e. the Darcy subproblem in the coupled Darcy-Stokes system (4).

4 Robust Preconditioners for the Darcy–Stokes system

In Example 2.1, we showed that the efficiency of the preconditioner (8) for the primal Darcy–Stokes problem (7) varied substantially with the material parameters even though the Stokes block and the Darcy block were preconditioned with appropriate preconditioners, and argued that the reason was a poor preconditioner at the interface.

In this section we demonstrate that robustness with respect to mesh resolution and variations in material parameters can be obtained by posing the Lagrange multiplier in properly weighted fractional spaces, namely the intersection space . No modifications of the velocity or pressure space norms will be required. Our analysis is closely related to Holter et al. , and based on creftypecap 1 along with an assumption of stability for the Stokes problem. We remark that although creftypecap 1 is motivated by the discrete problem, our analysis is carried out in a continuous setting.

Let for such that . We shall prove well-posedness of the coupled Darcy-Stokes problem (4) with spaces


We remark that in case intersects only the Dirichlet boundary the space needs to be modified to reflect as the appropriate trace space of . We refer to Holter et al. for a thorough discussion of the subject.

As a prerequisite for the coupled problem to be well-posed, we require that each subproblem is well-posed. For the Stokes subproblem the property has been demonstrated by numerical experiments in Holter et al. . Here we state the result without proof.

Assumption 2.

Let be such that , and . We define and the forms