Juliette Chabassier
email: juliette.chabassier@inria.fr

Sébastien Imperiale
email: sebastien.imperiale@inria.fr

Space/Time convergence analysis of a class of conservative schemes for linear wave equations

published or not. The documents may come

Space/Time convergence analysis of a class of conservative schemes for linear wave equations

Juliette Chabassier, Sebastien

Version française abrégée

On considère le problème de propagation d'onde suivant : pour une source donnée f ∈ C 0 ([0, T ], H) trouver u(t) ∈ V solution, pour tout t ∈ [0, T ], de ∂ tt u + Au = f dans V , u(0) = 0 dans V, ∂ t u(0) = 0 dans H, où V et H sont deux espaces séparables de Hilbert (la norme sur H est notée | • |) et V le dual de V , l'opérateur A est autoadjoint de V dans V . Considérons donnée une discrétisation en espace s'appuyant sur une famille d'espaces de dimension finie {V h } h>0 avec V h ⊂ V . Le problème semi-discret s'écrit alors : pour un terme source f h ∈ C 2 ([0, T ], V h ) trouver u h (t) ∈ V h , pour tout t ∈ [0, T ] satisfaisant

∂ tt u h + A h u h = f h , u h (0) = 0, ∂ t u h (0) = 0, dans V h ,
où A h est un opérateur de V h dans V h autoadjoint borné (non uniformément par rapport à h) résultant de la discrétisation en espace choisie. Nous supposerons qu'il existe une fonction positive δ(h) = o(h) telle que pour tout h > 0 on ait sup

t∈[0,T ] |u(t) -u h (t)| ≤ δ(h).
Nous étudions les schémas suivants : le θ-schéma (noté TS et dont le schéma saute-mouton explicite est un cas particulier), le schéma Saute Mouton Stabilisé (noté SLF ) et le (θ, ϕ)-schéma d'ordre quatre (noté TPS ). Soit ∆t > 0 le pas de temps et définissons t n := n∆t. Pour toute suite {v k h } k≥0 ⊂ V h , nous notons

[v n h ] ∆t 2 := v n+1 h -2v n h + v n-1 h ∆t 2 , {v n h } θ := θ v n+1 h + (1 -2θ)v n h + θ v n-1 h , n ≥ 1.
Nous cherchons une suite {u k h } k≥0 ⊂ V h telle que u n h approche u h (t n ). Les deux premiers termes sont supoosés donnés et pour tout n ≥ 1,

• TS : [u n h ] ∆t 2 + A h {u n h } θ = f h (t n ). • SLF : [u n h ] ∆t 2 + A h u n h + ∆t 2 16 A 2 h u n h = f h (t n ). • TPS : [u n h ] ∆t 2 + A h {u n h } θ + ∆t 2 θ - 1 12 A 2 h {u n h } ϕ = f h (t n ) + ∆t 2 ∂ tt 12 + θ - 1 12 A h f h (t n ).
Pour chaque schéma on peut définir une énergie discrète associée, une condition de type CFL assure sa positivité en tant que fonctionnelle : ∆t 2 ≤ α ρ(A h ) où α > 0 est positif, le plus grand possible, et dépend du schéma. Nous montrons que sous des hypothèses portant sur la discrtisation spatiale du problème ainsi que sous condition CFL, le résultat de convergence espace/temps suivant : il existe C > 0 indépendant de ∆t et h tel que

• θ-schéma (TS ) et schéma saute-mouton stabilisé (SLF ) sup t n ∈[0,T ] |u(t n ) -u n h | ≤ C ∆t 2 T 2 + δ(h) .
• (θ, ϕ)-schéma (TPS ) sup

t n ∈[0,T ] |u(t n ) -u n h | ≤ C ∆t 4 T 2 + δ(h) .
Notons qu'en particulier ces résultats de convergences sont valables lorsque ∆t 2 = α/ρ(A h ). Pour obtenir ce résultat nous nous ramenons à l'étude de polynômes dont les coefficients sont fonction des paramètres des schémas.

Abstract and semi-discrete wave propagation problem

In this work, we are interested in the simulation of linear wave propagation problems. The most simple example one could think of is given by the following problem: find

u(t) ∈ H 1 0 (Ω), for all t ∈ [0, T ] such that ∂ tt u -∆u = f, u(0) = 0, ∂ t u(0) = 0, in Ω, u = 0 on ∂Ω, (1)
To study wave propagation in a more abstract framework, following [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology -Volume 5 and 6 -Evolution Problems I and II[END_REF], chapter XVIII, we assume given separable Hilbert spaces H and V . The space H is equipped with the scalar product (•, •) H . The corresponding norms are denoted | • | and • respectively. Moreover we assume that V is dense in H, H is identified with its dual H and V is continuously embedded in H. Note that to solve (1) it is standard to choose H = L 2 (Ω) and V = H 1 0 (Ω). We assume given a continuous hermitian bilinear form a :

V × V → R that satisfies a(v, v) ≥ 0, a(v, v) + C o |v| 2 ≥ C a v 2 , ∀v ∈ V, (2)
where (C a , C o ) are real positive scalars. From this bilinear form we construct a self-adjoint positive unbounded operator A : V → V . It is defined by

A : u → Au such that Au, v = a(u, v), ∀v ∈ V.
where V denotes the dual of V and •, • denotes the duality product. We consider the following abstract wave propagation problems: for a given f ∈ C 0 ([0, T ], H) find u(t) ∈ V solution, for all t ∈ [0, T ], of

∂ tt u + Au = f in V , u(0) = 0 in V, ∂ t u(0) = 0 in H. (3)
Existence and uniqueness results for this problem are well known (see again [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology -Volume 5 and 6 -Evolution Problems I and II[END_REF], chapter XVIII): there exists a unique function u solution of (3), it satisfies

u ∈ C 0 ([0, T ], V ), ∂ t u ∈ C 0 ([0, T ], H).
We introduce the family of finite dimensional spaces {V h } h>0 with V h ⊂ V . As usual, the subscript h is devoted to tends to 0 and represents an approximation parameter of V h to V . For each h we define the operator A h as A h : V h → V h and

A h : u h → A h u h such that (A h u h , v h ) h = a h (u h , v h ), ∀v h ∈ V h ,
where the scalar product on V h denoted by (•, •) h as well as the continuous bilinear forms a h (•, •) represent some approximation of (•, •) H and a(•, •) respectively that account for instance for the use of quadrature formulae in the computation of integrals. Similarly we denote by | • | h the norm induced by the scalar product (•, •) h . We assume that there exists a positive constant C H such that

|v h | ≤ C H |v h | h , ∀v h ∈ V h
and we assume that a h is an hermitian bilinear form that satisfies the same property as in equation ( 2)

a h (v h , v h ) ≥ 0, a h (v h , v h ) + C o |v h | 2 h ≥ C a v h 2 , ∀v h ∈ V h .
Note that this assumption implies that the operator A h is self-adjoint and non-negative. Its spectrum, denoted Sp(A h ), is a set of finite number of real non-negative eigenvalues (eigenvalues of multiplicity greater than one are counted accordingly). To each eigenvalue λ we associate a real eigenvector w λ h defined by

(A h w λ h , v h ) h = λ (w λ h , v h ) h , ∀v h ∈ V h , |w λ h | h = 1, such that the family {w λ h } λ∈Sp(A h ) is an orthonormal basis of V h .
The semi-discrete equation we consider reads: for any given source term

f h ∈ C 2 ([0, T ], V h ) find u h (t) ∈ V h , for all t ∈ [0, T ] such that ∂ tt u h + A h u h = f h , u h (0) = 0, ∂ t u h (0) = 0, in V h . (4)
Existence and uniqueness results for this problem are direct consequences of the theory developed in infinite dimensional space (choose H = V = V h in the paragraph above): there exists a unique solution u h of (4), it satisfies u h ∈ C1 ([0, T ], V h ). We introduce the discrepancy error e h (t) = u(t) -u h (t) it depends on the approximation of the scalar product on H, on the bilinear form a h , on the approximation of the space V itself and on the approximation of the source term f . It is expected that e h vanishes when h → 0. Therefore we make the following assumption Hypothesis 1.1 Convergence of the semi-discrete problem. There exists a positive function δ(h) = o(h) such that for all h > 0 we have sup

t∈[0,T ] |e h (t)| ≤ δ(h).
Such a result is proved for continuous finite element approximation in [START_REF] Joly | Effective computational methods for wave propagation[END_REF] in the case a(•,

•) h ≡ a(•, •).
In [START_REF] Durufle | Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain[END_REF] the case of spectral elements is given (see [START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations, State-of-the-art surveys on computational mechanics[END_REF] for an introduction to spectral elements). For these convergence results to hold, the continuous solution must be regular enough in time and space. We do not detail the assumptions here. In what follows we specify the extra regularity in time and space we require for the semi-discrete solution u h in order to state the space/time convergence results, such regularity is a consequence of the space/time regularity of the source term.

Hypothesis 1.2 Stability of the semi-discrete problem. For a given couple (p, q) with p even, we assume that ∂ q t u h (t) ∈ V h for all t ∈ [0, T ] and there exists a constant C p,q such that for all h > 0 we have sup

t∈[0,T ] | A p/2 h ∂ q t u h (t) | h ≤ C p,q .

Analysis of a family of time discretizations

In the sequel we study the following two-steps conservative time discretizations: θ-Scheme (TS ) which is the family of conservative centered Newmark schemes, see [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology -Volume 5 and 6 -Evolution Problems I and II[END_REF], such as the classical leap-frog scheme, the Stabilized Leap-Frog scheme (SLF ) as introduced in [START_REF] Gilbert | Higher order time stepping for second order hyperbolic problems and optimal CFL conditions[END_REF] and the higher order (θ, ϕ)-Scheme (TPS ) developed in [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]. Let ∆t > 0 be the time step, a small parameter devoted to tend to zero and let define t n := n∆t. To shorten the writing we introduce the following notations, for a series

{v k h } k≥0 ⊂ V h , [v n h ] ∆t 2 := v n+1 h -2v n h + v n-1 h ∆t 2 , {v n h } θ := θ v n+1 h + (1 -2θ)v n h + θ v n-1 h , n ≥ 1.
We seek a series {u k h } k≥0 ⊂ V h such that u n h is an approximation of u h (t n ) solution of equation ( 4). The two first terms of this series are considered given 1 . The following terms of the series are computed using one of the schemes below, for n ≥ 1,

• TS : [u n h ] ∆t 2 + A h {u n h } θ = f h (t n ). • SLF : [u n h ] ∆t 2 + A h u n h + ∆t 2 16 A 2 h u n h = f h (t n ). • TPS : [u n h ] ∆t 2 + A h {u n h } θ + ∆t 2 θ - 1 12 A 2 h {u n h } ϕ = f h (t n ) + ∆t 2 ∂ tt 12 + θ - 1 12 A h f h (t n ).
By construction (SLF) is an explicit scheme, whereas (TS) is explicit when θ = 0 and (TPS) is explicit when (θ, ϕ) = (0, 0). In this latter case, the scheme corresponds to the modified equation approach presented in [START_REF] Shubin | A modified equation approach to constructing fourth-order methods for acoustic wave propagation[END_REF]. Notice that for the (TPS) the source term has to be modified in order to get adequate accuracy.

All these schemes can be written

P K (∆t 2 A h ) [u n h ] ∆t 2 + P P (∆t 2 A h ) A h {u n h } 1/4 = f n h , (5)
with f n h the source term and with the following definitions of the polynomial functions P K and P P

• TS :

P K (x) = 1 + θ - 1 4
x, P P (x) = 1.

• SLF :

P K (x) = 1 - x 4 +
x 2 64 , P P (x) = 1 -x 16 .

• TPS :

P K (x) = 1 + θ - 1 4 x + ϕ - 1 4 θ - 1 12 x 2 , P P (x) = 1 + θ - 1 12
x.

In the sequel we aim at establishing energy and error estimates for general schemes of the form (5) under some assumptions on P K and P P . Let the time discretization error be defined as

e n h := u h (t n ) -u n h , it satisfies P K (∆t 2 A h ) [e n h ] ∆t 2 + P P (∆t 2 A h ) A h {e n h } 1/4 = r n h . (6)
Assuming sufficient regularity of the semi-discrete solution, there exist t n-1 ≤ t n, , t n,♥ , t n,♣ ≤ t n+1 such that for each of the considered schemes, the remainder writes

• TS : r n h = ∆t 2 1 12 ∂ 4 t u h (t n, ) + θA h ∂ tt u h (t n,♥ ) . • SLF : r n h = ∆t 2 1 12 ∂ 4 t u h (t n, ) + 1 16 A 2 h u h (t n ) .
• TPS :

r n h = ∆t 4 1 360 ∂ 6 t u h (t n, ) + θ 12 A h ∂ 4 t u h (t n,♥ ) + ϕ θ - 1 12 A 2 h ∂ tt u h (t n,♣ ) .
This suggests that (TS) and (SLF) are second order accurate and (TPS) is fourth order accurate in time.

For all these schemes, a discrete energy is preserved (in the absence of sources

It can be shown that the total energy satisfies

E n+1/2 h -E n-1/2 h ∆t = r n h , e n+1 h -e n-1 h 2∆t h . (9)
A necessary condition for the schemes' stability is the non-negativity of the kinetic energy [START_REF] Shubin | A modified equation approach to constructing fourth-order methods for acoustic wave propagation[END_REF] and the potential energy [START_REF] Joly | Variational methods for time-dependent wave propagation problems[END_REF]. Such a condition is achieved if the following hypothesis holds Hypothesis 2.1 There exists α > 0 such that P K (x) and P P (x) are non-negative for all x ∈ [0, α].

For the schemes we consider we have

• TS (0 ≤ θ < 1/4) : α = 4/(1 -4θ). • TS (θ ≥ 1/4) : α = +∞. • SLF : α = 16.
• TPS (θ = 0, ϕ = 0) : α = 12. • TPS (θ ≥ 1/4, ϕ ≥ 1/4) : α = +∞.

• TPS ∀ (θ, ϕ) : see theorem 5. where ρ(A h ) is the largest eigenvalue of A h . Note that when α is finite (resp. infinite) the scheme is conditionally (resp. unconditionally) stable. The non-negativity of the energy E n+1/2 h is however not sufficient to ensure the stability of the schemes in H since it is only a semi-norm. In the sequel, the polynomial functions P P and P K will be assumed to satisfy the following hypothesis Hypothesis 2.2 We assume that there exists a partitioning of [0, α] into two disjoint subsets J K and J P such that there exist two strictly positive constants C K , C P , such that

C K ≤ P K (x), ∀x ∈ J K , C P ≤ x P P (x), ∀x ∈ J P .
Définition 2.3 We define for any v h ∈ V h , the projection operators

Π K (v h ) = ∆t 2 λ∈J K λ∈Sp(A h ) (v h , w λ h ) h w λ h and Π P (v h ) = ∆t 2 λ∈J P λ∈Sp(A h ) (v h , w λ h ) h w λ h
The operator Π K (resp. Π P ) corresponds to modal projection operators on subsets of V h for which the kinetic (resp. potential) part of the energy is a uniform upper bound of the H-norm with respect to h and ∆t. These ideas are specified in the following Lemma. Lemma 2.2 Assume Hypothesis 2.2 holds. Then, for any h > 0 and ∆t > 0 satisfying the CFL condition (10), and for any

v h ∈ V h , v h = Π K (v h ) + Π P (v h ) and | Π K (v h ) | 2 h ≤ C -1 K P K (∆t 2 A h )v h , v h h , | Π P (v h ) | 2 h ≤ C -1 P ∆t 2 A h P P (∆t 2 A h )v h , v h h . (11)
Then the following estimates can be written Lemma 2.3 Let the series {e n h } satisfy (6) and Hypothesis 2.2 holds. Then for all h > 0 and ∆t > 0 satisfying the CFL condition (10), and for all n ≥ 1,

E n+1/2 h ≤ E 1/2 h + √ 2 γ ∆t n =1 r h h (12) e n+1 h h ≤ √ 2 e 1 h h + 2 γ t n 2 E 1/2 h + 4 γ 2 ∆t 2 n =1 k=1 r k h h , (13)
where

γ = C -1/2 P + C -1/2 K /2.
It is important to notify here that γ is independent of h and ∆t, this is the key point to obtain uniform estimates in ∆t. Indeed as pointed out in Remark 10 of [START_REF] Joly | Variational methods for time-dependent wave propagation problems[END_REF] the classical stability proof for the (TS) blows up when the time step approaches its greatest admissible value in inequality (10). The constant C K and C P (and therefore γ) are computed for the considered schemes in appendix, we find γ ≤ √ 2 for all schemes. Proof 2.4 (of Lemma 2.3) From Cauchy-Schwartz inequality applied to [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF] we get

E n+1/2 h -E n-1/2 h ∆t ≤ |r n h | h e n+1 h -e n-1 h 2∆t h ≤ |r n h | h Π K e n+1 h -e n-1 h 2∆t h Ξ + |r n h | h Π P e n+1 h -e n-1 h 2∆t h Φ ,
where Π K and Π P are the projection operators given by definition 2.3. Now we use the two inequalities of (11) on respectively Ξ and Φ and then recognize the square root of the kinetic and potential parts of the energy. To shorten the notation we introduce the time finite difference of the error:

d n+1/2 h := (e n+1 h -e n h )/∆t. We have 2 Ξ ≤ Π K (d n+1/2 h ) h + Π K (d n-1/2 h ) h ≤ C -1/2 K P K (∆t 2 A h ) d n+1/2 h , d n+1/2 h h + C -1/2 K P K (∆t 2 A h ) d n-1/2 h , d n-1/2 h h ≤ C -1/2 K 2 E n+1/2 h + 2 E n-1/2 h
, note that the second inequality is obtained by taking

v h = d n+1/2 h and v h = d n-1/2 h
in the first inequality of (11). A similar procedure leads to

Φ ≤ C -1/2 P 2 E n+1/2 h + 2 E n-1/2 h . By denoting γ = C -1/2 P + C -1/2 K /2 we get from E n+1/2 h -E n-1/2 h ∆t ≤ γ |r n h | h 2 E n+1/2 h + 2 E n-1/
Ω ≤ | Π K (e n h ) | h + Π K (d n+1/2 h ) h ≤ | Π K (e n h ) | h + ∆t C -1/2 K P K (∆t 2 A h )d n+1/2 h , d n+1/2 h h ≤ | Π K (e n h ) | h + ∆t C -1/2 K 2 E n+1/2 h
The same procedure leads to

Υ ≤ |Π P (e n h )| h + 2 ∆t C -1/2 P 2 E n+1/2 h .
We use recursively the former inequalities down to n = 1 to get,

e n+1 h h ≤ Π K (e 1 h ) h + Π P (e 1 h ) h ≤ √ 2 |e 1 h | h +2 √ 2 ∆t γ n =1 E +1/2 h .
We can reuse identity (12) to get (13).

Space time convergence result

Let ε n h be the error between the continuous solution and the fully discrete solution. Note that we have u(t n ) -u n h = e h (t n ) + e n h , thanks to this decomposition, assumption 1.1 and lemma 2.3 we can state our final result. Theorem 3.1 Assume that e 0 h = e 1 h = 0 and assume that Hypothesis 2.2 holds and that Hypothesis 1.2 holds for every sufficiently large (p, q), then, for all h > 0 and all ∆t satisfying the CFL condition (2.1) the following uniform convergence results hold

• θ-scheme (TS) sup t n ∈[0,T ] |u(t n ) -u n h | ≤ C H ∆t 2 γ 2 T 2 C 0,4 6 + 2 θ C 2,2 + δ(h) .
• Stabilized leap-frog scheme (SLF)

sup t n ∈[0,T ] |u(t n ) -u n h | ≤ C H ∆t 2 γ 2 T 2 C 0,4 6 + C 4,0 8 + δ(h) .
• (θ, ϕ)-scheme (TPS)

sup t n ∈[0,T ] |u(t n ) -u n h | ≤ C H ∆t 4 γ 2 T 2 C 0,6 180 + θ C 2,4 6 + 2 ϕ θ - 1 12 C 4,2 + δ(h) .
Note that, in particular these convergence results hold if ∆t 2 = α/ρ(A h ).

Appendix: Expression of γ for the considered numerical schemes

• TS (0 ≤ θ < 1/4). We are looking for a partitioning J K ∪ J P of the interval [0, 4 (1-4θ) ] and two positive constant C K and C P such that

C K ≤ P K (x) = 1 + θ - 1 4 x, ∀x ∈ J K , C P ≤ x P P (x) = x, ∀x ∈ J P
We see that the first inequality cannot be fulfilled for x = 4 (1-4θ) ∈ J K and the second one for x = 0 ∈ J P . The proposed partitioning is the following

J K = 0, a 2 , C K = 1 + θ - 1 4 a 2 , J P = a 2 , 1 (1/4 -θ) , C P = a 2
where a can be chosen according to θ in order to minimize γ = C

-1/2 P + C -1/2 K /2.
The peculiar optimal value is given in Lemma 2.3 of [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF], it reads

a(θ) = 4 (1 -4θ) 2/3 + (1 -4θ) ; γ(θ) = 1 a(θ) + 1 4 -(1 -4θ)a(θ) 2
Especially if θ = 0 we get γ = √ 2.

• TS (θ ≥ 1/4). In this case the polynomial functions P K (x) and P P (x) are positive on the whole interval [0, +∞] and moreover P K (x) ≥ 1 ∀x ≥ 0. Therefore one can choose J K = [0, +∞] and J P = {∅}. As a consequence one can set C K = 1 and since C P has no influence and can be chosen formally equal to +∞. We get γ = 1/2.

• SLF. We are looking for a partitioning J K ∪ J P of the interval [0, 16] and two constant C K and C P such that The choice of a that minimizes the value of γ is a = 1 • TPS (θ = 0, ϕ = 0). The approach is similar to the one presented above. We are looking for a partitioning J K ∪ J P of the interval [0, 12] and two constant C K and C P such that • TPS (θ ≥ 1/4, ϕ ≥ 1/4). In the case of the (θ, ϕ)-scheme with θ ≥ 1/4 and ϕ ≥ 1/4, the polynomial functions P K (x) and P P (x) read

C K ≤ 1 - x 8
C K ≤ 1 - x 4 + x 2 48 , ∀x ∈ J K , C P ≤ x 1 - x 12
P K (x) = 1 + θ - 1 4 x + ϕ - 1 4 θ - 1 12
x 2 , P P (x) = 1 + θ -1 12

x.

The polynomial P K (x) and P P (x) are positive on the whole interval [0, +∞] and P K (x) ≥ 1 ∀x ≥ 0. Therefore one can choose J K = [0, +∞] and J P = {∅}. We get γ = 1/2.

  Out of symmetry with respect to x = 8, we look for J K = [0, 8(1 -a)[ ∪ ]8(1 + a), 16] and J P = [8(1 -a), 8(1 + a)]. Then C K = a 2 and C P = 4(1 -a 2 ). Therefore

√ 2

 2 and leads to γ = √ 2.

1 √ 3

 13 , ∀x ∈ J P We look for J K = [0, 6(1 -a)[ ∪ ]6(1 + a), 12] and J P = [6(1 -a), 6(1 + a)], this choice gives some values for C K and C P that can be optimized to minimize γ. The optimal choice is a = and γ = √ 2.

  ). This energy is computed at the intermediate time step t n+1/2 and is denoted E

							n+1/2 h	. It is the sum of a kinetic and potential energies,
	more precisely E h n+1/2	= E K n+1/2	+ E P n+1/2	. When considering equation (6) satisfied by the error, we have
			E K n+1/2	=	1 2	P K ∆t 2 A h	e n+1 h ∆t -e n h	,	e n+1 h ∆t -e n h	h	(7)
	and		E P n+1/2	=		1 2	P P (∆t 2 A h ) A h	e n+1 h 2 + e n h	,	e n+1 h 2 + e n h

h

.

  1 of [6].

	Lemma 2.1 Assume Hypothesis 2.1 holds. The energies E K n+1/2	and E	n+1/2 P	are non-negative for all
	n ≥ 0 if ∆t satisfies	∆t 2 ≤	α ρ(A h )	(CFL condition)	(10)

In practice, u 0 h and u 1 h are computed using the initial data and a sufficiently accurate one-step method (such as high-order Runge Kutta methods).


Jean-Baptiste Courbot
email: jb.courbot@unistra.fr

Emmanuel Monfrini

Vincent Mazet

Christophe Collet

ORIENTED TRIPLET MARKOV FIELDS FOR HYPERSPECTRAL IMAGE SEGMENTATION

Keywords: Triplet Markov Field, Bayesian Segmentation, Orientation Retrieving

Hyperspectral image processing benefits greatly from using spatial information. Markov field modeling is a well-known statistical model class for considering spatial relationships between sites of an image. Often, the model restricts to Hidden Markov Field, therefore cannot handle non-stationarities in the images. This paper presents a Triplet Markov Field model for hyperspectral image segmentation, allowing the joint retrieving of image classes and local orientations. Segmentation results on synthetic data validate the methods, and results on real astronomical data are presented.

INTRODUCTION

Cosmological theory predicts the existence of a "cosmic web", that is weak gas concentration, linking galaxies through elongated structures [START_REF] Van De Weygaert | The cosmic web: geometric analysis[END_REF]. These structures are assumed to trace galaxy formation. The Multi-Unit Spectroscopic Explorer (MUSE) instrument is designed to reveal such features in deep field hyperspectral images (HSI) [START_REF] Bacon | The MUSE 3D view of the Hubble Deep Field South[END_REF]. The purpose of this paper is to study the relevance of Triplet Markov Fields (TMF) [START_REF] Benboudjema | Unsupervised image segmentation using triplet Markov fields[END_REF] for segmenting oriented, faint structures in HSI. Indeed, modeling local orientation of elongated structures can be performed within TMF and could improve the segmentation.

Priors used for HSI segmentation are generally either spectral or spatial. For instance, spectral priors are used in supervised anomaly detection [START_REF] Matteoli | A tutorial overview of anomaly detection in hyperspectral images[END_REF]. On the other hand, spatial priors can translate into object-wise processing [START_REF] Meillier | Nonparametric Bayesian extraction of object configurations in massive data[END_REF], morphologybased processing [START_REF] Fauvel | Advances in spectralspatial classification of hyperspectral images[END_REF] or local priors, as in the case of Markov Fields. In this framework, numerous works were published for HSI processing, most of them using Hidden Markov Field (HMF) models [START_REF] Susan | Hyperspectral imagery: Clutter adaptation in anomaly detection[END_REF][START_REF] Rellier | Texture feature analysis using a Gauss-Markov model in hyperspectral image classification[END_REF][START_REF] Li | Spectralspatial hyperspectral image segmentation using subspace multinomial logistic regression and Markov random fields[END_REF][START_REF] Eches | Adaptive Markov random fields for joint unmixing and segmentation of hyperspectral images[END_REF][START_REF] Li | Hyperspectral image classification using Gaussian mixture models and Markov random fields[END_REF][START_REF] Xia | Spectral-spatial classification for hyperspectral data using rotation forests with local feature extraction and Markov random fields[END_REF]. Despite its common use, standard HMF modeling has significant lacks: especially, spatial noise correlation and non-stationarity are not modeled. Generalizing HMF, pairwise [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] and triplet [START_REF] Benboudjema | Unsupervised image segmentation using triplet Markov fields[END_REF] Markov field models handle these two aspects.

In addition, spatial priors can enhance directions, meaning that oriented object are searched for. To our knowledge, there exist only a few works reporting the search for orientations in HSI. Most of them are related to edge detection techniques, targeting for instance road network detection in remote sensing [START_REF] Kamandar | Linear feature extraction for hyperspectral images based on information theoretic learning[END_REF]. Varying local orientations in images can be seen as varying stationarities of the underlying stochastic processes. Hence, it is desirable to provide an accurate modeling in the context of HSI statistical segmentation.

This paper presents a TMF model for segmentation in HSI, recovering jointly classes of interest and local orientations in the image. The main objective is to enhance the segmentation of oriented structures in HSI. The novelties of the approach are both the use of TMF and the estimation of orientations in HSI.

This paper uses the following notations. S is the set of sites in the image. A set of random variables over the image is noted A = (A s ) s∈S and the corresponding realization set is a = (a s ) s∈S . We will note A Ns the set of random variables over the local neighborhood N s of the site s, and A S\s the random variable set over the image, except in the site s. In the remainder of this paper, Section 2 introduces the proposed model and corresponding methods, and Section 3 presents some experimental validations on HSI.

MODEL

Oriented Triplet Markov Fields

Triplet Markov Fields (TMF) [START_REF] Benboudjema | Unsupervised image segmentation using triplet Markov fields[END_REF] are Markov fields with three constituting processes: T = (Y, X, V). The process Y represents the observed HSI, containing the spectra Y s ∈ R Λ , with Λ bands per site s. The class process X takes its value from a finite set Ω x = {ω 1 , . . . ω K }. Similarly, the auxiliary process V takes values from Ω v = {ν 1 , . . . , ν L }.

TMF, while generalizing HMF modeling, still benefit from posterior Markovianity: (X, V) conditionally to Y is a Markov field [START_REF] Benboudjema | Unsupervised image segmentation using triplet Markov fields[END_REF]. This ensures that statistical segmentation is possible, as in HMF. The model verifies:

p(Y, X, V) ∝ p(Y|X, V)p(X, V) ∝ p(Y|X, V)p(X|V)p(V).
For Maximum Posterior Mode (MPM) segmentation [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF], posterior probabilities must be computed on each site. In our model, V is a Markov field. Besides, p(X|V) is a Markov field distribution, in which V rules X in each site independently.

Then, for each site s on the lattice S:

p(X s , V s |Y, X S\s , V S\s ) ∝ p(Y|X, V)p(X s |X Ns , V)p(V s |V Ns ) ∝ p(Y|X, V)p(X s |X Ns , V s )p(V s |V Ns ). (1)
Conditionally to X, Y and V are assumed independent. Furthermore, we assume that the Y s are independent conditionally to X. Finally, each Y s depends only on the corresponding X s :

p(Y|X, V) = p(Y|X) = s∈S p(Y s |X) = s∈S p(Y s |X s ).
V relies on a Ising/Potts potential (as used in [START_REF] Derin | Modeling and segmentation of noisy and textured images using Gibbs random fields[END_REF]), ruled by a positive, real-valued parameter α v :

∀s : p(V s |V Ns ) ∝ exp α v s ∈Ns 1 -2δ(V s , V s ) . (2)
We restrict to binary cliques, covering the site s and a neighboring site s within a 8-neighborhood.

Similarly, the distribution p(X s |X Ns , V s ), parametrized by α x|v ∈ R + , is expressed as:

∀s : p(X s |X Ns , V s ) ∝ exp α x|v s ∈Ns φ s (V s ) (1 -2δ(X s , X s )) ; (3)
where φ s is an orientation function, defining the clique potential weights used for X in the case of Oriented TMF (OTMF) Algorithm 1 SEM algorithm for OTMF Input: Observed realization y, stop criterion Output: Parameter estimation Θ SEM Initialization: Θ noise (0) , Θ field (0) while stop criterion not reached -iteration q do Simulate (x q , v q ) from p Θ q-1 (X, V|y). Estimate Θ noise (q) through MLE on (x q , y). Estimate Θ field (q) with DE, with (x q , v q ). end while accounting for the favored direction described in V s . The orientation function is:

φ s (V s ) =        | cos(V s )| if s ∈ {3, 7} | cos (V s -π/4) | if s ∈ {4, 0} | cos (V s -π/2) | if s ∈ {5, 1} | cos (V s -3π/4) | if s ∈ {6, 2} ;
where s denotes the index within the neighborhood N s (see Fig. 1).

To summarize, the distributions expressed in Eqs ( 2) and ( 3), together with the knowledge of the noise behavior, enable the Bayesian segmentation according to the MPM criterion [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF]. ∀s ∈ S:

(x s , vs ) MPM = arg max (ω,ν)∈Ωx×Ωv p(X s = ω, V s = v|y). (4)
There is no exact analytical expression for these estimations, however simulations are possible through Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]. Then, MPM estimates amount to choose for each site the most frequent Gibbs sampler realization.

Parameter Estimation

When performing unsupervised segmentation, the model parameters must be estimated. The parameter set Θ splits into field parameters, Θ fields = {α v , α x|v } and noise parameters, Θ noise . From a realization (x, v), the parameter set Θ fields is estimated by the method from Derin and Elliott (DE) [START_REF] Derin | Modeling and segmentation of noisy and textured images using Gibbs random fields[END_REF]. From complete data (y, x, v), Θ noise is obtained through the Maximum Likelihood Estimator (MLE). Complete data are however unavailable: this is why we use the Stochastic Expectation-Maximization (SEM) algorithm [START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF], whose principle is to mimic complete realization with iterative simulations. Let us note that it has been shown in [START_REF] José | An empirical comparison of EM, SEM and MCMC performance for problematic Gaussian mixture likelihoods[END_REF] [START_REF] Monfrini | Estimation de mélanges généralisés dans les arbres de Markov cachés, application à la segmentation des images de cartons d'orgue de barbarie[END_REF] that SEM provides in practice slightly better results than the Expectation-Maximization (EM) [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] or the Iterative Conditional Estimator (ICE) [START_REF] Delignon | Estimation of generalized mixtures and its application in image segmentation[END_REF]. The SEM algorithm is reported in Alg. 1.

The stop criterion is a similarity measure between consecutive estimations, to ensure that convergence is reached. Initial parameters are estimated through simple estimates: Kmeans clustering [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] provides an initial X, and the initial V is obtained through image gradient measure on the initial X.

EXPERIMENTS

Observation model

Detection aims at segmenting significant signal versus noise only. Therefore, Ω x contains two labels ω 0 = 0 and ω 1 = 1, denoting signal absence and presence, respectively. Furthermore, the spectral noise behavior is assumed identical with respect to the considered band, and correlated over three adjacent bands. The spectral covariance matrix is written as:

Σ = σ 2 ρ 1 ρ 2 0 0 ρ 1 ρ 2 0 0 ρ 2 ρ 1 0 0 ρ 2 ρ 1 σ 2                        
Assuming a zero-mean, additive Gaussian spectral noise implies that Θ noise = {µ, σ, ρ 1 , ρ 2 }. µ is the spectrum to detect, σ stands for the noise standard deviation, and ρ 1 , ρ 2 are the spectral noise correlation parameters. In this paper, we use two orthogonal directions for the auxiliary process: Ω v = {π/4, 3π/4} ; to ensure a proper discrimination between orientations. Note that we restrict to a simple instance, but that generalizing to more subdivisions is straightforward.

Images Three images are used for experiments:

-an image generated through OTMF; -a synthetic image generated conditionally to a fixed auxiliary process V, which is not a Markov field realization; -a real HSI from the MUSE instrument [START_REF] Bacon | The MUSE 3D view of the Hubble Deep Field South[END_REF], obviously without V nor X ground truths for such data. The synthetic images, of size 128 × 128 px, are presented in Figs. 2a and2b. We use a weak average signal spectrum, presenting a few non-zero coefficient (see Fig. 2c) over Λ = 20 spectral bands. Let us remark that we consider that the images are hyperspectral in the sense that they correspond to very fine and close bandwidths. Two cases are investigated: SNR = -5 dB and SNR = -10 dB. The SNR is defined as:

SNR = 10log 10 µ 2 2 /Tr(Σ) = 10log 10 µ 2 2 /Λσ 2 .
Finally, the correlation coefficients were set to ρ 1 = 0.5σ 2 and ρ 2 = 0.25σ 2 .

The MUSE HSI has 90×55 px over Λ = 50 spectral bands, and is illustrated in Fig. 3. Segmentation on this image should reveal faint, spectrally-sparse and spatially wide light emission (called the Lyman-alpha line) located in the vicinity of distant galaxies. In such deep field observation, these galaxies occupy no more than 5 to 10 pixels. Detection of an oriented structure in this image may be a clue for a filamentary inter-galactic structure. Fig. 2. Synthetic images, spectra and the correponding segmentation results.

Results

Results on parameter estimation and segmentation error rate are reported in Table 1. Segmentation results are also shown in Fig. 2. The segmentation results are compared with a classical HMF modeling, using also SEM for parameter estimation and MPM for segmentation. The estimated values for noise are very close to real parameters, with less than 1% difference in most cases. On the other hand, field parameters α x|v and α v are poorly estimated, with a strong variance (not shown). One can notice that this misestimation does not affect significantly segmentation result: in this model, their low influence on the posterior distribution from Eq. (1) makes it difficult to provide precise estimations (see e.g [START_REF] Salzenstein | Parameter estimation in hidden fuzzy Markov random fields and image segmentation[END_REF] for an instance of similar conclusions). ), when processing oriented, thin structures. Indeed, the isotropic potential used in HMF may make the corresponding segmentation fail in such cases. On the other hand, the MPM estimation on V (error rate: E OTMF v ) is not as efficient as on X. This is mainly because some cases (e.g. wide, homogeneous region) make the direction estimation hard to perform. Therefore, there is an intrinsic imprecision on the segmentation result. Let us also remark that the segmentation results are better on the first synthetic image than on the second, because the fixed V in the second case is not a Markov field realization.

Estimated parameters for the real MUSE HSI are also reported in Table 1, and the estimated signal spectrum is displayed in Fig. 4a. The SNR, estimated at -5.04 dB, is close to the SNR used in synthetic cases, which induces to believe that parameter estimation are of similar quality. One can especially notice that the correlation coefficients are significantly far from zero. This is due to MUSE data pre-processing, and hence such correlation must be considered when handling MUSE observations.

The segmentation results are reported in Fig. 4b. Ground truths are naturally not available. One can notice, however, that the detected region xMPM is significantly larger than what is expected from a single galaxy, and that the wide detected emission cannot be due to the instrument point spread function alone, because of its asymmetry. Concerning the estimated orientations vMPM , the result is difficult to interpret: we notice, however, that the border of the detected region and the corresponding orientations are coherent. Let us finally remark that spatial noise correlation could induce mis-estimations, because the proposed model assigns such correlation on the hidden processes X and V.

CONCLUSION

This paper presented a model for joint Bayesian segmentation of classes and orientations in an HSI. The method was validated in the case of signal detection on synthetic images, outperforming HMF, and applied to real astronomical HSI. Further work will include models and algorithms refinement, to handle more complex cases.

ACKNOWLEDGEMENT

The authors would like to thank Roland Bacon (France, CNRS-CRAL) for providing the MUSE data.

Fig. 1 .

 1 Fig. 1. Binary clique indexing and local orientation measurment within a local neighborhood.

  (a) First synthetic image, at SNR = -5 dB.

  Second synthetic image, at SNR = -10 dB.

  Instances of noisy synthetic spectra (dashed) and µ (plain).

Fig. 3 .

 3 Fig. 3. Astronomical HSI from MUSE (Im3).

  Estimated average spectrum µ (plain) and detected spectrum instance (dashed).

  xMPM vMPM (b) MPM estimations for the X and V processes.

Fig. 4 .

 4 Fig. 4. Estimation results on the MUSE HSI.

Table 1 .

 1 Real parameters, SEM estimations and segmentation error rates.

				α x|v	α v	σ 2	ρ1	ρ2	µ (MSE) E HMF x	E OTMF x	E OTMF v
	Im1	-5 dB -10 dB	Real Est. Real Est.	5.0 0.8194 0.8097 0.2131 0.1065 0.0532 20.0 0.2119 0.1059 0.0530 5.0 20.0 0.6700 0.3350 0.1675 0.6543 0.9304 0.6683 0.3335 0.1657	2.10 -4 9.10 -4	0.29% 3.20%	0.27% 0.78%	10.50% 14.09%
	Im2	-5 dB -10 dB	Real Est. Real Est.	8.0 0.9357 0.8641 0.2114 0.1058 0.0525 -0.2119 0.1059 0.0530 8.0 -0.6700 0.3350 0.1675 10.10 -4 5.10 -4 0.7493 0.8843 0.6699 0.3351 0.1672	0.40% 3.67%	0.27% 0.95%	17.41% 22.31%
	Im3	-	Est.	1.2527 0.5878 0.5856 0.1227 0.0154	-	-	-	-
	y spectral average			y at λ = 10
	y at λ = 25			y at λ = 45

339659-MUSICOS. This work was funded in part by the DSIM project under grant ANR-14-CE27-0005.

acknowledges support from the ERC advanced grant


Daian Yue

Vania Joloboff
email: vania.joloboff@inria.fr

Frédéric Mallet

Flexible Runtime Verification Based On Logical Clock Constraints

Keywords: Runtime verification, trace analysis, property specification, clock constraint specification language, simulation, debugging, model driven engineering, virtual prototyping

We present in this paper a method and tool for the verification of causal and temporal properties of embedded systems, by analyzing the trace streams resulting from virtual prototypes that combines simulated hardware and embedded software. The proposed method makes it possible to analyze different kinds of properties without rebuilding the simulation models. Logical clocks are used to identify relevant points to put observation probes and thus also reducing the trace streams size. We propose a property specification language, called PSML, and based on behavioral patterns that does not require knowledge of temporal logics. From a given PSML specification, simulation is instrumented to generate a trace and the code is dynamically loaded by the simulator. The resulting trace stream is analyzed by parallel automata generated from the specification. The experiments, developed over the SimSoC virtual prototyping framework, show flexibility, possibility of using multi-core platforms to parallelize simulation and verification, providing fast results.

I. INTRODUCTION

Virtual prototypes provide good solutions to model and verify systems combining hardware components with embedded software, hardware-in-the-loop, real-time operating systems, various devices and their drivers. Even though exhaustive verification techniques can address increasingly larger parts of the system, integration requires complementary techniques that consider the system as a whole. Virtual prototypes also provide a good framework to test alternative hardware architectures or develop the embedded software while the hardware architecture is not yet fully available.

We promote the use of a Model-based verification approach. The model provides a high-level abstraction to detect early problems, it contains a set of expected properties from the different components of the system. From this central model, several artifacts can be generated, including, in our case, verification modules to ensure that a component implementation satisfies the expected system properties. For known required properties, assertion-based verification techniques can be used [START_REF] Dahan | Combining system level modeling with assertion based verification[END_REF], [START_REF] Foster | Assertion-based verification: Industry myths to realities[END_REF], typically by running many simulations with varying parameters and varying interactions.

However, most often during the design phase the engineers discover a system failure that is not captured by the assertions as it was not anticipated. Investigating the cause of the failure may not be easy because it is a long chain of causality events that eventually lead to failure. To investigate the failure cause, simulation trace files are generated that record a lot of information which can be exploited with trace analysis tools. A trace stream usually consists of the sequence of mixed events, states and transition steps with variable values, annotated with time stamps. The engineers eventually discover the problem by successive iterations, typically searching in the trace stream at which point it deviates from the model. They also often discover new properties that should be verified but were not formulated as initial requirements, which can be added as new assertions. Trace analysis is thus very useful as long as the trace analysis tool is convenient enough to be handled by seasoned engineers and powerful enough to discover intricate properties.

Two issues in trace analysis are that firstly it may be impractical or not feasible to re-compile and re-build the simulator to check new properties; secondly that trace files have become huge in industrial applications (hundreds of Gigabytes) as systems have become fairly complex, although it is often the case that a lot of the trace stream data is actually irrelevant to the failure cause.

The purpose of this work is to make it easy to check for properties by quickly analyzing trace streams with a powerful tool able to verify causal and temporal properties, also able to reduce the size of trace streams. We are proposing a new approach for trace analysis that does not require rebuilding the simulator when adding new properties. Raw binary data are abstracted into a formal trace that can be used either dynamically as a run time verification tool to detect property violations (without slowing down the simulation) or statically by filtering out unnecessary data with respect to the properties.

The paper is structured as follows. The next section presents background work in the trace analysis and property specification areas. Section III presents our global approach, decomposed into several steps, based on a modeling approach and model transformation using Domain Specific Languages (DSL). Section IV presents the details of the tools function and an example drawn from a System on Chip simulator.

II. BACKGROUND AND RELATED WORK A. Related Work

Even though model-driven engineering approaches with model verification techniques are used, the engineers still face the problem of debugging simulation models [START_REF] Krahl | Debugging simulation models[END_REF] and runtime verification is still needed in most industrial applications [START_REF] Leucker | A brief account of runtime verification[END_REF]. Trace analysis has been a topic of work for some time in the simulation community. Several languages in the literature allow to express temporal properties, among which Linear Temporal Logics [START_REF] Pnueli | The Temporal Logic of Programs[END_REF]. LTL abstracts the system behavior as an infinite sequence of states and establishes either properties of invariance or eventuality. All the properties are expressed relative to the steps of execution, using modal temporal operators, and not relative to an external (physical or not) clock measure. Several efforts have been made to provide real-time extensions for temporal logics (RT-LTL and others).

Later, some efforts were made to renounce to the full expressiveness of some logic in favor of a representation more natural to the designers and more computationally tractable. This lead to the notion of Logic Of Constraints (LOC) initially proposed by Balarin and al. [START_REF] Balarin | Constraints specification at higher levels of abstraction[END_REF] and used for trace analysis by Chen and al. [START_REF] Chen | Automatic trace analysis for logic of constraints[END_REF]. This system allows to express functional and performance constraints containing event names, instances of events, index variables allowing to express generic formulas. For example, the formula t(Stimuli[i + 1]) -t(Stimuli[i]) < T denotes the property that two successive occurrences of the event Stimuli must occur within a given time frame T .

As SystemC is widely used by the embedded systems community, the specific VCD (Value Change Dump) trace format has been defined for signal tracing and many commercial or non commercial tools exist for tracing wave forms. A more flexible trace generator is proposed by CULT [START_REF] Hong | Cult: A unified framework for tracing and logging c-based designs[END_REF], but there is no associated tool to analyze the traces produced. A SystemC analyzer like [START_REF] Tabakov | A temporal language for systemc[END_REF] can be used to verify temporal properties holding in the SystemC processes, but it requires access to the SystemC modules and cannot check information originating from embedded software running over the simulated hardware. Although the experiments related in this article are based on SystemC, nothing in the approach requires SystemC, and our DSL makes it possible to check value changes on any variable, not restricted to signals.

The COSITA tool [START_REF] Khlif | Co-simulation trace analysis (cosita) tool for vehicle electronic architecture diagnosability analysis[END_REF] analyzes traces resulting from co-simulation between SystemC and Matlab models for automotive applications. The tool permits to notice differences between simulation and emulation to investigate the models but it does not include a property specification language. The Meta-Event Definition Language (MEDL) [START_REF] Bhargavan | Verisim: Formal Analysis of Network Simulations[END_REF] is used to verify properties in traces of packets in network simulation, using the Monitoring and Checking (MAC) framework [START_REF] Kim | Formally Specified Monitoring of Temporal Properties[END_REF]. It is based on a logic for events and conditions, which allows to consider instances of events by means of counters, similar to LOC.

Regarding property specification, the seminal paper from Dwyer and al [START_REF] Dwyer | Patterns in property specifications for finite-state verification[END_REF] has inspired many subsequent works. In particular the pattern specification language proposed by [START_REF] Konrad | Real-time specification patterns[END_REF] can capture causality and conditional properties with timing conditions, but does not seem to capture counting occurrences.

A similar though different approach from trace analysis consists in considering usage scenarios and either generate timed automata that can be model-checked such as [START_REF] Some | From scenarios to timed automata: Building specifications from users requirements[END_REF] or generate tests associated with verification of pre and post conditions [START_REF] Tsai | Rapid embedded system testing using verification patterns[END_REF]. The work from [START_REF] Tokarnia | Scenario patterns and trace-based temporal verification of reactive embedded systems[END_REF] has extended this mechanism to actually implement trace analysis to verify some temporal relations. Yet another approach is to replay execution of a trace over the model to investigate the example at stake [START_REF] Hegedüs | Replaying execution trace models for dynamic modeling languages[END_REF].

Regarding the tooling aspect, the Open Trace Format [START_REF] Eschweiler | Open trace format 2: The next generation of scalable trace formats and support libraries[END_REF] provides a flexible trace format that inspired our trace mapping work. Our mapping system is using the notion of meta model for defining a trace format introduced by [START_REF] Hamou-Lhadj | A metamodel for the compact but lossless exchange of execution traces[END_REF]. Our specification language is based on the Clock Constraint Specification Language (CCSL) [START_REF] Mallet | CCSL: specifying clock constraints with UML/Marte[END_REF]. CCSL is uncomparable to LTL and CTL (see [START_REF] Gascon | Logical time and temporal logics: Comparing UML MARTE/CCSL and PSL[END_REF]), it provides a set of safety patterns regarding both causal and temporal properties. It combines purely boolean expressions with unbounded counters, thus being more expressive than regular languages and adequate for properties targeted here. Although our tool does not use the software library from TimeSquare [START_REF] Deantoni | Timesquare: Treat your models with logical time[END_REF], it relies on the notion of logical clocks. In the sequel we use the word tick to mean an occurrence of a logical clock.

B. The SimSoC Framework

The work described here is based and integrated into SimSoC [START_REF] Helmstetter | SimSoC: A SystemC TLM integrated ISS for full system simulation[END_REF], an open-source simulation framework developed at Inria [START_REF]SimSoC open source software[END_REF]. The code generated by the tools described below is compiled and linked with the SimSoC Instruction Set Simulator to generate traces. Although our work is based on SimSoC, it is independent and can be re-targeted to any other simulator.

III. RUNTIME VERIFICATION A. Global Architecture

The global software architecture of our approach is shown on Figure 1. It is decomposed into four independent steps. The first step is the trace generation, based on a trace meta model so that all traces are conforming to some well defined model. The second step is a model transformation process, using a DSL called STML, in order to map the actual simulation binary trace data into an abstract trace of logical clocks. The third step is the property specification by means of another DSL called TPSL, an English-like pattern-oriented specification language. The final step is the property verification over the mapped trace, using code generated by the property specification compiler.

The individual steps, described in more details in the following sections, are motivated and structured as follows:

• We initially place a light constraint on the simulator: The simulator can generate arbitrary binary data into what is called the raw trace as long as each data item dumped into the trace can be modeled as a structure that can be described itself using the meta-model that is defined and explained below in sectionIII-B. Because the trace corresponds to such a model, it can be parsed and understood for extraction and transformation. Another constraint on the simulator implementation is that it must comply with the trace interface that we provide to produce trace items. • The runtime verification is based on the logical clock notion from CCSL. The verification process eventually operates on logical clocks. Thanks to the meta-model, a mapping process can be defined, which maps raw binary data into logical clock ticks. For example the fact that the variable value of a sensor has changed can be transformed into a specific clock tick.

• This mapping of raw data is achieved using a domain specific language (DSL) called STML (Simulation Trace Mapping Language), defined using the Eclipse Modeling Framework. The engineers can take advantage of STML to filter out useless data from the trace and significantly reduce the trace file size, or map complex data structures into clock ticks.

The STML compiler generates the mapping code as a dynamically loadable module implementing the trace interface. • After the raw trace has been mapped to a trace of logical clocks, it can be analyzed to verify the required properties, expressed in another DSL named TPSL (Trace Property Specification Language). TPSL has been designed to be a simple language usable by engineers that do not require knowledge of temporal logics or logical clock algebra. It has English keywords and simple common mathematical expressions. A TPSL module consists of an arbitrary number of properties.

The TPSL compiler compiles each module into a set of parallel automata corresponding to the properties. The resulting automata are applied to the mapped trace file and result into a PASS or FAIL status. In case of a FAIL an explanation is provided upon the failed expression in TPSL. Since the automata are parallel and no backtrack is ever operated, the verifier performance is linear with the trace file size.

A key notion in CCSL is the notion of simultaneity. In this framework, it is a parameter of the system. The simulator must provide a timestamp for all data item produced. These timestamps must be expressed in some unique real time clock, typically the hardware clock, the resolution of which is arbitrary. For verification of properties, a trace clock is used. The mapping transformer is parameterized by the trace clock frequency, which must be smaller than the real clock frequency. The trace emitter finally emits logical clock ticks annotated with tags related to the trace sampling frequency. Therefore, trace items that may have a different timestamp relative to the real time clock may translate into simultaneous clock ticks as observed by the trace clock, as they carry the same time stamp. The trace clock frequency should be chosen by the engineers depending on their use cases.

B. Trace Emission

To avoid recompiling and rebuilding the simulator when checking different properties, the trace generation is done inside a dynamically loaded module that is provided when the simulation starts as a parameter. The simulation must call this module using the provided interface. A key element of this interface is the trace_emit() function, an abstract polymorphic function that the simulation engineers must call from their model.

The interpretation of the raw trace data is defined by a Ecore/UML model. This model must comply with a meta-model that is essentially the following: each data item in the trace must carry a timestamp and belong to some class. Each class has a name and an arbitrary number of attributes. Each attribute has a name and a data type, which must be a basic value type, pointers and references are not supported. Embedded structures are not supported either at the moment, except fixed size arrays. These data items can represent anything, sensor values, states, transitions, events, transactions, no specific requirement is made.

A trace model following the meta model must be provided by the engineers constructing the simulator (using Eclipse Modeling Framework) who can output any piece of information described in the model. In our experiment, the simulator is the SimSoC simulator from Inria coded in SystemC, hence using C++ code. From this model, it is possible to derive trace generation code that captures and generates the simulation data into a well understood format, such as XML technology. For example, the trace model can describe a SystemC TLM transaction to have an initiator, a target, an operation, an address and argument data. The data types of those attributes are described using abstract type names. The code generator is parameterized to map these abstract types into the target programming language data types, also provided in the trace model.

A call to function trace_emit(transaction) generates a transaction data item in the trace whoseformat is the one derived from the trace model. As the trace generator is using code generated from the model, it can arrange the layout of the trace data as most convenient for later processing, while maintaining the integrity. For example the N th bit of a register can be translated into a boolean flag if deemed appropriate. The role of function trace_emit() is to decide how the trace data should be emitted.

C. Trace Mapping and Filtering

The next step in the process consists in mapping the trace events to logical clocks such as defined by CCSL [START_REF] Mallet | CCSL: specifying clock constraints with UML/Marte[END_REF]. The STML was devised to define such a mapping. A STML module takes as input two parameters, on one hand the trace model described above, on the other hand a set of clock definitions. This clock specification actually comes from the TPSL language described below. The mapping module then contains instructions that specify the mapping of the trace data into either nothing (the data is eliminated from the trace) or into one or more clock ticks.

A STML module consists of a set of rules, each one introduced by the keyword for. There must be one and exactly one rule for each trace item class described in the trace model passed as the first parameter. Each rule describes how the trace items of this class translate (or not) into the logical clock occurrences of the second parameter, using a syntax close to the C language for expressions.

In the sequel of this paper we illustrate the approach with a virtual prototype example that consists of a simulator of a platform with a multi-core Power architecture processor. Each core has its own MMU (Memory Management Unit), and is connected to various peripherals, including a OpenPIC multi-core interrupt controller and a transmit controller with a FIFO queue. The hardware is controlled by a real time software that has drivers monitoring the peripherals. The hardware FIFO queue receives items asynchronously from the software but it re-sends these items synchronously to the outside world. The hardware should (i) emit a F if oEmpty signal whenever the queue is empty (including at initialization phase), a F if oF ull signal when it is full, (ii) send the data items with limited jittering (iii) have a reasonable latency. On the other hand, the software should (i) enable interrupts and service them, and (ii) it should not put data into the queue when it is full.

In STML example on Fig. 2, only data items related to sending and receiving and related interrupts to the FIFO are mapped in order to analyze the FIFO handling behavior. All other trace data items are simply thrown away. Only the items related to the FIFO are mapped into clock ticks. Regarding interrupts, the trace items emitted by the interrupt controller translate into different clock ticks depending on the successive trace items as the interrupt controller is tracing both whether it receives a signal and whether the corresponding interrupt is enabled. In the latter case, the interrupt may stay pending (e.g. because there is a higher priority interrupt) or it is raised to some core. These different situations map to different clock ticks.

The STML language helps reducing the trace file size by either filtering out irrelevant data with regards to the properties to be verified or by combining complex binary data into a single clock tick. Another feature reducing the trace and facilitating property verification is the ability to emit clock ticks only when a value changes, eliminating the need to repeatedly store the same value. The STML compiler generates a dynamic linking library, loaded dynamically by the simulator. This dynamic library is called by trace_emit() interface and implements the specified mapping. Therefore one can dynamically generate as many traces as desired from the same simulator, without rebuilding the virtual prototype (which may not be an easy and quick task). Conversely, using static analysis, a huge file can be dumped using no filtering, but the library can be stored and reused to transform that huge file into as many clock mappings as desired. Any raw trace can be reduced to a stream of logical clock ticks. Since the STML compiler requires to handle all possible classes specified in the trace model, the user must specify for each trace item whether it translates into clock ticks or whether it is ignored. Users may have to specify a fairly large number of clocks to distinguish specific cases. When specifying TPSL properties as shown below, it is possible to either specify specific properties for specific clocks or to group them using the union operator.

D. The Property Specification Language

Following the same reasoning as Balarin et al. [START_REF] Balarin | Constraints specification at higher levels of abstraction[END_REF] or Drechler et al. [START_REF] Drechsler | Formal specification level: Towards verification-driven design based on natural language processing[END_REF], we advocate for a specification language close to the domain and easy to compute rather than as expressive as temporal logics but failing to get a wide adoption among engineers. Another domain-specific language, named TPSL, is dedicated to expressing properties on simulation traces. TPSL 1 also uses the Eclipse Modeling Framework for easy integration in development tools so that engineers can express properties that should be verified on traces. This language first defines the list of logical clocks (that are exported to be generated by the mapping language). We have captured in this language the properties that cover a large part of the CPS engineering applications:

• causality: the fact that when some event occurs or some state is reached, there should be a consequence in the future.

• delay: the fact that two events should not be spaced by more or less than a given (time) value.

• counting occurrences of events and comparing counts.

• specific well known properties such as jitter, latency and drift. The property specification language in its current version focuses mostly on exploring causality factors. It is our experience that many failures occur because some complex causality fails. Verifying a TPSL property boils down to checking expressions in the CCSL clock algebra. In this paper, we do not re-introduce the formal definition of CCSL algebra. The reader is referred to [START_REF] André | Syntax and Semantics of the Clock Constraint Specification Language (CCSL)[END_REF] to obtain the formal definitions. TPSL so far only uses a subset of these clock operators, namely:

• the union (noted + in CCSL) of two clocks, that ticks whenever one of them ticks.

• the intersection (noted * in CCSL) of two clocks, that ticks whenever both of them ticks simultaneously (as defined above in our context). • sampling (noted in CCSL) of a clock A by another clock B, which ticks whenever the second clock B ticks after some clock tick of A. • subclocking (noted ⊂ in CCSL) of clock by some frequency divider. A TPSL module must first declare the clock identifiers and the constant values that are used in subsequent properties. These declared clocks are called the basic clocks. The users may define additional new clocks from the basic ones using one of the following operators:

• X = A or B to define X as the union of the two clocks A and B.

• X = A and B to define X as the intersection.

• X = A by N to define X as the CCSL subclock of A corresponding to frequency divider N. After the clock definitions the users may define properties related to these clocks. The language has a few English keywords to construct properties. The property language considers only two relations from the clock algebra the alternation (CCSL notation ∼) of two clocks, that ticks alternatively, and the precedence of two clocks (CCSL notation ≺ or ).

1 not to be confused with PSL [START_REF]Property Specification Language (PSL)[END_REF] The properties follow a generic pattern, which is:

Scope Clock_Expr Relation Clock Constraints

The scope can be either always if the property must hold at any time (the default value), or never if something must never happen. A Clock Expression results either directly from clock operators or from pattern matching. TPSL allows some kind of regular expression pattern matching on clock ticks in clock expressions, using two operators noted < and <>. A < B means that a precondition of a relation is that one or more ticks of A were encountered since the property is active, followed by a tick of B. The expression A <> B is matched whenever A is followed by B or B is followed by A. These pattern matching expressions must not be confused with the precedence relation ≺ of CCSL. The relation following a clock expression is either causes or alternates to mean that two clocks must constantly alternate.

The causality relation is refined in TPSL because it is often necessary to distinguish in practice between necessary and sufficient condition. For example, in a train scenario, it is sufficient that the alarm is raised to stop the train, but the train may stop for other reasons. On the contrary, in a embedded software context, whenever an interrupt is raised, there should be a reason for this interrupt, otherwise it is spurious. TPSL uses the keywords 'causes' to express the sufficient condition and the condition becomes necessary by adding the '!' symbol, thus in TPSL, Alarm causes TrainStop and Signal causes! Interrupt.

Another distinction in causality is uniqueness. It may be the case that multiple ticks of some clock should result in one single consequence (several reasons can cause the alarm but the train stops only once) or conversely that each action will result into another action (each message sent should be received). The keyword each is used to specify the latter.

In addition to the relationship, there may be additional constraints on the property and we distinguish three of them:

• additional conditions that must be satisfied. These conditions may be related to the respective timing of the clock ticks, or the count of tick occurrences. The keyword satisfies expresses such conditions. It means that for the property to be verified, not only the clock ticks must obey the rules but there are additional constraints on these clocks, typically a time delay between some ticks. • suspending condition for the causality, introduced by keyword unless, to express that a property has to be verified unless something happened, for example a Cancel order (e.g. Signal causes Interrupt unless Disabled). The semantics are that the property is void if the unless clause did occur between the first and the last relevant clock tick in the trace. The unless clause makes it possible to deal with absence in a number of cases, using a property template of the form Something causes Error unless OK such as the coffee machine that should return the coin when the system goes wrong Coin causes ReturnCoin unless CoffeeReady • activation condition for the property, introduced by keyword if. The semantics are that the property holds only if the condition is true. This is different from unless clause as the predicate is not a clock expression but a Boolean expression on counting occurrences or time stamps. For example, writing into the FIFO may cause the FIFO to be full only if the previous items were not sent already. • Causality properties specifying each for some clock C have a slightly different behavior. No if clause can be specified then as it would be contradictory; and the unless clause, if present, only applies to the most recent occurrence of C, but it does not void the causality of previous occurrences.

IV. IMPLEMENTATION AND TESTS A. STML

STML is a DSL described in Eclipse Modeling Framework [START_REF] Steinberg | EMF: eclipse modeling framework[END_REF] using the Xtext tool and an abstract syntax model also using Ecore. In the current implementation we have chosen to express properties based solely on logical clocks. Therefore STML currently offers basically two choices: a raw trace data item can simply be ignored because it is not relevant to the properties under investigation, or it can be translated into one or more logical clock ticks. Many irrelevant data can thus be simply eliminated from the trace.

To further reduce trace size and articulate related subsequent data items in the trace, STML allows to reduce several successive trace data items into one logical clock event if they meet some condition. Emitting a logical clock tick can be conditional based on some variable value. For example, the P ending interrupt clock can be emitted when some particular bit of the interrupt register is set to 1, but not when it is set to 0. But it can be also conditional to the occurrence of subsequent trace data. A question is thus how long can possibly be that sequence. The STML mapping library does consider only these simultaneous events (relatively to the trace clock) to map them into a single logical clock event.

To that end the mapping algorithm constructs a queue of the trace data items ordered by their timestamp in the simulation clock until the next period of the trace clock (the queue size increases as the trace clock frequency decreases). It then considers the earliest data and checks the STML rule for that item. Either it is simply ignored, or it directly maps to a logical clock tick, or there is a conditional path leading to some other conclusion. That conditional path can be either based on variables extracted from the simulation data (e.g. the N th bit of a hardware register), or based on the eventuality of a later event. In which case that event is searched in the queue. When such a conditional path successfully leads to a clock emit, or to an ignore, all of the matching events are removed from the queue.

Hence, any trace data in STML is considered only once. The same data item cannot participate twice or more into generating a clock tick. However it may be that simultaneous trace data (with respect to the trace clock) generate simultaneous logical clock ticks. As a result, the trace mapping is reducing a potentially huge trace stream into a much smaller one, by eliminating irrelevant data, or by synthesizing complex data possibly expressed on many bytes of data into a single logical clock tick. The logical clock trace stream can then be analyzed with automata resulting from TPSL compiler.

The Java code generator of STML generates the C++ code of the dynamic mapping library. This dynamic library is loaded by the simulator that generate trace files (in our case SystemC models). One can change the STML rules to generate different clocks to check different properties without recompiling the simulator.

B. TPSL Implementation

TPSL is the Trace Property Specification Language. It is also described as a DSL in the Eclipse EMF framework. The TPSL compiler generates code for the Verifier, which consists mainly of the dispatcher and the parallel automata. The dispatcher dispatches basic clock ticks to the parallel automata. Since the TPSL compiler knows which clock operators and which automata depends on which clock ticks, it can directly dispatch the basic clock input to their destination. A clock tick may be dispatched to multiple automata as the rules are verified in parallel.

The automata generated by TPSL compiler are symbolic finite state transducers (SFTs) [START_REF] Van Noord | Finite state transducers with predicates and identities[END_REF], [START_REF] Veanes | Symbolic finite state transducers: Algorithms and applications[END_REF], which are an extension of traditional automata by allowing transitions to be labeled with arbitrary formulas in a specified theory, in our case Boolean expressions over clock ticks.

The alphabet for the transducers is the clock alphabet, but the transducers include predicates that operate over a finite set of clock ticks. When taking transitions, the transducers add information to the history of the automaton. The history is a summary of what happened during the evaluation of the property. It maintains the timestamps and the counters of the first and last tick of each kind involved in the automaton, or the absence thereof. The TPSL language has been designed precisely such that maintaining the history is sufficient to evaluate the transducers predicates and no other data is required. The history has enough information to guarantee the decidability of the predicate bearing on the unless or if or satisfies clauses. The history is reset to void each time an automaton enters its initial state.

By construction of TPSL, the automata all follow some particular template. The TPSL compiler has several such templates that it can instantiate with the particular clock expressions and predicates from the source code. This can be implemented using C++ template facility combined with a functional style since C++ 11 standard supports lambda expressions. Figure 3 Each template automaton has a specific behavior, but they all follow the same generic pattern: from the initial state (when it is active), each transducer reaches a checking state when the clock expression has been encountered. The clock expression is actually a dependent automaton, also generated from a template. The transition is taken in the property automaton when the clock expression automaton terminates. Based on the predicates result, the transducer will either return to the initial state, or stay in the same state, or move to the expecting state. It then expects some clock tick to occur to meet the required causality. If it does, it will transition to the satisfying state, and then, depending on the required condition to be satisfied (which also bears on the history) it returns to its initial state or fails. While in the expecting state, other ticks may occur, for example a tick that would make true the unless clause, or ticks that modify the count in the history. Such clock ticks are told to be Relevant Clocks. An automaton may take a transition when some such relevant tick occurs. The TPSL compiler knows what these relevant ticks are and can generate code that dispatches only such relevant ticks to each automaton, and generate appropriate tests in the automata to take the transition.

The each automaton template is particular. As shown on Fig. 4, it is a potentially infinite automaton that waits for N expected clocks after N causing clocks. It will return to initial state only when the counts of occurrences are equal.

C. The Verifier

The verifier takes as input a trace file and outputs either a success message or one or more failures in the properties verification. Some part of the verifier is pre-defined (manually coded), essentially a library of base functions, the automaton templates, and the skeleton of the main program. This predefined code is linked with the code generated by the TPSL compiler to form a program that can read trace streams. It works essentially as follows: each line of a trace file has a time stamp (in the trace clock) followed by all clock ticks that occured on that timestamp.

The verifier first calls a setup function that creates all of the parallel automata generated by the TPSL compiler and it activates those that are in the always or never scope. It then repeatedly reads the trace stream line by line and calls the dispatcher, which does dispatch the clock ticks to each individual active automata. During this process some automaton may become active or inactive, and some may report failure. At the end of the trace stream, the automata are checked again to know if they have reached their final state. For any property, if the end of the trace is reached while a transducer is still in expecting state, it is a failure.

D. Tests

We have instrumented the SimSoC PowerPC Instruction Set Simulator and several simulation modules to generate traces according to our model. The main experiment simulates a multi-core PowerPC processor connected to various peripherals, including a OpenPIC multi-core interrupt controller and a transmit controller with a FIFO queue, controlled by a real time monitor that has drivers monitoring the peripherals. The hardware FIFO queue receives items asynchronously from the software but it re-sends these items synchronously. The hardware should (i) emit a F if oEmpty signal whenever the queue is empty (including at initialization phase), a F if oF ull signal when it is full, (ii) send the data items with limited jittering (iii) have a reasonable latency. On the other hand, the software should (i) enable interrupts and service them, and (ii) it should not put data into the queue when it is full.

The embedded software is a test program that emits data to the outside world, using the drivers and receiving interrupts. The software must enable interrupts, service them and acknowledge the interrupts with the hardware. On a multi-core configuration like the one considered, this induces a cascade of acknowledgments and signals. The full simulation trace contains data regarding hardware states, transactions between the cores and the controllers and events happening consequently to software actions.

The TPSL code regarding our example is: In this case, the trace reduction is very significant, as we are only verifying properties related to the FIFO queue and consequently many other hardware states are simply ignored. The trace file emitted by the test on SimSoC for a duration of 27 seconds is 21.162 Megabytes, the size of the logical clock trace is only 1.737 Megabytes, a reduction of over 90%. The trace is analyzed by our tool in less than 1 second.

V. CONCLUSION

We have demonstrated the feasibility of runtime verification of system properties for investigating system failure from simulation trace analysis, by abstracting raw binary simulation trace files into logical clock trace files. Thanks to a DSL that we have defined based on a simple yet powerful trace meta model, the huge binary trace files can be reduced and mapped to much smaller abstract traces relating occurrences of logical clock ticks. A language has been defined to capture temporal properties, compiled into parallel automata executed by a runtime verifier.

Our goal is not to formally offer more expressive power than any variant of temporal logic. Our tool is meant to be a convenience tool for engineers that have limited knowledge of formal methods but still can express temporal and causal properties in a easy-to-use language that is supported by a user-friendly environment, namely the Eclipse platform. Fig. 5: TPSL code for our example Another advantage of our technique is that it does not require to recompile the simulator to change the mapping and verify different types of properties. Different mappings can be achieved and different properties verified from a single trace. This provides a mechanism for the engineers to explore several paths for investigating a failure cause.

The SimSoC implementation has been modified so that a new argument can specify the trace mapping module to be loaded. In such case, the trace stream is generated such that it can be pipelined with the verifier. The verification can occur at runtime, concurrently with the simulation (possibly using at least two of the cores on a multi-core platform, introducing no penalty on the runtime verification), somehow like an extension of an assertion verifier, except that the assertions can be modified and checked without recompiling the simulator. In addition, our system can also work as a static analyzer for trace data stored into files.

In the current implementation, the mapping language translates binary data into genuine clock ticks. We are considering in the future to associate attributes to the clock ticks, then the predicates in the transducers could also evaluate Boolean expressions over such attributes, which would make it more powerful.

Fig. 1 :

 1 Fig. 1: Global architecture

Fig. 2 :

 2 Fig. 2: STML example

Fig. 3 :

 3 Fig. 3: Automaton template with if and unless clauses

Fig. 4 :

 4 Fig. 4: Automaton template for each clause