Olivier Delestre
email: olivier.delestre@unice.fr

Frédéric Darboux

François James

Carine Lucas

Christian Laguerre

Stéphane Cordier

FullSWOF: A free software package for the simulation of shallow water flows

Keywords: Overland flow, modeling, shallow-water equations, software, open source

Numerical simulations of flows are required in numerous applications and are typically performed using shallow-water equations. Here, we describe the FullSWOF software, which is based on up-to-date finite volume methods and on well-balanced schemes used to solve these types of equations. The software consists of a set of open-source C++ codes and is freely available to the community, easy to use, and open for further development. Several features make FullSWOF particularly suitable for applications in hydrology: small water depths and wet-dry transitions are robustly addressed, rainfall and infiltration are incorporated, and data from grid-based digital topographies can be used directly. A detailed mathematical description is given here, and the capabilities of FullSWOF are illustrated based on analytic solutions and on datasets compiled from real cases. The codes, available in 1D and 2D versions, have been validated on a large set of benchmark cases, which are available together with the download information and documentation at

Introduction

The simulation of geophysical surface flows is required when analyzing a large variety of natural and man-made situations. For such purposes, sev-eral mathematical models are available. More complex models utilize the full system of Navier-Stokes equations; these equations can be solved by Gerris, for example [START_REF] Popinet | Quadtree-adaptative tsunami modelling[END_REF]. However, the parametrization of the model is very intricate, and people typically prefer simplified models that are chosen based on the characteristics of the simulation. For avalanche modeling (on large slopes), a simplification of the Navier-Stokes equations leads to the Savage-Hutter system [START_REF] Savage | The dynamics of avalanches of granular materials from initiation to runout. part i: Analysis[END_REF]. In the case of shallow flows, if the vertical motion of the fluid can be neglected, the shallow-water equations can be used. Additional simplifications may be needed. For example, if the shallow-water equations continue to be too complex, they can be replaced by the diffusive wave or kinematic wave equations [START_REF] Moussa | Approximation zones of the Saint-Venant equations for flood routing with overbank flow[END_REF][START_REF] Novak | Hydraulic modelling -An Introduction[END_REF]. In this hierarchy of mathematical models, the present paper considers the shallow-water equations and presents the FullSWOF program. The shallow-water equations were first formulated by Saint-Venant (Barré de [START_REF] De Saint-Venant | Théorie du mouvement nonpermanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF] in the study of floods and tides. These equations are now classical equations in the fields of flood forecasting, pollutant transport, dam breaks, tsunamis, soil erosion by overland flow, etc. Because no explicit solutions to the shallow-water equations are known for general situations, efficient and robust numerical simulations are required.

Over the last forty years, numerous codes that make use of various methods have been developed in this field. The MacCormack scheme has been widely used for scientific purposes (e.g., [START_REF] Zhang | Modeling of two-dimensional overland flow[END_REF][START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Fiedler | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF]. Although it is relatively easy to program and performs computations quickly, this scheme neither guarantees the positivity of water depths at the wet-dry transitions nor preserves steady states, i.e., it is not well balanced [START_REF] Lee | Simple and efficient solution of the shallow water equations with source terms[END_REF], which requires some work if these issues are to be addressed (e.g., [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Fiedler | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF]. In the various industrial codes used in engineering (e.g., CANOE [START_REF] Tanguy | [END_REF], HEC-RAS [START_REF] Brunner | HEC-RAS River Analysis System[END_REF], ISIS [START_REF] Halcrow | ISIS 2D quick start guide[END_REF]MIKE11 (DHI Software, 2009)), the flow equations, namely, the shallow-water equations, are often solved in a nonconservative form [START_REF] Novak | Hydraulic modelling -An Introduction[END_REF] with either the Preissmann scheme or the Abbott-Ionescu scheme, leading to inaccurate calculations for transcritical flows and hydraulic jumps.

While all these pieces of software have been used for scientific research, most of their source codes have not been made available to the community, raising a major concern about research reproducibility. Reproducibility is a key component of the scientific method and has received increased interest in recent years in the computational modeling community. Numerical simulations rarely ensure this essential property, leading to a low confidence in scientific results and undermining the advancement of knowledge [START_REF] Claerbout | Electronic documents give reproducible research a new meaning[END_REF][START_REF] Stodden | Setting the default to reproducible: Reproducibility in computational and experimental mathematics[END_REF]. While the free availability of a source code is not a sufficient condition to make numerical results reproducible, it is a clear necessity [START_REF] Peng | Reproducible research in computational science[END_REF]. We refer to [START_REF] Harvey | The relevance of open source to hydroinformatics[END_REF], which details the advantages of open-source codes in hydroinformatics and claims that computer software is the lifeblood of hydroinformatics.

Some programs have adopted a strategy consisting of proposing a set of free, open-source codes to solve the shallow-water equations (e.g., GeoClaw [START_REF] Berger | The GeoClaw software for depth-averaged flows with adaptive refinement[END_REF], Gerris [START_REF] Popinet | Quadtree-adaptative tsunami modelling[END_REF], Telemac-Mascaret1 and Dassflow2 ). In an explicit effort to facilitate reproducibility in water-flow modeling and simulations, and considering that "they are no valuable excuses not to make the code available" [START_REF] Barnes | Publish your computer code: it is good enough[END_REF], the present paper describes a free software package for solving the shallow-water equations: FullSWOF. The source code of FullSWOF has been made public under a license that grants the freedom to use, study, share, and modify the code. The use of a standard version of C++ (ANSI) helps in increasing accessibility to numerous users in education, science, and industry.

FullSWOF has been developed through a joint effort between mathematicians and hydrologists. Using a set of analytic solutions to the shallow-water equations detailed in [START_REF] Delestre | SWASHES: a compilation of Shallow-Water analytic solutions for hydraulic and environmental studies[END_REF], FullSWOF includes a validation procedure that guarantees reproducibility to the users and non-regression to the developers. This procedure is, together with the version control system, bug tracking, etc., part of the quality assurance of FullSWOF. These procedures also facilitate contributions by third parties and make their inclusion clearly identifiable. To make the code's use and development easier, a graphical user interface and both a one-dimensional version and a two-dimensional version have been released. FullSWOF can also be included in third-party software, as already realized in openLISEM version 1.67 [START_REF] Baartman | Exploring the role of rainfall variability and extreme events in long-term landscape development[END_REF].

Another specific feature of the implementation of FullSWOF is its modular architecture. Other open-source developments with a modular approach that couples hybrid numerical methods (namely, the finite element and finite volume methods) are presented in [START_REF] Kolditz | Numerical analysis of coupled hydrosystems based on an object-oriented compartment approach[END_REF] and the references therein. Based on the presentation of NOAH (Newcastle Object-oriented Advanced Hydroinformatics), [START_REF] Kutija | An object-oriented approach to the modelling of free-surface flows[END_REF] explained why an object-oriented paradigm provides software that has a high computational efficiency and that is easy to maintain and extend.

As in the GeoClaw and Gerris programs, FullSWOF makes use of finite volume methods, but specific features make it better oriented toward applications in hydrology, hence its name: Full S hallow W ater equations for Overland F low. FullSWOF makes use of the finite volume method, which is preferred to the finite difference method because it ensures both mass conservation and the positivity of water depths. A well-balanced scheme guarantees the preservation of steady states. Special attention has been paid to specific hydrological features: transitions between wet and dry areas, small water depths, and various friction models. The software program incorporates rainfall and infiltration and enables the direct use of digital topographic grids. Some of these features were announced at a conference held in 2012 [START_REF] Delestre | FullSWOF: a software for overland flow simulation[END_REF].

After presenting the physically based model (section 2) and its most important properties (section 3), the numerical methods are described in section 4. Then, we proceed with the description of the FullSWOF program in section 5. Finally, in section 6, comparisons are made with explicit solutions representative of a wide variety of flow conditions, and application to three real cases -a laboratory experiment, rain on a field plot, and a dam break -are reported.

2 Model: Shallow-Water Equations

General Settings

The system of shallow-water equations is a simplified model for a class of free boundary, incompressible Navier-Stokes flows that can occur not only in rivers or channels and in the ocean (tides, tsunamis, etc.) but also in overland flow. They are characterized by the fact that the water depth h(t, x, y) [L] is small with respect to the horizontal dimensions of the considered domain (Figure 1 and Hervouet (2007, chap. 2)). In this context, two main hypotheses are assumed. First, the fluid velocity is constant along the vertical direction; therefore, we can use the horizontal components of the vertically averaged velocity u(t, x, y) and v(t, x, y) [L/T] instead of the three-component Navier-Stokes velocity vector. Next, the pressure of the fluid is hydrostatic; therefore, after integration along the vertical direction z, the pressure field is given by p(t, x, y) = gh(t, x, y) 2 /2, where g is the gravity constant [L/T 2 ].

Under these assumptions, the averaged Navier-Stokes system can be rewritten as the following balance laws:

             ∂ t h + ∂ x (hu) + ∂ y (hv) = R -I, ∂ t (hu) + ∂ x hu 2 + gh 2 2 + ∂ y (huv) = gh(S 0x -S f x ), ∂ t (hv) + ∂ x (huv) + ∂ y hv 2 + gh 2 2 = gh(S 0y -S f y ).
(1)

The first equation is the exact integrated form of the incompressibility condition, and hence, it is a mass balance. The other two equations are momentum balances and involve forces such as gravity and friction. In particular, the gh 2 /2 term is the hydrostatic pressure. Let us now describe each term, recalling their physical dimensions.

1. z b is the topography [L]. Because we do not consider erosion, it is a given function of space, z b (x, y), and we classically denote by S 0x and S 0y the opposites of the slopes in the x and y directions, respectively:

S 0x = -∂ x z b (x, y) and S 0y = -∂ y z b (x, y). 2. R is the rain intensity [L/T]. It is a given function R(t, x, y) ≥ 0.
In the current versions of FullSWOF, we consider the rain to be uniform in space.

3. I is the infiltration rate [L/T]. This term, I(t, x, y) ≥ 0, is defined through the coupling with an infiltration model such as the bi-layer Green-Ampt model (section 2.3).

4. S f = (S f x , S f y ) is the friction force, which is, in general, a nonlinear function of the velocity and water depth (section 2.2).

We shall pay particular attention to the one-dimensional version of system (1) because, on the one hand, it has practical applications when the flow can be considered homogeneous or when the effect of the edges can be neglected (e.g., wide channels or flood propagation in river networks). On the other hand, its study gives a better insight into the complete twodimensional model from both theoretical and numerical perspectives. This one-dimensional system can be written as

     ∂ t h + ∂ x (hu) = R -I, ∂ t (hu) + ∂ x hu 2 + gh 2 2 = gh(S 0x -S f x ).
(2)

In both systems (1) and ( 2), the homogeneous part, that is, the left-hand side, is called the transport (or convection) operator. This part corresponds to the flow of an ideal fluid on a flat bottom, without friction, rain or infiltration. In the one-dimensional setting, this is exactly the model introduced by Barré de [START_REF] De Saint-Venant | Théorie du mouvement nonpermanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF]. This operator contains several important properties of the flow; hence, in sections 3 and 4, we perform a careful analysis of the homogeneous system in the following, considering both theoretical and numerical studies.

Friction Terms

Friction terms depend on the flow velocity. In the formulae below, u is the velocity vector u = (u, v), with | u| = √ u 2 + v 2 , and q is the discharge q = (hu, hv) = h u. In hydrological models, based on empirical considerations, friction laws can be categorized into two families. On the one hand, the Manning-Strickler friction law reads

S f = C f u| u| h 4/3 = C f q| q| h 10/3 , (3)
with , 1959). On the other hand, the laws of Darcy-Weisbach and Chézy can be written as

C f = n 2 or C f = 1/K 2 , where n is the Manning coefficient [L -1/3 T] and K is the Strickler coefficient [L -1/3 T -1 ] (Chow
S f = C f u| u| h = C f q| q| h 3 . ( 4
)
Taking C f = f /(8g), with f as a dimensionless coefficient, or

C f = 1/C 2 , C [L 1/2 T -1 ],
we obtain the Darcy-Weisbach and Chézy friction laws, respectively. Readers are referred to chapter 5 of [START_REF] Chow | Open-Channel Hydraulics[END_REF] for details about and examples of friction laws. Note that the friction force may depend on the space variable, especially in large domains, but this is not considered in the following.

Infiltration Model

Infiltration is computed at each cell using a Green-Ampt model [START_REF] Green | Studies on soil physics. part I.-the flow of air and water through soils[END_REF][START_REF] Mein | Modeling infiltration during a steady rain[END_REF]. The main idea is to assume that the water infiltrates as an advancing wetting front (at the depth

Z f = Z f (t))
from a fully saturated zone (with moisture content θ s ) to another zone with an initial water content θ i (Figure 2a). Following [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF], we implemented a bi-layer Green-Ampt model [START_REF] Hillel | Transient infiltration into crust-topped profiles[END_REF][START_REF] Delestre | Rain water overland flow on agricultural fields simulation[END_REF] in which the upper layer is characterized by its thickness Z c and its hydraulic conductivity K c , and the second layer has an infinite extension and a hydraulic conductivity K s (Figure 2b). Note that the 'c' and 's' subscripts stand for 'crust' and 'soil', respectively. Z c can be set equal to zero for cases where a single layer model is relevant. In the following, we assume that the infiltration parameters can vary in space but are uniform in time.

wetting front O h n surf θ s -θ i Z θ Z n f θ i θ s (a) Wetting front. K s Z c K c (
At each cell, when there is water at the surface, the infiltration capacity I C [L/T] at time t n is given by

I C (t n ) = I n C =                        K s 1 + h f -h n surf Z n f if Z c = 0, K c 1 + h f -h n surf Z n f if Z n f ≤ Z c , K n e 1 + h f -h n surf Z n f else.
(

) 5
In these equalities, K n e = K e (t n ) is the effective hydraulic conductivity at time t n :

K n e = Z n f Z n f -Z c K s + Z c K c = 1 1 K s 1 -Z c ∆θ V n inf + Z c ∆θ V n inf 1 K c .
The value of h f depends on the soil; it is sometimes denoted by Ψ in the literature and represents the suction head at the wetting front. The (positive) quantity -h n surf = -h surf (t n ) is the water depth at the surface of the cell that is available for infiltration at time t n . Finally, Z n f = Z f (t n ) is the gravity force over the water column at time t n and can be written as Z n f = V n inf /∆θ, where V n inf is the infiltrated volume at time t n and ∆θ = θ sθ i . To avoid an infinite infiltration rate initially (when the infiltrated volume is still equal to zero), we add a threshold to obtain the infiltration rate I n = min(I n C , i max ). Because the infiltrated volume cannot exceed the water depth, the volume is updated as follows:

V n+1 inf = V n inf + min(-h n surf , I n × ∆t)
Finally, the water depth is updated.

Properties

In this section, we recall several mathematical properties of the shallow-water model for both the 1D and 2D cases.

One-Dimensional Model

To emphasize the mathematical properties of the shallow-water model, we first rewrite the one-dimensional homogeneous equations using vectors:

∂ t W +∂ x F (W ) = 0, where W = h hu , F (W ) =   hu hu 2 + gh 2 2   , (6)
with F (W ) being the flux of the equation. The transport is more clearly evidenced in the following non-conservative form:

∂ t W + A(W )∂ x W = 0, A(W ) = F (W ) = 0 1 -u 2 + gh 2u ,
where A(W ) is the matrix of transport coefficients. More precisely, when h > 0, the matrix A(W ) turns out to be diagonalizable, with eigenvalues

λ 1 (W ) = u - √ gh < u + √ gh = λ 2 (W )
. This important property of having two real and distinct eigenvalues is called strict hyperbolicity (e.g., see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] and references therein for more details). The eigenvalues are the velocities of the surface waves of the fluid, which are the fundamental characteristics of the flow. Note that the eigenvalues coincide if h = 0, that is, for dry zones. In this case, the system is no longer hyperbolic, which induces difficulties at both the theoretical and numerical levels.

From these formulae we recover a useful classification of flows based on the relative values of the velocities of the fluid, u, and of the waves, √ gh. If |u| < √ gh, the characteristic velocities u -√ gh and u + √ gh have opposite signs, and information propagates upward as well as downward; the flow is said to be subcritical or fluvial. In contrast, when |u| > √ gh, all the information propagates downward, and the flow is said to be supercritical or torrential.

This classification has consequences for the numerical scheme. Because we have two unknowns h and u (or, equivalently, h and q = hu), a subcritical flow is determined by one upstream value and one downstream value, whereas a supercritical flow is completely determined by the two upstream values. Thus, for numerical simulations, we use only one of the two variables for a subcritical inflow/outflow boundary. For a supercritical inflow boundary, we have to set both variables, and for a supercritical outflow boundary, the Neumann free-boundary conditions are considered (e.g., [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF]. In this context, it is useful to be able to determine whether the flow is subcritical or supercritical. To this end, we can consider two quantities. The first quantity is the Froude number:

Fr = |u| √ gh .
The flow is subcritical or supercritical if Fr < 1 or Fr > 1, respectively.

A more visual criterion is obtained through the so-called critical depth h c , which can be written as

h c = |q| √ g 2/3
, for a given discharge q = hu. The flow is subcritical or supercritical if h > h c or h < h c , respectively.

Source Terms and Equilibria

When source terms (e.g., topography, rain or friction) are involved, other properties have to be considered, in particular, the occurrence of steadystate (or equilibrium) solutions (solutions that do not depend on time, i.e., ∂ t ≡ 0). This amounts to some balance between the flux and source terms:

∂ x (hu) = R -I, ∂ x hu 2 + gh 2 2 = gh(S 0x -S f x ).
Among these solutions, the homogeneous states are of particular interest, namely, the so-called parallel flows

∂ x (hu) = R -I = 0, gh(S 0x -S f x ) = 0,
and lakes or puddles at rest

u = 0, h + z = const.
These solutions are important for two reasons: on the one hand, specific numerical methods have to be designed to obtain these solutions; on the other hand, these models furnish several explicit solutions that can be used as test cases for numerical methods (see section 6 and Delestre et al. ( 2013)).

Two-Dimensional Model

The two-dimensional shallow-water system (1) can be written under the following conservative form:

∂ t U + ∂ x G(U ) + ∂ y H(U ) = S(U, t, x, y), (7)
where

U =   h hu hv   , G(U ) =     hu hu 2 + gh 2 2 huv     , H(U ) =     hv huv hv 2 + gh 2 2     and S(U, t, x, y) =   R -I gh(S 0x -S f x ) gh(S 0y -S f y )   .
If we denote by DG x and DH y the Jacobian matrices of the fluxes, namely,

DG x =   0 1 0 -u 2 + gh 2u 0 -uv v u   and DH y =   0 0 1 -uv v u -v 2 + gh 0 2v   , then system (7) reads ∂ t U + DG x ∂ x U + DH y ∂ y U = S(U, t, x, y).
The notion of hyperbolicity is defined here, following [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], by studying the flow along any direction. For any unit vector ξ = (ξ x , ξ y ) ∈ R 2 , the velocity of the flow in the ξ direction is, by definition, u ξ = ξ x u + ξ y v, and we define a type of directional derivative as DF (ξ) = ξ x DG x + ξ y DH y . Then, it can be verified that DF (ξ) has three eigenvalues:

λ 1 (ξ) = u ξ -gh, λ 2 (ξ) = u ξ and λ 3 (ξ) = u ξ + gh. (8)
Outside the dry zones, that is, when h > 0, the three eigenvalues satisfy λ 1 (ξ) < λ 2 (ξ) < λ 3 (ξ); therefore, the matrix DF (ξ) is diagonalizable for all ξ. Specifically, the system ( 7) is strictly hyperbolic. As in the onedimensional case, this property is no longer true inside the dry zones. source terms. First, we perform a convective step, where the homogeneous part of the system is solved by a finite volume strategy. This leads to a first-order-accurate scheme, the second-order accuracy in space being obtained by reconstruction techniques. Second, this scheme is coupled with the source terms, and specific methods are used to consider steady states. In case a second-order approximation in time is requested, Heun's method is used. These steps are described below in detail for the one-dimensional setting and are then extended to the two-dimensional case.

General Settings

Let us first introduce several notations that will be used throughout this section. First, consider the time discretization: let t 0 = 0 be the initial time; we fix a time step ∆t > 0, and for n ≥ 0, we set t n+1 = t n + ∆t. Next, the space discretization is defined by constant positive space steps ∆x in one dimension and (∆x, ∆y) in two dimensions, where a rectangular mesh is assumed. Finally, we compute piecewise constant approximations of the vectors 3 for details of the notations).

W (in 1D) or U (in 2D). More precisely, W n i and U n ij are constant approximations of W on [t n , t n+1 [×]x i-1/2 , x i+1/2 [ and of U on [t n , t n+1 [×]x i-1/2 , x i+1/2 [×]y j-1/2 , y j+1/2 [, respectively (see Figure
x i+1 x i-1 x i-1/2 x i x i+1/2 x ∆x t n+1 O t t n ∆t
(a) One-dimensional time and space discretization.

x i-1 x i-1/2 x i x i+1/2 x ∆x O x i+1
y ∆y y j-1/2 y j y j+1/2 (b) Two-dimensional space discretization for every time step. For each n ≥ 1, we compute the components of W n and U n using explicit schemes in time that take the general recursive form

W n+1 i = W n i -∆tΦ(W n ), U n+1 ij = U n ij -∆tΦ(U n ). (9)
Formula ( 9) is nothing more than an Euler scheme in time, where the function Φ describes the discretization of the space derivatives and of the source terms.

At this point, we briefly recall several fundamental properties required in an effective numerical scheme (for detailed definitions, see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]). First, the consistency error quantifies the ability of the scheme to mimic the equation. This error is obtained by replacing each occurrence of W n i in the formula (9) by the exact solution W (t n , x i ). If this error goes to 0 when ∆t and ∆x → 0, the scheme is consistent. If the error behaves as O(∆t q + ∆x p ), the scheme is of order q in time and p in space. Next, stability, a notion that can take various forms, is a concern. Usually, it is formulated by stating that a given norm of the numerical solution at time t n is controlled by the norm at time 0. Regarding the shallow-water equations, a particularly important stability property is the preservation of the positivity of the water depth, especially if wet-dry transitions and thin water layers are to be simulated.

As far as precision in time is concerned, note that formula (9) defines a first-order scheme in time, regardless of the definition of Φ. There are several ways to increase the accuracy of such schemes following standard methods to solve ordinary differential equations. Here, we choose Heun's method (also referred to as the modified Euler method), which is a prediction-correction method. In the first two formulae below, W * i and W * * i are the predicted values; the last formula is the correction step:

W * i = W n i -∆tΦ(W n ), W * * i = W * i -∆tΦ(W * ), W n+1 i = W n i + W * * i 2 .
We conclude this paragraph by emphasizing that when Φ is determined, the numerical scheme is complete. Therefore, the remainder of this section explains step by step how to build the function Φ, which encodes all the stability and accuracy properties of the scheme in the space variables.

Convective Step for the One-Dimensional Model

To solve system (2), the first step consists of using a finite volume strategy to solve the homogeneous system [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. The idea is to integrate the system of equations in each time-space cell, as depicted by the blue rectangle in Figure 3a. We obtain the following approximation formula

W n+1 i = W n i - ∆t ∆x (F n i+1/2 -F n i-1/2 ), (10)
where F n i+1/2 and F n i-1/2 are approximations of the flux on the edges x i+1/2 and x i-1/2 , respectively. At this stage, the finite volume scheme is not completely determined yet; we need to specify a formula to compute this flux approximation at each edge. Here, we choose to set, for all indexes i,

F n i+1/2 = F(W n i , W n i+1 )
, where F(W L , W R ) is called the numerical flux, a given function of the states on the left and right of the interface.

Among the numerical fluxes proposed in the literature, such as Rusanov, HLL, VFRoe-ncv, and kinetic (see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF] for explicit formulae and additional references), several fluxes are available in FullSWOF. In this paper, we choose to present the HLL flux (introduced in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]) because this flux offers a good compromise between simplicity, computing time and robustness, as shown numerically in [START_REF] Delestre | Rain water overland flow on agricultural fields simulation[END_REF]. This flux can be written as

F (W L , W R ) =        F (W L ) if 0 ≤ c 1 c 2 F (W L ) -c 1 F (W R ) c 2 -c 1 + c 1 c 2 c 2 -c 1 (W R -W L ) if c 1 < 0 < c 2 F (W R ) if c 2 ≤ 0 , (11)
with two parameters c 1 < c 2 that are the approximations of the slowest and fastest wave speeds. We refer to [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF] for further discussion on the wave speed estimates. In FullSWOF, we have implemented

c 1 = inf W =W L ,W R ( inf j∈{1,2} λ j (W )) and c 2 = sup W =W L ,W R ( sup j∈{1,2} λ j (W )),
where λ 1 (W ) = u -√ gh and λ 2 (W ) = u + √ gh are the eigenvalues of the one-dimensional model (section 3.1). This numerical flux is said to be upwind, i.e., it mimics the wave propagation. When the flow is supercritical (0 ≤ c 1 or c 2 ≤ 0), all information is going from upstream to downstream; thus, the numerical flux is calculated from upstream values only. In contrast, for subcritical flows (c 1 < 0 < c 2 ), information is coming from both upstream and downstream; thus, the flux is calculated using both upstream and downstream values.

As described above, except for certain numerical fluxes, the scheme is first-order accurate in space. To obtain second-order accuracy, we perform a linear reconstruction in space on the variables W = (h, hu), thus obtaining new variables on each interface i + 1/2, namely, W i+1/2-on the left and W i+1/2+ on the right. Using F(W i+1/2-, W i+1/2+ ) in the finite volume scheme instead of F(W n i , W n i+1 ) leads to a second-order approximation in space. Several of the existing formulae for the linear reconstruction (e.g., MUSCL and ENO, see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF]) are implemented in FullSWOF.

We conclude this section with an important remark about time discretization. Explicit schemes imply a control on the time step, which cannot be too large for a given space step. More precisely, for a 3-point scheme, as defined above (W n+1 i depends only on W n i-1 , W n i and W n i+1 ), the numerical speed of propagation is ∆x/∆t. To avoid any loss of information, this velocity has to be larger than any possible physical velocity:

C ∆x ∆t ≥ sup (t,x) |u(t, x)| + gh(t, x) ,
where the supremum is taken over the whole time-space domain of interest and C is a parameter depending on the dimension and on the order of the schemes that are considered (in 1D, at the first order, C = 1, and, at the second order, C = 0.5; in 2D, at the first order, C = 0.5, and, at the second order, C = 0.25). When this condition is violated, the scheme becomes unstable, and oscillations appear. This theoretical formulation of the limitation, known as the CFL condition (Courant, Friedrichs, Lewy) [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], is not useful because, due to the nonlinearity, the righthand side is difficult to estimate. Moreover, even if such an estimate is obtained, it can lead to an underestimated fixed time step when there are no large variations in h and u. Thus, the computation of the time step could be replaced by a sequence of variable time steps ∆t n according to the following rule:

∆t n = C∆x sup i {|u n-1 i | + gh n-1 i } . ( 12
)
However, when water depths converge toward zero, the time steps will tend to become infinite, which is not reasonable. Consequently, we add an upper bound to the computation, which leads to the following formulation (implemented in FullSWOF_2D):

∆t n = C min     ∆x, ∆x sup i {|u n-1 i | + gh n-1 i }     .

Source Terms

When source terms are involved, specific methods have to be introduced to obtain equilibrium states. Numerical schemes that preserve stationary equilibria are called equilibrium schemes, or well-balanced schemes [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equation[END_REF]. The first class of methods consists of splitting-type methods, where the transport equation and source term are solved independently. Another strategy consists of applying a specific reconstruction technique, thus modifying the flux computation. Apparently, the latter method is well adapted to the topography term but not quite to the friction term. Hence, in FullSWOF, the following two approaches are used: the topography is solved using a well-balanced scheme (section 4.3.1), and the friction is solved with a semi-implicit method (section 4.3.2).

The Hydrostatic Reconstruction

The hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF], as its name suggests, reconstructs new variables to be used in the numerical flux.

It is designed to obtain a well-balanced scheme in the sense that it preserves, at least, the steady state at rest (i.e., the hydrostatic equilibrium) as well as the positivity of the water depth. The hydrostatic reconstruction procedure is first applied to the original variables W n i , thus leading to a first-order scheme.

Here, we detail how to obtain a second-order version, which is used in FullSWOF. The second-order scheme consists of first performing the linear reconstruction (chosen for the convective step) not only to W = (h, hu), as before, but also to h + z b . Then, the hydrostatic reconstruction is applied to these modified variables. Note that this strategy introduces an artificial time dependence on the topography, which is in some sense reconstructed as well. This is mandatory for addressing equilibrium solutions and for preserving the positivity of h [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. The linear reconstruction gives the values (h i+1/2-, z i+1/2-, u i+1/2-) on the left of the interface i + 1/2 and (h i+1/2+ , z i+1/2+ , u i+1/2+ ) on the right of the interface. The final formula for the hydrostatic reconstruction can now be written as

               h i+1/2L = max h i+1/2-+ z i+1/2--max z i+1/2-, z i+1/2+ , 0 , W i+1/2L = h i+1/2L h i+1/2L u i+1/2- , h i+1/2R = max h i+1/2+ + z i+1/2+ -max z i+1/2-, z i+1/2+ , 0 , W i+1/2R = h i+1/2R h i+1/2R u i+1/2+ . (13)
A given space discretization may exhibit abnormal behaviors for some combinations of slopes and water depths [START_REF] Delestre | A limitation of the hydrostatic reconstruction technique for Shallow Water equations[END_REF]. Particularly obvious for the first-order scheme and on a coarse mesh, they disappear when refining the mesh and are hardly noticeable at second order.

When applying the hydrostatic reconstruction, formula (10) has to be modified to preserve the consistency of the scheme. This can be written as

W n+1 i = W n i -∆tΦ(W n ), with Φ(W n ) = 1 ∆x F n i+1/2L -F n i-1/2R -F c n i ,
where F n i+1/2L and F n i-1/2R are given by

                 F n i+1/2L = F(W n i+1/2L , W n i+1/2R ) +   0 g 2 h n i+1/2- 2 -h n i+1/2L 2   , F n i-1/2R = F(W n i-1/2L , W n i-1/2R ) +   0 g 2 h n i-1/2+ 2 -h n i-1/2R 2   ,
where W n i+1/2L and W n i+1/2R are computed using formula (13). The additional centered term F c n i is determined to preserve consistency and to ensure a well-balanced scheme [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]:

F c n i = 0 - g 2 h n i-1/2+ + h n i+1/2- z n i+1/2--z n i-1/2+
.

Friction

Friction may be addressed in a well-balanced scheme by first introducing the friction in the topography term and by then applying the hydrostatic reconstruction (section 4.3.1); this approach is therefore named "apparent topography". This method has been used, for example, to solve a shallowwater system with a Coriolis force [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF] and with Coulomb friction [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF][START_REF] Mangeney | Numerical modeling of self-channeling granular flows and of their levee-channel deposits[END_REF]. The main idea is to use the modified topography z app defined by z app = z bb, with ∂ x b = S f x . The use of this class of methods gives rise to a well-balanced scheme for friction as well as for the topography and thus computes equilibrium states correctly. However, the solutions are not completely satisfactory for transitory situations, as noted in [START_REF] Delestre | Simulation of rain-water overland-flow[END_REF] and Delestre and James (2010); a spurious peak appears at the wet-dry front before equilibrium is reached. Therefore, we utilize splitting methods. The explicit discretization, despite its simplicity, is not relevant for the type of problems we are interested in because this discretization leads to instabilities and overestimations of the velocity at wetdry interfaces [START_REF] Paquier | Modeling and simulation of the propagation of the dambreak wave[END_REF]. However, the fully implicit method has a high computational cost and cannot be generalized to the two-dimensional system due to the nonlinear form of the friction law. Another possibility consists of using a Strang time splitting to obtain the second order in time calculation; however, this leads to a more complicated algorithm, and moreover, there is no significant gain in accuracy compared to first-order methods [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. Finally, the best compromise between accuracy, stability and computational complexity lies in semi-implicit methods [START_REF] Fiedler | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF][START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. We choose the semi-implicit treatment proposed in [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF] not only because it preserves steady states at rest but also for its stability. At first order in time, after a convective step

W * i = W n i -∆tΦ(W n ), the value W n+1 i is given by W n+1 i = h n+1 i q n+1 i =    h * i q * i 1 + gn 2 ∆t |q n i | h n i h n+1 i 4/3 -1   
for the Manning friction law (3) and by

W n+1 i = h n+1 i q n+1 i =   h * i q * i 1 + ∆t f 8 |q n i | h n i h n+1 i -1  
for the Darcy-Weisbach friction law (4). Note the simplicity of the method, which gives an explicit value for W n+1 i .

Rain and Infiltration

Unlike the friction and topography source terms, rain and infiltration produce no particular numerical difficulties, such as steady-state or stability preservation. Moreover, addressing infiltration implicitly would complicate the integration of other infiltration models (such as Richard's or Darcy's models) into FullSWOF. For these reasons, we have chosen to treat the rain and infiltration terms explicitly.

Convective Step for the Two-Dimensional Model

In general, the 2D system can be treated in the same manner as the 1D system because the calculations are performed for each interface of each cell (and thus do not depend on the cell geometry). Digital elevation models (DEMs) (i.e., digital topographic maps) are mainly represented as structured grids. Some DEMs can be represented as vector-based triangular networks (TINs). Because TINs can easily be converted into DEMs, we chose to make developments dedicated to structured meshes (Figure 3b). To obtain a 2D version, we perform the linear reconstruction of U and z b and the hydrostatic reconstruction as in one dimension (equation ( 13)). We obtain U •L,• and U •R,• along the x direction, and we obtain U •,•L and U •,•R along the y direction.

Then, the two-dimensional finite volume scheme reads

U * i,j = U n i,j - ∆t ∆x G n i+1/2L,j -G n i-1/2R,j -Gc n i,j - ∆t ∆y H n i,j+1/2L -H n i,j-1/2R -Hc n i,j with                G n i+1/2L,j = G U n i+1/2L,j , U n i+1/2R,j + S n i+1/2L,j , G n i-1/2R,j = G U n i-1/2L,j , U n i-1/2R,j + S n i-1/2R,j , H n i,j+1/2L = H U n i,j+1/2L , U n i,j+1/2R + S n i,j+1/2L , H n i,j-1/2R = H U n i,j-1/2L , U n i,j-1/2R + S n i,j-1/2R
, where

S n i+1/2L,j =     0 g 2 h n i+1/2-,j 2 -h n i+1/2L,j 2 0     , S n i-1/2R,j =     0 g 2 h n i-1/2+,j 2 -h n i-1/2R,j 2 0     , S n i,j+1/2L =     0 0 g 2 h n i,j+1/2- 2 -h n i,j+1/2L 2     , S n i,j-1/2R =     0 0 g 2 h n i,j-1/2+ 2 -h n i,j-1/2R 2     .
G and H are the numerical fluxes used to obtain the solution of the homogeneous system. The first two components of G, G 1 and G 2 , as well as the first and third components of H, H 1 and H 3 , are computed as the components of F. In the HLL Riemann solver (equation ( 11)), being a two-wave numerical flux solver, this assumption is correct only for hyperbolic systems of two equations, such as the one-dimensional shallow-water equations. In two spatial dimensions, there are three equations and hence three eigenvalues (equation ( 8)). In this case, the HLL solver described in equation ( 11) is not sufficiently precise because it only involves two waves. To address this issue, Toro proposed an extension of HLL called the HLLC solver [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]. The solver we propose in FullSWOF is inspired by this reference and defines G 3 and H 2 as

G 3 U n i+1/2L,j , U n i+1/2R,j =    v n i+1/2L,j G 1 U n i+1/2L,j , U n i+1/2R,j if u n i+1/2L,j + u n i+1/2R,j > 0, v n i+1/2R,j G 1 U n i+1/2L,j , U n i+1/2R,j if u n i+1/2L,j + u n i+1/2R,j ≤ 0, H 2 U n i,j+1/2L , U n i,j+1/2R =    u n i,j+1/2L H 1 U n i,j+1/2L , U n i,j+1/2R if v n i,j+1/2L + v n i,j+1/2R, > 0, u n i,j+1/2R H 1 U n i,j+1/2L , U n i,j+1/2R if v n i,j+1/2L + v n i,j+1/2R, ≤ 0.

Description of the FullSWOF Software Package

FullSWOF stands for "Full Shallow-Water equations for Overland Flow". The names FullSWOF_1D and FullSWOF_2D define the one-dimensional and two-dimensional versions, respectively.

Having been designed to encourage research reproducibility, the source codes (in C++) for FullSWOF are available and can be downloaded from the websites https://sourcesup.renater.fr/projects/fullswof-1d/ for the one-dimensional version and from https://sourcesup.renater.fr/projects/ fullswof-2d/ for the two-dimensional version. The common optional graphical interface is named FullSWOF_UI ("UI" for User Interface), was developed in Java, and is available for download at https://sourcesup.renater. fr/projects/fullswof-ui/. Each piece of software is distributed under the CeCILL-V2 (GPL-compatible) free software license, which allows use of the software package without any limitation as to its field of application. For more details, we refer to the documentation on the websites.

The structure of the source code is designed to make future development easy, especially for new developers; for example, a new friction law can eas-ily be added to the libfriction library by creating a new friction file. The documentation for programmers is included directly in the C++ code using doxygen specific comments [START_REF] Van Heesch | Manual for version 1[END_REF]. This leads to the automatic extraction of the programmer's manual and simplifies the documenting task. Due to the hosting on a software forge, FullSWOF obtains the benefits of systems of version control, bug tracking, package release, etc.

The one-dimensional program is primarily designed, of course, to perform 1D flow simulations. It is also used as a development tool to test new numerical methods (e.g., fluxes) and to introduce new features into the models (such as new friction and infiltration laws). The new code can then be later integrated into the two-dimensional program; because the 2D mesh is a structured mesh, it is easy to adapt the code from 1D to 2D. In addition, a parallelized version (MPI) of FullSWOF_2D has been developed to run large test cases (such as the Malpasset dam break, section 6). Currently considered in an early stage of development, this version is available only in a dedicated branch of the version control system.

The FullSWOF software is designed to use point-wise-defined topographies and/or initial conditions. Currently, several topographies and initialization values are hard coded (e.g., a parabolic topography and the wet dam-break initial data). For other cases, input text files can be read. The friction law can be chosen between the Manning and Darcy-Weisbach laws. Several boundary conditions are available (e.g., wall, Neumann, periodic, and imposed depth). The rain is constant in space but can vary in time. Several numerical methods are implemented for the flux, the linear reconstruction and the order of the scheme. Based on a study of numerical methods for overland flows [START_REF] Delestre | Rain water overland flow on agricultural fields simulation[END_REF], default values have been selected for an HLLtype method for the flux, the MUSCL formula for the reconstruction and a second-order scheme. We refer to the documentation for more details.

Some classical benchmarks and analytic solutions from the literature are defined in FullSWOF. They have been chosen among those gathered in SWASHES, a compilation of shallow-water analytic solutions for hydraulic and environmental studies [START_REF] Delestre | SWASHES: a compilation of Shallow-Water analytic solutions for hydraulic and environmental studies[END_REF]. They are used to validate each new version of the code before its release, assuring the quality of the software. This validation is partially automated due to a dedicated script and is useful to both users (who can verify that their results do not differ from reference results -and do not depend on their compiler, operating system or hardware) and developers (who can verify that their changes do not cause a regression in the result quality). Additional test cases can be added to validate new conditions introduced in the code.

Numerical Illustrations

In this section, we illustrate some results from FullSWOF using a few classical test cases, using analytic solutions available in the literature, and using three real datasets. This is by no means an analysis of the performance of FullSWOF. We merely attempt to demonstrate the ability of FullSWOF to simulate a wide range of flow conditions. In particular, real datasets (section 6.3) are used without any filtering of the data (no smoothing, etc.). Note that the illustrated test cases and analytic solutions are part of the benchmark set included in the FullSWOF_1D or FullSWOF_2D codes. In the following, FullSWOF_1D version 1.01.00 (2013-05-17) and FullSWOF_2D version 1.04.08 (2013-11-07) have been used.

Classical Test Cases

We begin by running FullSWOF on three analytic cases: one steady-state and two classical transitory solutions. These examples correspond to the original shallow-water system with a spatially varying topography; there is no friction, no rain and no infiltration.

Lake at Rest with an Emerged Bump

This one-dimensional steady-state solution has been developed in [START_REF] Delestre | Rain water overland flow on agricultural fields simulation[END_REF] and in Delestre et al. (2013, §3.1.2) as a test case for the preservation of steady states and the boundary condition treatment. This solution is based on the topography provided by [START_REF] Goutal | Proceedings of the 2 nd workshop on dambreak wave simulation[END_REF]. The initial condition satisfies the hydrostatic equilibrium h + z = const and q = 0 m 2 s -1 .

(

The domain length is set to L = 25 m, with a topography given by

z(x) = 0.2 -0.05(x -10) 2 if 8 m < x < 12 m, 0 else.
The maximum water depth is smaller than the amplitude of the topography to simulate a lake at rest with an emerged bump (Figure 4a). In such a configuration, starting from the steady state, the velocity must remain null, and the water surface should stay flat. This is the exact behavior simulated by FullSWOF_1D, showing the interest in a well-balanced scheme (Figure 4). Figure 4: Lake at rest with an emerged bump: comparison of the analytic solution with the FullSWOF_1D results. Simulation with 500 cells at T = 100 s.

Dam Break on a Dry Domain without Friction

We turn now to a transitory one-dimensional case, namely, the analytic solution of a dam break on a dry domain without friction on a flat and horizontal topography (Delestre et al., 2013, §4.1.2). This case is known as Ritter's solution [START_REF] Ritter | Die Fortpflanzung der Wasserwellen[END_REF][START_REF] Hervouet | Hydrodynamics of free surface flows[END_REF]. The difficulties here are (a) the existence of a shock and (b) a wet-dry transition. This case also tests whether the scheme preserves the positivity of the water depth because this property is usually violated near the wetting front.

The initial condition for this configuration is the following Riemann problem

h(x) = h l = 0.005 m for 0 m ≤ x ≤ 5 m, h r = 0 m for 5 m < x ≤ 10 m, with u(x) = 0 m s -1 .
Initially, the free surface exhibits the following structure: starting from upstream, there is a constant water depth at rest (with h = 0.005 m) connected by a parabola to a dry zone downstream (with h = 0 m). The left extremity of the parabola moves upstream, while its right end slides downstream. Figure 5 displays this solution at time t = 6 s.

There is an overall agreement between the FullSWOF_1D results and the analytic solution; the code is able to represent the shock, locate and correctly treat the wet-dry transition, and preserve the positivity of the water depth. However, differences are observed at the two connections between the constant states and the parabola. These differences become smaller when the space step is decreased (results not shown).

Planar Surface Rotating in a Paraboloid

Finally, we present a transitory two-dimensional solution: a planar water surface rotating in a paraboloid (Delestre et al., 2013, §4.2.2). In the literature, this analytic solution is known as Thacker's 2D case, named after its author [START_REF] Thacker | Some exact solutions to the nonlinear shallow-water wave equations[END_REF]. The topography is a paraboloid of revolution, and the shoreline is a moving circle. The free surface exhibits periodic motion and remains planar in time. To visualize this case, wine being swirled in a glass can be considered (cross-section in Figure 6a). This is a solution with a variable slope (in space) for which the wet-dry transitions are constantly moving. Because there is no friction, the rotation should not decrease over time. Hence, this case tests both the ability of the schemes to simulate flows with comings and goings and the numerical diffusion of the scheme (which causes a damping in the water depth over time). The analytic solution at t = 0 s is taken as the initial condition for the computation. The results are considered after three periods. Overall, FullSWOF_2D produces results in good agreement with the analytic solution (Figure 6). There is no spurious point at the wet-dry transitions, exemplifying the capability of the well-balanced schemes to properly compute these common situations for natural surface flows. The slope of the water surface given by FullSWOF_2D is slightly lower than that of the analytic solution. In addition, the water flow rate is slightly lower than expected. This is a consequence of the numerical diffusion of the scheme. Several techniques can be applied to improve the results of this specific case, but they may not be relevant for the general case because they would increase the complexity of the method and the computational cost. Simulation with 100 cells in x and y at T = 13.4571 s.

MacDonald-Type Solutions

We turn now to a class of stationary solutions, more complex than the previous solutions (section 6.1) in the sense that the model is now complete, with friction and rain. These solutions are obtained by the procedure introduced by I. MacDonald: the water depth profile and the discharge are given, and the corresponding topography is then computed. In the original works, the Manning friction law was considered [START_REF] Macdonald | Analysis and computation of steady open channel flow[END_REF][START_REF] Macdonald | Analytic benchmark solutions for open-channel flows[END_REF], but the method allows many variants: other friction laws, rain, diffusion, etc. The three selected solutions show the ability of the scheme to address stationary states induced by the topography and by friction under a wide range of flow conditions.

Short Channel with a Smooth Transition and a Shock

The length of the channel is 100 m, and the discharge at steady state is q = 2 m 2 s -1 (Figure 7). The flow is subcritical both upstream and downstream (Delestre et al., 2013, §3.2.2). In the intermediate part, the flow is supercritical. Hence, the flow goes from subcritical to supercritical via a sonic point and then -through a shock (located at x = 200/3 ≈ 66.67 m) -becomes subcritical again (Figure 7a). The Manning friction coefficient n is equal to 0.0328 m -1/3 s. The run was performed at a resolution of 0.2 m.

The water flow rate simulated by FullSWOF_1D matches the analytic solution except at the shock (Figure 7b). A careful examination shows that the mismatch is limited to two cells. This mismatch is due to the product of approximations of two discontinuous functions (u and h) on a discretization grid. Further testing showed that the shock was better represented with finer grids, which is the expected behavior. Future developments -by us or by others, considering that FullSWOF is open source -could improve the results of this test case. x (m)

Rain on a Long Channel with a Supercritical Flow

This analytic case is similar to the previous case but includes rain (Delestre et al., 2013, §3.3.2). The channel length is set equal to 1000 m (Figure 8). Because the flow is supercritical along the entire channel, we consider a constant discharge and a constant water depth at inflow and a free outflow.

In the simulation, the channel is initially dry. There is no rain until 1500 s; after this time, the rain intensity is set to 0.001 m s -1 until the end of the simulation. The Darcy-Weisbach friction coefficient is f = 0.065, and the inflow discharge is q 0 = 2.5 m 2 s -1 (Figure 8).

The results from FullSWOF_1D are in good agreement with the analytic solution, which demonstrates that the inclusion of the rain source term in the model and its implementation in the code are satisfactory.

Pseudo-2D Channel with a Supercritical Flow

In this section, we consider a pseudo-2D shallow-water system in a channel with a varying width. By pseudo-2D, we mean an intermediate model between the one-dimensional and two-dimensional models. More precisely, the pseudo-2D shallow-water equations are obtained by averaging all quantities both in the vertical direction and over the width of the channel (y direction), i.e., perpendicular to the flow (recall that the two-dimensional shallow-water equations involve an average solely in the vertical direction).

The derivation is detailed in [START_REF] Goutal | A kinetic interpretation of the section-averaged Saint-Venant system for natural river hydraulics[END_REF]. In that paper, the authors showed that the new terms that appear in the pseudo-2D shallow-water system due to the y direction averaging procedure produce numerical difficulties that must be treated carefully. We perform a simulation of the full two-dimensional equations with a topography and discharge corresponding to a stationary state for the pseudo-2D system, which enables some comparisons.

We consider a 200-m-long channel with a rectangular cross-section. The width and slope of the channel depend on x (Figure 9a). The inflow rate is fixed at q = 20 m 3 s -1 , and the water depth is prescribed at outflow (Delestre et al., 2013, §3.5.2). The Manning coefficient is set to 0.03 m -1/3 s. The channel is initially dry, with a small puddle downstream (because of the outflow condition).

At steady state, FullSWOF_2D produces a depth profile similar to the analytic solution (Figure 9b). An exact match is not to be expected because the results of a 2D code are being compared with the analytic solution of a pseudo-2D case.

Real Datasets

Finally, we provide the results from FullSWOF using three real datasets. The first case is an original laboratory experiment, the second case comes from an experimental plot in Senegal, and the third case is the well-known Malpasset dam break. We do not pretend to provide new insights into these phenomena; we merely demonstrate the ability of FullSWOF to simulate real situations without data filtering. The refined calibration of the parameters and a comparison with other pieces of software will be the topic of future works.

Flow over a Corrugated Bottom

This first example is an experimental flow over an inclined channel having a corrugated bottom (Figure 10). The discharge is imposed upstream (measured value: 0.69 l m -1 s -1 ). The water depths are measured at steady state along a 55-cm-long profile with a 0.5-mm resolution using the device described in [START_REF] Legout | High spatial resolution mapping of surface velocities and depths for shallow overland flow[END_REF]. FullSWOF_1D is run using the measured topography profile as the input. The only calibrated parameter is the Manning friction coefficient n = 0.0127 m -1/3 s. The comparison between the computed and measured values shows that FullSWOF_1D is able to reproduce the qualitative behavior after the first bump (x > 0.06 m) (Figure 10b). FullSWOF_1D reproduces all the hydraulic jumps (for x > 0.13 m) and correctly locates them at the minima of the measured solution. However, the simulated solution exhibits shocks that are steeper than the measured shocks. This, as well as the poor restitution at the beginning of the channel (at approximately x = 0.05 m), is likely due to a shortcoming of the model. Indeed, in this range of water depths and velocities, the surface tension likely cannot be neglected. This is a possible extension of the FullSWOF software. The purpose of this section is to test FullSWOF_2D using a real system: a plot in Thiès, Senegal [START_REF] Tatard | Measurement and modelling of high-resolution flow-velocity data under simulated rainfall on a low-slope sandy soil[END_REF][START_REF] Mügler | Comparison of roughness models to simulate overland flow and tracer transport experiments under simulated rainfall at plot scale[END_REF], an experimental system instrumented by IRD. This system consists of an artificial rainfall of approximately 2 hours with an intensity of approximately 70 mm h -1 on a sandy-soil plot with an area of 4 × 10 m 2 . The plot has the classical configuration of Wooding's open book, with a 1% slope along the Oxand Oy-axes. The complete dataset is freely available at http://www.umr-lisah.fr/Thies_2004/. FullSWOF_2D is used with a variable time step (the CFL value is fixed at 0.4) and with the same parameters as in [START_REF] Tatard | Measurement and modelling of high-resolution flow-velocity data under simulated rainfall on a low-slope sandy soil[END_REF]: f = 0.26, h f = 0.06 m, θ sθ i = 0.12, K s = 4.4 10 -6 m s -1 and K c = 0. Qualitatively, the results (Figure 11) are similar to those obtained with other software packages [START_REF] Tatard | Measurement and modelling of high-resolution flow-velocity data under simulated rainfall on a low-slope sandy soil[END_REF][START_REF] Mügler | Comparison of roughness models to simulate overland flow and tracer transport experiments under simulated rainfall at plot scale[END_REF]. At this stage, these results illustrate the ability of FullSWOF_2D to simulate a dynamic flow without any filtering of the data, in contrast to the previously utilized software packages. A comparison will be performed in a more detailed study. This section is devoted to a large-scale test case. We implement the wellknown Malpasset dam break, which occurred in 1959 in the south of France [START_REF] Hervouet | A high resolution 2-D dam-break model using parallelization[END_REF], Hervouet (2007, p.281-288)). Because of its varying topography and complex geometry, this is a classical test for numerical methods and hydraulics software validation. For this case, we used the parallel version of FullSWOF_2D [START_REF] Cordier | FullSWOF_Paral: Comparison of two parallelization strategies (MPI and SkelGIS) on a software designed for hydrology applications[END_REF] to reduce the computational time to a few hours. The dimensions of the computing domain are L x = 17273.9 m with 1000 cells and L y = 8381.3 m with 486 cells. The topography comes from a map published before the disaster that was subsequently digitized. We consider the Manning law with n = 0.033 m -1/3 s, as advised in [START_REF] Hervouet | A high resolution 2-D dam-break model using parallelization[END_REF]. No calibration was performed for this simulation, which uses a total simulation time of T = 4000 s.

The results are presented in Figure 12 a-d at four time steps. A scaled physical experiment was built by the Laboratoire National d'Hydraulique in 1964 to study the dam-break flow (for more details, see [START_REF] Hervouet | A high resolution 2-D dam-break model using parallelization[END_REF] and [START_REF] Hervouet | Hydrodynamics of free surface flows[END_REF]). The maximum water level was recorded at nine gauges (labeled 6 to 14) during the physical experiment (Figure 12e). The maximum water elevations of the scaled experiment are well captured by FullSWOF_2D (Figure 12f).

Conclusions

The FullSWOF software project grew out of a long-term collaboration between mathematicians and hydrologists in Orléans (France). Several specific features make it particularly suitable for applications in hydrology: FullSWOF considers rain, infiltration and classical friction laws and can operate with DEMs. Moreover, the software package has demonstrated its ability to address a wide range of flow conditions. The latest releases of the three codes (FullSWOF_1D, FullSWOF_2D and FullSWOF_UI) are fairly stable and robust. They are freely available for users outside of the original team for industrial, scientific and educational purposes.

The project is still undergoing several improvements. Apart from the usual bug correction work, several improvements to the current codes (e.g., implementing non-homogeneous friction) are regularly produced, and a parallelized 2D version is under development. More conceptual modifications are under consideration, such as the introduction of erosion and sedimentation models. Another direction of improvement concerns the surface tension. Although these modifications currently originate mainly from the initial development team, we emphasize that external contributions are welcome.

Figure 1 :

 1 Figure 1: Notations for 2D shallow-water equations.

  b) Geometry of the bi-layer model.

Figure 2 :

 2 Figure 2: Notations for the Green-Ampt infiltration model.

Figure 3 :

 3 Figure 3: Discretization of time and space in FullSWOF.

Figure 5 :

 5 Figure 5: Dam break on a dry domain without friction: comparison of the analytic solution with the FullSWOF_1D results. Simulation with 500 cells at T = 6 s.

Figure 6 :

 6 Figure 6: Planar surface rotating in a paraboloid: comparison of the analytic solution with the FullSWOF_2D results. Cross-sections at y = 2 m. Simulation with 100 cells in x and y at T = 13.4571 s.

Figure 7 :

 7 Figure 7: Short channel with a smooth transition and a shock: comparison of the analytic solution with the FullSWOF_1D results. Simulation with 500 cells at T = 1500 s.

Figure 8 :

 8 Figure 8: Rain on a long channel with a supercritical flow: comparison of the analytic solution with the FullSWOF_1D results. Simulation with 500 cells in x at T = 3000 s.

  Water depths and topography. For FullSWOF_2D, the mean depth inside the channel is plotted.

Figure 9 :

 9 Figure 9: Pseudo-2D channel with a supercritical flow: comparison of the analytic solution with the FullSWOF_2D results. Simulation with 400 cells in x and 201 cells in y at T = 200 s.

  (a) Close-up photography of the experimental setup. The red line is the measured profile.

  Measured topography, experimental water surface and water surface computed by FullSWOF_1D (Simulation with 1104 cells at T = 10 s).

Figure 10 :

 10 Figure 10: Experimental flow over a corrugated bottom.

  6.3.3 The Malpasset Dam Break: a Large Dataset Demanding the Parallel Version of FullSWOF_2D

Figure 11 :

 11 Figure 11: Results of FullSWOF_2D for the Thiès field plot. Simulation with 160 cells in x and 200 cells in y at T = 7400 s.

Figure 12 :

 12 Figure 12: Malpasset dam break: evolution of the water flow simulated by FullSWOF_2D (parallel version) and comparison with experimental data.

http://www.opentelemac.org/

University of Toulouse, CNRS, INSA. http://www.math.ups-tlse.fr/DassFlow

Numerical MethodsIn this section, we detail the numerical methods that are shown to be well adapted to finding solutions of the shallow-water system and that are implemented in FullSWOF. To obtain the solutions, the shallow-water system is divided into two parts: the transport (i.e., convective) operator and the

Acknowledgments

This study is part of the ANR METHODE granted by the French National Agency for Research, ANR-07-BLAN-0232. The authors would like to thank Olivier Planchon (IRD) for the data used in section 6.3.2 for the rain over a field plot (Thiès, Senegal).

http://www.univ-orleans.fr/mapmo/


Van-Chan Ngo

Jean-Pierre Talpin

Thierry Gautier

Loïc Besnard

Paul Le Guernic

Modular translation validation of a full-sized synchronous compiler using off-the-shelf verification tools (abstract)

The aim of this presentation is to demonstrate a scalable, modular, refinable methodology to design, assess and improve the trustability of an existing (20 years old), large (500k lines of C), open source (Eclipse/Polarsys IWG project POP) code generation suite using off-the-shelf, open-source, SAT/SMT verification tools (Yices), by adapting and optimizing the translation validation principle introduced by Pnueli et al. in 1998. This methodology results from the ANR project VERISYNC, in which we aimed at revisiting Pnueli's seminal work on translation validation using off-the-shelf, up-to-date, verification technology. In face of the enormous task at hand, the verification of a compiler infrastructure comprising around 500 000 lines of C code, we devised to narrow down and isolate the problem to the very data-structures manipulated by the infrastructure at the successive steps of code generation, in order to both optimize the whole verification process and make the implementation of a working prototype at all doable. Our presentation outlines the successive steps of this endeavour, from clock synthesis, static scheduling to target code production.

INTRODUCTION

Synchronous programming languages like Signal, Lustre and Esterel propose a formal semantic framework to give high-level specification of safety-critical software in automotive and avionics systems [START_REF] Berry | The foundations of Esterel[END_REF][START_REF] Le Guernic | Programming real-time applications with Signal[END_REF][START_REF] Halbwachs | A synchronous language at work: the story of lustre[END_REF][START_REF] Jackson | A direct path to dependable software[END_REF]. Safety-critical systems are those systems whose failure could result in loss of life, or damage to the environment. They need to be validated to ensure that their specified safety properties are implemented correctly. Software validation is traditionally done by using testing techniques which, in the case of safety-critical systems, is not sufficient [START_REF] Moy | Testing or formal verification: Do-178c alternatives and industrial experience[END_REF]. Since synchronous languages are based on formal semantic models, they provide much higher level of abstraction, expressivity, and clarity at source level rather than once compiled into C code. That makes the application of formal methods much simpler to enforce safety properties. However, a synchronous compiler is still a large and complex program which often consists of hundreds of thousands lines of code, divided into numerous packages. Moreover, compiler modules often interact in sophisticated ways, and the design and implementation of a compiler are substantial engineering tasks. The compilation process involves many analyzes, program transformations and optimizations, some may introduce additional information or constrain the compiled program, some may refine its meaning and specialize its behavior to meet a specific safety or optimization goal, and all these compile-time decisions should additionally be traced.

STATE OF THE ART

Proving the correctness of a compiler can be based on the examination of the developed compiler's source code itself, meaning that a qualification process applies on the development of the compiler, the source of the compiler, and/or the compiler's output. Qualifying a compiler is rare because of the tremendous administrative effort involved. Qualification amounts to demonstrate compliance with all recommendations and objectives specified in the certification standards for safety-critical softwares: Do-178C and its European equivalent Ed-12 [START_REF]Software considerations in airborne systems and equipment certification[END_REF]. Although Do-178 has been successful in industry, the cost of complying with it is significant: the activities on verification it incurs may well cost seven times more than the development effort needed [START_REF] Rushby | New challenges in certification for aircraft software[END_REF]. A more traditional method is therefore to solely inspect or formally verify the compiler's output. This task requires less unitary effort, but has to be repeated every time target code is generated. For instance, Astrée [START_REF]Astrée[END_REF] [3] is a special-purpose static program analyzer based on abstract interpretation to verify the absence of run time errors in the C code generated from Scade programs. One last resort is to formally verify the correctness of the compiler itself. Certifying compilation [START_REF] Necula | Proof-carrying code[END_REF] attests that the generated object code satisfies the properties established on the source program by generating concrete evidences along the compilation into object code. Systematic compiler verification techniques use formal methods. For the purpose of compiler verification, there are two approaches to prove the software correctness.

Formal compiler verification consists of specifying the behavior of the compiler in a formal specification language and build a proof that the compiler satisfies behavioral equivalence or refinement. Formal verification can be done through many approaches. One such approach is deductive verification. It consists of providing deductive proofs that a system behaves in a certain way that is described in the specification, with the aid of either interactive theorem provers (such as Hol [START_REF] Gordon | Introduction to HOL: A theorem proving environment for higher order logic[END_REF], Isabelle [START_REF]The Isabelle proof assistant[END_REF], or Coq [START_REF]Coq proof assistant[END_REF]), or an automated theorem prover. Another approach is model checking [START_REF] Clarke | Automatic verification of finite-state concurrent systems using temporal logic specifications[END_REF][START_REF] Quielle | Specification and verification of concurrent systems in Cesar[END_REF]. It involves building an abstract model of the system and ensure it complies with specified requirements by exploring all its accessible states. Requirements are represented in temporal logics, such as Linear Temporal Logic (Ltl) or Computational Tree Logic (Ctl) and verification produces a confirmation that the system model conforms to requirements or a counterexample that can be used to locate and eliminate an error. Some techniques can be used to deal with the state explosion problem including abstract interpretation, symbolic simulation and abstract refinement [START_REF] Cousot | Refining model checking by abstract interpretation[END_REF]. A variant of model checking, Bounded model checking (Bmc) [START_REF] Clarke | Bounded model checking using satisfiability solving[END_REF], encodes the fact that potential executions of the system model do not conform to the specification in incremental fashion as propositional satisfiability formulas. The bounded number of evaluation steps is increased as long as the resulting propositional formula is satisfiable. Then a concrete counterexample can be extracted as a trace of system states leading to an error state in the system model. Inductive reasoning can also be used to prove that a system conforms its specification. Advances in Satisfiability Modulo Theories (Smt) have been useful for checking systems inductively. With SMT solvers, systems can be modelled efficiently, require fewer limitations to represent specifications and meet significant performances. [START_REF] Pnueli | Translation Validation[END_REF] as an approach to verify the correctness of translators (compilers, code generators). The main idea of translation validation is that instead of proving the correctness of a compiler, each of its individual translations (e.g., run of the code generator) is followed by a validation phase to check that the target program correctly implements the source specification. A translation validator consists of two objects:

• The Model builder is a simple module that formally represents the semantics of the source and target programs of the translator (e.g., labeled transition system, first-order logic formula).

• The Analyzer formalizes correctness as a refinement relation between the models of the source and target of the validator. The analyzer provides an automated proof method on the existence of the refinement between the formal models. If the analyzer succeeds, a proof script is created. If it doesn't, it generates a counter example, which can be decompiled to help spot the error. Translation validation does not modify or instrument the compiler. It treats it as a "black box" (as long as there is no error in it). It only considers the input program and its compiled result. Hence, it is not affected by updates and modifications to the compiler, as long as its data-structures remain the same. In general, the validator is much simpler and smaller than the compiler itself. Thus, the proof of correctness of the validator takes less effort than that of the compiler. Plus, verification is fully automated and scales to large programs.

SCALABLE TRANSLATION VALIDATION

Our approach is to scale translation validation not only in a modular fashion, by decomposing the problem into the successive transformations performed by the compiler on the intermediate representation of a program [START_REF] Ngo | Formal Verification of Compiler Transformations on Polychronous Equations[END_REF], but by narrowing it further to the actual datastructure that are being used to represent the transformation problem and the actual algorithmic operations performed on it. In all cases, we show that translation validation is amenable to simple SAT/SMT verification (the semantic inclusion of one data-structure into another) instead of the more general problem of simulation-based conformance-checking of the transformed program w.r.t. the input specification [START_REF] Ngo | Translation Validation for Synchronous Data-flow Specification in the SIGNAL Compiler[END_REF][START_REF] Ngo | Translation validation for clock transformations in a synchronous compiler[END_REF][START_REF] Ngo | Efficient deadlock detection for polychronous data-flow specifications[END_REF]. In the case of the synchronous data-flow language Signal, a 500k-lines big code generation infrastructure, the compilation process can be divided into three phases depicted in Figure 1 (top row). A source Signal program is the synchronous composition of discrete equations on signals, e.g. x :" y `1 | y " x$1 defines x by y `1 at all times and y by the value of x delayed by 1 evaluation tick: its previous value. Its compilation may be seen as a sequence of morphisms that refine and rewrite the source specification with information gained from analysis. C or Java code production is performed on ultimately transformed program, e.g., y " x; x :" y `1.

• Clock calculation. This stage determines the clock of all signals in the program and defines a Boolean abstraction of the program. The clock of a signal defines when the value of the signal shall be evaluated. A proof of semantic-preservation can be decomposed into the preservation of clock semantics at the clock calculation phase and that of data dependencies at the static scheduling phase. Valueequivalence of variables can be then checked at the code generation phase. Figure 1 shows the integration of this verification framework into the compilation process of the Signal compiler. For each phase, the validator takes the source program and its compiled counterpart, and constructs the corresponding formal model of the program. Then, it checks the existence of the refinement relation to prove semantic-preservation.

Preservation of clock models [24]

The first verification stage focuses on proving that all clock relations associated with signals in the source and transformed program are equivalent. A clock model is a first-order logic formula that characterizes the presence/absence status of all signals at all times in a Signal program. Given two clock models, a clock refinement relation is defined to express the semantic preservation of clock semantics. The existence of a clock refinement is defined as a satisfiability problem which is automatically and efficiently proved by an Smt solver.

Example Consider the Signal program DEC which counts through the output N from the value of input FB to 1, ZN := N$1 init 1 defined ZN as the previous value of N, N := FB default (ZN-1) assigns FB to N when the input FB is present, and ZN-1 otherwise, FB ^= when (ZN<=1) synchronizes FB to the condition (ZN<=1). The clock model of the source program is:

ΦpDECq " p x FB ô y ZN 1 ^Ą ZN 1 q ^p y ZN 1 ô y v 1 ă" ô x ZNq ^p y ZN 1 ñ p Ą ZN 1 " Ą v 1 ă" qq ^p x ZN ô p Nq ^p x ZN ñ p Ă ZN " m.N ^m.N 1 " Ñqq ^pm.N 0 " 1q ^p p N ô x FB_ y ZN 2 q^p p N ñ pp x FB^Ñ " FBq_p x FB^Ñ " Ą ZN 2 qqq ^p y ZN 2 ô x v 1 ´ô x ZNq ^p y ZN 2 ñ p Ą ZN 2 " Ă v 1
´qq Terms x, resp. x, represent the clock resp. value of a signal x.

The model of the transformed program DEC' is twice as large.

Checking it a correct refinement amounts to establish a variable mapping y X DEC z y X IO " αp z X DEC 1 z y X IO q between DEC and DEC' and delegate the checking validity of the formula pΦpDEC 1 q^y X DEC z y X IO " αp { X DEC_BASIC_TRA z y X IO q ñ ΦpDECqq, named ϕ, to the Smt solver under the logical context defined by the variable mapping and the following assertions:

p Ă ZN " Ą ZN c ^1 " 1q ñ ´p Ą v 1 ă" " Ą v 1c ă" q ^pĂ v 1 ´" Ă v 1c ´q3
.2 Preservation of data dependency [START_REF] Ngo | Efficient deadlock detection for polychronous data-flow specifications[END_REF] The goal of this stage is to prove that the existence of a data dependency between two signals in the source specification is a property preserved by the target program which, in addition, makes a sequential schedule of computations explicit.

d ZN1 = d ZN2 = c ZN = b N c FB = d ZN1 ^g ZN1 b N = c FB _ d ZN2 c FB FB N d ZN1 ZN ZN1 ZN2 c FB ¬ c FB ^d ZN2 c FB c ZN d ZN1 d ZN2 d ZN2 d ZN1 c FB c FB

Translation validation of SDDG

In this section, we adopt the translation validation approach to prove the correctness of the compiler in the static scheduling phase of the compilation process. Given two SDDGs, we first formalize the notion of "correct implementation" as a dependency refinement relation. This refinement expresses the semantic preservation of data dependencies. Then, we propose a method to implement our verification framework by the use of a SMT solver for checking the existence of the above refinement relation. Next, it finds all dependency paths from the signal FB to the signal N, and the cycles in the graphs SDDG(DEC) and SDDG(DEC_SEQ_TRA) to generate the formulas for checking the dependency refinement. In this example, we get the following formula:

FB1 c N [ FB1 c ZN3 c C_FB c ZN1 c [ ZN3 c ZN2 c ZN c C_CLK_12 c [ FB1 c ¬ [ FB1 c ^[ ZN3 c [ FB1 c [ FB1 c [ FB1 c \ C_FB c [ ZN1 c [ ZN3 c [ ZN2 c [ ZN3 c [ ZN3 c \ C_CLK_12 c
c FB ) ( [ FB1 c ^[ FB1 c )
In the third step, we delegate the checking validity of the formula to YICES solver. We get unsat when checking the satisfiability of

¬( c FB ) ( [ FB1 c ^[ FB1 c
)), which means the formula is valid. Therefore, we can conclude that SDDG(DEC_SEQ_TRA) vdep SDDG(DEC) Checking formula construction establishes the variable mapping much like for clock models. Next, it amounts to finding all cycles and dependency paths from FB to N to generate the formulas for checking the dependency refinement, the validity of which is delegated to the Yices solver.

x FB ñ p z FB1 c ^z FB1 c q 3.3 Value-equivalence of variables [START_REF] Ngo | Translation Validation for Synchronous Data-flow Specification in the SIGNAL Compiler[END_REF] This else N = N -1; w_DEC_N(N); DEC_step_finalize(); return TRUE; }

SDVG translation validation 159

In the first step, we shall compute the shared value-graph for both programs to represent the computation of all signals and their corresponding variables. This graph is depicted in Figure 6.34. Listing 6.11 Generated C code of DEC Note that in the C program, the variable N c ("c" is added as superscript for the C program variables, to distinguish them from the signals in the SIGNAL program) is always updated (line ( 8)). In lines ( 3) and ( 8), the references to the variable N c are the references to the last value of N c denoted by m.N c . The variable FB c which corresponds to the input signal FB is updated only when the variable C_FB c is true.

{ N, c ZN} _ { c FB} ^<= { f ZN} g m.N 1 {N} f { Ñ} f f FB ? {ZN} f {FB} f m.N c {C_FB c } <= {N c , f N c } f g FB c {FB c } f
The shared Sdvg are then normalized by employing a series of rewriting rules to merge nodes refering to the same value. This yields the following graph, proving that the generated step function conforms its specification 6.5 Discussion 161 conclude that the output signal N and its corresponding variable N c is equivalent since they point to the same node in the final normalized graph.

{ g m.N, ZN, f ZN} m.N c {C_FB c , c FB} <= 1 {N c , f N c , N, e N} f { g FB c } f FB ? {FB c , FB} f

Discussion

There is a wide range of works for value-graph representations of expression evaluations in a program. For example, in [? ], Weise et al. present a nice summary of the various types of value-graph. In our context, the value-graph is used to represent the computation of variables in both source program and its generated C code and identical structures are shared. We believe that this representation will reduce the required storage and make the normalizing process more efficient. Another remark is that the calculation of clocks as well as the special value, the absent value, are also represented in the shared graph.

Another related work which adopts the translation validation approach in verification of optimizations, Tristan et al. [? ], recently proposed a framework for translation validation of LLVM optimizer. For a function and its optimized counterpart, they compute a shared value-graph. The graph is normalized (the graph is reduced). After the normalizing, if the outputs of two functions are represented by the same sub-graph, they can safely conclude that both functions are equivalent.

On the other hand, Tate et al. [? ] proposed a framework for translation validation. Given a function in the input program and the corresponding optimized version of the function in the output program, they compute two value-graphs to represent the computations of the variables. Then they transform the graph by adding equivalent terms through a process called equality saturation. After the saturation, if both value-graphs are the same, they can conclude that the return value of two given functions are the same. However, for the translation validation purpose, our normalization process is more efficient and scalable since we

RELATED WORKS

The notion of translation validation was introduced in [START_REF] Pnueli | Translation Validation[END_REF][START_REF] Pnueli | Translation Validation: From Signal to C. Correct System Design Recent Insights and Advances[END_REF] by Pnueli et al., using a symbolic model of Synchronous Transition Systems (Sts) to represent both source and target programs. An Sts is a set of logic formulas which describes the functional and temporal constraints of the whole program and its generated C code. Bdds [START_REF] Bryant | Graph-based Algorithms for Boolean Function Manipulation[END_REF] implement the symbolic Sts models. Our approach improves over standard translation validation by explicitly capturing the clock semantics in the model, which also results in much smaller models. Another related work is the static analysis of Signal programs for efficient code generation [START_REF] Gamatié | Static Analysis of Synchronous Programs in Signal for Efficient Design of Multi-Clocked Embedded Systems[END_REF], where linear relations among clocks and values are represented by first-order logic formulas with the help of numerical interval abstraction techniques. The objective is to make generated code more efficient by detecting and removing the dead-code segments (e.g., segment of code to compute a data-flow which is always absent). They determine the existence of empty clocks, mutual exclusion of two or more clocks, or clock inclusions, by reasoning on the formal model using an Smt solver. Related works have also adopted the translation validation approach in verification of transformations, and optimizations. In [START_REF] Necula | Translation Validation for an Optimizing Compiler[END_REF], the translation validation is used to verify several common optimizations such as common subexpression elimination, register allocation, and loop inversion. The validator checks the existence of a simulation relation between two programs. Leroy et al. [START_REF]The CompCert Project[END_REF][START_REF] Leroy | Formal Certification of a Compiler Back-end, or Programming a Compiler with a Proof Assistant[END_REF] used this technique to develop the CompCert high-assurance C compiler. The programs before and after the transformations and optimizations of the compiler are represented in a common intermediate form, then the preservation of semantics is checked by using symbolic execution in the proof assistant Coq. Tristan et al. [START_REF] Tristan | Evaluating Value-graph Translation Validation for LLVM[END_REF] recently proposed a framework for translation validation of Llvm optimizer. For a function and its optimized counterpart, they compute a shared value-graph. The graph is normalized. If the outputs of two functions are represented by the same sub-graph, they can safely conclude that two functions are equivalent. We believe that our approach is more modular and efficient in both design time, space and time than these based on proof automation and simulation relations: by reducing each of the translation steps to these of the very data-structures subject to refinements (the clock hierarchy, the data-flow graph, the value graphs) considerably reduces the size of the refinement or simulation problem to solve, and that using off-the-shelf verifiers like Z3, Yices, SMTLib guarantees both speed and correctness.

CONCLUSION

We have presented a technique based on Smt solving to prove the preservation of clock semantics during the compilation of a synchronous data-flow compiler. Our approach focuses on the transformations performed by the compiler using the simplest structures to represent them: SAT/SMT formulas represent the refinement of clock models and the reinforcement of data-flow graphs, value graphs are used to represent the production of target code patterns from and a specification's syntax tree. This reduces the whole process of proving a refinement relation between the source specification and the generated code to a couple of SMT SATchecking on formulas of minimal size, and to a symbolic rewriting on a reduced graph to check value equivalence. Our validator does not modify or instrument the compiler. It treats it as a "black box" (as long as there is no error in it). It only considers an input program and its transformed result. Hence, it is not affected by an update an modification made to this or that compilation stage, as long as its principle and data-structure remains the same. The validator is much simpler and smaller than the compiler itself. Proving its correctness (the model builder, the verifier) would take a lot less effort than for the compiler as well. Verification is fully automated and scales to large programs very well by employing state-of-theart verification tools and by minimizing the representation of the problem to solve. For that purpose, we represent the desired program semantics using a scalable abstraction and we use efficient Smt libraries [START_REF]The Yices Sat-solver[END_REF] to achieve the expected goals: traceability and formal evidence. We believe that this approach provides a, both technically and economically, attractive alternative to developing a certified compiler. The individual modules designed in the context of this project are being integrated in the open-source environment of the Eclipse project POP with the Polarsys Industry Working Group [START_REF]Polychrony On Polarsys[END_REF].

Figure 1

 1 Figure 1 The translation validation for the Signal compiler

Fig. 5 . 5

 55 Fig. 5.5 The SDDG of DEC

5. 3 . 1

 31 Definition of correct implementation: Dependency refinement Let SDDG(A) = hNA, EA, IA, OA,CA, mN A , mE A i and SDDG(C) = hNC, EC, IC,OC,CC, mN C , mE C i, to which we refer respectively as the data dependency representations of the source program and its transformed counterpart produced by the SIGNAL compiler. Let x and y be two signals in both programs A and C. A dependency path from the signal x to the signal y in In the first step, we shall construct the synchronous data-flow dependency graphs which are depicted in Figure 5.5 and Figure 5.7. Two fresh variables FB1 and ZN3 are used to replace the expressions FB when C_FB and ZN2 when C_CLK_12. We omit the dependencies among the signals FB, c FB,C_FB c and \ C_FB c in the graph of the program DEC_SEQ_TRA.

FB

Fig. 5 .

 5 Fig. 5.7 The SDDG of DEC_SEQ_TRA In the second step, checking formula construction, our tool will establish the variable mapping d XDEC \ d XIO = a( \ XDEC_BASIC_TRA \ d XIO) as in Section 4.4.5.1.Next, it finds all dependency paths from the signal FB to the signal N, and the cycles in the graphs SDDG(DEC) and SDDG(DEC_SEQ_TRA) to generate the formulas for checking the dependency refinement. In this example, we get the following formula:

Fig. 6 .

 6 Fig. 6.34 The shared value-graph of DEC and DEC_step

Fig. 6 .

 6 Fig. 6.36 The final normalized graph of DEC and DEC_step

  Sequential C or Java code is directly generated from the structure of the clocked and scheduled Signal program.

• Static scheduling. Based on the clock information and the Boolean abstraction obtained at the first stage, the compiler constructs a Conditional Dependency Graph (Cdg), which represents the schedule of signals' evaluations. • Code generation.

  Along the way, the validator further checks the target program deadlock-free. Data dependencies among signals are represented by a Synchronous Data-flow Dependency Graph (Sddg). An Sddg is a labeled directed graph in which node are signals and clocks and edges represent dependencies between nodes. Edge are labeled by clocks. An edge clock tells when the dependency is effective: when its clock is present. The correctness of a schedule is formalized as a dependency refinement relation between the source and target Sddgs. It is implemented by Smt-checking the existence of the refinement relation.

	94	Translation validation of SDDG
	graph are encoded as the first-order formulas given as follows.
		Example The Sddgs of DEC and the scheduled DEC' are the following

  stage focuses on proving that every output signal in the source program and the corresponding variable in the generated C program are assigned the same values at all times. The defining equation of a signals and its C translation are represented by a shared valuegraphs, called Synchronous Data-flow Value-Graph (Sdvg). To prove that a signal and a variable have the same value with an Sdvg, we just need to check that they are represented by the same sub-graph, meaning that they point to the same graph node.Example The code generated from DEC consists of the following step function (left). The value graphs of DEC and DEC_step are constructed and composed (right).

	logical DEC_step() {
	C_FB = N <= 1;
	if (C_FB)
	if (!r_DEC_FB(&FB))
	return FALSE;
	if (C_FB) N = FB;