Nhung Pham

Philippe Helluy

Laurent Navoret

HYPERBOLIC APPROXIMATION OF THE FOURIER TRANSFORMED VLASOV EQUATION

published or not. The documents may come

Introduction

Solving the Vlasov-Poisson equation is challenging. Some popular methods for studying this equation are the Particle-In-Cell (PIC) method [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] or the semi-lagrangian approach [5].

In a previous work [START_REF] Helluy | Space-only hyperbolic approximation of the Vlasov equation[END_REF], we constructed a reduced Vlasov-Poisson model with a velocity basis expansion.

In this paper, we consider a Fourier velocity transformation of the Vlasov equation. We construct a reduced model where the unknown depends on space and time instead of the full phase-space variables. The reduced model is a linear hyperbolic system, with non-linear source terms. We present numerical results for classical plasma physics test cases.

Plasma mathematical model

In our work, we consider the one-dimensional Vlasov equation

∂ t f + v∂ x f + E∂ v f = 0, (1)
where the unknown distribution function f depends on the space variable x ∈ R/LZ, the velocity variable v ∈ R and the time variable t ∈ R + . The electric eld E depends on x and t and is the solution of the Poisson equation

∂ x E = -1 + ˆv f, ˆL x=0 E = 0. (2)
The equations (1-2) are supplemented by an initial condition

f (x, v, 0) = f 0 (x, v).
1 IRMA, 7 rue Descartes, 67000 Strasbourg, pham@math.unistra.fr 2 Inria TONUS, philippe.helluy@inria.fr 3 Inria TONUS, laurent.navoret@math.unistra.fr c EDP Sciences, SMAI 2013 ˆv f = 1.

For practical reasons, we will allow that f and E take complex values, however, of course, only the real parts are physically relevant.

We consider a Fourier transformation with respect to the velocity variable (we denote by

I = √ -1) φ(x, η, t) = ˆ+∞ v=-∞ f (x, v, t) exp(-Iηv)dv.
The Fourier velocity variable is denoted by η ∈ R.

∂ t φ + I∂ x ∂ η φ + IEηφ = 0. (3)
In addition, the Poisson equation becomes

∂ x E(x, t) = -1 + φ(x, 0, t). (4)
We call the new model the Vlasov-Fourier equation.

Discretization of the Vlasov-Fourier equation with respect to the Fourier velocity variable

We will perform a semi-discretization of (3) with respect to the variable η in order to obtain a rst order hyperbolic system set only in (x, t). We shall call this new system of equations the reduced Vlasov-Fourier model. We could expand the function φ on a basis of arbitrary functions depending on η. See for instance [START_REF] Bourdiec | Numerical solution of the Vlasov-Poisson system using generalized Hermite functions[END_REF] and included references. For practical numerical reasons, we decide to choose a classical Lagrange nite element interpolation basis, because it leads to a sparse matrix representation of the hyperbolic system.

Continuous interpolation by the nite element method

In practice, φ almost vanishes at the boundaries η → ±∞. We consider thus a truncated domain η ∈ [-η max , η max ] and the following dissipative boundary conditions at ±η max ∂ x φ(x, ±η max , t) ± Iγφ(x, ±η max , t) = 0, (5) in which γ ≥ 0.

We recall now how the nite element basis is constructed. We consider an arbitrary polynomial degree d. The reference element is dened by

Q = [-1, 1].
We dene the d + 1 reference nodes by

Ni = -1 + 2 i -1 d , i = 1 • • • d + 1.
We mesh the interval [-η max , η max ] with N nite elements

(Q i ) i=1•••N and nodes (N j ) j=1•••P . The total number of nodes in this interval is P = d • N + 1.
In practice, we suppose that the nodes are equally spaced in [-η max , η max ]

N j = -η max + 2η max dN (j -1).
We introduce a connectivity array for detecting that node N j is the k th local node of a given element

Q i j = connec(k, i) = k + (i -1)d, 1 ≤ k ≤ d + 1, 1 ≤ i ≤ N.
We also use the notation

N j = N k,i and then, element Q i has its support in the interval [N 1,i , N d+1,i ].
We construct a transformation τ i that maps element Q onto Q i . For this purpose we consider the Lagrange polynomials on Ê, dened by

Lk (η) = l =k η -Nl Nk -Nl . ( 6
)
The transformation is then given by

τ i (η) = d+1 k=1 Lk (η)N k,i . (7)
Because the nodes of the mesh are equally spaced in our application, the transformation τ i is linear. We construct the interpolation basis in such a way that each basis function ϕ j is associated to a node N j of the mesh and satises

ϕ j (N i ) = δ ij ,
where δ ij denotes the Kronecker symbol. We recall how to compute the basis function ϕ

j . Let η ∈ [-η max , η max ].
Necessarily, η belongs at least to one nite element Q i . Two cases are possible

(1) Node N j belongs to nite element Q i , i.e. ∃k,

N j = N k,i , then ϕ j (η) = Lk (η), where η = τ i (η). (8)
(2) Node N j does not belong to Q i , then ϕ j (η) = 0.

Application to Vlasov-Fourier discretization

We suppose that the function φ(x, η, t) is well approximated by an expansion on the basis

{ϕ j } j=1•••P φ(x, η, t) = P j=1 w j (x, t)ϕ j (η), (9)
we shall also use the convention of sum on repeated indices φ(x, η, t) = w j (x, t)ϕ j (η).

(10)

Because of the interpolation property of the basis{ϕ j } j=1

•••P ϕ i (N j ) = δ ij , we have φ(x, N i , t) = P j=1 w j (x, t)ϕ j (N i ) = w i (x, t).
Therefore, we can approximate the initial condition in the following way

w j (x, 0) = φ(x, N j , 0) = φ 0 (x, N j ).
Considering the equation (3) and the boundary condition (5) we can consider the following weak formula of the problem : nd φ(x, η, t) such that for all (continuous) test function ϕ(η) we have

ˆη ∂ t φϕ + ˆη I∂ x ∂ η φϕ + ˆη IEηφϕ - 1 2 ϕ(η max )I∂ x φ(•, η max , •) + 1 2 ϕ(-η max )I∂ x φ(•, -η max , •) + 1 2 ϕ(η max )γφ(•, η max , •) + 1 2 ϕ(-η max )γφ(•, -η max , •) = 0 (11)
This "semi-weak" formula is equivalent with the initial problem (3) supplemented with the boundary conditions (5).

Indeed, if φ is a solution of (3) with the conditions (5), it is evident that φ is also a solution of (11).

Reciprocally, if we suppose that φ is a solution of (11), because (11) is true for arbitrary test function. Thus for every function ϕ such that ϕ(-η max ) = ϕ(η max ) = 0, we obtain ˆη ∂ t φϕ + ˆη I∂ x ∂ η φϕ + ˆη IEηφϕ = 0 and thus (3).

Then taking test functions ϕ that do not vanish at η = ±η max we obtain the boundary conditions (5).

Morever we show now that the choice of boundary conditions ( 5) ensures the hyperbolicity of the system and ensures that the energy is not increasing.

We can introduce the following matrices of dimension

P × P M = ( ˆη ϕ i ϕ j ), A = I Ã, B = B E + D, (12)
where

Ãij = ˆη ϕ i ϕ j - 1 2 ϕ i (η max )ϕ j (η max ) + 1 2 ϕ i (-η max )ϕ j (-η max ), (B E ) ij = IE ˆη ηϕ i ϕ j , D ij = 1 2 γ(ϕ i ϕ j (η max ) + ϕ i ϕ j (-η max )).
We obtain the following equation

M ∂ t W + A∂ x W + BW = 0, (13)
in which W (x, t) is the complex vector of P components W = (w 1 , w 2 , ..., w P ) T .

Obviously, the mass matrix M is positive hermitian. An integration by parts in η shows that A is hermitian. Finally, B E is skew-hermitian and D is diagonal non-negative. It is then classical to prove that system (13) is hyperbolic (i.e. that M -1 A is diagonalizable with real eigenvalues [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]) and energy dissipative.

In practice, to compute the matrices M, A, B we use the Gauss-Legendre integration and sparse matrix representations.

Finite volume schemes

We describe now the numerical approximation. We assume that the spatial domain ]0, L[ is split into N x cells. The cell C i is the interval x i-1/2 , x i+1/2 , i = 1..N x . For practical reasons, we also consider two virtual cells C 0 and C Nx+1 for applying the periodic boundary condition. At the beginning of a time step, we copy the values of the cell C Nx to the cell C 0 , and the values of the cell C 1 to the cell C Nx+1 . The center of the cell C i is

x i = i x -x 2 .
The space step is x = L/N x . We also consider a sequence of times t n , n ∈ N, such that t 0 = 0 and t n = n t, where t satises the following CFL condition

t = α x η 2d , 0 < α ≤ 1.
We consider a nite volume approximation of (13). We denote by W i (t) a piecewise constant approximation of W in each cell W i (t) W (x, t), x ∈ C i .

We obtain the following semi-discrete (in space) approximation

M ∂ t W i = - F (W i , W i+1 ) -F (W i-1 , W i ) x -BW i .
where

(W L , W R ) → F (W L , W R ) denotes the numerical ux.
We then introduce a time discretization to compute

W n i W n i W (x, t n ), x ∈ C i .
We use a time second order scheme given by the following algorithm

M W n+1/2 i -W n i t/2 = - F (W n i , W n i+1 ) -F (W n i-1 , W n i ) x -BW n i , M W n+1 i -W n i t = - F (W n+1/2 i , W n+1/2 i+1 ) -F (W n+1/2 i-1 , W n+1/2 i ) x -BW n+1/2 i . ( 14
)
We consider several choices for the numerical ux F (W L , W R ). We consider the centered ux or a numerical ux with small numerical viscosity ("slightly upwinded ux"). The centered ux is given by

F (W L , W R ) = A W L + W R 2 ,
and the slightly upwinded ux with δ > 0

F (W L , W R ) = A W L + W R 2 - δ 2 (W R -W L ). (15)
In practice, for saving CPU time and memory, we use two subroutines for computations with the sparse matrices in skyline format. The rst one computes the product of a sparse matrix and a vector. The other one is used to solve a linear system (by the LU method).

Test cases

In this section, we will apply our numerical scheme, (14-15) and ( 12), on two test cases: the Landau damping and the two stream instability. The electric eld, solution of the Poisson equation ( 4), is computed with the FFT (Fast Fourier Transform) algorithm.

We will compute the electric energy dened by

E(t) = ˆL 0 E(x, t) 2 dx.
We are also interested in the distribution function in physical variable (x, v). The formula of the inverse Fourier transform reads

f (x, v, t) = 1 2π ˆ+∞ -∞ φ(x, η, t)e Iηv dη 1 2π ˆηmax -ηmax φ(x, η, t)e Iηv dη. (16)
We apply the rectangle method with oversampling for computing (16), in order to avoid Shannon aliasing. For this computation we can use a naive DFT computation instead of the FFT algorithm, because this step is applied only at the beginning and the end of the simulation.

In our numerical experiments, the discretization parameters are N = 40, d = 5, N x = 256, γ = 0. Numerical investigations for γ > 0 will be carried out in future work.

The Landau damping

In this test case, the initial distribution function and the initial electric eld are given by

f 0 (x, v) = (1 + ε cos(kx)) 1 √ 2π e -v 2 2 , E 0 (x) = ε k sin(kx),
where ε > 0 and k ∈ N * , and the domain size is L = 2π/k.

So, the initial distribution function in velocity Fourier variable writes For small ε, thanks to a linear approximation of the non-linear Vlasov-Poisson system, it is possible to compute an approximate analytical solution of the electric eld. The details of the computation are given in [START_REF]Sonnendrücker Approximation numerique des equations de Vlasov-Maxwell[END_REF].

In addition, the distribution function can be computed by a well-validated method, such as the PIC method. Let us take the values of parameter k = 0.2 and ε = 5 × 10 -2 . We compare our numerical results with the PIC results and also with the analytical solution.

We compare the distribution function and the electric energy of the Vlasov-Fourier method and of the PIC method (taken from [3]), for example at time t = 100 on Figure 1. We also compare the time evolution of the electric energy in the domain obtained by our method with the analytical solution (see Figure 2). The value of parameters for this test case are k = 0.2, ε = 5 × 10 -3 δ = 0.05 and v 0 = 3. The distribution function is plotted at times t = 25 and t = 50 in Figures 3 and4. We compare the PIC method and the Vlasov-Fourier method with the centered ux. At time t = 50, we remark small oscillations, which are maybe due to the fact that we have almost no upwind mechanism in the resolution of the transport equation. In order to remove these oscillations, we use the slightly upwinded ux (15) with parameter δ = 0.05 instead of the centered ux and we obtain the results of Figure 5 in which the oscillations have disappeared.

Figure 1 .

 1 Figure 1. The distribution function of the Landau damping test case at time t = 100. Left: Vlasov-Fourier method. Right: the PIC method.

Figure 2 .

 2 Figure 2. The electric energy of the Landau damping test case up to time t = 100, the green curve asymptotic is the analytical solution and the red curve energy is computed with the Vlasov-Fourier method.

  x, v)e -Iηv dv = (1 + ε cos(kx))e -η 2 2 .

Figure 3 .

 3 Figure 3. The distribution function of the two-stream test case at time t = 25. Left: Vlasov-Fourier method (with the centered ux). Right: the PIC method.

Figure 4 . 2 + 2 , 1 2 η 2 cos

 42212 Figure 4. The distribution function of the two-stream test case at time t = 50. Left: Vlasov-Fourier method (with the centered ux). Right: the PIC method.

  The distribution function φ(x, η, t) satises the Fourier trans-

	formed Vlasov equation [9]


Emmanuel Frenod
email: emmanuel.frenod@univ-ubs.fr

Ibrahima Faye
email: ibrahima.faye@uadb.edu.sn

Diaraf Seck
email: diaraf.seck@ucad.edu.sn

F S T Emmanuel Frénod

TWO-SCALE NUMERICAL SIMULATION OF SAND TRANSPORT PROBLEMS

Keywords: 2000 Mathematics Subject Classification. Primary: 35K65, 35B25, 35B10 ; Secondary: 92F05, 86A60 Homogenization, Asymptotic Analysis, Asymptotic Expansion, Long Time Behavior, Dune and Megaripple Morphodynamics, Modeling Coastal Zone Phenomena, Numerical Simulation

Two-Scale numerical simulation of sand transport problems

1.

Introduction. This paper deals with numerical simulations of sand transport problems. Its goal is to build a Two-Scale Numerical Method to simulate dynamics of dunes in tidal area. This paper enters a work program concerning the development of Two-Scale Numerical Methods to solve PDEs with oscillatory singular perturbations linked with physical phenomena. In Ailliot, Frénod and Monbet [START_REF] Aillot | Long term object drift in the ocean with tide and wind[END_REF], such a method is used to manage the tide oscillation for long term drift forecast of objects in coastal ocean waters. Frénod, Mouton and Sonnendrücker [START_REF] Frénod | Sonnendrücker Two-Scale numerical simulation of the weakly compressible 1D isentropic Euler equations[END_REF] made simulations of the 1D Euler equation using a Two-Scale Numerical Method. In Frénod, Salvarani and Sonnendrücker [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a Two-Scale PIC-method[END_REF], such a method is used to simulate a charged particle beam in a periodic focusing channel. Mouton [START_REF] Mouton | Approximation multi-échelles de l'équation de Vlasov[END_REF][START_REF] Mouton | Two-Scale semi-Lagrangian simulation of a charged particules beam in a periodic focusing channel[END_REF] developped a Two-Scale Semi Lagrangian Method for beam and plasma applications.

We consider the following model, valid for short-term dynamics of dunes, built and studied in [START_REF] Faye | Singularly perturbed degenerated parabolic equations and application to seabed morphodynamics in tided environment[END_REF]:

   ∂z ∂t - 1 ∇ • (A ∇z ) = 1 ∇ • C , z |t=0 = z 0 , (1.1)
where z = z (t, x) is the dimensionless seabed altitude. For a given T, t ∈ (0, T ) stands for the dimensionless time and x ∈ T 2 , T 2 being the two dimensional torus R 2 /Z 2 , stands for the dimensionless position and A , C are given by A (t, x) = A (t, x) + A 1 (t, x),

and

C (t, x) = C (t, x) + C 1 (t, x), (1.3)
where, for three positive constants a, b and c, U and M are the dimensionless water velocity and height.

A (t, x) = A(t
The small parameter involved in the model is the ratio between the main tide period 1 ω = 13 hours and an observation time which is about three months i.e. = 1 t ω = 1 200 . The following hypotheses on ga, gc, U and M given in (1.8) and (1.9) are technical assumptions and are needed to prove Theorem 1.1. Functions ga and gc are regular functions on R + and satisfy

                   ga ≥ gc ≥ 0, gc(0) = g c (0) = 0, ∃ d ≥ 0, sup u∈R + |ga(u)| + sup u∈R + |g a (u)| ≤ d, sup u∈R + |gc(u)| + sup u∈R + |g c (u)| ≤ d, ∃ U thr ≥ 0, ∃ G thr > 0, such that u ≥ U thr =⇒ ga(u) ≥ G thr . (1.8)
Functions U and M are regular and satisfy:

                                                 θ -→ (U , M) is periodic of period 1, |U |, | ∂U ∂t |, | ∂U ∂θ |, |∇U |, |M|, | ∂M ∂t |, | ∂M ∂θ |, |∇M| are bounded by d, ∀ (t, θ, x) ∈ R + × R × T 2 , |U (t, θ, x)| ≤ U thr =⇒ ∂U ∂t (t, θ, x) = 0, ∇U (t, θ, x) = 0, ∂M ∂t (t, θ, x) = 0, and ∇M(t, θ, x) = 0 , ∃ θα < θω ∈ [0, 1] such that ∀ θ ∈ [θα, θω] =⇒ |U (t, θ, x)| ≥ U thr .
(1.9)

To develop the Two-Scale Numerical Method, we use that in [START_REF] Faye | Singularly perturbed degenerated parabolic equations and application to seabed morphodynamics in tided environment[END_REF] we proved that under assumptions (1.8) and (1.9) the solution z of (1.1) exists, is unique and moreover asymptotically behaves, as → 0, the way given by the following theorem.

Theorem 1.1. Under assumptions (1.8) and (1.9), for any T, not depending on , the sequence (z ) of solutions to (1.1), with coefficients given by (1.2) coupled with (1.4) and (1.3), (1.5) and (1.6), Two-Scale converges to the profile

Z ∈ L ∞ ([0, T ], L ∞ # (R, L 2 (T 2 ))) solution to ∂Z ∂θ -∇ • ( A∇Z) = ∇ • C, (1.10)
where A and C are given by

A(t, θ, x) = a ga(|U (t, θ, x)|) and C(t, θ, x) = c gc(|U (t, θ, x)|) U (t, θ, x) |U (t, θ, x)| . (1.11)
Futhermore, if the supplementary assumption

U thr = 0, (1.12)
is done, we have

A(t, θ, x) ≥ G thr for any t, θ, x ∈ [0, T ] × R × T 2 , (1.13)
and, defining Z = Z (t, x) = Z(t, t , x), the following estimate holds for z -

Z z -Z L ∞ ([0,T ),L 2 (T 2 )) ≤ α, (1.14)
where α is a constant not depending on .

Because of assumptions (1.8) and (1.9),

A, C, A 1 , C 1 , A , A 1 , C
, and C 1 are regular and bounded.

(1.15)

2. Two-Scale Numerical Method Building. In this section, we develop the Two-Scale Numerical Method in order to approach the solution z of (1.1). The idea is to get a good approximation of z (t, x) seeing Theorem 1.1 content as z (t, x) ∼ Z(t, t , x).

The strategy is to consider a Fourier expansion of Z solution to (1.10). In this equation, t is only a parameter. The Fourier expansion of Z is given as follows:

Z(t, θ, x) = l,m,n Z l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) , (2.1)
where Z l,m,n (t), l = 0, 1, 2, . . ., m = 0, 1, 2, . . ., n = 0, 1, 2, . . . , are the unknown complex coefficients of the Fourier expansion of Z. Using (2.1), the Fourier expansion of ∂Z ∂θ is given by

∂Z ∂θ (t, θ, x) = l,m,n 2iπ l Z l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) . (2.2)
To obtain the system satisfied by the Fourier expansion (2.

+ l,m,n A l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) l,m,n 4π 2 (m 2 + n 2 )Z l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) = l,m,n C l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) , (2.8)
which gives after identification, the following algebraic system for (Z l,m,n ):

2iπ l Z l,m,n (t) - i,j,k 2iπ A grad i,j,k (t) • m -j n -k Z l-i,m-j,n-k (t) +4π 2 i,j,k A i,j,k (t)((m -j) 2 + (n -k) 2 )Z l-i,m-j,n-k (t) = C l,m,n (t).
(2.9)

In formula (2.1), the integers m, n and l vary from -∞ to +∞. But in practice, we will consider the truncated Fourier series of order P ∈ N defined by

Z P (t, θ, x) = 0≤l≤P,0≤m≤P,0≤n≤P
Z l,m,n (t) e 2iπ(lθ+mx 1 +nx 2 ) .

(2.10) Using (2.10), formula (2.9) becomes:

2iπ l Z l,m,n (t) - 0≤i≤P, 1≤j≤P, 0≤k≤P 2iπ A grad i,j,k (t) • m -j n -k Z l-i,m-j,n-k (t) +4π 2 0≤i≤P, 0≤j≤P, 0≤k≤P A i,j,k (t)((m -j) 2 + (n -k) 2 )Z l-i,m-j,n-k (t) = C l,m,n (t).
(2.11)

3. Convergence result.

Proof. of Theorem 1.1. For self-containedness, we recall the proof of Theorem 1.1. Firstly, we obtain an estimate leading to that z is bounded in L ∞ ((0, T ); L 2 (T 2 )). Secondly, defining test function ψ (t, x) = ψ(t, t , x) for any ψ(t, θ, x), regular with a compact support over [0, T ) × T 2 and 1-periodic in θ, multiplying (1.1) by ψ and integrating over [0, T ) × T 2 gives

T 2 T 0 ∂z ∂t ψ dtdx - 1 T 2 T 0 ∇ • (A ∇z )ψ dtdx = 1 T 2 T 0 ∇ • C ψ dtdx. (3.1)
Then integrating by parts in the first integral over [0, T ) and using the Green formula in T 2 in the second integral we have

- T 2 z 0 (x)ψ(0, 0, x)dx - T 2 T 0 ∂ψ ∂t z dtdx + 1 T 2 T 0 A ∇z ∇ψ dtdx = 1 T 2 T 0 ∇ • C ψ dtdx. (3.2)
Again using the Green formula in the third integral we obtain

- T 2 z 0 (x)ψ(0, 0, x) dx - T 2 T 0 ∂ψ ∂t z dtdx - 1 T 2 T 0 z ∇ • (A ∇ψ ) dtdx = 1 T 2 T 0 ∇ • C ψ dtdx. (3.3) But ∂ψ ∂t = ∂ψ ∂t + 1 ∂ψ ∂θ , (3.4)
where

∂ψ ∂t (t, x) = ∂ψ ∂t (t, t , x) and ∂ψ ∂θ (t, x) = ∂ψ ∂θ (t, t , x), (3.5)
then we have

T 2 T 0 z ∂ψ ∂t + 1 ∂ψ ∂θ + 1 ∇ • (A ∇ψ ) dxdt + 1 T 2 T 0 ∇ • C ψ dtdx = - T 2 z 0 (x)ψ(0, 0, x) dx. (3.6)
Using the Two-Scale convergence due to Nguetseng [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] and Allaire [START_REF] Allaire | Homogenization and Two-Scale convergence[END_REF] (see also Frénod Raviart and Sonnendrücker [START_REF] Frénod | Asymptotic expansion of the Vlasov equation in a large external magnetic field[END_REF]), since z is bounded in L ∞ ([0, T ), L 2 (T 2 )), there exists a profile Z(t, θ, x), periodic of period 1 with respect to θ, such that for all ψ(t, θ, x), regular with a compact support with respect to (t, x) and 1-periodic with respect to θ, we have

T 2 T 0 z ψ dtdx -→ T 2 T 0 1 0
Zψ dθdtdx, as tends to zero, (

for a subsequence extracted from (z ). Multiplying (3.6) by , passing to the limit as → 0 and using (3.7) we have

T 2 T 0 1 0 Z ∂ψ ∂θ dθdtdx + lim →0 T 2 T 0 z ∇ • (A ∇ψ ) dtdx = lim →0 T 2 T 0 C • ∇ψ dtdx, (3.8)
for an extracted subsequence. As A and C are bounded and ψ is a regular function, A ∇ψ and ∇ψ can be considered as test functions. Using (3.7) we have

T 2 T 0 z ∇ • (A ∇ψ )dtdx -→ T 2 T 0 1 0 Z∇ • ( A∇ψ) dθdtdx, (3.9)
and

T 2 T 0 C • ∇ψ dtdx Two-Scale converges to T 2 T 0 1 0 C • ∇ψ dθdtdx. ( 3
∂z ∂t - 1 ∇ • ( A ∇z ) = 1 ∇ • C + ∇ • ( A 1 ∇z ) + ∇ • C 1 . (3.11)
For Z , we have

∂Z ∂t = ∂Z ∂t + 1 ∂Z ∂θ , (3.12)
where

∂Z ∂t (t, x) = ∂Z ∂t (t, t , x) and ∂Z ∂θ (t, x) = ∂Z ∂θ (t, t , x). (3.13) Using (1.10), Z is solution to ∂Z ∂t - 1 ∇ • A ∇Z = 1 ∇ • C + ∂Z ∂t . ( 3
∂(z -Z ) ∂t - 1 ∇ • A ∇(z -Z ) = ∇ • C 1 + ∂Z ∂t + ∇ • ( A 1 ∇z ). (3.15)
Multiplying equation (3.15) by 1 and using the fact that z = z -Z + Z in the right hand side of equation

(3.15), z -Z is solution to: ∂ z -Z ∂t - 1 ∇ • ( A + A 1 )∇( z -Z ) = 1 ∇ • C 1 + ( ∂Z ∂t ) + ∇ • ( A 1 ∇Z ) . (3.16)
Our aim here is to prove that z -Z is bounded by a constant α not depending on . For this let us use that A , A 1 , C and C 1 are regular and bounded coefficients (see (1.15)) and that A ≥ G thr (see

(1.13)). Hence, ∇ • C 1 is bounded, ∇ • ( A 1 ∇Z ) is also bounded. Since Z is solution to (3.14), ∂Z ∂t satisfies the following equation ∂ ∂Z ∂t ∂θ -∇ • A∇ ∂Z ∂t = ∂∇ • C ∂t + ∇ • ∂ A ∂t ∇Z .
(3.17) Equation (3.17) is linear with regular and bounded coefficients. Using a result of Ladyzenskaja, Solonnikov and Ural'Ceva [START_REF] Ladyzenskaja | Linear and quasi-linear equations of parabolic type[END_REF], ∂Z ∂t is regular and bounded and so the coefficients of equations (3.16) are regular and bounded. Then, using the same arguments as in the proof of Theorem 1.1 in [START_REF] Faye | Singularly perturbed degenerated parabolic equations and application to seabed morphodynamics in tided environment[END_REF] we obtain that z -Z is bounded.

To determine the value of the constant α, we proceed in the same way as in the proof of Theorem 3.16 of [START_REF] Faye | Singularly perturbed degenerated parabolic equations and application to seabed morphodynamics in tided environment[END_REF]. Since

the coefficients A , A 1 , C and C 1 , ∇ • C 1 , ∇ • ( A 1 ∇Z
), and ∂Z ∂t are bounded by constants, let β denotes the maximum between all these constants. Then we use the same argument as in the proof of Theorems 1.1 and 3.16 and we get:

z -Z L ∞ ([0,T ),L 2 (T 2 )) ≤ z 0 (•) -Z(0, 0, •) 2 β + β 3 G thr + 2β T. (3.18)
Theorem 3.1. Let be a positive real, z be the solution to (1.1), Z P be the truncated Fourier series (defined by (2.10)) of Z solution to (1.10) and Z P defined by Z P (t, x) = Z P (t, t , x). Then, under assumptions (1.8), (1.9) and (1.12), z -Z P satisfies the following estimate:

z -Z P L ∞ ([0,T ),L 2 (T 2 )) ≤ z 0 (•) -Z(0, 0, •) 2 β + β 3 G thr + 2β T + f (P ), (3.19)
where f is a non-negative function of P not depending on and satisfying lim P →+∞ f (P ) = 0.

Proof. We can write :

z -Z P L ∞ ([0,T ),L 2 (T 2 )) = z -Z + Z -Z p L ∞ ([0,T ),L 2 (T 2 )) ≤ z -Z L ∞ ([0,T ),L 2 (T 2 )) + Z -Z p L ∞ ([0,T ),L 2 (T 2 )) . ( 3
z -Z L ∞ ([0,T ),L 2 (T 2 )) ≤ z 0 (•) -Z(0, 0, •) 2 β + β 3 G thr + 2βT. (3.21)
For the second term of (3.20), using classical results of Fourier series theory, since Z -Z P is nothing but the rest of the Fourier series of order P of Z and since Z is regular (because it is the solution of (1.10) which has regular coefficients), the non-negative function f satisfying lim P →+∞ f (P ) = 0 such that

Z -Zp L ∞ ([0,T ],L ∞ # (R,L 2 (T 2 ))) ≤ f (P ), (3.22)
exists. From this last inequality,

Z -Z p L ∞ ([0,T ),L 2 (T 2 )) ≤ f (P ), (3.23

Reference solution.

Having Fourier coefficients of Z on hand, we will do the same for function z (t, x) solution to (1.1) in order to compare it to the profile Z for a given , in a fixed time. The Fourier expansion of z is given by

z (t, x 1 , x 2 ) = m,n zm,n(t) e 2πi(mx 1 +nx 2 ) , (4.1)
where m = 0, 1, 2, . . . and n = 0, 1, 2, . . . , then the Fourier expansion of ∂z ∂t is

∂z ∂t = m,n żm,n(t) e 2πi(mx 1 +nx 2 ) . (4.2)
Using the same idea as in the Fourier expansion of Z, we obtain the following infinite system of Ordinary Differential Equations

∂zm,n ∂t (t) - 1 i,j 2iπA grad i,j (t) • m -i n -j z m-i,n-j (t) + 1 4π 2 i,j A i,j (t)((m -i) 2 + (n -j) 2 )z m-i,n-j (t) = 1 Cm,n(t), (4.3)
where A grad i,j (t), A i,j (t) and Cm,n(t) are respectively the Fourier coefficients of ∇A , A and ∇ • C . In the same way, the truncated Fourier series of order P ∈ N of z is given by

z P (t, x 1 , x 2 ) = P m,n=0 zm,n(t) e 2πi(mx 1 +nx 2 ) , (4.4)
which gives from (4.3) the following system Ordinary Differential Equations ∂zm,n ∂t (t) -

1 P i,j=0 2iπA grad i,j (t) • m -i n -j z m-i,n-j (t) + 1 4π 2 P i,j=0 A i,j (t)((m -i) 2 + (n -j) 2 )z m-i,n-j (t) = 1 Cm,n(t). (4.5)
In (4.5), we will use an initial condition zm,n(0, x). To solve (4.5) we use, for the discretization in time, a Runge-Kutta method (ode45).

Comparison Two-Scale Numerical Solution and reference solution.

In this paragraph, we consider the truncated solution z P (t, x 1 , x 2 ) and Z P (t, t , x 1 , x 2 ). The objective here is to compare for a fixed and a given time, the quantity |z P (t, x 1 , x 2 ) -Z P (t, t , x 1 , x 2 )| when the water velocity U is given.

4.2.1. Comparisons of z P (t, x) and Z P (t, t , x) with U given by (4.6). For the numerical simulations, concerning z , we take z 0 (x 1 , x 2 ) = cos 2πx 1 + cos 4πx 1 and z 0 (x 1 , x 2 ) = Z(0, 0, x 1 , x 2 ). In what concerns the water velocity field, we consider the function

U (t, θ, x 1 , x 2 ) = sin πx 1 sin 2πθ e 1 , (4.6
) where e 1 and e 2 are respectively the first and the second vector of the canonical basis of R 2 and x 1 , x 2 are the first and the second components of x.

In Figure 1 , we can see the space distribution of the first component of the velocity U for a given time t = 1 and for various values of θ: 0.3, 0.55, and 0.7. In Figure 2, we see, for a fixed point x = (x 1 , x 2 ), how the water velocity U (θ) evolves with respect to θ. In Figure 3, the θ-evolution of A(θ) is also given in various points (x 1 , x 2 ) ∈ R 2 . In this paragraph, we present numerical simulations in order to validate the Two-Scale convergence presented in Theorem 1.1. For a given , we compare Z P (t, t , x), where Z P is the Fourier expansion of order P of the solution to (1.10) and z P (t, x) the Fourier expansion of order P of the solution to the reference problem. The simulations presented are given for P = 4. The calculation of z P (t, x) implies knowledge of z 0 (x). For an initial condition z 0 (x) well prepared and equal to Z(0, 0, x), we obtain the results of Figure 4 and we remark that the results obtained are the same for z P (t, x) and Z P (t, t , x). t,•) and Z P (t, t , •), P = 4, at time t = 1, = 0.001, when U is given by (4.6) and when z 0 (•) = Z(0, 0, •). On the left z P (t, •), on the right Z P (t, t , •).

In practice, the solution Z P , P ∈ N evolves according to P. For the simulations, we made the value of the integer P vary and we saw that this variation is very low from P ≥ 6. To better show that Z P (t, t , x 1 , x 2 ) is close to the reference solution z P (t, x 1 , x 2 ), we plot and compare Z P (t, t , x 1 , 0) and z P (t, x 1 , 0), at different times t. In these comparisons the initial condition z 0 (x 1 , x 2 ) = cos 2πx 1 + cos 4πx 1 is different from Z(0, 0, x 1 , x 2 ). The results are shown in Figure 5 and Figure 6. We see in these figures that the solution z P (t, x) get closer and closer to Z P (t, t , x) with time of order . Z(t,t/ε,x 1 ,0) and z ε (t,x 1 ,0), t=2ε, ε=0.001

z ε Z Figure 5.
Comparison of z P (t, x 1 , 0) and Z P (t, t , x 1 , 0)), P = 4. On the left t = 0, in the middle t = and t = 2 on the right, = 0.001. Z(t,t/ε,x 1 ,0) and z ε (t,x 1 ,0), t=ε, ε=0.01

z ε Z 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -4 -3 -2 -1 0 1 2 3 4 x 10 -3
Z(t,t/ε,x 1 ,0) and z ε (t,x 1 ,0), t=2ε, ε=0.01

z ε Z Figure 6.
Comparison of z P (t, x 1 , 0) and Z P (t, t , x 1 , 0)). On the left t = 0, in the middle t = and t = 2 on the right, = 0.01.

So we can see from these figures that the solution Z of the Two-Scale limit problem is such that Z(t, t , •, •) is close to the solution z (t, •, •) of the reference problem. In the presently considered case where the initial condition for z is not Z(0, 0, •, •), we saw in Figure 5 and Figure 6 that z P tends to reach a steady state. This steady state is an oscillatory one in the sense that for large t, z P (t, •, •) behaves like Z P (t, t , •, •). This is illustrated by Figure 7 where z P (t, x 1 , 0) and Z P (t, t , x 1 , 0) are given for various value of t in a period of lenght . More precisely, in this figure we see that within a period of time of lenght , z P (t, •, •) and Z P (t, t , •, •) do not glue together completly. Nevertheless, despite this phenomenon which is linked with the fact that the Two-Scale approximation of z (t, •, •) by Z(t, t , •, •) is only of order 1 in , the two solutions re-glue well together at the end of the period. 4.2.2. Comparisons of z (t, x) and Z(t, t , x) with U is given by (4.7). In this subsection, we do the same as in the precedent one, but when the velocity fields U given by (4.7). The results are all identical to the precedent one i.e. the Two-Scale limit Z P (t, t , x 1 , x 2 ) is very close to the solution z P (t, x 1 , x 2 ) to the reference problem when P = 4. The initial condition z 0 (x 1 , x 2 ) = Z(0, 0, x 1 , x 2 ) and is the same as in the subsection 4.2.1. The results are given for = 0.1 and = 0.005 and for various time t. We notice that z comes very close to Z(t, t , x 1 , x 2 ) when is very small. We begin by giving the space distribution of U at various time and the θ-evolution of U and A. The second velocity fields is given by

U (t, θ, x 1 , x 2 ) = U (t, θ, x) =                                                            0 in [0, θ 1 ], θ-θ 1 θ 2 -θ 1 U thr e 2 in [θ 1 , θ 2 ], U thr e 2 + φ( θ-θ 2 θ 3 -θ 2 )ψ(t, x) in [θ 2 , θ 3 ], θ-θ 3 θ 4 -θ 3 U thr e 2 in [θ 3 , θ 4 ], 0 in [θ 4 , θ 5 ], θ-θ 5 θ 6 -θ 5 U thr e 2 in [θ 5 , θ 6 ], -U thr e 2 -φ( θ-θ 6 θ 7 -θ 6 )ψ(t, x) in [θ 6 , θ 7 ], -θ-θ 7 θ 8 -θ 7 U thr e 2 in [θ 7 , θ 8 ], 0 in [θ 8 , 1], (4.7)
where U thr > 0, φ is a regular positive function satisfying φ(s) = s(1 -s) and ψ(t, x 1 ) = (1 + sin π 30 t)(U thr e 2 + 1 10 (1 + sin 2πx 1 )e 1 ), θ i = i+1 10 , i = 1, . . . , 8. The θ-evolution of U , given by (4.7), is given in Figure 9 for various position in [0, 1] 2 . Function ga(u) = gc(u) = |u| 3 , a = c = 1 and M(t, θ, x) = 0 which yields a θ-evolution of A(θ) which is drawn for various positions in Figure 10. Using this, we compute Z P (t, t , x 1 , x 2 ) and z P (t, x) for P = 4. To compute z P (t, x) we take z 0 (x 1 , x 2 ) = cos 2πx 1 + cos 4πx 1 which is not Z(0, 0, x 1 , x 2 ). First we study the errors Z P (t, t , x 1 , x 2 ) -z P (t, x) at t = 1. This quantity decreases when decreases as illustrated in the following tabular. : Errors norm Z P (t, t , x 1 , x 2 ) -z P (t, x 1 , x 2 ), P = (4, 4), P = (4, 4, 4), t = 1.

The results given in this table show that, at time t = 1, z (t, x) is closer to Z(t, t , x) when is very small. These results validate the results obtained in Theorem 1.1. In Figures 11 and12, we present simulations at times t = 0.75 and t = 0.775. We see that Z P (t, t , x 1 , x 2 ) is close to z P (t, x 1 , x 2 ). The numerical results shown in these figures are made with = 0.1. In Figure 13 and 14, we do the same but for = 0.005. The numerical results show that z P (t, x) is also very close to Z P (t, t , x 1 , x 2 ). We remark that for = 0.1 and = 0.005, the solution z P (t, x) is very close to Z P (t, t , x). But the approximation z P (t, x) ∼ Z P (t, t , x) is very good when is very small. To show that z P is very close to Z P , we construct the same figures as previously but in dimension 2 with = 0.005 i.e. we construct z P (t, x 1 , 0) and Z P (t, t , x 1 , 0) for = 0.005 at time t = 0.775. This is given in Figure 15. Figure 11. Comparison of z P (t, x 1 , x 2 ) and Z P (t, t , x 1 , x 2 ), P = 4; t = 0.75, = 0.1, z 0 (x 1 , x 2 ) = cos 2πx 1 + cos 4πx 1 . On the left Z P (t, t , x 1 , x 2 ), on the right z P (t, x 1 , x 2 ).

IBRAHIMA FAYE, EMMANUEL FR ÉNOD AND DIARAF SECK

Figure 13. Comparison of z P (t, x 1 , x 2 ) and Z P (t, t , x 1 , x 2 ), P = 4; t = 0.75, = 0.005, z 0 (x 1 , x 2 ) = cos 2πx 1 + cos 4πx 1 . On the left Z P (t, t , x 1 , x 2 ), on the right z P (t, x 1 , x 2 ). of z P (t, x 1 , 0) and Z P (t, t , x 1 , 0), t = 0.775, = 0.005. On the left Z P (t, t , x 1 , 0), on the right z P (t, x 1 , 0).

The results in Figure 16 show that Z P and z P have the same behavior in the same period and Z P is very close to z P . We also notice that, despite the small shifts that occur during a period, the two solutions glue together.

. 10 )

 10 Passing to the limit as → 0 we obtain from (3.8) a weak formulation of the equation (1.10) satisfied by Z. Using (1.2) and (1.3) equation (1.1) becomes

  .14) Formulas (3.11) and (3.14) give

  .20) Using (3.18), the first term in the right hand side of (3.20) is bounded by

4 .

 4 ) follows and coupling this with (3.21) and (3.20) gives inequality(3.19). Numerical illustration of Theorem 3.1.

Figure 1 .

 1 Figure 1. Space distribution of the first component of U(1, 0.3, (x 1 , x 2 )), U(1, 0.55, (x 1 , x 2 )) and U(1, 0.7, (x 1 , x 2 )) when U is given by (4.6).

Figure 2 .Figure 3 .

 23 Figure 2. θ-evolution of U(θ, (1/2, 0)) and U(θ, (1/4, 0)) when U is given by (4.6)

Figure 4 .

 4 Figure 4. Comparison of z P (t, •) and Z P (t, t , •), P = 4, at time t = 1, = 0.001, when U is given by (4.6) and when z 0 (•) = Z(0, 0, •). On the left z P (t, •), on the right Z P (t, t , •).

3 Z

 3 (t,t/ε,x 1 ,0) and z ε (t,x 1 ,0), t=ε, ε=0.001

Figure 7 .

 7 Figure7. Evolution of Z P (t, t , x 1 , 0) in the top and z P (t, x 1 , 0) in the bottom, t = 1 + n 4 , n = 0, 1, 2, 3.

Figure 8 .

 8 Figure 8. Space distribution of the first component of U(1, 0.25, (x 1 , x 2 )), U(1, 0.275, (x 1 , x 2 )) and U(1, 0.75, (x 1 , x 2 )) when U is given by (4.7).

Figure 9 .

 9 Figure 9. θ-evolution of U(1, θ, (1, 0)), U(1, θ, (4, 0)) and U(1, θ, (1/3, 1/3)) when U is given by (4.7).

Figure 10 .

 10 Figure 10. θ-evolution of A(1, θ, (1, 0)), A(1, θ, (4, 0)) and A(1, θ, (1/3, 1/3)) when U is given by (4.7).

  Figure15. Comparison of z P (t, x 1 , 0) and Z P (t, t , x 1 , 0), t = 0.775, = 0.005. On the left Z P (t, t , x 1 , 0), on the right z P (t, x 1 , 0).

  Figure15. Comparison of z P (t, x 1 , 0) and Z P (t, t , x 1 , 0), t = 0.775, = 0.005. On the left Z P (t, t , x 1 , 0), on the right z P (t, x 1 , 0).

5 Z

 5 (t,t/ε,x 1 ,0), t=1+ε/2, ,t/ε,x 1 ,0), t=1+3ε/4, ε=0.005 1 z ε (t,x 1 ,0), t=1,ε=0.005 0.4

z ε (t,x 1 ,0), t=1+ε/4, ε=0.005


Christophe Lecoutre

Chavalit Likitvivatanavong

Roland H C Yap

A Path-Optimal GAC Algorithm for Table Constraints

A path-optimal GAC algorithm for table constraints

Introduction

Constraint propagation, which consists in calling iteratively filtering algorithms associated with constraints, is one of the most attractive features of Constraint Programming (CP). The various levels of possible filtering are typically described by properties of constraint networks called consistencies. Generalized Arc Consistency (GAC), which corresponds to the highest filtering level of variable domains when constraints are considered independently, is such a property. GAC is one very important lever to solve efficiently Constraint Satisfaction Problems (CSPs) because it identifies many more inconsistent values than limited consistency forms (e.g., those defined on domain bounds or dependent on the number of assigned variables), reducing search space considerably as a consequence.

Table constraints are defined in extension by explicitly listing all permitted combinations of values (positive tables) or all forbidden ones (negative tables). Table constraints naturally arise in many application areas such as configuration and databases, and besides they can be viewed as a universal mechanism for representing any constraints. Classical filtering algorithms (e.g. [START_REF] Bessiere | Arc consistency for general constraint networks: preliminary results[END_REF][START_REF] Bessière | An Optimal Coarse-Grained Arc Consistency Algorithm[END_REF][START_REF] Lecoutre | Generalized arc consistency for positive table constraints[END_REF][START_REF] Lhomme | A fast arc consistency algorithm for n-ary constraints[END_REF]) generally do not alter table constraints during backtrack search. Recent developments, however, suggested handling tables directly, which leads to faster algorithms [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF][START_REF] Ullmann | Partition search for non-binary constraint satisfaction[END_REF]. Alternatively, specially-constructed intermediate structures such as tries [START_REF] Gent | Data structures for generalised arc consistency for extensional constraints[END_REF] or Multi-valued Decision Diagrams (MDDs) [START_REF] Cheng | An MDD-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints[END_REF] have been proposed. In any case, the search space gets smaller as the tables or their equivalent structures are reduced.

Most GAC algorithms follow the same pattern: a value is proved to be consistent by producing a valid tuple containing that value (in the case of positive table constraints) or by producing evidence from auxiliary structures (a path from top to bottom in the case of MDDs). This is usually performed by traversing these structures and running tests on each chunk. Reducing the amount of traversal has long been the focus of a many works resulting in many optimization techniques.

Simple Tabular Reduction (STR) and its improvements [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF][START_REF] Ullmann | Partition search for non-binary constraint satisfaction[END_REF] fall into this category and have been shown to be among the best GAC algorithms for positive table constraints. The main idea of STR is to remove invalid tuples systematically from tables immediately.

In this paper, we introduce a new GAC algorithm called STR3. It is based on the same principle as STR but employs a different representation of table constraints. Indeed, domain values are the focal points: a set of tuple identifiers (row numbers) is associated with each value, indicating the different tuples (rows) where the value appears in. Figure 1 gives a ternary constraint example showing the standard table and our equivalent representation, i.e. both representations have the same semantics. We also introduce a novel technique for maintaining valid supports of domain values. Two kinds of data structures are involved: the first one for partitioning each set of rows (as in Figure 1b) into invalid and untested areas, and the other for keeping track of which value depends on which row as its valid support. The former requires maintenance while the latter is backtrack-free. We show that the synergy between these two structures leads to greater efficiency and provides us with extra valid supports with no additional cost when the search backtracks.

It is worthwhile to note that algorithms such as STR2 [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF] or mddc [START_REF] Cheng | An MDD-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints[END_REF] may suffer from multiple traversals of the same region when they are invoked successively. On the other hand, STR3, similarly to GAC4, is path-optimal: each element of a table is examined at most once along any path of the search tree. Importantly, STR3 is designed to be used directly during search, where maintaining consistency is needed and cost of backtracking should be minimized. GAC4, however, is to enforce consistency in a standalone context. While it is possible to convert GAC4 to MGAC4 (Maintaining GAC4), this is not simple and not been treated anywhere as far as we know. STR3 also allows some freedom during initialization, e.g. using STR2.

Our experiments show that STR3 is rather complementary to STR2. Where simple tabular reduction can eliminate so many tuples from the tables that they become largely empty, STR2 is faster than STR3. STR3, by contrast, outperforms STR2 when the average size of the tables during search is not too low. A finite constraint network P is a pair (X , C) where X is a finite set of n variables and C is a finite set of e constraints. D(X) represents the domain of X ∈ X , i.e., the set of values that can be assigned to X. During search, D c (X) denotes the current domain of X. If a ∈ D c (X), we say that a is present in D(X); otherwise a is absent from D(X). We use (X, a) to denote the value a ∈ D(X) (or simply a when the context is clear). Each constraint C ∈ C involves an ordered subset of variables of X , called its scope and denoted by scp(C), and a relation denoted by rel(C). For any r-ary constraint C with scp(C) = {X1, . . . , Xr}, rel(C) ⊆ r i=1 D(Xi) denotes the set of satisfying combinations of values for the variables in scp(C). A solution to a constraint network is an assignment of a value to each variable such that every constraint is satisfied. A constraint network is satisfiable iff at least one solution exists.

X Y Z 1 a f l 2 b f m 3 e g m 4 a f m 5 b g o 6 a h o 7 d h o 8 b i n 9 c j k (a) Standard table X Y Z a {1,4,6} f {1,2,4} k {9} b {2,5,8} g {3,5} l {1} c {9} h {6,7} m {2,3,4} d {7} i {8} n {8} e {3} j {9} o {5,6,7} (
For any r-tuple t = (a1, . . . , ar

) ∈ rel(C) such that scp(C) = {X1, . . . , Xr}, t[Xi] denotes ai. A tuple t ∈ rel(C) is valid iff t[X] ∈ D c (X) for each X ∈ scp(C). A tuple t ∈ rel(C) is a support of (X, a) on C iff t[X] = a. A value (X, a) is Generalized Arc-Consistent (GAC) on a constraint C involving X iff there exists a valid support t of (X, a) on C; (X, a) is GAC iff it is GAC on every constraint C involving X. A variable X is GAC iff D c (X) = ∅ and (X, a) is GAC for each a ∈ D c (X). A constraint network is GAC iff each of its variables is GAC.
We assume a total ordering for every and define pos(C, t) to be the position of the tuple t in that ordering. Given C ∈ C, X ∈ scp(C), and a ∈ D(X), we define row(C, X, a) to be the set {pos(C, t) | t ∈ rel(C) ∧ t[X] = a} -the set of indices (called rows) to each support of (X, a) on C. We say that row k

of constraint C is (in)valid iff the tuple t such that pos(C, t) = k is (in)valid. The set of invalid rows of C is denoted by inv(C) = {pos(C, t) | t ∈ rel(C) ∧ t is invalid }.

STR3

In this section, we introduce STR3, an algorithm based on simple tabular reduction for enforcing GAC on positive table constraints. A central operation for a GAC algorithm is to check whether a given value a ∈ D c (X) has a valid support on a table constraint C, which usually entails going through every tuple t in rel(C) and testing if t[X] = a and t[Y ] ∈ D c (Y ) for every Y ∈ scp(C) \ X. Checking each row (tuple) and column (tuple value) in this manner has been the cornerstone of many GAC algorithms. Optimizations efforts are often concentrated on reducing the amount of traversal, typically by skipping over irrelevant rows or columns of the tables [START_REF] Gent | Data structures for generalised arc consistency for extensional constraints[END_REF][START_REF] Lecoutre | Generalized arc consistency for positive table constraints[END_REF][START_REF] Lhomme | A fast arc consistency algorithm for n-ary constraints[END_REF].

In Simple Tabular Reduction (STR) [START_REF] Ullmann | Partition search for non-binary constraint satisfaction[END_REF], tables are dynamically maintained so that they contain valid tuples only. STR2 [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF] features two improvements over standard STR. When a tuple t is being inspected, t[X] is skipped over if it is known that every single value of D c (X) is already supported. Second, there is no need to check whether t[X] ∈ D c (X) if there has been no change to the domain of X since the last time STR2 was called (throughout this paper we actually refer to the variant called STR2+ in [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF]).

We propose STR3 which uses the same underlying principle as STR and STR2: we no longer examine a tuple once it has been recognized as invalid. But unlike STR and STR2, we do not explicitly discard the tuple from the table. Rather, we accomplish this indirectly by partitioning off invalid tuples in a different but equivalent representation. This allows us to avoid duplicated effort in re-establishing the consistency of values across the search tree as it commonly happens with conventional GAC algorithms. STR3, which is repair-based and fine-grained, works as follows:

• Every time a value (X, a) is deleted, STR3 is invoked for every constraint C involving X. row(C, X, a) is then merged into inv(C). • Because row(C, X, a) contains all supports of (X, a) on C, to verify whether a domain value (X, a) is GAC on C, STR3 needs only to test if row(C, X, a) \ inv(C) = ∅. • Each set row(C, X, a) is associated with a separator (can be thought of as a cursor) which partitions the set into two areas: one containing the tuples known to be invalid, the other the tuples yet to explored. The separator moves sequentially from one end of row(C, X, a) to the other in a fixed direction. As search progresses, the invalid area grows until it encompasses the whole set, at which point (X, a) has been proven to be not GAC on C.

• To check the validity of a row k in row(C, X, a), STR3 tests if k ∈ inv(C).
• Each row of the table is associated with a list of domain values, indicating that this row is a valid support for these values. Whenever the row becomes invalid, STR3 must look for a new valid support for every value in the list.

Implementation

Detailed operations of STR3 is given as pseudo-code in Figure 2. We first explain the data structures used in the algorithm:

• row(C,X,a) is implemented as an array. row(C,X,a).size is the number of elements of this array while row(C,X,a).curr is a number ranging from 0 to size -1, called the separator of row(C,X,a), indicating that row(C,X,a)[row(C,X,a).curr] corresponds to the last known valid support of (X, a) on C. For brevity, we shall use row(C,X,a

)[↑] to denote row(C,X,a)[row(C,X,a).curr].
The value of curr is maintained throughout the search. • inv(C) is implemented as a sparse set: inv(C).members gives the position of the last current element in inv(C) and inv(C).dense is the array containing all elements (see [START_REF] Briggs | An efficient representation for sparse sets[END_REF][START_REF] Cheng | An MDD-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints[END_REF] for details). In inv(C), we only need to keep at most a single copy of each value once the search starts. This is sound because the sets row(C, X, a) are fixed -only row(C, X, a).curr may change during search and must be restored when backtracking occurs. inv(C).members is also maintained throughout the search. • dep(C) is called the dependency list of C, implemented as an array of sets. dep(C) [k] is the set of values (X, a) such that the tuple in row k is a valid support of (X, a) on C; we say that (X, a) depends on row k. dep is not maintained during search.

Because STR3 can maintain GAC but does not establish it from scratch, a different GAC algorithm is needed before search (in the preprocessing stage). GACinit (Lines 1-6) is first called to remove all invalid tuples and to initialize all data structures. During search, STR3 (Lines 7-28) is called on a constraint C every time a value a is removed from the domain of a variable X ∈ scp(C). For each such value (X, a), every row in row(C, X, a) becomes invalid. STR3 then appends these rows to inv(C) if they are not already present (Line 9-11). Values that need new valid supports are later processed (Lines 14-27); we discuss this part of the algorithm in the next subsection. Upon backtracking, Functions restoreR and restoreI are called so as to restore values row(C, X, a).curr and inv(C).members through the use of the stacks stateR and stateI. Values are stored in these stacks at Lines 13 and 25 by calling Function save (Lines 29-31).

Synchronized Supports

Central to the implementation is the relationship between the separators and the dependency list. A present value (X, a) is GAC on C either because (X, a) ∈ dep(C) [v] for some row v / ∈ inv(C), or because row(C, X, a)[↑] = w for some row w / ∈ inv(C). Only one of the conditions is necessary for (X, a) to be GAC, and when both conditions are true, v does not have to be the same as w. We study the circumstances involving these two conditions and their values here.

When v = w, we say that the dependency list and the separators are synchronized at (X, a) (or that the supports of (X, a) are synchronized). When the search starts, GACinit initializes curr and dep so that each of them refers to an opposite end of row(C, X, a). Both are valid supports of (X, a), as any invalid row is removed at Line 2 during preprocessing. Only the elimination of row k would trigger the search for a new valid support for each (X, a)

∈ dep(C)[k].
When the separators and the dependency list are synchronized at a particular value, we are provided with a single valid support instead of two. In this case, the role of the dependency list is straightforward: it just mirrors what happens to the separators. When row k becomes invalid, we look for a new valid support for each value that depends on k according to the dependency list (Line 16). Potential supports are tested one by one against inv(C) (Lines 18-19). If no valid support is found the value is removed (Line 21). STR3 immediately fails when that value is the last one left in the domain (Line 22). Other-wise, the previous curr is recorded for backtrack purpose (Line 25) and the position where the support is found is set to be the value of curr (Line 26). dep(C) is updated accordingly at Line 27. The separators and the dependency list remain synchronized.

Given (X, a) ∈ dep(C)[k], when the search algorithm backtracks dep(C)[k] will hold on to its values while row(C, X, a)[↑] must revert back to its previous state if applicable. We may end up with a situation where the separators and the dependency list are no longer synchronized at (X, a). In such cases, tuples at row(C, X, a)[↑] and k are two distinct valid supports of (X, a) on C. We consider what happens inside STR3 when their validity later change: These relationships are summarized in Figure 3. The separators and the dependency list are comparable to watched literals [START_REF] Moskewisz | Chaff: Engineering an efficient SAT solver[END_REF] introduced for SAT. Significant differences are as follows. To begin with, dep is the only activation point, working as a primary valid support while curr serves as a possible backup; curr points to a supporting row that may or may not be valid. In contrast, there are always two watched literals for SAT, both functionally equivalent. curr is rigid and must maintain its value at all times while dep is not maintained. The two watched literals are unmaintained. dep and curr can be synchronized or unsynchronized depending on circumstances, in effect providing either a single support or two distinct supports whereas for SAT the two watched literals are always distinct where possible.

•

Related Works

STR3 is guided by deleted values, making it a fine-grained algorithm. Other fine-grained (G)AC algorithms have been proposed in the literature such as AC6 [START_REF] Bessiere | Arc consistency and arc consistency again[END_REF], AC7 [START_REF] Bessiere | Using constraint metaknowledge to reduce arc consistency computation[END_REF] and GAC4 [START_REF] Mohr | Good old discrete relaxation[END_REF]. All these algorithms use dependency lists. However, STR3 differs significantly from AC6, AC7, and GAC4 by the choice of the additional data struc-tures. The closest algorithm to STR3 is GAC4, which we will consider in more detail.

MGAC4 requires complicated management of dependency lists, which have to be implemented as doubly linked-lists along with additional structures to keep complexity cost down. Because a row position appears in more than one list, difficulty arises when backtrack occurs: MGAC4 has to be careful not to restore row positions that have been removed at shallower depths. This entails record keeping and thus increases overhead. STR3 avoids this problem by sequentially scanning the lists and cordoning off invalid members rather than performing random-access operations on any location. In a way, STR3 can be seen as a highly optimized version of MGAC4 (for which, we are not aware of any efficient implementation published in the literature) through the mechanism of simple tabular reduction.

Example

As an illustration, we consider the table constraint depicted in Figure 1a to demonstrate how the algorithm works. For each value (X, a), row(C, X, a) is given in Figure 1b. After GAC preprocessing, row(C, X, a).curr and dep(C) are initialized as shown in Fig- ure 4. For clarity, we synchronize the separators and the dependency list before the search starts, as opposed to the unsynchronized version in the actual code of GACinit. The symbol t indicates that the associated value is assigned as curr at time t.

Assume values h, i, and o are eliminated. During the execution of STR3, all rows involving these values are appended to inv(C). We update dep(C) [k] for each k ∈ inv(C). The result is shown in Figure 5. We denote the fact that inv(C).curr is assigned the value k at time t by placing ↑t at column k. A square box surrounding a domain value means that this value has been deleted.

Suppose now that we backtrack to t = 1. dep(C) is unaffected while row(C, X, a).curr and inv(C).members are rolled back to the ones at time t = 1. The result is shown in Figure 6. Next, suppose e, and n are eliminated. Rows 3 and 8 are added to inv(C). We will look at the changes to dep(C) in details (Figure 7). Value e in dep(C) [START_REF] Bessiere | Arc consistency for general constraint networks: preliminary results[END_REF] and value n in dep(C) [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF] are already deleted so they remain in their positions according to Line 16 of STR3. Value i is removed because it has no further valid support. We now consider values g and b. A valid support has to be found for g because it depends on row 3. Recall that g ∈ dep(C) [START_REF] Bessiere | Arc consistency for general constraint networks: preliminary results[END_REF] tells us that the tuple in row 3 is a valid support of g. As soon as row 3 becomes invalid, a new valid support must be found. Because row(C, Y, g)[↑] = 5 is a valid support, g is moved from dep(C) [START_REF] Bessiere | Arc consistency for general constraint networks: preliminary results[END_REF] to dep(C) [START_REF] Briggs | An efficient representation for sparse sets[END_REF]. On the other hand, while the invalidation of row(C, X, b)[↑] = 8 deprives b of a valid support, because b is not contained in dep(C) [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF] we do not need to look for a new valid support for b. In fact, b ∈ dep(C) [START_REF] Bessiere | Using constraint metaknowledge to reduce arc consistency computation[END_REF], so the process of seeking a new valid support for b is activated only when row 2 is removed.

Correctness and Complexity

Correctness is guaranteed by the two following invariants.

Invariant 1 Given C ∈ C, X ∈ scp(C), and a ∈ D c (X), no valid support of (X, a) on C exists in row(C, X, a)[k] for any k such that row(C, X, a).curr < k < row(C, X, a).size.

Proof The invariant holds when the search starts since GACinit has already eliminated all invalid tuples and curr is assigned the maximum value. Considering STR3, we see that curr is decreased only when the row it points to becomes invalid, in which case a new valid support is needed. If a new valid support is found, the previous value of curr is saved so that the search can restart from this point after backtrack. If no valid support is found, curr is left unchanged. Therefore, if a valid support exists in row(C, X, a)[k] for some k such that row(C, X, a).curr < k < row(C, X, a).size, this can only happen as a result of some backtracking from a particular depth. Because any change to curr is recorded in stateR, its previous value must be restored upon backtracking as well, meaning the value of curr must be at least k, contradicting our assumption.

2

Invariant 2 If a ∈ D c (X) and (X, a) ∈ dep(C)[k] then the tuple in row k is a valid support of (X, a) on C.
Proof The invariant holds right after GACinit. We now look at STR3, which is invoked when at least one value becomes absent.

From the code, we see that whenever row k is invalid, any (X, a) in dep(C)[k] will be moved to another dep(C)[j] when a different valid row j is found (Line 27). The invariants for dep(C)[k] and dep(C)[j] are then maintained. If no valid alternative is found, (X, a) becomes invalid. The invariant remains true because we only deal with present values. However, when backtracking occurs we have to re-examine the relationship between row(C, X, a).curr and dep(C).

If (X, a) switches from being absent to present after backtrack, the invariant for it remains true, because either (1) (X, a) is removed as consequence of the instantiation of X to some other value b = a, in which case the invariant is unaffected, or (2) chronological backtracking makes sure that the row (X, a) depended on most recently is restored as well.

An interesting situation happens when (X, a) is present before and after backtrack. In this case, the value of row(C, X, a).curr may be reverted. Assume the value of row(C, X, a)[↑] before the backtrack is k and after backtrack it is j (k < j). This means (X, a) is in 2 STR3 is designed to be incremental by being capable of eliminating repeated domain checks along the same path in the search tree. In our complexity analysis, we consider the worst-case accumulated cost along a single path of length m in the search tree involving a positive r-ary table constraint containing t tuples.

Theorem 1 The worst-case accumulated cost along a single path of length m in the search tree involving a positive r-ary table constraint containing t tuples is O(rt + m) for STR3.

Sketch: STR3's operations can be seen as belonging to two independent phases. First, invalid rows are collected incrementally. indices are duplicated r times in the representation, the collection cost for a single path is O(rt). Second, curr pointers are moved in one direction from one end to the other. In a single path of the search tree, this is equivalent to traversing each element of each tuple in the table once. Hence, the traversal cost is O(rt). Besides, each call to STR3 requires another fixed cost O(1) involving other miscellaneous operations, whose cost can be kept low thanks mainly to the various constant-time sparse set operations. The total cost is O(rt

+ rt + m) = O(rt + m). 2
Since there are rt elements in a table and m << rt in general, it follows that STR3 is path-optimal.

Observation 1 The accumulated cost along a single path of length m in the search tree involving a positive r-ary table constraint containing t tuples can be as much as O(rtm) for STR2.

Reasoning: Recall that STR2 improves over standard STR in two major ways. First, any (X, a) can be disregarded if D(X) is fully supported. Second, no validity check is necessary for (X, a) if it is known that there was no change to the domain of X since the last time STR2 is called. Because STR2 is also sensitive to ordering, we can build a table constraint and a search path such that (1) each call to STR2 involves a domain reduction of exactly one value on every domain, so that the second improvement is useless. (2) each call to STR2 results in exactly one tuple eliminated and this tuple is found at the end of the table. As a result, the cost is O(

m i=1 r(t -i)), which is O(rtm) when m << t. 2
On the other hand, each backtrack costs O(rd) in the worst-case for STR3, where d is the domain size, whereas it is O(r) for STR2.

Similarly, it can be shown that mddc [START_REF] Cheng | An MDD-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints[END_REF] or tries [START_REF] Gent | Data structures for generalised arc consistency for extensional constraints[END_REF] are not pathoptimal. Not counting space for the table representation itself, the space complexity of STR3 for a single table constraint is is O(rd+t) whereas the space complexity for STR2 is O(r).

Experimental Results

In order to show the practical interest of the approach we propose, we have implemented the algorithm STR3 and conducted an experimentation using a cluster of Xeon 3.0GHz with 1GiB of RAM under Linux. Because it has been shown that STR2 is state-of-theart on many series of instances [START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF], we have compared the respective behavior of STR3 and STR2 (we also include results from STR as a baseline). We have considered some series of instances 3 involving positive table constraints with arity greater than 2. We use MAC with the dom/ddeg variable ordering and lexico as value ordering heuristic, to solve all these instances; a time-out of 1, 200 seconds was set per instance. The two chosen heuristics guarantee that we explore the very same search tree regardless of filtering algorithm used.

Table 1 shows mean results (cpu time in seconds and memory usage in MiB) per series. Below the name of each series, we give the number p of instances solved by MAC with the three algorithms (under the form #p) as well as the average proportion q of remaining tuples (i.e., the ratio "size of the current table" to "size of the initial table") over all table constraints and over all nodes of the search tree (under the form avgP = q%). For the crossword series, we have discarded the numerous easy instances, those that require less than 3 seconds to be solved. A first observation is that STR3 requires on average two or three times more memory than STR2; more memory was expected, but this is much better than what worst-case complexity suggests. A second observation is that the results seem to vary widely. STR2 and STR3 are respectively the best approaches on different series: half-7-25 and rand-8-20 for STR2, and rand-5-12 and rand-10-60 for STR3. On series rand-3-20 and crosswords, STR3 offers a more limited benefit. What is interesting to note is that a certain correlation (if we discard renault4 ) between the value of avgP and the ranking of STR2 and STR3 is visible: the higher the value of avgP is, the more competitive STR3 becomes. The ratio cpu STR2 to cpu STR3 is plotted against avgS (average size of tables during search). Dots correspond to instances.

Intuitively, higher value of avgP also implies that there are fewer chances that the solver can reach deeper levels of the search tree. This value seems also to be related to unsatisfiability. To confirm this, Table 2 divides crossword instances according to both satisfiability and an avgP threshold (8%). Table 3 gives details on some representative instances. In particular, on the four crossword instances words and ogd in Table 3, the relationship between the efficiency of STR3 and the value of avgP is obvious. Finally, Figure 8 plots the relative efficiencies of STR2 against STR3 against the value of avgS (the average size of the tables during search), when considering the 658 instances of our experimental study. On the more difficult case where tables remain large (>1000 in Figure 8), STR3 can then be up to 3.6x faster and only 0.6x slower than STR2.

Summary and Future Work

We introduce STR3, a new GAC algorithm for table constraints that is competitive with STR2, a state-of-the art algorithm. Interestingly, STR3 is path-optimal, by being able to completely avoid unnecessary traversal of tables. We have shown that the performance of STR3 correlates to the average number of tuples remaining in the tables during search. The advantage of STR2 appears to depend largely on excessively high rates of table reduction (that is, very low avgP). As soon as the reduction rate drops below 90%, STR2 becomes much less effective. STR2 dominates on benchmarks tailor-made to suit table reduction algorithms -we have seen a surprising number of problems with tables that are virtually wiped out -but apart from these cases STR3 has the upper hand.

The results in this paper exhibit a clear and actionable boundary for which exploitation of both STR2 and STR3 is obvious: a problem instance from a class with low avgP should be solved with STR2; the others should be solved with STR3. To cope with problems of unknown quality, many possibilities exist. For example, one can imagine a simple hybrid algorithm that probes the problem first by tentatively solving it for a while with STR2 while collecting data on table ratios. If it turns out that avgP exceeds a certain threshold, the algorithm may restart and switch to STR3.

Figure 1 :

 1 Figure 1: Two representations of the same table constraint.

1

  GACinit(C: Constraint) 2 remove invalid tuples from rel(C) 3 inv(C) ← ∅ 4 foreach X ∈ scp(C) and a ∈ D c (X) do 5 row(C, X, a).curr ← row(C, X, a).size -1 6 dep(C)[row(C, X, a)[0]] ← {(X, a)} 7 STR3(C: Constraint, X: Variable, a : Value) 8 prevMembers ← inv(C).members 9 for k ← 0 to row(C, X, a).curr do 10 if row(C, X, a)[k] / ∈ inv(C) then 11 add row(C, X, a)[k] to inv(C) 12 if prevMembers = inv(C).members then return true 13 save(C, prevMembers, stateI) 14 foreach i ∈ {prevMembers + 1, . . . , inv(C).members} do 15 k ← inv(C).dense[i] 16 foreach (Y, b) ∈ dep(C)[k] such that b ∈ D c (Y ) do 17 p ← row(C, Y, b).curr 18 while p ≥ 0 and row(C, Y, b)[p] ∈ inv(C) do 19 p ← p -1 20 if p < 0 then 21 removeValue(Y, b) 22 if D c (Y ) = ∅ then return false 23 else 24 if p = row(C, Y, b).curr then 25 save((C, Y, b), row(C, Y, b).curr, stateR) 26 row(C, Y, b).curr ← p 27 move (Y, b) from dep(C)[k] to dep(C)[row(C, Y, b)[p]] 28 return true 29 save(key, newData, store) 30 if (key, oldData) / ∈ top(store) for any oldData then 31 insert (key, newData) to top(store) 32 restoreR() 33 list ← pop(stateR) 34 foreach ((C, X, a), k) ∈ list do row(C, X, a).curr ← k 35 restoreI() 36 list ← pop(stateI) 37 foreach (C, k) ∈ list do inv(C).members ← k 38 removeValue(X: Variable, a: Value) 39 remove a from D c (X) 40 add (X, a) to the propagation queue

Figure 2 :

 2 Figure 2: Algorithm STR3

  The tuple at row(C, X, a)[↑] becomes invalid while the tuple at row k remains valid. Because we seek a new valid support only when k is invalid (Line 16), nothing needs to be done. • The tuple at row(C, X, a)[↑] remains valid while the tuple at row k becomes invalid. Value (X, a) is simply updated (Line 27). There is no need to seek a new valid support, only verification (Line 18) is required. The dependency list and the separators are synchronized at (X, a) as a result. • The tuple at row(C, X, a)[↑] becomes invalid first, then the tuple at row k becomes invalid afterward. The search for a new valid support proceeds as usual. The dependency list and the separators are synchronized at (X, a) if the search succeeds.

Figure 3 :

 3 Figure 3: Transition diagram with respect to (X, a). When dep fails, STR3 is triggered and we end up with synchronized supports. STR3 is never called when curr fails.

  a

Figure 7 :

 7 Figure 7: Status at t = 3

  dep(C)[k] before backtrack. We consider dep(C)[k] and dep(C)[j] after backtrack. Because backtracking never invalidates tuple, the tuple in row k must still be valid after backtrack. Because dep(C) is not maintained, (X, a) remains in dep(C)[k]. Therefore, the invariant for dep(C)[k] is still true, although row(C, X, a)[↑] is no longer k. The invariants involving values in dep(C)[j] are unaffected. Next, consider what happens if the search goes forward when two distinct valid supports exist. That is, (X, a) ∈ dep[k] while row(C, X, a)[↑] = j = k for some value (X, a). If row k becomes invalid, we need to find a new valid support for (X, a). If there exists 0 < i ≤ j such that row i is valid we merely move (X, a) to dep(C)[i]. The invariants for dep(C)[k] and dep(C)[j] hold afterward. If no valid support is found, (X, a) remains in dep[k] and a becomes absent, making the invariant trivially true.

Figure 8 :

 8 Figure 8: The ratio cpu STR2 to cpu STR3 is plotted against avgS (average size of tables during search). Dots correspond to instances.

Table 1 :

 1 Mean cpu time (in seconds) to solve instances from different series (a time-out of 1, 200 seconds was set per instance) with MAC.

			STR	STR2	STR3
	crosswords	cpu	127	72	68
	(#82 -avgP=14.4%)	mem	25M	27M	96M
	renault	cpu	22	16.5	17.0
	(#47 -avgP=30.6%)	mem	31M	32M	65M
	rand-3-20	cpu	126	92	74
	(#49 -avgP=8.4%)	mem	25M	25M	40M
	rand-5-12	cpu	60	38	14
	(#50 -avgP=24.5%)	mem	137M	137M	382M
	half-7-25	cpu	230	125	516
	(#5 -avgP=0.9%)	mem	160M	160M	410M
	rand-8-20	cpu	17	16	235
	(#18 -avgP=0.2%)	mem	116M	116M	244M
	rand-10-60	cpu	465	236	125
	(#29 -avgP=23.0%)	mem	74M	76M	281M

Table 2 :

 2 Mean cpu time (in seconds) to solve Crossword instances (a time-out of 1, 200 seconds was set per instance) with MAC.

			STR	STR2	STR3
	Structured instances
	words-23-02	cpu	409	362	459
	(sat -avgP=2.5%)	mem	120M	123M	314M
	ogd-11-13	cpu	> 1, 200	1,017	784
	(unsat -avgP=11.9%)	mem		32M	177M
	ogd-11-14	cpu	884	412	253
	(unsat -avgP=20.8%)	mem	30M	30M	163M
	ogd-14-14	cpu	42.6	19.8	9.5
	(unsat -avgP=34.1%)	mem	20M	21M	112M
	renault-42	cpu	18.6	14.3	15.4
	(unsat -avgP=33.6%)	mem	34M	34M	60M
	Random instances
	rand-3-20-26	cpu	275	188	130
	(sat -avgP=7.7%)	mem	35M	35M	50M
	rand-5-12-26	cpu	52.8	35.7	11.8
	(unsat -avgP=25.5%)	mem	138M	138M	383M
	half-7-25-2	cpu	210	124	438
	(sat -avgP=1.2%)	mem	156M	156M	404M
	rand-8-20-8	cpu	93.5	54	525
	(sat -avgP=0.2%)	mem	115M	117M	273M
	rand-10-60-5	cpu	258	154	71
	(unsat -avgP=26.5%)	mem	123M	123M	335M

Table 3 :

 3 Detailed results on selected instances.

CRIL-CNRS, UMR 8188, Université d'artois, France, lecoutre@cril.fr

School of Computing, National University of Singapore, Republic of Singapore, likitchav@gmail.com, ryap@comp.nus.edu.sg

Available at http://www.cril.univ-artois.fr/CSC09.

In[START_REF] Lecoutre | STR2: Optimized simple tabular reduction for table constraints[END_REF] renault was also reported as one of the few series of benchmarks where there was little difference between STR2 and mddc.