Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Agnès Maurel
email: agnes.maurel@espci.fr

Jean-Jacques Marigo
email: marigo@lms.polytechnique.fr

Numerical modeling of the acoustic wave propagation across an homogenized rigid microstructure in the time domain

Keywords: homogenization, effective jump conditions, time-domain wave propagation, immersed interface method, ADER scheme

Homogenization of a thin micro-structure yields effective jump conditions that incorporate the geometrical features of the scatterers. These jump conditions apply across a thin but nonzero thickness interface whose interior is disregarded. This paper aims (i) to propose a numerical method able to handle the jump conditions in order to simulate the homogenized problem in the time domain, (ii) to inspect the validity of the homogenized problem when compared to the real one. For this purpose, we adapt the Explicit Simplified Interface Method originally developed for standard jump conditions across a zero-thickness interface. Doing so allows us to handle arbitrary-shaped interfaces on a Cartesian grid with the same efficiency and accuracy of the numerical scheme than those obtained in an homogeneous medium. Numerical experiments are performed to test the properties of the numerical method and to inspect the validity of the homogenization problem.

Introduction

The description of the interaction of waves with many scatterers of size much smaller than the wavelength is in principle simple since the scattering is weak, and several approximated methods can be applied, owing to a small parameter being the scattering strength. In the case of many scatterers located periodically, homogenization techniques are well adapted to handle the problem within a rigorous mathematical framework. The classical homogenization of massive media predicts that scatterers occupying a large area can be replaced by an equivalent homogeneous and in general anisotropic medium occupying the same area [START_REF] Cioranescu | An introduction to homogenization[END_REF]. The case of scatterers located periodically along a mean line Γ (Fig. 1(a)) has been less regarded but it is nowadays accepted that jump conditions across an equivalent interface have to be thought (Fig. 1(b)). Such jump conditions can be established using homogenization techniques which basically rely on the same ingredients than the classical homogenization of massive media, see e.g. [START_REF] Capdeville | Second-order homogenization of the elastic wave equation for non-periodic layered media[END_REF][START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for p-sv waves[END_REF][START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF][START_REF] Delourme | High-order asymptotics for the electromagnetic scattering by thin periodic layers[END_REF]; we call them interface homogenizations. When the numerical resolution of the problem is sought, such equivalent media are of particular interest since they avoid to mesh the very small scales which are essential in the real problem. Indeed, these small scales define the scales of variation of the evanescent near wavefield excited in the vicinity of the scatterers. This is particularly true for scatterers located along a line since the resulting film is dominated by the effect of the near field, or in other words by boundary layer effects. In the case of sound hard scatterers, such interface homogenization has been proposed recently in [START_REF] Marigo | Two scale homogenization to determine effective parameters of thin metallic structured films[END_REF][START_REF] Maurel | Homogenization of ultrathin metallodielectric structures leading to transmission conditions at an equivalent interface[END_REF][START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF]. These works follow those developed in solid mechanics [START_REF] Marigo | The effective behavior of elastic bodies containing microcracks or microholes localized on a surface[END_REF] and they are adapted to the wave equation in the harmonic regime [START_REF] Marigo | Two scale homogenization to determine effective parameters of thin metallic structured films[END_REF][START_REF] Maurel | Homogenization of ultrathin metallodielectric structures leading to transmission conditions at an equivalent interface[END_REF] and in the time domain [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF].

The goal of the present paper is twofold. First, we propose an accurate timedomain numerical scheme for the homogenized problem, thus incorporating the jump conditions across the equivalent thick interface. Next, we inspect the validity of the homogenized result; notably we exemplify the robustness of the jump conditions with respect to the ratio of the array spacing with the typical wavelength and with respect to a possible curvature of the line Γ. Concerning the proposed numerical scheme, specific aspects will be addressed:

• the capability to handle arbitrary-shaped interfaces on a Cartesian grid, without introducing spurious diffractions due to a naive stair-step discretization of the interfaces;

• the accuracy of the numerical scheme, despite the non-smoothness of the solution across the thick interface;

• the performances of the scheme in terms of computational cost; specifically, the scheme has to guaranty an additional cost in the homogenized problem (due to the treatment of the jump conditions) which is highly negligible compared with the computational cost in the real problem.

To do so, cartesian grid methods such as the Immersed Interface Method (IIM) constitute an attractive strategy. The IIM modifies the stencil at the grid nodes near the interfaces. It was originally developed for elliptic and parabolic equations [START_REF] Li | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF][START_REF] Li | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF][START_REF] Li | The Immersed Interface Method for the Navier-Stokes equations with singular forces[END_REF], to handle e.g. flows with irregular domains. Then it has been applied to linear hyperbolic problems in [START_REF] Zhang | The Immersed Interface Method for acoustic wave equations with discontinuous coefficients[END_REF]. An overview may be found in the review book [START_REF] Li | The Immersed Interface Method. Numerical Solution of PDEs Involving Interfaces and Irregular Domains[END_REF].

Alternative approaches have been developed in the context of wave propagation by one author of the present paper. The Explicit Simplified Interface Method (ESIM) does not modify explicitly the stencil, which may be a tedious task for a high-order scheme. On the contrary, the numerical values used by the scheme are modified, based on the jump conditions, the conservation laws and the local geometry. Generalization to different constitutive laws on both sides of the interface can be handled in a naturel way, contrary to the IIM. This approach has been applied to a broad class of problems, such as acoustics / elastodynamics [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF], seismology [START_REF] Lombard | Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves[END_REF], poroelasticity [START_REF] Chiavassa | Time domain numerical modeling of wave propagation in 2D heterogeneous porous media[END_REF]. Modeling of imperfect contacts between solids has also been investigated by this approach [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF].

The extension of the ESIM to the present homogenized problem involves two ingredients. First, it requires to implement a generalized version of the usual boundary conditions at an interface, able to involve both the field and its spatial derivatives. Second, it must handle an interface with a non-zero thickness: the values of the solution on both sides of the interface are linked together (via the jump conditions), whereas no field is computed inside the thin interface.

The paper is organized as follows. In section 2, the actual and the homogenized problems are presented. Section 3 details the numerical methodology: a fourth-order ADER scheme combined with the ESIM, this latter constituting the core of the work. Section 4 presents numerical experiments. Comparisons with exact solutions confirm the efficiency and the accuracy of the numerical modeling. Next, comparisons with direct simulations for a real microstructure confirm the second-order accuracy of the effective model. Lastly, some perspectives are drawn in section 5.

The real problem and its homogenized version

The real problem concerns the propagation of acoustic waves through sound hard scatterers located periodically onto a mean line Γ (Fig. 1(a)). Periodically means a constant spacing h between two scatterers defined by the arc length along Γ. The resulting curved array is surrounded by a fluid with mass density ρ and isentropic compressibility χ. In the fluid, the linearized Euler equations apply and Neumann boundary conditions for the pressure apply for sound-hard scatterers.

Denoting Ω the computational domain containing the array of scatterers, the real problem consists to solve in

Ω                    ρ ∂V ∂t = -∇p, (1a)
χ ∂p ∂t = -∇.V, (1b)
∇p .N = 0 on the boundaries of the scatterers,

with V = (v x , v y ) T the acoustic velocity and p the acoustic pressure (χ is often written in terms of the sound speed c: χ = (ρc 2 ) -1 ). The Neumann boundary condition ∇p .N = 0, with N the vector locally normal to the scatterer boundaries, accounts for a large mass density of the scatterers; it results a vanishing normal velocity at their boundaries.

If the central wavelength (or the minimum one) λ 0 imposed by the wave source is much larger than h, the real problem can be replaced by an equivalent homoge-nized problem, owing to the introduction of the small parameter

ε ≡ 2π h λ 0 . (2)
In [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF], such a homogenization has been proposed; the problem ends up with jump conditions for the pressure and for the normal velocity. In this reference, the derivation is performed for a line Γ being straight (Fig. 2(a)). Here, we heuristically extend this result to a curved line Γ, just by replacing locally the jump conditions expressed in Cartesian coordinates by jump conditions expressed in the local coordinates defined by the normal and tangent vectors to Γ (Fig. 2(b), and we shall discuss the validity of this extension in this paper). Next, in [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF] (see also the previous works of [START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF][START_REF] Delourme | High-order asymptotics for the electromagnetic scattering by thin periodic layers[END_REF]), the homogenized problem is shown to be associated with a satisfactory equation of energy conservation if the jump conditions are expressed across an "enlarged" interface of interior Ω e delimited by the two lines Γ -and Γ + . Doing so allows us to define a positive interface energy, thus guaranties notably the well-posedness in the time domain, which is necessary for our present purpose. The thickness of the interface is a priori arbitrary as soon as it guaranties a positive interface energy. In this paper, following [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF], it is taken equal to the thickness e of the array of Neumann scatterers. We restrict ourselves to scatterers symmetrical w.r.t. the direction perpendicular to Γ.

In this case, with (n, τ ) denoting the local normal and tangent unitary vectors to Γ, the jump conditions read

           p = B ∂p ∂n , v n = C 1 ∂v n ∂n + C 2 ∂v τ ∂τ , (3)
where v n = V.n, v τ = V.τ , and for any function f ,

∂ n f = ∇f.n, ∂ τ f = ∇f.τ .
Also, we define

f = f |Γ + -f |Γ -, and f = 1 2 f |Γ -+ f |Γ + (4)
being the jump of f across the homogenized interface and the mean value of f respectively (the values f |Γ -and f |Γ + are defined locally on Γ -and Γ + along the n direction). In the above expressions, B, C 1 and C 2 are interface parameters which depend on the shape of the scatterers only. As defined in (3), these parameters have the dimension of lengths. Finally, the homogenized problem consists in solving in

Ω\Ω e                  ρ ∂V ∂t = -∇p, (5a)
χ ∂p ∂t = -∇.V, (5b)
Jump conditions [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for p-sv waves[END_REF] across

Ω e . ( 5c
)
Similar non-trivial jump conditions have been obtained by various authors in related problems. See for instance [START_REF] Rohan | Homogenization of the acoustic transmission through a perforated layer[END_REF][START_REF] Dörlemann | Transmission conditions for the Helmholtz-equation in perforated domains[END_REF] for a study of acoustics in perforated domains. In the following section, the numerical resolution of the homogenized problem ( 5) is specifically addressed.

Numerical methods

The numerical scheme used to solved (5) is implemented using the first-order hyperbolic system

∂ ∂t u + A ∂ ∂x u + B ∂ ∂y u = 0, (6)
where ( 6) is deduced from (5) by setting

u =    v x v y p    , A =    0 0 1/ρ 0 0 0 1/χ 0 0    , B =    0 0 0 0 0 1/ρ 0 1/χ 0    . (7)
The ADER-r scheme is used [START_REF] Schwartzkopff | Fast high-order ADER schemes for linear hyperbolic equations[END_REF] to integrate [START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF]. It is an explicit and two-step finite-difference scheme of order r in both space and time; here we use r = 4, which amounts to a fourth-order Lax-Wendroff scheme [START_REF] Lörcher | Lax-Wendroff-type schemes of arbitrary order in several space dimensions[END_REF]. It is dispersive of order 4 and dissipative of order 6 [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF]. Finally, it is stable under the CFL condition β = c ∆t ∆x ≤ 1 (in two dimensions). The solution u is discretized in space on a uniform Cartesian grid with mesh sizes ∆x = ∆y, and in time with a time step ∆t, and we denote u n i,j the discretized numerical value of u at M (i, j) at (x = i∆x, y = j∆y) and at time t = n∆t. The calculation of u n i,j involves a stencil of (r + 1) 2 nodes centered at M (25 nodes for ADER 4), which is written formally as the time-marching

u n+1 i,j = H u n i 0 ,j 0 , {i 0 -i, j 0 -j} ∈ {-2, • • • , +2} , (8)
with H being deduced from ( 6) and [START_REF] Delourme | High-order asymptotics for the electromagnetic scattering by thin periodic layers[END_REF].

The jump conditions (3) of the homogenized problem are discretized by the scheme [START_REF] Dörlemann | Transmission conditions for the Helmholtz-equation in perforated domains[END_REF]. We introduce now a numerical method that incorporates these conditions in the numerical scheme, for a negligible additional cost.

Modified ADER scheme in the neighborhoods of Ω e

To solve the homogenized problem (5), the scheme (8) has to be adapted. Because no constitutive law is defined in the thick interface Ω e , the solution u in this region is not defined. Thus, we shall first distinguish so-called regular and irregular points M depending on the position of the stencil of M with respect to the interface. In the following, we denote Ω ± the subdomains of Ω above Γ + and below Γ -(Fig. 3), and

Ω + ∪ Ω -= Ω\Ω e .
When all the nodes P of the stencil at M fall in the same physical medium (Ω - or Ω + ), the point M is called a regular point. For these points, (8) can be used straightforwardly using u(M ) = u n i,j and u(P ) = u n i 0 ,j 0 being the discretized values of the continuous solution u; we call them direct values.

In the neighborhood of Ω e , it happens that the stencil of M crosses Γ -or Γ + (Fig. 3(a)). It results that nodes Q in the stencil of M belong to Ω e . Such points M are called irregular since u(M ) = u n i,j in (8) requires the values of u at the nodes Q, and u is not defined in Ω e . In the forthcoming discussion, we focus on an irregular node Q close to Γ -: it means that Q is used in the stencils of grid nodes in Ω -.

To account for the regular and irregular points, it is possible to adapt (8) in a simple way. This is done by attributing modified values u * (Q) for the points Q in Ω e , and the direct values otherwise. Specifically, the modified ADER 4 is modified as

             u n+1
i,j = H ũn i 0 ,j 0 , {i 0 -i, j 0 -j} ∈ {-2, . . . , +2} , ũn i 0 ,j 0 = u n i 0 ,j 0 , (x i 0 , y j 0 ) / ∈ Ω e , ũn i 0 ,j 0 = u n * i 0 ,j 0 , (x i 0 , y j 0 ) ∈ Ω e .

In [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF], the numerical values of the solution are affected to all the nodes in the computational domain (with direct values u in Ω -and Ω + and modified values u * in Ω e ). We shall see in the following section that the modified values (at the points Q) are expressed in terms of the direct values at points P ± in Ω ± (Fig. 3(b), the choice of the P ± is incidental at this stage). It follows that ( 9) is solved implicitly on the direct values only: the solution in Ω e is not questioned, as expected (Q does not appear as the center of a stencil in (9)).

Construction of the modified values u

* at the nodes Q in Ω e
As previously said, the solution in Ω e is undefined (u * (Q) does not exist). Only the jump conditions across Ω e make sense. Thus, we shall start by defining u * (Q) and this will be done using Q ± being the projections of Q on Γ ± along n, and using the jump conditions (3) which apply between Q + and Q -. Next, because Q ± are not nodes (except by casuality), the solutions at Q ± will be expressed in terms of two sets of direct values u(P ± ). This is the meaning of the nodes P ± , chosen in the vicinity of Q ± .

The construction of u * (Q) presented in the forthcoming section applies for any sets of P ± being nodes of Ω ± in the vicinities of Q ± . Fig. 3(b) illustrates our choice: P ± are the nodes of Ω ± in the disks D ± centered at Q ± and of radius d. In practice, d 3.5∆x and this is discussed further in Sec. 3.3.

The ingredients of the construction of u * (Q)

In this section, the construction of u * (Q) is presented formally. It will be detailed precisely in the Sections 3.2.2 to 3.2.5. Let us recall that Q is used for time-marching at points of Ω -, and thus is close to Γ -. Because the solution is not defined in Ω e , we start by defining u * (Q) as the smooth extension of the solution in Ω -. This is done using a Taylor expansion for Q in the neighborhood of Q -, written formally

u * (Q) = T(Q, Q -)w(Q -). (10)
In one dimension at order 1,

u * (Q) = u(Q -) + (x Q -x Q -)∂ x u(Q -); w(Q -) encapsulates the weights (u(Q -), ∂ x u(Q -)) while T(Q, Q -)
are the polynomial forms depending on both Q and Q -, here 1 and

(x Q -x Q -).
It is worth noting that [START_REF] Li | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF] introduces a disymmetry between Ω -and Ω + . Indeed, if u * (Q) is a smooth extension of the solution in Ω -, it cannot be a smooth extension of the solution in Ω + because of the jump conditions. These jump conditions apply between Q - and Q + but Q ± do not coincide with nodes. Thus, we use the sets of nodes P ± in the vicinities of Q ± , and the Taylor expansions

u(P ± ) = T(P ± , Q ± )w(Q ± ), (11)
which involves the direct values u(P ± ).

From [START_REF] Li | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF][START_REF] Li | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF], it is visible that u * (Q) can be expressed as a function of u(P ± ) if a relation between w(Q + ) and w(Q -) is established. This relation will obviously involve the jump conditions. Before doing so, the number of terms in w is reduced owing to high-order compatibility conditions, deduced from the initial condition ∇ ∧ v = 0 (coming from (5a) and valid in Ω ± ). The new vector v collects the reduced numbers of unknowns of w, with

w(Q ± ) = Gv(Q ± ), w(P ± ) = Gv(P ± ), (12)
and G is a constant matrix. Finally, the jump conditions are used to get a relation between v(Q + ) and v(Q -), written

v(Q + ) = Func v(Q -) . (13)
The relations [START_REF] Li | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF] to ( 13) allow to conclude. First, u(P ± ) can be expressed as a function of v(Q -) only. From ( 11), we have u(P -) = T(P -, Q -)Gv(Q -). For u(P + ), it starts the same, with u(P + ) = T(P + , Q + )Gv(Q + ) and ( 13) allows to conclude (and ( 13) is essential since it encapsulates the jump conditions). Using these relations for all the nodes P ± ∈ D ± and collecting the direct values u(P ± ) in a single vector

U = u(P -) |P -∈D -, u(P + ) |P + ∈D + , the formal relation U = Mv(Q -) can be inverted to get v(Q -) = M -1 u(P -) u(P + ) . (14)
Finally, ( 14) is injected in [START_REF] Li | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF], using [START_REF] Li | The Immersed Interface Method for the Navier-Stokes equations with singular forces[END_REF], to get the modified values u * (Q) as a function of the direct values u(P ± )

u * (Q) = T(Q, Q -) G M -1 u(P -) u(P + ) . ( 15
)
We shall now detail the steps in the construction of u * (Q) as implemented numerically for k-th order Taylor expansions.

The Taylor expansions, Eqs. (10)-(11)

Here, we simply specify the notations in ( 10) and ( 11) to get Taylor expansions at the order k. The matrix Next, we collect in a single vector w k the n W limit values of u and of its successive spatial derivatives up to the k-th order, at Q ± :

T k of k-th order expansions for Q near Q -reads T k (Q, Q -) = I 3 , • • • , 1 ! ( -m)! (x Q -x Q -) -m (y Q -y Q -) m I 3 , • • • , (y Q -y Q -) k k! I 3 , (16
w k (Q ± ) = lim Q ± ∈Ω ± u T , ..., ∂ ∂x -m ∂y m u T , ..., ∂ k ∂y k u T T , (17)
with = 0, ..., k and m = 0, ..., . The modified value u * (Q) is defined as a smooth extension of the solution in

Ω - u * (Q) = T k (Q, Q -) w k (Q -), (18)
and u * (Q) appears to depend on k. As previously said, once this definition is chosen, u * (Q) cannot be a smooth extension of the solution in Ω + because of the jump conditions.

The story is different for the Taylor extensions of the direct values u(P ± ) written crudely in [START_REF] Li | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF]. Because u(P ± ) are the discretized versions of the exact solution, their Taylor expansions have to be written as approximations, namely

u(P ± ) = T k (P ± , Q ± )w k (Q ± ) + O(∆x k+1 ). (19)
3.2.3. High-order compatibility condition, Eq. ( 12)

The equation (5a) provides a compatibility condition, ∇ ∧ v = 0, which tell us that the fluid is irrotational in Ω ± . Assuming sufficiently smooth solutions in Ω ± , this relation can be differentiated ( -1) times w.r.t. to x and y:

∂ v y ∂x -m-1 ∂y m+1 = ∂ v x ∂x -m ∂y m , ≥ 1, m = 0, • • • , -1, (20)
whose version with = 1, m = 0 is the originate condition. Doing so for = 1, . . . , k provides k(k + 1)/2 high-order compatibility conditions. The equations are valid everywhere outside Ω e , and in particular at Q ± . This allows to use vectors v k containing only the remaining independent derivatives, thus being reduced in size (to

n V = (k + 1)(k + 3) components). The v k are linked to w k by w k (Q ± ) = G k v k (Q ± ), (21)
where G k is a n W × n V matrix (an algorithm to compute G k can be found in [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF], see the appendix A in this reference). Now, [START_REF] Maurel | Homogenization of ultrathin metallodielectric structures leading to transmission conditions at an equivalent interface[END_REF] can be written as a function of v k (Q ± ), using [START_REF] Marigo | An interface model for homogenization of acoustic metafilms[END_REF],

u(P ± ) = T k (P ± , Q ± )G k v k (Q ± ) + O(∆x k+1 ). ( 22
)

High-order order jump conditions, (13)

To express the jump conditions between Q -and Q + , it is useful to come back to the construction of Γ ± and of Q ± (Fig. 4). Because Γ ± are deduced from Γ by the translations of ±e/2 along n, the vectors locally normal and tangent to Γ ± at Q ± are the same, and the same as the vectors (n, τ ), locally normal and tangent to Γ at Q Γ . In practice, we use a parametrization of Γ with s, with [(X(s), Y (s)] the coordinates of Q Γ . We denote n = (Y (s), -X (s)) and τ = (X (s), Y (s)), with prime denoting the derivative. From (3) and ( 4), the fields at Q ± are related by

           p(Q + ) -p(Q -) = B 2 ∂p ∂n (Q + ) + ∂p ∂n (Q -) , v n (Q + ) -v n (Q -) = C 1 2 ∂v n ∂n (Q + ) + ∂v n ∂n (Q -) + C 2 2 ∂v τ ∂τ (Q + ) + ∂v τ ∂τ (Q -) . (23)
The normal and tangential derivatives involved in [START_REF] Pendry | Mimicking surface plasmons with structured surfaces[END_REF] are expressed in terms of the parametric representation of Γ: from the components of (n, τ ):

v n = Y v 1 -X v 2 , v τ = X v 1 + Y v 2 , ∂v n ∂n = Y 2 ∂v 1 ∂x -X Y ∂v 2 ∂x + ∂v 1 ∂y + X 2 ∂v 2 ∂y , ∂v τ ∂τ = X 2 ∂v 1 ∂x + X Y ∂v 2 ∂x + ∂v 1 ∂y + Y 2 ∂v 2 ∂y , ∂p ∂n = Y ∂p ∂x -X ∂p ∂y . (24)
Injecting ( 24) into ( 23) leads to lengthy expressions, which are summed up by the matrix form

E 0 w 0 (Q + ) -w 0 (Q -) = 1 2 F 1 w 1 (Q + ) + w 1 (Q -) , (25)
where the matrices E 0 and F 1 depend on the geometry of the interface because of their dependance on (n, τ ). E 0 is a 2 × 3 matrix defined by

E 0 = Y -X 0 0 0 1 , (26)
and F 1 is a 2 × 9 matrix whose non-zero components are

F 1 (1, 4) = C 1 Y 2 + C 2 X 2 , F 1 (1, 5) = F 1 (1, 7) = -C 1 X Y + C 2 X Y , F 1 (1, 8) = C 1 X 2 + C 2 Y 2 , F 1 (2, 6) = BY , F 1 (2, 9) = -BX . (27)
Now, the goal is to obtain an extended version of the jump conditions involving derivatives of u up the k-th order, namely involving w k . For this purpose, ( 25) is differentiated w.r.t the time t and the parameter s, and using

         ∂ ∂t w = -A ∂ ∂x w --B ∂ ∂y w , ∂ ∂s = X (s) ∂ ∂x + Y (s) ∂ ∂y , (28)
where the first relation holds for any . Applying to [START_REF] Rohan | Homogenization of the acoustic transmission through a perforated layer[END_REF] the chain-rule

∂ ( -i) ∂t ( -i) ∂ i ∂s i E 0 w 0 (Q + ) -w 0 (Q -) - 1 2 F 1 w 1 (Q + ) + w 1 (Q -) = 0, (29)
where i = 0, . . . , , provides ( + 1) matrix relations, or 2( + 1) scalar relations, at each order . These relations involve spatial derivatives up to the order ( + 1) and the chain-rule is stopped at the order k. The n J = (k + 1)(k + 2) obtained relations for = 0 to k are collected in the matrix relation

E k w k (Q + ) -w k (Q -) = 1 2 F k+1 w k+1 (Q + ) + w k+1 (Q -) . ( 30
)
For instance, at k = 1 one gets

E 1 =      E 0 0 E 0 0 -E 0 A -E 0 B E 0 X E 0 Y E 0      . (31)
Computing the matrices E k and F k+1 for k > 1 is a tedious task, that has been automatized thanks to computer algebra tools. It is worth noting that they depend on the first k-th derivatives of X(s) and Y (s), thus they encapsulate informations on the interface shape well beyond its position only.

In (30), w k+1 involves spatial derivatives up to the order (k +1), and this order is not considered in the resolution at the order k. Considering F k+1 the restriction of F k+1 obtained by removing the last (n W + 1)-th column in the matrix F k+1 , (30) simplifies to

E + k w k (Q + ) = E - k w k (Q -), (32)
where 21) and (32) yields

E ± k ≡ E k ∓ 1 2 F k+1 are two n J × n W matrices. Combining (
S + k v k (Q + ) = S - k v k (Q -), (33)
with the n J × n V matrices S ± k = E ± k G k being smaller than the n J × n W E ± k . With n J = (k + 1)(k + 2)
and n V = (k + 1)(k + 3), the system (33) is underdetermined. It is inverted in the sense of the least-squares sense using Singular Value Decomposition, leading to

v k (Q + ) = S + k -1 S - k | K S + k   v k (Q -) Λ k   , (34)
where (S + k ) -1 is the least-square pseudo-inverse of S + k , K S + k is the matrix filled with the kernel of S + k , and Λ k is a set of the (n V -n J ) Lagrange multipliers which are the coordinates of v k (Q + ) onto the kernel. A singular value decomposition of S + k is used to build (S + k ) -1 and the kernel K S + k [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF].

3.2.5. Final step in the construction of modified values, ( 14)-( 15)

We want to express v k (Q -) as a function of the set of direct values u(P ± ) (we denote n ± the numbers of nodes P ± ). For the set of nodes P -in Ω -, ( 19) is simply re-written

u(P -) = T k (P -, Q -)G k (1 | 0)   v k (Q -) Λ k   + O(∆x k+1 ). (35)
In the above relation, 1 stands for the n V × n V identity matrix and 0 for the n V × (n V -n J ) zero matrix. For the set of nodes P + in Ω + , we use (34) in [START_REF] Marigo | The effective behavior of elastic bodies containing microcracks or microholes localized on a surface[END_REF], whence

u(P + ) = T k (P + , Q + )G k S + k -1 S - k | K S + k   v k (Q -) Λ k   + O(∆x k+1 ).
(36) Now, we collect in a single vector U the (n -+ n + ) vectors u(P ± )

U =              u(P - 1 ) . . . u(P - n -) u(P + 1 ) . . . u(P + n + )              , (37)
and use ( 35) and (36) written in the matrix form

U = M   v k (Q -) Λ k   + O(∆x k+1 ), (38)
where

M ≡                 T k (P - 1 , Q -)G k (1 | 0) . . . T k (P - n -, Q -)G k (1 | 0) T k (P + 1 , Q + )G k S + k -1 S - k | K S + k . . . T k (P + n+ , Q + )G k S + k -1 S - k | K S + k                 (39) is a (3(n -+ n + ) × (2 n V -n J ) matrix.
To ensure that the system (38) is overdetermined, the radius R D of the discs D ± has to be chosen in order to ensure that

γ(R D , k) ≡ 3 (n -+ n + ) (k + 1)(k + 4) ≥ 1, (40)
and this condition will be further discussed in Section 3.3. The condition (40) being ensured, (38) is solved using the least-squares inverse M -1 of M. The Lagrange multipliers Λ k have been incorporated in the construction of M, but they are not needed to build the modified value. Thus, they are removed using the

(n W -n V ) × 3 (n -+ n + ) restriction M -1 of M -1 , to get only v r (Q -) = M -1 U,
and finally using (18)

u * (Q) = T k (Q, Q -) G k M -1 U. (41)
The matrices

T k (P ± , Q ± ), G k and S + k -1 S - k | K S + k are involved in (41) through M -1 .

Comments and practical details

1. The described algorithm is applied to the irregular points in Ω ± . The sizes of the matrices involved are summarized in table 1. Since the jump conditions do not vary with time, the evaluation of the matrices in (41) is done during a preprocessing step. Only small matrix-vector products are therefore required at each time step. After optimization of the computer codes, this additional cost is made negligible, lower than 1% of the time-marching of the ADER 4 scheme.

n W = 3 (k + 1) (k + 2)/2 n V = (k + 1) (k + 3) < n W n J = (k + 1) (k + 2) < n V Quantity Size T k 3 × n W G k n W × n V S ± k n J × n V M 3 (n + + n -) × (2 n W -2 n V -n J ) M -1 3 (n + + n -) × (n W -n V )
Table 1: Quantities involved in the computation of the modified values (section 3.2).

2. The matrix M in (38) depends on the subcell positions of Q ± inside the mesh and on the jump conditions expressed between Q + and Q -, involving the local geometry and the curvature of Γ ± at Q ± . Consequently, all these insights are incorporated in the modified value (41), and thus also in the scheme. 3. The optimal order k depends on the order r of the scheme and on the jump conditions. Let us begin with the classical case of acoustics where the jump conditions do not involve spatial derivatives: for instance, v n = 0 and p = 0. In this case, taking k = r maintains a r-th order global accuracy [START_REF] Chiavassa | Time domain numerical modeling of wave propagation in 2D heterogeneous porous media[END_REF] (the criterion k = r -1 is even sufficient [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF]). In the non-classical case studied here, the jump conditions involve first-order spatial derivatives. After successive derivations, the higher-order terms are canceled in (32), which introduces a loss of accuracy. To maintain the r-th order convergence, the order of the ESIM must be increased of one unity: k = r + 1 ≡ 5 for the ADER 4 scheme.

Summary of the algorithm

The numerical strategy presented all along this section couples two numerical methods: a finite difference ADER scheme and the Explicit Simplified Interface Method. To clarify the interactions between these methods, the algorithm can be summarized as follows: Preprocessing

• Detection of irregular grid points

• At each irregular point M :

-detection of the orthogonal projections Q ± -computation of the matrices S ± k (33) -filling of M (39) -computation of extrapolation matrix in (41)

• Initialization of the computation at t = 0 Time iterations

• Computation of the modified values (41)

• Time-marching at the regular points (8)

• Time-marching at the irregular points (9)

Numerical experiments

In this section, we first validate the ESIM implemented on the thick interface as presented in the preceding section. This is done for a plane wave at normal and oblique incidences on the interface. These cases allow for analytical solutions of the homogenized problem, already validated in [START_REF] Marigo | Two scale homogenization to determine effective parameters of thin metallic structured films[END_REF][START_REF] Marigo | An interface model for homogenization of acoustic metafilms[END_REF] in the harmonic regime. The extension to the time domain being done by means of discrete inverse Fourier transforms, we do not repeat the validation by comparison of the solutions of the homogenized and of the direct problems. We inspect the accuracy of the numerical scheme varying the order k in the scheme, and discuss the convergence obtained with the estimated optimal value k = 5.

Next, we exemplify the capability of the homogenized problem (5) to mimic the real problem (1) in the time domain. The solutions of the direct and the homogenized problems are compared in the case of a source emitting a short pulse with a central frequency that we vary to inspect the robustness of the homogenized solution w.r.t. the small parameter ε, defined in (2) (and λ 0 will be defined latter). This is done for a straight and curved line array.

The following characteristics hold for all our simulations:

• The real problem is an array of 50 rectangular sound-hard scatterers in water. The periodicity of the array is h = 20 m with a filling ratio of the scatterers ϕ = 0.5; the thickness of the array is e = 20 m. For these dimensions of the scatterers, the interface coefficients entering in the jump conditions of the homogenized problem are

B = 44.412, C 1 = 10, C 2 = 8.338. (42)
(see Appendix A). For water, we use ρ = 1000 kg/m 3 and c = 1500 m/s.

• The numerical results have been performed in a domain of extension 1200 m× 1200 m, discretized on N x × N y points, with N x = N y . The time step follows from the CFL condition: ∆t = β∆x/c, and we used β = 0.95.

• Discrepancies between two solutions yielding the pressure fields p 1 and p 2 are given by the relative difference ||p 1 -p 2 ||/||p 1 || (and ||.|| refers to the discrete L 2 -norm)

• The temporal signals are built using h(t) being a combination of truncated sinusoids

h(t) =      4 m=1 a m sin(β m ω 0 t) if 0 < t < 1 f 0 , 0 otherwise, (43)
where β m = 2 m-1 , ω 0 = 2π f 0 ; the coefficients a m are: a 1 = 1, a 2 = -21/32, a 3 = 63/768, a 4 = -1/512, ensuring C 6 smoothness of h(t). The Fourier transform of (43) reads

ĥ(ω) = ω 0 2π 4 m=1 a m β m e 2iπω 0 /ω -1 ω 2 -β m ω 2 0 , (44)
with a maximum slightly greater than f 0 and a cut-off frequency f m at around 5f 0 . h(t) and ĥ(ω) are shown in Figs. 5 for f 0 = 10 Hz.

In the simulations, various central frequencies are considered: f 0 = 2.5, 5 and 10 Hz. The resulting values of ε calculated using k 0 = 2πf 0 /c are ε = 0.21, 0.41 and 0.83. However the spectral content of the source being large, the wave packet contains wavenumbers 5 times smaller than k 0 ; these small scales are associated to ε larger than 1 (up to 4).

Validation of the ESIM for large interface

To begin with, we examine how accurately the ESIM discretizes the homogenized jump conditions. To that aim, we consider a plane wave at oblique incidence θ I on the interface. This case allows for an exact, or say reference, solution of the homogenized problem in the frequency domain (see Appendix B) u(x, y, ω) = u I (x, y, ω) + R(ω)u R (x, y, ω) + T (ω)u T (x, y, ω),

with u I,R,T given by (B.2) and (R, T ) by (B.5). Afterwards the solution u(x, y, t) in the time domain (called reference homogenized solution in the following) is deduced by discrete inverse Fourier transform of u(x, y, ω).

Numerically, this solution has to be recovered by imposing, at each time step, the reference solution u(i∆x, j∆x, n∆t) on the 2 lines (i = 0, 1 and (N x -1), N x ) and the 2 columns (j = 0, 1 and (N y -1), N y ) at the edges of the computational domain. Once these boundary conditions have been imposed, the numerical scheme has to be able to produce the solution in the whole domain. This numerical solution is referred as the numerical homogenized solution in the following.

In the simulations, we used a temporal signal with Fourier dependence given by (44) at the central frequency f 0 = 10 Hz, resulting in a wave packet with the central and smallest wavelengths of 150 m and 30 m respectively. The homogenized interface mimicking the array of Neumann rectangles (with spacing 20 m, filling fraction 0.5 and thickness 20 m) is thus 20 m large, and associated to the interface parameters (42). The Fourier synthesis is done using N f = 256 modes with a uniform frequency step ∆f = 0.4 Hz around f 0 and yields the reference solution u imposed at the 2 points boundaries of the computational domain, as previously described. The computational domain is 1200 m× 1200 m large and it is discretized using ∆x = ∆y = 2 m (N x = N y = 600), thus from the CFL condition ∆t = 0.0013 s.

Plane wave at normal incidence

We first consider a normal incidence, for which u is independent of x. Fig. 6 shows the pressure field of the reference solution which is imposed in the whole domain at the initial time t = 0; this initial time has been chosen before the wave hits the interface (the profile of this one-dimensional field is also reported). The reference and the numerical homogenized solutions are then computed in time; the pressure fields of the numerical solution after 158 iterations (t = 0.2 s) is reported in Fig. 7(a). Fig. 7(b) reports the profiles along y of the two solutions. The discrepancy between both solutions is 0.2%, a very low error which can be attributable to the discretization. We conclude that the ESIM discretized correctly the jump conditions.

To further inspect the accuracy of the method, we report results varying the k-values and the spatial discretization. With the same ∆x = 2 m, we used the orders k = 1 to 5 to compute the numerical homogenized solution. Fig. 8(a) reports the resulting pressure profiles along y at t = 0.2 s (the profile for k = 5 is indiscernible of the one for k = 4 and it is not reported). Fig. 8(b) reports the errors between these profiles and their reference counterpart as a function of k. Next, we used k = 5 and various discretizations for N x = 300 up to N x = 4800. The error is reported in Fig. 9 as a function of N x . The order measured (being the slope of the curve) is close to 4 at this order as soon as N x = 300 (∆x = 4 m) and this corresponds to the best convergence that we can expect since the order 4 is the order of ADER 4 in homogeneous medium. It is worth noting that this order 4 is not found for k < 5, and such high order has a numerical cost. However, the difference in accuracy between k = 4 and k = 5 being very small, a good compromise in practice is to choose k = 4. In terms of the spatial resolution, we have said that the smallest wavelength is about 30 m; with an error less than 1% for N x = 600 (∆x = 2 m), we can estimate that 15 grid nodes per wavelength calculated for the smallest wavelength is a good criterion to fix the spatial resolution (and this conclusion holds for k = 4 as well).

Plane wave at oblique incidence and tilted interface

The same comparison between the reference and computed homogenized solutions is performed in the case of a wave at oblique incidence and, more importantly, using a tilted interface. As in the previous section, the reference homogenized solution is obtained by discrete inverse Fourier transform of (45) with θ I the angle between the incident wave and the tilted interface (see Appendix B). The case of a tilted interface allows us to inspect another aspect of the efficiency of the ESIM, namely its capacity to account for the interface shape with a subcell resolution on a Cartesian grid. In the present case, this means that the real slope of the interface is accounted for, instead of a crude stair-step discretization.

We considered the interface with a tilt angle of about 10 • with Ox and the incident wave packet makes an angle -30 • with Ox. Figs. 10 show the pressure fields of the numerical homogenized solution at the initial time (identical to the reference one) and after the wave packet has propagated (t = 0.2 s corresponding to 158 iterations). As expected, no spurious diffractions have appeared. More quantitatively, the pressure profile along the vertical slice (x = 740 m and y ∈ [100, 800] m) is reported in Fig. 11 together with the corresponding reference solution (the scattered wave packets (R and T) are visible on these profiles). The discrepancy between the two profiles is of 0.5% , as small as in the case of a non tilted interface. This accuracy could not be obtained with a piecewise constant approximation of the interface shape, and it confirms that the slope of interface is accurately accounted for in the numerical scheme.

Validation of the homogenized problem in the time domain

In the previous section, we have inspected the ability of the ESIM to properly account for the homogenized jump conditions (3). Another question is whether or not the homogenized problem ( 5) is a good approximation of the real one (1); this question is addressed now through two examples. First, we consider an array of scatterers along a straight line Γ, which corresponds to the configuration for which the jump conditions have been derived [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF]. Next, we consider an array located and numerical (open symbols) homogenized solutions; the agreement between both is 0.5% . onto a curved line Γ, for which we extended heuristically the jump conditions. This allows us to inspect the intuitive idea that such extension is possible for small curvature of Γ.

The real problem is solved numerically following the method presented in [START_REF] Lombard | Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves[END_REF]. The numerical method uses a ADER 4 scheme, and is able to accurately account for free boundaries at an interface by affecting so-called fictitious values of the solution inside the sound-hard scatterers. In the spirit, these fictitious values are the equivalent of the modified values presented in this paper, and used in the homogenized thick interface.

The case of a straight array of sound-hard scatterers

The array of rectangular scatterers is placed at y = 500 m along the x-axis. Each rectangle is 10 m large along x, with spacing h =20 m, and e =20 m thick along y. In the homogenized problem, the interface is e =20 m thick (the region y = [490, 510] m is not resolved), and associated to the interface parameters (42). We considered a source point at (x s = 600, y s = 620) m emitting the short pulse h(t) in (43), with central frequency f 0 . When not specified, we used a fine grid ∆x = 0.125 m for the real problem and a coarse grid ∆x = 2 m for the homogenized problem (in both cases, k = 4 has been considered). These are the mesh sizes required to get converged solutions (from ∆x = 2 m, reducing the mesh size to ∆x = 0.5 m produces about 20% variations of the solution of the real problem, while the solutions of the homogenized have already converged, with variations less than 0.1%).

The pressure fields p r computed in the real problem and p in the homogenized one are reported in Figs. 12 for central frequencies f 0 = 2.5, 5 and 10 Hz, at t = 0.38 s (corresponding profiles along the centerline x = 600 m are reported in Figs. [START_REF] Li | The Immersed Interface Method. Numerical Solution of PDEs Involving Interfaces and Irregular Domains[END_REF]). The discrepancy between the two fields is of about 5% for f 0 = 2.5 and 5 Hz, and it is of 10% for f 0 = 10 Hz, and these orders of magnitude are in agreement with those reported in the frequency regime, see [START_REF] Marigo | Two scale homogenization to determine effective parameters of thin metallic structured films[END_REF] (we measured the discrepancy between both fields outside the thick interface y ∈ [490, 510] m where p is not defined). Note that, at f 0 = 10 Hz, the central and smallest wavelengths are 150 and 30 m, leading to ε ∼ 0.8 -4, so overcoming the intuitive limit ε = 1 for the validity of the homogenization. We report in Fig. 13(d) the relative errors for increasing ε-value. Here, the error has been calculated on the part of the profiles corresponding to the transmitted wave packet y ∈ [0, 490] m (thus avoiding to cross the interface) and with ε = 2πf 0 h/c (an extra point at f 0 = 20 Hz has been added). The observed ε 2 scaling is consistant with the expected second-order accuracy of the interface homogenization model.

The computations have been performed using a fine grid ∆x = 0.125 m for the real problem and a coarse grid ∆x = 2 m for the homogenized problem; as previously said, we checked that these meshes are required to get converged solutions. This is not obvious at f 0 = 10 Hz where the source generates wavelengths of the same order of magnitude than the size of the array; thus, we could expect that the two problems, real and homogenized, require the same mesh size, but we observe that it is not the case. Heuristically, this can be explained as follows. The fine grid needed for the real problem resolves the smallest scale, and it turns out that this smallest scale is associated to the evanescent field, excited in the vicinity of the array (this small scale is visible in Fig. 12(a) for f 0 = 10 Hz). This means that the usual rule of say 15 nodes per wavelength has to apply to this near field scale and not only to the incident wavelength; however, the near field scale is not known a priori and it depends on how deeply the evanescent field is excited, so it is difficult to anticipate how fine has to be the grid (a discussion of this point is presented in [START_REF] Marigo | An interface model for homogenization of acoustic metafilms[END_REF]).

In the homogenized problem, the near field effect is encapsulated in the interface parameters and this is possible since the near field is essentially a static field. Thus, the mesh size is limited by the usual rule on the incident field only. This is confirmed here; with a minimum incident wavelength of 30 m, a mesh size ∆x = 2 m is sufficient to get a converged homogenized solution; next the validity of the field depends on how close the real near field is close to the static one, and this is lost progressively by increasing the frequency. Now, let us inspect how different is the story for the mesh size needed in the real problem. Fig. 14 shows a zoom of the profiles along the centerline x = 0 (y ∈ [50, 250] m) for decreasing mesh size ∆x. It is visible that the solution of the real problem continues to converge toward the solution of the homogenized problem (up to the error due to model) for ∆x well beyond ∆x = 2 m, and it has not converged before ∆x = 0.25 m. This means that the evanescent field experiences rapid variations with a typical exponential decrease over lengths of few meters (∆x = 0.25 m solves, as a rule of the thumb, variations over 4 meters of the evanescent field, to be compared to the 30 meters wavelength). . In conclusion, the gain in replacing the real problem by the homogenized one at low frequency is simply given by the ratio between the geometrical size of the array (h, e) on the wavelength. Besides, the smaller is the frequency, the better is the agreement between the homogenized and the real problem. It is obviously for these small ε values that the homogenized problem can mimic the real problem in the most efficient way.

The situation is more involved for intermediate frequencies (say ε of order unity). In this regime, the homogenized solution becomes less efficient to describe the real problem. Nevertheless, the gain in the numerical cost remains important; this is because the smallest scale to be resolved in the real problem is associated to the near field variation, and not anymore to the array size. Because higher frequency produces stronger scattering, the evanescent field may contain scales significantly smaller than the array size. For these higher frequencies, a compromise between the accuracy of the homogenized solution and the numerical gain has to be found, and this depends on the wanted precision. As an indication in the numerical gain at f 0 = 10 Hz, a computational time of 1 minute for the homogenized problem corresponds to a computational time of 10 hours in the real problem. The extra time needed in the real problem is not only due to the smaller mesh grid, but also to the smaller time step imposed by the CFL condition.

Variable homogenized interface

As previously said, the numerical implementation of the jump conditions (3) have been extended to the case of a curved line Γ. Below, we report results varying the curvatures. For small curvatures, this allows us to validate the numerical implementation of the jump conditions along a curved line; also, we inspect the error due to the model when increasing the curvature. For high curvatures, we expect the jump conditions to be modified. It is outside the scope of the present paper to derive such conditions, but let us estimate the maximum curvature below which we expect the jump conditions (3) to be unaffected. The parameters (B, C 1 , C 2 ) have been calculated in static problems to account for the boundary layers effects near the scatterers; if these boundary layers are significantly modified because of the local curvature, the parameters will be affected as well; it is worth noting that, if the case, the homogenized problem becomes more tricky since the parameters will vary along Γ (if the curvature varies, that is for any curve Γ except a circle). Fig. 15 shows the relative position of two -rectangular-scatterers for a local curvature 1/R, producing a minimum distance of δh between them (and δh = h for infinite R, or zero curvature). The modification of the boundary layers in the static problems involved to calculate (B, C 1 , C 2 ) are expected to depend of how δh/h is far from 1. Note also the condition δh > ϕh for no overlapping of the scatterers. It is easy to see that this leads to

1 ≥ δh h = R -e/2 R > ϕ. (46)
Now, in addition to the condition of no-overlapping which imposes e/R < 2(1ϕ), small values of e/R-values (or small 1 -δh/h) are expected to ensure the validity of (3), and this is what we shall inspect further.

To easily increase the curvature, we consider a curved line Γ in the form of a sinusoid with mean value at y m = 500 m and y Γ = y m + A cos 2π(x -x s )/D (x s = 600 m). We kept a fixed A = 10 m value and varied D from D = 1000 m and 125 m. The corresponding maximum local curvatures are given by κ = A(2π/D) 2 , leading to a minimum local radius of curvature R from R = 2500 m to 40 m. This smallest value of R corresponds to twice the minimum value imposed by non overlapping (46).

To begin with, we report in Figs. 16 and17(a-c) the wavefields and the profiles on the centerline x = 600 m in the real problem and in the homogenized problem for f 0 = 10 Hz. The reported time is t = 0.38 s, and the calculations have been performed in the same conditions as in Figs. 12- 12-13, it is visible that even a small curvature produces significant modification in the field pattern. This is particularly noticeable regarding the signal often referred as the "coda", which corresponds to the signal between the two main wavefronts directly transmitted and reflected by the interface. In the real problem, increasing the curvature of Γ enhances multiple scattering effect in the region of the scatterers, which feeds the coda region. Also noticeable is the fact that the homogenized solutions reasonably reproduce the main features of this coda region, although they are unable to reproduce its finest scales.

Finally, Fig. 17(d) shows the error as a function of (1-δh/h) where additional D values have been considered. The result is as expected: for δh h, the effect of the curvature is incidental and the error remains the same as for a straight line Γ; increasing further the curvature produces an increase in the error. In the presented case, this happens for (1-δh/h) above 10%, thus for a relative minimum distance between the scatterers 5 times smaller than the minimum distance of overlapping.

Conclusion

We have proposed a numerical method to implement jump conditions established within a homogenization approach. These jump conditions allow us to replace the real problem of the wave propagation though an array of sound hard scatterers by an equivalent homogenized problem, much simpler to handle numerically. The numerical method has been validated considering a scattering problem for which an explicit solution is available, and comparisons of the solutions in the homogenized problem and in the real problem have confirmed the ability of the homogenized problem to advantageously replace the real one. Advantageously is meant with a computational time highly reduced; in the cases reported here, the computational time is typically reduced by a factor 500. In fact, and as previously said, accounting for the jump conditions requires additional calculations being time independent. Thus, they are done once and for all independently of the wave source and before the time resolution is considered. Then, the time resolution is not more demanding than for a problem of wave propagation in a homogeneous medium.

In terms of the possible extensions of the numerical method, several directions seem of interest to us. The first is rather incremental. We have considered thicknesses of the scatterers being small but sufficiently large with respect to the expected mesh size; this is because we assumed that points of certain stencils may fall within the homogenized interface (which is not resolved). If the scatterer thickness passes below the mesh size, the region of the interface will always be contained between two nodes; this case has been considered already, notably for vanishing scatterer thickness [1] (the harmonic regime is considered in this reference) and it is strictly more simple to handle numerically. Nevertheless, because of the practical interest in such thin arrays (the gain in the numerical implementation is even more significant), this has to be done. The second extension is numerically not so demanding but it requires to adapt the homogenization approach. The jump conditions have been established for the scatterers being located along a straight line and we have inspected their possible extension to the case of a curved line. This has confirmed the intuitive idea that large local curvatures require a modified version of the jump conditions. It is worth noting that this would lead to a more involved homogenized problem since the curvature being defined locally along the mean line, the parameters entering in the jump conditions would be local as well. Nevertheless, and again in regard with the practical interest of such configurations, a generalization to scatterers located onto curved lines deserves interest.

Finally, several extensions concern the nature of the scatterers, and let us men-tion two situations which are not trivial extensions of the present work. The first concerns scatterers associated to Dirichlet boundary conditions: it is typically metallic arrays illuminated by a polarized electromagnetic wave in the far infrared regime. In this case, effective boundary conditions at each side of the interface have to be considered, rather than jump conditions [START_REF] Marigo | Two scale homogenization to determine effective parameters of thin metallic structured films[END_REF]. The second is a bolder extension of the present work. It concerns scatterers with material properties having high contrasts with respect to the surrounding medium, such that resonances inside the scatterers are possible. In this case, the jump conditions cannot be derived in the time domain. The calculations are done in the harmonic regime, revealing interface parameters being frequency dependent. Thus, the numerical implementation of these jump conditions in the time domain requires to handle these frequency dependent parameters. in [START_REF] Delourme | High-order asymptotics for the electromagnetic scattering by thin periodic layers[END_REF]. Below, we shall calculate the pressure p afterwards u = (v x , v y , p) T is deduced using [START_REF] Cioranescu | An introduction to homogenization[END_REF] written in the harmonic regime, with time dependence e -iωt , whence For this one dimensional problem, the pressure field (7) reads p(x, y, ω) = e iωx sin θ/c × e -iω(y-e/2) cos + Re iω(y-e/2) cos θ/c , x > e/2, T e -iω(y+e/2) cos θ/c , x < e/2, (B.2) and we want to determine the scattering coefficients (R, T ). It is sufficient to inject (B.2) in the jump conditions (3) using ( 5

Figure 1 :

 1 Figure 1: (a) The real problem involves a structured film composed of a periodic set of soundhard scatterers in a homogeneous fluid, (b) the equivalent homogenized problem involves a thick interface across which jump conditions apply. The wavefield, on which the jump conditions apply, is not represented.

Figure 2 :

 2 Figure 2: In the homogenized problem, the array of sound-hard scatterers located periodically along a mean line Γ is replaced by a domain Ω e . Ω e is delimited by the lines Γ -and Γ + across which the jump conditions (3) apply. (a) in [20], the jump conditions are expressed for a straight line Γ. (b) we extend these conditions to a curved line Γ; Γ -and Γ + are deduced from Γ by a translation of ±e/2 along the normal vector n to Γ; the normal and tangent vectors to Γ ± at the resulting points are the same than at the initial point on Γ.

Figure 3 :

 3 Figure 3: (a) Example of irregular point M ; in the stencil of M (grey region), Q (•) are in Ω e (the usual nodes, in Ω ± are denoted P (×)). (b) Construction of the nodes P ± (×) in Ω ± used to calculated the modified value u * (Q). With Q ± the projections of Q on Γ ± along the normal n, the P ± are the nodes in the half-disks D ± of centers Q ± (P ± belong to Ω -and Ω + ).

  ) with = 0, • • • , k and m = 0, • • • , and I 3 the 3 × 3 identity matrix. With ( + 1) polynomial forms at the order , T k is a 3×n W matrix, with n W = 3(k+1)(k+2)/2.

Figure 4 :

 4 Figure 4: Parametrization of Γ with the parameter s; the vectors (n, τ ) locally normal and tangent to Γ at Q Γ [X(s), Y (s)] are also normal and tangent to Γ ± at Q ± by construction.

4 .

 4 The simulations indicate that overestimation of γ in (40) ensures the stability of the ESIM. Numerical experiments have shown that d = 3.5 ∆x is a good choice. Typically, it gives n ± ≈ 10 and γ ≈ 4.

Figure 5 :

 5 Figure 5: (a) Temporal signal at a source point (43) for f 0 = 10 Hz, and (b) corresponding spectrum.

Figure 6 :

 6 Figure 6: (a) Pressure field of the incident wave (for the reference solution) imposed at the initial time of the simulation and (b) corresponding yprofile.

Figure 7 :

 7 Figure 7: (a) Pressure field of the numerical homogenized solution at t = 0.2 s (see text) and (b) Pressure profiles of the reference (black line) and numerical (open symbols) homogenized solutions; the relative error between both profiles is of 0.2%

Figure 8 :

 8 Figure 8: Influence of the order k in the computation of the numerical homogenized solution. (a) Zoom on the pressure profile of Fig. 7(b); the profiles for k = 1, 2, 3 and 4 are shown (the profiles for k = 4 and 5 are indiscernable). (b) Errors of the numerical solution compared to the reference solution (calculated for y ∈ [0, 1200]) as a function of the order k.

Figure 9 :

 9 Figure 9: Error as a function of the mesh size N x (for k = 5). The numerical values indicate the slopes between two successive points, the slope -4 is shown in dotted line.

Figure 10 :

 10 Figure 10: (a) Pressure field of the numerical homogenized solution (identical to the reference one) imposed at the initial time of the simulation, (b) Pressure field of the numerical homogenized solution at t = 0.2 s; the vertical slice (dotted lines) at x = 740 m, y ∈ [100, 800] m is used in Fig. 11.

Figure 11 :

 11 Figure 11: Pressure profiles along the slice in dotted line in Fig. 10(b), of the reference (black line) and numerical (open symbols) homogenized solutions; the agreement between both is 0.5% .

Figure 12 :

 12 Figure 12: Pressure fields p r (real problem with the array) and p (homogenized problem) computed with ADER 4 at t = 0.38 s for f 0 = 2.5, 5 and 10 Hz. f 0 is the central frequency of the signal h(t) in (43) imposed by the source point (indicated by the black point).

Figure 13 :

 13 Figure 13: (a-c) Pressure profiles along the centerline x = 600 m from Figs. 12, p r (blue symbols) and p (black lines). (d) Error |p r -p|/|p| calculated on the transmitted wave y ∈ [0490] m as a function of ε = 2πf 0 /c (an extra point at f 0 = 20 Hz has been considered); dotted grey line shows the ε 2 law.

Figure 14 :

 14 Figure 14: Zoom on the profiles of the solution in the real problem (symbols) compared to the homogenized solution (black plain line) for f = 10 Hz (same profile as in Fig. 13); in the real problem, ∆x has been reduced from 2 m (as used in the homogenized problem) to 0.25 m.

Figure 15 :

 15 Figure 15: Effect of the local curvature on the relative positions of two scatterers. R is the local radius of curvature, producing a smallest distance δh ≤ h between the two scatterers (and δh = h for 1/R = 0).

  13). Here, we have considered D = 250 m (R = 160 m), D = 160 m (R = 65 m) and D = 125 m (R = 40 m). Corresponding values of (1 -δh/h) are 0.06, 0.15 and 0.25 respectively, for a maximum allowed value of 0.5. By comparison of the results in Figs. 16-17 with those obtained with a straight line Γ in Figs.

(a) 1 -

 1 δh/h = 0.06(b) 1 -δh/h = 0.15 (c) 1δh/h = 0.25

Figure 16 :

 16 Figure 16: Wavefields for a curved mean line Γ. Expect the shape of Γ, the calculations are identical to those of Fig. 12 for f 0 = 10 Hz. On (a-c), the left panels show the wavefields p r in the real problem (with dotted white lines indicating Γ); the right panels the wavefields p in the homogenized problem (with the white region indicating the homogenized interface.

Figure 17 :

 17 Figure 17: (a-c) Pressure profiles along the centerline x = 600 m from Figs. 16, p r (blue symbols) and p (black lines). (d) Error |p r -p|/|p| as a function of (1 -δh/h), Eq. (46).

Figure B. 18 :

 18 Figure B.18: Incident plane wave (I) impacting a plane homogenized interface, leading to a reflected (R) and a transmitted (T) plane waves.

= ρ 2 C 1

 21 ), and setting the impedance Z = ρ c cos θ and the parameters (α, β)α cos 2 θ + C 2 sin 2 θ , β = R) (1 + iω β) = (1 -iω β) , (T + R) (Z -iω α) = (Z + iω α) . Zβ) (Z + iωα)(1 + iωβ) , T = Z -ω 2 αβ (Z + iωα)(1 + iωβ) , (B.5)

Appendix A. Interface parameters (B, C 1 , C 2 ) for rectangular sound-hard scatterers

The interface parameters involved in the jump conditions (3) have been calculated for rectangular sound-hard scatterers in [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF][START_REF] Marigo | An interface model for homogenization of acoustic metafilms[END_REF]. As they are written in (3), they have the dimension of lengths and they read

and, as a rule of the thumb for the last parameter

for rectangular scatterers of length ϕh, with spacing h and thickness e. Here ϕ is the filling fraction of the inclusion. A different definition is used in [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF], where it refers to the filling fraction of the matrix, so ϕ is replaced here by 1-ϕ. A more accurate estimation of C 2 can be obtained by solving a so-called elementary problem, and a simple script to do so is provided in [START_REF] Marigo | Homogenization models for thin rigid structured surfaces and films[END_REF].

Appendix B. Scattering of a plane wave at oblique incidence on a plane homogenized interface

We consider a plane wave at oblique incidence θ on the thick interface and the problem to solve is [START_REF] Cioranescu | An introduction to homogenization[END_REF]. We want to determine u(x, y, ω) in (45), and u is defined


Olivier Omedes
email: omedes@cadence.com

Michel Robert
email: michel.robert@lirmm.fr

Mohamed Ramdani
email: mohamed.ramdani@eseo.fr

Eseo Angers

Block constraints budgeting in timing-driven hierarchical flow

In this paper, we introduce a new block budgeting algorithm that speeds up timing closure in timing driven hierarchical flows. After a brief description of the addressed flow, block budgeting challenges are detailed. Then, we explain why existing budgeting approaches are not adapted to fulfil these challenges. A new block budgeting algorithm is proposed. In order to derive relevant block constraints, this algorithm analyzes the design flexibility. This Flexibility Aware Budgeting (FAB) approach is then compared to some previous ones. Experiments based on commercial EDA tools and real designs show up to 55 % reduction in hierarchical flow run time and lead to a good flow timing closure.

I. BLOCK BUDGETING CHALLENGES

Any physical synthesis solution has limitations on the size of circuits that can be handled in a single run. "Divide and Conquer" approaches have been introduced to overcome these limitations. In this kind of approach, large designs are sub-divided into smaller synthesizable sub-blocks. Fig. 1 depicts a typical timing-driven hierarchical flow. This flow starts with RTL synthesis and technology mapping. Resulting netlist is assumed to be too big or too complex to meet the specified performance. Thus, the design is partitioned into sub-blocks. These blocks are floorplanned upon the chip die. In fact, the floorplanning step includes design physical partitioning, blocks placement and inter-block net global routing. Then, blocks are optimized. Let's notice that partitioning technique not only permits to implement large designs, but also allows optimizing blocks concurrently, that can be precious to decrease time to market -this represents one more reason to use the hierarchical flow. Optimized blocks are reassembled at the top-level and top optimization is run. If chip constraints are not met, the whole process can be repeated. To implement sub-blocks, EDA tools need constraints. Blocks constraints are computed by the block budgeting step. The budgeting step derives blocks IO constraints from the chip constraints. In order to speed up the hierarchical flow timing closure, the budgeting process has, first to assign feasible blocks constraints, and, second, to ensure that if blocks implementation succeeded, all chip constraints will be met. These two conditions qualify budgets quality. Another key point of block budgeting step stays in its low resources consumption (budgeting that runs slower than flat optimization is not interesting). If various delay budgeting approaches have been proposed to drive cell placement or net routing [START_REF] Sarrafzadeh | A delay budgeting algorithm ensuring maximum flexibility in placement[END_REF] [2] [START_REF] Frankle | Iterative and adaptive slack allocation for performance-driven layout[END_REF], there is only a few academic or industrial works dealing with IO constraints computing. [START_REF] Kuo | Delay budgeting for a timing-closure-driven design method[END_REF] proposes an interresting RTL delay budgeting approach but reported that it needs an accurate area-time projection for each block; that is too much time consuming for our flow. [START_REF] Shi | Hybrid hierarchical timing closure methodology for a high performance and low power dsp[END_REF] also presents an IO budgeting method for the timing-driven hierarchical flow. Unfortunately, the proposed budgeting is only applied evolving from a previous design implementation. Thus, it is usable for a second flow iteration but not to compute IO budgets at a so early stage of the design process. In fact, there are two main issues to resolve when computing IO budgets at this design stage. First, at this stage, design is not optimized yet (just area mapped upon library standard cells). Consequently, current logic structure is likely to largely evolve; this logic flexibility has to be taken into account. Second, IO budgets are not only timing constraints. Indeed, in the same way a cell delay depends on the driven load, an output required time assertion, for example, should depend on the connected load. Thus, block budgeting has to compute both timing and electrical constraints.

In this paper, we propose a block budgeting method that addresses these two issues. Logic flexibility awareness is obtained using a logical effort cell modeling [START_REF] Sutherland | Logical Effort: Designing Fast CMOS Circuits[END_REF] and some simple restructuring or rebalancing algorithms inspired from [7] [8]. The logical effort theory is also used to better correlate timing constraints with block load specifications.

The following section is devoted to this Flexibility Aware Budgeting (FAB) description.

II. FLEXIBILITY AWARE BUDGETING

A. Principles

In timing-driven hierarchical flow, at budgeting stage, design is likely to be only area mapped (or not homogeneously optimized). Cells are just implementing RTL functions and are surely not adapted to the load they will have to drive. Moreover, most of time, their logical structure is balanced and it is not taking into account that, for example, one input is more timing critical than another one. Consequently, at budgeting stage, current structure delays fall far short of the design achievable performance; this initial design implementation is highly flexible. If we run delay budgeting using these initial delays and loads, we are quite sure to produce irrelevant IO constraints. To avoid that, the first FAB step aims at computing more relevant delay and load values for all design components. "more relevant" stands for delay and load values that take into account the optimization possibilities: buffering, restructuring, ... Then, timing analysis is performed. This timing analysis gives information about the amount of slack available for each timing paths -slack is defined as the difference between required and arrival times. To ensure timing closure, these slacks (positive or negative) must be distributed between timing paths components. This slack distribution is also called delay budgeting since it consists to increase or decrease components delays -also called budgets. From an optimization point of view, it's obvious that some components are harder to optimize and should, then, receive more budget. One way to favor these components is to use a weighted slack allocation algorithm. In FAB approach, we implement an IMP like algorithm [START_REF] Youssef | Timing constraints for correct performances[END_REF] -weights are detailed in the following. As we said, in our approach, budgets are not only delays but also load specifications. The very last FAB step computes the block IO constraints from the budgeted design.

B. Logical Effort Theory

In fact, the very first FAB step consists in standard cells characterization using the logical effort theory. The concept of logical effort has been proposed to model the load-independent gate delay:

d abs τ¡ d τ ¢ p £ f ¤ τ ¢ p £ g¡ h¤ (1)
This model divides the gate delay into two parts. First part is the load independent parasitic delay p (its major contribution is capacitance of transistors source / drain). Second part ( f g¡ h) is the effort delay. The logical effort g depends only on the gate topology and on its ability to produce output current. The electrical effort h (also called gain) is defined to be:

h C out ¥ C in
Where C out is the connected load value and C in , the gate input capacitance. Finally, τ is a technology dependent factor. All these values are easy to characterize for simple gates (inverter, nand, nor, xnor...). The logical effort theory also defines the branching effort b as the ratio of the total output capacitance to the capacitance which is on the path being analyzed. Given a string of N g monofanout (b = 1) gates, the path delay is:

D ∑ N g p i £ g i h i P £ ∑ N g g i h i
This string, or path, is also characterized by a path logical effort G and a path electrical effort H :

G ∏ g i H ∏ h i C in¦ path C out ¦ path
The GH product is also called path effort and noted F. [START_REF] Sutherland | Logical Effort: Designing Fast CMOS Circuits[END_REF] explains that D is minimized when F is equally spread between stages : f

Ng § F D N g f £ P (2)
Unfortunately, the number of stages N is very likely to evolve during optimization, and thus optimal effort f computation can be fairly ineffective. Indeed, on hard to optimize paths, optimization process will surely add some buffering components (inverters or buffers). Since, buffering components logical effort is defined as 1, we can write:

N N g £ N bu f f N § F D N ¡ f £ P £ ¢ N ¨N g ¤ © p bu f
N bu f and p bu f are respectively the number of added buffering components, and their unit parasitic delay. Differentiating with respect of N and setting the result to zero gives the following displaymath for f to minimize path delay:

f ¢ 1 ¨l n ¢ f ¤ ¤ £ p bu f 0 ( 3 )
(3) can be solved using numerical method (Newton-Raphson for example) and obtained f value permits to define an optimal electrical effort for each type of gate:

ĥ f ¥ g
This result is fundamental for our approach since ĥ will be used for delays and drive strengths initialization. Following table shows some experimental results for a 90 nm technology (τ 7¡ 2 ps).

¡

Buffer trees: a buffer tree is a set of interconnections and buffering components. In following algorithms, a simple mono-fanout net is considered like a particular buffer tree. FAB initialization process depends on a standard load value C re f that can be either computed during the library characterization or user specified. In practice, its value is similar to the 1X inverter input capacitance. FAB initialization starts with the less flexible design instances which are IPs (Intellectual Property Blocs or small custom-designed blocks) and registers (sequential elements). IPs and registers timing arc delays are set to their propagation delay values when outputs are driving C re f . Budgeted load C b of IPs and registers outputs is set to C re f . Budgeted load of IPs and registers inputs is set to their input capacitance value and is marked as fixed. Then, the flexibility aware initialization (FAI) step can start.

FAI is based upon two design transforms that are quite easy to anticipate. The first one consists in gates and cones fan-in ordering [START_REF] Paulin | Logic decomposition algorithms for the timing optimization of multi-level logic[END_REF] [START_REF] Singh | Timing optimisation of combinational logic[END_REF]. Fig. 2 illustrates such a fan-in ordering on a small cone. Ordering starts when budgeted arrival times (BAT) and loads at cone inputs have been computed. It begins by remapping F C using only inverter, NAND and XNOR gates (cf. Fig. 2-b) -NOR gate is voluntary not used because of its low driving capability. Then, the structure is unbalanced in order to get a minimal BAT at the cone output (cf. Fig. 2-c). Ensure: sink ¢ SinkOfCone ¢ C ¤ treated Remap F C (Fig. 2-b) and rebalance it (Fig. 2-c).

? / ? 0 / 1 0 / 2 Budgeted AT / Load 3 / 1 0 / 1 0 / 2 3 / 1 ? / ? a b 0 / 1 0 / 2 3 / 1 6 / 2 c 2 / 2 : = = h d Nand2 1 / 1 : = = h d Inv
d ¢ v £ sink¤ ¢ ∑ v¤ sink d C b ¢ sink¤ ¢ min v¥ π¦ § C ¨© C b ¢ v¤ © ∏ v¤ sink ĥ
The second transform is about buffering. It consists in 2 phases. The first phase computes a balance buffer tree to drive the whole B buffer tree fanout (number of leaves). Budgets are computed for all buffer tree branches in the same time. On the contrary, in the second phase, leaves are treated one after the other. This second phase aims at adding some budget to branches that are driving large fixed loads (such as IPs input capacitance). Some budgets are also added in order to do block interconnections buffering (interconnection capacitances are estimated from the initial floorplan). Algorithm 2 and fig. 3 These two transforms are combined and evaluated by traversing the graph in topological order (BFS) from primary inputs to primary outputs (this ensures that when treating node v, all π ¦ ¢ v¤ nodes have already been treated). It's important to remark that these transform are not really applied (netlist is not edited) but just evaluated to get some budget values. This allows a significant computing gain. For this reason, these transforms are qualified of "virtual". FAB initialization produces a graph which is annotated with timing and load initial budgets.

D. FAB allocation

Timing analysis on this budgeted graph shows paths with both negative and positive slacks. Negative slack identifies hard to optimize paths whereas positive one characterizes more flexible paths where budgets (both d and C b ) can be increased. As explained in FAB principles section, path slack S π is distributed among path arcs using the IMP algorithm. IMP is an iterative algorithm which adds a timing budget ∆d

Ensure: π

¢ src¤ treated. if Fanout ¢ B ¤ © C re f C b ¢ src¤ then bufferingDelay ¢ 0 tempC b ¢ C b ¢ src¤ ¥ Fanout ¢ B ¤ else effort ¢ Fanout ¢ B ¤ © C re f ¥ C b ¢ src¤ N b ¢ log ĥbu f ¢ e f f ort ¤ bufferingDelay ¢ N b ¢ g bu f © N b § e f f ort £ p bu f ¤ (cf. Eq. 2) tempC b ¢ C re f end if /* End of first phase */ for each v ¡ π ¢ src¤ do d ¢ src £ v¤ ¢ bufferingDelay if isFixed ¢ C b ¢ v¤ ¤ then effort ¢ C b ¢ v¤ ¥ tempC b if effort ¢ 1 then N b ¢ log ĥbu f ¢ effort¤ d ¢ src £ v¤ ¢ N b ¢ g bu f © N b § effort £ p bu f ¤ (cf. Eq. 2) end if else C b ¢ v¤ ¢ tempC b end if end for /* End of second phase */
Where w e is the e arc weight, Π e is the set of all paths traversing e and, S π and W π are respectively π path slack and total weight. A proof of IMP convergence to null slack on final budgeted timing graph is given in [START_REF] Youssef | Bounds on net delays for VLSI circuits[END_REF]. The main difficulty in slack allocation is to find out relevant weights for arcs. In FAB context, weights have to qualify arcs optimization hardness. Given the way initial arcs delays have been computed (using virtual restructuring and buffering), they give a good idea of optimization timing needs. Consequently, we use these initial delays as weights for the IMP allocation. IPs (or fixed blocks) and sequential elements receives a null weight since their flexibility is also nearly null; indeed IPs are impossible to remap and sequential elements are also fixed since we fixed the load they will have to drive. If this slack allocation is quite typical, FAB not only allocates timing budgets, it also modifies node budgeted capacitance according to the added delay budget. Indeed, if we add delay to an arc, it's obvious to think that the load it can drive is also increased. So, if we differentiate (1) with respect to C out , we obtain:

∂d ∂C out g C in £ ∆C out C in g ∆d
Consequently, FAB allocation is ran as described by algorithm 3. At the end of this allocation, we get a full budgeted timing graph with null slack everywhere and each node is also constrained by a maximal drivable capacitance.

W π d ¢ e¤ ¢ ∆d if isNull ¢ ∆d¤ then Weight ¢ e¤ ¢ 0 setFixed ¢ C b ¢ sink ¢ e¤ ¤ ¤ end if if isNotFixed ¢ C b ¢ sink ¢ e¤ ¤ ¤ then C b ¢ sink ¢ e¤ ¤ ¢ ∆d © C b ¢ source ¢ e¤ ¤ ¥ g e end if
¢ BAT ¢ p¤ ¤ set max input capacitance ¢ C b ¢ p¤ ¤ else set required time ¢ BAT ¢ p¤ ¤ set output load ¢ C b ¢ p¤ ¤ end if end for
Arrival Time and is computed by propagating budgeted delays. Algorithm 4 gives optimization constraints to facilitate block integration and top level timing closure. Indeed, for an input port, FAB ensures that if the max input capacitance is not exceeded, the port arrival time will be at least the one asserted. In the same way, block optimization have to implement block such that logic will be able to drive the output load with respect to the required time assertion.

III. VALIDATION

A. Experimental Setup and Metrics

To qualify block IO budgets quality, we need to remember two important points. First, block optimization is a very chaotic step for which small constraints changes are likely to produce very different block implementations. Second, block constraints set is quite large and constraints are highly linked between them. Hence, if we want to compare budgeting algorithms results just looking at blocks implementations or constraints, we are likely to fail. In fact, the most significant approach to evaluate budgeting methods is to use them in

B. Some other used approaches

In the results section, we compare FAB approach with two other used block budgeting algorithms. The first one is called critical-path budgeting (CPB). Its principle is to treat block IO ports one after other and, to adjust current port arrival time according to delays of the critical path which crosses it. So, block IO port p budgeted arrival time is computed as follow:

BAT ¢ p¤ CurrentArrivalTime ¢ p¤ £ Slack ¢ p¤ d le f t d le f t £ d right
Where d le f t and d right are respectively delays of the left and right part of the most critical path traversing p. This budgeted arrival time is then used in algorithm 4 to assert block constraints. Capacitance assertions are computed using current load value for outputs load and driving cell load limit for inputs max capacitance. Of course, due to completely not optimized design state, current values (arrival times, loads and driving cells) are surely irrelevant and thus, resulting IO budgets quality is doubtful. However, this approach has the advantage to be really fast (runtime wise).

The second approach, called IMP T, is based upon IMP algorithm but doesn't include a FAI like step. For reasons explained above, it would be ineffective to use current design delays as weights. However, it could be a good idea to set weights using the logical effort theory. Concerning loads assertions, a common way to avoid asserting huge values is to truncate them (limits can be specified by users). Fig. 5 shows truncation examples. This kind of truncation is part of IMP T approach.

C. Results

Benchmarks specifications are given in table I. As we can see in the "Initial Slack" column, initial designs slack falls 4 flow results when using either CPB or IMP T budgeting algorithms to drive one pass of hierarchical optimization. Runtime and memory wise, there is no doubt that the "divide and conquer" approach is faster and lighter than the flat one. Even if the timing closure is not obtained at the end of this first pass, the average runtime gain of about 66 % lets us with time to run a second iteration and to achieve this timing closure. The average memory need, which is about 35 % less than the one observed when using the flat flow, will not be increased in following iterations. Now, let's compare the CPB and IMP T timing results. The WBSlack column shows the worst block slack measured after blocks optimization. The -0.94 ns WBSlack obtained when using the CPB algorithm shows well that this approach assigns some irrelevant block constraints. On the contrary, blocks constraints derived by IMP T budgeting seem more relevant (or at less more feasible) since the average WBSlack is only -0.27 ns -drivers and loads truncation is very helpful for that. Of course, having more feasible blocks constraints is not sufficient to ensure top level timing closure; these constraints have also to allow a feasible top level optimization. However, if we look at FTSlack columns for both CPB and IMP T algorithms, we can conclude that IMP T computes better blocks constraints than CPB. Indeed, at the end of this first hierarchical optimization pass, the top level slack is nearer from timing closure when using the IMP T budgeting algorithm (IMP T average top slack is about -0.40 ns whereas CPB one is about -0.94 ns). Hence, IMP T budgeting leads to a better implementation to start a second iteration. Now, let's look at the results obtained when using our new budgeting approach (Table III). Resources wise, the average gain stays similar to previous budgeting algorithms. Concerning the timing closure, load and drive constraints derivation -in parallel with timing ones -as done in FAB seems to be more relevant than simple truncation. Indeed, FAB WBSlacks stand out better than IMP T ones (average gain is about 0.18 ns) and, as expected, this improves top-level optimization results (average gain is about 0.15 ns). When we look closer at why timing closure was not obtained at the end of this first pass, we found out that some drive sizing and some buffering transforms was hardly impossible to anticipate since cells were not placed yet at budgeting step. However, for all testcases, the top-level timing closure has been obtained by just running a quick top-level optimization after the first flow iteration in order to fix these drive / buffering problems -this is not always true with CPB and IMP T approaches for which deeper optimizations are needed (these results have not been reported in this paper). It's important to notice that, even with this additional top-level optimization, average runtime stays about 30 % less compared to flat flow runtime (memory gain is about 20 %).

IV. CONCLUSIONS

In this paper, we introduced a new budgeting approach that speeds up timing closure of a timing driven hierarchical flow. The initial budget allocation step of this new approach (FAI) allowed us to quickly identify paths that were "really" critical (following optimizations iterate on these paths). FAB IO constraints computing (BRun) is about 10 to 20 times slower than CPB and about the same than IMP T, but this runtime stays lower than 10 % of FRun (e.g., about 15 minutes for TEST5). In addition, if we look at the results obtained for 5 commercial designs, it is obvious that taking more time to derive relevant IO constraints leads to a better timing performance predictability and to a faster hierarchical flow timing closure.

In the future, by integrating FAB approach to a physical prototyping tool, we should be able to take into account some more physical parameters such as blocks congestion and aspect ratios. This could also be helpful to compute more relevant weights or budgets for blocks interconnections.

Fig. 1 .

 1 Fig. 1. Typical hierarchical flow

Fig. 2 .Algorithm 1

 21 Fig. 2. C Restructuring

Fig. 3 .

 3 Fig. 3. B Buffering

Algorithm 2

 2 Quick B Buffering Evaluation Require: src ¢ Source ¢ B ¤ already treated.

Fig. 4 .

 4 Fig. 4. Validation metrics

Fig. 5 .

 5 Fig. 5. Truncation of block input port driver and block output port load.

  ¤ which represents the circuit timing structure. V denotes the set of vertices or nodes -nodes are circuit and cells IO. E denotes the set of edges or arcs which represent timing relations between nodes. Each edge e is annotated with a delay value d ¢ e¤ which corresponds to its delay budget. Each node v is associated with a budgeted load C b Cones: a cone is a set of mono-fanout logic and is completely identified by a set of cone inputs and one cone output. Cone C implements a boolean function F C .

	¢
	v¤ which corresponds
	to the maximum drivable load of this node. π ¦ ¢ v¤ and π	¢ v¤
	respectively denote the set of v predecessors and successors.
	In the timing graph, we distinguish some interesting structures
	that are:
		Cell	g	ĥ	p	d ps
		Inverter	1.01	3.71	1.16	4.90
		Nand 2	1.41	2.62	2.41	6.10
		Xnor 2	1.62	2.31	12.59	16.33
		C. FAB Initialization process
		FAB walks on an acyclic directed timing graph G

¢

V E

¡

  typical timing-driven hierarchical flow. Then, metrics like optimization runtime, and slacks of blocks and top level can be used to appreciate budgeting results. Metrics we use are defined on fig.4. Of course, this flow and these metrics favor rapid to close methods but this is what we are looking for.

	Tests (Hierarchical) Frequencies : Nb Gates : 100k 5M 50 200 MHz	In	IP
	Techno : 250, 180, 130 nm
	Timing constraints	VHDL / Verilog			Out
	RTL Synthesis / Partitionning
			Budgeting	BRun
	-BRun : Budgeting Runtime Metrics Definition :	Concurrent Block Optimization	WBSlack	FMem FRun
	-FRun : Total Flow Runtime		Top Level
	-FMem : Total Flow Memory	Optimization
	-WBSlack : Worst Block Slack	Timing Closure ?
	-FTSlack : Final Top Slack		FTSlack
	DCIS 2004				-933 -	ISBN 2-9522971-0-X

a

  Table I also summarizes flat flow results. Table II details fig.

TABLE III NEW BLOCK BUDGETING RESULTS -FAB ALGORITHM.

 III

	Design		FAB Hier. Flow Results
	Name	FRun (s)	FMem (Mb)	WBSlack (ns)	FTSlack (ns)
	TEST1	800	107	0.00	-0.04
	TEST2	3300	456	-0.08	-0.33
	TEST3	2900	930	-0.03	-0.16
	TEST4	14000	690	-0.12	-0.25
	TEST5	11700	1347	-0.25	-0.47
	∆	-60 %	-31 %	-0.09 ns	-0.25 ns