Raphaèle Herbin
email: herbin@cmi.univ-mrs.fr

Florence Hubert
email: fhubert@cmi.univ-mrs.fr

Benchmark on Discretization Schemes for Anisotropic Diffusion Problems on General Grids

Keywords: Anisotropic medium, diffusion process, finite volume schemes, benchmark

We present here a number of test cases and meshes which were designed to form a benchmark for finite volume schemes and give a summary of some of the results which were presented by the participants to this benchmark. We address a two-dimensional anisotropic diffusion problem, which is discretized on general, possibly non-conforming meshes. In most cases, the diffusion tensor is taken to be anisotropic, and at times heterogeneous and/or discontinuous. The meshes are either triangular or quadrangular, and sometimes quite distorted. Several methods were tested, among which finite element, discontinous Galerkin, cell centred and vertex centred finite volume methods, discrete duality finite volume methods, mimetic methods. The results given by the participants to the benchmark range from the number of unknowns, the errors on the fluxes or the minimum and maximum values and energy, to the order of convergence (when available).

Introduction

The aim of this benchmark is to provide a number of test cases in order to test the properties (convergence, robustness...) of existing discretization schemes for anisotropic diffusion problems using general grids. In all test cases except test 8, the domain Ω is the unit square. The boundary of the domain is divided into ∂Ω = Γ D ∪Γ N where Dirichlet (resp. Neumann) boundary conditions are given on Γ D (resp. on Γ N ). The considered diffusion problem is formulated as:

-∇•(K∇u) = f on Ω, [1] u = ū on Γ D , [2] K∇u • n = g on Γ N , [3]
where K : Ω → IR 2×2 is the diffusion (or permeability) tensor, f the source term, ū and g the Dirichlet and Neumann boundary conditions, and n denotes the outward unit normal vector to Γ N . For each test case, we propose some meshes which are used for the discretization of Problem [1] by the various schemes. These meshes are depicted in section 15 below.

The participating schemes and teams

Even though the benchmark is associated with the FVCA5 conference, the call for submission was by no means restricted to finite volume schemes, and indeed several types of schemes were submitted:

Cell centred schemes

• CMPFA: Compact-stencil MPFA method for heterogeneous highly anisotropic second-order elliptic problems, by S. Mundal, D. A. Di Pietro and I. Aavatsmark, [?].

• FVHYB: A symmetric finite volume scheme for anisotropic heterogeneous second-order elliptic problems, by L. Agelas and D. A. Di Pietro, [?].

• FVSYM: Numerical results with two cell-centered finite volume schemes for heteregeneous anisotropic diffusion operators, by C. Le Potier, [?].

• SUSHI-P or SUSHI-NP (which are cell centred in their full barycentric version): SUSHI: A Scheme Using Stabilization and Hybrid Interfaces for anisotropic heterogeneous diffusion problems, by R. Eymard, T. Gallouët and R. Herbin, [?].

Control volume finite element schemes

• CVFE: Numerical simulation for the anisotropic benchmark by a vertex-centred finite volume method, by M. Afif and B. Amaziane, [?].

Discontinous Galerkin schemes

• DG-C: The Compact Discontinuous Galerkin Method for Elliptic Problems, by A. Dedner and R. Klöfkorn, [?].

• DG-W: A discontinuous Galerkin flux for anisotropic heterogeneous secondorder elliptic problems, by D. A. Di Pietro and A. Ern, [?].

Discrete duality finite volume schemes

• DDFV-BHU: The DDFV "discrete duality finite volumes" and m-DDFV schemes, by F. Boyer and F. Hubert, [?].

• DDFV-HER: Numerical experiments with the DDFV method, F. Hermeline, [?].

• DDFV-MNI: Some MPFA methods of DDFV type, by I. Moukouop Nguena and A. Njifenjou, [?].

• DDFV-OMN: Tests with the Discrete Duality Finite Volume method, by P. Omnes, [?].

Finite element schemes

• FEP1, FEP2, FEQ1, FEQ2: A Galerkin Finite Element Solution, by G. Ansanay-Alex, B. Piar and D. Vola, [?].

Lattice Boltzmann schemes

• LATTB Using Lattice Boltzmann scheme for anisotropic diffusion problems, by F. Dubois, P. Lallemand and M. M. Tekitek, [?].

Mixed or hybrid methods

• MFD-BLS: Mimetic finite difference method, by K. Lipnikov, [?].

• MFD-FHE: Numerical investigation of a mimetic finite difference method, by B. Flemisch and R. Helmig, [?].

• MFD-MAN: The mimetic finite difference method, by G. Manzini, [?].

• MFD-MAR: A mimetic finite difference method, S. Marnach, [?].

• MFE: Mixed finite element (Raviart Thomas) programmed by the benchmark organizers for lack of a submission.

• MFV: Use of the mixed finite volume method, by C. Chainais-Hillairet, J. Droniou and R. Eymard [?].

• SUSHI-P or SUSHI-NP (which are hybrid in their full hybrid version): SUSHI: A Scheme Using Stabilization and Hybrid Interfaces for anisotropic heterogeneous diffusion problems, by R. Eymard, T. Gallouët and R. Herbin, [?].

Nonlinear schemes

The schemes are nonlinear in order to ensure the positivity of the scheme (that is, if the right hand side is positive then the solution is positive) or the discrete maximum principle (that is, if the linear system stems from the discretization of an elliptic equation satisfying the maximum principle, then its solution is also bounded by the bounds of the continuous system).

• FVPMMD: Numerical results with two cell-centered finite volume schemes for heteregeneous anisotropic diffusion operators, by C. Le Potier, [?].

• NMFV: Nonlinear monotone finite volume method, by D. Svyatskiy, [?].

Of course, the above choice of categories is neither exhaustive nor unique. In fact, most of these categories intersect: schemes are not so easy to classify, and some schemes are known to be identical in special cases. We refer to the above-cited papers for the details of the schemes and their implementation. Our purpose is to give here a synthesis of the results presented by the participants and raise some questions on the future trends in this field of research.

Nature of the results and notations

Some of the output results are given for all tests, and others only for those tests for which we know of an analytical solution. Here are both lists.

For all tests

• nunkw number of unknowns • nnmat number of non-zero terms in the matrix • sumflux the discrete flux balance, that is: sumflux = flux0+ flux1 + fluy0 +fluy1sumf, where flux0, flux1, fluy0, fluy1 are the outward numerical fluxes at the boundaries x = 0, x = 1, y = 0 and y = 1 (for instance flux0 is an approximation of -x=0 K∇u • n ds) and sumf = K∈T |K|f (x K ) where x K denotes some point of the control volume K. Let us note that for vertex methods such as the finite element schemes FE or the CVFE scheme, the points x K are the vertices of the mesh. For most cell centred schemes, x K is chosen as the centre of gravity for most schemes. Some papers concerning DDFV schemes (DDFV-OMN, DDFV-BHU) report results for values at two different points, but we shall only present here the results with x K the centre of gravity of the primal mesh.

The residual sumflux is a measure of the global conservativity of the scheme; hence, it is expected to be of the order of the machine precision for all schemes which are based on a local conservation of fluxes. In fact, it was very useful in debugging codes; however, for the sake of conciseness, we choose not to display its value here since it is indeed of the machine precision order for FV type schemes and of the mesh size order for purely Galerkin finite element type schemes (FE).

• umin: value of the minimum of the approximate solution.

• umax: value of the maximum of the approximate solution.

• ener1, ener2, where ener1 and ener2 are approximations of the energy following the two expressions (which are identical in the continuous setting): ener1, ener2

E 1 = Ω K∇u • ∇udx, E 2 = ∂Ω K∇u • n udx [4]
Even though E 1 and E 2 converge to the same value on fine grids, there can be a noticeable difference between E 1 and E 2 on the coarsest meshes. Note also that we get an upper bound of the continuous value E 1 = E 2 by using a conformal finite element method and lower bound with the mixed finite element (Raviart-Thomas) method.

When the analytical solution is known and the mesh refined

Let us denote by u the exact solution, by T the mesh and by u T = (u K ) K∈T the piecewise constant approximate solution.

• erl2, relative discrete L 2 norm of the error, for instance:

erl2 =     K∈T |K|(u(x K ) -u K ) 2 K∈T |K|u(x K ) 2     1 2
.

• ergrad relative L 2 norm of the error on the gradient, if available

• ratiol2: for i ≥ 2, ratiol2(i) = -2 ln(erl2(i)) -ln(erl2(i -1)) ln(nunkw(i)) -ln(nunkw(i -1))
• ratiograd, for i ≥ 2, same formula as above with ergrad instead of erl2.

• erflx0,erflx1, erfly0,erfly1 relative error between flux0, flux1, fluy0, fluy1 and the corresponding flux of the exact solution:

erflx0 = flux0 + x=0 K∇u • n x=0 K∇u • n
(except for the fluxes at y=0 and y=1 for case test 2 -numerical locking -because these are zero).

• erflm L ∞ norm of the error on the meanvalue of the flux through the edges of the mesh, if available (give the definition of numerical flux

(K∇u • n) T ) erflm = max 1 |σ| σ (K∇u • n -(K∇u • n) T ) , σ edges of T .
• ocvl2 order of convergence of the method in the L 2 norm of the solution as defined by errl2 with respect to the mesh size:

ocvl2 = ln(erl2(imax)) -ln(erl2(imax -1)) ln(h(imax)) -ln(h(imax -1))
where h is the maximum of the diameter of the control volume.

• ocvgradl2 order of convergence of the method in the L 2 norm of the gradient as defined by ergradl2 with respect to the mesh size, same formula as above with ergrad instead of erl2.

Test 1: mild anisotropy

A homogeneous anisotropic tensor is first considered: K = 1.5 0.5 0.5 1.5 .

Test 1.1: smooth solution, mesh1 (triangular mesh) and mesh4 (distorted quadrangular mesh)

In this first test, we consider a very regular solution: u(x, y) = 16x(1-x)y(1-y), we set f = -∇ • (K∇u), and we consider Dirichlet boundary conditions: Γ D = ∂Ω, Γ N = ∅, and ū = u| ∂Ω . The schemes were first tested on a series of "regular" triangular meshes: mesh1, a picture of which is given in Figure 13. We give in Figure 1 the L 2 norm of the error versus the number of unknowns (upper left), the square root of the number of non zero matrix terms (upper right) and the step size h (lower left) for all schemes for meshes mesh1. As can be seen on Figure 1, all schemes show an order 2 convergence except for the two higher order schemes FEP2 and DG-C, which are order 3. The schemes DDFV-HER, FEP2 and DG-C (see Figure 2) have an order 2 convergence with respect to the L 2 norm of the gradient; all other schemes show an order 1 convergence.

Note that when calculating the error with respect to nunknw and nnmat, the order of convergence remains the same (as expected) but the relative position of the schemes vary, as may be seen on the graphs of Figure 1. The behaviour or the L 2 norm of the gradient with respect to the three previous variables is depicted in Figure 2.

The same test was then carried out on two distorted quadrangular meshes mesh4_1 and mesh4_2 depicted in Figure 15. From Table 1 below, we can see that almost all schemes satisfy the discrete maximum principle on this test and meshes. The DG schemes DG-C and DG-W provide small negative values while the nonlinear (positive) scheme NMFV exceeds the maximum value 1. Note that some schemes include the boundary conditions in the calculation of the minimum and maximum values and others not. In the second case, if these minimum and maximum values are not equal to the boundary conditions, they can be taken as a measure of the accuracy of the scheme.

Test 1.2: mesh1 (triangular mesh) and mesh3 (locally refined non-conforming rectangular mesh)

We consider here the following exact solution: u(x, y) = sin ((1 -x)(1 -y)) + (1 -x) 3 (1 -y) 2 , with the right hand side f = -∇ • (K∇u), boundary conditions (Γ D = ∂Ω, Γ N = ∅) and ū = u| ∂Ω . We use both the triangular mesh mesh1 and the locally refined non-conforming rectangular mesh mesh3, which is depicted in Figure 14. The L 2 norm of the error and the gradient with respect with nunknw, nnmat and h are given in Figures 3 and4 for the triangular mesh mesh1 and, in Figures 5 and6 for the locally refined non-conforming rectangular mesh mesh3.

As can be seen on Figure 3, the orders of convergence of the schemes are the same as in the case of Test 1.1. Similar results are observed on the locally refined non-conforming rectangular mesh mesh3. 2,4 and 6 that the L 2 norm of the gradient error is lower with the DDFV schemes than for the other second order schemes.

Note from Figures

Test 2: Numerical locking [?, ?]

In this test, the permeability tensor is defined by: K = 1 0 0 δ , with δ = 10 5 or 10 6 . The exact solution is taken to be u(x, y) = sin(2πx)e -2π √ 1/δy ; note that since δ is large, the solution is almost constant in the y variable. We consider Problem This test was mentioned in [?] as causing some numerical locking problems for finite volume and finite element schemes, and in fact, some schemes could not give the solution. We give in Tables 2 and3 the minimum and maximum values obtained with those schemes which passed the test, and on what grid they are obtained, along with the convergence order for the error of the solution and its gradient. The behaviour of the P 1 finite element scheme seems good except on the very first grid where, even though the maximum principle is satisfied, the maximum and minimum values are close to 0, and quite far from the exact ones; this is also the case for the schemes DG-W SUSHI-NP (but this latter scheme does not converge well on finer grids). The P 2 finite elements also give good results, even if the rate of convergence is a bit under the theoretical rate on the last grids. Note that the scheme DDFV-OMN satisfies the discrete maximum principle for both δ = 10 5 and δ = 10 6 . Most schemes seem to encounter some difficulties on the first two (for δ = 10 5 ) or three (for δ = ) coarsest meshes (locking effect). Note also the very similar behaviour of the schemes MFE, MFD-FHE, MFD-BLS, MFD-MAN, MFD-MAR and SUSHI-P on this test case.

Test 3: oblique flow

This test case represents a flow with boundary conditions such that the pressure driven flow "would like" to go from vertex (0,0) to vertex (1,1), but is impeded by a homogeneous anisotropic tensor with high permeability in a direction at 40 degrees from the horizontal and low permeability in the orthogonal direction. This test case is inspired by a talk given by I. Aavatsmark in Paris in December at GDR MOMAS. The permeability tensor is:

K = R θ 1 0 0 δ R -1 θ ,
where R θ is the rotation of angle θ = 40 degrees and δ = 10 The source term is equal to zero (f = 0). We use a uniform rectangular family of meshes mesh2_i. Dirichlet boundary conditions are considered, Γ D = ∂Ω with ū a continuous and piecewise linear function defined by:

ū(x, y) =        1 on ((0, .2) × {0.} ∪ {0.} × (0, .2) 0 on ((.8, 1.) × {1.} ∪ {1.} × (.8, 1.) 1 2 on ((.3, 1.) × {0} ∪ {0} × (.3, 1.) 1 2
on ((0., .7) × {1.} ∪ {1.} × (0., 0.7)

The solution features a Z across the y = x axis; it is depicted in Figure 7 (approximate solution obtained by the scheme DDFV-BHU on a very fine grid). The maximum principle is not always easy to verify for such a solution and is violated, on coarse and fine grid, by most DDFV implementations, the DG schemes and the scheme FVHYB (see Table 4). Note also that, as previously mentioned, the FEQ2 energy, equal to 2.42E -01, is an upper bound of the exact energy, while the MFE energy, equal to 2.41E -01, is a lower bound of this energy. On the coarse grid, the scheme which gives the closest energy to these values is LATTB (see Table 4). However, it may be seen in Table 5 that the scheme LATTB is not so good with respect to the fluxes, for which the scheme MFD-FHE seems to be the most accurate. From this table, we can notice that most of the schemes give a reasonable approximation of the outward flux on the boundary.

Test 4: vertical fault

The medium considered here is a pile of anisotropic layers with a fault in the middle, which leads to a discontinuity of the layers at x = .5. Each geological layer is meshed with one layer of discretization cells only. The domain Ω is decomposed as

Ω = Ω 1 ∪ Ω 2 , with Ω 2 = Ω \ Ω 1 , with Ω 1 = Ω 1 ∪ Ω r 1 ,
and It is described in Figure 8 where Ω 1 is in black and Ω 2 in white. The diffusion tensor K is anisotropic and heterogeneous, and is given by:

Ω
(0,0) (1,1) (0,1) (1,0) Figure 8.
K = α 0 0 β , with          α β = 10 2 10 on Ω 1 , α β = 10 -2 10 -3 on Ω 2
A Dirichlet boundary condition is imposed: ū(x, y) = 1 -x on Γ D = ∂Ω (so that Γ N = ∅) and the right hand side is assumed equal to zero f = 0. Note that the exact solution of this problem should stay between 0 and 1; it is therefore interesting to see whether the schemes respect these bounds (discrete maximum principle), especially on coarse meshes (which are used in oil industry, for instance). Problem [1] is discretized using the non-conforming rectangular mesh mesh5 (see Figure 16), the conforming 20 × 20 square mesh mesh5 reg and a reference mesh, for instance the 320 × 320 square mesh mesh5 ref . From Table 8, we see that the energies given by the FEQ1 scheme and the MFE scheme on the reference mesh are both equal to 43.2, which means that this is also the value of the continuous energy. On the fine mesh, all schemes give an approximate energy which is close to this value, and the results on the fluxes are also quite homogeneous. If we now look at what are the schemes that give the energy which is closest to 43.2 on the coarse mesh (Table 6), we find that the closest is given by DG-W (error of less than 0.5 %), and that the outward fluxes (Table 7) given by DG-W are also close to the reference ones (Table 9). However, neither the upper nor the lower bound is respected by the scheme. This illustrates the well-known fact that precision and robustness are difficult to obtain simultaneously. Other schemes (DDFV-BHU, DDFV-OMN, FVHYB, MFD-MAR, SUSHI-NP) give an energy and fluxes which are reasonably close to the reference ones (some schemes are better for the energy, other better for

Test 5: heterogeneous rotating anisotropy

This test is inspired from [?, ?]. The permeability tensor is a rotating anisotropic tensor:

K = 1 (x 2 + y 2 )
10 -3 x 2 + y 2 (10 -3 -1)xy (10 -3 -1)xy x 2 + 10 -3 y 2 , and we consider the following smooth exact solution:

u(x, y) = sin πx sin πy, f = -∇ • (K∇u).
Problem [1] is then solved with non-homogeneous boundary conditions: Γ D = ∂Ω, Γ N = ∅ and ū(x, y) = sin πx sin πy, and approximated on the uniform rectangular mesh2.

We give in Figure 9 the L 2 norm of the error versus the number of unknowns (upper left), the square root of the number of non-zero matrix terms (upper right) and the step size h (lower left) for all schemes for meshes mesh2. As can be seen in Figure 9, all schemes show an order 2 convergence except for the two higher order schemes FEP2 and DG-C, which are order 3. As for the L 2 norm of the gradient is concerned, all the schemes, even higher order schemes, show an order 2 convergence, with an error with respect to nunkw, nnmat and h minimized by the DDFV schemes (see Figure 10). All schemes converge when the mesh is refined, but some schemes violate the discrete maximum principle, as may be seen in Table 10. Note that the schemes MFE and FVHYB violate it even on the fine grids.

Test 6: oblique drain

This test case represents a situation which is encountered in underground flow engineering where an oblique drain consisting in a very permeable layer concentrates most part of the flow; this drain is meshed with only one layer of discretization cells. In the case of a pressure gradient driven transport, as often described in reservoir engineering, it seems important that the discretization cells consist in only one homogeneous material: numerical experiments show that otherwise the solution may be badly approximated. Here we consider the steady case, but wish to verify that the outward fluxes are as close as possible to the exact values for the meshes considered here, both for the conforming and non-conforming meshes.

The domain Ω is composed of 3 subdomains:

Ω 1 = {(x, y) ∈ Ω; φ 1 (x, y) < 0}, Ω 2 = {(x, y) ∈ Ω; φ 1 (x, y) > 0, φ 2 (x, y) < 0}, Ω 3 = {(x, y) ∈ Ω; φ 2 (x, y) > 0}, with φ 1 (x, y) = y -δ(x -.5) -.475 and φ 2 (x, y) = φ 1 (x, y) -0.05.
We take the slope of the drain δ = 0.2 and define the exact solution and the source term by: u(x, y) = -x -δy, on Ω, f = -∇(K∇u),

where the permeability tensor K is such that its principal axes are parallel and perpendicular to the drain:

K = R θ α 0 0 β R -1 θ ,
with θ such that δ = tan θ and:

α β = 10 2 10 on Ω 2 and α β = 1 10 -1 , on Ω 1 ∪ Ω 3 .
Problem [1] is considered with non-homogeneous Dirichlet boundary conditions: Γ D = ∂Ω, Γ N = ∅ and ū(x, y) = -x -δy, and approximated on the conforming mesh: mesh6 and the non-conforming mesh mesh7.

As the exact solution is affine on the whole domain, only a few schemes (CVFE, DDFV-HER, MFV, NMFV) are not exact on such a problem. Moreover, they all satisfy the discrete maximum principle. The interest of this case lies in the approximation of the interface fluxes, which could be used for instance for the approximation of a coupled system of diffusion and convection equations.

Test 7: oblique barrier

This test case is similar to the Test 6, except that we now have to deal with a barrier, and the aim is that the scheme should respect this barrier as well as the outward fluxes.

We take the same geometry as test 6 above, with the slope of the drain δ = 0.2.

We take the exact solution to be

u(x, y) =        -φ 1 (x, y) on Ω 1 , -φ 1 (x, y)/10 -2 on Ω 2 ,
-φ 2 (x, y) -0.05/10 -2 on Ω 3 , and f = -∇•(K∇u), where the permeability tensor K is heterogeneous and isotropic:

K = α 0 0 α , with α = 1 on Ω 1 ∪ Ω 3 , 10 -2 on Ω 2 .
We consider Problem [1] with Dirichlet boundary conditions: Γ D = ∂Ω, Γ N = ∅, and ū = u| ∂Ω , and discretize the problem using the mesh mesh6.

The exact solution of this problem being piecewise affine, most of the schemes are exact on this test; the non-exact schemes are: CVFE, DDFV-HER, DDFV-OMN, MFV, NMFV. A few schemes (namely some DDFV and DG schemes) have values which are less than the exact minimum (see Table 11).

Test 8: perturbed parallelograms [?]

This test case was given to us by I. Aavatsmark [?], and is meant to test the schemes for the violation of the maximum principle within the domain. The domain Ω is a parallelogram, which is represented in Figure 11. The parameters shown in Figure 11 are X = 1, Y = 1/30 and θ = 30 Note that the solution u of this problem should be a function with a maximum in cell (6, 6), decreasing smoothly to zero towards the boundary. If u shows internal oscillations or if u < 0, Hopf's first lemma is violated. In the first line of Tables 12 and13, we give the minimum and maximum values and the fluxes obtained with the scheme SUSHI-P for a 201 × 201 uniform grid which was chosen parallel to the axes. The behaviour of the corresponding approximate solution is represented in Figure 12. The minimum value umin given by the schemes ranges from -7.21E-02 (FVSYM) to 3.26E -06 (SUSHI-P), while the maximum value umax ranges from 6.77E -03 (SUSHI-P) to 3.99E -01 (FVPMMD).

There are only three schemes which remain positive, namely NMFV, FVPMMD and SUSHI-P. Both latter schemes, however, provide fluxes which differ from the expected fluxes (which are calculated from a fine grid discretization and given on the first line of Table 13). The only scheme which respects positive while also yielding the correct flux is NMFV, which is a nonlinear scheme. This reinforces the conjecture [?] that nonlinear schemes should be used if positiveness is sought, even in the case of linear problems,.

Test 9: anisotropy and wells [?]

Here, Ω is again the square unit domain Ω = (0, 1) × (0, 1). The medium is homogeneous and anisotropic with

K = M (-θ) 1 0 0 10 -3 M (θ), M (θ) = cos θ sin θ -sin θ cos θ , [5]
where θ = 67.5 • .

The source density f is zero in all cells. Problem [1] is considered with homogeneous Neumann boundary conditions (Γ N = ∂Ω) on a square uniform grid mesh9 with 11 × 11 cells. The pressure is fixed in two cells, approximating a sink and a source with fixed pressure:

u = 0 in cell (4, 6), u = 1 in cell (8, 6).
[6]

The solution u of this problem should satisfy u ∈ [0, 1]. This discrete problem could arise from a coarse discretization of a flow field with two wells, and where the pressure is given at the well, each well being represented by one grid cell. Of course, if the mesh were refined, then all the new fine cells that are included in the coarse mesh cells (4, 6) and (8, 6) should be set to 0 and 1 respectively .

However, the point here is not the convergence of the method as the mesh tends to 0, but rather its behaviour with respect to the discrete maximum principle on a coarse grid. In particular, we would like the scheme to give an approximate solution that: 1) stays between 0 and 1, 2) does not oscillate within the domain, 3) has no extremum on the outer (no-flow) boundary. On this test, most of the schemes which were tested violate the Hopf lemmas, as can be seen in Table 14. The SUSHI-NPscheme, although it does not show internal oscillations, has extreme values which are very much out of the bounds (-1.00 and 2.00) while the extremal values of most of the other schemes violate their bounds from less than 10 per cent. Most of the schemes for which the values are provided show oscillations of the solution. The only schemes which remain positive are the schemes FVPMMD, SUSHI-P and NMFV. This latter scheme, however, provides values which are greater than 1 (it is positive, but does not respect the discrete maximum principle) and shows some small internal oscillations. The schemes FVPMMD and SUSHI-P provide values which stay within the bounds 0 and 1, and show very little oscillations within the domain. However, the values provided by the schemes are quite different, since the values given by FVPMMD on the boundary range from roughly 0.3 to 0.6, while those given by SUSHI-P range from 0.02 to 0.98., values which are close to those given by the MFE scheme in Table 15.

Conclusion

This paper proposes a comparison of about 20 numerical approximations for a family of two-dimensional anisotropic diffusion problems. The nine tests presented here involve both a wide class of diffusion tensors (anisotropic and at time heterogeneous and/or discontinuous) and a wide class of (sometimes non-conforming) meshes.

A first remarkable feature of the benchmark is the relative homogeneity of the results which are obtained by the different methods. In some cases, the convergence curves are so close to one another that it is difficult to distinguish them. It is therefore difficult to rate the schemes, and we did not attempt to do so. The choice of a scheme is problem dependent. It may result from considerations on robustness (postiveness, maximum principle, etc.) and accuracy, but can also be largely influenced from other considerations such as existing codes, coupling problems or merely data structures.

As predicted by the theory, we observe that the higher order schemes FEP2, FEQ2 or DG-C feature a better order of convergence for the L 2 norm of the solution and also asymptotically for the L 2 norm of its gradient. Nevertheless, DDFV methods are able to provide particularly low errors of the gradient, which are even smaller than the higher order schemes on coarse meshes.

The robustness of the schemes with respect to the maximum principle is a great challenge. The results vary from one scheme to another, according to the meshes or to the anisotropy of the diffusion tensor. Cell centred schemes and the mimetic scheme generally seem a bit more robust than the DG or DDFV schemes. However, in more severe cases, it seems, as conjectured in [?], that nonlinear schemes like FVPMMD or NMFV should be used even for linear problems to ensure the positiveness of the schemes.

We expect to enrich this benchmark with other existing methods and with some other challenging problems in view to improve the existing methods or develop some new algorithms. We also hope that this benchmark will serve as a bridge between the different scheme developers: it has proven that the behaviour of the schemes is comparable on the cases which are studied here, and it would be interesting to compare them on a more mathematical basis. We finally plan to extend this benchmark to some coupling between such an anisotropic elliptic equation and a transport equation and why not to 3D anisotropic diffusion problems.

Acknowledgment

The present paper is a synthesis of the work of many people:

• we first warmly thank all participants for their active collaboration in this benchmark, and for letting us use their results in this synthesis,

• we are also very grateful to I. Aavatsmark, F. Boyer, R. Eymard, T. Gallouët, C. Le Potier, G. Manzini and R. Masson for their help in the benchmark design. This work was supported by GDR MOMAS, CNRS/PACEN.

The meshes

We give here the figures of the families of meshes which were used here, except for the uniform square meshes. The family of meshes mesh2 is a uniform square mesh, with size ranging from 3.54E-01 (mesh2_2) to 5.52E-03 (mesh2_7). The mesh mesh9 is a 11×11 uniform square mesh. For more details see http://www.latp. univ-mrs.fr/fvca5/, where the data files are also given. [DIP 08] DI PIETRO D., ERN A., "Benchmark for anisotropic problems. A discontinuous Galerkin flux for anisotropic heteroge-neous second-order elliptic problems", these proceedings, 2008.

[DUB 08] DUBOIS F., LALLEMAND P., TEKITEK M. M., "Benchmark for anisotropic problems. Using Lattice Boltzmann scheme for anisotropic diffusion problems", these proceedings, 2008.

[EYM 08] EYMARD R., GALLOUËT T., HERBIN R., "Benchmark for anisotropic problems. SUSHI: A Scheme Using Stabilization and Hybrid Interfaces for anisotropic heterogeneous diffusion problems", these proceedings, 2008.

Figure 1 .Figure 2 .

 12 Figure 1. L 2 norm of the solution error for test 1.1, versus nunknw (upper left), nnmat (upper right) and h (lower left) for the triangular meshes mesh1_i

  [1] with non-homogeneous Neumann boundary conditions (Γ D = ∅ and Γ N = ∂Ω), so that the right hand side is f = -∇ • (K∇u) and the Neumann boundary condition g = (K∇u • n)| ∂Ω . Furthermore, we shall ensure the uniqueness of the solution to [1] by enforcing the condition Ω u dx = 0. The meshes are the triangular meshes mesh1.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 3456 Figure 3. L 2 norm of the solution error for test 1.2, versus nunknw (upper left), nnmat (upper right) and h (lower left) for the triangular mesh1_i

Figure 7 .

 7 Figure 7. Approximate solution on a fine grid for Test 3, oblique flow

  1 = (0.; .5] × 4 k=0 [.05 + 2k × .1; .05 + (2k + 1) × .1) , Ω r 1 = (.5; 1) × 4 k=0 [2k × .1; (2k + 1) × .1) .

Figure 9 .Figure 10 .

 910 Figure 9. L 2 norm of the solution error for test 5, versus nunknw (upper left), nnmat (upper right) and h (lower left), for the rectangular mesh mesh2

Figure 11 .

 11 Figure 11. Parallelogram-shaped domain Ω showing the distances X and Y and the angle θ

Figure 12 .

 12 Figure 12. Approximate solution on a fine uniform grid for Test 8

Figure 13 .Figure 14 .Figure 15 .

 131415 Figure 13. Triangular mesh with acute angles: meshes mesh1_1 (left) and mesh1_4 (right) -Mesh size ranging from 2.50E-01 (mesh1_1) to 3.91E-03 (mesh1_7)

Figure 16 .Figure 17 .

 1617 Figure 16. Non-conforming regular rectangular mesh mesh5 -Mesh size 1.41E-01

Figure 18 .

 18 Figure 18. Perturbed parallelogram grid mesh8 with 11 × 11 cells -Mesh size 1.24E-01. To visualize the grid, on this picture the height is 3 times the real height of the grid described in the test case

Table 1 .

 1 Test

	DG-C	-2.33E-03 9.96E-01	-3.24E-04 9.99E-01
	DG-W	-7.90E-05 9.22E-01	-8.18E-06 9.66E-01
	FEQ1	0.00E+00 8.61E-01	0.00E+00 9.37E-01
	FEQ2	0.00E+00 9.99E-01	0.00E+00 1.00E+00
	FVHYB	2.14E-03 9.84E-01	7.16E-04 9.93E-01
	FVSYM	7.34E-03 9.59E-01	2.33E-03 9.89E-01
	MFD-BLS	8.54E-03 9.55E-01	2.44E-03 9.87E-01
	MFD-FHE	9.73E-03 9.45E-01	2.90E-03 9.83E-01
	MFD-MAN	6.64E-03 9.71E-01	1.50E-03 9.93E-01
	MFD-MAR	8.82E-03 9.60E-01	2.47E-03 9.88E-01
	MFV	1.08E-02 9.42E-01	3.34E-03 9.82E-01
	NMFV	1.30E-02 1.11E+00	3.61E-03 1.04E+00
	SUSHI-NP	7.64E-03 8.88E-01	2.33E-03 9.61E-01

1.1: Minimum and maximum of the approximate solutions for mesh 4_1 and mesh 4_2

Table 2 .

 2 Test

		, 1	1.04E+00, 3	1.09	/	6.36E+02
	CVFE	-1.01E-00, 2	1.01E-01, 2	2.00	1.00	1.55E+01
	DDFV-BHU	-9.27E-01, 1 1.17E+00, 1	1.76	1.21	7.51E+00
	DDFV-HER	-4.20E-01, 2 9.12E+00, 4	/	/	7.16E-03
	DDFV-OMN	-8.24E-01, 1	7.76E-01, 1	2.00	1.00	2.11E+00
	DG-W	-1.18E-01, 1	1.18E-01, 1	2.00	1.00	1.69E+01
	FEP1	-9.48E-03, 1	9.75E-03, 1	2.00	1.01	/
	FEP2	-9.56E-01, 1	9.56E-01, 1	2.97	2.00	/
	FVSYM	-1.76E+00, 2 1.80E+00, 2	2.38	1.47	7.29E+00
	MFD-BLS	-6.50E+00, 2 5.75E+00, 2	2.54	/	3.59E+01
	MFD-FHE	-6.50E+00, 2 5.75E+00, 2	2.54	1.51	3.59E+01
	MFD-MAN	-6.62E+00, 2 5.50E+00, 2	2.49	1.50	3.58E+01
	MFD-MAR	-6.50E+00, 2 5.75E+00, 2	2.53	/	3.59E+01
	MFE	-6.50E+00, 2 5.75E+00, 2	2.53	1.47	3.59E+01
	MFV	-6.50E+00, 2 5.75E+00, 2	2.41	1.51	3.58E+01
	SUSHI-P	-6.50E+00, 2 5.75E+00, 2	2.53	1.47	3.59E+01
	SUSHI-NP	-1.93E-02, 4	1.89E-02, 4	0.37	1.99	/

2, (numerical locking) δ = 10 5

Table 3 .

 3 Test

	-3 .

2 (numerical locking) δ = 10 6

Table 4 .

 4 The computational domain and approximate solution on a fine grid (320 × 320) for Test 3, the values of umin, umax and the energies

	Test 4, vertical fault

Table 5 .

 5 Test 3, the fluxes

Table 6 .

 6 Test

		umin	umax	ener1	eren
	CMPFA	4.55E-02	9.56E-01	/	/
	CVFE	0.00E+00 1.00E+00	45.9	1.04E-02
	DDFV-BHU	4.04E-02	9.61E-01	42.1	3.65E-02
	DDFV-HER	0.00E+00 1.00E+00	49.3	1.75E-01
	DDFV-OMN	4.04E-02	9.62E-01	42.2	3.65E-02
	DG-W	-3.34E-01 1.33E+00	43.5	1.38E-02
	FVHYB	4.52E-02	9.59E-01	41.4	6.12E-02
	MFD-BLS	3.31E-02	9.71E-01	33.9	7.93E-14
	MFD-MAN	2.84E-02	9.75E-01	31.4	1.16E-12
	MFD-MAR	4.03E-02	9.60E-01	41.1	1.30E-13
	MFV	4.93E-02	9.54E-01	49.9	4.21E-05
	NMFV	4.33E-02	9.58E-01	/	/
	SUSHI-NP	1.32E-03	9.99E-01	39.1	6.67E-02

4, non-conforming rectangular mesh: mesh5. The values of umin, umax and the energies

Table 7 .

 7 Test 4, non-conforming rectangular mesh: mesh5. The fluxes

	DDFV-OMN	-40.0 41.8 -1.81 9.08E-04
	DG-W		-43.1 45.3 -2.19 1.50E-03
	FVHYB		-44.3 46.3	0.49	1.55E-04
	MFD-BLS		-32.3 36.2 -3.94 1.22E-03
	MFD-MAN		-29.7 34.1 -4.37 1.01E-03
	MFD-MAR		-39.8 42.5 -2.68 9.95E-04
	MFV		-44.0 50.3 -8.03 1.72E+00
	NMFV		-43.2 44.5 -1.23 2.32E-04
	SUSHI-NP		-40.9 43.1 -2.21 6.94E-04
			umin	umax		ener1	eren
	CMPFA	1.32E-03	9.99E-01	/	/
	CVFE	0.00E+00 1.00E+00	43.3	6.25E-04
	DDFV-BHU	1.31E-03	9.98E-01	43.2	1.27E-03
	DDFV-HER	0.00E+00 1.00E+00	43.8	1.64E-02
	DDFV-MNI		00E+00	1.00E+00	43.8	6.23E-02
	DDFV-OMN	1.32E-03	9.99E-01	43.2	1.28E-03
	DG-C	-2.20E-06 1.00E+00	43.2	1.46E-04
	DG-W	-2.21E-09 1.00E+00	43.2	7.63E-04
	FEQ1	0.00E+00 1.00E+00	43.3	2.31E-03
	FEQ2	0.00E+00 1.00E+00	43.2	0.00E+00
	MFD-BLS	1.32E-03	9.99E-01	43.2	2.84E-12
	MFD-FHE	2.12E-02	1.00E+00	43.2	3.53E-04
	MFD-MAN	1.32E-03	9.99E-01	43.2	4.71E-14
	MFD-MAR	1.32E-03	9.99E-01	43.2	2.69E-12
	MFE	1.32E-03	9.99E-01	43.2	4.20E-04
	MFV	1.32E-03	9.99E-01	43.2	1.88E-05
	NMFV	1.32E-03	9.99E-01	43.2	5.92E-04
	SUSHI-NP	1.32E-03	9.99E-01	43.1	8.88E-04

Table 8 .

 8 Test 4, the square 320 × 320 mesh mesh5 ref . The values of umin, umax and the

		flux0 flux1 fluy0	fluy1
	CMPFA	-42.1 44.4 -2.33 7.97E-04
	CVFE	-42.2 44.5 -2.25 7.98E-04
	DDFV-BHU	-42.1 44.4 -2.33 7.94E-04
	DDFV-HER	-42.0 44.3 -2.35 7.97E-04
	DDFV-MNI	-39.9 42.6 -2.68 8.01E-04
	DDFV-OMN	-42.1 44.4 -2.33 7.97E-04
	DG-C	-42.1 44.5 -2.32 7.98E-04
	DG-W	-42.1 44.5 -2.32 7.96E-04
	FEQ1	-42.2 44.5 -2.16 7.96E-04
	FEQ2	-42.1 44.5 -2.32 7.96E-04
	MFD-BLS	-42.1 44.4 -2.33 7.96E-04
	MFD-FHE	-42.1 44.5 -2.47 7.98E-04
	MFD-MAN	-42.1 44.4 -2.33 7.97E-04
	MFD-MAR	-42.1 44.4 -2.33 7.96E-04
	MFE	-42.1 44.4 -2.33 -7.96E-04
	MFV	-42.1 44.4 -2.33 8.33E-04
	NMFV	-42.1 44.4 -2.33 7.97E-04
	SUSHI-NP	-42.1 44.4 -2.33 7.97E-04

energies the fluxes), while respecting the discrete maximum principle. Note that the MFD implementations MFD-BLS and MFD-MAN give surprisingly low values of both energy and fluxes in comparison with MFD-MAR.

Table 9 .

 9 Test 4, the square 320 × 320 mesh mesh5 ref . The fluxes

Table 10 .

 10 Test 5: minimum and maximum values on the coarsest mesh for the schemes which do not satisfy the maximum principle

Table 11 .

 11 Test 7, the conforming mesh: mesh6 ( = machine precision)

		erl2	ergrad	umin	umax
	CMPFA	1.23E-15	/	-5.54 5.37E-01
	CVFE	1.83E-05 3.86E-05 -5.57 5.75E-01
	DDFV-BHU	1.38E-14 5.86E-14 -5.575 5.75E-01
	DDFV-HER	6.53E-08	/	-5.575 5.75E-01
	DDFV-MNI	1.38E-15 3.48E-14 -5.58 5.75E-01
	DDFV-OMN	3.79E-08 4.51E-08 -5.58 5.75E-01
	DG-C	1.39E-13 1.88E-12 -5.58 5.75E-01
	DG-W	1.18E-14 1.89E-14 -5.58 5.75E-01
	FEQ1	3.79E-15 6.49E-14 -5.575 5.75E-01
	FEQ2	6.12E-15 8.19E-14 -5.575 5.75E-01
	FVHYB	1.24E-15 5.10E-14 -5.54 5.37E-01
	FVSYM			-5.53 5.37E-01
	MFD-BLS		/	-5.54 5.37E-01
	MFD-FHE	3.25E-15 4.73E-15 -5.54 5.37E-01
	MFD-MAN	1.39E-15 1.92E-15 -5.54 5.37E-01
	MFD-MAR	8.19E-13	/	-5.54 5.37E-01
	MFV	1.40E-08 4.69E-08 -5.54 5.37E-01
	NMFV	4.98E-03	/	-5.54 5.37E-01
	SUSHI-NP	1.30E-15 1.35E-14 -5.54 5.37E-01

• .

Table 12 .

 12 Minimum and maximum values for Test 8

		umin	umax	umin	umax
	Fine grid	1.07E-24 4.10E-01
	CMPFA	-2.31E-02 1.03E-01	FVHYB	-3.38E-02 1.12E-01
	CVFE	-1.23E-03 4.24E-02	FVSYM	-7.21E-02 1.52E-01
	DDFV-BHU	-1.25E-03 8.22E-02	FVPMMD	1.22E-09 3.99E-01
	DDFV-HER	-1.61E-03 8.99E-02	MFD-BLS	-1.03E-01 1.85E-01
	DDFV-MNI	-1.46E-03 6.69E-02	MFD-FHE	-6.54E-02 1.44E-01
	DDFV-OMN	-1.77E-03 8.36E-02	MFD-MAR	-2.62E-02 9.07E-02
	DG-C	-7.33E-03 1.05E-01	MFV	-8.08E-03 5.81E-02
	DG-W	-9.03E-03 6.57E-02	NMFV	3.05E-15 9.42E-02
	FEQ1	-4.17E-03 4.90E-02	SUSHI-NP	-1.19E-03 5.65E-02
	FEQ2	-5.07E-03 8.04E-02	SUSHI-P	3.26E-06 6.77E-03
				flux0	flux1	fluy0	fluy1
	Fine grid		5.46E-21	5.46E-21 5.00E-01 5.00E-01
	CVFE		-1.17E-05	2.63E-05 2.87E-01 5.54E-01
	DDFV-BHU	-5.814E-10 -3.35E-10 4.97E-01 5.02E-01
	DDFV-HER		2.45E-10	-1.83E-10 4.80E-01 5.11E-01
	DDFV-MNI	-4.37E-05	6.50E-05 5.07E-01 4.93E-01
	DDFV-OMN	-6.86E-10 -4.86E-10 4.98E-01 5.02E-01
	DG-C		-1.66E-07	1.41E-07 5.08E-01 4.92E-01
	DG-W			5.02E-01	0.00E+00 4.98E-01 0.00E+00
	FEQ1			5.51E-06	7.15E-05 5.46E-01 4.89E-01
	FEQ2		-5.52E-05	1.96E-05 4.98E-01 5.01E-01
	FVSYM			1.37E-04	-1.15E-04 4.96E-01 5.04E-01
	FVPMMD		1.76E-06	3.5E-06	4.55E-01 5.44E-01
	MFD-BLS	-5.14E-04 -3.13E-03 5.01E-01 5.03E-01
	MFD-FHE		4.48E-04	-4.08E-03 5.03E-01 5.01E-01
	MFD-MAR	-1.79E-02	2.61E-03 5.05E-01 5.10E-01
	MFV		-2.30E-02	4.95E-02 2.74E-01 6.99E-01
	NMFV			0.00E+00	0.00E+00 4.99E-01 5.01E-01
	SUSHI-NP		7.35E-04	1.29E-04 4.99E-01 5.00E-01
	SUSHI-P	-4.21E-02 -3.29E-02 5.38E-01 5.37E-01

Table 13 .

 13 Boundary fluxes for Test 8

Table 14 .

 14 Values of umin and umax for Test 9

	0.03 0.03 0.04 0.03 0.02 0.08 0.25 0.48 0.73 0.90 0.97
	0.03 0.03 0.03 0.02 0.03 0.14 0.34 0.59 0.81 0.94 0.97
	0.03 0.03 0.04 0.02 0.06 0.21 0.44 0.69 0.88 0.96 0.97
	0.03 0.03 0.04 0.02 0.10 0.30 0.53 0.78 0.94 0.97 0.97
	0.03 0.03 0.05 0.02 0.17 0.41 0.62 0.89 0.96 0.96 0.97
	0.03 0.04 0.04 0.00 0.28 0.50 0.72 1.00 0.96 0.96 0.97
	0.03 0.04 0.04 0.11 0.38 0.59 0.83 0.98 0.95 0.97 0.97
	0.03 0.03 0.06 0.22 0.47 0.70 0.90 0.98 0.96 0.97 0.97
	0.03 0.04 0.12 0.31 0.56 0.79 0.94 0.98 0.96 0.97 0.97
	0.03 0.06 0.19 0.41 0.66 0.86 0.97 0.98 0.97 0.97 0.97
	0.03 0.10 0.27 0.52 0.75 0.92 0.98 0.97 0.96 0.97 0.97

Table 15 .

 15 Values given by the MFE scheme for Test 9


L Rakotondrainibe

Y Kokar

G Zaharia

G El Zein

E Tanguy

H Li

B Charbonnier

Hybrid, Optical and Wireless Near-Gigabit Communications System

This paper presents the study and the realization of a hybrid 60 GHz wireless communications system. As the 60 GHz radio link operates only in a single-room configuration, an additional Radio over Fibre (RoF) link is used to ensure the communications in all the rooms of a residential environment. A single carrier architecture is adopted. The system uses low complexity baseband processing modules. A byte/frame synchronization technique is designed to provide a high value of the preamble detection probability and a very small value of the false alarm probability. Conventional RS (255, 239) encoder and decoder are used to correct errors in the transmission channel. Results of Bit Error Rate (BER) measurements are presented for various antennas configurations.

I. INTRODUCTION

The 60 GHz band, due to a large bandwidth is one of the most promising solutions to achieve a gigabit class for short distance high speed communications. Aspects including standardization, justification of using 60 GHz frequency, 60 GHz propagation, antennas and key system design issues were addressed in [START_REF] Smulders | On the design of Low-Cost 60 GHz Radios for Multigigabit-per-Second Transmission over Short Distances[END_REF][START_REF]ECMA International High Rate Short Range Wireless Communications[END_REF][START_REF] Young | An Overview of Multigigabit Wireless through Millimeter Wave Technology: Potentials and Technical Challenges[END_REF][START_REF] Funada | A design of single carrier based PHY for IEEE 802.15.3c standard[END_REF]. The characteristics of the 60 GHz band affect the system design. Problems such as power amplifier (PA) non-linearity and oscillator phase noise are more important for the circuits design. These effects should be taken into account in the overall communication system.

This paper proposes a hybrid optical/wireless system derived from simplified IEEE802.15.3c PHY layer proposal to ensure near 1 Gbps data rate on the air interface. The first system application in a point-to-point configuration is the highspeed file transfer. The system must operate in indoor, line-ofsight (LOS) domestic environments.

In sections II and III, the general design of the transmitter and the receiver are respectively described. In these sections, the intermediate frequency (IF) and radiofrequency (RF) blocks are first presented. Then, the baseband (BB) blocks are described. The byte/frame synchronization method is presented in baseband section. Measurement results are presented in section IV. Section V concludes the work. II. TRANSMITTER DESIGN Fig. 1 shows the block diagram of the transmitter (Tx). The system uses a single carrier architecture based on a Differential encoded Binary Phase Shift Keying (DBPSK) modulation. A differential encoder is used before the BPSK modulator to remove the phase ambiguity at the receiver (Rx), knowing that a differential demodulation is used.

A. IF-Tx and RF-Tx architecture

After channel coding and scrambling, the input data are differentially encoded using logic circuits (PECL). The differential encoder performs the delayed modulo two addition of the input data with the output. The obtained data are used to modulate an IF carrier generated by a 3.5 GHz phase locked oscillator (PLO) with a 70 MHz external reference. The IF modulated signal is fed into a band-pass filter (BPF) with 2 GHz bandwidth, and transmitted through a 300 meters fibre link. This IF signal is used to modulate directly the current of the Vertical Cavity Surface Emitting Laser (VCSEL) operating at 850 nm through a band-pass filter with a bandwidth of 2 GHz. The VCSEL input RF power must not exceed -3 dBm to avoid signal distortions. After transmission, the optical signal is converted to an electrical signal by a PIN diode and amplified. The overall RoF link has 0 dB gain and an 8 GHz bandwidth. This bandwidth could be increased if necessary by using a VCSEL and a photodetector of broader bandwidth.

Following the RoF link, the IF signal is sent to the RF block. This block is composed of a mixer, a frequency tripler, a PLO at 18.83 GHz and a band-pass filter (59-61 GHz). The local oscillator frequency is obtained with an 18.83 GHz PLO with the same 70 MHz reference and a frequency tripler. The phase noise of the 18.83 GHz PLO signal is about -110 dBc/Hz at 10 kHz off-carrier. The BPF prevents spill-over into adjacent channels and removes out-of-band spurious signals caused by the modulator operation. The 0 dBm obtained signal is fed into the horn antenna with a gain of 22.4 dBi and a half-power beamwidth (HPBW) of 12°.

B. BB-Tx architecture

The Gigabit-Ethernet (G.E.) interface is used to connect a home server to a wireless link with around 800 Mbps bit rate, as shown in Fig. 2. The Gigabit Media Independent Interface (GMII) is an interface between the Media Access Control (MAC) device and the physical layer (PHY). The interface defines speeds up to 1 Gbps, implemented using 8 bits data interface clocked synchronized at 125 MHz. However, this frequency is different from the clock sampling data (104.54 MHz) generated by a clock manager, as explained latter. These clocks asynchronism provide data packet loss. In order to avoid jitters transmission, a programmable circuit (FPGA) is used as part of buffers memory. The byte stream from the G.E. interface is transfered into the dual port FIFO memory (due to the RS encoder), as shown in Fig. 3. Otherwise, the data input must be paralellized by a serial-to-parallel (S/P) converter. The frame format consists of 4 preamble bytes, 239 data burst bytes, 16 check bytes and an additional dummy byte, as shown in Fig. 4. A known preamble is sent at the beginning of each frame to achieve the frame synchronization. The used preamble is a Pseudo-Noise (PN) sequence of 31 bits + 1 bit to provide 4 preamble bytes. A dummy byte should be added in the frame structure in order to make a code word of 256 bytes (a multiple of 4) useful for the scrambling operation. Hence, owing to the frame structure, two different clock frequencies f 1 and f 2 are used in the baseband blocks: The structure of the frame is obtained as follows: the input byte stream is written into the dual port FIFO memory at frequency f 1 . The FIFO memory has been set up to use two different frequencies for writing at f 1 and reading at f 2 . Therefore, reading can be started when the FIFO memory is half-full (to avoid overflow/underflow). The encoding control generates 4 preamble bytes and reads 239 data burst bytes (number of bytes used in the frame to be coded). Consequently, the minimum of the maximum crosscorrelation is obtained for several values of k. For each value of k, the cross-correlation between P and T is analyzed (for 1 ≤ i ≤ 8) in order to identify the best dummy byte which gives the lowest cross-correlation values. The best result is obtained for k = 64 which gives d = [0 0 0 0 0 0 1 0]. As shown in Fig. 6, the maximum cross-correlation is equal to 16. This value is obtained twice; all the other values are smaller. This optimal dummy byte gives a smaller value of the preamble false detection probability. The scrambler is a PN-sequence of 4 bytes (31 bits + 1 additional bit). This scrambler is chosen in order to provide the lowest cross-correlation values between the received scrambled data and the PHY preamble. This choice decreases the false alarm detection probability of the PHY preamble in the scrambled data. The obtained byte stream is finally serialized just before the differential encoder. III. RECEIVER DESIGN Fig. 7 shows the receiver system block diagram. To simplify the receiver architecture, a non-coherent demodulation is used.

A. IF-Rx and RF-Rx architecture

The input band-pass filter removes the out-of-band noise. The RF signal at the output of the filter is down-converted to an IF signal centered at 3.5 GHz and fed into a band-pass filter with a bandwidth of 2 GHz. An Automatic Gain Control (AGC) with a dynamic range of 20 dB is used to ensure a quasi-constant signal level at the demodulator input. A Low Noise Amplifier (LNA) with a gain of 40 dB is used to achieve sufficient gain. A simple differential demodulation enables the coded signal to be demodulated and decoded. Compared to a coherent demodulation, this method is less performing in additive white Gaussian noise (AWGN) channel. The significant impact on the system caused by the radio channel is the frequency selectivity which induces intersymbol interference (ISI). Nevertheless, the differential demodulation is more resistant to ISI effect. In fact, the differential demodulation, based on a mixer and a delay line (delay equal to the symbol duration Ts = 1.14 ns), compares the signal phase of two consecutive symbols. Due to the product of two consecutive symbols, the rate between the main lobe and the second lobe of the impulse response of the channel is increased. Following the loop, a Low-Pass Filter (LPF) with 1 GHz cut-off frequency removes the highfrequency components of the obtained signal. For better acquisition of the Clock and Data Recovery (CDR) circuit, long sequences of '0' or '1' must be avoided. Therefore, the use of a scrambler and a descrambler are compulsory.

B. BB-Rx architecture

Fig. 8 shows the block diagram of the BB-Rx architecture. The synchronized data after the CDR are converted into a byte stream (due to the RS decoder). The preamble detection is based on the cross-correlation between 32 successive received bits and the PHY preamble of 32 bits. Further, each correlator must analyze a 1-bit shifted sequence. Therefore, the preamble detection is performed with 32 + 7 = 39 bits (+ 7 because of different possible shifts of a byte). In all, there are 8 corrrelators of 32 bits in each correlators-bank, as shown in Fig. 9.

In order to minimize the preamble miss detection probability, two banks of correlators are used (taking into account the periodical repetition of the PHY preamble). The decision is made from 264 successive bytes (P1 + 256 data bytes + P2) stored in a register. Each value of the correlation between 32 successive received bits and the PHY preamble is compared to the threshold. The preamble is detected when the same correlator C k of each correlators-bank validated its presence.

The byte/frame synchronisation performance is characterized by the miss detection probability P m , the false alarm probability P f and the channel error probability p [START_REF] Urkowitz | Signal Theory and Random Processes[END_REF]. Fig. 10 shows the miss detection probability of the PHY preamble as a function of channel error probability p. The threshold  = 32 is equal to the maximum cross-correlation value obtained when the received preamble has no error. The byte/frame synchronisation also depends on the false alarm events. Occasionally, due to the random values of the binary transmitted data, the architecture shown in Fig. 9 can falsely declare that the preamble is detected. A false alarm is declared when 2 correlators C k of the two correlator-banks indicate the detection of the preamble within the data burst bytes. Fig. 11 shows the false alarm probability P F1 (with one correlators-bank) and P F2 (with two correlators-banks) as functions of p. We can see that for  = 28, P F2 = 10 -10 and P F1 = 10 -5 . Therefore, by means of two correlators-banks, the best tradeoff  = 28 gives the high detection probability of the preamble and a very small false alarm probability. The effect of p on the false alarm probability is insignificant since the random data bits "0" and "1" are considered equiprobable.

After the byte alignment and the preamble detection, the descrambler performs the modulo 2 addition between 4 successive received data and the scrambler sequence. The RS decoder processes the descrambled bytes and attempts to correct the eventual errors. The RS (255, 239) decoder can correct up to 8 erroneous bytes and operates at a high clock frequency (109.37 MHz). The byte stream obtained at the decoder output is written discontinuously into the dual port FIFO memory at the clock frequency f 2 . A clock frequency f 1 = 100.54 MHz read out continuously the data bytes stored in the FIFO memory. Finally, the byte stream is transmitted to the receiver Gigabit Ethernet interface, as shown in Fig. 12. IV. MEASUREMENT RESULTS An HP 8753D vector network analyzer was used to determine the frequency response and the impulse response of RF blocks (Tx + Rx) including the LOS propagation channel. The objective was to measure the system bandwidth and to estimate the multipath channel effects. Measurements were located in a corridor where the major part of the transmitted power is focused to the receiver. The RF-Tx and RF-Rx were placed at a height of 1.5 m. After measurement set-up and calibration, 2 GHz available bandwidth was measured, as shown in Fig. 13 A perfect system must have an impulse response with only one lobe. Fig. 14 presents the impulse response of the RF-Tx and RF-Rx blocks placed at 10 m distance. Few side lobes were obtained which are mainly due to RF components imperfections. A back-to-back test was realized by

P F2 P F1
exchanging the Tx and Rx antennas with a 45 dB attenuator but similar results were obtained. Moreover, in order to evaluate the transmission quality, an HP70841B pattern generator was used to generate random data at the transmitter and an HP 708842B error detector was used at the receiver. Fig. 15 indicates a high transmission quality of the received data (good eye opening at 875 Mbps rate). As shown by this eye pattern, the recovered data can be obtained by sampling at half-period in order to be noise resistant. Fig. 16 shows the BER measurements at 875 Mbps as a function of the Tx-Rx distance. In our experiments, four antennas were used: two horn antennas and two patch antennas. Each patch antenna has an 8 dBi gain with a 30° at -3 dB beamwidth. The use of high-directivity horn antennas gives a remarkable BER performance and greatly reduces the harmful effects of multipath propagation. This leads to the saving of link budget as high as 28.8 dB compared to patch antennas. However, between (1) and ( 2), the link budget ratio is around of 15 dB. The system with high-gain directive antennas is acceptable for point-to-point (PP) communications links, with minimal multi-path interference. But, the 60 GHz radio links are sensitive to shadowing due to high attenuation of non-line-of sight (NLOS) propagation. In addition, the Tx and Rx antenna have to be aligned, otherwise the beam-pointing errors will cause an significant degradation of the channel quality. For properly aligned antennas, if the direct path is blocked by moving objects, the communication can be completely lost. In order to improve the system reliability, the choice of a centralized transmitter, preferably located on the centre of the room ceiling must be considered. The height of the transmitter (with less directional antenna) can partly reduce the influence of the shadowing. The receiver could be equipped with more directional antenna, pointed toward the transmitter. If antennas beamwidth is large, equalization should be adopted to overcome multi-path interference while maintaining a high data rate. Futures works will provide BER (with channel coding) measurement in different environments.

V. CONCLUSION This paper presents the design and the implementation of a 60 GHz system for high data rate wireless communications. The proposed system provides a good trade-off between performance and complexity. An original method used for the byte synchronization is also described. This method allows a high preamble detection probability and a very small false detection probability. For 1 Gbps reliable communications within a large room, antennas must have a relatively high gain.

Increasing the data rate is achieved using higher order modulations such as QPSK. Equalization methods are still under study. The demonstrator will be further enhanced to prove the feasibility of wireless communications at data rates of several Gbps in different configurations, especially in non line-of-sight (NLOS) environments.

Fig. 1 .

 1 Fig. 1. 60 GHz wireless Gigabit Ethernet transmitter

Fig. 2 .

 2 Fig. 2. Transmitter Gigabit Ethernet interface

Fig. 3 .

 3 Fig. 3. BB-Tx architecture

Fig. 4 .

 4 Fig. 4. Frame structure

Fig. 5 .

 5 Fig. 5. Maximum cross-correlation between P and T for 0 ≤ k ≤ 255

Fig. 6 .

 6 Fig. 6. Cross-correlation between P and T, for k = 64

Fig. 7 .

 7 Fig. 7. 60 GHz Wireless Gigabit Ethernet receiver (Rx)

Fig. 8 .

 8 Fig. 8. BB-Rx architecture Fig. 9 shows the architecture used for preamble detection and byte alignment.

Fig. 9 .

 9 Fig. 9. Preamble detection based on periodicity and control alignment

Fig. 10 .

 10 Fig. 10. Miss detection probability of PHY Preamble

Fig. 11 .

 11 Fig. 11. False alarm probability of PHY Preamble

Fig. 12 .

 12 Fig. 12. Receiver Gigabit Ethernet interface

  .

Fig. 13 .

 13 Fig. 13. Frequency response of RF Tx-Rx (with horn antennas)

Fig. 14 .

 14 Fig. 14. Impulse response of RF blocks (using horn antennas)

Fig. 15 .

 15 Fig. 15. Eye diagram at 875 Mbps (10 m Tx-Rx distance)

ACKNOWLEDGMENT

This work is part of research project Techim@ges supported by the French "Media & Network Cluster" and the Comidom project supported by the "Région Bretagne".


Gilles Simon
email: gilles.simon@loria.fr

Immersive Image-Based Modeling of Polyhedral Scenes

Keywords: Augmented Reality, Image-Based Modeling, Construction at a Distance, 3D User Interfaces, Wearable Computers Index Terms: I.2.10 [Artificial Intelligence]: Vision and Scene Understanding-Modeling and recovery of physical attributes I.3.6 [Computer Graphics]: Methodology and Techniques-Interaction techniques

In this paper, we describe a purely image-based system that allows a user to interactively capture the 3D geometry of a polyhedral scene with the aid of its physical presence. A video camera is used as both an interaction and tracking device. The 3D user interface is intuitive to a non-expert and the mouseless control procedure makes the system particularly suitable for mobile devices such as PDAs and mobile phones. The efficiency and accuracy of the method are demonstrated on a polyhedral scene made of two house-like boxes.

INTRODUCTION

Acquiring the 3D geometry of arbitrary scenes has been a primary objective of both the computer vision and graphics communities for many decades. Applications are numerous in various domains such as construction, GIS and 3D maps, virtual tours, visual effects and AR. Existing modeling methods usually rely on two separate stages: first, some data about the scene (photographs, videos, laser measurements, ...) are acquired on-site. Then these data are treated off-line using some specific manipulations and algorithms. Unfortunately, this process can be time-consuming and tedious; moreover, there is no guarantee after the first stage that the required model is fully extractable from the acquired data and additional acquisitions are sometimes needed to supplement the missing parts.

This paper tries to bridge the gap between data acquisition and their exploitation. A pioneer work has been described in [START_REF] Piekarski | Tinmith-Metro: New Outdoor Techniques for Creating City Models with an Augmented Reality Wearable Computer[END_REF], where mobile AR users are enabled to interactively capture planar shaped objects with the aid of their physical presence. Users interact with the computer using a set of pinch gloves and the camera is tracked using an inertial sensor and a GPS. Our system differs from that work in that the only device required for both tracking and interacting is a handheld camera coupled with three keyboard keys (e.g. a mobile phone); this makes it suitable for indoor as well as outdoor use. A purely image-based interactive model building system has recently been presented in [START_REF] Bunnun | OutlinAR: an assisted interactive model building system with reduced computational effort[END_REF]. This work has several common aspects with our work, such as the use of a camera-mouse and model-based tracking. However, the modeling task is less constrained in our system in that it does not require any initial template to be placed in the scene, nor that the current model partly stays in the camera field of view while new structures are added. Moreover, defining a vertex in [START_REF] Bunnun | OutlinAR: an assisted interactive model building system with reduced computational effort[END_REF] is based on parallax motion which is not well appropriate for modeling objects far from the camera. Finally, our system can be seen as an immersive version of the widely used 3D drawing software SketchUp TM (http://sketchup.google.com). This tool combines some of the features of pencil-and-paper sketching and some of the features of CAD systems to provide a lightweight, gesture-based interface for 3D polyhedral modeling. In the latest releases of SketchUp, the user is able to align the world axes to match a photo perspective. With this done, he can create models using the photo as a direct reference; mouse strokes are converted into 3D-space using inverse ray intersections with the previously defined geometry (or the ground plane by default). In addition, a virtual camera can be manipulated in order to move around the 3D model under construction. All these principles have been taken up in our implementation, but with the important difference that we consider dynamic video images instead of static ones; moreover, the video camera is used as the interaction device instead of a mouse and virtual camera manipulations are replaced by real ones.

OVERVIEW OF THE SYSTEM

Scene modeling is done by alternating modeling interaction and camera moving sequences (Fig. 1):

Modeling mode. Pure rotations are applied to the camera (e.g. when using a head-mounted display, this amounts to look around the scene by only turning the head). Using the camera as an interaction device, the user is able both to calibrate the camera and to define the scene geometry. It must be noticed that this mode may be used alone, providing already interesting contributions to classical single view metrology, as (i) by turning the head the reachable field of view is much larger than the camera's field of view and (ii) mouse manipulations are replaced by camera manipulations, which can be of great interest in a mobile context.

Tracking mode. When at least one face of the scene has been described and the camera undergoes a general motion, a 6-DOF (degrees of freedom) camera tracking is performed, based on the available geometry [START_REF] Vigueras | Iterative multi-planar camera calibration: Improving stability using model selection[END_REF]. This allows the user to get closer to some parts of the scene or make some new faces visible before continuing modeling. In this mode, the geometry previously modeled has to stay (at least partly) visible in the camera field of view, which is not required in modeling mode. A pose recovery procedure based on SIFT features can be called upon user request or each time a tracking failure is detected by the system.

The system starts in modeling mode. Once a 3D geometry has been initialized, switches between modeling and tracking modes are done automatically by the system, depending on the motion applied to the camera (Akaike's motion model selection is used [START_REF] Vigueras | Iterative multi-planar camera calibration: Improving stability using model selection[END_REF]).

MODELING INTERACTIONS

A key idea in this work is that the user can interact with the images using the video camera itself. Indeed, let us consider that he/she wants to draw a line stroke on a video image between two physical points A and B of a real-world scene; there are two possible ways to do this: (i) the camera stays fixed and a mouse is used to move a cursor to the two respective endpoints; (ii) the camera is rotated so that the physical points apparently move to a fixed cursor, for instance at the center of the image. The second solution, which is used in our implementation, does not require mouse support. However, the position of point A has to be updated from frame two frame while point B is aimed at. Fortunately, camera rotations only induce homographic deformations of the image, that can be computed easily e.g. using keypoint matches. In both cases a click has to be done when the cursor (or the camera) is correctly positioned. In the second solution, mouse clicks can be replaced by key presses or any other input controls such as buttons, voice commands, etc. In our system, only three keys are used for all operations: one to "click" or "drag and drop" in the image, a second to cancel the current operation or request a pose recovery and a third to scroll the tools menu.

Before being able to model the scene, the user has to calibrate the camera: this is done by indicating two sets of horizontal parallel lines orthogonal to each other (Fig. 2, frame a). In addition, the user can drag and drop the origin of the world frame, providing a ray of possible 3D positions of that point. The depth of the origin is chosen so that the unit up vector has an arbitrary size in the image plane. The user can also change the scale of the scene by moving the extremity of the unit up vector. Once the camera parameters are known, the user can start modeling. New faces are instantiated in contact with existing faces or the ground plane by default. Contacts are guaranteed by "attracting" the clicked points to existing 3D points which appear in a specific color when close to the cursor: yellow for a face vertex, blue for the middle of an edge and red for the closest point on an edge. Faces under the cursor are also specifically outlined and new 3D vertices can be generated by using inverse ray intersection. In addition, most actions are guided using a line stroke, which itself can be attracted by the projected world axes.

Six different tools have been implemented which are summarized in Tab. 1. Fig. 2 illustrates how a simple house-like box can be modeled using the Add, Extrude, Cut and Move tools: a rectangular face is first created using two line strokes AB and BC: according to the priority rules described at bottom of Tab. 1, AB is instantiated in the ground plane and parallel to the green axis (frame b); then C is dragged in the plane orthogonal to AB (frame c) and finally such that BC is aligned with the blue axis (frame d); the new face is then extruded forming a box (frame e) whose top face is subdivided so that the cut line joins the middles of two opposite edges; this line is moved along the blue axis in order to form the roof (frame f). Table 1: Manipulation tools used in our prototype.

EVALUATION

A video is associated with this paper showing the system in action (http://www.loria.fr/˜gsimon/ismar09/). A Dell Precision M6300 laptop coupled with a simple Logitech webcam were used for this experiment. The system runs at video rate in standard mode and 2 to 6Hz when the recovery procedure is called. A simple polyhedral scene, made of two house-like boxes (Fig. 2), is modeled. Several modeling and tracking sequences alternate during a 6-minute working session (the video has been cut and accelerated 1.5×). Note that it can be easily distinguished when the system is running in modeling mode and when it is running in tracking mode, as the cross cursor disappears in the second case. Fig. 3 shows the errors obtained on the recovered 3D geometry. The largest relative error (12.3%, 9 mm) actually corresponds to a depth inaccuracy on d 6 ; getting closer to the related physical face before modeling it or making it appear more "fronto-parallel" would probably have led to better accuracy.

Figure 1 :

 1 Figure 1: The system alternates between modeling operations (only pure rotations of the camera are allowed) and 6-DOF camera tracking. Thick lines represent the modeled faces of the scene.

Figure 2 :

 2 Figure 2: Snapshots of the system in use.

  Add a rectangular face by clicking three points A, B,C ( * ) . The rectangle is generated in the plane ABC using A and C as diagonally opposite vertices.Extract texture and tracking features from the selected face.Extrude the selected face.Move the selected vertex, edge or face.Subdivide the selected face by joining two of its vertices.Delete the selected face. ( * ) The inverse ray intersection priorities are: A: selected vertex > selected face > ground B: selected vertex > projected world axis AB > plane(A) C: selected vertex > projected world axis BC > plane ⊥ [AB]

Figure 3 :

 3 Figure 3: Recovered/expected distances (in cm) and percentage errors obtained on the recovered geometry.

CONCLUSION

We regard this work as a proof-of-concept application. Its main flaw resides in that 6-DOF tracking is only possible when a face of the modeled scene is visible; this can make it difficult to go "behind" the modeled scene. To tackle this issue, two different approaches may be considered: the first one may consist in allowing the user to model different parts of the scene independently each other and merge these parts later using a specific tool; the second approach may be to automatically reconstruct 3D points out of the model faces during 6-DOF tracking and use these points to track the camera when no face is visible or in addition to the visible faces; this may be done within a SLAM framework.


Elvio Di Salvo

Ziad J Ajaltouni

Time Reversal Violation: Suggested Tests for Non-Leptonic Two-Body Decays

We discuss a theoretical condition for time reversal violation in weak non-leptonic two-body decays involving spinning and unstable particles in the final state. In some cases this condition looks particularly simple and suitable for tests of time reversal violation.

Introduction

People started to believe in Time Reversal Violation (TRV), since when the first CP violation was discovered, in 1964. This is a consequence of the CPT symmetry, derived under very mild assumptions and confirmed by stringent tests. However till now only one experiment has given evidence for direct TRV. We refer to the CPLEAR experiment [1], which was realized in 1998 and consisted of establishing that the K 0 → K 0 transition rate is different than the K 0 → K 0 one. Generally it is quite difficult to realize the inverse process to a weak decay and, above all, to compare the two decay rates with enough precision for discovering possible differences. An alternative way of coping with the problem is to look for theoretical consequences of Time Reversal Invariance (TRI) and to investigate whether these are contradicted by experiment. In particular, let us consider the two non-leptonic decays K → ππ and Λ → πp. Such a kind of decays involves strong Final State Interactions (FSI) between the decay products. In the cases considered, given the low mass of the parent resonances, FSI amount essentially to elastic scattering. Then, thanks to the Fermi-Watson theorem, TRI implies the following condition:

A = e 2iδ A * , ( 1
)
where A is the weak decay amplitude and δ is the strong phase-shift of elastic scattering. In other words, if one has

A = ae i(δ+φw ) , (2)
with a real and φ w = 0, π, we conclude that Time Reversal (TR) is violated. φ w is generally known as "weak" phase; for example, it could be related to the CP-violating phase which appears in the Cabibbo-Kobayashi-Maskawa (CKM) parametrization of the Yukawa coupling with three-quark generations. However also this strategy appears impervious, since we are faced with the difficulty of determining phases with sufficient precision. The main aim of our talk is to pose the question as to whether we may be luckier with decays of higher mass resonances, produced copiously at facilities like the LHC. In particular, we shall modify condition (1) for such decays and we shall suggest how to exploit data in order to test TRV.

In section 2 we define the observables which can be extracted in two-body decays with spinning decay products. In section 3 we illustrate how to modify the condition for TRV, when inelastic channels are open. Section 4 is devoted to suggesting tests for TRV. Lastly, a short conclusion is drawn in section 5.

Observables in two-body decays

We consider two-body decays involving beauty, and therefore detectable at LHCb, like Λ b → ΛJ/ψ and B → V V , V being a vector meson. We examine the observables that can be extracted from such decays. These can be defined, for example, in the so-called helicity frame. To this end we consider a canonical frame, such that the xy-plane coincides with the production plane of the parent resonance. If p is the momentum of this resonance in the canonical frame, the helicity frame is defined, in the rest frame of one of the two decay products, by three mutually orthogonal unit vectors:

e L = p |p| , e T = u × e L |u × e L | , e N = e T × e L , (3)
where u is the unit vector in the direction of the z-azis of the canonical frame. The observables -angular distribution, polarization components and polarization correlations -can be suitably defined starting from the density matrix of the decay products.

We have

I(θ, φ) = trρ, (4)
I(θ, φ)P i (θ, φ) = tr(ρs i ), (5)
I(θ, φ)P ij (θ, φ) = tr(ρs i s j ). ( 6
)
where ρ is the density matrix of the two decay products, s i the components of the spin operator, θ and φ are respectively the polar and azimuthal angle of the momentum of one of the decay products in the helicity frame and i, j run over the indices L, T, N of the unit vectors. It is worth noting that P ij = P i P j because of quantum effects.

As an example, we give the expressions of some of the above observables, relative to the decay Λ b → ΛJ/ψ [2], as functions of the rotationally invariant decay amplitudes, A λ 1 ,λ 2 , λ 1 and λ 2 being the helicities of the two decay products. We have

I(θ, φ) = 1 4π (|A 1/2,0 | 2 + |A -1/2,0 | 2 + |A 1/2,1 | 2 + |A -1/2,-1 | 2 ), (7)
I(θ, φ)P Λ L (θ, φ) = 1 4π (|A 1/2,0 | 2 -|A -1/2,0 | 2 + |A 1/2,1 | 2 -|A -1/2,-1 | 2 ), (8)
I(θ, φ)P T N (θ, φ) = 1 4π √ 2 Im(A -1/2,-1 A * 1/2,0 -A -1/2,-1 A * 1/2,0 ). ( 9
)
Really the above formulae hold in the case of unpolarized Λ b , which is generally not true, since spin-orbit coupling causes a transverse polarization in the production reaction. However, also with this simplifying assumption, it can be shown [2] that one obtains an over-determined linear system, whose unknowns are the products

A λ 1 ,λ 2 A * λ ′ 1 ,λ ′ 2
. Solving this system allows to determine all moduli of the amplitudes and all phases up to one, which we set conventionally to zero. It can be shown that this is a rather general result, valid also, e. g., for decays of the type 0 → 1 1.

A condition for TRV with inelastic channels open

Now we generalize eq. ( 1) to the case when more decay modes are present and therefore inelastic processes are admitted in the scattering between the decay products. In this case TRI implies[3]

A m = S mn A * n , (10)
where S is the scattering matrix and A m is the decay amplitude to a given eigenstate of the total angular momentum, denoted by m. The solution to this equation can be obtained by diagonalizing the S-matrix. To this end we recall a result by Suzuki [3], that is,

S mn = k O mk e 2iδ k O T kn , (11)
where O is an orthogonal matrix and the δ k are strong phase-shifts. Eqs. ( 10) and (11) yield

A m = n O mn a n e iδn , (12)
where a n are real amplitudes representing the effects of weak interactions on the decay process. Obviously complex values of one or more such amplitudes -that is, "weak" phases -would imply TRV, analogously to the case when only elastic scattering is allowed between the decay products. We remark that, in a local field theory -like the SM -, even one complex a n implies CP violation, owing to CPT symmetry. As told in the introduction, in the SM the phase of such a complex amplitude is related to the CKM parametrization.

Suggested tests for TRV

Unfortunately the result found in the previous section is generally not useful for detecting TRV. Indeed, first of all, the elements of the matrix O are not known: at best one can elaborate models, as in ref. 3. Secondly, the amplitudes A n may be determined, at best, up to a phase per decay mode. In fact, as shown in section 2, a decay to spinning and unstable particles allows to determine, through angular distribution, polarizations and polarization correlations [4,5,2], all products of the type

A λ 1 λ 2 A * λ ′ 1 λ ′ 2
. These products, in turn, allow to determine all moduli of the amplitudes and their phases relative to a given amplitude, taken as a reference. This is not sufficient to solve the linear system (12) with respect to the products a n e iδn . Therefore generally we cannot elaborate tests for TRV through the method described.

However relation ( 12) may be considerably simplified, provided we choose suitable decay modes. Let us consider, for example, the decay [2] Λ b → ΛJ/ψ.

(

) 13
or the following decays of B and B s , already studied, both theoretically [4,5,6,7] and experimentally [8,9,10]:

B + → (J/ψK * + ); (14) B 0 → (J/ψK * 0 ); (15)
B s → (J/ψφ), (J/ψK * 0 ). ( 16) Such decays do not involve isotopic spin, in order to avoid interference among different isospin amplitudes. We examine the scattering corresponding to FSI between the decay products in decays of the type just illustrated. Then the momentum of the "initial" hadrons in the center-of-mass system is 1.6 to 1.7 GeV /c, therefore the scattering involves at least 15 partial waves. Since the decays considered imply orbital angular momenta not greater than 2, we are faced essentially with central collisions. These give rise mainly to deep inelastic collisions, whose shadow constitutes the largest part of elastic scattering. Instead, anelastic reactions -that is, with excitation of one or both decay products to higher mass states -are suppressed, as well as spin-flip elastic scattering, since these processes occur only in peripheral collisions, so as to keep coherence with the initial state.

But deep inelastic collisions imply a complete loss of coherence with respect to the initial particles. Indeed, any such process includes a) an exchange of gluons between the two hadrons; b) hard collisions among partons; c) parton fragmentations and/or recombinations to hadrons; d) gluon exchange among final hadrons. Therefore the sum (12) consists of quite a lot of terms. Moreover it appears logical to assume for the phases in ( 12) a very rapid dependence on the index n (n = m), that is, a sudden variation of the phase-shift from state to state. On the contrary, there is no reason for the terms a n (weak amplitudes) and O mn (kinematic quantities) to depend so strongly on n. Then, for n = m, the sum can be approximated by the integral of the product of a smooth function times a rapidly varying phase. Therefore only the term with n = m survives in that sum, i. e.,

A m = O mm a m e iδm . ( 17
)
This result was obtained also by Wolfenstein [11] with a different line of reasoning.

Eq. ( 17) can be tested by comparing the decays considered with the CP -conjugate ones and assuming CP T symmetry. Indeed, under this assumption, a m differs from a m just by a phase, while O mm and δ m are CP -invariant, since they depend only on strong interactions. Here and in the following the barred quantities refer to the CP -conjugate process, obviously with opposite helicities, denoted synthetically by m. Then we have

|A m | = |A m |, (18)
We stress that this test is not mandatory for the method we are going to suggest for detecting TRV. In order to state tests for TRV, we define, preliminarily, a particular observable that we can extract from analyses of decays, that is, for a given decay mode,

Φ m = arg(A m ) -arg(A m 0 ). (19)
Here, as explained before, A m 0 is conventionally taken to be real. Then we define the following asymmetries:

A CP = Φ m -Φ m Φ m + Φ m , A C = Φ m -Φ m Φ m + Φ m , A P = Φ m -Φ m Φ m + Φ m . (20)
Here, again, the barred quantities refer to the phases of the C-conjugated amplitudes.

Conclusion

We have shown a theoretical condition which implies TRV. Moreover we have suggested tests to be applied to non-leptonic two-body decays involving spinning particles. The tests may be realized by means of analyses of angular distributions, polarizations and polarization correlations of the decay products. Given the wealth of b -b pairs which will be produced at LHC per year (10 12 ), such tests appear to be not so unrealistic. Furthermore, in principle, comparing the results of the suggested tests with the SM predictions could help detecting NP effects, since the condition found is independent of any specific model. Lastly we observe that the condition we require for obtaining evidence for TRV is opposite to the one demanded for detecting CP violations, which needs interference among amplitudes, and therefore, necessarily, more terms in the sum (12) [12].


Lahatra Rakotondrainibe

Yvan Kokar
email: yvan.kokar2@insa-rennes.fr

Gheorghe Zaharia
email: gheorghe.zaharia@insa-rennes.fr

Ghais El Zein
email: ghais.el-zein@insa-rennes.fr

Système de Communications Sans Fil Très Haut Débit à 60 GHz

Cet article présente l'étude et la réalisation à l'IETR d'un système de communications sans fil très haut débit à 60 GHz. Le système utilise une architecture monoporteuse simple. La structure du récepteur proposé est basée sur une démodulation différentielle permettant de réduire l'effet des interférences inter-symboles (IIS) et sur une unité de traitement du signal composée d'un bloc de synchronisation conjointe trame et octet et d'un décodeur correcteur d'erreurs RS (255, 239). La technique de synchronisation utilisée permet d'obtenir une forte probabilité de détection de préambule et une très faible probabilité de fausse alarme. Les premiers résultats de mesures ont montré une bonne qualité de transmission pour une liaison en visibilité directe avec des antennes directives.

Abstract -This paper presents the study and the realization at IETR of a high data rate 60 GHz wireless communications system. The system uses a simple single carrier architecture. The receiver architecture is based on a differential demodulation which minimizes the intersymbol interference (ISI) effect and a signal processing unit composed of a joint frame and byte synchronization block and a conventional RS (255, 239) decoder. The byte synchronization technique provides a high preamble detection probability and a very small value of the false detection probability. First measurement results show a good communication link quality in line of sight environments with directional antennas.

Introduction

La montée en fréquence, en particulier à 60 GHz, constitue une des solutions les plus prometteuses pour accroître le débit (quelques Gbit/s) des futurs réseaux locaux sans fil (WPAN). Des études montrant les applications à 60 GHz, la modélisation des canaux de propagation à 60 GHz et les antennes sont citées dans [START_REF] Smulders | On the design of Low-Cost 60 GHz Radios for Multigigabit-per-Second Transmission over Short Distances[END_REF][START_REF] Collonge | Wideband and Dynamic Characterization of the 60 GHz Indoor Radio Propagation-future Home WLAN Architectures[END_REF]. Récemment, les groupes de normalisation IEEE 802.15.3c, ECMA et VHT ont été formés pour normaliser les WPAN dans la bande de 60 GHz [START_REF] Funada | A design of single carrier based PHY for IEEE 802.15.3c standard[END_REF]. Différentes architectures ont été analysées afin de développer de nouveaux systèmes de communications millimétriques pour des applications multimédia.

Cette étude s'inscrit dans le cadre du projet « Techim@ges ». Ce projet exploratoire s'intéresse à l'étude et la réalisation de systèmes de radiocommunication très haut débit (autour de 1 Gbit/s) à 60 GHz pour de futurs réseaux domestiques. La première application visée est le téléchargement rapide de fichiers volumineux. Le système doit fonctionner en liaison directe dans un environnement domestique. A 60 GHz, l'atténuation provoquée par la traversée d'obstacles tels que les murs est très forte. Compte tenu que la liaison radio à 60 GHz est mono-pièce, une liaison par fibre optique est utilisée pour couvrir toutes les pièces d'un environnement domestique.

Les défauts liés aux blocs radiofréquence (RF) du système (imperfection des composants, bruit de phase, non-linéarité, etc.) peuvent engendrer une dégradation importante des performances. L'intérieur des bâtiments constitue aussi un canal à trajets multiples. Les nombreux obstacles présents (murs, cloisons, plafonds, mobilier) sont autant de surfaces plus ou moins réfléchissantes pour les ondes. L'existence de trajets multiples est la cause du phénomène d'évanouissements et des IIS nécessitant un traitement de signal performant en réception (synchronisation, égalisation, correction d'erreurs), surtout quand le débit est important.

Le suite de l'article contient 3 parties. La partie 2 décrit l'architecture générale du système. La bande de base et le traitement du signal associé sont présentés dans la partie 3. Enfin, la partie 4 traite des résultats de mesure.

(OL), ainsi qu'aux non-linéarités de l'amplificateur de puissance.

En émission, après les opérations de codage de canal et de brassage, le flux d'information est appliqué à l'entrée du codeur différentiel. Les symboles codés vont ensuite moduler une porteuse en FI à 3.5 GHz générée par un oscillateur verrouillé en phase (PLO). Le signal modulé en FI est filtré dans une bande de Pour chaque valeur k on calcule 8 valeurs de cette fonction de corrélation, en fonction du nombre i de bits de l'octet d considérés avant le préambule P amputé de ses derniers i bits. La méthode consiste à effectuer une corrélation entre chaque séquence de 32 bits successifs reçus et le préambule (4 octets). On décale ensuite d'un bit pour chaque corrélateur la séquence de données à analyser. En tout, pour la détection d'un préambule et la synchronisation octet, il y a 8 corrélateurs de 32 bits. Par conséquent, on effectue la détection d'un préambule à l'aide de 32 + 7 = 39 bits (+ 7 à cause des différents décalages possibles d'un octet). Pour obtenir une faible probabilité de non-detection de préambule et de fausse alarme [START_REF] Urkowitz | Signal Theory and Random Processes[END_REF], on tient compte de la périodicité des préambules. En détectant 2 préambules successifs, chaque valeur de la fonction de corrélation entre le préambule et la séquence des données reçues est comparée à un seuil à déterminer. La prise de décision doit se faire à partir de 264 octets successifs (préambule + données + préambule) stockés dans un registre. Dès que les corrélateurs de même rang C k des deux réseaux de corrélateurs indiquent la détection de préambule, l'opération est validée. Par conséquent, à partir de 8 + 7 = 15 bits successifs, on peut assurer l'alignement octet.

Dès que l'alignement octet est assuré et le débrassage effectué, le décodeur RS calcule pour chaque trame reçue une pluralité de syndromes d'erreurs représentatifs du nombre et des emplacements des erreurs détectées dans cette trame. En fonction de ces syndromes, le décodeur détermine les corrections à effectuer, dans la mesure où le nombre d'erreurs ne dépasse pas la capacité de correction du code (8 octets). Finalement, les données sont fournies sur 8 bits à l'interface G.E.

Résultats de mesures

Afin de valider le système, une transmission de données à 875 Mbit/s a été réalisée en visibilité directe. Selon la figure 10, le diagramme de l'oeil obtenu montre une bonne ouverture (pour une distance Tx-Rx de 10 m) traduisant une bonne qualité de transmission. Afin de déterminer la qualité de la liaison, des mesures de TEB ont été effectuées avec deux types d'antennes : antennes cornet décrites précédemment et antennes imprimées avec un gain de 8 dBi et un angle d'ouverture de 30°. Les résultats de TEB sont donnés par la figure 11. Ces résultats de mesure montrent qu'avec des antennes cornets, les performances du système sont améliorées par rapport à des antennes imprimées. Cependant, s'il y a un obstacle, la communication peut être coupée.

Conclusions

Ce papier a présenté la conception et la réalisation d'un système de communications sans fil très haut débit à 60 GHz. Les premiers résultats de mesure ont montré une bonne qualité de transmission en liaison visibilité directe avec des antennes directives. Avec des antennes moins directives, l'ajout d'une égalisation est essentiel pour combattre la sélectivité fréquentielle du canal. De même, il est prévu de mettre en oeuvre des modulations à plus grand nombre d'états pour améliorer l'efficacité de transmission.

2

 2 GHz puis transmis par une fibre optique d'environ 300 m. Ce signal FI est utilisé pour moduler directement le courant d'un laser émettant à 850 nm. Le signal optique est reconverti en signal électrique au moyen d'une diode PIN. A la sortie de la fibre optique, le signal est amplifié et transmis au bloc RF d'émission à 60 GHz. Ce bloc est composé d'un mélangeur, d'un tripleur de fréquence, d'un PLO à 18.83 GHz et d'un filtre passe-bande (59-61 GHz). Le bruit de phase du PLO à 18.83 GHz est de -110 dBc/Hz à 10 kHz de la porteuse. Le signal RF de puissance 0 dBm est transmis via une antenne cornet avec un gain de 22.4 dBi et un angle d'ouverture de 12°. La figure 2 présente le schéma bloc du récepteur.

Figure 2 :

 2 Figure 2 : Architecture du récepteur (Rx) à 60 GHz En réception, le signal RF est d'abord filtré pour limiter l'effet du bruit à l'entrée du récepteur. Le signal obtenu est transposé en FI à 3.5 GHz et filtré. Un contrôle automatique de gain (CAG) ayant une dynamique de 20 dB permet d'assurer un signal avec un niveau quasi-constant à l'entrée du démodulateur. On utilise une démodulation différentielle qui assure à la fois la démodulation et le décodage par transition. Comparée à une démodulation cohérente, la démodulation différentielle est une solution techniquement plus simple, mais moins performante dans un canal à Bruit Blanc Additif Gaussien (BBAG). Par contre, elle s'avère moins sensible à l'influence des trajets multiples et par conséquent aux interférences entre symboles qu'ils peuvent engendrer. Grâce au produit de deux symboles consécutifs, le rapport entre le lobe principal et le deuxième lobe de la réponse impulsionnelle du canal est élevé au carré. Par conséquent, les autres trajets sont fortement atténués par rapport au premier trajet. Un filtrage passe-bas (fréquence de coupure 1 GHz) permet de restituer le signal en bande de base. La récupération d'horloge se fait à partir des données reçues remises en forme. De longues séquences de 0 ou de 1 consécutifs sont à éviter

Figure 3 :Figure 4 :

 34 Figure 3 : Format de la trame et architecture de BB-Tx Le codage de canal est réalisé avec un codeur RS(255,239), travaillant en octet. Lorsque la source de données est en série, une conversion série/parallèle (S/P) est nécessaire. Il est possible d'accepter des données par octet fournies par une interface Gigabit Ethernet (G.E.). La trame à transmettre est formée de 4 octets de préambule, 239 octets de données utiles, 16 octets de contrôle et 1 extra-octet afin d'assurer une trame contenant un nombre d'octets (sauf le préambule) multiple de 8 (pour simplifier les opérations de brassage/débrassage). Comme la trame à transmettre après codage est de 260 octets, il est indispensable de disposer de deux fréquences f 1 et f 2, générées dans le bloc gestion d'horloge, telles que :1 2 1 2 F F f = =100.54 MHz, f = =109.37 MHz. 8 8

Figure 5 :

 5 Figure 5 : Evolution de la fonction de corrélation entre P et T La meilleure situation indiquée à la figure 5 est obtenue pour k = 64. On a : d = [0 0 0 0 0 0 1 0], donc seulement d(7) = 1, les autres bits sont tous égaux à 0. La figure 6 présente l'architecture de BB-Rx. En réception, le flux de données série à 875 Mbit/s en sortie de RHD est parallélisé sur 8 bits. Les procédés de détection du préambule et du contrôle d'alignement sont présentés dans la figure 7.

Figure 6 :

 6 Figure 6 : Architecture de BB-Tx

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7 : Détection du préambule et du contrôle d'alignement

Figure 10 :

 10 Figure 10 : Diagramme de l'oeil à 875 Mbit/s (Tx-Rx : 10 m)

Figure 11 :

 11 Figure 11: Résultats de mesure de TEB (sans codage)

Architecture générale du systèmeLa figure1présente le schéma bloc de l'émetteur. Figure 1 : Architecture de l'émetteur (Tx) à 60 GHz La modulation MDP2-D est simple et permet de résoudre le problème d'ambiguïté de phase en réception grâce à une démodulation différentielle. Comparée à l'OFDM, cette technique est moins sensible à l'instabilité et au bruit de phase de l'oscillateur local

Remerciements

Ce travail est réalisé dans le cadre du projet « Techim@ges » du Pôle de Compétitivité « Images & Réseaux » et du projet « Comidom », avec le soutien financier de la Région Bretagne.