Marc Honnorat
email: marc.honnorat@imag.fr>

Joel Marin
email: joel.marin@imag.fr>

Jerome Monnier
email: jerome.monnier@imag.fr>

Xijun Lai
email: laixijin@gmail.com>

Marc Honnorat -Joël

Marin -Jérôme Monnier -Xijun Lai

Joël Marin

Jérôme Monnier

DassFlow v1.0: a variational data assimilation software for 2D river flows

Keywords: Variational data assimilation, parameter identification, river hydraulics, flood Assimilation variationnelle de données, identification de paramètres, hydraulique fluviale, inondations

come

Introduction

The St-Venant equations (shallow water) can describe accurately river flows, in the main channel (1D geometry) or in the major bed (flood plain, 2D geometry). Nevertheless, in order to carry out a reliable numerical simulation, models require input parameters such as bed elevation, roughness coefficients in addition to initial and boundary conditions. Generally, these parameters are approximated (eg roughness coefficients), missing (eg inflow boundary conditions or initial conditions) or subjected to uncertainties. Moreover, the governing equations themselves may model ideal configurations, which not necessarily correspond to real configurations. Classicaly, expert-users of computational river flow softwares, see e.g. [START_REF] Edf | Telemac : numerical modelling system for free surface hydrodynamics, hydroinformatics[END_REF], must perform many simulations (hence spend time) in order to calibrate the numerical model related to a given configuration (trial-error tests). On the other hand, observations from in-situ gauge stations (eg water elevation) or remote sensing observations (satellite images, air photographs, video images) are very important sources of information on the flow state, but they are not qualitatively integrated into the simulation process. To improve the quality of the simulation, data assimilation methods combine in an optimal sense the information from the model and observated data. This allows to identify some parameter values consistent with reality, and /or initial conditions. Variational data assimilation method, see e.g. [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF], [START_REF] Bouttier | Data assimilation concepts and methods[END_REF], is based on the optimal control theory, [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF], and it aims to fuse the dynamical model and observations by minimizing a cost function which measures the discrepancy between the simulation results and physical measurements. This method is operational in meterology and since more recently in oceanography. In river hydraulics, variational data assimilation methods have been used for the identification of model parameters in 1D channels, [START_REF] Khatibi | Identification problem of open-channel friction parameters[END_REF][START_REF] Roux | Data assimilation applied to hydraulic parameter identification[END_REF], and in 2D, [START_REF] Atanov | Estimation of roughness profile in trapezoidal open channels[END_REF][START_REF] Yang | Assimilation de données variationnelle pour les problèmes de transport des sédiments en rivière[END_REF][START_REF] Mazauric | Assimilation de données pour les modèles d'hydraulique fluviale. Estimation de paramètres, analyse de sensibilité et décomposition de domaine[END_REF][START_REF] Ding | Identification of Manning's roughness coefficients in shallow water flows[END_REF][START_REF] Dimet | Prospects for the use of data assimilation for flood prediction[END_REF][START_REF] Bélanger | Data assimilation (4D-VAR) to forecast flood in shallowwaters with sediment erosion[END_REF][START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF][START_REF] Lai | Assimilation of spatial distributed water levels into a shallowwater flood model. part i: Method and toy test case[END_REF]. In real configurations, observations are available only in very small quantities. At best, water levels are measured at very few gauge stations (information very sparse in space); velocity measurements are even more scarce since they require complex human interventions.

DassFlow is a river hydraulics simulation software designed for variational data assimilation on simple hydraulic configurations for research purpose. The model is based on the bidimensional shallow-water equations, solved by the finite volume method using the HLLC approximate Riemann solver. The minimization procedure is based on the L-BFGS algorithm, [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF]), which requires the computation of the gradient of the cost function. Since the forward code is written in Fortran 95, the adjoint code is generated by the automatic differentiation tool Tapenade [28] developped by the Tropics project-team at INRIA Sophia-Antipolis. DassFlow software is interfaced with few free and commercial pre and post-processors (SIG tools, mesh generators, visualization tools), which allows to performs computations with real data.

In Section 2, we present the Shallow Water model. In section 3, we present the Variational Data Assimilation method (4D-Var). In Section 4, we detail the finite volume scheme, then the way the adjoint code is automatically generated. In Section 5, we compare some solutions of the forward model with either explicit solutions or real measures (benchmarks INRIA for the forward code). In Section 6, we consider twin experiments which show the capabilities of DassFlow (benchmarks for the data assimilation software). In appendix we detail characteristics boundary conditions and Riemann invariants.

The present research report treats of the mathematical and numerical features of Dassflow only. Concerning the code itself, we refer to the user manual [START_REF] Honnorat | DASSFLOW : USER GUIDE[END_REF] and the developer guide [START_REF] Honnorat | DASSFLOW : DEVELOPER GUIDE[END_REF].

2 The shallow water model (St Venant)

The equations

The river hydraulics model is based on the bidimensional shallow water equations in their conservative formulation. The state variables are the water depth h and the local discharge q = hu, where u = (u, v) T is the depth-averaged velocity vector. On a domain Ω and for a computational time interval [0, T ], the equations are:

           ∂ t h + div( q ) = 0 in Ω×] 0, T ] ∂ t q + div 1 h q ⊗ q + 1 2 g∇h 2 + gh∇z b + g n 2 q 2 h 7/3 q = 0 in Ω×] 0, T ] h(0) = h 0 , q(0) = q 0 + boundary conditions (see section 2.2) ( 1
)
where g is the magnitude of the gravity, z b the bed elevation, n the Manning roughness coefficient, h 0 and q 0 the initial conditions for the state variables. The quantity c = √ gh denotes the local wave celerity. It is important to notice that if the bed elevation is constant (∇z b = 0) and without the friction term, (1) is an hyperbolic system. Moreover, a transport equation related to a tracer Φ is considered:

∂ t hΦ + div( qΦ ) = s f (2)
with initial and boundary conditions. s f is the tracer source term.

Boundary conditions

We split the boundary Γ of the domain Ω in three parts (see figure 2.1). We denote by: Γ in the part where the flow is incoming (type INFLOW)

Γ out the part where the flow is outgoing (type OUTFLOW)

Γ wall the part where it is a closed boundary (type WALL) The user has a large choice of boundary conditions, depending on the boundary condition type.

Wall boundary conditions

On Γ wall , we consider a slip condition:

u • n | Γc (t) = 0 , ∂h ∂n | Γc (t) = 0 ∀ t ∈ ]0, T ]

Inflow conditions

Few boundary conditions are possible at inflow ie on Γ in .

Discharge imposed

We have:

(q • n) | Γin = -q in (t) , ∂h ∂n | Γin (t) = 0 ∀ t ∈ ]0, T ]

INRIA

Incoming characteristics

The characteristic variables of the system are:

w 1 = u • n + g h0 h, w 2 = u • τ and w 3 = u•n-g h0 h, associated to the eigenvalues λ 1 = u 0 •n+c, λ 2 = u 0 •τ and λ 3 = u 0 •n-c respectively. Let us recall that: c = √ gh 0 .
We detail more concerning the characteristics in appendix A. The characteristic w k is incoming if λ k < 0 and outgoing otherwise. The characteristics boundary conditions are implemented as follows:

. For each point of the boundary, we compute the three eigenvalues λ k .

. For each λ k , -if λ k ≥ 0, then w k is incoming; thus it is evaluated from external data, -if λ k < 0, then w k is outgoing; thus it is evaluated from the internal simulation domain.

. Then, the model variables can be retrieved as follows:

-u • n bnd = w1+w3 2 -u × n bnd = w 2 -h bnd = h0 g w1-w3 2

Outflow conditions

Few boundary conditions are possible at outflow ie on Γ out .

Homogeneous Neumann

∂h ∂n | Γout (t) = 0 , ∂q ∂n | Γout (t) = 0 ∀ t ∈ ]0, T ] Water elevation prescribed h | Γout (t) = z out (t) -z b | Γout , ∂(u•n+2c) ∂n | Γout (t) = 0 ∀ t ∈ ]0, T ]
This boundary condition should be applied only in case of sub-critical flow. We refer to appendix B for a justification of the relation

∂(u•n+2c) ∂n | Γout (t) = 0 RR n 6150

Incoming characteristics

This case is similar to those related to the inflow boundary Γ in .

Rating curve

A rating curve is a relation between the normal discharge q • n and the water height h : q • n = f (h). Such a relation can be used to specify outflow boundary conditions if combined with the relation:

∂(U+2 √ gh) ∂n = 0.
To this end, we solve the following non-linear system (see figure 2.2):

q • n = f (h) u • n + 2 √ gh = u • n + 2 √ gh int (3)
where u • n + 2 √ gh int is computed from the interior of the domain. The DassFlow software is designed for variational data assimilation (classicaly called 4D-var method). Variational data assimilation is based on the optimal control theory, see [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF], [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]. Briefly, this method consists to compute a control vector value k minimizing a cost function which measures the discrepancy between the computed variable and available data.

In our case the control vector k can include the initial condition, boundary conditions, manning coefficient and bed elevation. For example, if inflow discharge is prescribed and water elevation is prescribed, we get: k = (h 0 , q 0 , q in , z out , n, z b ) T

INRIA

If we impose incoming characteristics for both inflow and outflow boundaries, we get:

k = (h 0 , q 0 , [w in 1 , w in 2 , w in 3 ], [w out 1 , w out 2 , w out 3 ], z out , n, z b ) T
We set below the observation function H(k; h, q) and the cost function J(k) = H(k; h, q) which measures the discrepancy between the state of the system (computed variable) and observations. Then the optimal control problem we solve reads:

M in k J(k) (4)
where (h, q) is the solution of the forward model [START_REF] Atanov | Estimation of roughness profile in trapezoidal open channels[END_REF].

This optimization problem is solved numerically by a descent algorithm. Thus, we need to compute the gradient of the cost function. This is classicaly done by introducing the adjoint model.

Let us point out that initial conditions, boundary conditions, manning coefficients and bed elevation are only potential control variables. In practice, one manages to identify very few of them only at the same time.

Observations and cost function

Given observations, we define the cost function J which measures the discrepancy between the computed variable (state of the system) and the available data.

In the present version, we assume that water depth and discharge are available at some point and some time instant. In the forthcoming version of DassFlow, we will assume to have some extra observations such as trajectories at surface (lagrangian data).

Given such (eulerian) observations, we define the cost function as the sum of three terms.

The first term measures the discrepancy between observation data and computed state variable :

J obs (k) = 1 2 T 0 C h h(k, t) -h obs (t) 2 + T 0 C q q(k, t) -q obs (t) 2
where C h and C q are restriction operators.

The second term measures the discrepancy of the flux: The third term is a classical regularization term. For example, if we control inflow discharge, we set:

J f lux (k) = 1 2 T 0 G(k, t) -G obs (k, t) 2 RR n 6150 with G(k, t) -G obs (k, t) = ∆t |K i | 4 j=1 T -1 ij G 1 (U L , U R ) -G 1 ( ŨL , U R ) with ŨL = [h obs , h obs u, h obs v] T and G 1 (U L , U R )
J p (k) = 1 2 ∂ 2 tt q in 2
Finally, the cost function writes :

J(k) = α obs J obs (k) + α f lux J f lux (k) + α p J p (k) (5)
where α obs , α f lux and α p are scaling coefficients which must be setted by user.

Adjoint model

In order to compute efficiently all partial derivatives of the cost function J(k) with respect to the components of the control vector k, we introduce the adjoint model, [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]. [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF].

INRIA

Given the state of the system (h, q), solution of the forward model, we compute (h * , q * ), solution of the adjoint model:

                             ∂ t h * -1 h 2 (q • ∇) q * • q + gh div(q * ) -g q * • ∇z b + 7 3 g n 2 q h 7/3 q • q * = α obs C T h C h h -h obs + α f lux D k (G -G obs ) T (G(k, t) -G obs (k, t)) ∂ t q * + ∇h * + 1 h (q • ∇) q * + 1 h (∇q * ) T q -g n 2 q h 7/3 q * -g n 2 h 7/3 q (q⊗q)q * = α obs C T q C q q -q obs h * (T ) = 0 , q * (T ) = 0 . (6)
Let us point out that the adjoint model is reverse in time, and it includes the observations. Concerning its boundary conditions, for example in the case of inflow discharge prescribed at inflow and water elevation prescribed at outflow, we get:

C.L. : q * | Γin = 0 , (q * • n) | Γ wall = 0 , (q * • τ ) | Γout = 0 , (∂ n h * ) | Γin∪Γ wall = 0 , h * + 2(u • n)(q * • n) | Γout = 0 . (7)
Given the state (h, q) and the adjoint state (h * , q * ), one obtains the following expression of the partial derivatives of the cost function (independent of the derivative of the state):

∂j ∂h 0 (k) = -h * (0) , ∂j ∂q 0 (k) = -q * (0) , ∂j ∂z b (k) = - T 0 div gh(t)q * (t) dt , ∂j ∂n (k) = 2 g n T 0 u(t) h(t) -1 3 u(t), q * (t) Ω dt , ∂j ∂q in (k) = α p ∂ 2 tt T ∂ 2 tt q in -h * | Γq , ∂j ∂z out (k) = (q * • n) c 2 -(u • n) 2 | Γz . (8)
A single integration of the direct model ( 1) followed by a single integration backward in time of the adjoint model ( 6)-( 7) allow to compute all components of the gradient of the cost function (see figure 3.4).

The optimal control problem ( 4) is solved using a local descent algorithm, more precisely the L-BFGS algorithm (a quasi-Newton method), implemented in the M1QN3 routine, see [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF][START_REF] Gilbert | M1qn3 web site[END_REF]. Thus, these partial derivatives are used as input to the minimization algorithm M1QN3. The global optimization process is represented in Fig. 3.4.

In practice, we don't implement the adjoint model presented above but we compute the partial derivatives directly by differentiating the forward code using an automatic differentiation tool, see next section.

Forward code

Adjoint code

Numerical implementation

In this section, we describe on one hand the finite volume scheme solving the forward model and its implementation, on the other hand the automatic differentiation of the forward code in order to obtain the partial derivatives of J. For more details concerning the numerical implementation of the finite volume scheme, one can consult [START_REF] Honnorat | Assimilation de données Lagrangiennes pour la simulation numérique en hydraulique fluviale[END_REF], [START_REF] Honnorat | DASSFLOW : DEVELOPER GUIDE[END_REF].

INRIA

Forward code and numerical scheme

The forward (or direct) model ( 1) is discretized using by a finite volume method based on the approximate Riemann solver HLLC, see e.g. [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF]. The 2D mesh can be a mix of triangular and quadrangular cells.

The finite volume scheme

Equations ( 1) are rewritten as follows :

∂ t U + ∂ x (G(U )) + ∂ y (H(U )) -(S g (U ) + S f (U )) = 0, (x, y) ∈ Ω 2 , t ∈ (0, T ) where U = [h, q x , q y ] T , G(U ) = [q x , q 2 x /h + gh 2 /2, q x q y /h] T H(U ) = [q y , q x q y /h, q 2 y /h + gh 2 /2] T S g (U ) = [0, gh∂ x z b , gH∂ y z b ] T
z b is the topography term, S f is the 2D friction term (Manning law) and g is the magnitude of the gravity. For the finite volume K i of the mesh, we define the mean value of the state variable U :

U i = 1 |K i | Ki U dΩ , (9)
where |K i | denotes the surface of the cell. We integrate the equation over K i , we use the divergence theorem, we use the rotational invariance property of the equations [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF], then we obtain:

Ki ∂ t U dx + Ni j=1 Eij T -1 ij G(T ij U ) ds - Ki (S g (U ) + S f (U )) dx = 0 ( 10
)
where N i denotes the number of faces of the cell K i (3 or 4), E ij is the cell interface and T ij is the rotation of angle θ ij between the normal to E ij and the horizontal plane, Fig. 4.5:

T =   1 0 0 0 cos θ sin θ 0 -sin θ cos θ  
Using the forward Euler scheme for temporal discretization, the equation becomes:

U m+1 i = U m i - ∆t |K i | Ni j=1 T -1 ij G m (T ij U i ) + ∆t S m g (U i ) + ∆t S m+1 f (U i ) (11)
where m = 0, ..., T /∆t is the time index, ∆t is the time step used for the 2D model integration. S m+1 f denotes the friction term treated semi-implicitely in time.

The difficulty is to compute the numerical flux

G m (T ij U i ) = Eij G(T ij U ) ds.
The rotational invariance property of the equations allows to reduce the 2D shallowwater problem to a sum of 1D local Riemann problems. These 1D Riemann problems are solved using the HLLC solver which consists to approximate the 1D flux G, see e.g. [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF] and [START_REF] Monnier | Well-balanced hllc solvers for nonhomogeneous shallow-water equations with consistent intermediate wave speed[END_REF]. Moreover, the topography term S g is included in the HLLC solver (see section 4.1.2). Thus, [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF] writes :

U m+1 i = U m i - ∆t |K i | Ni j=1 T -1 ij Gm (T ij U i ) + ∆t S m+1 f (U i ) (12)
Only the friction term is treated outside the HLLC numerical flux. It is discretized semi-implicitely in time as follows:

S m+1 (U i ) =     0 -g n 2 u h 4/3 m (q x ) m+1 -g n 2 u h 4/3 m (q y ) m+1    
where n is the manning coefficient.

The local 1D Riemann problem

We detail below the computation of an approximation of the 1D flux G. This is done as follows. 2)

If V = [H, Q n , Q τ ]
T where (τ, n) denotes the tangential and normal unit vectors respectively, we have the local normal 2D flux:

G(V ) = [Q n , Q 2 n /H + gH 2 /2, Q n Q τ H] T
In the homogeneous case (S g = 0), the 2D finite volume solver consists to solve, for each edge, the following local 1D Riemann problem:

∂ t V + ∂ n G(V ) = 0 with V (x, 0) = V L if x n < 0 ; V (x, 0) = V R if x n > 0.
The normal flux G(V ) is computed using the HLLC solver, [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF]. Finally, we obtain the new solution U m+1 i using [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF].

In the non-homogeneous case (S g = 0), the topography term is included in the flux term, using the following modification of G 2 (V ) proposed by LeVeque [START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods : The quasi-steady wave-propagation algorithm[END_REF] :

G(V ) = G(V ) +   0 1 2 g(h L + h R )(z bR -z bL ) 0   4.1.

Stability condition

Since we use the forward Euler scheme in time, we do not have linear system to solve but the time step must respect a stability condition (CFL type). In the present case, this CFL stability condition is, see [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF]:

CFL = ∆t max ( u + c) min(d L,R ) ≤ 1 , (13)
where d L,R is the distance between the cell L and the center of the edge separing it from cell R.

Boundary conditions

In this section, we describe the way we impose the different boundary conditions described in section 2.2.

Let us consider a "boundary cell" L, ie it has no neighboor through its edge i, (see figure 4.6). In the context of a finite volume scheme, we have to evaluate the flux through the edge i. Thus, it becomes natural to define a "ghost-cell" R, for each boundary cell L. In order to impose a boundary condition (e.g. inflow dicharge, incoming characteristics etc), we set a particular value to U R = (h R , u R , v R ) T such that the boundary condition is imposed when solving the local Riemann problem through the edge i (ie by computing the numeric flux G(U L , U R )). Let us detail the values U R required for each type of boundary conditions.

L R

Wall conditions On Γ wall , we consider slip type conditions, then:

Homogeneous Neumann condition on h is imposed by setting

h R = h L Slip condition on u is imposed by setting u R = -u L and v R = v L
Inflow conditions On Γ in , we impose one of the two following boundary conditions.

. Inflow discharge

The inflow discharge is imposed (unit in m 3 s -1 ). This corresponds to the discharge through the whole boundary and it is distributed in each cell such that the normal velocity u L is constant along the boundary.

-Homogeneous Neumann on h are prescribed by setting h R = h L .

-The inflow discharge is prescribed by setting u R = -q in (t) A n , where A n is the boundary wet section area, and v R = v L .

.

Incoming characteristics

Once incoming and outgoing characteristics w 1 , w 2 and w 3 are computed in the ghost cell as described in appendix A, we set:

-h R = w1-w3 2c -u R = w1+w3 2hR -v R = w 2

INRIA

Outflow conditions On Γ out , we impose one of the four following boundary conditions.

. Water elevation

-We impose h R = z prescribed -(z b ) R , where (z b ) R is the bathymetry in the ghost- cell. -v R = v L -We impose (u R + 2 √ gh R ) = (u L + 2 √ gh L ), that is u R = (u L + 2 √ gh L -2 √ gh R ).
.

Homogeneous Neumann We impose h R = h L , u R = u L and v R = v L .
. Incoming characteristics Once incoming and outgoing characteristics w 1 , w 2 and w 3 are computed in the ghost cell as described in appendix A, we set:

-h R = w1-w3 2c -u R = w1+w3 2hR -v R = w 2
. Rating curve -As described in section 2.2.3, we solve:

q R = f (h R ) u R + 2 √ gh R = u L + 2 √ gh L (14)
-In addition we impose: v R = v L

Adjoint code

In practice, there exists three approaches to compute the adjoint state variable, then the gradient of the cost function.

1. The discretized continuous gradient can be obtained from the adjoint model ( 6)-( 7) discretized using an appropriate numerical scheme which is then implemented.

2. The discrete gradient can be obtained from the adjoint of the direct numerical scheme and its implementation.

3. The "computational gradient" is obtained from the differentiation of the forward code directly.

This last approach ensures a better consistency between the computed cost function (including all types of errors -errors of discretization, rounding errors, iterative algorithms etc-) and its gradient since it is the computed cost function which is differentiated. A large part of this extensive task can be automated using algorithmic differentiation, see [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF]. In the case of DassFlow, the direct code is written in Fortran 95 and it is derived using the automatic differentiation tool Tapenade [START_REF] Hascot | Tapenade 2.1 user's guide[END_REF].

How to use the adjoint code?

We describe how to define the direct code, then what is the response of the adjoint code automatically generated and finally how to use it.

Let K be the space of control variables and Y the space of the forward code response. In the case of DassFlow, we have :

k = (y 0 , q in , z out ,
dY = ∂M ∂k (k) • dk
The adjoint model is defined as the adjoint operator of the tangent model. This can be represented as follows: ∂M ∂k (k) * : Y ′ -→ K ′ . It takes dY * ∈ Y ′ an input variable and provides the adjoint variable dk * ∈ K ′ at output, see figure 4.9:

dk * = ∂M ∂k (k) * • dY *
Now, let us make the link between the adjoint code and the "computational gradient". • dY * , dk

INRIA dY dk ∂M ∂k (k)
K ′ ×K = dY * , ∂M ∂k • dk Y ′ ×Y (15)
or, using the relations presented above:

dk * , dk K ′ ×K = dY * , dY Y ′ ×Y . (16)
If we set dY * = (0, 1) T and by denoting the perturbation vector dk = (δy 0 , δn, δz b , δq in ) T , we obtain:

0 1 , dy * dJ * Y ′ ×Y =       δy * 0 δn * δz * b δq in * δz * out       ,       δy 0 δn δz b δq in δz out       K ′ ×K
Moreover, we have by definition:

dJ = ∂J ∂y0 (k) • δy 0 + ∂J ∂n (k) • δn + ∂J ∂z b (k) • δz b + ∂J ∂qin (k) • δq in + ∂J ∂zout (k)
• δz out Therefore, the adjoint variable dk * (output of the adjoint code with dY * = (0, 1) T ) corresponds to the partial derivatives of the cost function J:

∂J ∂y0 (k) = y * 0 ∂J ∂n (k) = n * ∂J ∂z b (k) = z * b ∂J ∂qin (k) = q * in ∂J ∂zout (k) = z * out RR n 6150
In summary, in order to compute the "computational gradient" (partial derivatives of the cost function J using differentiation of the forward code), first, one runs the direct code then one runs the adjoint code with dY * = (0, 1) T as input.

Automatic differentiation Tapenade [START_REF] Hascot | Tapenade 2.1 user's guide[END_REF] is an automatic differentiation tool for Fortran programs. It is devopped by the Tropics team [28] at INRIA Sophia-Antipolis. Tapenade works using source code transformation : it builds the tangent and/or adjoint code automatically from the direct code written in Fortran. One refers also to [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF] to get a detailed description of how automatic differentiation works.

Adjoint code validation

We describe below how we check the validity of the adjoint code. Classicaly, we check that it is actually the adjoint of the tangent linear code (scalar product test) and that it computes correctly the partial derivative of the cost function (gradient test).

Scalar product test

The objective of this test is to check if the adjoint code is actually the adjoint of the tangent linear code. In other words, we check the relation ( 15) :

• Given an arbitrary dk ∈ K, we compute using the tangent linear code : dY = ∂M ∂k • dk • Given an arbitrary dY * ∈ Y, we compute using the adjoint code :

dk * = ∂M ∂k * • dY *
• Then, we compute the following scalar products :

. sp 1 = dY * , dY Y

. sp 2 = dk * , dk K

• And we check if sp 1 = sp 2 or not.

Gradient test

The objective of this test is to check if the adjoint variables dX * computed by the adjoint code correspond to the partial derivatives of the cost function. The Taylor expansion of the cost function j at k for a small perturbation α δk (where α ∈ R + ) writes :

j(k + α δk) = j(k) + α ∂j ∂k (k) • δk + o α δk . ( 17
)

INRIA

We set:

I α = j(k + α δk) -j(k) α ∂j ∂k (k) • δk . ( 18
)
According to [START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF], one must have: lim

α→0 I α = 1.
The gradient test consists to check this property:

• For an arbitrairy k, we compute ∂j ∂k (k) with the adjoint code. • With the direct code, we compute j(k).

• For n = 0, . . . , N :

We compute α = 2 -n ; With the direct code, we compute j(k + α δk) ; We compute I α ;

• We check if lim α→0 I α = 1 or not. . When α is smaller, the approximation of the partial derivatives is not reliable anymore due to truncation errors, see e.g. [START_REF] Elizondo | Automatic-versus manual-differentiation for non-linear inverse modeling[END_REF].

5 Benchmarks of the forward code

Explicit solutions

Academical tests cases with the exact solution known are performed in order to check the finite volume solver (forward code). The first test is the water at rest with a variable topography. We check if the finite volume solver is "well balanced" or not, see e.g. [START_REF] Monnier | Well-balanced hllc solvers for nonhomogeneous shallow-water equations with consistent intermediate wave speed[END_REF]. That means that it must preserves at least water at rest solutions.

|I α -1| α #################################### ## TEST DU PRODUIT SCALAIRE ## #################################### Appel du

Steady solution and Manning-Strickler formula

This test case aims to check the well known Manning-Strickler formula that links water height, discharge and manning coefficient in the case of a steady and uniform state. In the present case (a rectangular and very large channel), the water height h can be computed by: h = q 3/5 n 3/5 L -3/5 i -3/10 (19

)
where q is the water discharge, n the Manning coefficient, L the canal length and i = sin(α).

This test case consists to compare the steady state solution obtained by DassFlow, to the expression [START_REF] Hostache | Assimilation of spatial distributed water levels into a shallow-water flood model. part ii: Flood event at moselle river[END_REF].

This test case allows to validate or not the balance between the topography term and the friction term.

Dambreak

This test case simulates a dam-break with a flat bottom and without friction. The domain dimension is 50m × 10m, the topography is constant. There is a discontinuity at x = 25m. The initial condition is the following : The simulation time is 3s and the time step is 0.01. We perform this test case with a regular structured mesh (with n x = 200 and n y = 20) and with an unstructured triangular mesh (based on the same mesh for boundaries). Results are shown in figure 5.13.

h(x, y) = 1 if x ≤ 25 0.1 if x > 25 see fig.

Real data test case: Moselle river

In this section, we present a forward run on real data and compare the results with available measures. The study concerns the flood event from Feb.25 1997, 12h00 to Mars.2 1997, 12h00 of Moselle river. The mesh has 2340 cells and 2430 nodes, and it is a mix of triangles and quadrangles, see fig. 5.15. It includes the topography defined at the center of cells and was generated using the ArcGIS software. Time step is dt = 2s. This flood event test case is studied in [START_REF] Hostache | Assimilation of spatial distributed water levels into a shallow-water flood model. part ii: Flood event at moselle river[END_REF] and data come from [START_REF] Puech | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF].

Boundary conditions are imposed using real measures obtained during the flood event. The hydrographs showing the inflow discharge prescribed at inflow boundary and the water elevation prescribed at outflow boundary are shown in figure 5. [START_REF] Hascot | Tapenade 2.1 user's guide[END_REF].

In order to validate the numerical solution, we compare it to measures available at the EDF gauge station located approximatively in the middle of the main channel, see fig. 5.16. concerning data assimilation experiments conducted on this configuration we refer to [START_REF] Hostache | Assimilation of spatial distributed water levels into a shallow-water flood model. part ii: Flood event at moselle river[END_REF]. 6 Benchmarks for the full code (data assimilation)

In this section, we present some twin experiments results. A "twin experiment" means that a first run of the direct model provides observations. Then control variables are changed (we set them as a "first guess") and we run DassFlow in the data assimilation mode in order to try to retreive the control variable values that generated the observations.

The inflow boundary is at x = 0, the outflow boundary at x = 30. Boundaries y = 0 and y = 4 are walls. We use a rectangular structured mesh of dimension 90 × 20. Bed roughness, defined by its Manning coefficient, is uniform (n = 0.025). We impose a constant discharge q in = 8 m 3 /s at x = 0 and a constant water height h out = 1.4 m at x = 30. We obtain a steady state solution after about 80 s of simulation. Figure 6.17 shows the water height of this steady state solution and the topography. From this steady state solution, we extract the forthcoming observations: h obs and u obs every 0.02 s during 20 s on each volume. The objective of this test case is to retreive the topography. The first guess used is a flat bottom. We use the datassim mode, see [START_REF] Honnorat | DASSFLOW : USER GUIDE[END_REF], with the following cost function :

j 1 (z b ) = 1 2 T 0 h(t) -h obs (t) 2 Ω + q(t) -q obs (t) 2 Ω dt , (21)
Figure 6.18 shows the cost function and the norm of its gradient normalized by its initial values, vs iterates (a) and the identified topography (b). We can notice that convergence is obtained and the reference topography is well retreived. Further similar assimilation experiments but using lagrangian floaters are conducted in [START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF].

Identification of inflow discharge

We consider a toy test case which includes many features of real cases (eg. Moselle river). The computational domain contains a main channel (river) and floodplains, see figures 6.19 and 6.20).

Again, the present test case is a twin experiment. At the inflow boundary, we set the inflow discharge shown if figure 6.21 (a) simulating a flood event.

Then we perform a forward run to generate observations at points 1 and 2 shown with black stars in figure 6.20(a). Then, we suppose that the inflow discharge is constant (4.95 m 3 s -1 ), and we try to retreive its real value by assimilating observations. We present in Figure 6.21 the identified inflow discharge for different experiments. In Fig. 6.21(a), observations are h and q at each cell and each time step. In Fig. 6.21(b), observations are h at point 1 and (h, q) at point 2, both at each time step. In Fig. 6.21(c), observations are h at point 1 only, but at each time step.

We can notice that the identified inflow discharge is good even with the observation of h at point 1 only. Also, we can notice that the end of the flood event is not well identified. This is the "blind period" phenomena: for example in case (c), the inflow discharge after 270 s can not be identified because the information from the inflow boundary did not reach yet the gauge station.

Identification of Manning coefficients

We consider the same toy test case presented above but we seek to identify manning coefficients in given land-use. Manning coefficients are defined using five different land uses as follows, see figure 6.22:

In the main channel : n = 0.025 (gravelly main channel)

Right to the main channel : n = 0.066 (flood plain with bushes)

Left to the main channel : n = 0.04 (flood plain with little vegetation)

Near gauging station number 1 : n = 0.03 (pasture, farmland)

Near outflow boundary condition : n = 0.10 (urban)

Again, this test case is a twin experiment: a forward run generates the observations. For the present test case, the observations are water height at any time at gauge station 1, see figure 6.20(a). As first guess we consider a constant Manning (n = 0.010). The data assimilation process allows to recover perfectly the values used to generate the observations: Extra numerical experiments with real data (Moselle river, see [START_REF] Hostache | Assimilation of spatial distributed water levels into a shallow-water flood model. part ii: Flood event at moselle river[END_REF]), show the efficiency of the present approach.

Area

INRIA

A Characteristics boundary conditions

In this section, we explain how we define the open boundary conditions based on the theory of characteristics. We refer to [START_REF] Debreu | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF] and [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF]. Let us consider (1) near the boundary, in the non-conservative form and linearized around a mean value (h 0 , u 0 , v 0 ), with a flat bottom and without friction :

       ∂u ∂t + u 0 ∂u ∂x + v 0 ∂u ∂y g ∂h ∂x = 0 ∂v ∂t + u 0 ∂v ∂x + v 0 ∂v ∂y + g ∂h ∂y = 0 ∂h ∂t + u 0 ∂h ∂x + v 0 ∂h ∂y + h 0 ∂u ∂x + ∂v ∂y = 0 (22)
In a matrix form, this gives:

U t + A 1 U x + A 2 U y = 0 (23) with U = [h, u, v] T , A 1 =   u 0 0 g 0 u 0 0 h 0 0 u 0   , and A 2 =   v 0 0 0 0 v 0 g 0 h 0 v 0   Let n = [n 1 , n 2 ]
T and τ be respectively the normal and the tangent vector to the boundary (see Fig If we neglect the variations along x τ , (24) becomes a system of transport equations of w k at speed λ k in the normal direction n. Given an open boundary, w k is incoming if λ k < 0 and outgoing otherwise.

A.23). The matrix A = n 1 A 1 + n 2 A 2 has 3 eigenvectors: w 1 = u • n + g h0 h, w 2 = u • τ and w 3 = u • n -
Then, the characteristics boundary conditions are implemented as follows:

For each point of the boundary, we compute the eigenvalues λ 1 , λ 2 and λ 3 .

For each λ k , -if λ k < 0, then w k is incoming and must be specified from external data.

-if λ k ≥ 0, then w k is outgoing and must be computed from internal data.

Then the state variables can be recovered from characteristics values by:

-u • n bnd = w1+w3 2 -u • τ bnd = w 2 -h bnd = h0 g w1-w3 2
The remaining question is to compute the outgoing characteristics using internal data. Since the finite volume solver is explicit, at current time step n, we need values of the state variable in ghosts cells at time step n -1, see section 4.1.4. This is done using a simple linear extrapolation in the normal direction of the normal, see figure A.23.

B Riemann invariants and boundary conditions

This section details why the relation ∂ ∂n (u • n + 2c) = 0 is imposed if either water elevation or a rating curve is prescribed. We refer to ([2], section 1.3).

Let us consider the 1D shallow water equations with flat topography and without source term: The idea is to transform this equation in the form :

U t + B U x = 0 (
d dt (U + ω) = 0 ( 26
)
where ω is a quantity to be determined. Since gA b is not constant, we have to solve :

dω = gb A dh (27)
In case of a rectangular channel, we obtain: w = 2 √ gh. This means that (U + 2 √ gh) and (U -2 √ gh) are conserved along characteristic curves of direction λ 1 and λ 2 respectively.

For an outflow boundary (u • n > 0), λ 1 > 0 and (U + 2 √ gh) must be computed from the interior of the domain. Since the present finite volume solver is explicit in time, we need values of the state variable on the boundary at time n -1 when computing those of time n. Then, we use a simple extrapolation in the direction of the normal to compute the quantity (U + 2 √ gh) on the boundary. Using a 0-order extrapolation, we find the well-known condition :

∂(U + 2 √ gh) ∂n = 0
Link with the characteristics boundary conditions RR n 6150

Figure 2 . 1 :

 21 Figure 2.1: The simulation domain and its boundaries

Figure 2 . 2 : 0 3

 220 Figure 2.2: Outflow BC: Rating curve and U + 2 √ gh = 0

Figure 3 . 3 :

 33 Figure 3.3: Flux through all edges of cell L

Figure 3 . 4 :

 34 Figure 3.4: Variational data assimilation process

Figure 4 . 5 :

 45 Figure 4.5: Adjacent cells and use of the rotational invariance property of the equations

Figure 4 . 6 :

 46 Figure 4.6: Internal cell L at boundary and ghost cell R

Figure 4 . 7 :

 47 Figure 4.7: Representation of the direct model. The tangent model becomes ∂M ∂k (k) : K -→ Y. It takes as input variable a perturbation of the control vector dk ∈ K, then it gives the variation dY ∈ Y as output variable, see figure 4.8) :

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Representation of the tangent model.

Figure 4 .

 4 Figure 4.10 (b) shows an example of the scalar product test.

Figure 4 .

 4 Figure 4.10 (a) shows a result of the gradient test. |I α -1| is plotted against α in logarithmic scale. The convergence is good until α > 10 -7. When α is smaller, the approximation of the partial derivatives is not reliable anymore due to truncation errors, see e.g.[START_REF] Elizondo | Automatic-versus manual-differentiation for non-linear inverse modeling[END_REF].

Figure 4 . 10 :

 410 Figure 4.10: Adjoint code validation

Figure 5 . 11 :

 511 Figure 5.11: Steady state solution for the Manning test case

Figure 5 . 12 :

 512 Figure 5.12: Initial condition on h for the dambreak test case

Figure 5 . 13 :

 513 Figure 5.13: Dam break test case (all y values are superposed). a) and b) on a regular structured mesh, c) and d) on an unstructured triangular mesh. a) and c) represent h, b) and d) represent u.

Figure 5 . 14 :

 514 Figure 5.14: Moselle river. a) Inflow discharge imposed. b) Elevation prescribed at outflow boundary.

Figure 5 . 15 :

 515 Figure 5.15: Moselle river. Mesh and water depth computed on Feb.26 1997 14h00

Figure 5 . 16 :

 516 Figure 5.16: Moselle river. Comparison of the measured and computed hydrographs at EDF gauge station

Figure 6 . 17 :

 617 Figure 6.17: Identification of the topography. Topography and steady state elevation.

Figure 6 . 18 :

 618 Figure 6.18: Value of the cost function and of the norm of its gradient, normalized by their initial values (a) and the identified topography (b)

Figure 6 . 19 :Figure 6 . 20 :

 619620 Figure 6.19: Toy test case mesh (a) and bathymetry (b)

Figure 6 . 21 :

 621 Figure 6.21: In green : reference inflow discharge. In red : identified inflow discharge. (a): Observation of (h, q) everywhere; (b): Observation of h at point 1 and (h, q) at point 2. (c): Observation of h at point 1 only.

Figure 6 .

 6 Figure 6.22: Manning coefficients

g+ λ 3 ∂w3∂xn + u 0 • τ ∂w3 ∂xτ -c ∂v ∂xτ = 0 ∂w2∂t + λ 2 ∂w2∂xn + u 0 • τ ∂w2 ∂xτ + c 2 ∂

 3022 h0 h. They are associated to the eigenvaluesλ 1 = u 0 • n + c, λ 2 = u 0 • τ and λ 3 = u 0 • n -c (c = √ gh 0 )respectively. These eigenvectors w 1 , w 2 and w 3 are the so-called characteristic variables.We rewrite[START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods : The quasi-steady wave-propagation algorithm[END_REF], using w 1 , w 2 and w 3 :

Figure A. 23 :

 23 Figure A.23: Computation of characteristic variables in a ghost cell using a linear extrapolation of the values in cells (black points).

+ λ 1 ∂

 1 25) INRIA with U = [h, U ] T , B = U A/b g U , and b(x) is the canal width. The matrix B has 2 eigenvalues λ 1 = U + g A b and λ 2 = U -g A b . The corresponding left eigenvectors are Λ 1 = g, ∂x represents the Lagrangian derivative in the direction of the curve of slope λ 1 .

  code linaire tangent...

		Appel du code adjoint...
		<Xd,Xb> =	-682.277083033688
		<Yd,Yb> =	-682.277082428555
		relative error : -8.869324203484993E-010
	(a) Gradient test	(b) Scalar product test

RR n 6150

Honnorat & Marin & Monnier & Lai

One can notice that the transported quantities obtained are different from those of appendix A. This can be explained by the use of linearized equations in appendix A.

Let us linearize [START_REF] Puech | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]

We get:

By setting in [START_REF] Puech | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF], we get:

Assuming h ′ and U ′ are small, we find the expression of section A :


Frédéric Giroire
email: frederic.giroire@cnrs.fr

Frédéric Havet

Joanna Moulierac

Compressing Two-dimensional Routing Tables with Order ⋆

Keywords: routing, compact routing tables, software defined networks, complexity, approximation algorithms

A communication in a network is a pair of nodes (s, t). The node s is called the source source and t the destination. A routing of a communication (s, t) is a path in the network from s to t. A routing of a communication set is a union of routings of its communications.

At each node, there is a set X of communications whose routing path goes through this node. The node needs to find for each communication (s, t) in X , the port that the routing path of (s, t) uses to leave it. An easy way of doing it is to store the list of all triples (s, t, k), where (s, t) ∈ X and k is the port used by the (s, t)-path to leave the node.

However, such a list might be very large. Motivated by routing in telecommunication network using Software Defined Network (SDN) technologies, we consider the problem of compacting this list using aggregation rules. Hence, in addition, we can use some additional triples, called * -triples. As an example, a t-destination triple ( * , t, p), means that every communication with destination t leaves on port p.

We carry out in this work a study of the problem complexity, providing results of NP-completeness, of Fixed-Parameter Tractability and approximation algorithms.

Introduction

Motivation. Software Defined Technology (SDN), e.g. OpenFlow [START_REF] Mckeown | Openflow: Enabling innovation in campus networks[END_REF] is a new promising approach to operate telecommunication networks. Its promise is to allow dynamic routing decisions by decoupling the control plane (the system making decisions) from the data plane (which forwards the packets). This way, a centralized controller receives the data monitored in the system (e.g. load, delay, ...) and then, based on this information, computes appropriate routing decisions, e.g. to improve energy efficiency [START_REF] Giroire | Minimizing routing energy consumption: from theoretical to practical results[END_REF][START_REF] Giroire | Minimization of network power consumption with redundancy elimination[END_REF]. Each time a new flow arrives, the router contacts the controller and waits for the decisions to be pushed into its forwarding table. The routing tables thus are populated with flow-based rules with header informations (source IP, destination IP, ...) → exit port.

However, SDN hardware uses specific memory, e.g. TCAM memory [START_REF] Kang | Optimizing the "one big switch abstraction" in software-defined networks[END_REF][START_REF] Kanizo | Palette: Distributing tables in softwaredefined networks[END_REF], which is very expensive and of small size. Thus, the number of entries of the routing tables is limited to only a few thousands [START_REF] Narayanan | Macroflows and microflows: Enabling rapid network innovation through a split sdn data plane[END_REF][START_REF] Stephens | Past: Scalable ethernet for data centers[END_REF] and grows linearly with the number of flows passing through a router, causing a problem of scalability. It is thus an important area of research to obtain routing using only a limited number of rules per router. [START_REF] Cohen | On the effect of forwarding table size on SDN network utilization[END_REF] studies the problem of choosing routing with a limited number of entry per router using linear programming. Another way to compact the forwarding tables is to use aggregation rules. With such an aggregation, we can set routing entries such as "(*,destination)→ port" or "(source,*)→ port" or also a default entry such as "(*,*)→ port". For example, [START_REF] Giroire | Optimizing Rule Placement in Software-Defined Networks for Energy-aware Routing[END_REF] studies how to use default ports to reduce the size of routing tables. In this work, we consider the problem of compressing a routing table using aggregation rules.

We consider here two-dimensional routing tables in which the routing decision is not done exclusively on the destination IP address, but on the source and destination IP addresses. Indeed, the commonly implemented destinationbased routing has its limitations, especially in delivering quality of service which is a goal of SDN paradigm. One suggested remedy is to base the routing decision on additional fields in the packet header. One of the most important field is the source host. For instance, this would permit selective routing to provide a high bandwidth connection between two different sites of a company. Such refined forwarding is part of the next generation Internet design, and falls within the broader scope of layer four packet classification, where packets are routed using arbitrary fields of the packet header [START_REF] Hari | Detecting and resolving packet filter conflicts[END_REF], [START_REF] Eppstein | Internet packet filter management and rectangle geometry[END_REF], [START_REF] Lakshman | High-speed policy-based packet forwarding using efficient multi-dimensional range matching[END_REF], [START_REF] Buddhikot | Space decomposition techniques for fast layer-4 switching[END_REF]. Routers capable of packet classification can implement many advanced services, such as firewall access control, Virtual Private Networks, and quality of service routing, which are all promises of the SDN paradigm.

Modeling. A communication in a network is a pair of nodes (s, t). The node s is called the source source and t the destination. We use the source and destination fields in our examples, although our ideas apply to any two prefix fields in Internet protocol networks. A communication set is a set of distinct communications, i.e. two communications might have the same source or the same destination, but they cannot have both same source and same destination. A routing of a communication (s, t) is a path in the network from s to t. A routing of a communication set is a union of routings of its communications.

At each node, there is a set X of communications whose routing path goes through this node. The node needs to find for each communication (s, t) in X , the port that the routing path of (s, t) uses to leave it. An easy way of doing it is to store the list of all triples (s, t, k), where (s, t) ∈ X and k is the port used by the (s, t)-path to leave the node. Such triples are called communication triples.

However, such a list might be very large. So we want to reduce it as much as possible using the * symbol. Hence, in addition, we can use some additional triples, called * -triples. There are two kinds of * -triples:

• t-destination triple ( * , t, p), meaning that every communication with destination t leaves on port p.

• s-source triple (s, * , p), meaning that every communication with source t leaves on port p.

A routing list is an ordered list T 1 , . . . , T r of triples (either communication, or source, or destination ones). A communication is then assigned the port of the first triple in the list, that applies to it. It is crucial to remark that using * -triples introduces an order of the rules in the routing list.

Let C be a set of communication triples. A routing list R emulates C if each communication of C is assigned the same port by C and R. Observe that R may route more communications than C. For example, if the port of all triples of C have source s and port p, then the singleton list made of the source triple (s, * , p) emulates C, even if there is not a triple in C for all communications.

Problems. The problem is then to find the shortest routing list that emulates C. We denote by rmin(C) the minimum number of triples in a routing list emulating C.

Routing List:

Input: A set C of communication triples and an integer r. Question: rmin(C) ≤ r?

The number of saved triples is sav(C) = |C| -rmin(C). The complementary problem to Routing List is the following.

List Reduction

Input: A set C of communication triples and an integer z. Question: sav(C) ≥ z?

Contributions. In this work, we study the complexity of the above problems. We provide NP-completeness results and an approximation algorithm. Due to lack of space, some proofs are omitted. Full proofs can be found in [START_REF] Giroire | Compressing Two-dimensional Routing Tables with Order[END_REF].

Our work answers an open question of [START_REF] Suri | Compressing two-dimensional routing tables[END_REF]. Similarly to us, the authors consider the problem of determining a compact routing table using aggregation rules that has the same behavior as the original routing table. The difference with our problem is that their goal is to find what they call a conflict-free routing table in which the rules can be taken in any order. On the contrary, as noted above, the order is crucial in our problems.

Results

We show that Routing List and List Reduction are NP-complete by reduction from Feedback Arc Set problem, one of Karp's 21 NP-complete problems [START_REF] Karp | Reducibility among combinatorial problems[END_REF].

Theorem 2.1 Routing List and List Reduction are NP-complete, even if there are only two ports.

Because of the NP-hardness of the problem it is interesting to get approximation algorithms. We show that sav can be 2-approximated in polynomial time using the simple heuristic described below:

The Direction-based Heuristic Let C be a set of communication triples with destination set S and destination set T . Set n = |S| and m = |T |. For any port p, let C(p) be the set of triples of C with port p. For a source s (resp. destination t) and a port p, let C(s, p) (resp. C(t, p)) be the set of triples of C with source s (resp. destination t) and port p. For every source s, let M(s) := max p |C(s, p)| be the maximum number of triples in C with source s and same port, and for every destination t, let M(t) := max p |C(t, p)| be the maximum number of triples in C with destination t and same port. Set

Z -(C) = s∈S (M(s) -1) = s∈S M(s) -n and Z | (C) = t∈T (M(t) -1) = t∈T M(t) -m.
Any routing list emulating C yields an upper bound on rmin(C). One such list can be obtained by routing source by source. One source s after another we route all triples of C with source s. This can be done by using the triple (s, * , p) for p a port such that there are M(s) triples with source s and port p after all triples with source s and port distinct from p. Doing so, we save M(s) -1 triples when routing the triples with source s. Hence, we obtain a routing list of size |C| -Z -(C). Such a list is called a source-based routing list.

Proceeding similarly according to the destinations, we obtain a routing list, called destination-based of size |C| -Z | (C).

Setting Proof. Let C be a set of communication triples with destination set S and destination set T . Set n = |S| and m = |T |. We can order S and T by decreasing order according to the function

Z(C) = max{Z -(C), Z | (C)},
M. That is M(s 1 ) ≥ M(s 2 ) ≥ • • • ≥ M(s n ) and M(t 1 ) ≥ M(t 2 ) ≥ • • • ≥ M(t m ).
The directed {0, 1, . . . , n} × {0, 1, . . . , m}-grid, denoted by G n,m is the digraph defined by

V (G n,m ) = {(i, j) | 0 ≤ i ≤ n and 0 ≤ j ≤ m} A(G n,m ) = A h (G n,m ) ∪ A v (G n,m ), where A h (G n,m ) = {a h (i, j) | 1 ≤ i ≤ n and 0 ≤ j ≤ m}, with a h (i, j) = ((i -1, j), (i, j))
, is the set of horizontal arcs, and A v (G n,m ) = {a v (i, j) | 0 ≤ i ≤ n and 1 ≤ j ≤ m}, with a v (i, j) = ((i, j -1), (i, j)), is the set of vertical arcs.

An edge-weighted digraph is a pair (G, w) where G is a digraph and w a weight function on its arcs, that is a function from A(D) into R. The length of a path U in an edge-weighted digraph (G, w) is the sum of the weights of its arcs: w(U) = a∈A(U ) w(a).

Let w C be the weight function defined on A(G n,m ) by w C (a h (i, j)) = min{M(s i ), m-j}-1 and w C (a v (i, j)) = min{M(t j ), n-i}-1. Let W (C) be the maximum length of a path in (G n,m , w C ). Observe that it is attained by a path from (0, 0) to one of the sides {n} × {0, 1, . . . , n} and {0, 1, . . . , n} × {n} because the weight of an arc is negative if and only if it is in

{a h (i, m) | 1 ≤ i ≤ n} ∪ {a v (n, j) | 1 ≤ j ≤ m}.
We shall prove that

sav(C) ≤ W(C). (2)
Let R be a routing list that emulates C. One can show (see [START_REF] Giroire | Compressing Two-dimensional Routing Tables with Order[END_REF]) that we can freely rearrange the triples of R, so that R is canonical. A routing list is canonical if it is the concatenation of sublist B 1 , . . . , B q , called blocks having the following properties for every 1 ≤ ℓ ≤ q:

(i) in B ℓ , there is a unique * -triple and it is the last one;

(ii) if the * -triple of B ℓ is an s-source triple (resp. t-destination triple), then all triples of B ℓ have source s (resp. destination t).

For 1 ≤ k ≤ q, let T ℓ k be the * -triple of B k , let r k be the number of triples of C routed by T ℓ k , and let i k (resp. j k ) be the number of source triples (resp. destination triples) with index at most ℓ k . We have |R| = |C|-Σ with Σ = q k=1 (r k -1).

Let K s (resp. K t ) be the set of indices k such that T ℓ k is a source (resp. destination) triple. For k ∈ K s , s ′ k be the source of T ℓ k , and for k ∈ K t , t ′ k be the destination of T ℓ k .

Assume k ∈ K s . When routing according to R, before considering the block B k , there are at most mj k triples of C with source s ′ k and at most

M(s ′ k ) triples of C with source s ′ k . Therefore, r k ≤ min{M(s ′ k ), m -j k }. Similarly, if k ∈ K t , then r k ≤ min{M(t ′ k ), n -i k }. Thus Σ ≤ k∈Ks (min{M(s ′ k ), m -j k } -1) + k∈Kt (min{M(t ′ k ), n -i k } -1).
Set u 0 = (0, 0) and u k = (i k , j k ) for 1 ≤ k ≤ q. Note that i k + j k = k for all 0 ≤ k ≤ q. The sequence U = u 0 , u 1 , . . . , u q can be seen as a path in the directed grid G n,m . Observe that if k ∈ K s then u k-1 u k is a horizontal arc and w C (u k-1 u k ) = min{M(s i k ), mj k } -1, and if k ∈ K t then u k-1 u k is a vertical arc and w C (u k-1 u k ) = min{M(t j k ), ni k } -1. Hence w C (U) = k∈Ks (min{M(s i k ), mj k } -1) + k∈Kt (min{M(t j k ), ni k } -1). Now since M(s 1 ) ≥ M(s 2 ) ≥ • • • ≥ M(s n ) and j 1 ≤ j 2 ≤ • • • ≤ j q , we have k∈Ks (min{M(s ′ k ), mj k } -1) ≤ k∈Ks (min{M(s i k ), mj k } -1). Similarly, k∈Kt (min{M(t ′ k ), ni k } -1) ≤ k∈Kt (min{M(t j k ), ni k } -1). Hence Σ ≤ w C (U). Thus Σ ≤ W (C), which implies [START_REF] Cohen | On the effect of forwarding table size on SDN network utilization[END_REF].

One can check (see [START_REF] Giroire | Compressing Two-dimensional Routing Tables with Order[END_REF]) that W (C) ≤ 2Z(C), yielding the theorem. ✷

Conclusion

We studied the complexity of Routing List and of List Reduction. We provide results of NP-completeness and an approximation algorithm.

We leave several questions as open problems:

• The Destination-based Heuristic is a 2-approximation for List Reduction. However, the heuristic often returns a routing list that saves more than a half of sav(C) triples. This leads to think that the approximation of 2 is not best possible.

  we have sav(C) ≥ Z(C) and rmin(C) ≤ |C| -Z(C). (1) The algorithm consisting in computing a source-based routing list and a destination-based routing list and taking the shortest of the two, is called the Direction-based Heuristic. It provides a routing list emulating C of size Z(C). The Direction-based Heuristic is a 2-approximation for List Reduction.

Theorem 2 . 2

 22 The Destination-based Heuristic is a 2-approximation for List Reduction.

Problem 3.1 What is the best approximation ratio for List Reduction? • We proved the NP-Completeness of Routing List by a reduction to Feedback Arc Set. Knowing that Feedback Arc Set is APX-complete, Problem 3.2 Is Routing List also APX-complete? For which approximation ratio?

⋆ This work has been partially supported by ANR project Stint under reference ANR-13-BS02-0007, ANR program Investments for the Future under reference ANR-11-LABX-0031-01, ANR VISE, CNRS-FUNCAP project GAIATO, the associated Inria team AlDyNet, the project ECOS-Sud Chile.