Laurence Halpern
email: halpern@math.univ-paris13.fr

Jérémie Szeftel
email: jszeftel@math.princeton.edu.

Nonlinear Nonoverlapping Schwarz Waveform Relaxation for Semilinear Wave Propagation

We introduce a non-overlapping variant of the Schwarz waveform relaxation algorithm for semilinear wave propagation in one dimension. Using the theory of absorbing boundary conditions, we derive a new nonlinear algorithm. We show that the algorithm is well-posed and we prove its convergence by energy estimates and a Galerkin method. We then introduce an explicit scheme. We prove the convergence of the discrete algorithm with suitable assumptions on the nonlinearity. We finally illustrate our analysis with numerical experiments.

Introduction

Schwarz waveform relaxation is a new class of algorithms for domain decomposition in the frame of time dependant partial differential equations. They are well-adapted to evolution problems, designed to solve the equations separately on each spatial subdomain on the whole time interval , exchanging informations on the space-time boundary of the subdomains, overlapping or not [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF]. In particular for wave equations, it is of great importance, due to numerical dispersion, to be able to handle local time and space meshes, and this is allowed by the present method. We presented the method for the linear wave equation in [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF] and [START_REF] Gander | Absorbing boundary conditions for the wave equation and parallel computing[END_REF] . When using overlapping subdomains and "classical" Schwarz waveform relaxation -by a Dirichlet exchange of informations on the boundary-the so defined algorithm converges in a finite number of iterations, inversely proportional to the size of the overlap, which can be penalizing. We introduced optimized transmission conditions, relying on the theory of absorbing boundary conditions, which improve drastically the convergence of the algorithm.

Very little has been done so far about nonlinear Schwarz algorithms. An analysis of the classical Schwarz waveform relaxation algorithm was performed in [START_REF] Gander | Overlapping Schwarz waveform relaxation for convection dominated nonlinear conservation laws[END_REF] for a conservation law. The goal of this paper is to define new Schwarz waveform relaxation algorithms for the semilinear wave equation. We introduce two nonoverlapping algorithms. The first one referred to as linear uses the absorbing boundary condition of the linear problem whereas the second one referred to as nonlinear uses the nonlinear absorbing boundary conditions designed by J. Szeftel in [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF].

In Section 2, we introduce the definitions of the algorithms.

In Section 3 , we prove the algorithms to be well-posed. For the precise analysis, we use a fixed point algorithm with regularity estimates on a linear problem.

In Section 4, we prove the convergence of the algorithms. The proof is an extension of a clever trick in [START_REF] Lions | On the Schwarz alternating method[END_REF], already used for linear algorithms, either hyperbolic or parabolic (see [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF]). However the nonlinearity requires a very fine analysis.

In Section 5, we design discrete Schwarz waveform relaxation algorithms. In each subdomain, the interior scheme is the usual leapfrog scheme for the linear part, with a downwinding in time for the nonlinear part. The exchange of informations on the boundary is naturally taken into account by a finite volume strategy. In Section 6 we study the convergence of the algorithms, by discrete energy estimates.

As it is always the case for nonlinear problems, the well-posedness and convergence results hold only locally in time. Therefore numerical experiments are very important to bypass the limitations of the theory. We present the results in Section 7, showing in particular that our nonlinear algorithm gives optimal results within a large class of algorithms. Remark Due to the complexity of the mathematical theory, we restrain ourselves to the one dimensional case. The multidimensional study contains additional difficulties due to the geometry and should be the heart of a forthcoming paper.

Problem Description

We consider the second order semilinear wave equation in one dimension,

(∂ 2 t -∂ 2 x ) U = f (U, ∂ t U , ∂ x U ) (1)
on the domain R × (0, T ) with initial conditions U (•, 0) = p, ∂ t U (•, 0) = q.

Absorbing Boundary Conditions for the Semilinear Wave Equation

The question of absorbing boundary conditions arises when one wants to make computations on an unbounded domain: a bounded computational domain is introduced, on the boundary of which boundary conditions must be prescribed. These boundary conditions must be absorbing to the waves leaving the domain. A whole strategy has been designed by Engquist and Majda for linear problems with variable coefficients, using pseudo-differential operators [START_REF] Engquist | Radiation boundary conditions for acoustic and elastic calculations[END_REF]. Recently it has been extended to nonlinear operators by J. Szeftel, in particular for the semilinear wave equation [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF], using the paradifferential calculus of [START_REF] Bony | Calcul symbolique et propagation des singularités pour les équations aux dérivées partielles non linéaires[END_REF] and [START_REF] Sablé-Tougeron | Régularité microlocale pour des problèmes aux limites non linéaires[END_REF]. We introduce a family of operators

B ± (g ± )u = ∂ t u ± ∂ x u + g ± (u), (2)
for C ∞ functions g ± such that g ± (0) = 0. The linear absorbing boundary operators are given by g ± = 0.

In the case where f (u, u t , u x ) = f 1 (u) + f 2 (u)u t + f 3 (u)u x with f j in C ∞ (R), 1 ≤ j ≤ 3, and f 1 (0) = 0, the following nonlinear boundary operators are given in [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF]:

g + (u) := - 1 2 u 0 (f 2 -f 3 )(ξ)dξ, g -(u) := - 1 2 u 0 (f 2 + f 3 )(ξ)dξ. (3)
We replace the problem on the domain R by a boundary value problem in Ω 0 = (a, b):

(∂ 2 t -∂ 2 x )ū = f (ū, ∂ t ū, ∂ x ū) in Ω 0 × (0, T ), B -(g -)ū(a, •) = 0, B + (g + )ū(b, •) = 0, (4)
with initial values p and q. Such boundary conditions give well-posed initial boundary value problems, and are absorbing provided the intial data be compactly supported in Ω 0 , see [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF]. Following the strategy in [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF], we use such absorbing operators for domain decomposition.

A General Non-Overlapping Schwarz Waveform Relaxation Algorithm

We decompose the domain (a, b) into I non overlapping subdomains Ω i = (a i , a i+1 ), a j < a i for j < i and a 1 = a, a I+1 = b, and we introduce a general non overlapping Schwarz waveform relaxation algorithm. An initial guess {h ±,0 i } 1≤i≤I+1 is given. For k ≥ 1, one step of the algorithm is

  (∂ 2 t -∂ 2 x )u k i = f (u k i , ∂ t u k i , ∂ x u k i ) in Ω i × (0, T ), u k i (•, 0) = p, ∂ t u k i (•, 0) = q in Ω i , B -(g -)u k i (a i , •) = h -, k-1 i , B + (g + )u k i (a i+1 , •) = h +, k-1 i in (0, T ), h -, k i = B -(g -)u k i-1 (a i , •), h +, k i = B + (g + )u k i+1 (a i+1 , •) in (0, T ), (5)
where B ± are given in [START_REF] Chemin | Mécanique des fluides[END_REF]. For ease of notations, we defined here h ±, k 1 = 0 and h ±, k I+1 = 0, so that the index i in (5) ranges from i = 1, 2, . . . , I. In the sequel, we call linear transmission condition the choice g ± = 0 and nonlinear transmission condition the choice (3). For the classical linear homogeneous wave equation, it has been proved in [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF] that the algorithm converges optimally if T is small enough (which means in two iterations, independently of the number of subdomains), and the transmission operators B ± are given by B ± = ∂ t ± ∂ x . This behavior is due to the finite speed of propagation, together with the fact that these operators are the exact Dirichlet Neumann operators in this case. In the nonlinear case, the propagation still takes place with the finite speed, but we can use only approximate Dirichlet Neumann operators. Therefore the classical Schwarz algorithm with overlap is still convergent, and for our nonoverlapping nonlinear algorithms, we will use energy estimates.

Well-posedness For The Subproblems

The study of the nonlinear problem relies on an iterative linear scheme. Therefore a first step for the definition of the algorithm is the study of the nonhomogeneous initial boundary value problem for a general domain Ω = (a -, a + ),

(∂ 2 t -∂ 2 x )u = f (u, ∂ t u, ∂ x u) in Ω × (0, T ), B -(g -)u(a -, •) = h -, B + (g + )u(a + , •) = h + , (6)
with initial values p and q. We will use for j ≤ 2 the spaces

V j (Ω, T ) = {u ∈ L ∞ (0, T ; L 2 (Ω)), ∂ α u ∈ L ∞ (0, T ; L 2 (Ω)), |α| ≤ j}. (7)
In formula [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF], α is a 2-index in N 2 , the first coordinate in α stands for the time, and the second one stands for the space, so for instance

∂ α u = ∂ tx u for α = (1, 1). V j (Ω, T ) is equipped with the norm u Vj(Ω,T ) = max |α|≤j ∂ α u L ∞ (0,T ;L 2 (Ω)) .
Theorem 3.1 Let p in H 2 (Ω) and q in H 1 (Ω). There exists a time T * such that for any T ≤ T * , for h ± in H 1 (0, T ) with the compatibility conditions

h ± (0) = q(a ± ) ± p ′ (a ± ) + g ± (p(a ± )), (8)
(6) has a unique solution u in V 2 (Ω, T ), with ∂ t u(a ± , •) and ∂ x u(a ± , •) in H 1 (0, T ). Furthermore there exists a positive real number C * such that

u 2 V2(Ω,T ) + ± |α|=1 ∂ α u(a ± , •) 2 H 1 (0,T ) ≤ C * ( p 2 H 2 (Ω) + q 2 H 1 (Ω) + ± h ± 2 H 1 (0,T ) ), (9)
where T * and C * depend on the data p, q, f, g ± , h ± .

This result has been first proved in [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF] with homogeneous boundary conditions (i.e. h ± = 0). The additional difficulty comes from the boundary conditions, and we give here the main steps of the proof. It relies on the construction of a sequence of linear problems of the form:

∂ 2 t ũ -∂ 2 x ũ + ũ = F in Ω × (0, T ), (∂ t ũ -∂ x ũ)(a -, •) = H -, (∂ t ũ + ∂ x ũ)(a + , •) = H + . ( 10
)
Proposition 3.2 Let p in H 2 (Ω) and q in H 1 (Ω). For any positive time T , let F in H 1 ((0, T ) × Ω), and H ± in H 1 (0, T ) with the compatibility conditions

H ± (0) = q(a ± ) ± p ′ (a ± ). ( 11
)
Then, [START_REF] Szeftel | Absorbing boundary conditions for non linear partial differential equations[END_REF] with initial data p and q has a unique solution ũ in V 2 (Ω, T ), with ∂ t ũ(a ± , •) and ∂ x ũ(a ± , •) in H 1 (0, T ). Moreover we have the following bounds on the solution

ũ 2 V2(Ω,T ) + ± ∂ t ũ(a ± , •) 2 H 1 (0,T ) + ∂ x ũ(a ± , •) 2 H 1 (0,T ) ≤ C 1 e T F 2 H 1 (0,T ;L 2 (Ω)) + F (•, 0) 2 L 2 (Ω) + ± H ± 2 H 1 (0,T ) + p 2 H 2 (Ω) + q 2 H 1 (Ω) , (12)
where C 1 is a universal constant.

Proof We start with the a priori estimates. We multiply (10) by ∂ t ũ and integrate by parts in Ω:

1 2 d dt ∂ t ũ(•, t) 2 L 2 (Ω) + ũ(•, t) 2 H 1 (Ω) + (∂ t ũ(a -, t)) 2 + (∂ t ũ(a + , t)) 2 = (F (•, t), ∂ t ũ(•, t)) L 2 (Ω) + H -(t)∂ t ũ(a -, t) + H + (t)∂ t ũ(a + , t).
Using the Cauchy-Schwarz inequality on the right hand side, together with the inequality αβ ≤ 1 2 α 2 + 1 2 β 2 for all α, β ∈ R, and finally integrating in time, we obtain

∂ t ũ(•, t) 2 L 2 (Ω) + ũ(•, t) 2 H 1 (Ω) + t 0 [(∂ t ũ(a -, s)) 2 + (∂ t ũ(a + , s)) 2 ] ds ≤ t 0 [ ∂ t ũ(•, s) 2 L 2 (Ω) ds + t 0 F (•, s) 2 L 2 (Ω) ds + q 2 L 2 (Ω) + p 2 H 1 (Ω) + t 0 [(H -(s)) 2 + (H + (s)) 2 ] ds.
By Gronwall Lemma, we deduce that

∂ t ũ(•, t) 2 L 2 (Ω) + ũ(•, t) 2 H 1 (Ω) + t 0 [(∂ t ũ(a -, s)) 2 + (∂ t ũ(a + , s)) 2 ] ds ≤ e T F 2 L 2 ((0,T )×Ω) + q 2 L 2 (Ω) + p 2 H 1 (Ω) + ± H ± 2 L 2 (0,T ) , which gives max |α|≤1 ∂ α ũ 2 L ∞ (0,T ;L 2 (Ω)) + ∂ t ũ(a -, •) 2 L 2 (0,T ) + ∂ t ũ(a + , •) 2 L 2 (0,T ) ≤ e T F 2 L 2 ((0,T )×Ω) + q 2 L 2 (Ω) + p 2 H 1 (Ω) + ± H ± 2 L 2 (0,T ) . (13)
Differentiating in time in [START_REF] Szeftel | Absorbing boundary conditions for non linear partial differential equations[END_REF], we now apply (13) to ∂ t ũ, and obtain:

max |α|≤1 ∂ α ∂ t ũ 2 L ∞ (0,T ;L 2 (Ω)) + ∂ 2 t ũ(a -, •) 2 L 2 (0,T ) + ∂ 2 t ũ(a + , •) 2 L 2 (0,T ) ≤ e T ( ∂ t F 2 L 2 ((0,T )×Ω) + ∂ 2 t u(•, 0) 2 L 2 (Ω) + q 2 H 1 (Ω) + ± ∂ t H ± 2 L 2 (0,T ) ). ( 14
)
We must estimate ∂ 2 t ũ(•, 0) in the righthand side of (14). We multiply (10) by ∂ 2 t ũ, integrate in space and evaluate at time 0:

∂ 2 t ũ(•, 0) 2 L 2 (Ω) + (∂ x p, ∂ x ∂ 2 t ũ(•, 0)) + (p, ∂ 2 t ũ(•, 0)) + q(a -)∂ 2 t ũ(a -, 0) + q(a + )∂ 2 t ũ(a + , 0) = H -(0)∂ 2 t ũ(a -, 0) + H + (0)∂ 2 t ũ(a + , 0) + (F (•, 0), ∂ 2 t ũ(•, 0)).
We integrate by parts in the second term, and rewrite the equality as

∂ 2 t ũ(•, 0) 2 L 2 (Ω) -(∂ 2 x ũ(•, 0), ∂ 2 t ũ(•, 0)) + (p, ∂ 2 t ũ(•, 0)) = (H -(0) -q(a -) + ∂ x p(a -))∂ 2 t ũ(a -, 0) + (H + (0) -q(a + ) -∂ x p(a + ))∂ 2 t ũ(a + , 0) + (F (•, 0), ∂ 2 t ũ(•, 0)).
The boundary terms on the right-hand side vanish by the compatibility conditions, and we get

∂ 2 t ũ(•, 0) 2 L 2 (Ω) = (∂ 2 x p -p + F (•, 0), ∂ 2 t ũ(•, 0)).
Using the Cauchy-Schwarz Lemma, we obtain

∂ 2 t ũ(•, 0) L 2 (Ω) ≤ ∂ 2 x p -p + F (•, 0) L 2 (Ω) .
We replace the term ∂ 2 t ũ(•, 0) L 2 (Ω) in (13), and we deduce the second a priori estimate:

max |α|≤1 ∂ α ∂ t ũ 2 L ∞ (0,T ;L 2 (Ω)) + ∂ 2 t ũ(a -, •) 2 L 2 (0,T ) + ∂ 2 t ũ(a + , •) 2 L 2 (0,T ) ≤ 3 e T ∂ t F 2 L 2 ((0,T )×Ω) + F (•, 0) 2 L 2 (Ω) + p 2 H 2 (Ω) + q 2 H 1 (Ω) + ± ∂ t H ± 2 L 2 (0,T ) . ( 15
)
We still need to estimate the mixed derivatives ∂ xx ũ in the interior and ∂ xt ũ on the boundaries. We use the equation, which gives in the interior

∂ xx ũ L ∞ (0,T ;L 2 (Ω)) ≤ ∂ tt ũ L ∞ (0,T ;L 2 (Ω)) + ũ L ∞ (0,T ;L 2 (Ω)) + F L ∞ (0,T ;L 2 (Ω)) .
We now introduce the inequality

F 2 L ∞ (0,T ;L 2 (Ω)) ≤ 2( F (•, 0) 2 L 2 (Ω) + T ∂ t F 2 L 2 (0,T ;L 2 (Ω)) ),
and by ( 15) and (13), we get, using that e T ≥ 1 and e T ≥ T ,

∂ xx ũ 2 L ∞ (0,T ;L 2 (Ω)) ≤ 15e T F 2 H 1 (0,T ;L 2 (Ω)) + F (•, 0) 2 L 2 (Ω) + p 2 H 2 (Ω) + q 2 H 1 (Ω) + ± H ± 2 H 1 (0,T ) .
As for the boundary term, we get for instance on the left boundary

∂ tx ũ(a -, •) L 2 (0,T ) ≤ ∂ tt ũ(a -, •) L 2 (0,T ) + ∂ t H - L 2 (0,T ) .
Squaring the inequality, and adding the term coming from the right boundary leads to

± ∂ tx ũ(a ± , •) 2 L 2 (0,T ) ≤ 2 ± ∂ tt ũ(a ± , •) 2 L 2 (0,T ) + ∂ t H ± 2 L 2 (0,T ) ,
which provides the last estimate announced in the proposition. The well-posedness is then derived in a standard way by the Galerkin method.

The solution ū of the nonlinear subdomain problem is now defined through an iterative scheme. The initial guess is ū0 = p. At step k, ūk being known, we define

F (w) = w + f (w, ∂ t w, ∂ x w), G ± (w) = g ± (w(a ± , •)), H ± (w) = h ± -G ± (w). ( 16
)
ūk+1 is the solution of the linear initial boundary value problem [START_REF] Szeftel | Absorbing boundary conditions for non linear partial differential equations[END_REF] with data

f k = F (ū k ), g ± k = G ± (ū k ), h ± k = H(ū k )
, and initial data p and q. The proof of convergence for the sequence ūk is written in details in [START_REF] Szeftel | A non linear approach to absorbing boundary conditions for the semilinear wave equation[END_REF]. The uniqueness follows from the result: Lemma 3.3 There exists a real positive increasing function θ such that, for any time T , for any v in

V 2 (Ω, T ), ∂ α F (v) is in V 1 (Ω, T ), G ± (v) and H(v) are in L ∞ (0, T ). Moreover, for v 1 , v 2 in V 2 (Ω, T ), we have H ± (v 1 ) -H ± (v 2 ) L ∞ (0,T ) ≤ θ( v 1 V2(Ω,T ) + v 2 V2(Ω,T ) ) v 1 -v 2 V1(Ω,T ) , F (v 1 ) -F(v 2 ) V1(Ω,T ) ≤ θ( v 1 V2(Ω,T ) + v 2 V2(Ω,T ) ) v 1 -v 2 V2(Ω,T ) , (17)
As a consequence, we have the well-posedness of problem (4).

Corollary 3.4 Let p in H 2 0 (Ω) and q in H 1 0 (Ω). There exists a time T * 0 such that for any T ≤ T * 0 , (4) has a unique solution ū in V 2 (Ω, T ), with ∂ α ū(a, •) and ∂ α ū(b, •) in H 1 (0, T ) for |α| = 1. Furthermore there exists a positive real number C * depending only on the size of Ω such that

ū 2 V2(Ω,T ) + |α|=1 ∂ α ū(a, •) 2 H 1 (0,T ) |α|=1

Convergence of The Algorithm

We now study the convergence of the Schwarz waveform Relaxation Algorithm [START_REF] Gander | Absorbing boundary conditions for the wave equation and parallel computing[END_REF]. In order to define the algorithm, we need a regularity result:

Proposition 4.1 For any ǫ, 0 < ǫ < 1, V 2 (Ω, T ) ⊂ C 0 (0, T ; H 2-ǫ (Ω)) ∩ C 1 (0, T ; H 1-ǫ (Ω)).
Proof By using extension operators in time and space, it suffices to prove the result in R × R. We make use of the Littlewood-Paley theory (see for example [START_REF] Chemin | Mécanique des fluides[END_REF]). In particular, there exists ϕ and χ two tempered distributions on R, with ϕ supported in (-8/3, -3/4) ∪ (3/4, 8/3), χ supported in (-4/3, 4/3), and

χ(ξ) + q≥0 ϕ(2 -q ξ) = 1, ∀ξ ∈ R.
We define the dyadic projectors ∆ q by their action on a function u,

∆ -1 u = χ(D)u, ∆ q u = ϕ(2 -q D)u for q ≥ 0, (18)
where D = -i∂. These operators give an equivalent norm in H s (R),

|u| s = q≥-1 2 2qs ∆ q u 2 1 2
.

They can also be used to define the Zygmund spaces

C r * = {u ∈ S ′ , ||| u||| r = sup q≥-1 2 qr ∆ q u 2 < +∞}.
C r * co ïncides with the usual Hölder space when r is not an integer. For any positive r, we know that W r,∞ , the space of functions in L ∞ with derivatives of order up to r in

L ∞ , is included in C r * . Therefore we have V 2 (R, R) ⊂ C 0 * (R, H 2 (R)) ∩ C 1 * (R, H 1 (R)) ∩ C 2 * (R, L 2 (R)
). We need an interpolation lemma. Lemma 4.2 For any positive α, β, a, b, for any θ, 0 ≤ θ ≤ 1,

C α * (R, H a (R)) ∩ C β * (R, H b (R)) ⊂ C θα+(1-θ)β * (R, H θa+(1-θ)b (R)).
Applying the lemma with successively (α, β, a, b) = (0, 1, 2, 1) and (α, β, a, b) = (1, 2, 1, 0), we find for any θ, θ ′ in (0, 1),

V 2 (R, R) ⊂ C 1-θ * (R, H 1+θ (R)) ∩ C 2-θ ′ * (R, H θ ′ (R)).
Since for any ǫ > 0 we have C ǫ * ⊂ C 0 and C 1+ǫ * ⊂ C 1 , this concludes the proof of Proposition 4.1.

Proof of Lemma 4.2 It relies on the convexity of the exponential function.

u 2 C θα+(1-θ)β * (R,H θa+(1-θ)b (R)) = sup j≥-1 2 2j(θα+(1-θ)β) k≥-1 2 2k(θa+(1-θ)b ∆ t j ∆ x k u 2
where ∆ t j (resp. ∆ x j ) is the Littlewood Paley operator acting in the time (resp. space) variable.

k≥-1

2 2k(θa+(1-θ)b ∆ t j ∆ x k u 2 = k≥-1 2 2ka ∆ t j ∆ x k u 2 θ 2 2kb ∆ t j ∆ x k u 2 (1-θ) ≤ k≥-1 2 2ka ∆ t j ∆ x k u 2 θ k≥-1 2 2kb ∆ t j ∆ x k u 2 1-θ .
Therefore we have

u 2 C θα+(1-θ)β * (R,H θa+(1-θ)b (R)) ≤ sup j≥-1 2 2jα k≥-1 2 2ka ∆ t j ∆ x k u 2 θ sup j≥-1 2 2jβ k≥-1 2 2kb ∆ t j ∆ x k u 2 1-θ which writes u 2 C θα+(1-θ)β * (R,H θa+(1-θ)b (R)) ≤ u 2 C α * (R,H a (R)) θ u 2 C β * (R,H b (R)) 1-θ . Theorem 4.3 Let p in H 2 0 (Ω) and q in H 1 0 (Ω).
There exists a time T 1 ≤ T * 0 such that for any T ≤ T 1 , for any initial guess h ± i in H 1 (0, T ) with the compatibility conditions

h + i (0) = q(a i+1 )+p ′ (a i+1 )+g + (p(a i+1 )) and h - i (0) = q(a i ) -p ′ (a i ) + g -(p(a i ))
, the algorithm ( 5) is defined and converges in ∪ i V 2 (Ω i , T ) to the solution ū of (4).

Proof We first prove that the algorithm is well-defined : with the assumptions on h ± i in the theorem, we know by Theorem 3.1 that [START_REF] Gander | Absorbing boundary conditions for the wave equation and parallel computing[END_REF]

defines in each Ω i a u 1 i in V 2 (Ω i , T ), with ∂ t u 1 i (a i , •), ∂ t u 1 i (a i+1 , •), ∂ x u 1 i (a i , •) and ∂ x u 1 i (a i+1 , •) in H 1 (0, T ) for T ≤ T k i . Furthermore, by Lemma 3.3, B -(g -)u 1 i-1 (a i , •) and B + (g + )u 1 i+1 (a i+1 , •) are in H 1 (0, T ).
As for the compatibility conditions, we have

B -(g -)u 1 i-1 (a i , 0) = lim t→0 ∂ t u 1 i-1 (a i , t) -∂ x u 1 i-1 (a i , t) + g -(u 1 i-1 (a i , t))
and by Proposition 4.1, we can pass to the limit and get

B -(g -)u 1 i-1 (a i , 0) = q(a i ) -p ′ (a i ) + g -(p(a i )).
This, together with the same regularity result on a i+1 , permits the recursion.

We define for T ≤ min(T * 0 , min i (T * k i )), for k ≥ 1, the quantities (with u i = ū/ Ωi ) for 1 ≤ i, j ≤ I, j = i or j = i -1,

e k i = u k i -u i , f k i = F (u k i ) -F(u i ) h k,- i,j = g -(u k j (a i )) -g -(u j (a i )), h k,+ i,j = g + (u k j (a i+1 )) -g + (u j (a i+1 )).
The operators F is defined in (16). The error e k i in Ω i at iteration k is a solution of

(∂ 2 t -∂ 2 x )e k i + e k i = f k i in Ω i × (0, T ) (19) (∂ t -∂ x )e k i + h k,- i,i = (∂ t -∂ x )e k-1 i-1 + h k-1,- i,i-1 on {a i } × (0, T ) (20)
(∂ t + ∂ x )e k i + h k,+ i,i = (∂ t + ∂ x )e k-1 i+1 + h k-1,+ i,i+1 on {a i+1 } × (0, T ) (21)
with vanishing initial values and e k 0 ≡ 0, e k I+2 ≡ 0, h k,- 1,0 ≡ 0, h k,+ I,I+1 = 0 . In order to get a new energy estimate in Ω i , we multiply (19) by ∂ t e k i and integrate by parts:

d dt E Ωi (e k i ) -[∂ t e k i ∂ x e k i (a i+1 , •) -∂ t e k i ∂ x e k i (a i , •)] = (f k i , ∂ t e k i ) (22) with E Ω (u) = 1 2 ( ∂ t u 2 L 2 (Ω) + ∂ x u 2 L 2 (Ω) + u 2 L 2 (Ω)
). We rewrite the boundary terms using the boundary operators:

∂ t e k i ∂ x e k i (a i+1 , •) = 1 4 ((∂ t + ∂ x )e k i (a i+1 , •) + h k,+ i,i ) 2 - 1 4 ((∂ t -∂ x )e k i (a i+1 , •) + h k,- i+1,i ) 2 + R k i,i+1 , -∂ t e k i ∂ x e k i (a i , •) = 1 4 ((∂ t -∂ x )e k i (a i , •) + h k,- i,i ) 2 - 1 4 ((∂ t + ∂ x )e k i (a i , •) + h k,+ i-1,i ) 2 + R k i,i-1 . (23)
The remainders R k i,i+1 and R k i,i-1 will be evaluated later. We insert (23) into ( 22), and obtain

d dt E Ωi (e k i ) + 1 4 ((∂ t -∂ x )e k i (a i+1 , •) + h k,- i+1,i ) 2 + 1 4 ((∂ t + ∂ x )e k i (a i , •) + h k,+ i-1,i ) 2 = 1 4 ((∂ t + ∂ x )e k i (a i+1 , •) + h k,+ i,i ) 2 + 1 4 ((∂ t -∂ x )e k i (a i , •) + h k,- i,i ) 2 + R k i,i+1 + R k i,i-1 + (f k i , ∂ t e k i ).
Using the transmission conditions (20), (21), we get

d dt E Ωi (e k i ) + 1 4 ((∂ t -∂ x )e k i (a i+1 , •) + h k,- i+1,i ) 2 + 1 4 ((∂ t + ∂ x )e k i (a i , •) + h k,+ i-1,i ) 2 = 1 4 ((∂ t + ∂ x )e k-1 i+1 (a i+1 , •) + h k-1,+ i,i+1 ) 2 + 1 4 ((∂ t -∂ x )e k-1 i-1 (a i , •) + h k-1,- i,i-1 ) 2 + R k i,i+1 + R k i,i-1 + (f k i , ∂ t e k i ). ( 24
)
We sum (24) on the indexes i, 1 ≤ i ≤ I and integrate in time. We translate the domain indexes in the right-hand side. Defining

E ∂Ωi (e k i ) = 1 4 [((∂ t -∂ x )e k i (a i+1 , •) + h k,- i+1,i ) 2 + ((∂ t + ∂ x )e k i (a i , •) + h k,+ i-1,i ) 2 ],
we get, since the initial data vanish,

I i=1 E Ωi (e k i )(t) + t 0 I i=1 E ∂Ωi (e k i )(s)ds ≤ t 0 I i=1 E ∂Ωi (e k-1 i )(s)ds + t 0 I i=1 (R k i,i+1 + R k i,i-1 )(s)ds + t 0 I i=1 (f k i , ∂ t e k i )(s)ds. (25)
Differentiating the equation and the transmission conditions in time yields the bound on ∂ t e k i :

I i=1 E Ωi (∂ t e k i )(t) + t 0 I i=1 E ∂Ωi (∂ t e k i )(s)ds ≤ t 0 I i=1 E ∂Ωi (∂ t e k-1 i )(s)ds + t 0 I i=1 ( Rk i,i+1 + Rk i,i-1 )(s)ds + t 0 I i=1 (∂ t f k+1 i , ∂ tt e k i )(s)ds. ( 26
)
We now estimate the remainders. We start with R k i,i+1 (ignoring the superscript k):

R i,i+1 = 1 4 [-h + i,i (2(∂ t e i + ∂ x e i )(a i+1 , •) + h + i,i ) + h - i+1,i (2(∂ t e i -∂ x e i )(a i+1 , •) + h - i+1,i )],
and we get a bound on the integral of R i,i+1 :

t 0 R i,i+1 (s)ds ≤ 3 4 ( h + i,i 2
H 1 2 (0,t) + h - i+1,i 2
H 1 2 (0,t) ) + 1 2 ( ∂ t e i (a i+1 , •) 2 H -1 2 (0,t) + ∂ x e i (a i+1 , •) 2 H -1 2 (0,t)
).

We can treat R i,i-1 , Ri,i-1 , and Ri,i-1 the same way and obtain

t 0 R i,i-1 (s)ds ≤ 3 4 ( h - i,i 2
H 1 2 (0,t) + h + i-1,i 2 H 1 2 (0,t) ) + 1 2 ( ∂ t e i (a i , •) 2 H -1 2 (0,t) + ∂ x e i (a i , •) 2 H -1 2 (0,t) ), t 0 Ri,i+1 (s)ds ≤ 3 4 ( ∂ t h + i,i 2
H 1 2 (0,t) + ∂ t h - i+1,i 2
H 1 2 (0,t) ) + 1 2 ( ∂ 2 t e i (a i+1 , •) 2 H -1 2 (0,t) + ∂ xt e i (a i+1 , •) 2 H -1 2 (0,t) ), t 0 Ri,i-1 (s)ds ≤ 3 4 ( ∂ t h - i,i 2
H 1 2 (0,t) + ∂ t h + i-1,i 2 H 1 2 (0,t) ) + 1 2 ( ∂ 2 t e i (a i , •) 2 H -1 2 (0,t) + ∂ xt e i (a i , •) 2 H -1 2 (0,t)
).

At point a i for instance, by the Trace Theorem, there is a constant C 3 independent of T , such that for any α with |α| = 1, we have

∂ α e i (a i , •) H -1 2 (0,t) ≤ ∂ α e i (a i , •) L 2 (0,t) ≤ C 3 ∂ α e i H 1 (Ωi×(0,t)) , ∂ α ∂ t e i (a i , •) H -1 2 (0,t) ≤ ∂ α e i (a i , •) H 1 2 (0,t) ≤ C 3 ∂ α e i H 1 (Ωi×(0,t)) ,
which gives our first bounds on the remainders:

t 0 R i,i+1 (s)ds ≤ 1 2 ( h + i,i 2
H 1 2 (0,t) + h - i+1,i 2
H 1 2 (0,t) ) + C 2 3 2 e i 2 H 2 (Ωi×(0,t)) , t 0 Ri,i+1 (s)ds ≤ 1 2 ( ∂ t h + i,i 2
H 1 2 (0,t) + ∂ t h - i+1,i 2
H 1 2 (0,t) ) + C 2 3 2 e i 2 H 2 (Ωi×(0,t)) .
We now insert the previous estimates in ( 25) and (26). By Cauchy-Schwarz inequality we get

I i=1 E Ωi (e k i )(t) + t 0 I i=1 E ∂Ωi (e k i )(s)ds ≤ t 0 I i=1 E ∂Ωi (e k-1 i )(s)ds + C 2 3 + 1 2 I i=1 e k i 2
H 2 (Ωi×(0,t)) + 1 2 I i=1 t 0 f k i 2 L 2 (Ωi) (s)ds + 3 4 I i=1 ( h k,+ i,i 2
H 1 2 (0,t) + h k,- i+1,i 2
H 1 2 (0,t) ) + 3 4 I i=1 ( h k,- i,i 2
H 1 2 (0,t) + h k,+ i-1,i 2 H 1 2 (0,t)
), ( 27)

I i=1 E Ωi (∂ t e k i )(t) + t 0 I i=1 E ∂Ωi (∂ t e k i )(s)ds ≤ t 0 I i=1 E ∂Ωi (∂ t e k-1 i )(s)ds + C 2 3 + 1 2 I i=1 e k i 2 H 2 (Ωi×(0,t)) + 1 2 I i=1 t 0 ∂ t f k i 2 L 2 (Ωi) (s)ds + 3 4 I i=1 ( ∂ t h k,+ i,i 2
H 1 2 (0,t) + ∂ t h k,- i+1,i 2
H 1 2 (0,t) ) + 3 4 I i=1 ( ∂ t h k,- i,i 2
H 1 2 (0,t) + ∂ t h k,+ i-1,i 2 H 1 2 (0,t) ) (28)
Adding ( 27) and (28), we can write

I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i ))(t) + t 0 I i=1 (E ∂Ωi (e k i ) + E ∂Ωi (∂ t e k i ))(s)ds ≤ t 0 I i=1 (E ∂Ωi (e k-1 i ) + E ∂Ωi (∂ t e k-1 i ))(s)ds + (C 2 3 + 1) I i=1 e k i 2 H 2 (Ωi×(0,t)) + 3 4 I i=1 ( h k,+ i,i 2
H 1 2 (0,t) + h k,+ i+1,i 2
H 1 2 (0,t) + h k,- i,i 2
H 1 2 (0,t) + h k,+ i-1,i 2
H 1 2 (0,t) ) + 3 4 I i=1 ( ∂ t h k,+ i,i 2 H 1 2 (0,t) + ∂ t h k,+ i+1,i 2 H 1 2 (0,t) + ∂ t h k,- i,i 2 H 1 2 (0,t) + ∂ t h k,+ i-1,i 2 H 1 2 (0,t) ) + I i=1 1 2 ( f k i 2 L 2 (Ω×(0,T )) + ∂ t f k i 2 L 2 (Ω×(0,T )) ).
We now estimate the quantities involving the h k,± i,j . Refining the results in Lemma 3.3, we have a real positive increasing function θ 2 , such that

h k,± i,j 2 H 1 2 (0,t) + ∂ t h k,± i,j 2 H 1 2 (0,t) ≤ θ 2 2 ( |α|≤2 α =(0,2) ( ∂ α u i 2 L 2 (Ω×(0,T )) + ∂ α u k i 2 L 2 (Ω×(0,T )) )) |α|≤2 α =(0,2) ∂ α e k i 2 L 2 (Ω×(0,T ))
which gives

I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i ))(t) + t 0 I i=1 (E ∂Ωi (e k i ) + E ∂Ωi (∂ t e k i ))(s)ds ≤ t 0 I i=1 (E ∂Ωi (e k-1 i ) + E ∂Ωi (∂ t e k-1 i ))(s)ds + I i=1 θ 3 ( |α|≤2 α =(0,2) ( ∂ α u i 2 L 2 (Ω×(0,T )) + ∂ α u k i 2 L 2 (Ω×(0,T )) )) e k i 2 H 2 (Ωi×(0,t)) + I i=1 1 2 ( f k i 2 L 2 (Ω×(0,T )) + ∂ t f k i 2 L 2 (Ω×(0,T )) ), (29)
with θ 3 = 3θ 2 2 + C 2 3 + 1. We now evaluate the terms in the right-hand side. We first note that

e k i 2 H 2 (Ωi×(0,t)) ≤ t 0 (E Ωi (e k i ) + E Ωi (∂ t e k i ) + ∂ xx e k i 2 L 2 (Ωi) )
∂ xx e k i 2 L 2 (Ωi×(0,t)) ≤ 3( ∂ tt e k i 2 L 2 (Ωi×(0,t)) + e k i 2 L 2 (Ωi×(0,t)) + f k i 2 L 2 (Ωi×(0,t)) ),
from which we deduce

e k i 2 H 2 (Ωi×(0,t)) ≤ 4 t 0 (E Ωi (e k i ) + E Ωi (∂ t e k i ))ds + 3 f k i 2 L 2 (Ωi×(0,t)) .
There remains only in (29)

I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i ))(t) + t 0 I i=1 (E ∂Ωi (e k i ) + E ∂Ωi (∂ t e k i ))(s)ds ≤ t 0 I i=1 (E ∂Ωi (e k-1 i ) + E ∂Ωi (∂ t e k-1 i ))(s)ds + I i=1 4θ 3 ( |α|≤2 α =(0,2) ( ∂ α u i 2 L 2 (Ω×(0,T )) + ∂ α u k i 2 L 2 (Ω×(0,T )) )) t 0 (E Ωi (e k i ) + E Ωi (∂ t e k i ))ds + I i=1 (3θ 3 ( |α|≤2 α =(0,2) ( ∂ α u i 2 L 2 (Ω×(0,T )) + ∂ α u k i 2 L 2 (Ω×(0,T )) )) + 1 2 )( f k i 2 L 2 (Ω×(0,T )) + ∂ t f k i 2 L 2 (Ω×(0,T )) ).
(30) Again, as in Lemma 3.3, there exists a positive increasing function θ 4 such that

f k i (t, •) 2 L 2 (Ω) + ∂ t f k i (t, •) 2 L 2 (Ω) ≤ θ 2 4 ( |α|≤2 α =(0,2) ( ∂ α u i 2 L 2 (Ω) + ∂ α u k i 2 L 2 (Ω) ) |α|≤2 α =(0,2) ∂ α e k i 2 L 2 (Ω) ,
the latter sum means that no term ∂ xx are present, therefore we can bound the sum by twice the energy. Furthermore we know that the energy of ū is bounded on the interval (0, T ). Thus there exists a new positive increasing function θ 5 , depending on u, such that

f k i (t, •) 2 L 2 (Ω) + ∂ t f k i (t, •) 2 L 2 (Ω) ≤ (θ 2 5 (E Ωi (e k i ) + E Ωi (∂ t e k i ))(E Ωi (e k i ) + E Ωi (∂ t e k i )))(t). ( 31
)
10

We insert (31) into (30), and get (with a new function θ 6 )

I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i ))(t) + t 0 I i=1 (E ∂Ωi (e k i ) + E ∂Ωi (∂ t e k i ))(s)ds ≤ t 0 I i=1 (E ∂Ωi (e k-1 i ) + E ∂Ωi (∂ t e k-1 i ))(s)ds + I i=1 t 0 (θ 6 (E Ωi (e k i ) + E Ωi (∂ t e k i )))(E Ωi (e k i ) + E Ωi (∂ t e k i ))(s)ds. (32)
For clarity we define

E k int = I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i )), E k b = I i=1 (E ∂Ωi (e k i ) + E ∂Ωi (∂ t e k i )).
and we can rewrite (32) as

E k int (t) + t 0 E k b (s)ds ≤ t 0 E k-1 b (s)ds + t 0 θ 6 (E k int (s))E k int (s)ds (33)
Summing in k, we define E K int = K k=1 E k int , and we have

E K int (t) + t 0 E K b (s)ds ≤ t 0 E 0 b (s)ds + t 0 θ 6 ( E K int (s)) E K int (s)ds.
Let now C > 0. If t tends to 0, t 0 E 1 b (s)ds + tCθ 6 (C) tends to 0. Therefore there exists a T 1 such that

T1 0 E 1 b (s)ds + T 1 Cθ 6 (C) = C
, and so we have for t ≤ T 1 ,

E K int (t) ≤ C.
We conclude that u k i exists on the time interval (0, T 1 ), and that I i=1 (E Ωi (e k i ) + E Ωi (∂ t e k i )) tends to 0 when k tends to infinity: the sequence u k i converges to ū on (0, T 1 ) in each subdomain in the norm of energy.

A Finite Volume Discretization

We use here a finite volumes scheme, which has been described in [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF] for the linear one-dimensional wave equation, and extended to the non linear boundary value problems in the frame of absorbing boundary conditions in [START_REF] Szeftel | Absorbing boundary conditions for non linear partial differential equations[END_REF]. We restrict ourselves to uniform meshes in time and space.

Discretization of the Subdomain Problem (6)

The domain Ω × (0, T ) is meshed by a rectangular grid, with uniform mesh sizes ∆x and ∆t. There are J + 1 points in space with ∆x = (a + -a -)/J, and N + 1 points in time, with ∆t = T /N . We denote the numerical approximation to u(a -+ j∆x, n∆t) by U (j, n). We introduce the notations:

D + t U (j, n) = U (j, n + 1) -U (j, n) ∆t , D - t U (j, n) = U (j, n) -U (j, n -1) ∆t , D + x U (j, n) = U (j + 1, n) -U (j, n) ∆x , D - x U (j, n) = U (j, n) -U (j -1, n) ∆x , D 0 x U (j, n) = U (j + 1, n) -U (j -1, n) 2∆x , D 0 t U (j, n) = U (j, n + 1) -U (j, n -1) 2∆t , D -- t U (j, n) = 3 U (j, n) -4U (j, n -1) + U (j, n -2) 2∆t , D - * t U (j, n) = D -- t U (j, n) for n ≥ 2, D - t U (j, n) for n = 1. ( 34
)
The last finite derivative in (34) is a second order approximation of ∂ t u, to be used in the nonlinear term, in order to design an explicit scheme. The scheme in the interior writes

D + t D - t -D + x D - x U (j, n) -f (U (j, n), D - * t U (j, n), D 0 x U (j, n)) = 0 in [ |1, J -1| ] × [ |1, N -1| ]. (35)
We define the discrete initial value as

P (j) = p(a -+ j∆x), Q(j) = q(a -+ j∆x),
and we obtain the initial scheme

(D + t - ∆t 2 D + x D - x )U (j, 0) = Q(j) + ∆t 2 f (P (j), Q(j), D 0 x Q(j)), in [ |1, J -1| ]. ( 36
)
For the boundary conditions, we define the discrete boundary operators as:

B -(f, g -) U (0, n) = (D 0 t -D + x + ∆x 2 D + t D - t )U (0, n) - ∆x 2 f (U (0, n), D - * t U (0, n), D + x U (0, n)) + g -(U (0, n)), B -(f, g -) U (0, 0) = (D + t -D + x + ∆x ∆t D + t ) U (0, 0) - ∆x ∆t Q(0) - ∆x 2 f (P (0), Q i (0), D + x P (0)) + g -(P (0)), (37)
B + (f, g + ) U (J, n) = (D 0 t + D - x + ∆x 2 D + t D - t )U (J, n) - ∆x 2 f (U (J, n), D - * t U i (J, n), D - x U (J, n)) + g + (U (J, n)). B + (f, g + ) U (J, 0) = (D + t + D - x + ∆x ∆t D + t )U (J, 0) - ∆x ∆t Q(J) - ∆x 2 f (P (J), Q i (J), D - x P (J)) + g + (P (J)). ( 38
)
We define the boundary data as

H ± (n) = 1 ∆t tn+∆t/2 tn-∆t/2 h ± (τ )dτ, 1 ≤ n ≤ N, H ± (0) = 2 ∆t ∆t/2 0 h ± (τ )dτ (39)
The discretization of Problem ( 6) is now given by ( 35), (36), with boundary conditions

B -(f, g -)U (0, •) = H -, B + (f, g + )U (J, •) = H + in [ |0, n ≤ N | ], ( 40
)
where the discrete boundary operators B ± are given in (37), (38).

Our numerical computations indicate that this scheme is second order both in space and time.

The Discrete Schwarz Waveform Relaxation Algorithm

The equation is now discretized on each subdomain Ω i × (0, T ), i = 1, . . . , I separately, using an uniform mesh with sizes ∆x and ∆t. There are J i +1 points in space and N +1 grid points in time in subdomain Ω i , with ∆x = (a i+1 -a i )/J i and ∆t = T /N . We denote the numerical approximation to u k i (a i + j∆x, n∆t) on Ω i at iteration step k by U k i (j, n). The problem in domain Ω i is now defined through boundary data H ±, k-1 i , coming from the neighboring subdomains Ω i±1 .Therefore, we define the extraction operator from domain Ω i to his neighbours as:

B + (f, g + ) U i (0, n) = (D 0 t + D + x - ∆x 2 D + t D - t )U i (0, n) + ∆x 2 f (U i (0, n), D - * t U i (0, n), D + x U i (0, n)) + g + (U i (0, n)), B + (f, g + ) U i (0, 0) = (D + t + D + x - ∆x ∆t D + t ) U i (0, 0) + ∆x ∆t Q i (0) + ∆x 2 f (P i (0), Q i (0), D + x P i (0)) + g + (P i (0)). (41) B -(f, g -)U i (J i , n) = (D 0 t -D - x - ∆x 2 D + t D - t )U i (J i , n) + ∆x 2 f (U i (J i , n), D - * t U i (J i , n), D - x U i (J i , n)) + g -(U i (J i , n)), B -(f, g -) U i (J i , 0) = (D + t -D - x - ∆x ∆t D + t ) U i (J i , 0) + ∆x ∆t Q i (J i ) + ∆x 2 f (P i (J i ), Q i (J i ), D - x P i (J i )) + g -(P i (J i )). ( 42
)
The discrete Schwarz waveform relaxation algorithm on subdomains Ω i , i = 1, . . . , I is defined as follows.

An initial guess {H 0,± i } 1≤i≤I is given. For k ≥ 1, we solve

       D + t D - t -D + x D - x U k i = f (U k i , D - * t U k i , D 0 x U k i ) in [ |1, J i -1| ] × [ |1, N | ], U k i (•, 0) = P i , (D + t - ∆t 2 D + x D - x )U k i (•, 0) = Q i + ∆t 2 f (P i , Q i , D 0 x Q i ) in [ |1, J i -1| ], B -(f, g -)U k i (0, •) = H -, k-1 i , B + (f, g + )U k i (J i , •) = H +, k-1 i in [ |0, n ≤ N | ] H -, k i = B -(f, g -)U k i-1 (J i-1 , •) , H +, k i = B + (f, g + )U k i+1 (0, •) in [ |0, n ≤ N | ]. (43)
As in the continuous algorithm, we set U k 0 ≡ 0 and U k I+2 ≡ 0. We denote by Ū the discrete approximation of problem (4), obtained by solving (35,36,40) on Ω = (a, b) with J = I i=1 J i intervals of length ∆x, and H ± = 0. Each subproblem is an explicit scheme, thus has a unique solution. Therefore the Schwarz waveform Relaxation Algorithm is well defined. If it converges, the limit in each subdomain is denoted by V i . It satisfies the same scheme as Ū at initial time, in the interior and on the exterior boundaries. At point a i it satisfies for any n ≥ 0

B -(f, g -)V i (0, n) = B -(f, g -)V i-1 (J i-1 , n) , B + (f, g + )V i-1 (J i-1 , n) = B + (f, g + )V i (0, n). (44)
Theorem 5.1 Suppose that f is affine in the third variable ∂ x u. Then, if the discrete algorithm converges, it converges to the discrete approximation Ū of problem (4).

Proof Since p is in H 2 (Ω) and q is in H 1 (Ω), there are both continuous and we have for any i V i (0, 0) =

P i (0) = V i-1 (J i-1 , 0) = P i-1 (J i-1 ) = p(a i ), and Q i-1 (J i-1 ) = Q i (0) = q(a i ).
We write the transmission conditions (44) for n = 0. The nonlinear terms containing g ± on both sides cancel out, and we have

(D + t -D + x + ∆x ∆t D + t )V i (0, 0) - ∆x ∆t q(a i ) - ∆x 2 f (p(a i ), q(a i ), D + x p(a i )) = (D + t -D - x - ∆x ∆t D + t )V i-1 (J i-1 , 0) + ∆x ∆t q(a i ) + ∆x 2 f (p(a i ), q(a i ), D - x p(a i )), (45)
(D + t + D + x - ∆x ∆t D + t )V i (0, 0) + ∆x ∆t q(a i ) + ∆x 2 f (p(a i ), q(a i ), D + x p(a i )) = (D + t + D - x + ∆x ∆t D + t )V i-1 (J i-1 , 0) - ∆x ∆t q(a i ) - ∆x 2 f (p(a i ), q(a i )), D - x p(a i )). (46)
Adding ( 45) and (46) yields

D + t V i (0, 0) = D + t V i-1 (J i-1 , 0), and hence V i (0, 1) = V i-1 (J i-1 , 1). We define now Ṽ (j, n) for 0 ≤ j ≤ J as Ṽ ( j, n) = V i (j, n) if j∆x = a i + j∆x, with 1 ≤ j ≤ J i . With the assumption on f , since D 0 x = (D + x + D - x )/2
, we can rewrite (45) as

(D + x -D - x -2 ∆x ∆t D + t ) Ṽ ( ji , 0) + 2 ∆x ∆t q(a i ) + ∆xf (p(a i ), q(a i ), D 0 x p(a i )) = 0,
with ji = J 1 + . . . + J i , which, multiplying by -∆t 2∆x , proves that Ṽ is solution of (36) at any point, and therefore U and Ṽ coincide at time 0 and 1. A simple recursion with the explicit schemes now proves that, for any n, for any i, V i (0, n) = V i-1 (J i-1 , n), and therefore U = Ṽ .

Remark 5.2

The assumption on f in Theorem 5.1 is fulfilled when

f (u, u t , u x ) = f 1 (u) + f 2 (u)u t + f 3 (u)u x .

Convergence of the Discrete Algorithm

According to Theorem 5.1, we suppose that f is affine in ∂ x u. We introduce the linear transmission operators defined by T ± = B ± (0, 0), T ± = B ± (0, 0)

We note Ūi = Ū / Ωi . With these notations, the error

Ū k i = U k i -Ūi is solution of the linear problem D + t D - t -D + x D - x + 1 Ū k i = F k i , in [ |1, J i -1| ] × [ |1, N | ] (47) (48)
with the initial value Ū k i (j, 1) = Ū k i (j, 0) = 0 and the transmission conditions

T -Ū k i (0, •) + G -, k i = T -Ū k-1 i-1 (J i-1 , •) + G -, k-1 i-1 , in [ |0, N | ], T + Ū k i (J i , •) + G +, k i = T + Ū k-1 i+1 (0, •) + G +, k-1 i+1 , in [ |0, N | ]. ( 49
)
The remainders are given by

F k i = f (U k i , D - * t U k i , D 0 x U k i ) -f ( Ūi , D - * t U k i , D 0 x Ūi ) + Ū k i for n ≥ 1, G -, k i = - ∆x 2 f k i (0, •) + g -(U k i (0, •)) -g( Ūi (0, •)), G -, k i-1 = ∆x 2 f k i-1 (J i-1 , •) + g -( Ū k i-1 (J i-1 , •)) -g -( Ūi-1 (J i-1 , •)), G +, k i = - ∆x 2 f k i (J i , •) + g + ( Ū k i (J i , •)) -g + ( Ūi (J i , •)), G +, k i+1 = ∆x 2 f k i+1 (0, •) + g + ( Ū k i+1 (0, •)) -g + ( Ūi+1 (0, •)), (50)
and

G ±, k i (0) = 0, G ±, k i (0) = 0.
For n = 0, the centered derivative in x are replaced in the expression of F k i by a forward or backward derivative.

We define now a discrete energy as follows. We consider sequences of the form V = {V (j)} 0≤j≤J in R J+1 , and we define a bilinear form on R J+1 by

a ∆ (V, W ) = ∆x( J j=1 D - x (V )(j) • D - x (W )(j) + J-1 j=1 V (j)W (j)). ( 51
)
For a mesh function V of time and space, we define

E K (V )(n) = ∆x 2 J-1 j=1 (D - t V (j, n)) 2 + (D - t V (j, n + 1)) 2 , E P (V )(n) = a ∆ (V (•, n), V (•, n -1)), E = E K + E P . (52)
The quantity E K is a discrete kinetic energy. It is less evident to identify E P as discrete potential energy.

The following lemma gives a lower bound for E under a CFL condition, and hence shows that E is then indeed an energy. The proof is classical ( [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF]) and is omitted here. Lemma 6.1 For any n ≥ 1, we have

E(V )(n) ≥ 1 - ∆t 2 ∆x 2 - ∆t 2 4 E K (V )(n). ( 53
)
Hence, under the CFL condition

∆t 2 ∆x 2 + ∆t 2 4 < 1, (54)
E is bounded from below by an energy.

The following energy estimate is obtained by a discrete integration by parts: Lemma 6.2 For any V solution of

D + t D - t -D + x D - x + 1 V = F, 1 ≤ j ≤ J -1, 1 ≤ n ≤ N, (55)
we have for any n ≥ 1,

E(V )(n) -E(V )(n -1) + ∆t 2 [( T + V (0, n)) 2 + ( T -V (J, n)) 2 ] = ∆t 2 [(T -V (0, n)) 2 + (T + V (J, n)) 2 ] + 2∆t∆x J j=1 F (j, n)D 0 t V (j, n), (56)
and for n = 0,

E K (U )(0) + E(U )(0) + ∆t 4 [( T + U (0, 0)) 2 + ( T -U (J + 1, 0)) 2 ] = ∆t 4 [(T -U (0, 0)) 2 + (T + U (J, 0)) 2 ] + a ∆ (P, P ) + 2∆x J j=1 Q(j)D + t U (j, 0). ( 57
)
These estimates are obtained by multiplying (55) with D 0 t V (j, n) and integrating by parts [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF]. We now state the main result of this section. Theorem 6.3 Suppose that f is affine with respect to ∂ x u. Defining the quantities

R k i (n) = G +, k i (n)( G +, k i (n) + 2 T + Ū k i (0, n)) + G -, k i (n)( G -, k 1 (n) + 2 T -Ū k i i , n)) -G -, k i (n)(G -, k i (n) + 2T -Ū k i (0, n)) -G +, k i (n)(G +, k i (n) + 2T + Ū k i (J i , n)), (58)
and assuming that there exists a positive constant M such that for any iteration number K, any domain Ω i and any discrete time n, the following estimate holds

2∆x J j=1 F k i (j, n)D 0 t Ū k i (j, n) + 1 2 R k i (n) ≤ M E( Ū k i )(n), (59)
then, for ∆t sufficiently small, the discrete Schwarz algorithm converges in the energy norm.

Proof We apply (56,57) to Ū k i . Since the initial data vanish, every term in (57) vanish. Thus E( Ū k i )(0) = 0, and we rewrite (56) as

E( Ū k i )(n) -E( Ū k i )(n -1) + ∆t 2 [( T + Ū k i (0, n) + G +, k i (n)) 2 + ( T -Ū k i (J i , n) + G -, k i (n)) 2 ] = ∆t 2 [(T -Ū k i (0, n) + G -, k i (n)) 2 + (T + Ū k i (J i , n) + G +, k i (n)) 2 ] + 2∆t∆x Ji j=1 f k i (j, n)D 0 t Ū k i (j, n) + ∆t 2 R k i (n). (60)
We now insert the transmission conditions (49), translate the indices in the righthand side, and add the contributions of all subdomains. We define a total internal energy and a total boundary energy as

E k I (n) = I i=1 E( Ū k i )(n), E k B (n) = ∆t 2 i [( T + Ū k i (0, n) + G +, k i (n)) 2 + ( T -Ū k i (J i , n) + G -, k i (n)) 2 ].
With these notations we can write

E k I (n) -E k I (n -1) + E k B (n) ≤ E k-1 B (n) + I i=1 [2∆t∆x Ji j=1 f k i (j, n)D 0 t Ū k i (j, n) + ∆t 2 R k i (n)]. (61)
We now sum up (61) for 1 ≤ k ≤ K, and define ÊK

I (n) = K k=1 E k I (n): ÊK I (n) -ÊK I (n -1) + E K B (n) ≤ E 0 B (n) + ∆t K k=1 I i=1 [2∆x Ji j=1 f k i (j, n)D 0 t Ū k i (j, n) + 1 2 R k i (n)]. (62)
Under assumption (59) we deduce that

ÊK I (n) -ÊK I (n -1) ≤ F 0 (n) + M ∆t ÊK I (n). ( 63
)
The recursive inequality is easy to solve. For ∆t sufficiently small, M ∆t < 1, and since ÊK (0) = 0, we get

ÊK I (n) ≤ e Mn∆t n p=1 F 0 (p) ≤ e MT N p=1 F 0 (p). ( 64
)
This proves that sup n≤N ÊK I (n) is bounded as K tends to infinity. Therefore we have ∀n ≤ N,

I i=1 E(U k i -Ū )(n) → 0 as k → +∞ (65)
which concludes the proof.

We are able to prove the assumption (59) in the case of the linear transmission conditions. Corollary 6.4 Suppose that f is affine with respect to ∂ x u. For ∆t sufficiently small, there exists a time T such that the discrete Schwarz waveform relaxation algorithm (43) with linear transmission conditions (i.e. g ± = 0) converges to the discrete approximation Ū of problem ( 4) on (0, T ).

Proof Here the remainder R k i (n) reduces to

R k i (n) = 2∆x(-f k i (0, n)D 0 t Ū k i (0, n) + f k i (J i , n)D 0 t Ū k i (J i , n)), (66)
and the estimate (59) amounts to proving that for any k ≤ K, n ≤ N , ∆x Ji j=0

f k i (j, n)D 0 t Ū k i (j, n) ≤ M 1 E( Ū k i )(n). ( 67
)
If f is globally Lipschitz in both variables, this is merely an application of the discrete Cauchy-Schwarz lemma.

7 Numerical Results

Remarks on Overlapping versus Nonoverlapping Schwarz Waveform Relaxation Algorithms and Variants

The original Schwarz algorithm uses overlapping domains (domains

Ω i = (a i , b i ) for 1 ≤ i ≤ I, with a i < a i+1 < b i < b i+1 )
, with an exchange of Dirichlet data on the boundary. It is known for elliptic problems to converge, but the smaller the overlap, the slower the convergence [START_REF] Lions | On the Schwarz alternating method[END_REF]. Due to the finite speed of propagation, it converges in a finite number of iterations, given by N = ⌈cT /L⌉ where c is the wave speed, and L the size of the overlap (see [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF] for the linear wave equation). Therefore the convergence can be very slow. In the linear case, using absorbing boundary conditions instead of the Dirichlet transmission, even without overlap, improves drastically the convergence, giving in one dimension a number of iterations equal to 2 for some T [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF]. In the nonlinear case, we have theoretical convergence results on the linear algorithm for sufficiently small T . However, we will see that the latter performs very well on a

The Non linear Nonoverlapping Schwarz Algorithms

We will see in this section that our strategy greatly improves the performances of the Classical Schwarz algorithm. Note that, whereas the classical algorithm only converges in the presence of an overlap, and the smaller the overlap, the slowlier the convergence, our algorithms are run without overlap, Here the residual is given by:

res k = (G +, k+1 -G +, k )(a 1 , •) 2 L 2 + (G -, k+1 -G -, k 2 )(a 1 , •) 2 L 2 1/2
.

We start with f = u 3 . In this case the transmission operators are the same and linear, since g ± = 0. We have proved in Corollary 6.4 that there exists a final time T for which the discrete algorithm is convergent.

In the forthcoming computations, the theoretical and numerical data are the same as before. Figure 2 plots the convergence history for various mesh sizes. The computation is stopped as soon as the residual reaches 0.5 10 -7 . We see that the algorithm converges very rapidly, independently of the mesh size.

We consider now the case f = u 2 u x , and nonlinear transmission conditions, i.e. g ± (u) = ±u 3 /6. In Figure 3, we plot the convergence history for various mesh sizes. The computation is stopped as soon as the residual reaches 0.5 10 -7 . Both the linear and the nonlinear transmission conditions behaves very well. The convergence takes place in 5 iterations with the linear transmission condition, in 4 iterations with the nonlinear one.

We also can vary the nonlinearities in the transmission conditions, we use a real parameter δ, and g ± = ±δu 3 . The nonlinear strategy corresponds to δ = 1/6, whereas the linear one is obtained for δ = 0. We draw in Figure 4 the error curves after 3 iterations for the same initial values as before, the mesh sizes are ∆x = 1/100, ∆t = 1/120. We observe that the nonlinear strategy corresponds precisely to the optimal numerical value of the parameter δ, validating the high frequency approach. We have carried out the same computations in the case f = u 2 u t . We do not display the results here since they are very similar.

Conclusion

We have presented a linear and a nonlinear Schwarz waveform relaxation algorithm without overlap for the semilinear wave equation. On the continuous level, we proved the convergence for sufficiently small time intervals. We designed a discrete algorithm, in such a way that, if convergent, the algorithm converges to the discrete solution in the whole domain, which we prove when the nonlinearity is affine in ∂ x u. In that case we proved the convergence for the linear transmission condition. Numerical experiments highlight the fast convergence to the discrete full domain solution in a large time domain in both linear

Figure 1 :

 1 Figure 1: Variation of the residual res k (left) and the error E k (right) as a function of the iteration number k.

Figure 2 :

 2 Figure 2: f = u 3 . Left residual, right: error. * : ∆x = 1/100, ∆t = 1/120, o: ∆x = 1/200, ∆t = 1/240, ♦: ∆x = 1/200, ∆t = 1/240.

Figure 3 :

 3 Figure 3: f = u 2 u x . Top: residual, bottom: error. Solid: linear, dash: nonlinear.

Figure 4 :

 4 Figure 4: f = u 2 u x . Variations of the error as a function of the nonlinearity coefficient δ

Acknowledgments

The authors are very grateful to Pr. Martin Gander from Genève University whose matlab scripts were a basis for the present work.

The second author is partially supported by NSF Grant DMS-0504720

large time interval, and that the nonlinear algorithm performs even better.

Our experiments concern the space domain Ω = [0, 4], simulating R with linear absorbing boundary conditions at each boundary. The time interval is [0,[START_REF] Chemin | Mécanique des fluides[END_REF]. Ω is divided in two subdomains. The initial value is supported in (0, 2), with p(x) = x 3 (2 -x) 3 , the initial velocity is q(x) = 3x 2 (2 -x) 2 (x -1). This is a good test since the solution is supported in the first subdomain at t = 0 and escapes in the second domain before the end of the computation.

Note that the Schwarz algorithm can be viewed as a fixed point algorithm applied to the interface problems

). In all cases, the stopping criterion in the algorithm will be on the residual res k for A. We also compute the exact discrete solution in Ω, and measure the discrete global error E k in L ∞ (0, T, L 2 (Ω)).

The Classical Overlapping Schwarz Algorithm

In this case, Ω 1 = (0, 2 + L) and Ω 2 = (2, 4). The stopping criterion pertains to the residual for the interface problem:

.

We run the computation until the residual res k is equal to zero. The theoretical minimal number of iterations for the discrete algorithm is

We start with the nonlinear term f (u, u t , u x ) = u 3 . Table 1 gives the number of iterations N comp needed to achieve convergence (i.e the error is zero), together with N th . On the left we choose the overlap equal to 8 grid points, and vary ∆t = 1/120 and ∆x = 1/100 to fulfill the CFL condition. On the right we fix ∆t = 1/120 and ∆x = 1/100, and vary the overlap L. We show on Figure 1 the convergence history of the classical Schwarz algorithm for various values of the mesh size and an overlap equal to eight gridpoints.We check in each case that the error E k vanishes together with the residual. Furthermore we can see that the error decays very slowly for many iterations, and reaches zero in a few iterations, independently of the mesh size (three or four in all cases). The behaviour is very similar to what happens for the linear wave equation : only the finite speed of propagation produces convergence, which takes place when the signal has left the domain. The algorithm behaves similarly for other nonlinearities. and nonlinear strategies, without overlap. Furthermore, we have shown that our nonlinear transmission conditions give optimal results within a large class of transmission conditions.


Michel Bergmann

Angelo Iollo

Bioinspired Swimming Simulations

Keywords: Bioinspired Swimming, Level Set, Skeleton, Optimal transportation, Cartesian mesh, Penalization Nage bio-inspirée, fonctions level set, squelettisation, transport optimal, maillage Cartesien, penalisation

We present a method to simulate the flow past bioinspired swimmers starting from pictures of an actual fish. The overall approach requires i) a skeleton graph generation to get a level-set function from pictures; ii) optimal transportation to obtain the velocity on the body surface; iii) flow simulations realized thanks to a cartesian method based on penalization. We illustrate this paradigm by simulating the swimming of a mackerel fish.

Introduction

Autonomous swimming robots can address in the future needs such as search and rescue, environmental monitoring and emergency response in rapidly unfolding scenarios. Fishes have the hability to evolve in complex and unpredictable environments that may abruptly change. They are capable of maneuvering in ways that can hardly be achieved by present engineering devices. In this sense, the design of such devices can significantly benefit from bioinspired principles.

However, quantitative evaluation of the geometrical and dynamical characteristics of swimming in experiments is challenging [START_REF] Triantafyllou | Hydrodynamics of flishlike swimming[END_REF][START_REF] Wen | Understanding undulatory locomotion in fishes using an inertiacompen sated flapping foil robotic device[END_REF]. In particular, precise control of the kinematics and dynamics of the experimental model is difficult; harder still is the accurate measurement of forces and power of a freely swimming model. In this regard, numerical models are well suited for such investigation.

Three-dimensional numerical simulation of swimming has only become viable in the last decade [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF][START_REF] Anup | A new mathematical formulation and fast algorithm for fully resolved simulation of self-propulsion[END_REF][START_REF] Curet | A versatile implicit iterative approach for fully resolved simulation of self-propulsion[END_REF][START_REF] Hieber | An immersed boundary method for smoothed particle hydrodynamics of self-propelled swimmers[END_REF][START_REF] Gazzola | Simulations of single and multiple swimmers with non-divergence free deforming geometries[END_REF][START_REF] Bergmann | Effect of caudal fin flexibility on the propulsive efficiency of a fish-like swimmer[END_REF] thanks to the progress of both numerical schemes and computers architectures. Nonetheless, many of the studies presented in the literature typically deal with prescribed swimmer geometries. Only few start from actual fish geometries, like for example Patankar's group [START_REF] Anup | A new mathematical formulation and fast algorithm for fully resolved simulation of self-propulsion[END_REF][START_REF] Curet | A versatile implicit iterative approach for fully resolved simulation of self-propulsion[END_REF][START_REF] Bhalla | A unified mathematical framework and an adaptive numerical method for fluid-structure interaction with rigid, deforming, and elastic bodies[END_REF]] and Mittal's group [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF][START_REF] Bozkurttas | Low-dimensional models and performance scaling of a highly deformable fish pectoral fin[END_REF][START_REF] Dong | Computational modelling and analysis of the hydrodynamics of a highly deformable fish pectoral fin[END_REF] .

The aim of this paper is to model a three-dimensional swimmer starting from actual fish pictures. The first step is to build the three-dimensional fish profile starting from two-dimensional data retrieved from one or several subsequent pictures of a fish at rest. To this end, a skeleton technique [START_REF] Galusinski | Skeleton and level set for channel construction and flow simulation[END_REF] is used to reconstruct the three-dimensional level set function describing the body surface. Then the skeleton can be deformed using an appropriate swimming law to obtain a sequence of level set functions corresponding to snapshots of the body surface uniformly taken at different instants. Otherwise, if a sequence of pictures of swimming is available, the body surface can be reconstructed with the same technique for each snapshot. The level set functions give at each instant an implicit description of the surface which is quite different from usual swimming simulations where the deforming body is fully described by in a Lagrangian framework.

Thanks to skeleton deformation or to a sequence of actual pictures we typically can reconstruct only a small portion of the snapshots necessary to simulate a swimming stroke, since the time scale of the simulation is significantly smaller than the time step between two subsequent reconstructed three-dimensional snapshots. Moreover, the surface deformation velocity is required to set the boundary conditions of the flow problem. For this reason it is necessary to build intermediate level set functions and to compute the deformation velocity field between subsequent reconstructed fish snapshots. Optimal transportation [START_REF] Villani | Optimal Transport, old and new[END_REF] is well suited to achieve this goal providing an objective model to compute intermediate geometries and deformation velocities.

Optimal transportation is an old topic first introduced by Monge in 1781 [START_REF] Monge | Memoire sur la theorie des des deblais et des remblais[END_REF]. It consists in a distance minimization problem that allows to define a mapping between an initial and a final mass distribution. This problem has attracted considerable attention these last years especially from a numerical point of view [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF][START_REF] Iollo | A lagrangian scheme for the solution of the optimal mass transfer problem[END_REF][START_REF] Saumier | An efficient numerical algorithm for the l2 optimal transport problem with periodic densities[END_REF][START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF][START_REF] Hug | Multi-physics optimal transportation and image interpolation[END_REF]. Here we exploit the scheme introduced in [START_REF] Bouharguane | Numerical solution of the monge-kantorovich problem by density lift-up continuation[END_REF] to define a non-linear interpolation mapping between geometries at different instants and to compute the corresponding deformation velocities.

In the following we describe the numerical scheme, the swimmer kinematic modeling and show validations with examples of applications.

Ω f (χ = 0) Ω s (χ = 1) Γ s u s u, p ∂Ω

Numerical modeling

Flow problem

The flow configuration is given in figure 1 where Ω s is the domain defined by the body, Ω f is the domain filled by the fluid, and Ω = Ω f ∪ Ω s is the whole domain. The boundary between Ω f and Ω s is noted Γ s .

The flow is modeled using the incompressible Navier-Stokes equations written in the fluid domain Ω f as follows

ρ ∂u ∂t + (u • ∇)u = -∇p + ∇ • 2µD(u) + ρg in Ω f , (1a)
∇ • u = 0 in Ω f . (1b)
where

D(u) = ∇u+∇ T u 2
, u is the velocity field, p is the pressure field, µ is the dynamic viscosity, ρ is the density and g is the gravity. Unsteady boundary conditions are imposed on the body boundary

u(x, t) = u s (x, t) on Γ s , (1c)
and on the external boundary

u(x, t) = u ext (x, t) on ∂Ω. ( 1d
)
The dynamics of the boundary Γ s that represents the swimmer will be described later.

Classically, system (1) is discretized using a body fitted mesh such that the boundary condition on the body boundary can be easily and accurately taken into account. However, these kind of meshes are difficult to generate in complex unsteady geometries in three-dimensions as large mesh deformations are required when the body is displaced or is changing topology.

Inria

In this study we will use a fixed cartesian grid. The body boundary is defined and tracked thanks to a level set function [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] and the boundary condition on Γ s is taken in account in the Navier-Stokes equation using a penalty term [START_REF] Angot | A penalization method to take into account obstacles in a incompressible flow[END_REF].

We thus consider the penalized Navier-Stokes equations in the whole domain Ω (fluid and body):

∂u ∂t + (u • ∇)u = - 1 ρ ∇p + 1 ρ ∇ • 2µD(u) + g + χλ(u s -u) in Ω, (2a)
∇ • u = 0 in Ω, (2b)
where boundary conditions on the body boundary are implicitly imposed through the penalty term χλ(u su). In the penalty term χ is the characteristic function (χ(x) = 1 if x ∈ Ω s and χ(x) = 0 elsewhere) and λ is a large penalty parameter (we chose λ = 10 8 ). In what follows, the characteristic fonction will be computed from the signed distance fonction ϕ, so that χ = H(ϕ) where H denotes the Heaviside function. A second order penalization method can easily be obtained thanks to the signed distance function [START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF].

It has been proven [START_REF] Van Yen | Approximation of the laplace and stokes operators with dirichlet boundary conditions through volume penalization: A spectral viewpoint[END_REF] that the solution of the system (2) tends to the solution of the system (1) as 1/ √ λ. The system (2) is solved in time using a projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des equations de navier-stokes par la methode des pas fractionnaires ii[END_REF] and is discretized in space on a uniform cartesian mesh using second order centered finite differences for the diffusive term and third-order upwind finite differences for the convective term.

The resolution of system (2) requires knowledge of both the unsteady characteristic function χ and the body velocity u s . The following sections will explain how to compute these quantities.

For a self propelled body the body velocity can be decomposed such that

u s (x, t) = u(x, t) + u θ (x, t) + u(x, t) (3)
where u(x, t), u(x, t) and u θ (x, t) are the deformation, the translation and the rotation velocities, respectively. While the deformation velocity is imposed (by muscles), the translation and rotation velocities are the results of the hydrodynamic forces and torques applied by the fluid onto the body. These latter are computed integrating the Newtons laws. Once the body geometry is obtained (or imposed), we can deform the backbone (the midline) following swimming laws observed in nature. Considering that the backbone of the steady body is 0 ≤ x ≤ ℓ, y = 0 and z = 0, we can apply a swimming law y(x, t) that remains the backbone length to be constant. The swimming law observed in nature for most of fishes is [START_REF] Barrett | Drag reduction in fish-like locomotion[END_REF]:

y(x, t) = a(x) sin(2π(x/λ + f t)), (4a)
where the envelop is

a(x) = c 0 + c 1 x + c 2 x 2 . ( 4b
)
This swimming law is defined with a constant phase velocity c p = λf , where λ and f are the wave length and the frequency of the oscillations. The parameters c 0 , c 1 and c 2 are adjusted to reach a maximal tail amplitude A/2 that is an important parameter for the swimming [START_REF] Lighthill | Aquatic animal propulsion of high hydrodynamical efficiency[END_REF]. It has been shown [START_REF] Videler | Fast continuous swimming of two pelagic predators, saithe (pollachius virens) and mackerel (scomber scombrus): a kinematic analysis[END_REF][START_REF] Yu | Three-dimensional numerical simulation of hydrodynamic interactions between pectoral-fin vortices and body undulation in a swimming fish[END_REF] that carangiform swimming mode is close to c 0 = 0.02, c 1 = -0.12 and c 2 = 0.2 for a unit length fish with 0 ≤ x ≤ 1. These parameters can be rescaled when modifying the fish length.

Generation of level set functions from pictures: the skeleton approach

Thanks to the scheme presented above it was possible to perform numerical simulations of fishlike swimmers [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF][START_REF] Bergmann | Effect of caudal fin flexibility on the propulsive efficiency of a fish-like swimmer[END_REF] or jellyfishes [START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF]. The swimmer geometries were prescribed to mimic realistic fishes or jellyfishes. The motion of these swimmers was described in a purely Lagrangian way by tracking markers distributed on the swimmer boundary. The objective is now to simulate bioinspired swimmers using geometries and possibly deformation laws obtained by pictures or sequences of pictures. As an example we consider the definition of a swimmer starting from the mackerel fish presented in figure 2 (left column) [START_REF] Hovnanian | Methode de Frontieres Immergees pour la Mecanique des Fluides[END_REF].

As mentioned in §2, our numerical approach requires the knowledge of the level set function (distance function) associated with a body to get a second order scheme in space. A characteristic function (mask function) is sufficient for a first order penalty method [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF].

The process is synthesized in figure 2. We start with pictures of a mackerel fish. Top, bottom and side views are necessary. To rebuild fish in 3D from only three pictures (front, side, top), an underlying assumption is made: the whole cross-section of the fish along length is geometrically similar to the fish on front view (front picture). A dilatation takes place in width and another takes place in height. In the tail region of the fish, two connected components can appear on the cross section, each connected component will be geometrically similar to the front view. The scaling factors will be extracted from the thickness of the fish (from the side view) and the width (from the top view), according to the longitudinal position of the transverse cut. Segmentations are first performed to extract the boundaries (black lines in center column). The distance function (right column) can be computed by solving Eikonal equation. Then the skeleton is computed (red lines in the center column) as being the center points of the inscribed circles. The natural tool to identify this information is the centerline or skeleton of an object as recalled in [START_REF] Galusinski | Skeleton and level set for channel construction and flow simulation[END_REF]. It identifies the locus of points equidistant from further points of the boundary. The special feature here is that the relevant skeleton is not constructed from the Euclidean distance. In order to identify the thicknesses and the widths of the fish, the distance is measured in a fixed direction (depending on the height or along the width) of the side view or top view. The skeleton Inria is then simplified compared to the skeleton associated with the Euclidean distance. Thus, the top view will have a skeleton made of a line from the head to the tail of the fish, with a single point of the skeleton on each transverse cut. The skeleton of the profile view for the vertical distance will consist of a central line splitting in the fish's tail. Finally, for each transverse cut of the profile, the number of connected components of the profile view is the number (1 or 2) of points of the skeleton of this cut. This property is true thanks to the choice of the distance and is not true with Euclidian distance. pictures has been performed so that the fish length is the same for top and side views, and height is the same for top and side views. A simple 3D reconstruction of the fish could then take place, but the goal is to deform the fish, which will have the effect of losing the notion of directional distance, more complex work is to be implemented. First, the skeleton of the top view (for the vertical distance) S v is distorted by law ( ref eq: swimlaw). the concept of directional distance is no longer suitable for reconstruction, then the chosen approach is to distort the classic skeleton associated with Euclidean distance, relying on distorted S v . This new skeleton allows to reconstruct, for the Euclidean distance, the deformed fish seen from top view. Then, the 3D reconstruction of the fish can be implemented by harnessing both concepts skeleton as shown in figure 4 (right).

It is thus possible to generate a series of level set functions corresponding to different time t for the swimming law [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]. In what follows, we considered the periodical swimming law (4) and uniformly extract N snapshots for the level set function over one period (one swimming stroke). Finally, the level set functions have been modified (transported) with a rigid velocity so that we do not add any linear or angular translations during deformation. Linear and angular translations have only to be the results of the hydrodynamic forces and torques.

Note that in what follows, the level set are chosen to be the sign distance functions (positive inside the body).

Computation of intermediate level-sets and deformation velocities: optimal transportation

As mentioned, we have access to N (here N = 20) snapshots {ϕ k } k=0,...,N -1 of the level set function over a swimming period T . Here, k denotes the time t = k T i/N for the i th swimming period. The number of available snapshots is of course inessential, but it is typically small compared to the number of snapshots necessary for a numerical simulation. It is then necessary to build missing level set function corresponding to the Naiver-Stokes time discretization. Moreover, it is also necessary to compute the fish velocity for the penalty method. This velocity can be decomposed into a deformation velocity, to be computed from the subsequent N level set functions, and a rigid velocity computed from hydrodynamics effects. The deformation velocity is derived from an optimal transportation model.

Optimal transportation is a model to compute a mapping between two densities. The model assumes that the mass is transported along straight lines with constant initial velocity. The initial velocity is found such that to overall distance measured in a given norm is minimum. The formulation of the L 2 optimal transportation problem is illustrated with respect to generic densities ̺ 0 and ̺ 1 . The main step is to compute the velocity field transforming ̺ 0 into ̺ 1 . Then it is possible to compute any desired intermediate densities for time t = r T i/N with 0 < r < 1.

Let ̺ 0 (ξ), ̺ 1 (x) be two smooth enough non-negative scalar density functions with compact support Ω 0 and Ω 1 , where ξ, x ∈ R d and d is the space dimension. We assume that

Ω0 ̺ 0 (ξ)dξ = Ω1 ̺ 1 (x)dx.
Let X : Ω 0 → Ω 1 a smooth one-to-one map such that X(ξ) realizes the transfer of ̺ 0 onto ̺ 1 , i.e., a map that satisfies the following Jacobian equation:

̺ 0 (ξ) = det(∇X(ξ))̺ 1 (X(ξ)).
This equation is under determined with respect to X(ξ) and a solution is selected among all possible maps by introducing the following L 2 Kantorovich-Wasserstein distance:

inf X Ω0 ̺ 0 (ξ)|X(ξ) -ξ| 2 dξ.
The L 2 Monge-Kantorovich problem (MKP) corresponds to finding a map X o such that this infimum is achieved. It has been proved that this problem admits a unique solution [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal Transport, old and new[END_REF], which is the gradient of a.e. convex function Ψ o : Ω 0 → R:

X o (ξ) = ∇Ψ o (ξ).
The main idea of solution relies on the assumption that we dispose of a mapping at iteration n that is a perturbation of the optimal mapping X o (ξ) [START_REF] Bouharguane | Numerical solution of the monge-kantorovich problem by density lift-up continuation[END_REF]. We derive a linear perturbation equation that is used to iteratively improve the initial guess. Let us assume that the mapping obtained at iteration n is

X n (ξ) = ∇ ξ Ψ o + ∇ ξ Ψ n ǫ .
Inria where the error Ψ n ǫ satisfies Ψ n ǫ 2 ≈ ǫ. We define ̺ n 0 (ξ) as the initial density at iteration n that mapped by X n (ξ) gives the exact final density ̺ 1 (x). Then, taking a first-order Taylor expansion, we have

̺ n 0 (ξ) := ̺ 1 (X n (ξ)) det (∇ ξ X n (ξ)) = ̺ 1 (X o ) det (∇ ξ X o ) =̺0(ξ) + det (∇ ξ X o )̺ 1 (X o )Tr (∇ ξ X o ) -1 ∇ 2 ξ Ψ n ǫ + det (∇ ξ X o )∇ x ̺ 1 (X o ) • ∇ ξ Ψ n ǫ + o(ǫ),
where Tr denotes the matrix trace operator. At first order in ǫ, we have

̺ n 0 (ξ) -̺ 0 (ξ) det (∇ ξ X o ) ≈ ̺ 1 (X o )Tr (∇ ξ X o ) -1 ∇ 2 ξ Ψ n ǫ + ∇ x ̺ 1 (X o ) • ∇ ξ Ψ n ǫ (5)
The mapping update is then found by the iteration

X n+1 = X n -α ∇Ψ n ǫ ,
that converges to X o as a geometric series, for α ∈ [0, 1], in the linearized approximation [START_REF] Loeper | Numerical solution of the monge-ampere equation by a newton's algorithm[END_REF][START_REF] Saumier | An efficient numerical algorithm for the l2 optimal transport problem with periodic densities[END_REF].

If the optimal solution and the approximated mapping are a perturbation of identity, then equation ( 5) reduces to

̺ n 0 (ξ) -̺ 0 (ξ) ≈ ∇ ξ • (̺ 1 (X n (ξ))∇ ξ Ψ n ǫ ) (6)
that is equivalent to a semidiscretization in time of the continuity equation written in Eulerian form.

Choice of the density functions

Since the level set functions corresponding to the N fish snapshots are not strictly positive, they are not admissible candidates for optimal transportation. A regularized mask function (Heaviside of the level set function) is also not satisfying since the weight of the thin caudal fin in the mapping is very limited and the mapping tends to smear it in comparison to the whole body.

We thus considered Gaussian density functions defined from the level sets {ϕ k } k=0,...,N -1 , with a support that is approximatively equal to ℓ/10 (ℓ being the length of the fish) with maximum near the interface but inside the body. In particular, we take

̺ 0 = G(ξ, k) = exp (ϕ k (ξ, t)-ϕ b ) 2 γ
+ε, with γ = 2•10 -5 , ϕ b = 1.5•10 -3 and the regularization parameter ε = 10 -3 . In a similar way, we take ̺ 1 = G(ξ, k + 1). Using these definitions of the densities to be transported, each point of the body surface has a comparable weight in the transportation. We denote by G 0 the value of the Gaussian function corresponding to the swimmer surface. The density G(ξ, k) is presented in figure 5 where the external white line is the body boundary G 0 .

Computation of intermediate images

Starting from the optimal mapping, the computation of intermediate images is straightforward. Let X(ξ, τ ) = ξ + τ (X o -ξ) and We have ̺(ξ, 0) = ̺ 0 (ξ), ̺(ξ, 1) = ̺ 1 (X o (ξ)) and for 0 < τ < 1 we obtain the intermediate images. As for the deformation velocity, given the time step ∆t s between two images (here ∆t s = T /20), the deformation velocity from ̺ 0 to ̺ 1 is

̺(ξ, τ ) = ̺ 0 (ξ) det (∇ ξ X(ξ, τ ))
u 0 = X o /∆t s (7)
The index 0 highlights the fact that this velocity has to be mapped to the actual position of the swimmer in space, since the level set snapshots are taken with respect to a center of mass fixed in space. Figure 6 shows the reconstruction the intermediate density for τ = 0.5 between densities ̺ 0 and ̺ 1 .

The same approach can be applied to the whole series of level sets {ϕ i } N i=0 , thus computing the corresponding deformation velocities and intermediate densities. Starting from these densities we can compute a characteristic function and a sign distance function thanks to the fast marching approach [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] or by solving an Hamilton Jacobi equation [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Luddens | Enablers for high-order level set methods in fluid mechanics[END_REF]. In what follows, ϕ 0 (ξ, t) are the level set function of the deformed profiles with fixed center of mass and u 0 (ξ, t) is the deformation velocity. Since we are interested in a self propelled swimmer, the position of the swimmer is determined thanks to the level set function ϕ s (ξ, t) that is an image of the level set function ϕ 0 (ξ, t) obtained by the transformation characterizing the rigid motion (translation plus rotation), X r (ξ, t). We have ϕ s (ξ, t) = ϕ 0 (X r (ξ, t)).

The characteristic function χ required in the penalty method is

χ(ξ, t) = H(ϕ s (ξ, t) -G 0 ) see §4.
1 for the definition of the constant G 0 . We can then impose χ(ξ, t) = 1 inside the whole body, and compute the signed distance function for second order penalty method [START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF].

In the same way, the velocity of the swimmer is u 0 (ξ, t) = u(ξ, t) + u θ (ξ, t) + u 0 (ξ, t), where u 0 (ξ, t) is computed from [START_REF] Bergmann | Effect of caudal fin flexibility on the propulsive efficiency of a fish-like swimmer[END_REF]. Since the position of the swimmer is modified, the velocity has to act onto the modified position and we have:

u s (ξ, t) = u 0 (X r (ξ, t)).
The transformation X r (ξ, t) as well as the translation velocity u(ξ, t) and rotation velocity u θ (ξ, t) are evaluated using the Newton's laws computing the forces and the torques exerted by the fluid onto the body.

Computation of the forces and the torques

The forces and the torques exerted by the fluid onto the body could be computed integrating the stress tensor T(u, p) = -pI + µ(∇u + ∇u T ) on the body boundary. Let n be the unit outward normal vector to the interface Γ, and r = xx G with x G the mass center, the force and the torques write:

F = - Γs T(u, p) n dx, (8a)
M = - Γs r ∧ T(u, p) n dx. (8b)
However, since the body boundary is defined in an implicit way as being the zero level of the sign distance function, integration of equations ( 8) is thus not straightforward. Thus we turn to moment balance as follows.

We consider an arbitrarily subdomain Ω f delimited by the boundaries Γ s and Γ f surrounding the body with domain Ω s such that the forces and torques write [START_REF] Noca | On the evaluation of time-dependent fluid-dynamic forces on bluff bodies[END_REF]:

F = - d dt Ω f (t) u dV + Γ f (t) (T + (u -u f ) ⊗ u) n dS - Γs(t) ((u -u s ) ⊗ u) n dS. (9a) M = - d dt Ω f (t) r i ∧ u dV + Γ f (t) r ∧ (T + (u -u f ) ⊗ u) n dS - Γs(t) r ∧ ((u -u s ) ⊗ u) n dS. (9b)
RR n°8874

where u s and u f are respectively the velocities of boundaries Γ s (thus, the body velocity) and Γ f . Note that in many applications (like for swimmers), the integral over the Γ s boundary vanish due to the fact that we have u = u s .

To simplify the computation, it is possible to chose Ω f = Ω f and thus Γ f = ∂Ω or smaller regular boxes that fit the mesh grid.

6 Validation and Numerical results

Validation

The numerical solver has been already intensively validated for several bodies where the geometry is defined by markers and the body followed in a Lagrangian way [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF][START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF][START_REF] Bergmann | Effect of caudal fin flexibility on the propulsive efficiency of a fish-like swimmer[END_REF]. Our goal is then to validate the proposed method using optimal transportation comparing it with a well defined test case. All the following numerical simulations are performed using cartesian mesh with ∆x = ∆y = ∆z = ℓ/150 on 256 CPUs.

In this section we do not consider an actual sequence of fish snapshots. We consider a fish model where the geometry is known and defined by markers. The geometry is however chosen to fit as best as possible the mackerel fish, except the lunate tail. Figure 7 presents the fish geometry for this test case as well as the markers (note that we actually use two times more markers in each directions). Our goal is then to compare the two numerical models for the fish just defined:

• Case I ("Lagrangian"): the level set functions and the deformation velocities are given at each Navier-Stokes time step by following the Lagrangian markers on the body surface;

• Case II ("Eulerian"): we suppose that we only know 20 snapshots of the level set function over one stroke period. We have thus to compute all the missing level set as well as the deformation velocity using optimal transportation.

In both cases the length fish is ℓ = 10 cm and the swimming law (4) is defined with frequency f = 4 Hz and c 0 = -0.002, c 1 = -0.12, c 2 = 2 and λ = ℓ. These parameters are close to Inria those observed in nature [START_REF] Yu | Three-dimensional numerical simulation of hydrodynamic interactions between pectoral-fin vortices and body undulation in a swimming fish[END_REF]. A comparison of the deformation velocity fields obtained in cases I ("exact") and II (approximated using optimal transportation) is presented in figure 8. The swimming law is defined in the plan (X, Z), and we present the velocities U and W where the velocity U is the velocity in the swimming direction. A very good agreement is obtained for the velocity W : the Eulerian deformation velocity field and the Lagrangian one continuously match at the body surface. Somme difference are however visible for the U component. Due to the fish configuration, the velocity U can be considered to be the velocity that is tangential to the body. Since optimal transportation is irrotational by definition, it is not able to model any skin rotation. The tangential velocity can not be computed if no other constraints are taken into account, such as done for example in [START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF][START_REF] Hug | Multi-physics optimal transportation and image interpolation[END_REF] using a rigid motion constraint or a divergence-free constraint.

(a) Velocity in the swimming direction (b) Velocity normal to the swimming direction Figure 8: comparison of the deformation velocity fields obtained in cases I ("exact", presented inside the fish) and II (approximated using optimal transportation, represented outside the fish).

Theses Lagrangian and the Eulerian velocities are then used to simulate the self propelled swimming in the both cases. These self propelled velocities are presented in figure 9. As expected, the swimming velocity W shows good agreement between cases I and II. The velocity U is influenced by the differences observed for the deformation velocity U in figure 8. However, considering all the modeling hypothesis, the Eulerian forward swimming velocity correctly approaches the actual Lagrangian forward velocity from which it has been derived.

Mackerel fish swimming from pictures

We have applied the method described above to recover the geometry of a Mackerel, to deform its skeleton by a swimming law and to simulate swimming. Figure 10 shows the wake flow generating by the fish built from pictures. We chose the same parameters as before for the swimming law, i.e. the fish length is10 cm and the swimming frequency is f = 4 Hz. The swimming velocity of the mackerel fish is slightly higher than the previous test cases (u ≈ -0.15 m/s). This difference can be explained by a different shape, and moreover by the lunate tail for the real fish.

Conclusions

In this work we have described a procedure to define a bioinspired geometry for a swimmer starting either from one image of a fish or from a sequence of images. Thanks to this procedure it is possible to recover an Eulerian deformation velocity on the base of a limited number of modeling hypothesis, namely those at the base of L 2 optimal transportation. We have shown an in silico example where we compare the swimming velocity of a swimmer whose geometry and deformation law is completely determined a priori, to that obtained by the identification method we propose. The results in terms of deformation velocity show that the irrotational hypothesis at the base of optimal transportation induces an error in the tangential velocity of the body. This error coupled to the discretization inaccuracies due to the fully Eulerian treatment of the geometry obtained from Lagrangian data reflects in an overestimation of the actual velocity of the swimmer of less than 10%. As an example of an actual application we have modeled a

Figure 1 :

 1 Figure 1: Sketch of the general flow configuration.

Figure 2 :

 2 Figure 2: Skeleton steps: pictures (left), segmentation with skeleton (center) and associated level set functions (right).

  Figure 3: Re-centering (top view).

Figure 4 :

 4 Figure 4: Three-dimensional profiles. Undeformed (left) and deformed (right).

Figure 5 :

 5 Figure 5: Gaussian density used in the plan z = 0. The external white contour corresponds to the fish boundary.

Figure 6 :

 6 Figure 6: Construction of an image for τ = 0.5∆t s (red) between ̺ 0 (green) and ̺ 1 (blue) in the plan z = 0.

Inria 5 Dynamics of the swimmer 5 . 1

 551 Computation of the swimmer velocity and position

Figure 7 :

 7 Figure 7: Geometry profile used to validate the method.

Figure 9 :

 9 Figure 9: Comparison of the self propelled velocities obtained in cases I and II.

Figure 10 :

 10 Figure 10: Swimming simulation of the mackerel fish starting form pictures.

RR n°8874

Inria

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

The authors thank C. Galusinski and C. Nguyen for contributions on fish construction. This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the ?Investments for the future? Programme IdEx Bordeaux-CPU (ANR-10-IDEX-03-02). The simulations presented in this paper were carried out using the PLAFRIM experimental parallel testbed, being developed un-der the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d?Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/).