Bernard Di Martino
email: dimartin@univ-corse.fr

Catherine Giacomoni

Jean-Martin Martin Paoli

Pierre Simonnet

Simulation of a viscous fluid spreading by a bidimensional shallow water model

Keywords: Shallow water equations, Free boundary problem, ALE discretization, Flooding and drying AMS Subject Classification : 35Q30, 65M60, 76D03

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The flooding and drying of a fluid on ground is a problem which has several applications in fluid mechanics such as coastal engineering, artificial lake filling, river overflow or sea intrusion on ground due to a tsunami. These works involves the modelling of the physical process near the triple contact line between the liquid (the water), the gas (the atmosphere) and the solid (the ground). The numerical simulation of this problem is very complex, due to the domain deformation and the triple contact line movement. Some numerical approaches to solve this problem have been previously explored. For example, O. Bokhove [START_REF] Bokhove | Flooding and drying in discontinuous galerkin finite-element discretizations as shallow water equations. Part 1: One dimension[END_REF] uses a conservative form of the shallow water equations associated to the discontinuous galerkin discretization. D. Yuan et al. [START_REF] Yuan | Simulating moving boundary using a linked groundwater and surface water flow model[END_REF] propose to use a non conservative form of the shallow water equations in total flow rate formulation and uses a finite difference scheme for the numerical computation.

In this paper, we use a two dimensional viscous shallow water model in eight velocity formulation. This model, called Stokes-like by P.L. Lions ([6], section 8.3) is, in some sense, intermediate between the semi-stationary model and the full model of compressible isentropic Navier-Stokes model. This model can be obtained by integrating vertically the so-called primitives equations of the ocean [START_REF] Lions | On the equation of the large scale ocean[END_REF] with some hypothesis on the viscous terms. We have chosen this model because we can prove the existence of a solution under certain considerations as the smoothness of the initial conditions and an acceptable hypothesis on a boundary operator [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF].

After the presentation of the model we describe the numerical method based on the demonstration of the convergence proposed in [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF]. Finally, we present some numerical examples in two idealized domains and we give some perspectives for other studies.

Mathematical model

We assume that we know a continuous function H from R 2 to R that represents the topography (the level of the ground with respect to a reference level).

We note Ω t (an open simply connected set of R 2 ) the horizontal domain occupied by the fluid at time t and η(t, x, y) the elevation of the fluid compared to a horizontal zero level. We set h(t, x, y) the eight of the water column, h(t, x, y) = H(x, y) + η(t, x, y).

The shallow water equations are based on a depth integration of an incompressible fluid conservation laws in a free surface-three dimensional domain. Governing equations for u and h can be obtained in the usual way ( [START_REF] Bosseur | Resolution by Galerkin method with a special basis of a geophysical flow in open sea : a Calvi's bay simulation[END_REF] for example) and the two-dimensional system can be written as follows [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF][START_REF] Martino | Peybernes Simulation of a spilled oil slick with a shallow water model with free boundary[END_REF][START_REF] Flori | Orenga Convergence of a Characteristic-Galerkin scheme for a shallow water problem[END_REF]:

∂u ∂t + (u • ∇)u + g∇η -ν∆u + C d |u|u = f 1 , in t∈[0,T ] {t} × Ω t (1)
∂η ∂t + div(uh) = 0, in t∈[0,T ] {t} × Ω t ( 2
)
where g is the gravity constant, ν the viscosity coefficient, C d a drag coefficient. f 1 represents external forcing (for example the wind stress). These equations are completed by initial conditions (u(t = 0) = u 0 , η(t = 0) = η 0 ) and boundary conditions. In the following, we assume that the boundary of the real domain is non-vertical. Hence the variation of the surface elevation at the boundary imposes a modification of the horizontal domain occupied by the fluid.

1.2 Boundary conditions 1.2.1 In a three dimensional domain

In many geophysical applications, the governing equations are solved in a domain [0, T ] × Ω 3 t whereby Ω 3 t = {(x, y, z) ∈ R 3 /(x, y) ∈ Ω t , -H(x, y) < z < η(t, x, y)} is a three dimensional moving domain. We denote by N = (n 1 , n 2 , n 3 ) the exterior unit normal to this domain, and we set N h = (n 1 , n 2 ) as its horizontal component. Moreover we denote V = (v 1 , v 2 , v 3 ) as the velocity of the fluid and v h = (v 1 , v 2 ) as its horizontal component.

In the following, we assume that the ground of the sea does not have a vertical wall. Then we have n 3 = 0. The impermeability boundary condition on the ground for the three dimensional problem imposes V (x, y, -H(x, y)) • N = 0. Then, we have v h • N h = -v 3 n 3 on the ground and, since n 3 = 0, we obtain We recall that we are interested in the effect of the overflow and we need to consider the problem in a moving domain that represents the domain actually occupied by the fluid. Then we need to impose a condition on the moving boundary. At the triple contact line (between air, ground and water, see figure 1), the horizontal velocity at the surface (v h (x, y, η, t)) is equal to the velocity at the bottom (v h (x, y, -H, t)). Similarly, we have v 3 (x, y, η, t) = v 3 (x, y, -H, t), coherent with the fact that the lagrangian derivative of h is equal to zero on this boundary.

v 3 = -v h • ∇H.
This property is not verified if, during its displacement, a part of the triple contact line touches a vertical wall. In this case, the boundary cannot be represented by the line h = 0. This question is very interesting but not treated in this paper.

In a two dimensional domain

The mean velocity u used in the equations 1 and 2 is related to the previous three dimensional velocity by u

= (u 1 , u 2 ) = 1 h η -H (v 1 , v 2
)dz where h = 0. The two dimensional domain Ω t is the projection of the three dimensional domain Ω 3 t on the horizontal plane.

For the depth integrated model, we denote by γ 0 the boundary of Ω 0 (at initial time). Assuming that γ 0 is smooth enough, we define the deformed boundary as follows :

γ t = {x ′ = x + d(x, t), x ∈ γ 0 }, where d represents the horizontal displacement d(x, t) = χ(0, t, x)-
x and where χ(s, t, x) denotes the Lagrangian flow, i.e. the position at time t of the particle located at x at time s. With these notations, we have the natural condition ∂d(χ(0, t, x), t)/∂t = U(x, t) on γ t , where U(x, t) = u(x + d(x, t), t).

The continuity equation ( 2) can be written as follows:

h(χ(x 0 , 0, t), t) = h(x 0 , 0) exp - t 0 div u . (3)
So, if at the initial time, h(x 0 , 0) = 0 then at all times h(χ(x 0 , 0, t), t) = 0. Moreover, if h(x 0 , 0) = 0 then h(χ(x 0 , 0, t), 0) = 0 at all times. We conclude that the boundary (and only the points of this boundary) at the initial time can define the boundary at all times. For example, if at initial time the domain is simply connected, the domain remains simply connected at all times t 0.

We characterize the boundary motion γ t by a condition on the normal component of the fluid stress tensor σ on the boundary

σ(x + d(x, t), t) • n(x + d(x, t), t)|detJ|(x, t) = A ∂U(x, t) ∂t , (4)
where A is a differential operator defined on γ 0 taking into account the boundary forcing on the fluid (see [START_REF] Beale | Large-time Regularity of Viscous Surface Waves[END_REF] for example) and J is the jacobian matrix associated to the transformation x → χ(x, 0, t). In [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF], the authors show that if the operator A is a Laplace-Beltrami operator (

γ 0 A(v)v = γ 0 A 1/2 (v)A 1/2 (v) = ||v|| 2 H 2 (γ 0 )
) and if the data are small enough, the problem has a solution. The condition on A is necessary to obtain the smoothness of the boundary and allows the use of the classical Sobolev injections [START_REF] Solonikov | Unsteady motion of a finite mass of fluid bounded by a free surface[END_REF][START_REF] Solonikov | The solvability of the problem of the motion of a viscous incompressible fluid that is bounded by a free surface[END_REF]. More precisely, in this paper, D(A) = H p (γ 0 ) with p = 2 and currently, we are not able to prove an existence result for p < 2. But numerically, it is possible to obtain physical acceptable results with less smoothness on this operator (for example in [START_REF] Martino | Peybernes Simulation of a spilled oil slick with a shallow water model with free boundary[END_REF] we have obtained results by taking A ≡ 0).

It is not easy to characterize the operator A. It needs a good physical interpretation and give sufficient mathematical smoothness to prove the existence of a solution of our problem. Physically, this operator represents a friction condition in the triple contact line (between the fluid, the solid domain and the air). In the literature, a large number of physical approaches in the case of a wetting film of fluid or for a droplet [START_REF] De Gennes | Wetting: Statics and Dynamics[END_REF] are found. Unfortunately, we cannot use directly these results. Firstly, we do not work at the microscopic scale (mainly governed by the Van Der Waals forces) and we need to take into account the effect of these phenomena at large scale. Secondly, the main effect studied for those phenomenon is in the vertical direction (for example for the capillary effect) but here, we use a depth integrated model.

2 The numerical method

The ALE equations

We use a numerical method based on the weak formulation of the problem in the domain Ω t depending on time. The ALE method that here we use, has been used to solve the Navier-Stokes equations in a moving domain [START_REF] Maury | Characteristics ALE method for the Unsteady 3D Navier Stokes Equations with a free surface[END_REF].

We assume that at the initial time t 0 , the fluid domain is covered by a regular mesh (for example, an usual finite element mesh). We also assume that on the boundary, each point of this mesh has the velocity of the fluid. The velocity of the internal points can be chosen arbitrarily, the only condition is to conserve the smoothness of the mesh at all times.

Then, at any time t ∈ [0, T ], the velocity c t of the moving mesh can be defined by solving the following problem:

∆c t = 0 in Ω t (5) c t = u(x, t) on γ t ( 6
)
and we consider the following mapping

c : Q → R 2 (x, t) → c(x, t) = c t (x),
where

Q = ∪ t∈[0,T ] (Ω t × {t}).
The relation between the different domains Ω t i is given by the mapping

C(•, t 1 , t 2 ) : Ω t 1 → Ω t 2 x 1 → x 2 = C(x 1 , t 1 , t 2 )
where C(., t 1 , t 2 ) is the characteristic curve from (x 1 , t 1 ) to (x 2 , t 2 ) corresponding to the velocity c t in the space-time domain Q. For each time τ , we denote by u τ and h τ the ALE velocity and the ALE thickness:

u τ (x, t) = u(C(x, τ, t), t), h τ (x, t) = h(C(x, τ, t), t). (7)
With these notations, we can write the ALE formulation of the problem

∂u τ ∂t + (u τ -c τ )∇u τ -µ∆u τ + C d u τ |u τ | + g∇h τ = O(t -τ ), (8)
∂h τ ∂t + (u τ -c τ )∇h τ + h τ div u τ = O(t -τ ). (9)
In what follows, we will use a first order time discretization to solve this problem (see [START_REF] Fourestey | Une Méthode des Caractéristiques d'ordre deux sur Maillages Mobiles pour la Résolution des équations de Navier-Stokes incompressibles par éléments finis[END_REF] for more details about second order schemes).

Time discretization

Let M 1, we note ∆t = T M as the time step. For 1 n M, 1 m M, we set Ω n = Ω t n with boundary γ n and

g n m (x) = g t m (x, t n ) for x ∈ Ω m , g n (x) = g n n (x) = g(x, t n ) for x ∈ Ω n .
In order to ensure the positivity of h t n , the continuity equation is renormalised as follows

∂ ∂t log h t n + (u t n -c t n )∇ log h t n + div u t n = 0. ( 10
)
With this renormalisation, we do not have conservation of the mass, but when we follow the evolution of the mass during a simulation (e.g. in the first test presented in the following section), only small variations of the mass quantity are observed.

We denote by U n (respectively H n and C n ) the approximation of the exact solution u t n (respectively h t n and c t n ).

Then, we note Ũn (x) = U n (X n (x, t n )) and Hn (x) = H n (X n (x, t n )) with X n (x, •) the characteristic curve solution of:

     ∂X n ∂t (x, t) = (U n -C n )(X n (x, t)) X n (x, t n+1 ) = x.
In our study, we approximate the foot of the characteristic by

X n (x, t n ) ≃ x-(U n -C n )(x)∆t. We set U n+1 n (resp. H n+1 n
) the approximation of u t n (x, t n+1 ) (resp. h t n (x, t n+1 ). Approximating the Lagrangian derivative by a first order Euler scheme, and using a linearised drag operator, we obtain the following:

U n+1 n + ∆t g∇H n+1 n + C d U n+1 n | Ũn | -ν ∆U n+1 n = Ũn , in Ω n (11) log H n+1 n + ∆t div U n+1 n = log Hn in Ω n , ( 12
) div U n+1 n = A ∂ Ũn ∂t on γ n , (13) curl
U n+1 n = H n+1 n = 0 on γ n . (14)
Implicitly, taking into account the equation 3, we have H n+1 n = 0 on γ n . This is an approximation at the first order in time of the initial shallow water problem.

Remark 2.1

This problem is solved by a fixed point technic. A first approximation of 12 is computed assuming U n+1 n = Ũn and the approximation of H n+1 n is then used in 11 to obtain a first approximation of U n+1 n . We repeat this operation in order to obtain convergence of the system. Prove of convergence can be fount in [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF].

We cannot easily increase the order of this scheme because the ALE formulation is only of the first order.

The previous problem is solved on Ω n by using the spatial discretization proposed in the following section. After, we need to solve the problem given by equations 5 and 6 in order to compute the mesh velocity. Each point p k with coordinates x k of the mesh is then moved with the first order approximation

x n+1 k = x n k + ∆tC n (x n k , t n ).

Spatial discretization

The spatial discretization of the previous problem is based on the Finite Element Galerkin method. In [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF], an approach based on the Galerkin method with a special basis is proposed, but this approach cannot be applied easily if A ≡ 0.

We note {φ k } the set of finite element functions of H 1 (Ω n ). The weak formulation of our approximated problem is

Ω n U n+1 n (x)φ k (x) + ν ∆t Ω n ∇U n+1 n (x)∇φ k (x) -g ∆t Ω n H n+1 n (x)div φ k (x) + C d ∆t Ω n U n+1 n (x)|U n n (x)|φ k (x) + γt σ(X + d(X, t), t) • n(X + d(X, t), t)γ(φ k (x)) = Ω n Ũn (x)φ k (x) (15) log H n+1 n (x) + ∆t div U n+1 n (x) = log Hn (x) in Ω n (16)
where the operator γ is the trace operator from Ω n to γ n .

Then the equations ( 4) and [START_REF] Solonikov | The solvability of the problem of the motion of a viscous incompressible fluid that is bounded by a free surface[END_REF] give

Ω n U n+1 n (x)φ k (x) + ν ∆t Ω n ∇U n+1 n (x)∇φ k (x) -g ∆t Ω n H n+1 n (x)div φ k (x) + C d Ω n U n+1 n (x)|U n n (x)|φ k (x) + γ 0 A 1/2 ∂ Ũ ∂t A 1/2 γ( φk (x 0 )) dγ 0 = Ω n Ũn (x)φ k (x) (17)
where g(x) = g(χ(x, 0, t)).

We then use a first order discretization of the partial time derivative:

γ 0 A 1/2 ∂ Ũ ∂t A 1/2 γ( φk (x 0 )) dγ 0 ≃ 1 ∆t γ 0 A 1/2 Ũn+1 n (x) A 1/2 γ( φk (x 0 )) dγ 0 - 1 ∆t γ 0 A 1/2 Ũn (x) A 1/2 γ( φk (x 0 )) dγ 0 . ( 18
)
The first part of the right hand of this equation needs to be included on the left hand of the global problem and the second part on the right hand.

The global weak formulation is:

Ω n U n+1 n (x)φ k (x) + ν ∆t Ω n ∇U n+1 n (x)∇φ k (x) -g ∆t Ω n H n+1 n (x)div φ k (x) Ω n U n+1 n (x)|U n n (x)|φ k (x) + 1 ∆t γ 0 A 1/2 Ũn+1 n (x) A 1/2 γ( φk (x 0 )) dγ 0 = Ω n Ũn (x)φ k (x) + 1 ∆t γ 0 A 1/2 Ũn (x) A 1/2 γ( φk (x 0 )) dγ 0 , ( 19
) log H n+1 n (x) + ∆t div U n+1 n (x) = log Hn (x). ( 20
) Remark 2.2 On γ t , U(l, t) = k a k φ k (l, t)
where l is a point of γ t and the sum is computed on all finite element functions. Then, on the initial boundary γ 0 , we have the decomposition Ũ(l 0 , t) = k a k φk (l 0 , t) with the same a k .

Remark 2.3 Formally, we can use two kinds of boundary conditions:

(1) A condition of normal displacement of the boundary. We write u • n = ∂d ∂t where d is the normal displacement. With this condition, the well posed weak formulation associated to the diffusion operator is Ω div udiv φ

+ Ω curl ucurl φ -γ div uφ • n + γ curl uφ • α(n) where α(u 1 , u 2 ) = (-u 2 , u 1 )
and curl u = ∂u 2 ∂y -∂u 1 ∂x . The boundary condition is then imposed on hdiv u and curl u since we do not want to impose a condition on u • n or u • α(n).

(2) A condition on the displacement of the boundary u = ∂d ∂t where d is a vector. With this condition, the well posed weak formulation associated to the diffusion operator is Ω ∇u∇φγ ∇uφ•n and the boundary condition is then imposed on ∇u (more precisely, we need to take into account the trace tensor including condition on h in the boundary, but we assume h ≡ 0 on this boundary).

Remark 2. [START_REF] De Gennes | Wetting: Statics and Dynamics[END_REF] We can note that the friction term C d u|u| is necessary to stabilize the flow. Indeed, if we use the following decomposition

(L 2 (Ω n )) 2 = ∇H 1 0 (Ω n ) ⊕ Curl H 1 0 (Ω n ) ⊕ ∇H(Ω n ), ( 21
)
where H is the intersection of H 1 and the space of harmonic functions, we can see that the functions of ∇H(Ω n ) are not controlled by the laplacian operator. We present in the following section a case with C d = 0 where we do not have damping of the flow.

First numerical test

The previous method is tested in order to simulate the behavior of a fluid in a simplified domain. We assume that the domain is axisymmetric using the function H(r) : r → ar 2 + 1. At the initial time t = 0, the fluid at rest touches the wall (h(a, t = 0) = 0) (see figure 2). Physically, our initial domain Ω 0 is a disc with a radius of 130 meters. This domain is meshed with triangles (420 triangles for the mesh M1, 470 for the mesh M2 and 1002 for the mesh M3).

x a 0 H(x) η(x,t) (a) (b)
The height of the column of water at the centre is 1 meter (since r = 0 at the centre of our domain). Fluid viscosity is 0.01m 2 .s -1 and gravity coefficient g is assumed equal to 1m.s -2 in order to amplify the elevation.

We apply a forcing at the surface of this fluid. This forcing is usually taken into account by applying the continuity of the horizontal stress tensor on the surface. So, for the three dimensional model,

∂u h ∂z = CW a |W a |
where W a is the wind velocity (at ten meters over the flow for the ocean for example) and C is a "drag coefficient" assumed to be constant. Using the vertical integration of the vertical operator ∂ ∂z ν z ∂u h ∂z , and assuming u h (η) ≃ u, we obtain a forcing condition on the mean velocity f ∝ Wa|Wa| H . We use here

     f = 1 if 0 < t < 20 s f = 0 if t > 20 s . ( 22
)
Hence, we observe an oscillation of the free surface of the water plan. At each oscillation, a part of the water moves on ground and a part of the ground is uncovered by this water (see figure 3). In the last part of the remark 2.4, we indicate that with C d = 0 a part of the flow component is not diffusive. More precisely, our solution is only a gradient. If we do not take into account boundary friction effects or bottom friction effects, a part of the fluid is not diffusive. To observe this effect, we plot (on figure 4) the level of a boundary point according to time. We can observe an accumulation of energy on In figure 6, we plot the time evolution of the kinetic energy of the flow for the different meshes (M1, M2 and M3). We only have a little difference for all these meshes but a very expensive cost for the mesh M3. After the forcing phase (20 seconds), we can observe the transformation of the kinetic energy in potential energy and conversely. When the kinetic energy vanishes, potential energy is maximal. If the drag coefficient is assumed to be equal to zero, the amplitude of the oscillation does not decrease.

Finally, figure 7 represents the evolution of the mass of fluid compared to the initial mass. Even if we do not have a rigorous conservation of the mass, due to the renormalisation of the mass equation, the variation of the mass is very little.

Second numerical test

In this second numerical experiment, we use a more complex domain to test our numerical method. We assume that the domain is axisymmetric using the function H(r) : r → ar 3 + br 2 + cr + d. b, c and d are chosen in order to at the distance 1/a of the centre of the fluid domain, the fluid touches the wall (see figure 8). We use the same external forcing as for the first experiment (22). We present some results of the simulation in figure 9. Recalling that the shallow water equations are based in the continuum mechanic, it cannot cut this domain with a finite energy. Hence, even if the thickness of the layer of water is very low, we conserve a thin film of water between all the parts of the domain occupied by the fluid. We recall that the computation is only made in the moving two-dimensional domain Ω t . Mathematical study proves that the main difficulty is to conserve the smoothness of the boundary. This result can be observed in this simulation because, at the final time, we have contact between two parts of the boundary (see last part of the figure 9). Since we do not use a "good" boundary operator, we do not have sufficient smoothness on the boundary. We observe then a change of connectedness and a hole appears in the fluid domain.

Conclusion

The presented work suggests that capillarity effect needs to be incorporated into the evolution equations, and more precisely in the triple contact line. Theoretical results presented in [START_REF] Flori | An existence result for a free boiundary shallow water model using a agrangian Scheme[END_REF] give sufficient smoothness for the capillarity term used in equation 4.

The main difficulty is to describe the capillarity effects on the depth integrated model and in a dynamical system. A possible approach is to test some boundary operator and to compare numerical and experimental result. This approach will be proposed in future studies.

Fig. 1 .

 1 Fig. 1. Vertical section of the three dimensional domain

Fig. 2 .

 2 Fig. 2. Simplified domain -(a) at rest. (b) at an arbitrarily time.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Oscillation of the water

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Oscillations of the water with a drag term

Fig. 8 .

 8 Fig. 8. Mesh of the ground and initial position of the water But, due to the specific form of the topography, a part of the fluid goes to the external crown.

Fig. 9 .

 9 Fig. 9. Evolution of the fluid domain


M Manguoglu

M Sathe

O Schenk

Ahmed Sameh

Ahmed Sameh A Scalable

Parallel Preconditioned

Sparse Linear

A Scalable Parallel Preconditioned Sparse Linear System Solver

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

List of authors: M. Manguoglu 1 M. Sathe O. Schenk 2 A. Sameh 3

Achieving high parallel scalability of sparse linear system solvers implemented on large-scale computing platforms comprised of tens of thousands of multicore processors is a task that offers many challenges. Towards achieving such a solver, we extract an effective banded (sparse within the band) preconditioner for use with an outer Krylov subspace method. The extraction of such banded preconditioners is made possible via a fast parallel reordering scheme based on the computation of the Fiedler vector of the weighted Laplacian of the corresponding sparse matrix under consideration. Solving systems involving these banded preconditioners (narrowbanded, generalized-banded, or wide-banded) is accomplished via members of the SPIKE family of solvers. The resulting hybrid parallel sparse system solver, P-SPIKE (for Pardiso-Spike), is made quite scalable via a specialized version of the sparse direct solver Pardiso. We show that our solver is quite competitive in performance and robustness compared to current parallel direct solvers, as well as Krylov subspace methods preconditioned via: (i) approximate LUfactorization, and (ii) algebraic multigrid schemes.