Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Advances in the simulation of viscoplastic fluid flows using interior-point methods

Keywords: yield stress fluids, viscoplasticity, Bingham model, interior point method, conic programming

We present a primal-dual interior point algorithm for the resolution of steady-state viscoplastic fluid flows formulated as a conic optimization problem. We give a complete description of the algorithm including some advanced aspects such as a predictor-corrector and scaling scheme to improve its efficiency. Our interior-point approach is shown to be largely more efficient than Augmented Lagrangian (AL) approaches which are traditionally used to solve such problems. In particular, the interior-point approach is roughly 5 times faster than the modern accelerated version of AL algorithms. The yield surfaces are shown to be accurately predicted and various examples ranging from channel flows to three-dimensional flows through a porous medium demonstrate its efficiency.

Introduction

Viscoplastic or yield-stress fluids are materials which behave like a solid below a critical yield stress and flow like a viscous fluid for stresses higher than this threshold. Such materials are now commonly found in civil engineering, petroleum or food industries. Although many models have been proposed to describe the behavior of such materials, the Bingham model is the most well-known and very simple as the material is considered to be rigid below a yield criterion described by von Mises criterion and is linearly viscous above the yield stress. Such a simple model is enough to reproduce interesting physics such as the existence of stopped and rigid-body plug regions in flows or the return to rest in a finite time. For recent reviews on yield stress fluids, we refer to [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF][START_REF] Coussot | Bingham's heritage[END_REF].

Regarding their numerical simulation, viscoplastic fluids are very challenging because of the separation between flowing and non flowing regimes. This non-smoothness in the expression of the constitutive relation necessitates to resort to more advanced techniques than those traditionally used for Newtonian fluids for instance. Regularized models have first been proposed to replace the non-smooth viscoplastic constitutive law by a smooth purely viscous model [START_REF] Bercovier | A finite-element method for incompressible non-Newtonian flows[END_REF][START_REF] Papanastasiou | Flows of materials with yield[END_REF]. Augmented Lagrangian (AL) approaches have then emerged as an interesting alternative to the use of regularized models to solve viscoplastic fluid flows, it is now one of the most popular method to solve such problems [START_REF] Fortin | Méthodes de lagrangien augmenté: applications à la résolution numérique de problèmes aux limites[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Saramito | An adaptive finite element method for viscoplastic fluid flows in pipes[END_REF][START_REF] Dean | On the numerical simulation of Bingham visco-plastic flow: old and new results[END_REF][START_REF] Glowinski | On the numerical simulation of viscoplastic fluid flow[END_REF]. However, it still suffers from a slow convergence rate so that three-dimensional simulations are still extremely expensive. Recently, it has been pointed out that such AL schemes can be accelerated quite easily, yielding a faster convergence rate [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF]. For a recent review on the simulation of viscoplastic fluids and various extensions, we refer to [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF] and references herein.

In [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF], it has first been proposed to use interior-point algorithms to solve the velocity minimum principle satisfied by the solution to a viscoplastic fluid flow. It follows the general trend that has been observed in the field of computational limit analysis/yield design theory in which interior-point methods (IPM) have been largely and successfully used for estimating the collapse load of various structures (continuum, beams, plates, shells, etc.) made of different materials (steel, masonry, soils, etc.). It is interesting to note that this field experienced a similar evolution regarding numerical techniques used to solve the rigid-perfectly plastic problem, which can be explained by the strong similarities shared with viscoplastic fluids. Regularization methods have indeed first been proposed [START_REF] Friaâ | Le matériau de Norton-Hoff généralisé et ses applications en analyse limite[END_REF][START_REF] Guennouni | Calcul à la rupture, régularisation de Norton-Hoff et Lagrangien augmenté[END_REF] until Augmented Lagrangian techniques became more appropriate [START_REF] Fortin | Méthodes de lagrangien augmenté: applications à la résolution numérique de problèmes aux limites[END_REF][START_REF] Vicente Da Silva | A non-linear programming method approach for upper bound limit analysis[END_REF][START_REF] Carlier | A projection approach to the numerical analysis of limit load problems[END_REF]. Since the development of the IPM in the mathematical programming community [START_REF] Wright | Primal-dual interior-point methods[END_REF], linear and non-linear convex optimization problems can now be efficiently solved and it became the state-of-the-art method in the field of computational limit analysis [START_REF] Andersen | Computing limit loads by minimizing a sum of norms[END_REF][START_REF] Lyamin | Upper bound limit analysis using linear finite elements and non-linear programming[END_REF][START_REF] Krabbenhoft | A general non-linear optimization algorithm for lower bound limit analysis[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and secondorder cone programming[END_REF].

More precisely, interior-point algorithms are very efficient for second-order cone programming (SOCP), a specific class of problems which can be reformulated using second-order cone constraints. It encompasses linear programming problems and many non-linear problems, such as quadratic programming (see [START_REF] Lobo | Applications of second-order cone programming[END_REF] for different illustrative applications). In the field of limit analysis, many yield criteria can be reformulated using such SOCP constraints and, in particular, the von Mises criterion. In [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF], it has been shown that the quadratic viscous term for the Bingham model could also be easily reformulated so that the associated optimization problems also falls in the SOCP framework. It was therefore proposed to use the industrial interior-point software Mosek [START_REF] Mosek | The Mosek optimization software[END_REF] as a black-box solver for computing the viscoplastic flows. To be solved by a dedicated solver such as Mosek, SOCP problems often need to be formulated in a standard format. For instance, the objective function can only be linear, which requires reformulating the quadratic viscous term using second-order cone constraints. Another requirement of the standard format is that a given optimization variable can belong to one cone only. For these reasons, the resulting optimization problems present a much larger number of variables or constraints compared to the initial formulation, leading to higher computational costs. The extension to a Herschel-Bulkley power-law behavior also introduces similar problems.

In the present paper, we follow the general idea of the IPM, as usually implemented in IPM solvers, to derive an efficient algorithm for computing viscoplastic flows. The work presented here will serve many purposes:

• give all the necessary information to implement a home-made interior point algorithm and share the associated code for a larger diffusion of IPM in the field of viscoplastic fluids (part of the codes used in this paper can be found at https://doi.org/10. 5281/zenodo.1038519 [START_REF] Bleyer | Viscoplastic flows : supplementary code for "advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF]);

• further assess its efficiency compared to Augmented Lagrangian algorithms;

• further validate the solutions computed with the IPM, in particular regarding the accuracy of yield surfaces;

• provide an efficient way to resolve the limitations introduced by standard formats of black-box solvers.

To answer all these goals, the manuscript will first rapidly recall Augmented Lagrangian approaches and their accelerated versions in section 2, since they will be used for numerical comparison with IPM results. In section 3, the SOCP framework and the general idea of the IPM will be given. Advanced aspects of the IPM will be detailed in section 4 which will lead to a much more efficient algorithm. To ease the notation, the presentation of the first sections will be made in the context of antiplane problems with a Bingham model, 2D and 3D problems as well as extensions to other constitutive models will be discussed in section 5. Finally, numerical results and illustrative applications are presented in section 6.

Augmented Lagrangian approach for viscoplastic fluids

In this section, we will first recall the Augmented Lagrangian formulation which is now widely used for computing viscoplastic flows along with accelerated versions of the corresponding optimization algorithms.

Augmented Lagrangian approach

Let us consider the following minimum principle satisfied by the steady-state velocity field u(x) for the antiplane flow through a cylindrical channel of section Ω lying in the Oxy-frame and subjected to a pressure gradient f (x):

min u∈V Ω η 2 ∇u 2 + τ 0 ∇u - Ω f u (1)
where η represents the viscosity and τ 0 the yield stress of the Bingham model. We assume no-slip boundary conditions on some part ∂Ω u of the domain boundary and V denotes the appropriate functional space for u.

By introducing the strain rate d as an additional optimization variable, the non-smoothness of the yield stress term and the derivative arising from the gradient operator are decoupled in the objective function and treated as a linear equality constraint:

min u,d Ω η 2 d 2 + τ 0 d - Ω f u (2a) s.t. ∇u = d (2b)
The previous formulation fits in the general framework of minimizing the sum of two convex functions G(x) + H(y) constrained by a linear equation of the form Ax + By = b:

min x,y G(x) + H(y) (3a) s.t. Ax + By = b (3b)
In our case, with x = d and y = u, we have G(x) = Ω g(x) with g(x) = η 2 x 2 + τ 0 x and H(y) = -f, y +ι V (y) with ι V being the indicator function of V, while A = Id, B = -∇ and b = 0. Two common methods for solving problems of the form (3) are the Alternating Direction Method of Multipliers (ADMM) and the Alternating Minimization Algorithm (AMA) which are Augmented Lagrangian approaches combined with an Uzawa-like algorithm. Both of them perform two separate minimization steps, once in x and once in y, before updating the Lagrangian multiplier associated with the linear constraint: ADMM

x k+1 = arg min x G(x) -λ k , Ax + r 2 Ax + By k -b 2 (4a)
y k+1 = arg min y H(y) -λ k , By + r 2 Ax k+1 + By -b 2 (4b)
λ k+1 = λ k + r(b -Ax k+1 -By k+1 ) (4c)
Whereas for AMA, the first step is performed without the quadratic regularization term if G is strictly convex (which is the case here since η > 0):

AMA x k+1 = arg min x G(x) -λ k , Ax (5a)
y k+1 = arg min y H(y) -λ k , By + r 2 Ax k+1 + By -b 2 (5b) λ k+1 = λ k + r(b -Ax k+1 -By k+1 ) (5c)
where in both cases r > 0 is the augmentation parameter of the AL approach.

Identifying the Lagrange multiplier λ k with a stress vector σ k , the first step amounts to the following update of the strain rate variable:

d k+1 =    σ k η 1 - τ 0 σ k if σ k ≥ τ 0 0 otherwise (6)
with

σ k = σ k + r∇u k , η = η + r for ADMM σ k = σ k , η = η for AMA (7)
The second step amounts, for both algorithms, to solving the following Stokes problem: Find u k+1 ∈ V such that:

Ω r∇u k+1 • ∇v = Ω f v - Ω (σ k -rd k+1 ) • ∇v ∀v ∈ V (8)
Finally, the Lagrange multiplier update reduces, for both algorithms, to:

σ k+1 = σ k + r(∇u k+1 -d k+1 ) (9)
Note that it is possible to invert the first two steps in both ADMM and AMA yielding slightly different algorithms with very similar efficiency.

Accelerated versions of Augmented Lagrangian methods

Both ADMM and AMA algorithms require solving the same Stokes problem at each iteration and perform local updates only. The cost per iteration is hence quite small since the Stokes stiffness matrix can be factorized once and for all at the beginning of the iterations. Despite this nice feature, the number of iterations required to reach convergence has often been reported to be very large, the theoretical convergence rate of the residuals being in O(1/k).

However, it is possible to accelerate the convergence of both algorithms to reach an optimal O(1/k 2 ) rate of the residuals [START_REF] Goldstein | Fast alternating direction optimization methods[END_REF]. Most importantly, these accelerated versions are obtained at a completely negligible cost, except for keeping in memory the fields from the previous iteration. Following Nesterov's predictor-corrector scheme [START_REF] Goldstein | Fast alternating direction optimization methods[END_REF][START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o (1/k2)[END_REF], the previous steps ( 6)-( 9) are now replaced by:

d k+1 =    σ k η 1 - τ 0 σ k if σ k ≥ τ 0 0 otherwise (10)
with

σ k = σ k + r∇u k , η = η + r for accelerated-ADMM σ k = σ k , η = η for accelerated-AMA (11)
Find u pred k+1 ∈ V such that:

Ω r∇u pred k+1 • ∇v = Ω f v - Ω (σ k -rd k+1 ) • ∇v ∀v ∈ V (12)
σ pred k+1 = σ k + r(d k+1 -∇u k+1 ) (13)
t k+1 = 1 + 1 + 4t 2 k /2 (14a)
u k+1 = u pred k+1 + t k -1 t k+1 (u pred k+1 -u pred k ) ( 14b
)
σ k+1 = σ pred k+1 + t k -1 t k+1 (σ pred k+1 -σ pred k ) ( 14c
)
where t 0 = 1 and u pred 0 = u 0 , σ pred 0 = σ 0 for the initialization.

It is clear that the first predictor steps are exactly the same as before. The only difference comes from the corrector step ( 14) which changes the previous predictions for the next iteration using some inertial terms accounting for the values obtained at the last iteration. Thus, contrary to the unaccelerated scheme, we must keep in memory the predictors of the last iteration k to perform the correction at iteration k + 1. We also observe that an additional scalar parameter t has been introduced which corresponds to an adaptively chosen relaxation step. We recover the unaccelerated version by fixing t k = 1.

Remark:. Convergence of the accelerated ADMM version is ensured only if both functions G and H are strongly convex [START_REF] Goldstein | Fast alternating direction optimization methods[END_REF]. In the present case, H is only weakly convex and a restart procedure would need to be implemented to ensure convergence. This has not been done here since we always observed convergence of the accelerated ADMM algorithm, at least for antiplane simulations. Convergence of the accelerated AMA algorithm is however always ensured [START_REF] Goldstein | Fast alternating direction optimization methods[END_REF].

The use of such accelerated schemes in the context of viscoplastic fluid flows has recently been proposed in [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF]. We share the point of view of this paper to use such accelerated variants in the future when implementing Augmented Lagrangian algorithms because of the small amount of work needed to accelerate an already existing AL code.

Primal-dual interior point method for viscoplastic fluids

We now describe the primal-dual interior point approach proposed for computing viscoplastic flows with a much higher efficiency than AL approaches. This section concentrates on the main ideas and most important ingredients of the method. More advanced aspects to obtain an efficient algorithm will be discussed in section 4.

Second-order cone programming formulation and optimality conditions

The interior-point approach can be applied to any linear or non-linear convex optimization problem but has been shown to be particularly efficient for solving specific classes of convex programming: linear programming (LP), second-order cone programming1 (SOCP) and semi-definite programming (SDP). It usually requires reformulating the considered problem as the minimization of a linear function under specific convex constraints: linear equality/inequality constraints for LP, second-order cone constraints for SOCP and semi-definite constraints for SDP. As shown in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF], the viscoplastic problem can be rewritten into a second-order cone program by introducing additional optimization variables and conic constraints to fit the standard formulation required by industrial IP softwares like Mosek. However, we will take advantage of the particular structure of the problem, namely an objective function which can be decomposed into a smooth (viscous) and a non-smooth (plastic) part to reformulate only the non-smooth part using second-order cone constraints and keep the non-linear but smooth part of the objective function as such. This leads to an important saving on the number of auxiliary variables compared to what has been proposed in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF].

By characterizing the subgradient of the objective function in [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF], its solution u satisfies the following conditions [START_REF] Dean | On the numerical simulation of Bingham visco-plastic flow: old and new results[END_REF][START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF]:

Find (u, t, λ) such that:

Ω (η∇u + τ 0 λ) • ∇v = Ω f v ∀v ∈ V (15a) ∇u = tλ , t(1 -λ ) = 0, t ≥ 0, λ ≤ 1 (15b)
which expresses the mechanical equilibrium for a stress state of the form σ = η∇u + τ 0 λ, that is a sum of a viscous part and a stress state τ 0 λ which is inside the yield criterion or at its boundary. Besides, the strain rate is associated to the plastic part by the flow rule ∇u = tλ with either t = 0 if λ < 1 or λ = 1 if t > 0.

Conditions (15b) can also be equivalently reformulated as the following primal-dual second order cone constraints:

∇u = tλ, t -∇u • λ = 0, ∇u ≤ t, λ ≤ 1 ⇐⇒ (t, ∇u) • (1, -λ) = (0, 0), (t, ∇u) ∈ Q 3 , (1, -λ) ∈ Q 3 ( 16
)
where the circular product is defined in Appendix A and with Q n denoting the quadratic Lorentz cone of dimension n. The Karush-Kuhn-Tucker (KKT) optimality conditions can now be rewritten in the form of linear and conic constraints: Find (u, t, λ) such that:

Ω (η∇u + τ 0 λ) • ∇v = Ω f v ∀v ∈ V (17a) (t, ∇u) • (1, -λ) = (0, 0) (17b) (t, ∇u) ∈ Q 3 , (1, -λ) ∈ Q 3 (17c)
Remark:. The generalization of the above derivation to Herschel-Bulkley fluids is immediate, one just has to replace the viscous quadratic term in (1) by the corresponding power-law potential, transforming the linear variational problem (17a) into a non-linear one.

General idea

The idea of the interior-point method (IPM) is to find a solution to (17) by following, more or less closely, a curve called the central path given by (u(µ), t(µ), λ(µ)) and parametrized by a barrier parameter µ ≥ 0. The central path is defined as the unique solution to the following perturbation of the KKT system (17) [START_REF] Megiddo | Pathways to the optimal set in linear programming[END_REF][START_REF] Alizadeh | Second-order cone programming[END_REF]:

Find (u, t, λ) such that :

Ω (η∇u + τ 0 λ) • ∇v = Ω f v ∀v ∈ V (18a) (t, ∇u) • (1, -λ) = (µ, 0) (18b) (t, ∇u) ∈ Q 3 , (1, -λ) ∈ Q 3 (18c)
For µ = 0, we have exactly [START_REF] Wright | Primal-dual interior-point methods[END_REF]. From the perturbed equation (18b), we obtain for the solution to (18) that:

t(µ) = µ 2 + 1 2 µ 2 + 4 ∇u 2 > ∇u (19a) λ(µ) = ∇u µ 2 + 1 2 µ 2 + 4 ∇u 2 -1 (19b)
so that the stress state is given by σ

= η∇u + τ 0 ∇u µ 2 + 1 2 µ 2 + 4 ∇u 2 -1
. As a consequence the solution of (18) for a given µ > 0 can be interpreted as the solution of a regularized viscoplastic behavior, as pointed out recently in [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF]. However, we will later discuss how the proposed method is fundamentally different, at least from the user point of view, from traditional regularized approaches.

A nice property of the central path is that it defines a continuous set of strictly feasible points which are well-centred, i.e. far from the boundary of the feasible region (except when approaching the optimum for very small µ) so that large descent steps can be taken when minimizing the objective function from points located near the central path (see Figure 1). Indeed, when computing descent directions from the unperturbed KKT equations ( 17), the iterates usually attain the feasible region boundary very rapidly so that only very small steps can be taken close to the boundary to preserve feasibility of the subsequent iterates, yielding a small convergence rate. If one is able to compute solutions to [START_REF] Andersen | Computing limit loads by minimizing a sum of norms[END_REF], the central path defines a preferential way to converge to the solution of the viscoplastic problem by driving µ to 0.

From the KKT conditions (18b), different strategies for finding a solution to [START_REF] Andersen | Computing limit loads by minimizing a sum of norms[END_REF] using (19b), working only on the primal variable u. The auxiliary variable t plays no role anymore and the second-order cone constraints are automatically satisfied. In this case, the interior-point method shares strong similarities with the barrier method and regularization approaches. We will not detail primal interior point methods here since their practical performance is not very good compared to primal-dual methods which will be described later. For such an approach, λ and t remain as unknowns and ( 18) is then as difficult to solve as [START_REF] Wright | Primal-dual interior-point methods[END_REF]. However, as will be made clear in section 3.4, we will never aim at finding points exactly on the central path but just use it as a guide towards the solution.

Finite element discretization

Before going into the general steps of the primal-dual IPM, let us first introduce the finite-element discretization of problem [START_REF] Wright | Primal-dual interior-point methods[END_REF]. We consider a linear interpolation of the fluid velocity field u on triangular finite elements. The strain rate vector, the auxiliary variable t associated with its norm as well as the dual normalized stress vector λ are all interpolated as piecewise constants on the same triangular mesh which will enable to satisfy everywhere all the constraints of ( 17). We adopt here the same mixed stress-velocity finite element approximation as traditionally used in Augmented Lagrangian approaches. Denoting the corresponding vectors of unknowns using bold symbols, the discrete version of the perturbed KKT conditions are given by: Find (u, t, d, λ) such that:

Ku + τ 0 B T λ = f (20a) Bu -d = 0 (20b) (t i , d i ) • (1, -λ i ) = µe ∀i = 1, . . . , N (20c) (t i , d i ) ∈ Q 3 , (1, -λ i ) ∈ Q 3 (20d)
where e = (1, 0, 0), B = ∇N u denotes the velocity shape function derivatives and K is the stiffness matrix arising from the viscosity term. As is classically done in standard interiorpoint conic solvers, an additional strain rate variable d has been introduced to avoid any linear relationship inside the cones so that they are all decoupled from each other whereas coupling is restored by the additional linear constraint Bud = 0. The conic constraints and complementarity equations are written for all elements i = 1, . . . , N with N being the total number of elements.

Primal-dual interior point algorithm

We are now ready to describe the different steps of the primal-dual IPM. As mentioned earlier, we do not need to find exactly points lying on the central path but will only try to remain close to it. Hence, at each iteration k, a Newton step on the perturbed system [START_REF] Krabbenhoft | A general non-linear optimization algorithm for lower bound limit analysis[END_REF] is computed towards the central path for a fixed value of µ k . From this new point, the barrier parameter is decreased by some amount and the process is repeated until the different residuals fall below a prescribed tolerance. This method generates a sequence of iterates (u k , t k , d k , λ k ) which need to remain strictly feasible regarding the conic constraints i.e. d k i < t k and λ k i < 1.

Solving the perturbed KKT system

Let us consider a given iteration k at the beginning of which a point z k = (u k , t k , d k , λ k ) satisfying the conic constraints is known and a value for the barrier parameter µ (k) has been chosen.The next iterate is obtained by computing an approximation of the point on the central path for µ = µ (k) . This approximation is obtained by using one iteration of the Newton method. Hence, a descent direction ∆z = (∆u, ∆t, ∆d, ∆λ) is computed by solving the linearization of the perturbed KKT system around z k :

A∆z =     K 0 0 τ 0 B T B 0 -I 0 0 I -C k λ -C k d 0 -diag(λ k ) I -C k t            ∆u ∆t ∆d ∆λ        =        r k d r k p r k c r k c        = r k (21)
where

C k λ = diag(. . . , (λ k i ) T , . . .), C k d = diag(. . . , (d k i ) T , . . .) and C k t = diag . . . , t k i t k i , . . .
arise from the linearization of (20c).

The vector r k corresponds to a vector of residuals for the different linearized equations of the KKT system:

r k d = f -Ku k -τ 0 B T λ k is the dual residual, r k p = d k -Bu k is the primal residual and r k c = µ (k) -t k + C k λ d k , r k c = C k t λ k -d k
are the centering residuals. Due to the linearity of the first two equations of [START_REF] Krabbenhoft | A general non-linear optimization algorithm for lower bound limit analysis[END_REF], it is possible that the primal and dual residuals are always zero if the initial point z 0 = (u 0 , t 0 , d 0 , λ 0 ) is chosen so that r 0 p = r 0 d = 0. If finding an initial point satisfying the linear constraints for a general interiorpoint problem may be difficult, in our case a simple initialization using the solutions of the corresponding Stokes flow and λ 0 = 0 satisfies these conditions. The case where r k p = r k d = 0 for all iterations is usually referred to as a feasible interior-point method in contrast to an infeasible interior-point method where the initial point does not satisfy the primal and dual residual equations. However, one has to recall that, in any case, each iterate will always need to be feasible with respect to the conic constraints. Although it may seem more efficient to use a feasible approach, numerical experiments show that infeasible approaches are almost as efficient so that the overall performance is quite insensitive to the choice of the initial point.

Step-length computation

The next iterate is given by z k+1 = z k + α∆z where α ∈ (0, 1] is a step-length which is chosen so that z k+1 still satisfies all conic constraints. In practice, one can adopt for instance a simple line-search procedure. A fraction-to-boundary rule is used to avoid that the next iterates lie exactly on the boundary of the feasible domain. In practice, this rule consists in choosing α = θα max with θ = 0.99 -0.995.

The primal and dual residuals are then reduced by a factor 1 -α. The efficiency of the interior-point approach is due to the fact that this step-length will in general remain quite large i.e. close to 1, because the iterates tend to follow the central path which is located far from the boundaries of the conic-feasible region. On the contrary, a traditional descent direction along the gradient of the objective function will reach the feasible region boundaries at early stages of the iterations and the subsequent step lengths will be very small in order to maintain feasibility, exhibiting then only little progress.

Choice of the barrier parameter

A key point of the interior-point method lies in the choice of the barrier parameter µ for each iteration. In order to solve the initial unperturbed problem, the method aims at decreasing this parameter from one iteration to another. At the end of iteration k, a feasible point z k+1 has been computed and a common choice for µ (k+1) makes use of the value of the complementarity gap

g k+1 i = t k+1 i -λ k+1 i • d k+1 i for each point i: µ (k+1) = γ 1 N N i=1 g k+1 i = γg k+1 (22)
i.e. µ (k+1) is equal to the global average complementarity gap g k+1 times a scalar parameter γ ∈ (0, 1] which will control the amount of centering for the next iteration. The idea of this centering parameter is that for γ = 0 the KKT equations of the next iteration will aim at reducing the complementarity gap but maybe to the detriment of going away from the central path. On the contrary, a value of γ = 1 means that the complementarity gap will not be reduced (see Appendix B) but the next iterate will be closer to the central path so that a larger step may be taken at the next iteration. We will see in the next section that this centering parameter can be chosen in an adaptive manner.

Initialization and convergence criteria

For the initialization of the algorithm, we chose to take u = 0 (except when non-zero imposed velocity boundary conditions are present, see example 6.2), d = λ = 0 and t = 1. This choice ensures that the initial point is strictly feasible regarding both conic constraints (t i , d i ) ∈ Q 3 and (1, λ i ) ∈ Q 3 . Besides, it also satisfies the strain rate-velocity relation Bu = d. As a consequence, the primal residual r p will always be zero during the iterative process.

Another possible initialization would be to use the solution u Stokes of the corresponding Stokes problem as a first guess for u and Bu Stokes as a first guess for d. In this case, t has to be chosen to be strictly greater than Bu Stokes,i everywhere. We did not see a significant improvement in the efficiency of the algorithm by choosing this initialization so we decided to keep the first choice.

For the chosen initialization, the initial complementarity gap is equal to 1 everywhere. As the residual vectors and the complementarity gap are both reduced at the same rate, we stopped the iterative process when the following conditions are met:

g k < and (r k p , r k d ) < (23)
with a tolerance of typically = 10 -8 . A point satisfying these conditions will then be considered as an optimal solution to the initial unperturbed problem [START_REF] Wright | Primal-dual interior-point methods[END_REF].

The algorithm is also stopped if the iteration step size becomes too small or if a maximum number of iterations (e.g. N iter,max = 200) is exceeded meaning that the algorithm failed to find an optimal solution. In all our numerical tests, this situation has never been encountered.

A simple version of the primal-dual interior point algorithm is summarized by Algorithm 1.

Advanced aspects of the primal-dual interior-point method

Predictor-corrector scheme and adaptive choice of the barrier parameter

One major difficulty in using the previous simple version of the primal-dual IPM lies in the optimal choice of the centering parameter γ. As mentioned before, the choice of γ represents a trade-off between improving centering and reducing the complementarity gap. One major improvement in the efficiency of primal-dual IPM has been realized by Mehrotra [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF] who proposed an adaptive choice of the centering parameter by using a predictor-corrector scheme. This scheme takes advantage of the fact that factorizing the KKT matrix A of equation ( 21) is the most demanding step for each iteration. Once this factorization has been found, it is quite inexpensive to solve a second linear system with the same matrix.

Algorithm 1 A simple primal-dual interior point algorithm 1: Initialization (k = 0): choose z 0 = (0, 1, 0, 0, 0), γ ∈ (0, 1] and > 0 2: (r 0 p , r 0 d ) = (f , 0)

3: g 0 = 1 N N i=1 (t 0 i -λ 0 i • d 0 i ) = 1 4: µ (0) = γg 0 5: while (r k p , r k d ) > or g k > do 6:
form residual vector r k and Newton matrix A if k > N iter,max or α ≤ then break "Algorithm failed."

13:

end if 14: end while More precisely, the predictor-corrector scheme uses the fact that:

(t + i , d + i ) • (1, -λ + i ) = (t i , d i ) • (1, -λ i ) +(∆t i -λ T i ∆d i -d T i ∆λ i , ∆d i -t i ∆λ i -∆t i λ i ) -(∆d T i ∆λ i , ∆t i ∆λ i )
with the notation z + = z + ∆z. When linearizing this quantity to form the Newton step of the perturbed system, the last quadratic term (∆d T i ∆λ i , ∆t i ∆λ i ) is neglected. Mehrotra suggested estimating this term from a first pure Newton step on the unperturbed system. Hence, ( 21) is first solved with µ (k) = 0, the corresponding right-hand side being given by:

r k d = f -Ku k -τ 0 B T λ k (24a) r k p = d k -Bu k (24b) r k c = -t k + C k λ d k (24c) r k c = C t k λ k -d k (24d)
This first step is called the affine step and its solution is denoted by ∆z a = (∆u a , ∆t a , ∆d a , ∆λ a ). We denote by α a max the maximum step length possible for this search direction. The centering parameter γ is then estimated from this maximum step length using the following heuristic [START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF]:

γ = (1 -α a max ) • min 0.5; (1 -α a max ) 2 (25)
This heuristic enables to choose a centering parameter which is small if a large affine step can be taken (α a max close to 1) so that the adopted step is close to a pure Newton step. Otherwise, centering is improved by choosing a larger value of γ, yielding a more conservative direction towards the central path. The final search direction is then computed as the solution ∆z of [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and secondorder cone programming[END_REF] with the following residuals:

r k d = (1 -γ) f -Ku k -τ 0 B T λ k (26a) r k p = (1 -γ) d k -Bu k (26b) r k c = µ (k) -t k + C k λ d k (26c) r k c = C t k λ k -d k + C k ∆t a ∆λ a (26d)
with µ (k) = γg k and where the quadratic correction C k ∆t a ∆λ a has been added to the centering residuals based on the affine step prediction (the quadratic correction term of (26c) is in fact zero, see Appendix B). We remark that a factor (1 -γ) has been added to the primal and dual residuals, this is motivated by the fact that the new primal and dual residuals corresponding to the point z + α∆z are then given by (r k+1 p ,

r k+1 d ) = (1 -α(1 -γ))(r k p , r k d )
and that the new average complementarity gap is also given by g k+1 = (1 -α(1 -γ))g k (see again Appendix B). This ensures that each quantity is decreased at the same rate which is beneficial for the algorithm efficiency.

Reduced linear system

As previously mentioned, the main computation cost of the IPM is the resolution of the Newton linear system (21) (the factorization of A in case of a predictor-corrector scheme) at each iteration. It can be observed from the structure of system (21) that A is a sparse non-symmetric matrix and is quite large compared to the size of the stiffness matrix K of a classical Stokes problem. However, it is possible to reduce the previous system to a smaller system, of the size of K, to solve (21) more efficiently, although the obtained matrix will be non-symmetric matrix if condensing matrix A directly. A specific rescaling procedure has, therefore, to be adopted to obtain a symmetric reduced system but also to ensure convergence of the interior-point algorithm.

To do so, let us first rewrite the complementarity conditions as follows:

(t i , d i ) • (1, -λ i ) = X i s i = X i S i e = S i X i e = S i x i (27a)
where

x i = t i d i (27b)
s i = 1 -λ i (27c)
X i = mat(x i ) = t i (d i ) T d i t i I (27d)
S i = mat(s i ) = 1 -(λ i ) T -λ i I (27e)
which are standard notations for conic programming problems. Denoting X = diag(. . . , X i , . . .) and S = diag(. . . , S i , . . .) the corresponding global block-diagonal matrices, equation (18b) can be rewritten as XS e = µ e where e is the global concatenation of N e vectors. The linearization of (27a) is then rewritten as:

S∆x + X∆s = S ∆t ∆d + X 0 -∆λ (28)
We now introduce a symmetrization of the complementarity equations and use one of the most famous called the Nesterov-Todd scaling, described in Appendix C, which instead of considering the consistency equation x i •s i = 0 replaces it with v i •v i = 0 where v i is a scaled point depending on x i and s i and is such that

(x i ) T s i = v i 2 (but with x i • s i = v i • v i ).
In particular, there exists a symmetric matrix F i such that v i = (F i ) -1 s i = F i x i . The linearized complementarity equation ( 28) is then replaced by the following:

VF∆x + VF -1 ∆s = R k c = µ e -V 2 e ( 29
)
where V (resp. F) is the block matrix made of each

V i = mat(v i ) (resp. F i ).
We now rewrite the linearized KKT system (21) as follows:

A∆z =   K 0 τ 0 B T J B J 0 0 VF VF -1      ∆u ∆x ∆s    =    r d r p R c    (30)
where we dropped the iteration counter k and with J = diag(. . . , [0 -I 2 ], . . .) so that Js = -λ. Since V is an invertible matrix for a strictly feasible primal-dual point (x, s) (see Appendix C), the last two rows can be easily condensed, using only pointwise computations at the finite element level, to form the following final reduced system (see Appendix D):

K + τ 0 B T (F -2 ) -1 B A ∆u = r d + τ 0 B T (F -2 ) -1 (r p + F -1 V -1 R c ) = r (31a)
Therefore, the reduced linear system form amounts only to solving the linear system (31a) A ∆u = r which is still sparse but is now much smaller and symmetric positive definite since (F -2 ) -1 is (see (C.18)). The corresponding matrix A can be interpreted as the sum of a purely viscous contribution with constant viscosity K and a locally variable contribution due to the yield stress term τ 0 B T (F -2 ) -1 B.

Besides, Mehrotra's predictor-corrector scheme can still be applied with this reduced form since the linear system matrix A, and therefore A , does not change between the affine and the corrector step. One just has to adapt the corresponding values of the reduced residual vectors for each step. The quadratic correction term for the second step of Mehrotra's scheme is here given by F i ∆x a i • F -1 i ∆s a i where ∆x a , ∆s a are the directions solutions to the affine step.

Finally, let us highlight that computing the reduced residual vectors, the reduced matrix A or the back-substitutions for the auxiliary variables involves the computation of F -1 , V -1 or (F -2 ) -1 which are given by closed-form formulae in Appendix C.

Summary of the final algorithm

We now summarize the different steps of the final primal-dual interior point algorithm (Algorithm 2) incorporating Mehrotra's predictor-corrector scheme and a reduced linear system form.

Algorithm 2 Predictor-corrector primal-dual interior point algorithm 1: Initialization (k = 0): choose z 0 = (0, 1, 0, 0, 0) and > 0 2: (r 0 p , r 0 d ) = (f , 0)

3: g 0 = 1 N N i=1 (t 0 i -λ 0 i • d 0 i ) = 1 4: while (r k p , r k d ) > or g k > do 5:
compute scaling matrices F, V from the current point z k 6:

form reduced Newton matrix A from (31a)

7:

form affine residual vector r a corresponding to the reduced version of (24)

8:

solve affine step A ∆u a = r a 9:

back-substitute for (∆t a , ∆d a , ∆λ a ) from (D.1b),(D.3b)

10:
compute maximum affine step length α a max 11:

estimate centering parameter γ through (25)

12:

set µ (k) = γg k 13:

form corrector residual vector r c corresponding to the reduced version of (26) if k > N iter,max or α ≤ then break "Algorithm failed."

20:

end if

21:

k ←-k + 1 22: end while

Numerical implementation

The numerical implementation of the proposed IPM algorithm for viscoplastic fluid flows has been realized using FEniCS [START_REF] Logg | Automated solution of differential equations by the finite element method: The FEniCS book[END_REF] as regards the finite-element discretization and building of the KKT linear systems. Thanks to the flexibility offered by FEniCS in the choice of the linear system solver, we used the direct solver MUMPS [START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF][START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF] which enables to solve large sparse linear systems. Computations have been realized on a Intel Core i7-4600U (2.1 GHz).

Comments on some extensions

Two and three-dimensional flows

Let us quickly consider here the case of two or three-dimensional flows, the overall steps of the IPM are exactly the same as for the antiplane case with only small differences mainly due to the incompressibility constraint. Indeed, the solution of a steady viscoplastic flow in a two or three-dimensional domain Ω can be obtained from the following minimum principle:

min u Ω ηd : d + τ 0 2d : d - Ω f • u (32a) s.t. div u = 0 (32b)
where d = ∇ s u = (∇u + ∇u T )/2 is the strain rate tensor. Following the same procedure as in Section (3), this minimum principle gives rise to the following optimality conditions: Find (u, p, t, λ) such that:

Ω (2η∇ s u + τ 0 λ) : ∇ s v + p div v = Ω f • v ∀v ∈ V (33a) div u = 0 (33b) (t, d) • (1, -λ) = (0, 0) (33c) (t, d) ∈ Q m , (1, -λ) ∈ Q m (33d)
with m = 4 in 2D, m = 7 in 3D 2 and where we introduced a scaled vectorial representation of the strain rate tensor components d ij as follows: [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF] and similarly for the tensor λ. The discretized KKT conditions of the previous optimality conditions are given by: Find (u, p, t, d, λ) such that: where D represents the divergence operator and p is the pressure-like Lagrange multiplier associated with the incompressibility condition. We can then observe that the top-left block of the Newton matrix A corresponding to the viscosity stiffness matrix K is now replaced by a mixed displacement-pressure system of the form K D T D 0 with a new primal residual vector r p associated with the incompressibility constraint. Similarly, for the reduced linear system formulation described in section 4.2, the procedure is exactly the same except that the pressure variable p cannot be eliminated so that the final system has the following symmetric mixed form:

                                   d = Bu =      √ 2d xx √ 2d yy 2d xy      in 2D d = Bu =                    √ 2d xx √ 2d yy √ 2d zz 2d xy 2d xz 2d yz                    in 3D (
Ku + D T p + τ 0 B T λ = f (35a) Du = 0 (35b) Bu -d = 0 (35c) (t i , d i ) • (1, -λ i ) = 0 ∀i (35d) (t i , d i ) ∈ Q m , (1, -λ i ) ∈ Q m (35e)
K + τ 0 B T (F -2 ) -1 B D T D 0 ∆u ∆p = r r p (36)
Finally, let us also mention that a traditional mixed Taylor-Hood P 2 -P 1 approximation for (u, p) is chosen while all strain/stress-like variables (t, d, λ) are discretized using discontinuous linear elements. As a consequence each conic constraint is verified at the three vertices of each triangular element which is enough to satisfy it everywhere from convexity arguments.

Extension to other viscoplastic models

It was previously mentioned that second-order cone programs need to be formulated in a specific standard format to be solved by the interior-point solver Mosek. In particular, the objective function has to be reformulated as a linear objective function by introducing additional variables and second-order cone constraints. If this reformulation is quite straightforward in the case of a quadratic viscous potential for a Bingham fluid, it is less obvious in the Herschel-Bulkley case. A specific procedure has been described in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF] to tackle the Herschel-Bulkley case with a rational power-law exponent. In addition to increasing the number of optimization variables, one drawback of such an approach is that the form of the reformulation is strongly dependent on the value of the exponent and has to be performed on a case-by-case basis.

In the proposed formulation, the quadratic viscous part of the potential has not been reformulated using cone constraints but kept as such in the objective function. It then appears through its constant Hessian matrix K when forming the linearized KKT system. Therefore, we propose to deal with the Herschel-Bulkley case in the exact same way. The K stiffness matrix will then be replaced, at each iteration, by the current tangent stiffness matrix related to the Herschel-Bulkley power-law behavior. The iterative process will then combine a traditional Newton method on the nonlinear viscous term as well as on the conic constraints. We do not expect a strong increase in the number of iterations to converge compared to the Bingham case since the power-law viscous behavior is smooth but this remains to be tested in the future.

More generally, the same splitting strategy could be used to any yield stress fluid constitutive models (Casson fluids for example) which can be split into a smooth and yield stress non-smooth part.

Illustrative applications

Channel flow in an eccentric annulus

As a first validation of the proposed approach, we consider the problem of the antiplane flow in a cylindrical channel, the section of which consists of an eccentric annulus [START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF][START_REF] Wachs | Numerical simulation of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier/fictitious domain and augmented Lagrangian methods[END_REF]. The outer circle has a normalized radius R e = 1 while we fixed the inner circle radius to R i = 0.4 (Figure 2). The center of the inner circle is located at a distance δ from that of the outer circle on the middle plane y = 0. For symmetry reasons, only the upper half of the section y ≥ 0 has been modelled. The flow is driven by a uniform unit pressure gradient f = 1, the fluid viscosity is also taken as unity η = 1 and a no-slip velocity condition is assumed on both circle boundaries.

Except stated otherwise, we used the predictor-corrector interior-point algorithm with a reduced linear system (Algorithm 2). To assess its efficiency, its computational performance has been compared to both a standard AL approach and its accelerated version as described in section 2, the following results have been obtained using the ADMM algorithm, although the performance of AMA is very similar. Choosing the augmentation parameter r is a delicate point of AL algorithms and an optimal value is often problem-dependent. In [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF], a link between such algorithms and accelerated proximal algorithms has been drawn provided that r is chosen as the Lipschitz constant of G(x) in (3), which corresponds to r = η in the present case. This choice also ensures convergence of the accelerated algorithm. For this reason, this choice has been retained in the subsequent simulations for both standard and accelerated variants. Since the IPM is a second-order Newton-like approach, the number of iterations is expected to be orders of magnitude lower than first-order approaches like AL but also more expensive to compute. For this reason, the comparison between these approaches will be made in terms of CPU time. As regards both AL approaches, one attractive feature is that at each iteration, the same Stokes-like problem needs to be solved with different right-hand side Lagrangian approach (N iter > 1000: did not converge to the prescribed tolerance), its accelerated version (N iter = 986) and the proposed interior-point method (N iter = 16). The accelerated AL exhibits a faster convergence rate than the standard AL but is still outperformed in terms of computational cost by the proposed IPM approach.

members. In our FEniCS implementation, we fully exploit this feature by factorizing the corresponding matrix once and reusing it at each iteration, similarly to what is done for the predictor-corrector IP algorithm. The cost of one AL iteration is therefore much smaller. The residual measure for AL is given by ∇u k -d k whereas for the IPM the evolution of the average complementarity gap g k is reported. The results are represented in Figure 3 for the annulus problem (δ = 0.04) with a fine mesh (N e = 66, 077 elements) and for a yield stress τ 0 = 0.1. It can clearly be observed that the accelerated AL algorithm exhibits a faster convergence than the standard AL algorithm. It reaches the tolerance of 10 -8 in 986 iterations whereas the standard AL algorithm fails to reduce the residual below the prescribed tolerance in less than 1000 iterations. However, the IP algorithm converges to the optimal solution in less than 30s i.e. approximately 4 times faster than the accelerated AL. The average cost of one AL (either standard or accelerated) iteration is around 0.12s against 1.77s for the IPM. CPU times for other mesh sizes have also been reported in The influence of the mesh size has been investigated in Figure 4 and shows that its influence with respect to the number of iterations required to reach the prescribed tolerance is quite weak for the IPM. In particular, it remains of the order of a few tens of iterations as generally reported in practice for state-of-the-art interior-point solvers. As regards AL, the mesh size has some influence on the standard version whereas it is more difficult to assess on the accelerated version since it is not a descent algorithm. We also verify that both AL algorithms reached their theoretical convergence rate.

In Figure 5, the simple variant of the IPM with a constant centering parameter γ (Algorithm 1) has been compared to Algorithm 2. For the comparison to be fair, steps 6 and 7 in Algorithm 1 are modified to solve the reduced linear system A ∆u = r as described in section 4.2 instead of the global KKT system, we then only investigate the influence of the adaptive choice of the centering parameter through the predictor-corrector scheme. It can clearly be observed that the predictor-corrector variant converges with much fewer iterations (16 against 50 to 100 depending on γ) but also in much less computation time than Algorithm 1. The average computing cost of one iteration for Algorithm 1 is around 1.3 s for this specific example whereas the cost of one iteration of Algorithm 2 is around 1.77 s. Hence, the cost of the corrector step in Algorithm 2 induces an overhead of only 36% per iteration but with an important saving in the number of iterations. Besides, Figure 5 also shows that an optimal value of γ may be difficult to find since the step length of each iteration can be very small for small values of γ despite a fast reduction of the complementarity gap. Conversely for larger values of γ, the complementarity gap is slow to decrease despite large step lengths. On the contrary, the adaptive heuristic (25) works very well in general with a fast reduction of the gap and sufficiently large step lengths so that Algorithm 2 does not require fine tuning of the algorithm parameters.

In [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF], it was mentioned that "replicating the yield surface positions found in the computations of [START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF] is a worthy challenge for the hardened regulariser! ". We therefore checked that the IPM is able to replicate such yield surfaces as shown in Figure 6 (see the results reported in [START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF][START_REF] Wachs | Numerical simulation of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier/fictitious domain and augmented Lagrangian methods[END_REF]), notice in particular the sharp re-entrant corner which are difficult to obtain with regularized models. These yield surfaces are identified as the regions for which σ ≤ τ 0 where we recall that σ = η∇u + τ 0 λ, thus combining information from both the primal variable u and the dual variable λ. At convergence, the strain rate norm inside the yield surface is below 10 -8 which is the order of the prescribed convergence tolerance. One could then argue that the strain rate is never zero in unyielded zones but it could easily be post-processed by projecting the strain rate to zero in unyielded regions. The difference between these two strain rate states would induce only a minor difference in the residual norms and would still satisfy the convergence tolerance. Hence, from a user perspective, the IPM does not seem to exhibit any of the typical difficulties of regularized approaches as the complementarity gap (and hence the barrier parameter µ) has been driven to zero, within the prescribed numerical tolerance.

For the δ = -0.15 case, we computed the flow rate Q = Ω udS of the optimal velocity field. The results are reported in Figure 7 in terms of a non-dimensional flow rate Q * = Q/Q 0 where Q 0 corresponds to the flow rate of a Newtonian fluid in a concentric annulus δ = 0 (see [START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF]) and for varying pressure gradient f and yield stress τ 0 . In the Newtonian case, we obtain a linear relationship as expected whereas for non-zero yield stress, no flow occurs below a critical pressure gradient, the value of which increases with the yield stress. The computed flow rate Q * in these situations is not strictly zero but is also around the order of the convergence tolerance. For high pressure gradients, the response is again linear with the same slope as for the Newtonian case.

Two-dimensional flow: the lid-driven cavity

We consider here the typical benchmark of the lid-driven cavity problem: a unit-square cavity with an applied velocity u = e x on the top surface y = 1, the other boundaries being subjected to a no-slip velocity condition. We define the Bingham number as Bi = τ 0 L/ηU , with here U = L = 1, as the non-dimensional number measuring the influence of yield stress effects over viscous effects. The initial pressure field for the starting point of the interiorpoint algorithm has been chosen to be zero whereas the velocity field is initialized to zero for all degrees of freedom except those located on the top surface on which the imposed velocity is prescribed. During the iterative process, the Newton directions ∆u are then sought as corrections kinematically admissible to zero i.e. ∆u i = 0 everywhere on the cavity boundary 3 . For this initialization, the strain rate primal residual r p is non zero since the initial zero strain rate is not compatible with this initial velocity field. In fact, all residuals are here non zero at the beginning.

We first report the convergence history of the different residuals with a structured discretization of element size h = 1/50 (Figure 8). It can be observed that all residuals, either the primal (r p for the strain rate and r p for the pressure), the dual r d or the average complementarity gap g, converge to zero at the same rate as previously remarked. When increasing the Bingham number to a high value, we observe a smaller overall convergence rate, thus needing a larger number of iterations to reach the prescribed tolerance. However, this increase remains moderate, the IPM needing only a few tens of iterations to converge. The same trend is observed when refining the mesh. On the contrary, the number of iterations for AL algorithms seems more dependent on the increase of the Bingham number.

During the iteration process, we observe that the iso-stress contour lines converge to a well-defined location. Therefore, no problem regarding the lack of point-wise convergence for stresses seem to appear with the IPM contrary to what has been observed for regularized models [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF][START_REF] Putz | On the lubrication paradox and the use of regularisation methods for lubrication flows[END_REF]. The location of the yielded regions are reported in Figure 9 for different Bingham numbers and are consistent with results from the literature (see for instance [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF][START_REF] Syrakos | Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method[END_REF]).

Three-dimensional flow in a porous medium

As a final illustration, we investigate a three-dimensional flow around a fixed bodycentered cubic stacking of rigid spheres representing a porous medium. The unit cell of this periodic structure is represented in Figure 10. Taking into account the various symmetries, only one eighth of the model is considered with appropriate boundary conditions. The flow is driven by a uniform downwards pressure gradient f = -f e z and no-slip velocity conditions are assumed on the sphere boundaries. The yield stress is here τ 0 = 10, the viscosity is still η = 1 and the sphere radius has been taken as R = 0.4.

The computed velocity field along with unyielded points have been represented in Figure 11 for f = 140 and in Figure 12 for f = 300. The first case corresponds to a situation for which the pressure gradient is slightly above the critical value for which flow occurs whereas in the second case, the pressure gradient is more than twice its critical value. It can be noted that, as expected, the flow is mainly concentrated along paths which are far from the spheres and that many points near the obstacles are still unyielded. For the larger pressure gradient, the flow spreads in the central part of the computational domain around the obstacles and less unyielded points are observed. When looking at the flow along a 2D slice of equation x = y, in the most tortuous part of the unit cell, we also observe a secondary flow at the bottom of the computational domain for the larger pressure gradient. The evolution of the downwards flow rate Q z = Ω -u z dV over the computational domain has been represented in Figure 13 as well as the horizontal flow rate Q x = Ω u x dV (also equal to Q y because of symmetry). We obtain the same overall behavior as in Figure 7 for the annulus channel, that is a non-flowing part below a critical pressure gradient and an almost linear response at high velocity gradients. As regards computing times, for a mesh made of approximately 50 000 tetrahedra, the complete optimization required about 15-20 minutes. For a mesh made of 110 000 tetrahedra, the computing time increased to slightly more than 2 hours.

Conclusions and perspectives

In this work, we tackled the problem of computing viscoplastic fluid flows by solving the associated velocity minimum principle, reformulated using second-order cone constraints.

29

We presented an interior-point algorithm to efficiently solve the associated convex optimization problem, without relying on black-box commercial solvers as proposed in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF]. In particular, it enables to reduce the number of auxiliary variables needed to formulate the problem in a standard form. It also paves the way to tackling more complex viscoplastic models by taking directly into account the smooth non-linear viscous part in the objective problem, which was not possible in the commercial software used previously. The detailed steps of the primal-dual IPM have been described and more advanced aspects such as the predictor-corrector scheme and the reduced linear system form ensuring a good efficiency have been discussed. The illustration of the algorithm on various numerical examples ranging from channel flows to 2D and 3D flows enabled to show that:

• the proposed IPM algorithm is very efficient in terms of total computing time compared to standard as well as accelerated AL algorithms. The typical number of iterations is moderate as usually reported for such methods. It seems also less sensitive to mesh size and Bingham number than AL approaches;

• although it is difficult to compare it precisely, the total computing time of our homemade IPM algorithm seems also quite comparable to what we experienced using the commercial IPM solver Mosek on similar problems;

• the yield surfaces obtained using this method are in agreement with those obtained using AL approaches and results from the literature;

• although the IPM can be interpreted as a Newton method on a regularized problem with a continuation of the regularization parameter towards zero, we did not observe any of the reported problems of regularization approaches with the IPM. From a user perspective, the obtained solution is close to the one corresponding to the original unperturbed problem within the prescribed numerical tolerance.

Although the efficiency of IPM in the context of viscoplastic flows is very encouraging, some work is still needed to reach large-scale three-dimensional simulations. The main bottleneck that needs to be tackled first to achieve this objective concerns the linear solver for the Newton system. Although MUMPS is an efficient direct solver, which enabled to solve moderate 3D problems on a standard computer, and can be used in parallel, the memory cost of the LU decomposition will be prohibitive for very large systems. The use of iterative solvers instead of direct solvers for the Newton system will however need the development of good preconditioners, the main problem being that the condition number of A deteriorates when approaching the solution (x and s are approaching the boundary of the second-order cone). Recent works regarding preconditioning techniques for viscoplastic flows can be found in [START_REF] Aposporidis | Multigrid preconditioning of the non-regularized augmented Bingham fluid problem[END_REF][START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF].

In the present work, only steady-state flows have been considered. The computation of time-dependent problems is, in theory, quite immediate as the time-discretization adds a quadratic "mass" term in the objective function and changes the loading vector by taking into account inertial terms from the last iteration. Time-dependent problems therefore require solving successively similar problems. However, it is known that IPM cannot easily take advantage of a previously known solution as a good initial guess ("warm-start") as it is too close from the feasible region boundary and far from the new central path. Finding good warm-starting strategies for IPM is still an active research topic [START_REF] Yildirim | Warm-start strategies in interior-point methods for linear programming[END_REF][START_REF] Gondzio | Reoptimization with the primal-dual interior point method[END_REF][START_REF] John | Implementation of warm-start strategies in interior-point methods for linear programming in fixed dimension[END_REF].

Finally, adaptive mesh refinement strategies have been demonstrated to increase the accuracy of the computed solutions by resolving the transition between flowing and nonflowing regions [START_REF] Saramito | An adaptive finite element method for viscoplastic fluid flows in pipes[END_REF][START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF]. The coupling of IPM with auto-adaptive meshes should not pose any difficulty.

If K i are Lorentz cones, the complementary slackness condition x T s = 0 is equivalent to x i • s i = 0 ∀i where :

x • s =
x 0 s 0 + xT s x 0 s + s 0 x (A.4) with x = (x 0 , x), the notation representing the tail of a vector i.e. the vector obtained by removing the first element (here noted x 0 ).

Appendix B. Results on the solutions to the Newton equations

We consider the solutions to the following Newton equations with a residual vector parametrized by γ, κ ∈ [0; 1]:

    K 0 0 τ 0 B T B 0 -I 0 0 I -C k λ -C k d 0 -diag(λ k ) I -C k t            ∆u ∆t ∆d ∆λ        =        κ(f -Ku k -τ 0 B T λ k ) κ(d k -Bu k ) γg k -t k + C k λ d k C k t λ k -d k        (B.1)
where we recall that g k is the average complementarity gap given by g k = 1 where z + = z k + α∆z and g + is the new average complementarity gap associated with z + .

Proof. The first two relations are easily obtained due to the linearity of the associated residuals. Besides, (B.2d) follows directly from (B.2c) and the third line of (B.1). Finally the proof to (B.2c) can be easily adapted from the one found in [START_REF] Xu | A simplified homogeneous and self-dual linear programming algorithm and its implementation[END_REF].

Let us remark that if κ = 1, (B.1) reduces to [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and secondorder cone programming[END_REF]. From the results of the previous lemma, we observe that the primal and dual residuals decrease by 1 -α whereas the average complementarity gap decreases by 1 -α(1 -γ(1 -α)). This justifies the fact that the complementarity gap is not reduced a lot if γ is close to 1 (see the discussion in section 3.4). Hence, a primal-dual IPM aims at using a small value of γ but, at the same time, small values of γ also usually lead to small step lengths α (see Figure 5). Similarly, choosing κ = 1 -γ yields that both residuals and complementarity gap are reduced at the same 1 -α(1 -γ) rate so that the algorithm is globally convergent if a strictly positive step α can be taken at each iteration.

In the predictor-corrector framework described in section 4, the first affine step corresponds to (B.1) with γ = 0 and κ = 1. Hence, the associated directions satisfy (∆d a i ) T ∆λ a i = 0 for all i and a purely affine step of length α a max would reduce residuals and complementarity gap by 1 -α a max . This justifies the choice of γ close to 1 if α a max is small for the corrector step in order to improve centering. Finally, because (∆d a i ) T ∆λ a i = 0, adding the quadratic term correction to the residuals of [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o (1/k2)[END_REF] preserves the results of the previous lemma for the corrector step with κ = 1 -γ.

Appendix C. Nesterov-Todd scaling for the reduced linear system This section gives general results on Jordan algebra and the Nesterov-Todd scaling for second-order cones. Further details can be found in [START_REF] Alizadeh | Second-order cone programming[END_REF][START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF][START_REF] Nesterov | Self-scaled barriers and interior-point methods for convex programming[END_REF][START_REF] Yamashita | A primal-dual interior point method for nonlinear optimization over secondorder cones[END_REF].

For a vector v = (v 0 , v) ∈ Q m+1 , we define:

mat(v) = v 0 vT v v 0 I m ∈ R m+1×m+1 (C.1) det(v) = det (mat(v)) = (v 0 ) 2 -v 2 (C.2) v = Qv = (v 0 , -v) (C.3) (mat(v)) -1 = 1 det(v)   v 0 -v T -v det(v)I m + vv T v 0   (C.4)
with I m being the identity matrix of size m and Q = diag(1, -I m ). Note that the last relation is valid only if v is located strictly inside Q m+1 , i.e. v ∈ int(Q m+1 ), which is equivalent to det(v) > 0.

Let us consider here the complementary slackness condition:

x • s = x 0 s 0 + xT s x 0 s + s 0 x = Xs = Sx = XSe (C.5)

where X = mat(x), S = mat(s) and e = (1, 0). We further assume that x, s ∈ int(Q m+1 ). It can then be shown that there exists a unique symmetric matrix F (depending on x and s) such that: Fx = x = s = F -1 s and x T s = x T s (C.6)

The scaled point v = x = s satisfies the following properties:

x ∈ Q m ⇔ v ∈ Q m (C.7) x ∈ int(Q m ) ⇔ v ∈ int(Q m ) (C.8)
where denotes the block matrix/vector obtained by keeping only the tail of each corresponding block (e.g. ∆s = -∆λ). This last system can the be further condensed into:

K + τ 0 B T (F -2 ) -1 B A ∆u = r d + τ 0 B T (F -2 ) -1 (r p + F -1 V -1 R c ) = r (D.3a) ∆λ = (F -2 ) -1 B∆u -r p -F -1 V -1 R c (D.3b)

Figure 1 :

 1 Figure1: General idea of the interior point method. The central path defines a set of well-centered points located far from the feasible region boundaries. A descent step rapidly ends near the feasible region boundary. By using steps which follow the central path, larger steps can be taken to avoid the feasible region boundary.

7 :solve A∆z = r k 8 :maximum step length α 9 :

 789 compute update variables z k+1 = z k + α∆z 10: compute new barrier parameter µ (k+1) = γg k+1 11: k ←-k + 1 12:

14 : 16 :maximum step length α 17 :

 141617 solve corrector step A ∆u = r c (reuse factorization of A ) 15: back-substitute for (∆t, ∆d, ∆λ) from (D.1b),(D.3b) compute update variables z k+1 = z k + α∆z 18: compute new complementarity gap g k+1 19:

Figure 2 :

 2 Figure 2: Geometry of the eccentric annulus channel section.

Figure 3 :

 3 Figure 3: Comparison of CPU time for the annulus problem (δ = 0.04) between the standard AugmentedLagrangian approach (N iter > 1000: did not converge to the prescribed tolerance), its accelerated version (N iter = 986) and the proposed interior-point method (N iter = 16). The accelerated AL exhibits a faster convergence rate than the standard AL but is still outperformed in terms of computational cost by the proposed IPM approach.

Figure 4 :

 4 Figure 4: The mesh size N e does not really change the number of iterations for convergence for IPM as well as AL. The accelerated AL version exhibits a O(1/k 2 ) convergence rate compared to the standard O(1/k) rate. The IPM requires, in practice, only a few tens of iterations to converge.

Figure 5 :

 5 Figure5: Left: Evolution of the complementarity gap g for the annulus problem (δ = 0.04) for the simple IPM algorithm with a constant centering parameter γ and for the adaptive predictor-corrector variant. The number of iterations and the total computing times are much lower with the adaptive variant than for the constant γ variant. Right: Step length of each iteration for the corresponding cases.

Figure 6 :

 6 Figure 6: Optimal velocity field and unyielded regions (in gray) for the annulus problem with δ = 0.04 (left) and δ = -0.15 (right). The unyielded region shape and location are consistent with those obtained with AL approaches and results from [35, 36].

Figure 7 :

 7 Figure 7: Non-dimensional flow rate Q * = Q/Q 0 as a function of the applied pressure gradient f for different yield stress τ 0 . No flow occurs below a critical pressure gradient whereas the response becomes linear at high pressure gradients with the same slope as the Newtonian case τ 0 = 0.

  (a) Bi = 2 (b) Bi = 200

Figure 8 :

 8 Figure 8: Convergence of the different IPM residuals for the lid-driven cavity problem (h = 1/50) and two different Bingham numbers. Each type of residual clearly converges at the same rate. The number of iterations to convergence slightly increases for high Bingham numbers.

Figure 9 :

 9 Figure 9: Velocity streamlines and unyielded elements (in gray) for the lid-driven cavity problem (h = 1/200) with varying Bingham numbers.

Figure 10 :

 10 Figure 10: Three-dimensional unit cell of the porous medium (R = 0.4) along with the computational domain taking into account the various symmetries.

Figure 11 :

 11 Figure 11: Velocity field for a pressure gradient f = 140 slightly above the critical value (red dots denote unyielded regions).

  2D slice in x = y plane

Figure 12 :

 12 Figure 12: Velocity field for a pressure gradient f = 300 much higher then the critical value (red dots denote unyielded regions).

Figure 13 :

 13 Figure 13: Vertical and horizontal flow rates over the computational domain as a function of the pressure gradient.

  then have the following result: Lemma Appendix B.1. The directions ∆z solutions to (B.1) are such that:f -Ku + -τ 0 B T λ + = (1 -ακ)(f -Ku k -τ 0 B T λ k ) (B.2a) d + -Bu + = (1 -ακ)(d k -Bu k ) (B.2b) ∆d T i ∆λ i = κ(1 -γ -κ)N g k (B.2c) g + = (1 -α(1 -γ) + α 2 κ(1 -γ -κ))g k (B.2d)

Table 1 .

 1 It can be observed that the IPM is 4 to 6 times faster than the accelerated AL for all mesh sizes.

	Mesh size	Accelerated AL (s) IPM (s) Speed-up factor
	N e = 4 092	13.9	2.2	6.3
	N e = 16 492	34.3	6.4	5.4
	N e = 66 077	117.6	28.4	4.1
	N e = 264 230	758.6	120.3	6.3

Table 1 :

 1 CPU times for different meshes and speed-up factor of the IPM over accelerated AL. Standard AL never converged in less than 1000 iterations.

which encompasses quadratic programming

Note that it would be possible instead to add the imposed velocity condition u = e x as an additional linear constraint to the optimization problem, the residual of which would also converge to zero during the IPM iterative process.

Appendix A. Quadratic second-order cone programming Let us consider here the following quadratic second-order cone program:

where P is symmetric positive definite and K is a tensorial product of cones x ∈ K ⇔ x i ∈ K i with K i being either:

• the full space of dimension n i : R n i (i.e. x i is a vector of free variables)

• the positive orthant: R n i

+

• the Lorentz second-order cone:

This problem is not a standard SOCP program because of the quadratic objective term, it can however be reformulated as a problem involving a linear objective function only by introducing auxiliary variables and additional conic constraints [START_REF] Lobo | Applications of second-order cone programming[END_REF][START_REF] Alizadeh | Second-order cone programming[END_REF].

The Lagrangian L of the previous program is given by:

and the corresponding optimality conditions by [START_REF] Boyd | Convex optimization[END_REF]:

where K * is the tensorial product of the dual cones

The linearization of the complementarity equation is thus rewritten around point v as follows:

In the case of second-order cones, this transformation is explicitly given by the following symmetric matrix:

(C.14)

We then have the following results:

From the last relation, we deduce that F -2 = θ -2 (I m + 2 w wT ) where is the tail i.e. the bottom-right m × m block from matrix . Applying the Sherman-Morrison formula, we obtain that:

which is easily checked to be symmetric positive definite.

Appendix D. Condensation of the KKT system

From system (30), elimination of the last row gives:

Now, since ∆s i = 0 -∆λ i for all i, the previous system is equivalent to:


Emmanuel Schang
email: emmanuel.schang@univ-orleans.fr

Extended Projections in a Guadeloupean TAG Grammar

Introduction

This paper presents a TAG (meta-)grammar of Guadeloupean Creole (GC). Guadeloupean is a French-based creole language spoken on the island of Guadeloupe (French West Indies). While sharing most of its lexicon with French, GC differs from French in its grammar (see [START_REF] Bernabé | Fondal-natal[END_REF] for a comprehensive description, and [START_REF] Damoiseau | Syntaxe créole comparée[END_REF] for a comparative approach). In particular, GC has preverbal Tense-Aspect Markers while French has inflectional suffixes for Tense. I will show that these TMAs can be correctly described as extended projections of verbs. They are generated as inflected forms of a verbal lexeme (i.e. in morphology rather than in syntax).

First, I present succinctly the key concepts of a TAG grammar to the reader. Next, section 2 describes the Tense-Aspect markers (TMAs) in GC and section 3 presents the way they are incorporated in the GC metagrammar1 . Section 4 briefly presents how the same approach can be applied to the nominal domain.

Before addressing the question of Tense and Aspect markers, let me briefly present the key notions of TAG.

Tree-Adjoining Grammar (TAG) allows two operations (for a complete description, see [START_REF] Joshi | Tree-Adjoining Grammars[END_REF], from which I take the following definitions):

-" Substitution takes only place on non-terminal nodes of the frontier of a tree. [...] By convention, the nodes on which substitution is allowed are marked by a down arrow (↓). When substitution occurs on a node n, the node is replaced by the tree to be substituted. When a node is marked for substitution, only trees derived from initial trees can be substituted for it." (Joshi and Schabes, 1997, p.4) -" Adjoining2 builds a new tree from an auxiliary tree β and a tree α (α is any tree, initial, auxiliary or derived). Let α be a tree containing a nonsubstitution node n labeled by X and let β be an auxiliary tree whose root node is also labeled by X. The resulting tree γ, obtained by adjoining β to α at node n is built as follow:

• the sub-tree of α dominated by n, call it t, is excised, leaving a copy of n behind. • the auxiliary tree β is attached at the copy of n and its root node is identified with the copy of n. • the sub-tree t is attached to the foot node of β and the root node of t (i.e. n) is identified with the foot node of β. " (Joshi and Schabes, 1997, p.4) As proposed in [START_REF] Frank | Phrase Structure Composition and Syntactic Dependencies[END_REF], Substitution and Adjunction are supposed to be universal operations. The differences between languages can only reside in the shape of the elementary trees. The way elementary trees are built is thus a crucial matter, and yet the locus of divergences between several TAG grammars. It is thus necessary to explain the principles that govern the building of the elementary trees.

The major reference for a TAG French grammar is [START_REF] Abeillé | Une Grammaire électronique du Français[END_REF] 3 which presents the following (linguistic) principles of elementary trees well-formedness:

Lexical Anchoring: An elementary tree must have (at least) one non-empty lexical head. Predicate-Argument Co-occurrence: A predicate elementary tree must have a node for each of its arguments. Semantic Anchoring: A syntactic elementary tree must correspond to a (non-empty) semantic element. Compositionality Principle: An elementary tree corresponds to one and only one semantic unit.

Furthermore, I adopt the Conditions on Elementary Tree Minimality (CETM) (Frank, 2002, 54) :

CETM: The syntactic heads in an elementary tree and their projections must form an extended projection of a single lexical head.

This leads me to adapt [START_REF] Grimshaw | Locality and extended projection[END_REF]'s definition of head and projection to the TAG framework and say that an elementary tree is a maximal projection of one (or several) lexical head(s) within which the categorial features are shared. This opens discussions about what is functional and what is lexical. In the case of prepositions, this is a difficult matter for which the debate remains open (see [START_REF] Cinque | The syntax of adjectives: a comparative study[END_REF]).

2 Tense-Aspect Markers 2.1 A Brief Description of Tense-Aspect Markers in GC Creole languages are known to make use of independent4 markers to express Tense and Aspect (see [START_REF] Winford | Creole Languages. The Oxford Handbook of Tense and Aspect[END_REF] for a synthesis) usually gathered under the label TMA markers (Tense, Mood and Aspect markers) or TAM (Tense and Aspect Markers). Guadeloupean behaves like other creoles w.r.t. TMAs. The inflectional morphology found in French verbs is replaced in GC by independent morphemes. For comprehensive analyses, the reader may refer to [START_REF] Mccrindle | Temps, mode et aspect, les creoles des Caraibes a base lexicale franc ¸aise[END_REF] and [START_REF] Pfänder | Describing São Tomense Using a Tree-Adjoining Meta-Grammar[END_REF]. In Table 1, I will use the description of [START_REF] Vaillant | A layered grammar model: Using tree-adjoining grammars to build a common syntactic kernel for related dialects[END_REF], taking as a paradigmatic example the verb dansé 'to dance'5 6 . [START_REF] Vaillant | Une grammaire formelle du créole martiniquais pour la génération automatique[END_REF] also notes that ké ka and té ké ka are attested, although rare7 . The precise semantic value and uses of these TMAs are beyond the scope of this paper and the reader may refer to [START_REF] Pfänder | Describing São Tomense Using a Tree-Adjoining Meta-Grammar[END_REF] for a complete description. The following lines are just intended to present the useful key points for a non-specialist reader.

Bare Verbs: As in many languages, bare verbs in GC are used to express the (past) perfective (or preterite) with dynamic processes (as in (1-a)) and express the present tense with stative verbs8 , as in (1-b).

(1) a. Jan Jean month 'He has a lot of money at the end of each month.' (example from (Delumeau, 2006, 117))

With non-stative predicates, the favored interpretation is Imperfective:

(5) Jan Jean ka ASP manjé. eat 'Jean is eating' ké triggers a future reading:

(6) Jan Jean ké ASP manjé. eat 'Jean will eat'

Projecting TMAs

TMAs in GC have already been described in the TAG framework by [START_REF] Vaillant | A layered grammar model: Using tree-adjoining grammars to build a common syntactic kernel for related dialects[END_REF]. This analysis differs from the one I present here on several aspects: First, it does not rely on a metagrammar, a point which matters since the implementation I present here relies heavily on the concept of metagrammar. Second, it calls upon adjunction to describe the TMA markers, an approach I will question here.

In this section, I will provide arguments in support of integrating TMA markers into the elementary trees as extended projection of a verbal head9 , instead of using adjunction as a way to integrate TMAs (as auxiliary trees) into the structure. Naturally, it is well known that every grammar using substitution can be rewritten using adjunction only. Formal considerations are therefore of little help in this discussion as far as the derived tree is concerned, but if derivation trees as intended to reflect the meaning of the sentence, the way the trees combine is a crucial matter. I base my argumentation on linguistic arguments and propose below a series of tests which show that TMAs behave differently from verbs.

Several tests have been proposed for Romance languages (in [START_REF] Abeillé | Les prédicats complexes dans les langues romanes. Les langues romanes[END_REF] that cannot be used for GC. These are based on clitic-climbing and infinitival constructions which do not exist as such in GC and are therefore of little help. Instead, I will base my arguments on coordination and cleft structures.

Coordination: While coordination can apply to lexical items almost without restriction (see [START_REF] Bernabé | Fondal-natal[END_REF] for a review of coordination in GC), TMA markers cannot be coordinated on the same position: This test highlights the difference between (still) real periphrastic elements and TMAs. For instance, it is usually accepted that the form kay (ka+ay) 'Imperfective+go' (example ( 13)) is one of the TMA (see [START_REF] Vaillant | A layered grammar model: Using tree-adjoining grammars to build a common syntactic kernel for related dialects[END_REF] and [START_REF] Bernabé | Fondal-natal[END_REF](Bernabé, , 1035))). Then, the verb ay 'to go' will be inserted using adjunction, as in [START_REF] Abeillé | Une Grammaire électronique du Français[END_REF] for aspectual verbs and auxiliaries in French. It anchors its own tree and is stored in the Lexicon, contrarily to the TMAs.

The negative marker is also excluded from the clefted position (15), lending weight to the hypothesis that the clefted element is below any verbal functional projections [ V P V (NP)].

( But the behavior of vlé, here, is atypical, since the negative marker follows the verb10 . From these tests, I conclude that TMA markers are functional elements, just like inflectional affixes in French are. The only -but significant-difference is that these markers are not expressed via synthesis but via a particular form of periphrasis (a situation where two or more words express the meaning expected for a single word). Syntactic intermediate projections are available between Tense and Aspect and allow some (rare) adverbs to be inserted between the TMA markers, as in ( 17 vin. come 'Pierre was already coming.' (from [START_REF] Bernabé | Fondal-natal[END_REF](Bernabé, , 1060)) Thus, I treat TMA markers as functional elements that form an extended projection of the lexical head, in conformity with the CETM.

This work is in accordance with many recent works attempting to reconsider the traditional division between morphology and syntax. For instance, [START_REF] Brown | Defining 'periphrasis': key notions[END_REF] try to avoid the binary division of labor between morphology and syntax and claim that periphrasis are both. I will follow here some key assumptions presented in the Paradigm Function Morphology (PFM) framework [START_REF] Stump | Inflectional morphology: a theory of paradigm structure[END_REF]; Bonami and Stump (prep)) which are highly compatible with a lexicalist approach of syntax. This separates clearly the TAG grammar I propose here from the standard GB analyses (in the wake of [START_REF] Chomsky | Lectures on the Theory of Government and Binding[END_REF]) which are non-lexicalist. In particular, I adopt the idea that "some syntactic constructions express the pairing of a lexeme with a morphosyntactic set, and should thus be incorporated into inflectional paradigms as inflectional periphrases" (Bonami and Stump,prep,p.20).

In the cases of the GC's TMAs, it is reasonable to claim that the TMAs and the lexical verb realize different forms of the same lexeme. [START_REF] Bonami | The phrase-structural diversity of periphrasis: a lexicalist account[END_REF] note that most of the morphological descriptions of functional periphrasis are rudimentary in their description of the periphrasis syntactic structures and that their syntactic behavior differs from one language to another (see also [START_REF] Abeillé | Les prédicats complexes dans les langues romanes. Les langues romanes[END_REF]. I suggest that TMAs are in part similar to proclitics, an analysis which has also been proposed in the HPSG framework by [START_REF] Henri | The Morphology of TMA in Mauritian and Creole[END_REF], but still allow syntactic nodes for adjunction. The present study -as well as [START_REF] Henri | The Morphology of TMA in Mauritian and Creole[END_REF] and [START_REF] Pfänder | Describing São Tomense Using a Tree-Adjoining Meta-Grammar[END_REF] for other creole languages -suggest a finer-grained characterization of the TMAs from a typological point of view.

3 TMAs in the TAG Metagrammar

eXtensible Meta-Grammar

XMG 11 is a declarative language for specifying tree-based grammars at a meta-level [START_REF] Crabbé | Metagrammar Redux[END_REF]; [START_REF] Crabbé | XMG : eXtensible MetaGrammar[END_REF]). XMG allows the linguist to capture generalizations on his grammar by defining tree fragments (Classes) that can combine via feature unification (conjunctive / disjunctive combinations of fragments). Once defined, Classes can be reused in distinct contexts, allowing elegant generalizations. A core grammar is described by fragments of elementary trees and these fragments combine to form the expanded grammar which is made of elementary trees 12 . I will use these concepts to describe the TMA part of the core grammar of GC (Crabbé and Duchier, 2004, p.1).

Fragments of Functional Projections

From the preceding section, I take for granted that TMAs are not stored in the Lexicon (they don't anchor any tree properly) but are co-anchors of the elementary tree associated with verbs 13 .

For instance, (19) illustrates the structure of ( 18).

(18) Jan Jean té ant ka imperf manjé eat 'Jean was eating.' 11 I am using XMG-2 version. 12 See [START_REF] Crabbé | Représentation informatique de grammaires d'arbres fortement lexicalisées : le cas de la grammaire d'arbres adjoints[END_REF] for a large metagrammar of French using XMG. 13 For a similar approach on Santomense's Tense and Aspect markers, see [START_REF] Pfänder | Describing São Tomense Using a Tree-Adjoining Meta-Grammar[END_REF] . 19), S is a projection of V, the maximal functional stretching of the verb.

In XMG's framework, the structure ( 19) is broken down into four pieces (i.e. classes) each containing minimal information. These Classes are listed below.

-CanS ub ject: for the External Argument of the verb. It is described in (20-a) .

-Imper f : as a projection of the Imperfective (Aspect) marker. It is described in (20-b). -T ensed: as a projection of Tense.

-Intransitive verb: the minimal projection of V. It is described in (20-d).

(20) a. S V [pro j:T |Imper f |V] NP↓ b.

[pro j:Imper f ]

V [pro j:V] ka c.

V [pro j:T ] V [pro j:Imper f |V] té d.

V [pro j:V] V Thus, ( 19) is built up from the following conjunction of Classes:

CanS ub ject ∧ Intransitive ∧ Imper f ∧ T ensed
As in [START_REF] Pfänder | Describing São Tomense Using a Tree-Adjoining Meta-Grammar[END_REF] for Santomense, the feature proj(ection) is used to rule out invalid combinations in the output elementary tree.14

From the conjunction of classes given above, the result of the metagrammar compilation are elementary trees for intransitive verbs, inflected with Tense and Aspect, as shown in (21)15 .

(21)

The same mechanism is used to describe verb families. For instance, the class of Intransitive verbs inflected for TMAs can be described as a class IntransV gathering the classes of the inflected forms for intransitive verbs16 :

BareV | Imper f V | ProspImper f V | T ensProspImper f V | T ensedV | ProspV | T ensImper f V | T ensProspV
As expected from the morphological approach I defend here, the TMAs do not appear as adjuncts but are co-anchors of the verb. ( 22) shows the derivation tree17 for the sentence Jan té ka manjé 'Jean was eating'. (22) α1-manjé[Ant; Imper f ] α2-Jan

The benefit we have here with regards to adjunction is that the semantic interpretation of the sequence N0 té ka V is directly derivable from the features on the verb (see Table 1). There is no need to postulate ambiguous TMA markers triggering various interpretations or different syntactic positions for the same marker, as is commonly proposed in creole languages descriptions (see [START_REF] Winford | Creole Languages. The Oxford Handbook of Tense and Aspect[END_REF]), i.e. postulating several distinct ka or Ø markers to derive the correct interpretation. On the contrary, it is compatible with the PFM idea that words (or multiword expressions) realize the morphosyntactic property set of the lexeme they are member of. Then, té ka manjé is the realization of the morphosyntactic set {Anterior, imper f ective} for the lexeme manjé.

Extended Projections in the Nominal Domain

The same methodology developed for the TMAs is adopted in the Nominal domain, incorporating Definite, Demonstrative and Plural markers as co-anchors of a lexical Noun, as in ( 23), where the Plural marker sé precedes the N and the Definite marker follows the N, as illustrated in (24). Note that in GC, the Plural requires the Definite, *sé timoun being agrammatical. The genitive preposition a18 'of' in (25), which only has a functional role akin to a Genitive case marker is also treated as a co-anchor.

In (25), the NP [a Lelette] is adjoined to the head noun kaz.

(25) [kaz house [a Genitive

Lelette]] Lelette 'Lelette's house'

The genitival form of the noun (as head of the genitive NP) is an auxiliary tree (for adjunction) containing the co-achor a, as shown in (26). The comparison of GC with a Matiniké Creole (spoken on the island of Martinique, closely related to GC and diachronically linked to GC) lends weight to this approach, since in Matiniké, the same NP as in ( 25) would be kaz Lelette, without preposition. The difference between these creoles relies only in the fact that Matiniké has a covert (silent) preposition for the same structure as (26).

To summarize, the inflected forms of a noun can be represented as a conjunction of classes:

Demonstrative -N | Plural -N | De f inite -N | Bare -N | Genitive -N
Note that while Plural, Demonstrative, Definite and Genitive appear as syntactic projections in GC, the same functional elements are concatenated at the word level in other languages.

Conclusion

I have shown that Tense and Aspect markers of Guadeloupean Creole are functional projections that can be described as co-anchors of a lexical head (extended projections). Their response to the tests of coordination and cleft structures shows that TMAs are not as free as other periphrastic elements (such as the verb ay 'to go'). Their place is somewhere between concatenative morphology at the word level and 'free' periphrastic elements, which I take to be adjoined elements. In the nominal domain, I have suggested that the definite, demonstrative, plural and genitive markers are also functional elements of the same sort. These functional elements are integrated in the elementary trees that form the GC grammar as co-anchors of the lexical item. I have demonstrated how these elements form fragments of (elementary) trees and how they combine to form the expanded grammar. To do this, I have used the XMG formalism for metagrammar. Since the combination of the fragments constitute inflected forms of a lexeme (as the compound tenses are still members of the verbal paradigm, see [START_REF] Ackerman | Paradigms and periphrastic expression: a study in realization-based lexicalism[END_REF]), the building of Elementary trees is as much a morphological as a syntactic operation. It thus casts a new light on Creole languages which are commonly thought to have little (or even no) morphology.

  Sophie.' (and not: * Jean loved Sophie)Tense: The anterior marker of GC is té. When combined with non-stative verbs, té (ant) provides a perfective interpretation: '(At this time,) Jean was loving Sophie' Aspect: GC Aspect markers are ka and ké for Imperfective and Prospective respectively. When ka is combined with stative verbs, the reading must be Iterative, as in:

  goes up, up, up!' (from(Bernabé, 1983, 545))(9) shows that the repetition of the predicate is the norm, since coordination of TMA is blocked, and (10) shows that verb coordination can appear below TMAs. (against him in the past), I won (this time) and I will win!' dancing and singing' Predicate Cleft: Predicate cleft is a frequent construction in GC (while impossible in French). In (11), the lexical verb is clefted without the TMA marker. The same sentence with clefted TMA is ungrammatical, as in (12). climbing to arrive where we are going' (intensive meaning) (12) a. *sé ka monté nou ka monté... [Aspectual marker prohibited] b. *sé té monté nou té monté...[Tense marker prohibited]

  going to come' But ay 'go' differs from the TMAs w.r.t this test, as shown in (14going to come.' (stress on the movement)

  is the conjunction of the classes:Plural | De f inite | Nounwhere Plural can only appear if Definite is present.

Table 1 .

 1 Tense and Aspect markers.

	Value	Form
	Accomplished/Aoristic	dansé
	Unaccomplished / Present	ka dansé
	Frequentative	ka dansé
	Progressive	ka dansé
	Future	ké dansé
	Unaccomplished Future (seldom)	ké ka dansé
	Accomplished past (pluperfect)	té dansé
	Unaccomplished past	té ka dansé
	Irrealis (Past)	té ké dansé
	Irrealis unaccomplished (extremely rare) té ké ka dansé
	Conditional / Optative	té dansé

I would like to thank Simon Petitjean, Yannick Parmentier and Denys Duchier for their precious help. I would also like to express my gratitude to the anonymous reviewers for their precious comments.

2 I will prefer the term adjunction in the remainder of the paper.

For English, see the XTAG project(XTAG Research Group (2001))

In a sense that will be clarified later.

There is another marker kay, which is sometimes analysed as part of this system, however, as I will show later, it can't be included among the TMAs.

This table shows the main uses of the TMA markers to give the reader a quick overview of the TMAs, but, as it will be claimed below, the interpretation of a TMA sequence is highly dependent on the context.

I leave aside the marker laka which can plausibly be analysed as la (locative) and ka.

And more generally with non-stative predicates, such as adjectival predicates.

This can be extended to predicates in general, since there are non-verbal predicates in GC.

It might illustrate a case of lexical integration of pa into the verbal stem.

In the values associated with feature proj, "|" refers to disjunction.

Which is a picture of the elementary tree as displayed by XMG.

The abbreviations are: Imperf(ective), Prosp(ective), Tens(ed).

Where α is the conventional label for an elementary tree.

And its allomorph an when preceding a nasal.