Liang Li

Stéphane Lanteri

Ronan Perrussel
email: perrussel@laplace.univ-tlse.fr

A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equations

Keywords: Maxwell's equations, time-harmonic, discontinuous Galerkin, hybridizable method, domain decomposition, Schwarz algorithm équations de Maxwell, régime harmonique, Galerkin discontinu, méthode hybride, décomposition de domaine, algorithme de Schwarz

A Schwarz-type domain decomposition method is presented for the solution of the system of 3d time-harmonic Maxwell equations. We introduce a hybridizable discontinuous Galerkin (HDG) scheme for the discretization of the problem based on a tetrahedrization of the computational domain. The discrete system of the HDG method on each subdomain is solved by an optimized sparse direct (LU factorization) solver. The solution of the interface system in the domain decomposition framework is accelerated by a Krylov subspace method. The formulation and the implementation of the resulting DD-HDG (Domain Decomposed-Hybridizable Discontinuous Galerkin) method are detailed. Numerical results show that the resulting DD-HDG solution strategy has an optimal convergence rate and can save both CPU time and memory cost compared to a classical upwind flux-based DD-DG (Domain Decomposed-Discontinuous Galerkin) approach.

Introduction

During the last 10 years, discontinuous Galerkin (DG) methods have been extensively considered for obtaining approximate solution of Maxwell's equations, see [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF][START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3d heterogeneous Maxwell equations on unstructured meshes[END_REF][START_REF] Hesthaven | Nodal high-order methods on unstructured grids: I. time-domain solution of maxwell's equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin appproximation of the Maxwell's equations[END_REF]. Thanks to the discontinuity of the approximation, this kind of methods has many advantages, such as adaptivity to complex geometries through the use of unstructured possibly non-conforming meshes, easily obtained high order accuracy, hp-adaptivity and natural parallelism [START_REF] Hesthaven | Nodal discontinuous Galerkin methodsalgorithms, analysis, and applications[END_REF]. However, despite these advantages, DG methods have one main drawback particularly sensitive for stationary problems: the number of globally coupled degrees of freedom (DOFs) is much greater than the number of DOFs required by conforming finite element methods for the same accuracy. Consequently, DG methods are expensive in terms of both CPU time and memory consumption, especially for time-harmonic problems [START_REF] Dolean | Solution of the time-harmonic Maxwell equations using discontinuous Galerkin methods[END_REF]. Hybridization of DG methods [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] is devoted to address this issue while keeping all the advantages of DG methods. The design of such a hybridizable discontinuous Galerkin (HDG) method for the discretization of the system of 3d time-harmonic Maxwell's equations is one of the main objectives of the present work.

HDG methods introduce an additional hybrid variable on the faces of the elements, on which the definition of the local (element-wise) solutions is based. A so-called conservativity condition is imposed on the numerical trace, whose definition involved the hybrid variable, at the interface between neighboring elements [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF][START_REF] Li | A hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF][START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF]. As a result, HDG methods produce a linear system in terms of the DOFs of the additional hybrid variable only. In this way, the number of globally coupled DOFs is reduced. The local values of the electromagnetic fields can be obtained by solving local problems element-by-element. In this work, for the system of 3d timeharmonic Maxwell's equations, we propose a HDG formulation taking the tangential component of the magnetic field as the hybrid variable. We show that the reduced system of the hybrid variable has a wave-equation-like characterization, and the tangential components of the numerical traces for both electric and magnetic fields are single-valued. Moreover, numerical results seems to indicate that the approximate solutions for both E and H have optimal convergence orders.

Though the HDG method results in a smaller linear system than the one associated to a classical upwind flux-based DG method, the size of this system is often too large to be solved by a direct solver as soon as one consider realistic 3d problems. In addition, for very largescale propagation problems, exploiting a multi-processor system is a mandatory path to reduce the solution time and have access to the required memory capacity. Unfortunately, it is also difficult to solve the resulting linear system with standard preconditioned iterative methods [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF]. Therefore, as a second objective of the present work, we propose here to design a hybrid iterative-direct solution strategy for the HDG system, based on domain decomposition principles [START_REF] Smith | Domain decomposition and parallel multilevel methods for elliptic partial differential equations[END_REF]. Our previous works on domain decomposition methods for the solution of the timeharmonic Maxwell equations show promising results of Schwarz-type algorithms through many test problems in 2d and 3d [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF][START_REF] Dolean | Optimized schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a discontinuous galerkin method[END_REF]. A classical Schwarz algorithm which exchanges impedance data between subdomains is studied [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF]. We consider here a similar algorithm and study its adaptation to the HDG discretization framework. We show that the interface system of the Schwarz algorithm is similar to the one obtain with the upwind flux-based DG discretization.

For solving the local sparse linear system within each subdomain in the Schwarz iteration, we adopt a mixed precision arithmetic methodology in order to further reduce the memory requirements and to tackle large problems. A multifrontal sparse direct solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] is used to compute the LU factorization of the coefficient matrix in single precision, and then the L and U factors are employed as preconditioners for a Krylov subspace method which works on double precision DOFs vectors, or an iterative refinement procedure based on these two factors L. Li et al. is performed [BDK + 08] to retrieve double precision accuracy.

This paper is organized as follows. Section 2 defines the considered boundary value problem for the 3d time-harmonic Maxwell's equations and gives some notations. A Schwarz solution algorithm is also stated at the continuous level in this section. In Section 3 we propose a HDG formulation for the one-domain problem and discuss the properties of the reduced system. Section 4 presents the discretization of the domain decomposition algorithm. Numerical examples are given in section 5 with two objectives in mind: first, using a simple (model) problem, the proposed DD-HDG solution strategy is validated and its convergence properties are assessed; second, by considering more complex problems, we study the overall efficiency of the DD-HDG solution strategy and show how much one can gain in terms of CPU time and memory consumption with the new discretization. Finally, we draw some conclusions and state future works in section 6.

2 Problem statement and notations

Boundary value problem

The system of 3d time-harmonic Maxwell's equations is considered

         iωε r E -curl H = -J, in Ω, iωµ r H + curl E = 0, in Ω, n × E = 0, on Γ m , n × E + n × (n × H) = n × E inc + n × (n × H inc ), on Γ a , (1)
where i is the imaginary unit, ω is the angular frequency, ε r and µ r are the relative permittivity and permeability, n is the outward unit normal vector, E and H are the electric and magnetic fields, J is a known current source, and (E inc , H inc ) is a given incident electromagnetic wave. The boundary of the computational domain Ω is ∂Ω = Γ m ∪ Γ a with Γ m ∩ Γ a = ∅. The boundary condition on Γ m indicates a metallic boundary condition (also called perfect electric conductor condition), while the second relation states a Silver-Müller (first order absorbing boundary) condition on Γ a .

Notations

We consider a simplicial mesh T h (consisting of tetrahedra) of the computational domain Ω. We denote by F I h the union of all interior interfaces of T h , by F B h the union of all the boundary interfaces of T h , and by

F h = F I h ∪ F B h . For an interface F = K + ∩ K -∈ F I h
, let v ± be the traces of v on F from the interior of K ± . On this face, we define mean (average) values { • } and jumps • ( • only consider in Section 4) as follows

   {v} F = 1 2 (v + + v -), v F = n + × v + + n -× v -, v F = v + -v -,
where n ± denotes the outward unit norm vector to K ± . For the boundary faces these expressions are modified as (assuming

F = ∂K + ∩ ∂Ω) {v} F = v + , v F = n + × v + , v F = 0,
Inria since we assume v is single valued on the boundaries. Let P p (D) denote the space of polynomial functions of degree at most p on a domain D. For any element K ∈ T h , let V p (K) ≡ (P p (K)) 3 . The discontinuous finite element spaces are then introduced as

V p h = v ∈ (L 2 (Ω)) 3 | v| K ∈ V p (K), ∀K ∈ T h , where L 2 (Ω)
is the space of square integrable functions on the domain Ω. We introduce a traced finite element space

M p h = η ∈ (L 2 (F h )) 3 | η| F ∈ (P p (F )) 3 , (η • n)| F = 0, ∀F ∈ F h . (2)
Note that M p h consists of vector-valued functions whose normal component is zero on any face F ∈ F h . For two vectorial functions u and v in

(L 2 (D)) 3 , we denote (u, v) D = D u • v dx provided D is a domain in R 3 , and we denote u, v F = F u • v ds if F is a two-dimensional face.
Accordingly, for the mesh T h we have

( • , • ) T h = K∈T h ( • , • ) K , • , • ∂T h = K∈T h • , • ∂K , • , • F h = F ∈F h • , • F , • , • Γa = F ∈F h ∩Γa • , • F .
We set

v t = -n × (n × v) and v n = n(n • v),
where v t and v n denote the tangential and normal components of v on a face of unit normal vector n, and v = v t + v n .

Continuous Schwarz algorithm

To simplify the presentation we assume a two-subdomain decomposition of the domain Ω into overlapping or non-overlapping subdomains Ω 1 and Ω 2 . We set Γ 1,2 = ∂Ω 1 ∩ Ω 2 and Γ 2,1 = ∂Ω 2 ∩ Ω 1 . By n i,j (i, j ∈ {1, 2}), we denote the outward unit normal vector to the interface Γ i,j .

The classical Schwarz algorithm in this setting is an iterative method where the (n + 1)-th iterate is defined from the n-th iterate by solving the subdomain problems

           iωε r E (1),n+1 -curl H (1),n+1 = -J in Ω 1 , iωµ r H (1),n+1 + curl E (1),n+1 = 0 in Ω 1 , B n1,2 (E (1),n+1 , H (1),n+1 ) = B n1,2 (E (2),n , H (2),n ), + Boundary conditions on ∂Ω 1 ∩ ∂Ω,            iωε r E (2),n+1 -curl H (2),n+1 = -J in Ω 2 , iωµ r H (2),n+1 + curl E (2),n+1 = 0 in Ω 2 , B n2,1 (E (2),n+1 , H (2),n+1 ) = B n2,1 (E (1),n , H (1),n ), + Boundary conditions on ∂Ω 2 ∩ ∂Ω, (3)
where (E (i),n , H (i),n ) denotes the electromagnetic field in the subdomain Ω i and this local electromagnetic field should converge towards the true electromagnetic field in Ω i as n goes to infinity. Moreover, B n denotes the transmission operator which, in the case of the classical Schwarz algorithm for the system of Maxwell's equations [START_REF] Dolean | Solution of the time-harmonic Maxwell equations using discontinuous Galerkin methods[END_REF] writes as

B n (E, H) = n × E Z r + n × (n × H) = n × E Z r -H t , (4)
L. Li et al.

with Z r = µ r ε r being the impedance. We will comment later on the definition of the transmission operator in terms of the impedance condition.

3 Discretization of the one-domain problem

For sake of simplicity, we omit the volume source term J in what follows but it can be straightforwardly added.

Principles

The discontinuous Galerkin (DG) method seeks an approximate solution

(E h , H h ) in the space V p h × V p h that should satisfy for all K in T h (iωε r E h , v) K -(curl H h , v) K = 0, ∀v ∈ V p (K), (iωµ r H h , v) K + (curl E h , v) K = 0, ∀v ∈ V p (K). (5)
Applying appropriate Green's formulas and replacing the boundary terms by numerical traces E h and H h , we obtain

(iωε r E h , v) K -(H h , curl v) K + H h , n × v ∂K = 0, ∀v ∈ V p (K), (iωµ r H h , v) K + (E h , curl v) K -E h , n × v ∂K = 0, ∀v ∈ V p (K). (6)
It is straightforward to verify that n × v = n × v t therefore we can rewrite (6) using numerical traces defined in terms of the tangential components E t h and H t h

(iωε r E h , v) K -(H h , curl v) K + H t h , n × v ∂K = 0, ∀v ∈ V p (K), (iωµ r H h , v) K + (E h , curl v) K -E t h , n × v ∂K = 0, ∀v ∈ V p (K). (7)
In a classical DG method, we directly formulate the numerical traces to couple the traces of the local fields E h and H h between the neighboring elements

H h = {H h } + α E E h , E h = {E h } + α H H h , (8)
where α E and α H are positive penalty parameters. In the setting of a HDG method, we introduce a hybrid variable Λ h that we define here as

Λ h := H t h , ∀F ∈ F h . (9)
We now want to formulate the local fields in K through (7) assuming that Λ h is known on all the faces of an element K. In order to achieve this, we consider a numerical trace E t h of the form

E t h = E t h + τ K n × (Λ h -H t h ) on ∂K, (10)
where τ K is a local stabilization parameter. Note that we have used the fact that n×H t h = n×H h . Adding the contributions of (7) over all elements and enforcing the continuity of the tangential component of E h , we can formulate a problem which is to find

(E h , H h , Λ h ) ∈ V p h × V p h × M p h such that for all v ∈ V p h and η ∈ M p h        (iωε r E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = 0, (iωµ r H h , v) T h + (E h , curl v) T h -E t h , n × v ∂T h = 0, E t h , η F h -Λ h , η Γa = g inc , η Γa , (11)

Inria

where the last equation is called the conservativity condition. With the definition (10) of E t h , we employ again a Green formula with the second equation of (11), in order to get the following formulation, for all v ∈ V p h and η

∈ M p h      (iωε r E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = 0, (iωµ r H h , v) T h + (curl E h , v) T h + τ n × (H h -Λ h ), n × v ∂T h = 0, n × E h , η ∂T h + τ (H t h -Λ h ), η ∂T h -Λ h , η Γa = g inc , η Γa . (12)
Note that we have used

n × v t = n × v and n × (n × v t ) = -v t , (13)
to obtain (12). Indeed, the last relation of (11) together with the definition (10

) of E t h yields n × E t h , η ∂T h + τ n × n × (H t h -Λ h ) , η ∂T h -Λ h , η Γa = g inc , η Γa ,
which can be rewritten, using (13),

n × E h , η ∂T h -τ (Λ h -H t h ), η ∂T h -Λ h , η Γa = g inc , η Γa .
Moreover, we note that the last relation of (11) for a boundary face F on Γ a is equivalent to, using the fact that

E t h = n × E t h , n × E t h , η F -H t h , η F = g inc , η F , that is n × E h + n × (n × H h ) , η F = g inc , η F ,
which is nothing else than the Silver-Müller boundary condition in (1) with the numerical traces E h and H h .

Remark 3.1 The equivalent numerical scheme for the Laplacian is called LDG-H [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF]. It also coincides with the scheme called HDGII in [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF].

Characterization of the reduced problem

For any η ∈ M h , we denote by (E η h , H η h ) the approximate solutions whose restriction to an element K of T h is the solution to the local problem given by, for all v ∈ V p h (K)

(iωε r E η h , v) K -(H η h , curl v) K + η, n × v ∂K = 0, (iωµ r H η h , v) K + (curl E η h , v) K + τ n × (H η h -η), n × v ∂K = 0, (14)
We rewrite the reduced problem as

a h (Λ h , η) = b h (η), ∀η ∈ M p h , (15)
with

a h (Λ h , η) = n × E Λ h , η F h -Λ h , η Γa , b h (η) = g inc , η Γa .

RR n°8251

In order to give an explicit formulation for a h (Λ h , η) and analyze its properties, we follow the approach adopted in 2d [START_REF] Li | A hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF]. Summing the contributions of ( 14) over all the elements of T h , we obtain the following formulations by recalling the definition of •

(iωε r E η h , v) T h -(H η h , curl v) T h + η, n × v F h = 0, (iωµ r H η h , v) T h + (curl E η h , v) T h + τ n × (n × (η -H η h )) , v F h = 0. (16)
Note that the first relation of ( 16) is obtained by taking the summation in the complex conjugation of the first relation of ( 14). The sesquilinear form in (15) can now be obtained

a h (Λ h , η) = n × E Λ h h , η F h -Λ h , η Γa = n × E Λ h h , η F h + n × ( E Λ h h -E Λ h h ) , η F h -Λ h , η Γa = (iωε r E Λ h h , E η h ) T h +(curl E Λ h h , H η h ) T h + n × ( E Λ h h -E Λ h h ) , η F h -Λ h , η Γa ,
by the first relation of ( 16), taking v = E Λ h h , and

a h (Λ h , η) = (iωε r E Λ h h , E η h ) T h -(iωµ r H Λ h h , H η h ) T h + τ n × (n × (Λ -H Λ h h )) , (η -H η h ) F h -Λ h , η Γa
by the second relation of ( 16) and considering the definition (10) of E h . Finally we obtain

a h (Λ h , η) = (iωε r E Λ h h , E η h ) T h -(iωµ r H Λ h h , H η h ) T h -τ n × (Λ h -H Λ h h ), n × (η -H η h ) ∂T h -Λ h , η Γa . (17)
When ε r and µ r are real-valued, we can infer that the corresponding coefficient matrix K is complex-valued symmetric and all the eigenvalues lie in the left half-plane of the complex plane, because of the third and fourth terms in (17). Moreover, the first two terms define the imaginary part of K if ε r , µ r and τ are real-valued. Note that K can be straightforwardly computed in a classical finite-element way. This matrix K has usual characteristics of the discretization of the wave equation: it is symmetric but indefinite as soon as ω is sufficiently large. A probably more classical formulation is obtained if we multiply (17) by iω in order to obtain

iωa h (Λ h , η) = (µ r (-iωH Λ h h ), (-iωH η h )) T h -ω 2 (ε r E Λ h h , E η h ) T h -iω τ (Λ h -H Λ h h ), (η -H η h ) ∂T h -iω Λ h , η Γa .

About the numerical fluxes

The conservativity condition holds on all the interior faces, which means E t h , η F I h = 0. From the choice of approximation spaces and if the approximation order p is uniform, we can infer that E t h = 0 on every interior face of a conforming mesh. Substituting E h with the expression (10), we have

E t h + τ n × (Λ h -H t h ) = E t h + τ K + H t+ h + τ K -H t- h -(τ K + + τ K -)Λ h = 0 on F I h .

Inria

Solving for Λ h , we obtain (assuming

τ K + + τ K -= 0) H t h = Λ h = 1 τ K + + τ K - (τ K + H t+ h + τ K -H t- h ) + 1 τ K + + τ K - E t h on F. (18)
Substituting ( 18) into (10), and taking the tangential components on both sides, we have

E t h = E t+ h = E t- h = 1 τ K + + τ K - (τ K -E t+ h + τ K + E t- h ) + τ K + τ K - τ K + + τ K - H t h on F. (19)
Thus, the tangential components of the numerical traces for both E h and H h fields are singlevalued.

Remark 3.2 If we take τ uniformly equal to 1, the numerical traces coincide with the upwind flux-based DG formulation considered in [START_REF] Dolean | Solution of the time-harmonic Maxwell equations using discontinuous Galerkin methods[END_REF]. Moreover, equipped only with Silver-Müller boundary conditions, both approaches provide the same discrete solution. For a more general upwind flux-based DG formulation, i.e. adapted to inhomogeneous media as proposed in [START_REF] Hesthaven | Nodal high-order methods on unstructured grids: I. time-domain solution of maxwell's equations[END_REF], the numerical traces coincide with the HDG method by the simple choice τ = Z r .

Assuming the interpolation order is p on every element, we obtain the number of the globally coupled degrees of freedom Upwind flux-based DG formulation

: (p + 1)(p + 2)(p + 3)N e , HDG formulation : (p + 1)(p + 2)N f ,
with N e be the number of elements and N f the number of faces. For a simplicial mesh, N f ≈ 2N e , thus the ratio of the globally coupled degrees of freedom of HDG method over DG method is roughly 2/(p + 3).

Discretization of the multi-domain problem

In this section we formulate a HDG-based discrete variant of the Schwarz domain decomposition algorithm introduced in section 2.3. We first make some preliminary remarks concerning the definition of the interface (or transmission) conditions.

Preamble

As a preliminary remark, note that an objective of the HDG scheme is that the hybrid variable and the numerical trace satisfy on each interior face of the computational domain the following conditions n × E h = 0 and

Λ h = 0. ( 20
) Recall that n × E h = n × E t h and E t h = -n × (n × E h ). In a given subdomain Ω i , the hybrid variable Λ (i)
h is by definition single-valued on each face and thus for any face in the interior of the subdomain, the condition Λ (i) h = 0 is automatically satisfied. However at the interface between two subdomains Ω i and Ω j , the hybrid variables is a priori double-valued and the condition

Λ h = Λ (i) h -Λ (j)
h = 0 has to be explicitly enforced. Consequently, on any face of an interface between two subdomains Ω 1 and Ω 2 , we have to enforce (20) but we can equivalently consider

n × E h + z 1 Λ h = 0 and n × E h + z 2 Λ h = 0, (21)
L. Li et al.

with z 1 and z 2 such that z 1 -z 2 = 0 (more complex conditions could also be considered especially in the framework of optimized Schwarz [START_REF] Dolean | Optimized schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a discontinuous galerkin method[END_REF], but the only important thing is that ( 21) is equivalent to (20) on the interface between subdomains). In the case of the Schwarz algorithm (3) based on the transmission operator in terms of the impedance condition (4), and using the definitions E h and H t h (and again the fact that n × E t h = n × E h ), we want to enforce

           n 12 × (E (1) h + τ n 12 × (Λ (1) h -H t, (1) h
)) -Z (1) r Λ

(1)

h = n 12 × (E (2) h + τ n 21 × (Λ (2) h -H t,(2) h )) -Z (1) r Λ (2) h , n 21 × (E (2) h + τ n 21 × (Λ (2) h -H t,(2) h )) -Z (2)
r Λ

(2)

h = n 21 × (E (1) h + τ n 12 × (Λ (1) h -H t,(1) h )) -Z (2) r Λ (1) h , (22)
or after simplifications

               n 12 × E (1) h + τ (H t,(1) h -Λ (1) h ) -Z (1)
r Λ

(1)

h = n 12 × E (2) h -τ (H t,(2) h -Λ (2) h ) -Z (1) r Λ (2) h , n 21 × E (2) h + τ (H t,(2) h -Λ (2) h ) -Z (2) r Λ (2) h = n 21 × E (1) h -τ (H t,(1) h -Λ (1) h ) -Z (2) r Λ (1) h . (23) If we choose as definition Λ h = Λ (1) h -Λ (2)
h and remark that n 12 = -n 21 , the conditions (23) are equivalent to

n × E h -Z (1) r Λ h = 0 and n × E h + Z (2) r Λ h = 0. (24)
It is, as soon as Z

(1) r + Z

(2) r = 0, equivalent to (20) on the interface between subdomains Ω 1 and Ω 2 . In the Schwarz iterative algorithm at iteration (n + 1), the transmission conditions will translate as

           n 12 × E (1),n+1 h + τ (H t,(1),n+1 h -Λ (1),n+1 h ) -Z (1) r Λ (1),n+1 h = n 12 × E (2),n h -τ (H t,(2),n h -Λ (2),n h ) -Z (1) r Λ (2),n h , n 21 × E (2),n+1 h + τ (H t,(2),n+1 h -Λ (2),n+1 h ) -Z (2) r Λ (2),n+1 h = n 21 × E (1),n h -τ (H t,(1),n h -Λ (1),n h ) -Z (2) r Λ (1),n h . (25)
Moreover, in the HDG scheme, interface conditions are explicitly enforced independently of the volume condition which is more convenient for the formulation of the domain decomposition algorithm in the discrete case compared to a classical DG scheme.

HDG formulation of the multi-domain problem

Let us now assume that the domain Ω is decomposed into N s subdomains Ω = Ns i=1 Ω i . A superscript i indicates that some notations are related to the subdomain Ω i . Thus,

T (i) h , V (i) h , ∂T (i) h and M (i)
h can be defined from those of T h , V h , ∂T h and M h . We also define Γ

(i) m = Γ m ∩∂Ω i and Γ (i) a = Γ a ∩ ∂Ω i .
Finally, Γ i,j will denote the interface ∂Ω i ∩ Ω j . Discretization of the Schwarz algorithm (3) by the HDG scheme (12) integrates conditions (25) and leads to the following problem which is to find (E

(i),n+1 h , H (i),n+1 h , Λ (i),n+1 h ) (n = 1, 2, • • • ) Inria until convergence, satisfying for all v ∈ V (i) h and η ∈ M (i) h ,                                                                  (iωε r E (i),n+1 h , v) T (i) h -(H (i),n+1 h , curl v) T (i) h + Λ (i),n+1 h , n × v ∂T (i) h = 0, (iωµ r H (i),n+1 h , v) T (i) h + (curl E (i),n+1 h , v) T (i) h + τ n × (H (i),n+1 h -Λ (i),n+1 h ), n × v ∂T (i) h = 0, n × E (i),n+1 h , η ∂T (i) h + τ H t,(i),n+1 h -Λ (i),n+1 h , η ∂T (i) h -Λ (i),n+1 h , η Γ (i) a - j Z (i) r Λ (i),n+1 h , η Γi,j = j n × E (j),n h , η Γi,j - j τ H t,(j),n h -Λ (j),n h , η Γi,j + g inc , η Γ (i) a - j Z (i) r Λ (j),n h , η Γi,j . (26)
For simplicity, we have dropped the superscript p which represents the interpolation order.

Formulation of the interface system

For sake of simplicity we consider a problem with two domains Ω 1 and Ω 2 . The reduced linear system for the solution on the whole domain Ω, issued from the sesquilinear form (17), can take the following form

   K (1) ii 0 K (1) ig 0 K (2) ii K (2) ig K (1) gi K (2) gi K (1) gg + K (2) gg       Λ (1) h,i Λ (2) h,i Λ h,g    =    b (1) h,i b (2) h,i b (1) h,g + b (2) h,g    , (27)
where

K (l) ii K (l) ig K (l) gi K (l) gg
, denotes the matrix assembled on Ω l for l in {1, 2}. The subscript g indicates that the concerned DOFs are on Γ 1,2 , and the subscript i that DOFs are in Ω 1 or Ω 2 but not on Γ 1,2 . Following the lines of the work [START_REF] Magoulès | Optimal discrete transmission conditions for a nonoveralpping domain decomposition method for the Helmholtz equation[END_REF], we can deduce the following proposition.

Proposition 4.1 Let A (1) and A (2) be any matrices such that A (1) + A (2) is invertible. There

RR n°8251

L. Li et al.

is only one pair of Lagrange multipliers (Σ (1) , Σ (2) ) such that the coupled problem

K (1) ii K (1) ig K (1) gi K (1) gg + A (1) Λ (1) h,i Λ (1) h,g = b (1) h,i b (1) h,g + Σ (1) , K (2)
ii

K (2) ig K (2) gi K (2) gg + A (2) Λ (2) h,i Λ (2) h,g = b (2) h,i b (2)
h,g + Σ (2) ,

Σ (1) + Σ (2) -(A (1) + A (2) )Λ (1) h,g = 0, Σ (1) + Σ (2) -(A (1) + A (2) )Λ (2) h,g = 0, ( 28
)
is equivalent to problem (27).
From (28), we are led to solve the following linear system for the interface variables Σ (1) and Σ (2) I I -

A (12) (S (2) + A (2) ) -1 I -A (12) (S (1) + A (1) ) -1 I Σ (1) Σ (2) = A (12) (S (2) + A (2) ) -1 c (2) h,g A (12) (S (1) + A (1) ) -1 c (1) h,g , (29)
with A (12) = A (1) + A (2) and where

S (l) = K (l) gg -K (l) gi (K (l) ii ) -1 K (l)
ig and c

(l) h,g = b (l) h,g -K (l) gi (K (l) ii ) -1 b (l) h,i , l ∈ {1, 2}.
In the present case, the matrix A (l) is assembled by computing the following term

(A (l) ) i,j = -Z (l) r η j , η i Γ1,2 ,
where (η i ) i is a basis of M p h restricted to Γ 1,2 . Thus A (1) + A (2) is invertible if and only if Z

(1) r + Z

(2) r = 0, which is obviously the case here. The interface system (29) is then solved by a Krylov subspace method which is discussed in the following section.

Algorithmic aspects

Matrix-vector product

As mentioned above, one has to solve the interface system (29) by a Krylov subspace method. For this, we need to know how to efficiently perform the following matrix-vector product

I I -A (12) (S (2) + A (2) ) -1 I -A (12) (S (1) + A (1) ) -1 I Σ (1) Σ (2) = W (1)
W (2) .

In practice, this could be done as follows 1. For l = 1, 2, solve independently

K (l) ii K (l) ig K (l) gi K (l) gg + A (l) Λ (l) h,i Λ (l) h,g = 0 Σ (l) .
Note that these systems are solved only for obtaining the hybrid variable in each subdomain.

Inria 2. For l = 1, 2, set Y (l) = Λ (l)
h,g and compute independently X (l) = Σ (l) -A (l) Y (l) .

3. Communicate X (1) , Y (1) to subdomain 2 and X (2) , Y (2) to subdomain 1.

Compute independently

W (1) = Σ (1) + X (2) -A (1) Y (2) , W (2) = Σ (2) + X (1) -A (2) Y (1) .

Definition of the right-hand side

Before starting the Krylov subspace method, we need to define the right hand-side of the interface system (29). This can be achieved as follows 1. For l = 1, 2, solve independently

K (l) ii K (l) ig K (l) gi K (l) gg + A (l) Λ (l) h,i Λ (l) h,g = b (l) h,i b (l) h,g , 2. For l = 1, 2, set Y (l) = Λ (l) h,g = (S (l) + A (l) ) -1 (b (l) h,g -K (l) gi (K (l) ii ) -1 b (l) h,i ) and compute independently X (l) = A (l) Y (l) .
3. Communicate X (1) , Y (1) to subdomain 2 and X (2) , Y (2) to subdomain 1.

4.

Compute independently the part of the right hand-side F corresponding to each subdomain

F (1) = X (2) + A (1) Y (2) , F (2) = X (1) + A (2) Y (1) .

Numerical and performance results

In this section we present numerical examples with two objectives in mind: first, using a simple (model) problem, the proposed DD-HDG solution strategy is validated and its convergence properties are assessed; second, by considering more complex problems, we study the overall efficiency of the DD-HDG solution strategy and show how much one can gain in terms of CPU time and memory consumption with the new discretization. Comparisons are made with the DD-DG solution strategy presented in [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF].

Solution strategies

The interface system (28) is solved by an unpreconditioned BiCGStab(l) [START_REF] Sleijpen | BiCGStab(l) for linear equations involving unsymmetric matrics with complex spectrum[END_REF] Krylov subspace method with l = 6. The threshold of the stopping criteria for the BiCGStab(l) iteration is set to be 10 -6 . Each iteration of the Krylov subspace method requires several matrix-vector products with the interface matrix of system (28). As described in subsection 4.4.1. each matrix-vector product is translated into the solution of the subdomain problems. A mixed precision arithmetic method is employed for the solution to the discretized reduced system. A LU factorization is computed and stored by MUMPS [ADL00] (a multifrontal sparse direct solver) in single precision arithmetic (32 bit). With this factorization we can consider three local solution strategies • DD-bicgl, a BiCGStab(l) (with l = 1) method which works on double precision vectors, preconditioned by the above LU factorization. The linear threshold is set to ε i = 10 -6 .

• DD-gmres, a GMRES(m) [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] (with m = 10) method which works on double precision vectors with the same stopping criteria as DD-bicgl, preconditioned by the above LU factorization.

• DD-itref, an iterative refinement procedure used to recover double precision arithmetic [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF].

Otherwise stated, the DD-itref strategy is used by default. Moreover the numerical simulations have been performed using a Fortran 90 computer code implementing the HDG and upwind flux-based DG methods in the 3d case using a P 1 approximation of the components of the electromagnetic field within each tetrahedron of the mesh T h . The parallel implementation of the DD-DG and DD-HDG solution strategies is based on a SPMD strategy combining a partitioning of the computational domain with a message passing programming based on the MPI standard. The adjacency graph associated to the tetrahedral mesh is partitioned using the MeTiS tool [START_REF] Karypis | Parallel multilevel k-way partition scheme for irregular graphs[END_REF]. In the following tables, the number of subdomains in the domain partitioning is denoted by N s which is also the number of parallel processes (i.e. each subdomain is mapped to one processing unit). Moreover, "#iter" refers to the required number of iterations of the Krylov method used for the solution of the interface system (28), "CPU (min/max)" gives the minimum and maximum values of the per-process CPU time and "REAL" is the total elapsed time. Otherwise stated, parallel simulations with the DD-DG and DD-HDG solution strategies have been performed on a cluster operated by the INRIA Sophia Antipolis -Méditerranée research center which consists of AMD Opteron 6174 nodes (48 × 2.2GHz cores per node and 128GB per node) interconnected by a high performance Infiniband network.

Plane wave propagation in vacuum

We first consider a model problem consisting in the propagation of a plane wave in the vacuum. The computational domain is chosen to be the unit cube Ω =] -0.5; 0.5[ 3 and the Silver Müller absorbing boundary condition is imposed on the whole ∂Ω. The electromagnetic parameters ε r and µ r are set to be 1 everywhere, and the angular frequency is ω = 2π. Finally, the penalty parameter τ is set to be 1. A series of regular tetrahedral meshes are employed, which divide the unit cube into (N x -1) × (N y -1) × (N z -1) little cubes, where N x , N y and N z are the number of grid points on the edges of the unit cube, and then each little cube is divided into 6 tetrahedrons. Table 1 summarizes the error and convergence behaviors of both the HDG-P 1 and HDG-P 2 methods. For these results we have employed a Matlab implementation of the HDG-P 1 and HDG-P 2 methods with nodal Lagrange basis expansions on a tetrahedral element. We also note that these results were presented for the first time in [START_REF] Li | A hybridizable discontinuous Galerkin method for solving 3D time-harmonic Maxwell's equations[END_REF] which was our very preliminary contribution related to the design of HDG methods for the solution of the 3d time-harmonic Maxwell equations. The corresponding simulations are run in sequential mode and the native sparse direct solver of Matlab is used. In Table 1, "Mesh size" denotes the edge length of the tetrahedrons on the edge of the unit cube (we set N x = N y = N z in the construction of the tetrahedral meshes). We observe that the asymptotic convergence orders of the approximate solutions for both E and H are optimal. We now switch to results obtained using the Fortran 90 computer code implementing the HDG and upwind flux-based DG methods combined to the Schwarz domain decomposition solution algorithm. Table 2 summarizes the characteristics of the two finer meshes used to test the performances of the domain decomposition solver. In Table 2, N = N x = N y = N z denotes the number of the equally distributed grid points on each edge of the unit cube; N v , N t and N f respectively denote the number of vertices, tetrahedra and faces of the meshes. In Table 3 of the Schwarz method with different discretization schemes. Using the HDG method and the upwind flux-based DG method we have identical accuracy of the approximate solutions on the same mesh with the same domain partitioning. Again, the HDG method outperforms the upwind flux-based DG method.

Scattering of a plane wave by a PEC cube

The test problem considered here consists in the scattering of a plane wave by a perfectly electrically conducting (PEC) cube of side length C = 1/3m centered at the origin. The artificial boundary on which the Silver-Müller absorbing condition applies is defined by a unitary cube centered at the origin. The frequency of the incident plane wave is ω = 6π and its polarization is such that We make us of a non-uniform tetrahedral mesh consisting of 179,311 vertices and 1,005,932 tetrahedra (see Figure 1 left for a view of the triangulation in a selected cut plane). Figure 1 right shows the contour lines of the real part of the E y component. Performance figures are summarized in Tables 5 and6. We note that the maximum RAM requirement for computing and storing the L and U factors in the case of the DD-HDG strategy is about 4.2 times lower than that of the DD-DG strategy for the 96 subdomain decomposition. However we also observe that, for both strategies, there is a rather large difference between the minimum and maximum RAM requirements (44.2 % increase for the DD-DG strategy and 33.6 % increase for the DD-HDG strategy) despite a perfectly balanced decomposition of the mesh elements using the MeTiS tool. When switching from 96 to 192 subdomains, the maximum RAM requirement of the DD-HDG strategy is reduced by a factor approaching 2.5 (while the maximum of the per subdomain factorization time decreases by a factor 3.4). Finally, an almost perfect parallel speedup is achieved for the DD-HDG strategy which is, in the case of 96 subdomains, 8.

k =   k x 0 0   , E =   0 E y 0   and H =   0 0 H z   . L.

Scattering of a plane wave by a coated PEC cube

The test problem here is similar to these described in Subsection 5.3 except that the PEC cube is coated by a dielectric material. We make use of a uniform tetrahedral mesh consisting of 131,922 vertices and 744,000 tetrahedra. In order to model the coating of the cube we consider that a thickness of elements corresponds to a material with r = 4.0. Contour lines of the E x and E y components in a selected cut plane are visualized in Figure 2 (bottom figures). Performance figures are summarized in Tables 7 and8. We observe a superlinear parallel speedup of the factorization phase when switching from 64 to 128 subdomains in the DD-HDG solution strategy. However this is not the case for the solution phase. Indeed a very poor acceleration is obtained which is explained here by an unfavorable computation to communication ratio in the case N s = 128.

X Y Z REY 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1 -1.1 -1.2

Scattering of a plane wave by an aircraft

The problem under consideration is the simulation of the scattering of a plane wave by an aircraft geometry. We consider several values of the angular frequency of the incident field: ω = 4/3π, L. Li ω = 8/3π and ω = 16/3π (in the following we often refer to the corresponding frequency F with respective values 200, 400 and 800MHz). The computational domain is artificially bounded by a parallelepiped box on which the Silver-Müller condition is imposed. The underlying tetrahedral mesh consists of 619,759 vertices and 3,483,551 tetrahedra. This problem has been solved using the HDG-P 1 method for the discretization of the subdomain problems in the Schwarz algorithm. Indeed, it has not been possible to use the DG-P 1 method because of the memory capacity required for storing the L and U factors associated to the factorization of the subdomain problems in that case, at least for the number of subdomains considered here i.e. 288. The triangulation of the surface of the aircraft is shown on the top part of Figure 3 while the contour lines of |E| on the aircraft surface for an incident plane wave with frequency F = 400MHz are shown in Figure 4. The convergence histories of the DD-HDG solution strategy for a decomposition of the mesh in 288 subdomains are shown in Figure 5 for each of the considered values of the frequency F .

Performance results are summarized in Tables 9 and10. The first of these tables is concerned with the factorization of the subdomain problems. Although the partitioning of the tetrahedral mesh using the MeTiS tool [START_REF] Karypis | Parallel multilevel k-way partition scheme for irregular graphs[END_REF] yields perfectly balanced subdomains in terms of the number of elements, we again note that there is a large difference between the minimum and maximum memory requirements for constructing and storing the L and U factors. These results show that, as the frequency of the incident wave increases, the interface system (28) is solved in less iterations of the BiCGStab(l) method. This problem which involves a total of 83,605,224 unknowns for the components of the electromagnetic field (E, H) 1 , and 42,343,638 unknowns for the components of the hybrid variable Λ h 2 is solved in about 8min on 288 cores when F = 200MHz, and in less than 2min when F = 800MHz.

Exposure of head tissues to a plane wave

We now present the application of the proposed numerical methodology to the simulation of a time-harmonic electromagnetic wave propagation problem in an irregularly shaped and hetero-1 This number is the total number of degrees of freedom for the electromagnetic field which is obtained when discretizing the problem using the UF-DG-P 1 method and a single domain, and is thus computed as 6 × 4 × #elements (i.e. 6 components for the electromagnetic field and 4 degrees of freedom per component, for each tetrahedron).

2 This number is the total number of degrees of freedom for the hybrid variable which is obtained when discretizing the problem using the HDG-P 1 method and a single domain, and is thus computed as 2 × 3 × #faces (i.e. 2 components for the hybrid variable and 3 degrees of freedom per component, for each face). geneous medium. The problem under consideration is concerned with the propagation of a plane wave in realistic geometrical models of head tissues. We refer to [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF] for a detailed description of the problem setting and in particular of the geometrical modeling approach considered for constructing realistic tetrahedral meshes of the tissues starting from medical images. Surface meshes of the skin and the skull are visualized in Figure 6. The underlying tetrahedral mesh that has been used for the simulations contains 309,599 vertices and 1,853,832 elements. The frequency of the incident plane wave is F = 1800MHz and its polarization is such that:

k =   k x 0 0   , E =   0 0 E z   and H =   0 H y 0   .
Albeit this propagation problem clearly involves irregularly shaped domains and non-uniform tetrahedral meshes, it is yet a simplified configuration with regards to the simulations usually used in numerical dosimetry studies of human exposition to mobile phone radiation [START_REF] Bernardi | Specific absorption rate and temperature increases in the head of a cellular phone user[END_REF], for two reasons: a mobile phone geometrical model has not been taken into account in the present simulation setting and, the electromagnetic parameters of the materials are set to artificial values for the purpose of exemplifying the characteristics of the propagation of the plane wave in the head tissues (null conductivity, ε r = 4.0 for the brain, ε r = 6.5 for the cerebrospinal fluid, ε r = 1.5 for the skull and ε r = 4.0 for the skin).

Parallel simulations have been performed on the PlaFRIM (Plateforme Fédérative pour la Recherche en Informatique et Mathématiques) system operated by the INRIA Bordeaux -Sud-Ouest research center. The cluster that we have used consists of nodes with Nehalem Intel Xeon X5550 processors (8 × 2.66GHz cores per node and 24GB per node), interconnected by a high performance Infiniband network. Contour lines of the real part of E z are visualized in Figures 7 and8. Performance results are summarized in Tables 11 and12. Because of the heterogeneity of the propagation medium, the required number of iterations for the Krylov methods to obtain the convergence of the interface system (28) are notably higher than those observed for the previous test problems. Similar remarks can be made for what concerns the memory consumption for Inria construction and storing the L and U factors of the subdomain problems. for the 256 subdomain decomposition of the mesh the maximum of the per subdomain memory consumption for the DD-HDG solution strategy is about 4 times lower than that of the DD-DG solution strategy. We observe a linear parallel speedup between the 128 and 256 subdomains simulations with the DD-HDG strategy, and in the case of 256 subdomains, the computing time of the DD-HDG solution strategy is about 4.9 times lower than that of the DD-DG solution strategy.

Scattering of multiple spheres in a waveguide

We conclude this series of numerical results by considering a configuration to assess the potentialities of the proposed DD-HDG solution strategy in treating very large-size problems. The proposed test problem is a L-shape waveguide-like structure including a set of 8 perfectly conducting spheres aligned with the geometry (see Figure 9). The spheres are identical with radius R = 0.06m. The distance between a sphere and the lateral walls is set to 0.1m. All the faces of the L-shape structure are assumed perfectly conducting at the exception of the face defined by x = -0.42m and y = -0.58m. The tetrahedral mesh that has been used for the simulations Inria contains 802,532 vertices and 4,495,769 elements. The corresponding total numbers of degrees of freedom for the components of the electromagnetic field (E, H) and for the components of the hybrid variable Λ h are respectively 107,898,456 and 54,635,364. The frequency of the incident plane wave is F = 10GHz and its polarization is similar to these proposed in the human head example. This problem has been solved using the HDG-P 1 method for the discretization of the subdomain problems in the Schwarz algorithm. Indeed, it has not been possible to use the DG-P 1 method because of the memory capacity required for storing the L and U factors associated to the factorization of the subdomain problems in that case, at least for the number of subdomains considered here i.e. 288. Performance results are summarized in Tables 13 and14. We could have expected that the high frequency of the incident wave would have facilitated the solution of the interface system (28) if we refer to the convergence histories obtained for the test problem of the scattering by an aircraft (see Figure 5), however the waveguide-like structure of the considered geometry seems to stiffen the problem. Despite this difficulty, the proposed DD-HDG solution strategy is able to solve the problem in about 52min (including the factorization time) which is a rather encouraging performance result for such a large system size.

Conclusions

We have presented a hybrid iterative/direct solution strategy for the large, sparse and complex coefficients algebraic systems resulting from the discretization of the three-dimensional timeharmonic Maxwell equations. This solution strategy relies on a Schwarz-type domain decomposition method making use of natural transmission conditions at the interface between neighboring subdomains. Within each subdomain, the discretization of the time-harmonic Maxwell equations relies on a new hybridizable discontinuous Galerkin scheme formulated on an unstructured tetrahedral mesh. As a result, the proposed numerical methodology is particularly well suited to the simulation of wave propagation problems in irregularly shaped media. Moreover, the local nature of the hybridizable discontinuous Galerkin formulation allows for a natural treatment of heterogeneous media. Moreover, the hybridizable discontinuous Galerkin scheme allows for a dratsic reduction of the number of globally coupled degrees of freedom (within each sudomain) as compared to a classical discontinuous Galerkin scheme such as the one adopted in our previous study [START_REF] Dolea | A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods[END_REF]. Numerical results have shown that the resulting DD-HDG solution strategy has an optimal convergence rate and can save both CPU time and memory cost compared to a classical upwind flux-based DD-DG approach. Our future efforts will be towards the improvement of the numerical efficiency of the Schwarz-type algorithm adopted in this study thanks to the design of discrete optimized interface conditions [START_REF] Dolean | Optimized Schwarz methods for Maxwell equations[END_REF] in the framework of a hybridizable discontinuous Galerkin formulation.

Figure 2 :

 2 Figure 2: Scattering of a plane wave by a coated PEC cube. Contour lines of the real part of E x (bottom left) and E y (bottom right) in a selected cut plane (x = 0.5).

Inria

Figure 3 :

 3 Figure 3: Scattering of a plane wave by an aircraft. Triangulation of the aircraft surface.

Figure 4 :

 4 Figure 4: Scattering of a plane wave by an aircraft. Contour lines of |E| on the aircraft surface (F=400 MHz).

Figure 5 :

 5 Figure 5: Scattering of a plane wave by an aircraft. Convergence history for 288 subdomains using the DD-HDG solution strategy.

Figure 6 :Figure 7 :

 67 Figure 6: Exposure of head tissues to a plane wave. Triangulations of the skin and the skull.

Figure 9 :Figure 10 :

 910 Figure 9: Scattering of multiple spheres in a waveguide. Geometrical setting

Table 1 :

 1 , we compare the CPU time and memory consumption required for the storage of the LU factors using Propagation of a plane wave in vacuum: numerical convergence of the HDG-P 1 and HDG-P 2 methods.

	Inria

Table 2 :

 2 Propagation of a plane wave in vacuum: characteristics of tetrahedral meshes.

	Mesh N s	CPU (min/max)	RAM (min/max)
			HDG	UF-DG	HDG	UF-DG
	M1	16 2.7s/2.8s	26.1s/26.7s	39MB/40MB 157MB/157MB
	-	32 0.8s/0.8s	5.8s/5.9s	15MB/15MB	60MB/60MB
	M2	32 7.3s/8.0s 108.7s/112.1s 81MB/85MB 329MB/329MB

the HDG-P 1 method (denoted by HDG) and the upwind flux-based DG-P 1 method (denoted by UF-DG). We observe notable gains resulting from the new HDG discretization in both CPU time and memory requirement. Thanks to the reduction of the number of degrees of freedom, with the HDG method we can deal with larger problems. In Table

4

, we compare the performances

Table 3 :

 3

Propagation of a plane wave in vacuum: comparisons between DD-HDG and upwind flux-based DD-DG solution strategies based on computation times and memory requirement for storing the L and U factors.

Table 4 :

 4 Li et al. Propagation of a plane wave in vacuum: comparisons between DD-HDG and upwind flux-based DD-DG solution strategies based on memory consumption and computing times.

	Mesh N s Local	#iter	CPU (min/max)	REAL
			solver		HDG	UF-DG	HDG UF-DG
	M1	16 DD-bicgl	7	38.9s/39.2s 103.2s/104.9s 40.3s	107.1s
		-	DD-gmres	7	20.5s/20.8s	53.9s/54.2s	21.4s	55.3s
		-	DD-itref	7	13.2s/13.5s	34.9s/35.3s	13.9s	36.2s
		32 DD-bicgl	9	24.4s/24.8s	54.3s/55.2s	25.6s	55.6s
		-	DD-gmres	9	13.0s/13.3s	28.5s/28.7s	13.7s	29.5s
		-	DD-itref	9	7.3s/7.6s	18.2s/18.4s	7.9s	19.0s
	M2	32 DD-bicgl	9	89.8s/91.2s 261.5s/267.7s 94.0s	272.8s
		-	DD-gmres	9	47.2s/47.9s 137.0s/137.6s 49.5s	140.3s
		-	DD-itref	9	32.4s/33.8s	87.2s/88.4s	35.2s	90.4s

  Figure 1: Scattering of a plane wave by a PEC cube. Discretization mesh (left) and contour lines of the real part of E y (right) in a selected cut plane (x = 0.5).

	Method N s CPU (min/max) RAM (min/max)
	UF-DG	96	743s/746s	3616MB/5161MB
	HDG	96		47s/48s	906MB/1211MB
	-	192		13s/14s	349MB/488MB
						5 faster than the
	DD-DG strategy.
			Z			Z
					REY
			X	Y	-0.1	X	Y
					-0.15
					-0.2
					-0.25
					-0.3
					-0.35
					-0.4
					-0.45
					-0.5
					-0.55
					-0.6
					-0.65
					-0.7
					-0.75
					-0.8
					-0.85
					-0.9
					-0.95
					-1
					-1.05
					-1.1
					-1.15
					-1.2
					-1.25
					-1.3
						Inria

Table 5 :

 5 Scattering of a plane wave by a PEC cube. Computing times and memory requirement for storing the L and U factors.

	Method N s #iter CPU (min/max) TOTAL
	UF-DG	96	14	1117s/1120s	1131s
	HDG	96	15	128s/130s	133s
	-	192	16	64s/65s	67s

Table 6 :

 6 Scattering of a plane wave by a PEC cube. Computing times for the solution.

Table 7 :

 7 et al. Scattering of a plane wave by a coated PEC cube. Computing times and memory requirement for storing the L and U factors.

	Method N s CPU (min/max) RAM (min/max)
	UF-DG 128	94s/144s	1115MB/1628MB
	HDG	64	44s/71s		1126MB/1479MB
	-	128	13s/19s		436MB/588MB
	Method N s #iter CPU (min/max) TOTAL
	UF-DG 128	24	572s/643s	661s
	HDG	64	21	249s/274s	284s
	-	128	23	195s/216s	229s

Table 8 :

 8 Scattering of a plane wave by a coated PEC cube. Computing times for the solution.

Table 9 :

 9 Scattering of a plane wave by an aircraft. Computing times and memory requirement for storing the L and U factors.

	Frequency Method N s #iter CPU (min/max) TOTAL
	200MHz	HDG	288	17	269s/455s	484s
	400MHz	-	-	10	143s/232s	243s
	800MHz	-	-	4	59s/96s	98s

Table 10 :

 10 Scattering of a plane wave by an aircraft. Computing times for the solution.

	RR n°8251

  L. Li et al.

		Z			Z
		Y	X		Y	X
				Y
				Z	X
	Figure 8: Exposure head tissues to a plane wave. Contour lines of real part of E z in selected
	planes.
	Method N s CPU (min/max) RAM (min/max)
	UF-DG 256	207s/210s	2174MB/3178MB
	HDG	128	48s/49s	1435MB/2004MB
	-	256	14s/15s	566MB/784MB

Table 11 :

 11 Exposure of head tissues to a plane wave. Computing times and memory requirement for storing the L and U factors.

	Method N s #iter CPU (min/max) TOTAL
	UF-DG 256	68	760s/764s	766s
	HDG	128	63	317s/319s	321s
	-	256	71	156s/157s	158s

Table 12 :

 12 Exposure of head tissues to a plane wave. Computing times for the solution.

Table 13 :

 13 Scattering of multiple spheres in a waveguide. Computing times and memory requirement for storing the L and U factors.

	Method N s #iter CPU (min/max) TOTAL
	HDG	288	63	2514s/2669s	2714s

Table 14 :

 14 Scattering of multiple spheres in a waveguide. Computing times for the resolution.

RR n°8251


Influence of Inertial Stimulus on Visuo-Vestibular Cues Conflict for Lateral Dynamics at Driving Simulators

Keywords: Driving simulator, Lateral dynamics, Visuo-vestibular interaction, Virtual reality, Sensory cues conflict

Introduction

Maintaining a sustainable reality of the represented dynamics is a very difficult and sometimes even an impossible issue by using driving simulators. The major leading reasons of this problem are the constrained workspace of the driving simulator and whether a motion platform exists integrated with the driving simulator. The first driving simulators were fixed-base and the simulation was principally only realized by the visual stimulus [START_REF] Stratulat | Does tilt/translation ratio affect perception of deceleration in driving simulators?[END_REF][START_REF] Bertin | Visuo-vestibular interaction in the reconstruction of travelled trajectories[END_REF] to constitute the self-motion perception. This perception is based on the displacement of visual scene flow on the retina referring to the information about the velocity, direction of the motion and the relative distances [START_REF] Bremmer | Linear vestibular self-motion signals in monkey medial superior temporal area[END_REF].

For the static simulator case, illusory self-motion 'vection' often occurs because the driver is stationary and the visual scenario is moving [START_REF] Lepecq | Galvanic vestibular stimulation modifies vection paths in healthy subjects[END_REF][START_REF] Kolasinski | Simulator sickness in virtual environments[END_REF][START_REF] Berthoz | Perception of linear horizontal self-motion induced by peripheral vision (linearvection) basic characteristics and visualvestibular interactions[END_REF][START_REF] Dizio | Perceived self-motion elicited by postrotary head tilts in a varying gravitoinertial force background[END_REF][START_REF] Draper | The adaptive effects of virtual interfaces: vestibulo-ocular reflex and simulator sickness[END_REF][START_REF] Hettinger | Vection and simulator sickness[END_REF][START_REF] Hettinger | Virtually induced motion sickness in virtual environments[END_REF][START_REF] Hettinger | Handbook of virtual environments Design implementation and applications[END_REF][START_REF] Mccauley | Cybersickness-perception of self-motion in virtual environments[END_REF].

It is obvious that inertial restitution addresses a significant role to maintain a developed fidelity of the driver behaviours on driving simulators. The dynamic simulators are being used since the mid 1960's (Stewart platform) [START_REF] Stewart | A platform with six degrees of freedom[END_REF] firstly for the flight simulators, then the use has spread to the automotive applications [START_REF] Haro | Multi-structure approach for complex systems identification[END_REF][START_REF] Aykent | Study of the influence of different washout algorithms on simulator sickness for a driving simulation task[END_REF][START_REF]MTi and MTx User Manual and Technical Documentation[END_REF][START_REF]Statistics Toolbox User's Guide[END_REF][START_REF] Berger | Simulating believable forward accelerations on a stewart motion platform[END_REF][START_REF] Berger | Simulating believable forward accelerations on a stewart motion platform[END_REF][START_REF] Oman | Motion sickness: a synthesis and evaluation of the sensory conflict theory[END_REF]. The utilization scope diversifies from driver training to research purposes such as vehicle dynamics control, advanced driving assistance systems (ADAS) [START_REF] Meywerk | Einfluss der Fahrdynamikregelung auf die Sicherheit von N1-Fahrzeugen bei unterschiedlichen Beladungszuständen[END_REF].

Subjects prefer verbally the dynamic platform rather than the static case [START_REF] Colombet | Motion cueing: What is the impact on the driver's behavior[END_REF][START_REF] Parrish | Comparison of a linear and a non linear washout for motion simulator utilizing objective and subjective data from transport landing approaches[END_REF][START_REF] Reid | Simulator Motion-Drive algorithms: a designer's perspective[END_REF][START_REF] Hall | The need for platform motion in modern piloted flight training simulators[END_REF]. Driving simulation sickness was assessed between dynamic and static simulators in some studies [START_REF] Colombet | Motion cueing: What is the impact on the driver's behavior[END_REF][START_REF] Curry | Kennedy ssq results: fixed-vs. motion-based FORD simulators[END_REF][START_REF] Watson | A Synthesis of Simulator Sickness Studies Conducted in a High-Fidelity Driving Simulator[END_REF]. It was declared that, simulation sickness was coincided lower when using dynamic simulators rather than static simulators respectively [START_REF] Draper | The adaptive effects of virtual interfaces: vestibulo-ocular reflex and simulator sickness[END_REF][START_REF] Curry | Kennedy ssq results: fixed-vs. motion-based FORD simulators[END_REF][START_REF] Watson | A Synthesis of Simulator Sickness Studies Conducted in a High-Fidelity Driving Simulator[END_REF]. [START_REF] Siegler | Sensorimotor integration in a driving simulator: contributions of motion cueing in elementary driving tasks[END_REF] stated that if the motion platform was activated, the bias in reaching increased levels of decelerations was reduced strongly comparing to inactivated platform case for a braking maneuver [START_REF] Siegler | Sensorimotor integration in a driving simulator: contributions of motion cueing in elementary driving tasks[END_REF] which was an indicator for avoidance of visuo-vestibular cues conflict. Berger et al. investigated the believability of the forward acceleration on a Stewart motion platform [START_REF] Berger | Simulating believable forward accelerations on a stewart motion platform[END_REF][START_REF] Berger | Simulating believable forward accelerations on a stewart motion platform[END_REF].

However, there is a lack of publications on reality of lateral dynamics in terms of multi sensory levels (vehicle model: visual lateral acceleration a yv =a y_veh , motion platform: inertial lateral acceleration, human head level: vestibular lateral acceleration ay_vest=ay, sensed cues). Because of that fact, the reality of lateral dynamics in absence and in presence of the motion platform that yield to two conditions was surveyed in this article. This paper surveys if there is any correlation between the visuovestibular level accelerations in case of static and dynamic simulators. Visuo-vestibular level accelerations stand for the real-time registered and measured visual and vestibular level accelerations.

Thus it aims to prove the reduction of simulator sickness by using motion platform in terms of visuo-vestibular level acceleration cues proximity (positive correlation).

This paper explains the effect of having an inertial stimulus (motion platform) for driving simulators on proximity to the reality for the sensed lateral dynamics with respect to the measurements and the perceptual fidelity using a questionnaire technique. To assess this objectively, the vestibular and vehicle level lateral accelerations (a y,sensed =ay_ vest, a yv =a y_veh ) were saved by using a motion tracking sensor and SCANeR studio software respectively. A confidence interval of 95% was chosen to test the correlations (Pearson's correlation) and to fit models for the distributions of the visual-vestibular lateral accelerations with the multiple linear regression between the conditions of static (N=16) and dynamic (N=21) platform cases in terms of visuo-vestibular level lateral accelerations for the group of subjects (N=37). The results showed that the dynamic platform provides a higher lateral dynamics reality (positive correlation with an incidence of 90.48% for N=21) compared to the static configuration (negative correlation with an incidence of 50% for N=16) from Pearson's correlation and a better fitted model and a lower visuo-vestibular cues' conflict for the dynamic (R 2 =0.429, the model is positive sloped, N=21) condition comparing to the static one (R 2 =0.072, the model is negative sloped, N=16) from the multiple linear regression models. A two-tailed Mann Whitney U test yielded that the U computed (2139)>U expected (1300.5) as p<0.0001, there was a significant difference between the sensed lateral acclerations for the static and dynamic platform cases. Disorientation related perception had positive correlations with the vestibular sensed lateral accelerations for the static condition whereas they were negatively correlated in the dynamic case. As conclusion, the dynamic platform presented a reduced level of motion sickness depending on the sensory conflict theory and the perception fidelity studies approved that dizziness was found to have a significant positive correlation with the vestibular level measured lateral acceleration in the static platform (r=0.293, p=0.037<0.05).

Here, the visual level acceleration refers to the lateral acceleration values registered from the vehicle model that moves in the visual environment for the driving simulator. Whereas the vestibular level acceleration represents the subjects' head lateral accelerations connected to the right ear by using a head phone. Also, perception fidelity was investigated by correlating the vestibular level sensed lateral accelerations with the subjective impression which was measured by using the proposed questionnaire.

Materials and Methods

The research method tackled in this study was to compare the motion platform's contribution on visuo-vestibular acceleration conflicts.

Furthermore the subjects were asked to neither steer the steering wheel nor use the gas, brake pedals and the gear.

The throughout experiment phases were realized as in "virtual driver" mode in the driving simulation software SCANeR studio (Figure 1).

In the data analysis part, the superposition principle of motion was used for evaluating the sensed lateral dynamics at head (vestibular) level (Figures 2 and5).

The measured lateral acceleration in vestibular level is calculated by equation [START_REF] Stratulat | Does tilt/translation ratio affect perception of deceleration in driving simulators?[END_REF] where was measured from the participants' right ear levels for the same driven scenario for the static and dynamic platforms via using the sensor in figure 5.

Dynamic driving simulator

This research work was performed under the dynamic as well as static operations of the SAAM (Simulateur Automobile Arts et Métiers) driving simulator (Figure 2). The dynamic driving simulator SAAM involves a 6 DOF (degree of freedom) motion system (Figure 3 and table 1) [START_REF] Aykent | Study of the influence of different washout algorithms on simulator sickness for a driving simulation task[END_REF]. It is operated on a RENAULT Twingo 2 cabin with the original control instruments (gas, brake pedals, steering wheel). The visual system is realized by a 150° cylindrical view (Figure 2). With the driving cabin of the simulator, the multi-level measuring techniques are available: vehicle model and motion platform dynamics levels real-time data acquisition via SCANeR studio driving simulation software, vestibular level dynamics real-time data acquisition via XSens motion tracker, arm and neck muscles dynamics measurement via Biopac EMG (electromyography) device, human's center of gravity displacements measuring equipment Technoconcept to check postural stability [START_REF] Stratulat | Does tilt/translation ratio affect perception of deceleration in driving simulators?[END_REF][START_REF] Aykent | The Influence of the feedback control of the hexapod platform of the SAAM dynamic driving simulator on neuromuscular dynamics of the drivers[END_REF][START_REF] Aykent | A LQR Washout Algorithm for a Driving Simulator Equipped with a Hexapod Platform: The Relationship of Neuromuscular Dynamics with the Sensed Illness Rating[END_REF]. Figure 2 illustrates the SAAM moving-base driving simulator. It could be operated as with static or dynamic platform by switching the "motion platform" module off and on respectively. As seen in the figure, in general there are three dynamical systems of the SAAM driving simulator. They are vehicle dynamics, motion platform dynamics (motion system) and human vestibular dynamics (proprioceptive system). By manipulating or controlling the vehicle dynamics that moves in the vision system and the motion platform dynamics via motion drive algorithms, their effect on human vestibular dynamics can be compared.

In this article, the effect of having an inertial stimulus (motion platform is active and passive separately for the same driving scenario in figures 6 and 7) on human vestibular dynamics and the perception fidelity were discussed [START_REF] Haro | Multi-structure approach for complex systems identification[END_REF]. a y,sensed measured to obtain the vestibular level lateral accelerations of the subjects and a yv , which was given in figure 4, registered from the vehicle model driven in real-time via the driving simulation software for the same driven scenario for the static and dynamic platforms via using the sensor in figure 5.

Vestibular level data acquisition

In order to save the acceleration data from the vestibular level, a motion tracking sensor was used (Figure 5). This motion tracker can measure the data such as the roll, pitch, yaw angles and rates as well as the accelerations in X, Y and Z. The data are calibrated due to three dimensional quaternion orientation. The sampling rate for the data registration during the sensor measurements was 20 Hz. For the calibrated data acquisition, the alignment reset was chosen which simply combined the object and the heading reset at a single instant in time. This had the advantage that all co-ordinate systems could be aligned with a single action [START_REF]MTi and MTx User Manual and Technical Documentation[END_REF]. The details about the XSens motion tracking sensor are given in [START_REF]MTi and MTx User Manual and Technical Documentation[END_REF]. a y_vest and a y,sensed are equal to each other and they represent the measurements at the participants' ears as in figure 6 relying on the superposition principle of the translational and the rotational motions based on equation 1.

Vehicle level data acquisition

Vehicle level data registered by SCANeR studio software can be splitted as; command data (steering wheel angle, gas, brake pedal input, etc.), motion platform level (translational and angular accelerations of the hexapod platform),vehicle level data (vehicle dynamics, engine, etc.), frequential analysis of the motion platform and vehicle levels (by using FFT (Fast Fourier Transform)).

Protocol

Two conditions were driven by using "virtual driver". This was a kind of driving of which the driving simulator was driven with "autopilot". In order to achieve this goal, we chose and used a driver handler option that was already appointed as "follow a speed target specified from command data".

The experiment protocol involved two phases of the driving situations as static and dynamic platform conditions on a country road (Figure 6). Figure 6 also depicts the X-Y trajectory which was attempted in the experiment phases. The whole experimental phase was completed with a constant velocity of 60 km/h in 126 seconds (Figure 7).

Subjects

The experimental procedure was done for static and dynamic platform cases. 37 subjects (N=37, 29 males and 8 females) participated in experiments.

For the static platform condition, 16 subjects (N=16, 11 male and 5 female participants) aged (mean: 33.44 years, SD: 7.66 years) and with driving licence experience (mean : 15.03 years, SD: 7.14 years) participated in this phase (SD: standard deviation).

For the dynamic platform condition, 21 subjects (N=21, 18 male and 3 female participants) aged (mean: 31.62 years, SD: 7.33 years) and with driving licence experience (mean: 13.07 years, SD: 6.90 years) participated in this phase (SD: standard deviation).

Data Analysis

The effect of having an inertial stimulus (motion platform) on proximity to the reality was discussed here for the sensed lateral dynamics regarding driving simulators.

In order to assess this, the vestibular and vehicle level lateral accelerations (a y ) were collected by using a motion tracking sensor (Figure 5) and SCANeR studio software respectively.

Pearson's correlation was computed [START_REF] Berger | Simulating believable forward accelerations on a stewart motion platform[END_REF] between the conditions of static and dynamic driving simulator situations in order to assess the visuo-vestibular sensory conflict levels. According this; if the lateral acceleration at vehicle level is negatively correlated to the lateral acceleration at vestibular level, it yields less realistic driving simulation session. And if they are positively correlated with each other, it shows a convergence to the reality.

A matlab code example used for the Pearson's correlation is given below: * [rvisves, pvisves] =corrcoef (ay_vest_classic_filt, ay_veh_static) Code ( * ) shows a matlab code to maintain the correlation coefficients between the vestibular and vehicle level lateral accelerations to establish (Table 3). This process was applied to every single subject and to the both cases. ay_vest_classic_filt: vestibular level lateral acceleration which was registered at 20 Hz with classical motion cueing algorithm (dynamic platform) ay_veh_static: vehicle level lateral acceleration which was registered at 20 Hz rvisves: correlation coefficient between the vestibular and vehicle level lateral accelerations pvisves: probability (p-value) Multiple linear regressions were assigned to evaluate the sensory cue conflict in both situations (Figures 1 and2) and their difference was compared with a two-tailed Mann-Whitney U test, which is a nonparametric hypothesis test, by using XLSTAT statistics software.

Finally, the subjective evaluations (questionnaire) and objective measurements from the head level (a y_vest =a y_sensed ) were checked in terms of Pearson's correlation coefficient in order to discuss the perception fidelity for the both cases.

Results and Discussion

In this section, the associations of lateral accelerations on the vehicle and the vestibular levels (Figure 8) and their level of significance were discussed as of having and not having the motion platform during the driving simulator operations. Figure 8 explained briefly the impact of the inertial stimulus (motion platform) as an example for a subject. The blue curve illustrated the vestibularly sensed lateral acceleration from the experiments which could also be computed from equation 1. The red curve was depicting the lateral acceleration of the vehicle center of gravity and it was same for the both cases.

Depending on these graphs, it was seen that the lateral acceleration sensed by the vestibular organ indicated a high disparity from the vehicle lateral acceleration between the curvatures (45-105s), during the second curvature turn (105-110s) and the end sections (110-126s) (Figures 678) for the static platform condition whereas there was a close match of visuo-vestibular cues for the dynamic platform case.

Moreover between 65-85s. and on the second curvature turn, the high frequent motion was coincided in a higher severity at the vestibular level of the participants for the static platform compared to the dynamic condition (Figure 8).

The visuo-vestibular lateral accelerations' gap reduction during the second curvature turn (105-110s) was sourcing from the onset cueing in general. The close visuo-vestibular lateral acceleration fit in the end sections (110-126s) (Figures 678) arose by the tilt coordination and the time delays which were integrated in real-time as seen in figure 4, tables 1 and 2 for the dynamic platform condition whereas there was a mismatch of visuo-vestibular cues for the static platform case.

Moreover between 65-85s and on the second curvature turn, the high frequent motion was coincided in a higher severity at the vestibular level of the participants for the static platform compared to the dynamic condition (Figure 8). The visuo-vestibular lateral accelerations' gap reduction in this period originated from the presence of onset cueing for the dynamic platform case.

Table 3 indicated the impact of the motion platform on the reality of the lateral acceleration. Due to this table, 8 out of the 16 subjects had a negative correlation in terms of 'vestibular lateral accelerationvehicle CG lateral acceleration' who participated in the condition of static platform. Those negative correlations emphasized the sensory cue conflict (50% of the sensory conflict incidence).

Due to this table, 8 out of the 16 subjects had a positive correlation (50% of the sensory conflict avoidance incidence) for the static platform. It meant that; apart from this subject, the rest of them yielded a sensory conflict (negative correlation between vestibular lateral acceleration and vehicle CG lateral acceleration) and divergence from the reality of the lateral dynamics represented on the driving simulator.

For the dynamic case, 19 out of the 21 subjects showed positive correlation. It gave an incidence of 90.48% for avoiding the sensory cue conflict based motion sickness. For the dynamic platform situation, there was only 2 people who indicated a negative correlation. It can be interpreted as a visuo-vestibular sensory conflict with an incidence of 9.52%.

The negative correlations between visual and vestibular lateral accelerations in the static case were based on the higher level of the lateral head movements which could be realized to compensate the lack of the motion platform influence. In other words, it allowed the gap to increase between visual and vestibular level lateral accelerations when the vehicle level lateral accelerations were the identical for the boths cases (Figure 8 and table 3).

Figure 9 illustrated the reality of the lateral dynamics depending on the effect of the inertial cues (motion platform) which were given in figure 1. Figure 9 was also a filled contour mapping for the data yielded in table 3.

These graphs were depicted with respect to the Pearson's correlation coefficients (r). When the correlation coefficient r was near to "1 (red colour)", it meant that the driving simulation was representing the lateral acceleration more real. In other words, it presented less contradicting visuo-vestibular cues as from lateral acceleration in accordance with sensory cue conflict theory. In contrast; when the r was closer to "-1 (blue colour)", it indicated that the lateral dynamics was represented less realistically. Therefore, higher values of discrepancy were resulted between the lateral accelerations values measured from vehicle (visual cue) and vestibular cue levels.

Code ( ** ) depicted how figure 9 was created in MATLAB in order to check visuo-vestibular accelerations conflict for the static and dynamic operations of the platform by utilizing the registered data from the SCANeR software. A multiple linear regression with a confidence interval of 95%, where the input is the lateral acceleration at visual cues (vehicle lateral acceleration) and the output is the lateral acceleration at vestibular cues (sensed lateral acceleration, Equation 1, was modelled for the static and the dynamic platform cases. Equation 2gives the linear regression model for the static platform condition. Figure 10 shows the cue conflict comparison of the both cases for the groups of the subjects. Blue dots indicate the real data for the visual and vestibular acceleration cues which were saved from the right ear level of the subjects and the vehicle model driven. Black continuous lines depict the linear regression models. Grey continuous and dashed lines illustrate the confidence interval of 95% for the mean and the observed values respectively.

According to figure 10, it can be seen that the real data were distributed far from the fitted model line (the coefficient of determination R 2 =0.072) for the static platform. On the other hand, the real data were more closely gathered to the estimated model line (the coefficient of determination R 2 =0.429) for the dynamic platform. Furthermore, it is resulted that the visual lateral accelerations (realtime vehicle model in SCANeR studio software) are getting closer (the slope of the model is positive: inclined to right hand side) to the vestibular lateral accelerations (real-time XSens record as a module in SCANeR studio software) for the dynamic situation.

Inversely, the visuo-vestibular lateral accelerations are getting far away from each other (the inclination of the model is negative: inclined to left hand side) which mean that the visuo-vestibular cue conflict are increasing.

Lastly, it is yielded that as the U computed (2139)>U expected (1300.5) and the computed p-value (p<0.0001) is lower than the significance level alpha=0.05, one should reject the null hypothesis H0, and accept the alternative hypothesis Ha after applying the two-tailed Mann-Whitney U test to check the difference significancy between the static and dynamic situations of the vestibular level sensed lateral accelerations where the tested hypotheses were given below: H0: The difference of location between the samples from the static (N=16) and the dynamic (N=21) cases is equal to 0.

Ha: The difference of location between the samples from the static (N=16) and the dynamic (N=21) cases is different from 0.

In this research, also the perception fidelity was compared for the two conditions. In order to achieve this goal, Pearson's correlation method was benefited. Three questions were given to the subjects and they were asked to fill in depending on what they perceived during each driving session. The questions regarding the motion sickness (disorientation related) were composed of as the following: Q1-Were you at the point to vomit? Q2-Have you felt nausea? Q3-Have you felt dizziness?

The answers were scaled in equal ten parts as from 1: too little to 10: too strong.

Due to table 4, it can be seen that the a y_vest_sta (vestibular level lateral acceleration at static motion platform) was positively correlated (r=0.143, p=0.318 and r= 0.056, p=0.695 between the sensed lateral acceleration and the propensity to vomit, feeling nausea respectively) with the disorientation related subjective evaluations whereas the a y_ vest_dyn (vestibular level lateral acceleration at dynamic motion platform) was negatively correlated (r=-0.114, p=0.424; r=-0.119, p=0.407 and r=-0.143, p=0.318 between the sensed lateral acceleration and the propensity to vomit, feeling nausea, feeling dizziness respectively) with the subjective perceptions. In general, it shows that the perception proximity to the measured vestibular level lateral accelerations (a y_vest ) was coincided for the static case rather than the dynamic one. Because the visuo-vestibular sensorial cue conflict increased at the static hexapod platform compared to the dynamic one (Table 4 and figures 8-10). In particular, merely the dizziness feeling had a significant positive correlation (r=0.293, p=0.037<0.05) with the sensed lateral acceleration at the static condition.

Conclusion and Future Work

The reality of the represented lateral dynamics on driving Values in bold are different from 0 with a significance level alpha=0.05 Q1 s , Q2 s , Q3 s were the answers given by the participants to the questions of Q1, Q2 and Q3 at the static platform experiments. Q1 d , Q2 d , Q3 d were the answers given by the participants to the questions of Q1, Q2 and Q3 at the dynamic platform experiments. simulators was discussed in this paper. After having completed these experimental phases, due to the Pearson's correlations, it was proved that the dynamic platform provided a closer lateral dynamics representation between real-time vehicle model (visual cues) and realtime vestibular cues levels. It can be concluded that having dynamic platform represented a higher lateral dynamics reality in terms of data acquisition and measurements. For the dynamic case, there was a positive correlation with an incidence of 90.48% for N=21, whereas a negative correlation with an incidence of 50% for N=16 was yielded for the static one.

The multiple linear regression model resulted as a better fit and a positive slope for the visual-vestibular lateral accelerations for the dynamic (R 2 =0.429, N=21) platform, whereas it was yielded as a weak fit and a negative slope for the visual-vestibular accelerations for the static (R 2 =0.072, N=16) platform.

The two-tailed Mann-Whitney U test proved that there was a significant difference between the static and the dynamic cases in terms of the vestibular level lateral accelerations (U computed (2139)>U expected (1300.5), p<0.0001).

The perception fidelity illuminated that there was a significant positive correlation (r=0.293, p=0.037<0.05) between the vestibular level accelerations and feeling dizziness for the static platform simulator whereas there was no significant correlation between the vestibular level accelerations and the disorientation perception.

As prospective researches, the reality of having motion cueing algorithms with feedback control will be surveyed in terms of multisensory level dynamics approach (neuromuscular cues (EMGelectromyography), visual cues (vehicle level), vestibular cues (head level), inertial cues (motion platform level)) and with respect to postural stability.

  translational acceleration at vestibular level (m/s 2 ) ϕ : Roll angle at vestibular level (°) g : Gravitational acceleration (m/s 2 )

Figure 1 :Figure 2 :Figure 3 :+

 123 Figure 1: Visuo-vestibular interaction on dynamic driving simulator.

Figure 4 :

 4 Figure 4: Motion cueing algorithm.

Figure 5 :

 5 Figure 5: Vestibular level data acquisition.

Figure 5

 5 Figure 5 describes the motion cueing algorithm used for the dynamic platform condition in this study. The motion cueing algorithm was included in the SCANeR studio driving simulation software via DLL plugin in order to accomplish the real-time driving experiments with the participations of the subjects.

Figure 6 :Figure 7 :

 67 Figure 6: Country road testing land and vehicle CG's (centre of gravity) X-Y trajectory during the experiments.

Figure 8 : 1 Figure 9 :

 819 Figure 8: Vehicle-vestibular level lateral accelerations at static and dynamic conditions.

)Figure 10 :

 10 Figure 10: Multiple linear regression models of visuo-vestibular lateral accelerations for static (N=16) and dynamic (N=21) platforms.

Table 1 and

 1 Table2give the limits and the classical motion cueing algorithm parameters respectively. DLL plug-in integrated in the software took the constraints of the simulator platform and the classical motion cueing algorithm parameters into account for the dynamic condition.

	DOF	Displacement	Velocity	Acceleration
	Pitch	± 22 deg	± 30 deg/s	± 500 deg/s 2
	Roll	± 21 deg	± 30 deg/s	± 500 deg/s 2
	Yaw	± 22 deg	± 40 deg/s	± 400 deg/s 2
	Heave	± 0.18 m	± 0.30 m/s	± 0.5 g
	Surge	± 0.25 m	± 0.5 m/s	± 0.6 g
	Sway	± 0.25 m	± 0.5 m/s	± 0.6 g

Table 1 :

 1 Limits of each degree of freedom (DOF) of the dynamic driving simulator SAAM.

Vehicle level lateral acceleration (m/S 2 ) Vehicle level lateral acceleration (m/S 2 ) Vestibular level lateral acceleration (m/S 2 ) Vestibular level lateral acceleration (m/S 2

	Active	Model	Active	Model
	Conf. interval (Mean 95%)	Conf. interval (Obs. 95%)	Conf. interval (Mean 95%)	Conf. interval (Obs. 95%)

Table 4 :

 4 Effect of inertial condition on lateral dynamics reality in absence of Visuo-vestibular restitution.

Symbol

Longitudinal Lateral Roll Pitch Yaw


Ghislain Landry

Tsafack Chetsa

Laurent Lefévre
email: laurent.lefevre@ens-lyon.fr

Jean-Patrick Gelas
email: jean-patrick.gelas@ens-lyon.fr

DTN based Management Framework for Green On/Off Networks

Keywords: efficient networks, on/off networks, proactive switch on/off, DTN, quality of service

The increasing cost of powering high performance networking infrastructure has led to the proposal of various energy saving schemes. The On/Off technique, being the most common energy saving scheme, consists of powering down partially or entirely a network infrastructure for energy saving purposes. Despite their capability to achieve great energy savings, On/Off networks experience high packet-loss rates due to the absence of reliability on packet delivery. Moreover, they cannot guarantee any response time to user applications. This paper presents the design and implementation of MFO2N; Experimental results show a correlation between offered quality of service and overall network power consumption, revealing that a tradeoff should be made 1 .

I. INTRODUCTION

Initiatives conducted by research laboratories call for introduction of Green On/Off networks capable of switching off links (either redundant or not) or nodes for power saving purposes. On/Off networks can be divided into two categories: [START_REF] Orgerie | Energy Efficient Bandwidth Reservation for Bulk data Transfer in Dedicated wired Networks[END_REF] Unplanned On/Off networks whose communication opportunities are completely or partially unknown. (2) Planned On/Off networks whose communication opportunities are planned and well known. Our work focuses on bulk data transport dedicated networks with advance bandwidth reservation features. The network is designed to switch off links and nodes not in use by any data transfer to achieve power saving. In this context, to avoid affecting the utilisation of the network as well as degrading network performance, users should be able to send bandwidth reservation requests irrespective of the state of the network. In addition, users should be able to receive a response within a reasonable time. This paper introduces the Management Framework for On/Off Network (MFO 2 N); a framework providing energy-efficient management of information transfer with quality of service in planned On/Off networks. The goal is to minimise the total power consumption of a dedicated network infrastructure by switching off unused nodes and links while making network unavailability transparent to users. For this purpose, we combined the DTN Store and Forward message switching with a proactive switch On/Off mechanism to service users request which may either (1) arrive at a network node where the link to the next hop towards its destination is down (or the next hop itself) or ( 2) timed out in a DTN router's forwarding queue. We present preliminary results that permit us to estimate the power cost of providing a Green On/Off network with management tools. This paper is organized as follows. Section II presents the network management model description. Section III describes the Management Framework for ON/OFF networks. Section IV analyzes evaluation results. Finally section V concludes this paper and proposes some future works.

II. NETWORK MANAGEMENT MODEL DESCRIPTION

The management model we propose associate each network node and link with an on/off agenda [START_REF] Orgerie | Energy Efficient Bandwidth Reservation for Bulk data Transfer in Dedicated wired Networks[END_REF]. The On/Off agenda, local to each node, is actually a chronologically ordered list of On/Off events describing the future network usage. Thus, given a network node N, N maintains an agenda whose elements are of the form (t i , N, ACTION), and (t i ,l, ACTION) where:

• t i is the time when the action ACTION starts.

• l designate a network link connected to N.

• There are two types of actions: On for switch on, and Off for switch off. We also associate each management request with a lifespan or time to live (TTL) fixed by the users and representing the maximum number of seconds after which the user expects to receive a response to its requests. We then leverage the DTN Store and Forward message switching paradigm to guarantee the delivery of users requests. To guarantee on demand delivery time, we combined the DTN Store and Forward message switching with the proactive switch-on mechanism which consists of making links and/or nodes available to service management information information transport.

III. MANAGEMENT FRAMEWORK FOR ON/OFF NETWORKS (MFO 2 N)

We present in this section the Management Framework for On/Off Networks (MFO 2 N) an implementation of the network management model described in Section II. The main objective of the proposed framework is to minimise the extra energy cost

ION db View Network events manager

Agenda manager

Links manager control module

Switching table

Virtual requests manager requests manager of bandwidth reservation requests transport while guaranteeing on-demand response times to users as well as reliable transport of requests. MFO 2 N lays on the Interplanetary Overlay Networking (ION) [START_REF] Jpl D-48259 | Interplanetary Overlay Network (Design and Operation)[END_REF] software, an implementation of the DTN architecture to which we have brought substantial changes. MFO 2 N consists of several functional elements (see Figure 1) whose description is not given in this paper due to lack of space.

IV. EVALUATION

The evaluation of MFO 2 N is performed on a 10 node network topology (see Figure 2) deployed on the large scale experimental platform Grid5000 [START_REF] Cappello | Grid'5000: A large scale, reconfigurable, controlable and monitorable grid platform[END_REF] via ION. We consider several agenda scenarios of a few minutes each in which end nodes, i.e, nodes 1, 2, 3 and 4, send and/or receive bandwidth reservation requests. We then inject in the network after every 10 seconds one bandwidth reservation request whose source and destination nodes are chosen randomly. We investigate the impact of the offered quality of service expressed essentially in terms of requests response time on the overall network power consumption. Figure 3 depicts the cumulative on time per network link for the various on-demand requests response times denoted as TTL in that figure. As expected, requests with short response delay are more power consuming. Figure 3 successfully deliver user data, mainly bandwidth reservation requests on a network which is down while guaranteeing reasonable response times.

V. CONCLUSION

The paper introduced the design of MFO 2 N. We run MFO 2 N Grid5000 platform and showed through many requests scenarios that the designed framework can ensure the transport of management information in Green On/Off networks. We then provided estimates of the extra energy cost of the proposed framework. Collected results exhibit a correlation between the overall network power consumption and the offered quality of service revealing that a trade-off should be made. This work is ongoing, and obtained results may depend on the DTN under-layer communication software. Future work envisage more complex network topologies and the implementation of a Web based management platform.

Figure 1 :Figure 2 :

 12 Figure 1: MFO 2 N architectural model.

Figure 3 :

 3 Figure 3: On time variation with respect to requests response times (TTL).

  also show the ability of MFO 2 N to

		250
										TTL-0
										TTL-20
		200								TTL-30
	Time in second	100 150
		50
		0
		( 1 ,5 )	( 2 ,7 )	( 3 ,8 )	( 4 ,1 0 )	( 5 ,6 )	( 6 ,7 )	( 7 ,8 )	( 8 ,9 )	( 9 ,1 0 )	( 1 0 ,5 )
						Links Ids

This work is supported by the Hemera INRIA Large Wingspan project. Experiments of this paper were performed on the Grid'5000 platform, an initiative from the French Ministry of Research through the ACI GRID incentive action, INRIA, CNRS and RENATER and other contributing partners (http://www.grid5000.fr)

ACKNOWLEDGEMENT

This work is supported by the Hemera INRIA Large Wingspan project. The authors would like to acknowledge the help of Anne-Cécile Orgerie.