Numerical Modeling of the Leak through Semipermeable Walls for 2D/3D Stokes Flow: Experimental Scalability of Dual Algorithms

: The paper deals with the Stokes ﬂow subject to the threshold leak boundary conditions in two and three space dimensions. The velocity–pressure formulation leads to the inequality type problem that is approximated by the P1-bubble/P1 mixed ﬁnite elements. The resulting algebraic system is nonsmooth. It is solved by the path-following variant of the interior point method, and by the active-set implementation of the semi-smooth Newton method. Inner linear systems are solved by the preconditioned conjugate gradient method. Numerical experiments illustrate scalability of the algorithms. The novelty of this work consists in applying dual strategies for solving the problem.


Introduction
The no-slip condition is the standard boundary condition in mathematical fluid flow models. It expresses the fact that a fluid adheres to the surface S, i.e., the velocity u = 0 on S. Nevertheless, in many situations a motion of the fluid on S such as a slip along impermeable S or a leak through semipermeable S can be observed. The Navier slip condition is the simplest one [1]: σ t = −ku t , where k > 0 is the adhesive coefficient, σ t , u t denotes the shear stress, and the tangential component of u, respectively. From its form we see that the fluid starts to slip whenever σ t = 0. However, this is not the case, in general. For example, a water drop on an inclined plane which is coated by a non-wetting (hydrophobic) material (teflon, e.g.,) slips only if the angle inclination reaches a specific value. Similarly, a leak of a fluid through semipermeable S may occur only if the magnitude of the stress vector σ (or some of its components) on S attains a certain critical value. The true slip or leak boundary conditions should reflect this threshold behavior. To this end one can use tools of convex analysis and write them in a compact form as inclusions for set-valued mappings, defined by the subgradient of appropriate convex functionals [2]. The overview of different threshold slip laws is presented and discussed in [3]. Since "slip" in fluid mechanics is a synonym for "friction" in solid mechanics (in particular in contact mechanics), it is not surprising that friction and slip models use the same terminology. The weak form of mathematical models under threshold slip or leak conditions leads to a variational inequality type problem, intricacy of which depends on the used law. The simplest one is the threshold slip and leak condition of the Tresca type defined by: (slip) |σ t | ≤ g, u t = 0 then σ t = −g u t |u t | on S, and (leak) |σ n | ≤ g, u n = 0 then σ n = −g u n |u n | on S, (2) respectively, where g is an a-priori given threshold bound, | · | stands for the absolute value of a real number or the Euclidean norm in (1) for 3D problems. The Stokes system under (1) and (2) has been studied by Fujita [4,5] and by Roux [6] who considered more general slip models, in which the bound g may depend on |u t |. Despite the fact that the mathematical analysis of this class of problems began relatively late (compared with friction problems of contact mechanics), a lot of theoretical papers on this topic has appeared during past twenty years. Results have been extended to other fluid models and slip conditions such as the steady and unsteady Stokes and Navier-Stokes equations, Coulomb's type friction law [7][8][9][10][11] and the references therein. The Signorini unilateral boundary conditions [12]: u n + g ≥ 0, σ n + σ * ≥ 0, (u n + g)(σ n + σ * ) = 0 on S, with given g, σ * is another type of conditions which leads to an inequality formulation. The authors used (3) in [13,14] as the artificial outflow conditions prescribed on a part of the boundary to model the blood flow in large arteries.
One of difficulties we face in numerical solution of variational inequalities is the fact that problems are generally nonsmooth and consequently, appropriate computational methods have to be used [15,16]. For the Stokes or Navier-Stokes system with slip or leak conditions nonsmoothness is caused by the presence of the nondifferentiable functional j whose subgradient defines (1) or (2). In the case of the Signorini conditions nonsmoothness is due to the kinematical constraint u n + g ≥ 0 on S which imposes additional restrictions on the set of admissible velocity fields. The velocity-pressure formulation of such problems is discretized using a pair of finite element spaces satisfying the Babuska-Brezzi condition [17,18]. One way to handle the nondifferentiable term j, the kinematical constraint in (3), is to use a regularization of j and a smooth penalization of the constraint, respectively in order to transform originally the nonsmooth problem into a sequence of smooth ones, which can be solved by standard methods. The regularization of j has been used in [19], and the penalty method in [13]. Another and the most frequent way of releasing constraints and the nonsmooth character of mathematical models is based on a dualization approach. In the case of the Stokes and Navier-Stokes system, this leads to the weak velocity-pressuretype formulation which is enhanced by the additional Lagrange multipliers. The resulting weak formulation is formally the same as the Karush-Kuhn-Tucker conditions in saddle point problems. The majority of papers devoted to the convergence/error analysis of this class of problems uses just the dualization approach [20][21][22][23][24][25][26][27]. Moreover, there exists a wide range of numerical methods for efficient computational realization of the resulting algebraic models [16].
The present paper is devoted to computational aspects of the Stokes problem with threshold leak conditions. We use the dual strategy, i.e., the velocity component is eliminated from an appropriate mixed finite element discretization of 2D and 3D problems. Our aim is to develop efficient algorithms for solving resulting algebraic systems which are satisfied by the dual variables, namely by the discrete pressure, shear and the normal stress on S, in our case. Let us note that this approach is frequently used in contact problems of solid mechanics. However, its simple transfer to problems of fluid mechanics is not possible due to the presence of the incompressibility condition prescribed in the whole computational domain. The modified path-following (PF) variant of the interior point method has been used in [27] for solving 2D Stokes system with threshold slip conditions, and in [28] for 3D problems. In addition, the latter paper uses also the semi-smooth Newton (SSN) and compares it with the PF method. Although numerical experiments demonstrated the scalable behavior of both approaches, the SSN method turned out to be more efficient because of its simpler implementation. Similar numerical tests and comparisons of SSN and PF will be done in the present paper. It is hard to predict the computational efficiency in advance, since the discrete threshold leak conditions in 2D and 3D problems always lead to simple box constraints unlike the separable spherical constraints in the case of the threshold slip in 3D. Another important feature is the fact that the dual algebraic formulations involve relatively a small number of constrained unknowns. This excludes the efficient use of some types of algorithms such as strictly feasible minimization methods based on an active-set strategy [29]. It has been found in [27] for 2D problems that they are less efficient than the PF method. The inefficiency is more significant for large-scale 3D problems. The main benefit of this work is an application of dual strategies for solving the problem. The methods presented in this paper may be easily modified for parallelization based on the FETI domain decomposition technique [30].
The paper is organized as follows. In Section 2 we present the classical and several weak formulations of the problem. The mixed finite element approximation based on the P1bubble/P1 finite element pair is introduced in Section 3. The resulting algebraic problems are starting point for constructing of algorithms in Section 4. Finally, Section 5 contains results of numerical experiments, including the problems with non-unique solutions. Concluding remarks are summarized in Section 6.
The following notation is used in the paper: H k (Ω), k ≥ 0 integer, stands for the Sobolev space of functions defined in Ω which are together with their generalized derivatives up to order k square integrable in a domain Ω ⊂ R d (H 0 (Ω) = L 2 (Ω)). The scalar product in L 2 (Ω) will be denoted by ( , ). L 2 + (Ω) denotes the set of non-negative functions from L 2 (Ω). If γ is a part of the boundary of Ω then H 1 2 (γ) is the trace space on γ of functions belonging to H 1 (Ω). The space of p × q matrices is denoted by R p×q , R p := R p×1 is the space of p-dimensional vectors, and R p + denotes the non-negative orthant of R p . Bold characters are used for vectors and matrices. By 0 we denote the zero matrix, or vector. Further I ∈ R p×p stands for the identity matrix. The scalar product of two vectors a, b ∈ R p , is denoted by a · b. If A = (a ij ) ∈ R p×q , B = (b ij ) ∈ R p×q are two matrices then A : B = a ij b ij (summation convention is used). Caligraphic symbols will be used for index sets, for instance: N = {1, 2, . . . , n} and A, I ⊆ N . If N ∈ R p×q and A = ∅, then N A ∈ R |A|×q stands for the matrix given by the rows of N whose indices belong to A. If A ∈ R p×p is symmetric, positive definite with the smallest, and largest eigenvalues 0 < α min ≤ α max , respectively, then the spectral condition number of A is cond(A) = α max /α min . The indicator matrix of a subset S ⊆ N is the diagonal matrix D(S ) = diag(s 1 , s 2 , . . . , s p ) ∈ R p×p , where s i = 1 for i ∈ S and s i = 0 if i ∈ S. Finally, the symbol · stands for the Euclidean norm in R p , p ≥ 2.

Setting of the Problem
The aim of this section is to recall the classical formulation of the Stokes system with leak boundary conditions and to present several weak formulations of this problem together with main existence/uniqueness results.
Let Ω ⊂ R d , d ∈ {2, 3} be a bounded domain with the Lipschitz boundary ∂Ω which is decomposed into two nonempty, nonoverlapping parts Γ and S both open in ∂Ω. The classical formulation of the problem reads as follows: find the velocity vector u : Ω → R d and the pressure p : Ω → R satisfying the following system of differential equations and boundary conditions: in Ω, u = 0 on Γ, u t = 0 on S, |σ n | ≤ g on S, u n = 0 ⇒ |σ n | = g & u n σ n ≤ 0 on S.
Here µ > 0 is the viscosity of the fluid, f : Ω → R d , g : S → R + denote an external force, and a nonnegative leak threshold, respectively. Further Du = 1 2 (∇u + (∇u) ) is the symmetric part of the gradient of u, n stands for the unit, outward normal vector to ∂Ω, u n = u · n, u t = u − u n n is the normal, and the tangential component of u on ∂Ω, respectively. Similarly, σ n = 2µDun · n − p is the normal component of the stress vector σ = 2µDun − pn on ∂Ω. It is worth noticing that unlike the shear stress σ t = σ − σ n n, the normal component σ n depends explicitly on the pressure p. The boundary conditions on S express the fact that there is no slip along S and the fluid can escape from Ω only if the threshold g is attained, i.e., |σ n | = g. In addition, the sign of σ n is opposite to the one of u n . Remark 1. The last two conditions in (4) can be written in the following equivalent way: where the symbol ∂ stands for the subgradient of convex functions.
To present variational formulations of (4) we shall need the following function spaces and forms: By the weak velocity-pressure formulation of (4) we call a problem of finding a pair (u, p) ∈ V(Ω) × L 2 (Ω) such that From Green's formula with an appropriate choice of test functions v one can verify that the classical and weak velocity-pressure formulations are equivalent for sufficiently smooth solutions (u, p) ( [5]).
Restricting to test functions v ∈ V div (Ω) we arrive at the weak velocity formulation of (4): which is equivalent to the following minimization problem: where Problem (8) is the variational inequality of the second kind and it has a unique solution u [15].
The following existence/uniqueness result has been established in [4,5]. Theorem 1. Problem (6) has a solution (u, p) for any f ∈ L 2 (Ω) d and g ∈ L 2 + (S). The velocity vector u is unique. If there is a set ω ⊆ S, meas d−1 ω > 0 such that |σ n | = g on ω then the pressure p is unique, as well, otherwise p is determined up to constants c satisfying |σ n + c| ≤ g a.e. on S.

Remark 3.
It has been shown in [5] that, if p is not unique, then the range of the admissible constants c satisfying |σ n + c| ≤ g a.e. on S is the interval k 1 , k 2 , where k 1 = sup S (σ n − g), k 2 = inf S (σ n + g), provided that the solution (u, p) is smooth enough.
The last formulation presented in this section is the variant of the velocity-pressure formulation (6) in which the Lagrange multipliers are used to release the no-slip condition u t = 0 on S and to regularize the nonsmooth term j. To this end we shall need the following trace spaces: and their duals X t , X n . If Ω is the domain with a smooth boundary (C 1,1 , e.g.,), then X t = (H 1 2 (S)) d−1 and X n = H 1 2 (S).

Remark 4.
If the solution (u, p) is sufficiently smooth, then the normal and shear components of σ on S have been defined pointwisely at the beginning of this section. If no regularity assumptions are imposed on (u, p) then σ n and σ t have to be understood as elements of X n and X t , respectively, defined by where the , stands for the duality pairings. It is easy to see that the right hand sides of (9) and (10) depend only on v t , and v n on S, respectively, since (u, p) ∈ V div (Ω) × L 2 (Ω) solves the Stokes system in a weak sense. Moreover, σ n belongs to L 2 (S) and it satisfies the last two conditions on S in (4). Indeed, The duality pairing σ n , v n is defined by the L 2 (S)-scalar product in this case.
be a closed, convex subset of L 2 (S). To derive the new formulation we use the saddle-point approach. It holds: It is known (see [2]) that the necessary and sufficient condition for (u * , p * , τ * , ν * ) to be a saddle-point of L on Q is the satisfaction of the Karush-Kuhn-Tucker conditions: It is easy to see that (u * , p * ) ∈ V div (Ω) × L 2 (Ω) satisfies the Stokes system and the boundary conditions on ∂Ω, i.e., (u * , p * ) solves (6). Moreover, τ * = σ * t , ν * = σ * n is the corresponding shear, and normal stress on S, respectively. The opposite assertion is also (11), where σ t , σ n is the corresponding shear, and normal stress on S, respectively.
Proof is straightforward. For more details we refer to [31].

Remark 5.
In the discrete counterpart of the problem we shall use the following more general form j κ instead of j: Then the weak formulation (6) is modified as follows: find a pair (u, The function κ : S → R + characterizes pore opening. The respective leak condition reads as Therefore the substitution σ n = σ n + κu n transforms the leak condition (13) formally into (5) so that all results of this section remain valid also for problem (12).

Mixed Finite Element Method and Algebraic Formulations
In this section we present the algebraic form of the problem, which will be used in the subsequent parts of the paper for numerical solution. It is based on the discretization of the velocity-pressure formulation (12) by a mixed finite element method. Next we shall suppose that Ω is a polygonal (d = 2), or polyhedral (d = 3) domain. If not, then Ω is approximated by them. Let V h ⊂ W h , Q h be finite element approximations of V(Ω), W(Ω), and L 2 (Ω), respectively, such that the bilinear form b satisfies the Babuska-Brezzi condition: where β > 0 does not depend on the discretization parameter h. The mixed finite element approximation of (12) reads as follows: In computations we use the P1-bubble/P1 finite element pair [32] on a regular partition T h of Ω into triangles, tetrahedras for d = 2, and 3, respectively. On T h we define the finite element spaces: where P 1 (T), B(T) are the spaces of polynomials of degree one and of bubble functions of degree three (d = 2) and four (d = 3) on T ∈ T h , respectively.
The algebraic counterpart of (14) reads as follows: where the absolute value is understood componentwisely, D κ = diag(κ), g, κ ∈ R n s + with n s being the number of nodes belonging to S \ Γ. Here, n p , n u stand for the number of the finite element nodes of T h in Ω, and Ω \ Γ, respectively, while n q denotes the number of the triangles/tetrahedras in T h . Further, A ∈ R d(n u +n q )×d(n u +n q ) is a symmetric, positive definite diffusion matrix, B ∈ R n p ×d(n u +n q ) is a full row-rank divergence matrix, and b ∈ R d(n u +n q ) is a vector of nodal forces. The unit outward, normal vector at the i-th node of S \ Γ defines nonzero entries in the i-th row of the matrix N ∈ R n s ×d(n u +n q ) . For d = 2, the tangential vector defines nonzero entries in the i-th row of the matrix T ∈ R n s ×2(n u +n q ) . For d = 3, there are two mutually orthonormal tangential vectors that define nonzero entries in the i-th rows of the matrices T 1 , T 2 ∈ R n s ×3(n u +n q ) . In this case, that is the local index set of the nodes belonging to S \ Γ. The algebraic version j h κ of j introduced above is the result of numerical integration, assuming that g, κ ∈ C(S). For d = 3 we use the following integration formula on the triangular element τ ∈ T h|S with the vertices where |τ| = meas τ, x τ is the centre of gravity of τ, and (Nv) i is the i-th component of Nv.
Summing up (17) over all τ ∈ T h|S we arrive at: where g = (g 1 , . . . , g n s ) T and n i is the number of the triangles τ i j ∈ T h|S sharing x i ∈ S \ Γ as the common vertex. In the same way we get: where κ = (κ 1 , . . . , κ n s ) T . We proceed analogously for d = 2 when τ ∈ T h|S is a line segment, n i ≤ 2, and From this and (16) Thus problem (15) can be written in the following equivalent form: where and It is easy to show that the first component u in (20) solves the discrete velocity formulation being the algebraic version of (8): where Formulation (22) is not suited for direct computations, as the constraints in V B can be hardly handled for large-scale problems. Moreover, the function J is nondifferentiable caused by the term j h . To overcome these difficulties, we will use the dual formulation of (22) and derive the discrete counterpart of (11).
We introduce three Lagrange multipliers: µ n ∈ R n s to regularize j h , µ t ∈ R (d−1)n s , q ∈ R n p to release the discrete stick, and incompressibility conditions, respectively and denote µ = (µ T n , µ T t , q T ) T . Let where µ n,i is the i-th component of µ n ∈ R n s . The term j h can be written as follows: Thus: min Until now, the velocity vector v ∈ R d(n u +n q ) incorporates dn q bubble components. These components are usually eliminated before the computational process [33]. In our case, we perform this elimination in a saddle-point formulation for L # , which leads to the reduced Lagrangian L : To simplify notation here and in what follows, we use the same symbols for the corresponding matrices and vectors before and after the elimination of the bubble components. The dimensions of the reduced matrices are: A κ ∈ R dn u ×dn u , N ∈ R n s ×dn u , T ∈ R (d−1)n s ×dn u , B ∈ R n p ×dn u , and b ∈ R dn u . Note that these matrices preserve the same properties as before the elimination, especially, the expression (21) remains valid with A ∈ R dn u ×dn u being symmetric, positive definite. The presence of the symmetric, positive semidefinite matrix E ∈ R n p ×n p and of the vector c ∈ R n p is due to this elimination. Note that E has the defect one and the eigenvector whose all components are equal to 1 corresponds to zero eigenvalue.
The new saddle-point formulation of (22) reads as follows: or, equivalently, where λ = (λ T n , λ T t , p T ) T .

Remark 6.
Here, λ n , λ t , p are the discrete counterparts of − σ n (see Remark 5), −σ t , and p, respectively. The change of the sign at σ n and σ t is for convenience of readers. One can see at a glance that the dual Hessian will be positive definite, as C introduced bellow is assembled from matrices without any negative signs.
Now we eliminate the unknown u. Denoting C = (N T , T T , B T ) T , the first equation in (24) Inserting u into the first inequality in (23), we arrive at the dual formulation of (22): where S : R dn s +n p → R, In what it follows we will assume that F κ is non-singular. The case Γ N = ∅ will be discussed in Example 3.
It should be noted that (25) is more convenient for numerical solution unlike (22) as the function S is differentiable (quadratic) and the feasible set X is defined by the box constraints. The solution to (25) may be computed by an appropriate algorithm of the constrained minimization. Another way how to solve our problem originates from (24) but with the algebraic leak boundary conditions expressed by projections on appropriate convex sets. The resulting system of nonsmooth algebraic equations can be solved by the nonsmooth Newton type method.

Algorithms
The aim of this section is to present main ideas of the path-following and the semismooth Newton methods that turned out to be highly efficient for solving frictional contact problems of solid mechanics [34][35][36]. These algorithms have been adapted for the Stokes system with the threshold slip boundary conditions in 2D [27,37] and 3D [28]. The situation with the threshold leak conditions is different since only box constraints are imposed on the individual components of λ n in 2D, as well as in 3D problems. In view of this fact the same implementation of the path-following method can be used in both cases. Also the semismooth Newton method works only with projections on compact intervals unlike the threshold slip boundary conditions in 3D in which case the components of λ t are pairwise subject to spherical constraints and projections on circles are needed. Consequently, the finite termination property for the semismooth Newton method can be expected also for 3D problems.

Path-Following Method
Let L : R dn s +n p × R 2n s + → R be the Lagrangian to (25): where ν := (ν T 1 , ν T 2 ) T ∈ R 2n s , ν ≥ 0, is the Lagrange multiplier vector releasing the box constraints appearing in X. Let z k := −∇ ν k L(µ, ν), k = 1, 2, z := (z T 1 , z T 2 ) T ∈ R 2n s be the new variables and define the function H : R (d+4)n s +n p → R (d+4)n s +n p , where ω = (µ T , ν T , z T ) T ∈ R (d+4)n s +n p , M = diag(ν), Z = diag(z), and e ∈ R n s is the vector whose all components are equal to 1. The solution λ to (25) is the first component of the vectorω = (λ T ,ν T ,z T ) T , which satisfies since (27) is equivalent to the Karush-Khun-Tucker conditions for (λ,ν) to be a saddle-point of (26) on R dn s +n p × R 2n s + . To derive the path-following algorithm, we replace (27) by the following perturbed problem: Find where τ ∈ R + . Solutions ω τ to (28) define a curve C(τ) in R (d+4)n s +n p called the central path. This curve approachesω when τ tends to zero. We combine the damped Newton method used for solving the equation in (28) with an appropriate change of τ which guarantees that the iterations belong to a neighbourhood N (c 1 , c 2 ) of C(τ) defined by: where c 1 ∈ (0, 1], c 2 ≥ 1, and ϑ := ϑ(ω) = ν T z/(2n s ). In the k-th iteration, we replace τ := τ (k) by the product of ϑ (k) := ϑ(ω (k) ) with the centering parameter c (k) chosen as in [38]. The algorithm uses also the Armijo-type condition (30) ensuring that the sequence {ϑ (k) } is monotonically decreasing. By J(ω) in (29), we denote the Jacobian matrix of H at ω.
The bounds imposed on the parameters mentioned in the initialization section follow from the convergence analysis presented in [34]. The computational efficiency depends on the way how the inner linear systems (29) are solved. The Jacobian matrix is non-symmetric and indefinite with the following block structure: Eliminating the 2nd and 3rd component of ∆ω (k+1) , we get the reduced linear system for ∆µ (k+1) with the Schur complement As µ (k) > 0, z (k) > 0, the matrix S (k) κ is symmetric, positive definite and the reduced linear system can be solved by the conjugate gradient method. In order to guarantee its convergence, we use the preconditioner: κ depend on the k iteration, the eigenvalue bounds of the preconditioned matrix (P (k) κ are positive and do not depend on k. Consequently, the spectral condition number cond is bounded by (see [34]): In computations, we approximate D F κ replacing A −1 κ in F κ by diag (A κ ) −1 , i.e., The inequality (32) remains valid also for this choice od D F κ . The conjugate gradient method (CGM) used in the k-th step of ALGORITHM PF is initialized and terminated adaptively. The initial CGM iteration is taken as the computed result in the previous iteration. The CGM iterations are terminated if the relative residuum is less than the bound, depending on the precision achieved in the outer loop err (k−1) and determined by: where 0 < r tol < 1, 0 < c fact < 1, err (−1) = 1, and tol (−1) = r tol /c fact .

Semismooth Newton Method
Let P [a,b] : R → [a, b] be the projection of R on the interval [a, b], a ≤ b, defined by the max-function: The last two lines in (24) representing the leak condition can be equivalently written as: where ρ i > 0, i ∈ N , are arbitrary but fixed parameters. We introduce the new variable s n ∈ R n s whose components s n,i = κ i (Nu) i + λ n,i approximate −σ n (x i ) at the nodes x i ∈ γ S \ γ D , i ∈ N . Then (34) takes the form: We distinguish two cases: The equations in (35) can be reformulated using (33). As we prefer the vector notation, we split N into two subsets: and introduce the sub-vectors s n, N 0 , s n, N + , g N 0 , g N + and the sub-matrices N N 0 , N N + , D κ + = diag(κ N + ). Then (35) is equivalent to: where the max-function is understood component-wisely. The max-function can be efficiently replaced by indicator matrices. First we define inactive/active sets for the second arguments of the max-functions in (36): where G : R d(n u +n s )+n p → R d(n u +n s )+n p is defined at using also (21). The standard differentiation rules lead to the following (generalized) Jacobian matrix of G at y: The Equation (37) will be solved by the Newton-type iterations: This iterative process generates a sequence {y (k) } starting from y (0) as an initial approximation. The right hand-side in (38) reads as follows: From this it is easy to show that s (k+1) n, and that the remaining components of the new iteration y (k+1) in (38) solve the following linear system: To solve the linear systems (39), we use the Schur complement S (k) to A-block defined by: where The right hand-sides of the Schur complement linear systems are given by Note that the dependence on the k iteration is through the active set A 0 . We arrive at the implementation of (38), in which the iterations are performed only with the last four components of y assembled in the vectorλ = (s T n, N 0 , s T n, N + , λ T t , p T ) T . ALGORITHM SSN: Letλ (0) ∈ R dn s +n p , ε ≥ 0, ρ > 0, and set k := 0.
and solve the Schur complement linear system for the remaining components ofλ (k+1) assembled in the vector inλ (k+1) .
The algorithm is the dual variant of the semismooth Newton method with the superlinear convergence rate established in [39] provided that the initial iteration is sufficiently close to the solution. The overall computational efficiency depends on the way how the (inner) linear systems are solved. Although the size of S (k) depends on the cardinality of A 0 (i.e., on the k iteration), it is easily seen that these matrices are always positive definite. Hence, the linear systems (41) can be solved by the conjugate gradient method. This method in the k-th step of ALGORITHM SSN is initialized and terminated adaptively using the same ideas as in ALGORITHM PF. Note that the preconditioning of the conjugate gradient iterations is needed. This follows from the presence of D −1 κ + in D (k) and the fact that κ i , i ∈ N + , depend on the mesh norms (on the area of the respective triangles in (18) for d = 3 and the length of line segments in (19) for d = 2). Hence, cond (S (k) ) tends to infinity, when the finite element mesh norm approaches zero. As S (k) is ill-conditioned through its diagonal, we use the diagonal preconditioner: where D F (k) = diag(F (k) ) and D D (k) = diag(0, D −1 κ + D(I + + ∪ I − + ), 0, diag(E)). From spectral analysis presented in [34], it follows that although S (k) and P (k) depend on κ i , i ∈ N + , the eigenvalue bounds of the preconditioned matrix (P (k) ) −1 S (k) are positive and independent of κ i , i ∈ N + . Consequently, cond ((P (k) ) −1 S (k) ) ≤ cond (D F (k) )cond (F (k) ). (42) In computations, we approximate D F (k) using analogous idea as in ALGORITHM PF.

Numerical Experiments
The computations were performed by the supercomputer Salomon at IT4I VŠB-TUO [40]. The Salomon cluster consists of 1009 compute nodes. Each node is a powerful x86-64 computer with Intel Xeon E5-2680v3 processors equipped with 24 cores and at least 128 GB RAM. All codes are implemented in Matlab R2020a. The velocity component is eliminated in both algorithms by solving auxiliary linear systems involving A κ or A with the preliminary Cholesky factorization realized by the Matlab function chol. We use ALGORITHM PF with c 1 = 10 −3 , c 2 = 10 9 , c min = 10 −12 , c max = 0.5, c 3 = 10 −2 , ε = 10 −3 , r tol = 0.9, c fact = 0.9. These values turned out to be optimal, as it follows from the tests in [34]. ALGORITHM SSN uses r tol = 0.01, c fact = 0.5, and ε = 10 −3 . The termination tolerance ε leads to the relative residua of order 10 −5 . In the tables below the numbers n it , n F of the outer iterations, and the matrix-vector multiplications by F κ or F (k) , respectively, are monitored. Note that n F determines overall complexities of computations. The partitions T h of Ω are generated by Iso2mesh toolbox and ANSYS software [41,42]. As we have already mentioned, the finite element spaces use P1-bubble/P1 element pairs on T h . The resulting mesh will be characterized by values of the parameters n u , n p , n s introduced in Section 3. The stiffness matrices are assembled by the vectorized code [33,43]. To guarantee the uniqueness of the solution to (4) we will consider that ∂Ω is decomposed into three non-overlapping and non-empty parts Γ, S, and Γ N open in ∂Ω (see Remark 2). Nevertheless, one example and comments when Γ N = ∅ will be presented, as well. Finally note that all physical quantities are considered in the SI system so that the units of g and κ are Pa and Pa.s.m −1 , respectively.
Note that u exp satisfies the homogeneous Dirichlet boundary condition on ∂Ω. Data are chosen in such a way that both, leak and no leak zones appear on S. Figure 1 (right) presents the streamlines of the computed solution. One can see that the leak on S generates the fluid suction on Γ N . Figure 2 shows the pressure and velocity distribution in Ω. The magnitudes of the normal velocity |u n | and of the normal stress |σ n | on S are seen in Figure 3 (left), while | σ n | is depicted in Figure 3 (right). In the tables below the values n it , n F , and CPU time (in seconds) for different finite element meshes with increasing n u , n p , and n s , and different values of g are compared. The computational complexities for g = 15 are shown in Table 1. The analogous characteristics for g = 0.1 (leak everywhere) and g = 100 (no leak) are summarized in Tables 2 and 3, respectively. In view of the values of n F , one can see that both algorithms are scalable. The computations without preconditioning are considerably less efficient as it is seen from Table 4.      Let Ω = (0, 1) 3 . To construct the partition T h we first cut Ω into small cubes and then each cube is split into five tetrahedras, see Figure 4. The partition of ∂Ω into Γ, S, and Γ N is defined as follows: Again, u exp satisfies the homogeneous Dirichlet boundary condition on ∂Ω and data are chosen in such a way that both, leak and no leak zones appear on S. Figure 5 shows the distribution of the velocity field on S and Γ N . One can see the leak zone in the middle of S and the fluid suction on γ left . The same suction effect appears on γ right . The normal stress |σ n | and | σ n | on S are seen in Figure 6. In Tables 5-8 we summarize analogous characteristics as in Example 1. The conclusions are also analogous, i.e., both algorithms turn out to be scalable also in 3D and the preconditioner plays an important role. Notice that the numbers n F of the matrix-vector multiplications are lower compared to 2D case.       . We change the partition of ∂Ω from Example 2 as follows: Γ N = ∅, Γ = γ left ∪ γ right ∪ γ back ∪ γ front ∪ γ top , and S = γ bottom . Let us remind that, if Γ N = ∅, the pressure might be non-unique depending on the choice of g.
Next we suppose that g is constant and κ = 0 on S. To establish the critical value of g we use the Stokes system with the homogeneous Dirichlet data prescribed on the whole boundary ∂Ω. Assuming that its solution is smooth enough, g crit = (sup S (σ n ) − inf S (σ n ))/2, where σ n is the corresponding normal stress. However since the whole procedure is applied to the finite element discretization, the regularity assumption is superfluous. For given physical data and the partition T h characterized by the parameters n u /n p /n s = 38, 088/15, 625/529 the value g crit = 18.31. It is easy to show that if g < g crit , the problem has a unique solution with both, leak and no-leak zones on S. If g > g crit , the solution is not unique and only no-leak zone appears on S. The behavior of the algorithms for different g is summarized in Table 9. One can see that both algorithms are stable in all cases. The table presents also the values sup S (σ g n ) and inf S (σ g n ) computed for the respective g. One can easily check that (sup S (σ g n ) − inf S (σ g n ))/2 = g crit for g ≥ g crit , while (sup S (σ g n ) − inf S (σ g n ))/2 = g if g < g crit . The distribution of the velocity field on S for g = 16.31 is seen from Figure 7. Note that the matrix C used in (25), especially its part (N T , B T ) T , has not the full row-rank due to Γ N = ∅. Consequently, F κ and F (k) for A 0 = N 0 are singular with the defect one and the eigenvector whose all components are equal to 1 corresponds to the zero eigenvalue. The matrix S (k) κ of the inner linear systems in ALGORITHM PF is regularized by adding the diagonal matrix J 12 M (k) (Z (k) ) −1 J T 12 to F κ as seen from (31) (there are non-zero diagonal entries on the level of N). Analogously, the matrix S (k) in ALGORITHM SSN is regularized by adding D (k) to F (k) (see (41)) and by the choice of the active set A 0 = N 0 . These heuristic arguments show that the matrices of the solved linear systems can be regularized by the construction of the algorithms. Finally note that the condition number bounds of the preconditioned Schur complements (32) and (42) are no longer valid. Nevertheless, the preconditioners can still be used and play the same important role as in the previous examples.

Example 4. (3D branched tube).
Let Ω be the branched tube as in Figure 8, where the partition T h and the partition of the boundary ∂Ω = S ∪ γ in ∪ γ out1 ∪ γ out2 is also depicted. This tube may represent a final part of the capillary in the human body. On Γ = γ in we prescribe the input velocity u |Γ,1 = a − a/r 2 (y 2 + z 2 ), u |Γ,2 = u |Γ,3 = 0 with a = 5 · 10 −4 , r = 5.94 · 10 −4 while Γ N = γ out1 ∪ γ out2 consists of two outputs of the tube with the natural outflow condition σ N = 0. The leak condition is prescribed on S = ∂Ω \ Γ ∪ Γ N . Further f = 0 in Ω, µ = 1/2, and κ = 30. The partition T h is characterized by the following values of the parameters: n u = 27, 408, n p = 9178, and n s = 4079. The behavior of ALGORITHM PF and ALGORITHM SSN for different g is reported in Table 10. The distribution of the velocity and the pressure for three selected g is depicted in Figure 9. Almost all liquid leak through the wall at the vicinity of γ in for the smallest value g = 10 can be observed, see Figure 9 (top). The leak through the wall is minimal for g = 50, see Figure 9 (bottom). The pressure is also growing for higher values of g. The zoomed parts of the velocity fields near γ in and γ out2 are depicted in Figure 10 (the velocity distribution inside the tube looks like the Poiseuille distribution). We used the reorthogonalization technique of the descent vectors inside the CGM solver [44] that increases considerably efficiency of computations; compare Table 10 with Table 11.

Conclusions and Comments
We have proposed two algorithms for solving the Stokes flow with the threshold leak boundary condition, both based on dual strategies. Unlike the threshold slip boundary conditions [28], the algorithms are identical in 2D and 3D. ALGORITHM PF is based on the path-following variant of the interior point method that generates strictly feasible iterations. This algorithm combines the dumped Newton method in the outer loop with the preconditioned conjugate gradient method used for solving the inner linear systems. The precision control of the inner loop is driven adaptively with respect to the accuracy achieved in the outer level. The crucial ingredient of the algorithm is its diagonal preconditioner that improves ill-conditioning of the system matrices in the later iterations of the outer loop.
ALGORITHM SSN uses the active-set implementation of the semi-smooth Newton method adapted for our type of boundary conditions. The inner linear systems are solved using the same ideas as in ALGORITHM PF. In this case, the diagonal preconditioner improves ill-conditioning caused by small finite element mesh norms. Experimental scalability of both algorithms is indicated by numerical experiments. Indeed, one can conclude that ALGORITHM SSN is more efficient than ALGORITHM PF. Problems with multiple solutions can be also solved since the singular dual Hessian is implicitly regularized by the construction of the algorithms.