Ibtihel Ben Gharbia

Joëlle Ferzly
email: joelle.ferzly@ifpen.fr

Martin Vohralík
email: martin.vohralik@inria.fr

Soleiman Yousef
email: soleiman.yousef@ifpen.fr

Semismooth and smoothing Newton methods for nonlinear systems with complementarity constraints: adaptivity and inexact resolution *

Keywords: nonlinear complementarity constraints, semismooth smoothing Newton methods, interior-point method, a posteriori error estimate, adaptivity, stopping criteria

We consider nonlinear algebraic systems with complementarity constraints stemming from numerical discretizations of nonlinear complementarity problems. The particularity is that they are non-differentiable, so that classical linearization schemes like the Newton method cannot be applied directly. To approximate the solution of such nonlinear systems, an iterative linearization algorithm like the semismooth Newton-min or an interior-point algorithm can be used. Alternatively, the non-differentiable nonlinearity can be smoothed, which allows a direct application of the Newton method. Corresponding linear systems can be solved only approximately using an iterative linear algebraic solver, leading to inexact approaches. In this work, we design a general framework to systematically steer these different ingredients. We first derive an a posteriori error estimate given by the norm of the considered system's residual. We then, relying on smoothing, design a simple strategy of tightening the smoothing parameter. We finally distinguish the smoothing, linearization, and algebraic error components, which enables us to formulate an adaptive algorithm where the linear and nonlinear solvers are stopped when the corresponding error components do not affect significantly the overall error. Numerical experiments indicate that the proposed algorithm, possibly in combination with the GMRES algebraic solver, ensures important savings in terms of the number of iterations and execution time. It appears rather promising in comparison with the other methods, namely since its performance seems remarkably stable over a range of academic and industrial problems.

Introduction

Consider a system of algebraic equations with complementarity constraints written in the following form: Find a vector X ∈ R n such that

EX = F , (1.1a) K(X) ≥ 0, G(X) ≥ 0, K(X) • G(X) = 0, (1.1b)
where for two integers n > 1 and 0 < m < n, E ∈ R n-m,n is a matrix, K : R n → R m and G : R n → R m are (linear) operators, and F ∈ R n-m is a given vector. The first line (1.1a) typically represents the discretization of a linear partial differential equation. The second line (1.1b) then represents the complementarity constraints. It states that the vectors K(X) and G(X) have nonnegative components and are orthogonal.

Complementarity problems have important applications in many fields: economics, engineering, operations research, nonlinear analysis... In the literature, many theoretical results and numerical methods have been proposed

to solve problem (1.1), see for example the books of Facchinei and Pang [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF], Ito and Kunisch [START_REF] Ito | Lagrange multiplier approach to variational problems and applications[END_REF], Ulbrich [START_REF] Ulbrich | Semismooth Newton methods for variational inequalities and constrained optimization problems in function spaces[END_REF],

Bonnans et al. [START_REF] Bonnans | Numerical optimization[END_REF], and the study of Aganagić [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF].

By means of so-called C-functions (C for complementarity), see [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF], the complementarity constraints (1.1b) can be rewritten as a system of equations C(X) = 0, where C : R n → R n is nonlinear and semismooth. We then obtain the following equivalent formulation of problem (1.1): Find X ∈ R n such that EX = F , C(X) = 0.

(1.2)

A direct application of the standard Newton method to (1.2) is, however, impeded by the fact that C(X) is not

differentiable. An introduction of the Clarke differential [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] allows to give a weaker differentiability meaning and leads to the class of semismooth Newton methods, with reputedly good convergence properties [START_REF] Munson | The semismooth algorithm for large scale complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Gharbia | Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P-matrix[END_REF][START_REF] Gharbia | An algorithmic characterization of P-matricity[END_REF][START_REF] Gharbia | An algorithmic characterization of P-matricity II: adjustments, refinements, and validation[END_REF][START_REF] Dussault | A lower bound on the iterative complexity of the Harker and Pang globalization technique of the Newton-min algorithm for solving the linear complementarity problem[END_REF][START_REF] Dussault | Polyhedral Newton-min algorithms for complementarity problems[END_REF]. These methods are in certain cases equivalent to primal-dual active set strategies, see Hintermüller et al. [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF]. Moreover, in [START_REF] Xiao | A regularized semi-smooth Newton method with projection steps for composite convex programs[END_REF], a regularized semismooth Newton method combined with a hyperplane projection technique was proposed.

Augmented Lagrangian method is one of the commonly used algorithms for constrained optimization, see, e.g., [START_REF] Ito | Augmented lagrangian methods for nonsmooth, convex optimization in hilbert spaces[END_REF] and the references therein. It seeks a solution by replacing the original constrained problem by a series of unconstrained problems and add to the objective function a penalty term, and another term designed to mimic a Lagrange multiplier.

An additional technique, often used in a function space setting, consists in introducing a proper regularization, motivated by the augmented Lagrangian method. It allows to apply an infinite-dimensional semismooth Newton method for the solution of the regularized problem, see, e.g., [START_REF] Ulbrich | Semismooth Newton methods for variational inequalities and constrained optimization problems in function spaces[END_REF]. In the present context, this leads to replacing the complementarity conditions (1.1b) by K(X) + min{0, -K(X) + γG(X)} = 0, for a parameter γ > 0. This method can be combined with a path-following strategy to update the regularization parameter γ, see for instance [START_REF] Stadler | Path-following and augmented Lagrangian methods for contact problems in linear elasticity[END_REF][START_REF] Hintermüller | Path-following methods for a class of constrained minimization problems in function space[END_REF][START_REF] Stadler | Semismooth Newton and augmented Lagrangian methods for a simplified friction problem[END_REF].

Another important class of methods for constrained optimization problems of the form (1.1) is formed by interiorpoint methods. These methods consist in generating a sequence in the feasible region K(X) ≥ 0 and G(X) ≥ 0, under the assumption of knowing a feasible initial point. We refer to the work of Wright [START_REF] Wright | The interior-point revolution in optimization: history, recent developments, and lasting consequences[END_REF], Bellavia et al. [START_REF] Bellavia | An affine scaling trust-region approach to bound-constrained nonlinear systems[END_REF],

and the references therein for a review.

Lastly, an additional notable method is the smoothing Newton method. The main idea of this approach is to approximate the semismooth (non-differentiable) function C from (1.2) by a smooth (differentiable) function that depends on a smoothing parameter. The problem is reformulated as a sequence of regularized smooth equations that can be solved by applying the standard Newton method, and where one drives the smoothing parameter down to zero, cf. [START_REF] Rui | A smoothing inexact Newton method for nonlinear complementarity problems[END_REF][START_REF] Qi | Smoothing functions and smoothing Newton method for complementarity and variational inequality problems[END_REF][START_REF] Qi | A smoothing Newton method for general nonlinear complementarity problems[END_REF] and the references therein.

In this work, we design a general framework to systematically steer the above different ingredients. Our main philosophy is adaptive smoothing (regularization). For µ j > 0, let a smoothed function C µ j (•), satisfy C µ j (X) -C(X) → 0 as µ j → 0, for X ∈ R n . The smoothing parameter µ j is reduced at each smoothing iteration j ≥ 1.

Thus, problem (1.1), or equivalently (1.2), can be reformulated as a system of smooth (differentiable) equations written in the form: Find X j ∈ R n such that

EX j = F , C µ j (X j ) = 0. (1.3)
Hence, Newton-type methods can be applied to solve system (1.3), yielding, at each linearization step k ≥ 1, a linear system

A j,k-1 µ j X j,k = B j,k-1 µ j , (1.4)
where A j,k-1 µ j ∈ R n,n is a matrix and B j,k-1 µ j ∈ R n is a vector.

Solving (1.4) with a direct method may be very expensive. A popular approach is to solve it approximately by applying only a few steps of an iterative algebraic solver. Such inexact approaches can be found in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF][START_REF] Martínez | Inexact Newton methods for solving nonsmooth equations[END_REF] for semismooth Newton methods, in [START_REF] Rui | A smoothing inexact Newton method for nonlinear complementarity problems[END_REF][START_REF] Ge | A smoothing inexact Newton method for variational inequalities with nonlinear constraints[END_REF] for smoothing Newton methods, in [START_REF] Zhao | A Newton-CG augmented Lagrangian method for semidefinite programming[END_REF] for augmented Lagrangian methods, and in [START_REF] Bellavia | Inexact interior-point method[END_REF] for interior-point methods. In the algorithms introduced therein, the iterations of different solvers are stopped according to a fixed maximal number of iterations, the Euclidean norm of the residual vector, or other parameters-dependant stopping criteria. In this work, the a posteriori estimate constitute a distinctive element at the heart of the proposed smoothing method. Importantly, it ensures the desired balance between each source of error at any resolution step, unlike existing approaches based on classical stopping criteria.

At each linear algebraic step i ≥ 1 for (1.4), one in particular obtains X j,k,i ∈ R n such that

A j,k-1 µ j X j,k,i = B j,k-1 µ j -R j,k,i alg ,
where R j,k,i alg ∈ R n is the algebraic residual vector of (1.4).

Our principal aim is to reduce the computational cost of the numerical resolution of (1.1) by employing an adaptive strategy based on a posteriori error estimates. There is a well-developed literature on a posteriori error estimates and mesh adaptivity for partial differential equations, see for instance the books of Ainsworth and Oden [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF], Repin [START_REF] Repin | A posteriori estimates for partial differential equations[END_REF], and Nochetto et al. [START_REF] Nochetto | Theory of adaptive finite element methods: an introduction[END_REF]. For variational inequalities, we can mention the contributions of Repin [START_REF] Repin | Functional a posteriori estimates for elliptic variational inequalities[END_REF],

Ben Belgacem et al. [START_REF] Ben Belgacem | On the unilateral contact between membranes. Part 2: a posteriori analysis and numerical experiments[END_REF], Bürg and Schröder [START_REF] Bürg | A posteriori error control of hp-finite elements for variational inequalities of the first and second kind[END_REF], and Dabaghi et al. [START_REF] Dabaghi | Adaptive inexact semismooth Newton methods for the contact problem between two membranes[END_REF]. Although smoothing Newton approaches have been widely studied, to the best of our knowledge, almost no work has been done to this day on a posteriori error estimates and adaptivity for solvers applied to discrete problems of the form (1.1).

We first derive an upper bound on the norm of the residual of system (1.2), given by

R(X j,k,i ) := F -EX j,k,i -C(X j,k,i ) .
Then, decomposing R(X j,k,i ), we distinguish the different error components. This leads to an a posteriori control of the form

R(X j,k,i ) r ≤ η j,k,i = η j,k,i sm + η j,k,i lin + η j,k,i alg . (1.5)
Here, η j,k,i is a fully computable upper bound that holds true at any smoothing (regularization) step j, linearization step k, and algebraic solver step i, whereas the role of the estimators η j,k,i sm , η j,k,i lin , and η j,k,i alg is to identify the smoothing, linearization, and algebraic components of the error. This error bound allows to define adaptive stopping criteria for the nonlinear and linear algebraic solvers, in the spirit of [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Dabaghi | Adaptive inexact semismooth Newton methods for the contact problem between two membranes[END_REF], and the references therein. These criteria, as well as a simple way to tighten the smoothing parameter µ j , are incorporated in a three-level adaptive algorithm. In contrast to common approaches, where the termination requires reaching a fixed threshold, the particularity of this adaptive algorithm is that the iterations are stopped when the error component of the concerned solver is smaller than the total error, up to a desired fraction. The efficiency of the proposed adaptive algorithm for (inexact) smoothing Newton methods and (inexact) interior-point methods is showcased numerically on practical problems.

It is relevant to mention that this work is extented in [START_REF] Ben Gharbia | Adaptive inexact smoothing Newton method for a nonconforming discretization of a variational inequality[END_REF], where the present approach is applied to a system of PDEs with complementarity constraints in infinite-dimensional space. In particular, taking into account the discretization error allows to adaptively steer the smoothing in system (1.3). Although we do not address mesh adaptivity in our work, we underline that a posteriori estimators are an important tool for adaptive mesh refinement strategies, see, e.g., [START_REF] Di Pietro | An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media[END_REF] and the references therein. Consequently, algorithms based on the previous criteria ensure significant computational gains in terms of total number of iterations and mesh cells.

Our manuscript is organized as follows. In Section 2, we recall a semismooth Newton method based on an equivalent reformulation of the complementarity constraints in the form (1.2). Section 3 is devoted to introduce our adaptive inexact smoothing Newton method based on the reformulation as a system of smooth equations as in (1.3). We establish here the a posteriori error estimates (1.5) and propose an adaptive algorithm with a posteriori stopping criteria. We survey a nonparametric interior-point method in Section 4, and introduce its adaptive version in Section 5. Finally, a detailed numerical study is presented in Sections 6 and 7.

Semismooth Newton method

The purpose of this section is to briefly recall the semismooth Newton method to approximate the solution of the nonlinear system of equations (1.1), see, e.g., [START_REF] Munson | The semismooth algorithm for large scale complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Dabaghi | Adaptive inexact semismooth Newton methods for the contact problem between two membranes[END_REF]. The complementarity constraints represented by (1.1b) as algebraic inequalities are here rewritten as non-differentiable algebraic equalities, using a complementarity function

(C-function). A function C : R m × R m - → R m , m ≥ 1, is called a C-function if C(x, y) = 0 ⇐⇒ x ≥ 0, y ≥ 0, x • y = 0 ∀(x, y) ∈ R m × R m .
A variety of C-functions can be found in the literature, see, e.g., [START_REF] Sun | On NCP-functions[END_REF][START_REF] Galántai | Properties and construction of NCP functions[END_REF]. We give as examples the minimum (min) function and the Fischer-Burmeister (F-B) function: for l = 1, . . . , m, Cmin (x, y)

l := (min(x, y)) l = (x l + y l )/2 -|x l -y l |/2, (2.1)
CFB (x, y)

l := x 2 l + y 2 l -(x l + y l ). (2.2)
In general, the C-functions are not Fréchet differentiable. The min and the Fischer-Burmeister functions are, for example, differentiable everywhere except in x = y and (0, 0), respectively. Let us introduce a function

C : R n → R m defined as C(X) := C (K(X), G(X))
, where C : R m × R m → R m is any C-function. By using this reformulation in (1.1b), it is obvious that problem (1.1) can be equivalently rewritten as: Find a vector X ∈ R n , such that

EX =F , (2.3a) C(X) = 0. (2.3b)
Next, we detail the semismooth Newton linearization. Let an initial vector X 0 ∈ R n be given. At the step

k ≥ 1, one looks for X k ∈ R n such that A k-1 X k = B k-1 , (2.4)
where the square matrix A k-1 ∈ R n,n and the right-hand side vector B k-1 ∈ R n are given by

A k-1 := E J C (X k-1 ) , B k-1 := F J C (X k-1 )X k-1 -C(X k-1 )
.

(2.5)

Note that the Jacobian corresponding to (2.3a) is constant and equal to E since it is linear. The semismooth nonlinearity occurs in the second line (2.3b): the notation J C in (2.5) stands for the Jacobian matrix in the sense of Clarke of the function C, cf. [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF]. To give an example, consider the semismooth min function (2.1) and define the matrices K and G ∈ R m,n respectively by K := [∇K(X)] and G := [∇G(X)] . Then the l th row of the Jacobian matrix in the sense of Clarke J C is either given by the

l th row of K, if (K(X k-1 )) l ≤ (G(X k-1 )) l , or by the l th row of G, if (G(X k-1 )) l < (K(X k-1 )) l .
We will need below the total residual vector of problem (2.3), defined by

R(V ) := F -EV -C(V ) , V ∈ R n . (2.6)
In this context, the relative norm of a vector V ∈ R n is given by

||V || r := ||V || / R(X 0 ) , where • is the L 2 -norm.

Adaptive inexact smoothing Newton method

In this section we introduce our adaptive inexact smoothing Newton method. Based on a posteriori error estimators, adaptive stopping criteria are formulated to conceive an adaptive iterative algorithm.

Smoothing of the C-functions

The key of our developments is to smooth the non-differentiable equation formulation (2.3b) of the complementarity constraints (1.1b) with the help of a smooth (i.e. continuously differentiable) function. This smoothing allows us to approximately transform the nonsmooth nonlinear system (2.3) to a smooth system of nonlinear equations to be solved by using the standard Newton method.

Let µ > 0 be a (small) smoothing parameter. We construct an approximation function Cµ :

R m × R m → R m of a C-function C such that Cµ (•, •) is of class C 1 on R m × R m and satisfies C(x, y) -Cµ (x, y) → 0 as µ → 0 for all (x, y) ∈ R m × R m .
For example, for l = 1, . . . , m, a possible smoothing of the min and the Fischer-Burmeister functions (2.1) and (2.2) can be

Cminµ (x, y) l = x l + y l 2 - |x -y| µ l 2 , with (|z| µ ) l = z 2 l + µ 2 , (3.1)
CFBµ (x, y) We define the function C µ : R n → R m as C µ (X) := Cµ (K(X), G(X)) , where Cµ : R m × R m → R m is any smoothed C-function of at least class C 1 . This allows to approximate problem (1.1) or (1.2) by a system of smooth equations: Find a vector X ∈ R n , such that

l = µ 2 + x 2 l + y 2 l -(x l + y l ), (3.2
EX = F , C µ (X) = 0. (3.3)
Thus, Newton-type methods can be applied to solve the system of nonlinear algebraic equations (3.3).

Fixing µ 1 > 0, we now describe an iterative method for solving problem (2.3). At the beginning of each smoothing iteration (outer iteration) denoted hereafter by j ≥ 1, an initial guess X j ∈ R n is given, and a smoothing parameter µ j is determined; µ j will be driven down to zero. Then some iterative nonlinear solver like the Newton method is employed to solve the smoothed problem written in the form: Find X j ∈ R n such that

EX j = F , C µ j (X j ) = 0.
(3.4)

Newton linearization of the nonlinear algebraic system

In what follows, we detail the Newton method employed to solve problem (3.4) at a fixed outer smoothing step j ≥ 1. Given an initial vector X j,0 (typically X j,0 = X j-1 ), Newton's algorithm generates a sequence (X j,k ) k≥1 with X j,k ∈ R n given by the following system of linear algebraic equations

A j,k-1 µ j X j,k = B j,k-1 µ j , (3.5)
where the Jacobian matrix A j,k-1 µ j ∈ R n,n and the right-hand side vector B j,k-1 µ j ∈ R n are defined by

A j,k-1 µ j := E J C µ j (X j,k-1 ) , B j,k-1 µ j := F J C µ j (X j,k-1 )X j,k-1 -C µ j (X j,k-1 ) , (3.6)
with J C µ j (X j,k-1 ) the Jacobian matrix of the smooth function C µ j at X j,k-1 .

Inexact solution of the linear algebraic system

The linearized system (3.5) may not be solved exactly, since the use of a direct method may be expensive. For this reason, we consider in this work also an inexact resolution. For a fixed smoothing step j ≥ 1, a fixed Newton step k ≥ 1, and an initial guess X j,k,0 (typically X j,k,0 = X j,k-1 ), only a few steps of an iterative linear algebraic solver can be applied to find an approximate solution to (3.5), yielding, on step i ≥ 1, an approximation X j,k,i to X j,k .

This satisfies (3.5) up to the residual vector given by

B j,k-1 µ j -A j,k-1 µ j X j,k,i . (3.7)
Define now the linearization function C j,k-1 µ j

: R n → R m of C µ j at smoothing step j and Newton step k as

C j,k-1 µ j (V ) := C µ j (X j,k-1 ) + J C µ j (X j,k-1 )(V -X j,k-1 ) ∀ V ∈ R n . (3.8)
This allows us to write the algebraic residual vector for V ∈ R n as

R AISN alg (V ) := B j,k-1 µ j -A j,k-1 µ j V = F -EV -C j,k-1 µ j (V )
.

(3.9)

An upper bound for the norm of the residual

We consider the total residual vector of problem (2.3) given in (2.6). By adding and subtracting C µ j (X j,k,i ) and its linearization C j,k-1 µ j (X j,k,i ) given by (3.8), the total residual vector can be decomposed as follows:

R(X j,k,i ) = F -EX j,k,i -C(X j,k,i ) ± C µ j (X j,k,i ) ± C j,k-1 µ j (X j,k,i ) = 0 C µ j (X j,k,i ) -C(X j,k,i ) smoothing + 0 C j,k-1 µ j (X j,k,i ) -C µ j (X j,k,i ) linearization + F -EX j,k,i -C j,k-1 µ j (X j,k,i ) algebraic .
It is reasonable to get these three terms. Indeed, the first one reflects the error due to the approximation of the semismooth function C by the smoothed function C µ j . The second term is related to the linearization of the nonlinear smooth problem (3.4). Taking into account that the resolution of the smooth linearized problem (3.5) is possibly done "inexactly", the remaining term represents the error of the inexact algebraic resolution. By the triangle inequality, the relative norm of R(X j,k,i ) is thus bounded by the smoothing, linearization, and algebraic estimators respectively defined as

η j,k,i sm,AISN := C µ j (X j,k,i ) -C(X j,k,i ) r , (3.10a)
η j,k,i lin,AISN := C j,k-1 µ j (X j,k,i ) -C µ j (X j,k,i ) r , (3.10b)
η j,k,i alg,AISN := F -EX j,k,i 2 r + C j,k-1 µ j (X j,k,i ) 2 r 1 2
.

(3.10c)

Note that η j,k,i alg,AISN is exactly equal to the relative norm of R AISN alg (X j,k,i ) given by (3.9). From these developments we conclude:

Theorem 3.1 Let X j,k,i ∈ R n arise from an inexact solve of (3.5). We have

R(X j,k,i ) r ≤ η j,k,i AISN := η j,k,i sm,AISN + η j,k,i lin,AISN + η j,k,i alg,AISN .

Adaptive inexact smoothing Newton algorithm

Theorem 3.1 motivates the following. Let two real parameters α lin and α alg be given in ]0, 1], representing the desired relative size of the algebraic and linearization errors, and let ε > 0 be a given desired tolerance for the total error. The stopping criteria for the linearization, algebraic, and smoothing steps, with the bars denoting the stopping indices, are respectively set as

η j,k,i alg,AISN < α alg η j,k,i lin,AISN , (3.11a)
η j,k,i lin,AISN < α lin η j,k,i sm,AISN , (3.11b) R(X j,k,i ) r < ε. (3.11c)
The first criterion (3.11a) for the algebraic iterative solver expresses that there is no need to continue with the algebraic steps when the linearization error becomes dominant. Similarly, the second one (3.11b) aims at stopping the linearization iterations when the linearization error does not substantially contribute to the smoothing error.

Finally, the termination criterion for the smoothing steps (3.11c) is of the standard type, that is when we stop the entire procedure, when the relative norm of the total residual vector lies below the desired tolerance ε.

The entire method is described by the following adaptive algorithm, which drives the smoothing parameter µ j to zero as µ j := αµ j-1 at each smoothing iteration. Other common empirical ways to progressively reduce µ j can be found, e.g., in [START_REF] Vu | A new approach for solving nonlinear algebraic systems with complementarity conditions. Application to compositional multiphase equilibrium problems[END_REF]. The adaptive inexact smoothing Newton algorithm is the following:

Algorithm 1 Adaptive inexact smoothing Newton algorithm

Initialization

Choose a tolerance ε > 0 and parameters α ∈ ]0, 1[ and α lin , α alg ∈ ]0, 1]. Fix µ 1 > 0 and an initial approximation X 0 ∈ R n . Set j := 1.

Smoothing loop

2.1 Set X j,0 := X 0 as an initial guess for the nonlinear solver. Set k := 1.

Newton linearization loop

2.2.1 From X j,k-1 define A j,k-1 µ j ∈ R n,n and B j,k-1 µ j ∈ R n by (3.6).

2.2.2

Consider the problem of finding a solution X j,k to

A j,k-1 µ j X j,k = B j,k-1 µ j . ( 3.12)
2.2.3 Set X j,k,0 := X j,k-1 as initial guess for the iterative algebraic solver. Set i := 1.

Algebraic solver loop i)

Starting from X j,k-1 , perform a step of the iterative algebraic solver for the solution of (3.12), yielding, on step i an approximation X j,k,i to X j,k satisfying

A j,k-1 µ j X j,k,i = B j,k-1 µ j -R AISN alg (X j,k,i ).
ii) Compute the estimators given in (3.10).

iii) If η j,k,i alg,AISN < α alg η j,k,i lin,AISN , set i := i and stop. If not, set i := i + 1 and go to i). 2.2.5 If η j,k,i lin,AISN < α lin η j,k,i sm,AISN , set k := k and stop. If not, set k := k + 1 and go to 2.2.1. 2.3 If R(X j,k,i ) r < ε, set j := j and stop.

If not, set j := j + 1, X j,0 := X j-1,k,i , and µ j := αµ j-1 . Then set k := 1 and go to 2.2.1.

Nonparametric interior-point method

Now we employ a nonparametric interior-point method to problem (1.1). More precisely, we consider the method introduced in [START_REF] Vu | Numerical resolution of algebraic systems with complementarity conditions[END_REF] where a systematic strategy is used to steer the sequence of smoothing parameters towards zero.

We introduce a vector µ = µ1 ∈ R m , where µ > 0 is the smoothing parameter and 1 ∈ R m is the vector with all components equal to 1. The original nonsmooth problem (1.1) is replaced by a smoothed problem written in the form: Find X ∈ R n such that

EX = F , (4.1a) K(X) ≥ 0, G(X) ≥ 0, K(X)G(X) = µ, (4.1b)
where

[(K(X)G(X)] m = [K(X)] m [G(X)] m .
In order to properly adjust the sequence of smoothing parameters, the smoothing parameter µ is treated as an unknown, by introducing the following new equation into system (4.1)

θµ + µ 2 = 0, (4.2)
where θ is a small positive real parameter, chosen once and for all. This equation prevents µ from rushing to zero in just one iteration, and ensures quadratic convergence, see [START_REF] Vu | Numerical resolution of algebraic systems with complementarity conditions[END_REF]. The unknown of system (4.1) is now the enlarged vector X = (X, µ) T ∈ R n+1 . We are thus brought back to applying the standard Newton method to a smooth problem.

Let X 0 ∈ R n such that K(X 0 ) ≥ 0 and G(X 0 ) ≥ 0 be given. To update the iterate X k-1 , we compute a search direction denoted by

d k = [d k X , d k µ ] ∈ R n+1 , where d k X ∈ R n and d k µ ∈ R.
Then, to preserve positivity of K(X k ) and G(X k ) at each step of the nonlinear solver, a truncation of the Newton direction d k is performed so that the corresponding update satisfies

K(X k-1 + κ k d k X ) ≥ 0 and G(X k-1 + κ k d k X ) ≥ 0 for some κ k ∈ ]0, 1], as
close to 1 as possible. After this, we can set

X k := X k-1 + κ k d k .
Recall that our goal is to make µ equal to 0 in the limit while ensuring the positivity of the updated iterate. Another choice for the additional equation (4.2) added to system (4.1) was developed and introduced in a recent work, see [START_REF] Vu | A new approach for solving nonlinear algebraic systems with complementarity conditions. Application to compositional multiphase equilibrium problems[END_REF]Section 3]. The proposed equation does not require to truncate the Newton direction, and couples µ and X in a tighter way.

We rewrite system (4.1) as an enlarged system of n + 1 equations

EX = F , K(X)G(X) -µ = 0, θµ + µ 2 = 0. (4.3)

Adaptive inexact interior-point method

We present in this section our adaptive inexact version of the nonparametric interior point method of Section 4. In contrast to Section 4, we consider, however, µ > 0 as a parameter, and not as an unknown. At each smoothing step j ≥ 1, we may solve the system of smoothing equations written as: Find X j ∈ R n such that K(X j ) ≥ 0, G(X j ) ≥ 0, and

EX j = F , (5.1a)
H µ j (X j ) := K(X j )G(X j ) -µ j = 0. (5.1b)
The values of µ j are gradually decreased at each smoothing iteration, creating a sequence of suitable µ j converging to zero.

Newton linearization of the nonlinear algebraic system

Let X 0 ∈ R n such that K(X 0 ) ≥ 0 and G(X 0 ) ≥ 0 be given. At each smoothing iteration j ≥ 1 and each linearization step k ≥ 1, starting with an initial approximation X j,0 such that K(X j,0 ) ≥ 0 and G(X j,0 ) ≥ 0 (typically X j,0 = X j-1 ), we try to approach the solution of problem (5.1) by finding X j,k ∈ R n such that

A j,k-1 µ j X j,k = B j,k-1 µ j , ( 5.2)
where the Jacobian matrix A j,k-1 µ j ∈ R n,n and the right-hand side vector B j,k-1 µ j ∈ R n are defined by

A j,k-1 µ j := E J H µ j (X j,k-1 ) , B j,k-1 µ j := F J H µ j (X j,k-1 )X j,k-1 -H µ j (X j,k-1 ) , (5.3)
with J H µ j the Jacobian matrix of H µ j . To ensure the positivity of the complementarity constraints, we then define the direction d j,k := X j,k -X j,k-1 ∈ R n and find κ j,k ∈ ]0, 1] such that

K(X j,k-1 + κ j,k d j,k ) ≥ 0 and G(X j,k-1 + κ j,k d j,k ) ≥ 0.

Inexact solution of the linear algebraic system

For a fixed smoothing iteration j ≥ 1, a fixed Newton step k ≥ 1, and an initial guess X j,k,0 (typically X j,k,0 = X j,k-1 ), an iterative algebraic solver can be applied to approach the solution of (5.2), yielding, on step i ≥ 1, an approximation X j,k,i to X j,k . This satisfies (5.2) up to a residual vector defined by

B j,k-1 µ j -A j,k-1 µ j X j,k,i . (5.4)
Introduce the linearization H j,k-1

µ j : R n → R m of H µ j (•) such that for V ∈ R n , H j,k-1 µ j (V ) := H µ j (X j,k-1 ) + J H µ j (X j,k-1 )(V -X j,k-1
).

(5.5) Using (5.5), the algebraic residual vector can be written as follows

R AIIP alg (V ) := B j,k-1 µ j -A j,k-1 µ j V = F -EV -H j,k-1 µ j (V ) , V ∈ R n . (5.6)
We now define the function

H : R n → R m by H(V ) := K(V )G(V ), V ∈ R n . (5.7)
and the total residual vector associated to the adaptive inexact interior-point method by

R AIIP (V ) := F -EV -H(V ) , V ∈ R n . (5.8)
Here again, the relative norm of a given vector V ∈ R n is given by ||V || r := ||V || / R AIIP (X 0 ) .

An upper bound for the norm of the residual

In the same spirit as in Section 3.4, we decompose at each smoothing step j ≥ 1, each linearization step k ≥ 1, and each algebraic step i ≥ 1 the total residual vector given by (5.8)

R AIIP (X j,k,i ) = 0 H µ j (X j,k,i ) -H(X j,k,i ) smoothing + 0 H j,k-1 µ j (X j,k,i ) -H µ j (X j,k,i ) linearization + F -EX j,k,i -H j,k-1 µ j (X j,k,i ) algebraic .
We then define the smoothing, linearization, and algebraic estimators by

η j,k,i sm,AIIP := H µ j (X j,k,i ) -H(X j,k,i ) r = µ j r , (5.9a) η j,k,i lin,AIIP := H j,k-1 µ j (X j,k,i ) -H µ j (X j,k,i ) r , (5.9b)
η j,k,i alg,AIIP := F -EX j,k,i 2 r + H j,k-1 µ j (X j,k,i ) 2 r 1 2 .
(5.9c)

Then we have an upper bound for the norm R AIIP (X j,k,i ) r :

Theorem 5.1 Let X j,k,i ∈ R n be the approximation of X given by an iterative algebraic solver. Then we have R AIIP (X j,k,i ) r ≤ η j,k,i AIIP := η j,k,i sm,AIIP + η j,k,i lin,AIIP + η j,k,i alg,AIIP .

Adaptive inexact interior-point algorithm

Our proposed adaptive inexact interior-point algorithm implements adaptive stopping criteria formulated using the error component estimators given by (5.9) is as follows:

Algorithm 2 Adaptive inexact interior-point algorithm

Initialization

Choose a tolerance ε > 0 and parameters α ∈ ]0, 1[ and α lin , α alg ∈ ]0, 1]. Fix µ 1 > 0 and an initial vector X 0 ∈ R n such that K(X 0 ) ≥ 0 and G(X 0 ) ≥ 0. Set j := 1.

Smoothing loop

2.1 Set X j,0 := X 0 as an initial guess for the linearization loop and k := 1.

2.2 Interior-point linearization loop

2.2.1 From X j,k-1 define A j,k-1 µ j ∈ R n,n and B j,k-1 µ j ∈ R n by (5.3).

2.2.2

Consider the problem of finding X j,k ∈ R n such that

A j,k-1 µ j X j,k = B j,k-1 µ j .
(5.10) 2.2.3 Set X j,k,0 := X j,k-1 as initial guess for the iterative algebraic solver. Set i := 1.

Algebraic solver loop i)

Starting from X j,k-1 perform a step of the iterative algebraic solver for (5.10), yielding, at step i ≥ 1, a vector X j,k,i ∈ R n such that

A j,k-1 µ j X j,k = B j,k-1 µ j -R AIIP alg (X j,k,i ).
ii) Set d j,k,i := X j,k -X j,k-1 and compute κ j,k,i ∈ ]0, 1] such that K(X j,k-1 + κ j,k,i d j,k,i ) ≥ 0 and G(X j,k-1 + κ j,k,i d j,k,i ) ≥ 0.

Then set X j,k,i := X j,k-1 + κ j,k,i d j,k,i . iii) Compute the estimators given by (5.9). iv) If η j,k,i alg,AIIP < α alg η j,k,i lin,AIIP , set i := i and stop. If not, set i := i + 1 and go to i). 2.2.5 If η j,k,i lin,AIIP < α lin η j,k,i sm,AIIP , set k := k and stop. If not, set k := k + 1 and go to 2.2.1.

2.3

If R AIIP (X j,k,i ) r < ε, set j := j and stop. If not, set j := j + 1, X j,0 := X j-1,k,i , and µ j := αµ j-1 . Then set k := 1 and go to 2.2.1.

Numerical experiments: Problem of contact between two membranes

This section reports some numerical illustrations obtained using the algorithms previously presented. We consider here the model problem of contact between two membranes.

Problem statement

Let Ω = (a, b) be a one-dimensional domain. The problem reads: Find u 1 , u 2 , and λ such that

           -µ 1 ∆u 1 -λ = f 1 in Ω, -µ 2 ∆u 2 + λ = f 2 in Ω, (u 1 -u 2 )λ = 0, u 1 -u 2 ≥ 0, λ ≥ 0 in Ω, u 1 = g on ∂Ω, u 2 = 0 on ∂Ω, (6.1)
where u 1 and u 2 represent the vertical displacements of the two membranes and λ is a Lagrange multiplier characterizing the action of the second membrane on the first one, -λ being the reaction. The constant parameters µ 1 , µ 2 > 0 correspond to the tension of each membrane, whereas f 1 , f 2 ∈ L 2 (Ω) are given external forces. The boundary condition prescribed by a constant g > 0 ensures that, on the boundary ∂Ω, the first membrane is above the second one. The third line of (6.1) represents the linear complementarity conditions which serve to distinguish two different physical situations: either the membranes are separated (u 1 > u 2 and λ = 0), or they are in contact (u 1 = u 2 and λ > 0). We discretize this problem by the finite volume method. The corresponding discretization can be written under the form of problem (1.1).

Test problem setting

Following [START_REF] Ben Belgacem | On the unilateral contact between membranes. Part 2: a posteriori analysis and numerical experiments[END_REF], we set Ω = (-1, 1) and consider the following analytical solution for x ∈ Ω

u 1 (x) := g(2x 2 -1), u 2 (x) := 2g(1 -x 2 )(2x 2 -1) if x < -1 √ 2 or x > 1 √ 2 , g(2x 2 -1) otherwise, λ(x) := 0 if x < -1 √ 2 or x > 1 √ 2 , 2g otherwise.
This triple is the solution of (6.1) for the data f 1 and f 2 given by

f 1 (x) := -4g if x < -1 √ 2 or x > 1 √ 2 , -6g otherwise,
and

f 2 (x) := -12g(1 -4x 2 ) if x < -1 √ 2 or x > 1 √ 2 , -2g
otherwise.

Throughout the computational experiments, the parameters µ 1 and µ 2 are set to 1 and the boundary condition g for the first membrane is taken equal to 0.1. Let N be the number of mesh elements. The initial guess X 0 ∈ R 3N has its first N components equal to g and its other components equal to zero for the semismooth and smoothing Newton methods. For the nonparametric interior-point method (resp. the adaptive interior-point method), the initialization is given by

X 0 = [0.1 0 0.5 0.05] T ∈ R 3N +1 (resp. X 0 = [0.1 0 0.5] T ∈ R 3N
). All the simulations are performed in MATLAB. We consider N = 25000 elements, leading to the matrix A of size n = 75000.

Semismooth Newton method

We start by presenting the numerical results of the semismooth Newton method described in Section 2, using the F-B function (2.2). The stopping criterion is on the total residual vector (2.6)

R(X k ) r < 10 -8 . (6.2)
To achieve this stopping criterion, 527 semismooth Newton-F-B iterations (CPU time: 68.9s) and 2232 Newtonmin iterations (CPU time: 338.9s) are needed. Figure 2 represents the evolution of R(X k ) r as a function of the semismooth Newton-F-B iterations. We can see that it decreases slowly during iterations, then the convergence gets extremely fast at the end. Newton iteration Relative norm of the total residual vector R(X k ) r 2), stopping criterion (6.

2)] Relative norm of the total residual vector (2.6) as a function of semismooth Newton iterations.

Adaptive smoothing Newton method

We now test the adaptive smoothing Newton method, denoted by ASN, with the smoothed F-B function (3.2). This consists in employing the method presented in Section 3, summarized in Algorithm 1, but with an exact resolution of the nonlinear system (3.5). The linearization and smoothing estimators are respectively defined by

η j,k lin,ASN := C µ j (X j,k ) r , (6.3a)
η j,k sm,ASN := C µ j (X j,k ) -C(X j,k ) r , (6.3b)
and the total estimator by η j,k ASN := η j,k sm,ASN + η j,k lin,ASN .

First, we analyze the performance of the adaptive stopping criterion based on the estimators for stopping the linearization steps. We compare it with the classical approach in where the linearization is continued until the relative norm of the linearization estimator becomes smaller than a threshold taken as 10 -8 , i.e., Classical stopping criterion: η j,k lin,ASN < 10 -8 , (

Adaptive stopping criterion: η j,k lin,ASN < α lin η j,k sm,ASN . (6.5)

We set µ 1 = 1, ε = 10 -8 , α lin = 1, and α = 0.1 in Algorithm 1. Figure 3 depicts the evolution of the estimators and the relative norm of the total residual vector R(X j,k ) given in (2.6) as a function of the smoothing Newton-F-B iterations, at a specific smoothing iteration j = 1 (µ 1 = 1), left, and j = 3 (µ 3 = 10 -2 ), right. We can observe from Figure 3, left, that, as expected, the smoothing estimator and R(X j,k ) r stagnates after few steps, since here the smoothing parameter µ 1 is equal to 1, whereas the linearization estimator steadily decreases. If we consider the classical stopping criterion (6.4), the linearization will only be stopped at step k = 8. On the other hand, with our adaptive stopping criterion (6.5), only one iteration is necessary. Clearly after a few linearization steps, the linearization estimator no longer affects significantly the smoothing estimator, and we can economize many useless iterations.

Next, we provide in Table 1 the results obtained using the adaptive stopping criterion (6.5) to stop the nonlinear solver. We terminate the smoothing iterations using the relative norm of the total residual vector (2.6) R(X j,k ) r < 10 -8 . (6.6)

We present the cumulated number of Newton iterations Niter, the estimators (6.3), and the relative norm of the total residual vector (2.6) at each smoothing step j. In terms of numbers, 10 smoothing iterations and 36 cumulated Newton iterations (CPU time: 6.9s) are needed to achieve the stopping criterion (6.6). From Table 1, one can see that for each value of µ j , the Newton iterations are stopped according to (6.5). R(X j,k ) r decreases until lying below 10 -8 . Figure 4 displays the curve of the estimators as a function of cumulated Newton iterations and smoothing iterations, as well as the relative norm of the total residual vector as a function of smoothing iterations. The improvement of the performance with respect to the semismooth Newton-F-B method of Section 6.3 is spectacular. 2), classical and adaptive stopping criteria (6.4) and (6.5)] Relative norm of the total residual vector (2.6) and estimators (6.3) as a function of Newton iterations k, at a specific smoothing iteration j = 1 (µ 1 = 1), left, and at j = 3 (µ 3 = 10 -2 ), right.

Smoothing iteration

Relative norm of the total residual vector R(X j,k ) r 2), adaptive stopping criterion (6.5)] Estimators (6.3) as a function of cumulated Newton iterations (left). Estimators (6.3) (middle) and relative norm of the total residual vector (2.6) (right) as a function of smoothing iterations j at convergence of the linearization solver. With the intention to compare the proposed method to existing methods, we complete the semismooth Newton method by a path-following strategy to solve problem (1.1), following [START_REF] Zhang | Path-following and semismooth Newton methods for the variational inequality arising from two membranes problem[END_REF]. For the sake of brevity, we shall not detail this here. The following test compares the semismooth Newton method (SSN) and the semismooth Newton method with path-following (SSN-pf) in which the linearization is stopped when the criterion (6.2) is satisfied, to the adaptive smoothing Newton method, using the smoothed min and F-B functions (3.1) and (3.2) and the stopping criteria (6.5) and (6.2) respectively for the linearization and smoothing iterations. We compare the number of cumulated linearization iterations and the global CPU time of the simulation for the different strategies. The results are displayed in Figure 5. They confirm the expected reduction of the computational cost of the numerical resolution with our adaptive approaches. Actually, we notice that the semismooth Newton method with pathfollowing (red curve) and the adaptive smoothing Newton method (purple and dark blue curves) require significantly fewer cumulated Newton iterations and time to converge, in comparison with the semismooth Newton method (green and orange curves). Therefore, employing the path-following strategy or the adaptive strategy based on a posteriori error estimates enables to save many unnecessary additional iterations, and yield much better results than the pure semismooth Newton method. We note that, using the adaptive smoothing Newton method, one obtains similar computational results using both the smoothed F-B or the smoothed min function.

µ j Niter η j,k lin,ASN η j,k sm,ASN R(X j,k ) r 1e+00 1

Adaptive inexact smoothing Newton method

We focus in this section on the adaptive inexact Newton method introduced in Section 3 and investigate the performance of Algorithm 1 using the smoothed F-B function (3.2) together with the restarted GMRES method.

Typically, we use a fixed restart parameter equal to 300. The behavior of the adaptive smoothing solvers can be improved dramatically by using good preconditioners. Here, we merely use an ILU preconditioner to speed-up the GMRES solver. For other possibilities for preconditioners, we refer to, e.g., [START_REF] Lacroix | Iterative solution methods for modeling multiphase flow in porous media fully implicitly[END_REF] and the references therein. To point out the efficiency of the adaptivity, we test two approaches. First, we stop the algebraic iterations using the standard GMRES stopping criterion on the relative residual given by

R rel := M 2 \(M 1 \(B j,k-1 µ j -A j,k-1 µ j X j,k,i )) M 2 \(M 1 \(B j,k-1 µ j -A j,k-1 µ j X j,k-1 )) ≤ 10 -10 , (6.7)
where M 1 and M 2 are the preconditioner matrices. Second, we incorporate the adaptive stopping criteria (3.11a) for the algebraic solver in Algorithm 1. We set the parameters µ 1 = 1, ε = 10 -5 , α alg = 10 -3 , α lin = 1, and α = 0.1. GMRES iterations are needed to achieve the standard stopping criterion (6.7), whereas in the adaptive resolution case, only 10 GMRES iterations are required to satisfy the adaptive stopping criterion (3.11a). In this case, we can avoid many unnecessary iterations. One can also see from the right part of Figure 6, for j = 3 and k = 1, that the overall gain in terms of algebraic iterations obtained using our stopping criteria is quite significant.

Smoothing iteration

Estimators

η j,k,i sm,AISN η j,k,i lin,AISN η j,k,i alg,AISN 10

Cumulated GMRES iteration

Estimators part of Figure 7, we can observe the behavior of the estimators at the end of the algebraic iterations, during the linearization iterations. We present 8 curves, each one corresponding to a specific value of µ j . We can see that at each smoothing iteration j, the smoothing estimator η j,k,i sm,AISN stagnates after about two iterations. The linearization estimator η j,k,i lin,AISN decreases until becoming smaller than the smoothing estimator, satisfying the stopping criterion (3.11b). Finally, the detected behavior in terms of all smoothing iterations j, linearization iterations k, and algebraic solver iterations i is presented in Figure 7, right, for j ≤ 2. The overall results are collected in Table 2. We present in particular the number of linearization and cumulated algebraic iterations per smoothing step j, Niter and Giter respectively, as well as the estimators (3.10) and the relative norm of the total residual vector (2.6) at the end of each smoothing step j. Using the adaptive stopping criteria (3.11), 8 smoothing iterations, 39 cumulated Newton iterations, and 5999 cumulated GMRES iterations are needed to ensure convergence. Figure 8 illustrates the performance of the adaptive inexact smoothing Newton method. It represents the ratio between: 1) the number of algebraic iterations (left) and the CPU time (right) using the classical GMRES stopping criterion (6.7) and 2)

η j,k,i sm,AISN η j,k,i lin,AISN η j,k,i alg,AISN j = 1, k = 1 j = 2, k = 1 j = 2, k = 2 j = 2, k = 3 j = 2, k = 4
the number of algebraic iterations and the CPU time using the adaptive stopping criterion (3.11a) for GMRES, as a function of the number of elements. For larger systems, 20-times fewer iterations and 18-times faster execution time are achieved.

µ j Niter Giter η j,k,i lin,AISN η j,k,i sm,AISN η j,k,i alg,AISN R(X j,k,i ) r 1e+00 1

Nonparametric interior-point method

We consider here the nonparametric interior-point approach of Section 4, where the dimension of the corresponding problem is n = 3N + 1. The value of the constant θ in the additional equation (4.2) is 10 -1 . The stopping criterion is on the relative norm of the linearization residual vector

R IP (X k ) r := R IP (X k ) / R IP (X 0 ) < 10 -8 , (6.8)
with

R IP (X k ) :=   F -EX k µ -K(X k )G(X k ) -θµ k -(µ k ) 2   .
Using this method, 19 Newton iterations (CPU time: 6.7s) are needed to reach the end of the simulation. Figure 9 shows that R IP (X k ) r decreases during the Newton interior-point iterations until satisfying the stopping criterion (6.8).

Adaptive interior-point method

Next, we consider the adaptive interior-point method, which is the method presented in Section 5, Algorithm 2 without applying an algebraic iterative solver to approximate the solution of the linear system (5.2). In this case, we can define the linearization and smoothing estimators respectively by η j,k lin,AIP := H µ j (X j,k ) r , (6.9a) Newton iteration Relative norm of the residual vector where H µ j (•) is defined in (5.1b), and the total estimator by η j,k AIP := η j,k sm,AIP + η j,k lin,AIP . Recall from (5.8) the definition of the total residual vector for V ∈ R n as

η j,k sm,AIP := µ j r , (6.9b)
R IP (X k ) r
R AIP (V ) := F -EV -H(V ) , (6.10)
where H(•) is defined in (5.7). The adaptive stopping criterion

η j,k lin,AIP < α lin η j,k sm,AIP (6.11)
is used to stop the nonlinear solver and a criterion on the relative norm of the total residual vector is applied to stop the smoothing iterations R AIP (X j,k ) r < 10 -8 . (6.12)

The initial smoothing vector is µ 1 = [1, . . . , 1] T ∈ R N and α lin = 1. Concerning the update of the smoothing parameter µ, we set α = 10 -1 . Table 3 summarizes the results. To achieve the stopping criterion (6.12), 11

smoothing iterations and 20 cumulated Newton iterations are needed (CPU time: 5.0s). In Figure 10, we plot the estimators (6.9) as a function of the cumulated Newton iterations (left), the smoothing iterations (middle), and the relative norm of the residual vector as a function of the smoothing iterations (right). The behavior of R AIP (X j,k ) r in Figure 10 appears a bit different from its behavior in Figure 4. This is related to the fact that the relative norm of the total residual given by (2.6) includes C(X) in the adaptive smoothing Newton method, whereas in this adaptive interior-point method, the relative norm of the total residual given by (6.10) includes K(X)G(X). Smoothing iteration Relative norm of the linearization residual vector R AIP (X j,k ) r

µ j Niter η j,k lin,AIP η j,k sm,AIP R AIP (X j,k ) r 1e+00 2
Figure 10: [Adaptive interior-point method] Estimators (6.9) as a function of cumulated Newton iterations (left). Estimators (6.9) (middle) and relative norm of the total residual vector (6.10) (right) as a function of smoothing iterations j at convergence of the linearization solver.

Adaptive inexact interior-point method

Let us now present the numerical results of the adaptive inexact interior-point method, detailed in Section 5. We employ Algorithm 2 with the GMRES algebraic solver and an ILU preconditioner. The parameters in Algorithm 2 are set as µ 1 = [1, . . . , 1] T ∈ R N , ε = 10 -5 , α alg = 1, α lin = 1, and α = 0.1. The restart parameter of restarted GMRES is chosen equal to 300. From Table 4, we can see that the method converged after 8 smoothing iterations, 20 cumulated linearization iterations, and 760 cumulated GMRES iterations. Figure 11, left, displays the curves of the estimators (5.9) as a function of the smoothing iteration. One can see that the estimators satisfy the adaptive stopping criteria incorporated in Algorithm 2. In Figure 11, right, the estimators are shown as a function of cumulated Newton iterations, at convergence of the linear solver.

µ j Niter Giter η j,k,i lin,AIIP η j,k,i sm,AIIP η j,k,i alg,AIIP R AIIP (X 4: [Adaptive inexact interior-point method, Algorithm 2] Number of cumulated Newton and GMRES iterations, estimators (5.9), and relative norm of the total residual vector (5.8) at each smoothing iteration j, at convergence of the algebraic and linearization solvers.

Comparison of the methods

This section is devoted to compare the semismooth Newton method (SSN), semismooth Newton method with pathfollowing (SSN-pf), nonparametric interior-point method (IP), adaptive smoothing Newton method (ASN), and adaptive interior-point method (AIP). For this purpose, we introduce a unified residual vector, for independent of the way the nonlinear complementarity constraints are reformulated. The stopping criterion of the nonlinear solver for the classical methods (SSN, SSN-pf, IP) is on the relative unified residual R unif (X k ) r lying below 10 -8 . Regarding the adaptive methods (ASN, AIP), to stop the nonlinear solver, we use the adaptive stopping criteria given respectively in (6.5) and (6.11). To stop the smoothing iterations, R unif (X j,k ) r is requested to become smaller than 10 -8 . 2)), semismooth Newton method with a path-following strategy, nonparametric interior-point method, adaptive interior-point method, and adaptive smoothing Newton method (smoothed F-B function (3.2))] Number of cumulated Newton iterations (left) and CPU time (right) as a function of the number of mesh elements, employing a stopping criterion on the relative norm of the unified residual vector (6.13).

V ∈ R n R unif (V ) :=     F -EV min(0, K(V )) min(0, G(V )) K(V ) • G(V )     , (6.
In Figure 12, we plot the cumulated number of the Newton iterations (left) and the CPU time (right) required by each method, as a function of the number of mesh elements. It is clearly seen that the semismooth Newton method (green curve) is typically more costly, both in terms of the required number of iterations and the CPU time, in comparison with the other methods. Precisely, we can observe an important gain between the semismooth Newton method (green curve) and the adaptive smoothing Newton method (purple curve). Moreover, as we can remark from the red curve, the combination of a path-following strategy to the semismooth Newton method seems to be efficient. Finally, one does not see a remarkable difference between the results of the nonparametric interior-point method (cyan curve) and the adaptive interior-point method (black curve) in this test case.

Adaptive smoothing Newton algorithm

Let ε > 0 be the desired relative tolerance, α lin ∈ ]0, 1] be the desired relative size of the linearization error, and α ∈ ]0, 1[ the smoothing decrease parameter. The unsteady adaptive smoothing Newton algorithm reads as follows:

Algorithm 3 Unsteady adaptive smoothing Newton algorithm Initialization: Fix ε > 0, α ∈ ]0, 1[, and α lin ∈ ]0, 1]. Set ν := 1 and t ν := 0. Choose X ν,0 ∈ R n . Time loop 1. Fix µ jν > 0 and set j := 1.

Smoothing loop

2.1 Set X ν,j,0 := X ν,0 and k := 1.

Newton linearization loop

2.2.1 From X ν,j,k-1 define A ν,j,k-1 µ jν ∈ R n,n and B ν,j,k-1 µ jν ∈ R n given by (7.6).
2.2.2 Find X ν,j,k solution to the linear system

A ν,j,k-1 µ jν X ν,j,k = B ν,j,k-
max η ν,j,k sm,ASN , R(X ν,j,k ) r < ε. (7.8)
If satisfied, set j := j and stop. If not, set j := j + 1, X ν,j,0 := X ν,j-1,k , and µ jν := αµ (j-1)ν . Then set k := 1 and go to 2.2.1.

If ν = N t , stop. If not, set ν := ν + 1, j = 1, X ν,j,0 := X ν-1,j , and t ν := τ ν + t ν-1 . Then set µ jν = µ j ν-1 , k = 1, and go to 2.2.1.

Description of Algorithm 3

For the first time step ν = 1, starting with an initial approximation X ν,0 ∈ R n and an initial smoothing parameter µ ν,1 > 0, we solve the smoothed nonlinear system (7.5) by the Newton linearization solver, and decrease the smoothing parameter µ ν,j at each smoothing step j, until the stopping criterion (7.8) on the smoothing estimator or the relative norm of the total residual vector is satisfied at step j. Then, we continue the time loop, for 2 ≤ ν ≤ N t , starting for j = 1 with X ν,j,0 := X ν-1,j and µ jν := µ j ν-1 , until satisfying the stopping criterion (7.8).

Numerical results

We consider a homogeneous porous medium in one dimension, supposed to be horizontal with length 2m, and a uniform spatial mesh with N = 1000 elements. The final time of simulation is t F = 100s, and the time step is constant τ ν = 10s. We assume that the medium is initially saturated with liquid, S l = 1, and containing no hydrogen, χ l h = 0, on which we impose an injection of gas (hydrogen), constant in time, in the first cell of the mesh.

The initial conditions are S l,ν=0 = 1, P l,ν=0 = 10 6 Pa, and χ l,ν=0 h = 0.

Semismooth Newton method. We begin by employing the semismooth Newton method presented in Section 2, with the min function (2.1) to solve the nonlinear system (7.2). On each time step ν ≥ 1, we request the relative norm of the total residual vector R(X ν,k ) given by (7.3) to drop below 10 -4 .

In Figure 13, the evolution of R(X ν,k ) r is shown at each time step. 31 cumulated Newton iterations are needed. Cumulated Newton iteration Relative norm of the total residual vector R(X ν,k ) r From Figure 14, one can see that at the first time step ν = 1 and at each smoothing step j ≤ 4, the linearization estimator decreases until lying below the smoothing estimator. The smoothing iterations are thus stopped in the first possibility according to the stopping criterion (7.7). On the other hand, at the 5 th smoothing step, η ν,j,k lin,ASN is smaller than the fixed tolerance but not smaller than η ν,j,k sm,ASN . Even after additional Newton iterations at this smoothing step, we will have the same observation. This justifies the modification applied in the adaptive stopping criterion (7.7). In Figure 14, right, we report the estimators and R(X 1,j,k ) r as a function of cumulated Newton iteration for ν = 1. The stopping criterion (7.8) is satisfied after 5 cumulated smoothing iterations, and 10 cumulated Newton iterations. Then, as presented in Table 5, starting at the second time step (ν = 2) with µ jν = 10 -5 , the smoothing parameter does not decrease since the stopping criterion (7.8) is satisfied at each time step after one smoothing step only. To reach the end of the simulation, 9 cumulated smoothing steps and 31 cumulated linearization steps are needed.

As a conclusion, the results confirm the expected behavior of Algorithm 3 featuring an adaptive stopping criterion for the nonlinear solver. In this case, though, the stopping criteria in the adative smoothing Newton method do not bring the number of iterations down since the semismooth Newton method already behaves very well here.

Conclusion and outlook

In this work, we have considered nonlinear algebraic systems with inequalities in a form of complementarity constraints. We have considered some existing methods, like the semismooth Newton method, possibly combined with a path-following strategy, or a nonparametric interior-point method. Our goal was to propose a systematic way to drive such methods with adaptive stopping criteria and possibly inexact algebraic solvers. We have achieved this by a reformulation of the complementarity constraints using a smoothed function and a posteriori error estimate that enables to distinguish the different error components. Numerical experiments confirmed that the proposed adaptive approaches yield significant computational savings compared to some standard approaches from literature. Moreover, their numerical performance seems to be notably good across a range of test problems. In [START_REF] Ben Gharbia | Adaptive inexact smoothing Newton method for a nonconforming discretization of a variational inequality[END_REF], we also take into account the discretization error of the considered problem, enabling to adaptively stop the outer smoothing loop in Algorithm 1, and employ the method to solve more involved problems.

Figure 1 :

 1 Figure 1: Left: Absolute value function | • | and smoothed absolute value function | • | µ . Right: Fischer-Burmeister function CFB (•) and smoothed Fischer-Burmeister function CFBµ (•), for different values of the smoothing parameter µ.

Figure 2 :

 2 Figure 2: [Semismooth Newton method, F-B function (2.2), stopping criterion (6.2)] Relative norm of the total residual vector (2.6) as a function of semismooth Newton iterations.

Figure 3 :

 3 Figure 3: [Adaptive smoothing Newton method, smoothed F-B function (3.2), classical and adaptive stopping criteria (6.4) and (6.5)] Relative norm of the total residual vector (2.6) and estimators (6.3) as a function of Newton iterations k, at a specific smoothing iteration j = 1 (µ 1 = 1), left, and at j = 3 (µ 3 = 10 -2 ), right.

Figure 4 :

 4 Figure 4: [Adaptive smoothing Newton method, smoothed F-B function (3.2), adaptive stopping criterion (6.5)] Estimators (6.3) as a function of cumulated Newton iterations (left). Estimators (6.3) (middle) and relative norm of the total residual vector (2.6) (right) as a function of smoothing iterations j at convergence of the linearization solver.

Figure 5 :

 5 Figure 5: [Semismooth Newton method (with and without a path-following strategy) and adaptive smoothing method] Cumulated number of Newton iterations (left) and CPU time (right) as a function of the number of mesh elements.

Figure 6

 6 Figure 6 depicts the evolution of the algebraic and linearization estimators and the GMRES relative residual during the algebraic resolution, for specific smoothing step j and linearization step k. For j = 2 and k = 2, we see that 22

Figure 7 ,Figure 6 :

 76 Figure 7, left, shows the evolution of the estimators during smoothing iterations, at convergence of the nonlinear and linear solvers. As expected, the estimators decrease when µ decreases at each smoothing step. In the middle

Figure 7 :

 7 Figure 7: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Estimators (3.10) as a function of smoothing iterations j at convergence of the algebraic and linearization solvers, left. Estimators as a function of cumulated Newton iterations at convergence of the algebraic solver, middle. Estimators as a function of cumulated GMRES iterations during the first two smoothing iterations (j = 1 and j = 2), right.

Number of mesh elements ( n 3 )Figure 8 :

 38 Figure 8: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Ratio between: the number of algebraic iterations (left) and CPU time (right) needed by the classical stopping criterion (6.7) to converge to the number and time needed by the adaptive stopping criterion (3.11a), as a function of the number of mesh elements.

Figure 9 :

 9 Figure 9: [Nonparametric interior-point method, stopping criterion (6.8)] Relative norm of the linearization residual vector (6.8) as a function of Newton iterations.

13

 13

  )

Figure 11 :

 11 Figure 11: [Adaptive inexact interior-point method, Algorithm 2] Estimators (5.9) as a function of smoothing iterations j at convergence of the algebraic and linearization solvers (left). Estimators as a function of cumulated Newton iterations k at convergence of the algebraic solver (right).

Figure 12 :

 12 Figure12: [Semismooth Newton method (F-B function (2.2)), semismooth Newton method with a path-following strategy, nonparametric interior-point method, adaptive interior-point method, and adaptive smoothing Newton method (smoothed F-B function (3.2))] Number of cumulated Newton iterations (left) and CPU time (right) as a function of the number of mesh elements, employing a stopping criterion on the relative norm of the unified residual vector (6.13).

Figure 13 : 2 1

 1321 Figure 13: [Semismooth Newton method, min function (2.1)] Relative norm of the total residual vector (7.3) as a function of cumulated Newton iterations along the time steps ν.

Figure 14 :

 14 Figure 14: [Adaptive smoothing Newton method, smoothed min function (3.1), Algorithm 3] Estimators (6.3) and relative norm of the total residual vector (7.3) at the first time step ν = 1 as a function of smoothing iterations j, at convergence of the linearization solver (ν = 1 fixed, j varies, k = k), left, and of cumulated Newton iterations, right, (ν = 1 fixed, j and k vary).

Table 2 :

 2 [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Number of Newton iterations and cumulated GMRES iterations, estimators (3.10), and relative norm of the total residual vector (2.6) at each smoothing iteration j, at convergence of the algebraic and linearization solvers.

			8	2.16e+03 4.24e+03 1.80e+00	2.19e+03
	1e-01	4	34	5.95e+01 2.31e+02 1.89e-02	1.80e+02
	1e-02	3	391	1.54e+01 1.73e+01 1.41e-02	6.75e+00
	1e-03	4	198	5.04e-01	8.16e-01	4.60e-04	5.95e-01
	1e-04	10	796	7.99e-03	3.53e-02	5.58e-06	3.43e-02
	1e-05	10	684	8.54e-04	2.12e-03	7.61e-07	1.94e-03
	1e-06	4	513	9.03e-05	1.48e-04	7.42e-08	1.05e-04
	1e-07	3	3375	6.04e-06	8.14e-06	4.27e-09	4.26e-06

Table 3 :

 3 [Adaptive interior-point method] Number of Newton iterations, estimators (6.9), and relative norm of the total residual vector (6.10) at each smoothing step j, at convergence of the linearization solver.

			1.11e+01 2.00e+01	3.00e+01
	1e-01	2	1.24e+00 2.00e+00	3.20e+00
	1e-02	2	1.15e-01	2.00e-01	3.11e-01
	1e-03	2	6.51e-03	2.00e-02	2.43e-02
	1e-04	2	3.38e-04	2.00e-03	2.14e-03
	1e-05	1	1.58e-04	2.00e-04	2.82e-04
	1e-06	2	3.67e-06	2.00e-05	2.10e-05
	1e-07	2	1.00e-07	2.00e-06	2.02e-06
	1e-08	1	1.86e-07	2.00e-07	3.84e-07
	1e-09	2	9.33e-10	2.00e-08	2.01e-08
	1e-10	2	2.55e-11	2.00e-09	2.00e-09

  j,k,i )

							r
	1e+00	3	7	1.15e+01 2.00e+01 3.36e+00	5.59e+00
	1e-01	2	12	5.44e-01 2.00e+00 1.78e-02	2.00e+00
	1e-02	3	20	9.75e-02	2.00e-01	2.80e-02	2.04e-01
	1e-03	3	29	4.82e-03	2.00e-02	1.74e-03	2.01e-02
	1e-04	3	56	2.52e-04	2.00e-03	2.19e-04	2.01e-03
	1e-05	2	62	1.77e-04	2.00e-04	1.08e-04	2.29e-04
	1e-06	2	110	1.49e-05	2.00e-05	1.42e-05	2.46e-05
	1e-07	2	464	1.34e-06	2.00e-06	1.16e-06	2.31e-06
	Table

  Compute the estimators and check the stopping criterion for the nonlinear solver

		µ jν	1	.
	2.2.3 η ν,j,k lin,ASN < α lin η ν,j,k sm,ASN	or	η ν,j,k lin,ASN < ε .	(7.7)

If satisfied, set k := k and stop. If not, set k := k + 1 and go to 2.2.1.

2.3

Check the stopping criterion for the smoothing iterations in the form:

Table 5 :

 5 ν,j,k ) [Adaptive smoothing Newton method, smoothed min function (3.1), Algorithm 3] Relative norm of the total residual vector (7.3) and estimators (6.3) at each time step ν, at convergence of the linearization solver.

	r

* This project has partly received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).

7 Numerical experiments: Two-phase flow with phase transition

The second model problem that we consider in our numerical tests is a two-phase flow model (liquid-gas) with phase transition in porous media following [START_REF] Gharbia | Résolution de problèmes de complémentarité. Application à un écoulement diphasique dans un milieu poreux[END_REF][START_REF] Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF][START_REF] Gharbia | Study of compositional multiphase flow formulation using complementarity conditions[END_REF]. Each of the liquid phase, denoted by l, and the gas phase, denoted by g, is composed of two components, water and hydrogen, denoted respectively by w and h.

Problem statement

The problem at hand can be formulated as a system of nonlinear partial differential equations with nonlinear complementarity constraints at each time step τ ν . Let T h be the spatial mesh, we denote respectively by S ν K , P ν K , and χ ν K the discrete elementwise unknowns approximating the values of the saturation S l , the pressure P l , and the molar fraction of hydrogen in the liquid phase χ l h in the element K ∈ T h and on time step 1 ≤ ν ≤ N t . Let N be the number of elements in the mesh T h . If one introduces the appropriate nonlinear function H ν c,K : R 3N → R, c ∈ {w, h}, and suitable functions F K : R 3 → R and G K : R 3 → R, the discrete problem written elementwise consists in finding

, where n = 3N, and

The formulation (7.1) allows to model the transition from a single-phase flow to a two-phase flow during the appearance and disappearance of the gas phase and vice versa. As an example, a detailed finite volume discretization can be found in [7, Section 3.2]. The first 2N lines of system (7.1) can be written globally as

where H n : R 3N → R 2N is defined elementwise by (7.1a).

Considering a C-function

. This leads us to apply a semismooth Newton method to find a solution for problem (7.1) written as

The total residual vector R(V ) of problem (7.2) is thus given by

Adaptive smoothing Newton method

We introduce a function

where C ν µ is a smoothed C-function. Line (7.1b) can be approximated as a smoothed nonlinear equation C ν µ (X ν ) = 0, making it possible to apply the standard Newton method to solve the resulting nonlinear system in the form: Find X ν,j ∈ R 3N at each time step ν,

At each time step 1 ≤ ν ≤ N t , each smoothing step j ≥ 1, and each linearization step k ≥ 1, fixing X ν,j,0 ∈ R n , we try to approach the solution of problem (7.4) by finding a solution X ν,j,k ∈ R n such that

where the Jacobian matrix A ν,j,k-1 µ jν ∈ R n,n and the right-hand side vector B ν,j,k-1 µ jν ∈ R n are defined by

B ν,j,k-1 µ jν := J H ν (X ν,j,k-1 )X ν,j,k-1 -H ν (X ν,j,k-1 ) J C ν µ jν (X ν,j,k-1 )X ν,j,k-1 -C ν µ jν (X ν,j,k-1 ) ,

with J H ν (X ν,j,k-1 ) and J C ν µ jν (X ν,j,k-1 ) the Jacobian matrices of the function H ν and the smoothed function C ν µ jν , respectively, at the point X ν,j,k-1 obtained by a Newton linearization.