Angel Duran

Denys Dutykh
email: denys.dutykh@ucd.ie

Dimitrios Mitsotakis

On the Galilean invariance of some dispersive wave equations

Keywords: water waves, Galilean invariance, Boussinesq equations, Peregrine system, BBM equation, dispersive waves, solitary waves 2010 Mathematics Subject Classification. 76B25 (primary), 76B07, 65M70 (secondary)

Surface water waves in ideal fluids have been typically modeled by asymptotic approximations of the full Euler equations. Some of these simplified models lose relevant properties of the full water wave problem. One of these properties is the Galilean symmetry, i.e. the invariance under Galilean transformations. In this paper, a mechanism to incorporate Galilean invariance in classical water wave models is proposed. The technique is applied to the Benajmin-Bona-Mahony (BBM) equation and the Peregrine (classical Boussinesq) system, leading to the corresponding Galilean invariant versions of these models. Some properties of the new equations are presented, with special emphasis on the computation and interaction of solitary wave solutions. A comparison with the Euler equations demonstrates the relevance of the Galilean invariance in the description of water waves.

Introduction

The purpose of this paper is to investigate the relevance of the Galilean invariance in the description of water waves. A mechanism to incorporate invariance under Galilean transformations in some classical approximate models without this property is introduced. The corresponding Galilean invariant versions of two models, the BBM equation and the Peregrine system (also known to as the 'classical' Boussinesq system), are formulated. These new versions are compared numerically with their non-invariant counterparts, with some other classical models and with the Euler equations. The comparison is focused on the existence and the dynamics of solitary waves.

A natural argument in mathematical modeling is the inheritance of the physical properties of the phenomenon under consideration through the introduction of mathematical devices. In the case of the water wave theory, approximations to the Euler equations lead to some mathematical models where some fundamental properties of the original problem can be lost. This is relevant in the case of symmetries. For example, when asymptotic expansions around the still water level are performed, the invariance under vertical translations can be broken and the derived model is valid only in this particular frame of reference. Dispersive wave models possessing the property of invariance under vertical translations have been shown to be particularly robust for the simulation of the long wave runup, cf. [START_REF] Dutykh | Dispersive wave runup on non-uniform shores[END_REF][START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF]. A second symmetry, on which this paper is focused on, is related to the universality of mechanical laws in all inertial frames of reference. The Galilean invariance (or Galilean relativity) is one of the fundamental properties of any mathematical model arising in classical mechanics. This principle was empirically established by Galileo Galilei 55 years before the formulation of Newton's laws of mechanics in 1687, [52]. Nowadays, it is common to speak about this principle in terms of a symmetry of the governing equations. Hereinbelow, the Galilean invariance property is referred to the fact that the governing equations are invariant under a Galilean boost transformation, in the sense that this kind of transformations preserves the space of the solutions of the respective problem, [54]. For example, the particular form of the Galilean boost will depend on the fact that if the KdV (or BBM) equation is written in terms of the free surface elevation or the horizontal velocity variable. Nevertheless, in both cases one can easily check whether a model admits Galilean invariance as a symmetry or not.

The complete water wave problem possesses naturally this property. (For a systematic study of the symmetries and the conservation laws of the full water wave formulation we refer to the work of [START_REF] Benjamin | Hamiltonian structure, symmetries and conservation laws for water waves[END_REF], [START_REF] Benjamin | Hamiltonian structure, symmetries and conservation laws for water waves[END_REF], and in particular §4.1). Nevertheless, numerous dispersive wave equations are not invariant under the Galilean transformation. This issue was already addressed in a previous study by [START_REF] Christov | An energy-consistent dispersive shallow-water model[END_REF], [18]. We note that some fully nonlinear approximations such as the Nonlinear Shallow Water Equations (NSWE), [23], the improved Shallow Water Equations, [START_REF] Dutykh | Shallow water equations for large bathymetry variations[END_REF], and the Serre-Green-Naghdi equations (sometimes also referred to as the Su-Gardner equations), [64,[START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF]35,36,[START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF], are invariant under the vertical translation and the Galilean boost. On the other hand most of the Boussinesq models are not Galilean invariant, cf. [START_REF] Bona | A model for the two-way propagation of water waves in a channel[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]7].

Some consequences of the presence of symmetries in the models are well known. They include, for example, the Hamiltonian formulation and the generation of conserved quantities. In the case of Galilean invariance, it is worth mentioning its influence in the existence and stability of periodic traveling wave solutions (see e. g., [START_REF] Pava | Stability of cnoidal waves[END_REF]14]). It is also noted that the idea of exploiting symmetries of continuous equations has already been shown very beneficial in improving the behavior of underlying numerical discretizations, [START_REF] Kim | Invariantization of numerical schemes using moving frames[END_REF]41,16,17]. However, to our knowledge, the practical implications of the loss of the Galilean symmetry are not sufficiently known. In the present study we try to shed some light on this issue and its influence on the approximate dispersive wave models. Some other properties, such as mass, momentum and energy conservation, are at least equally important. These questions have been already addressed in the literature and, consequently, are not central to our study.

The paper is organized as follows. In Section 2 some classical models of water wave theory are reviewed and the invariant counterparts of the BBM equation and the classical Peregrine system are derived. The methodology of this derivation can be extended to other non invariant systems. These new equations are not attempted to provide new asymptotic models, but to investigate the implications of the introduction of the Galilean symmetry by comparison with the original models, other classical equations and the full Euler system. To this end, in Section 3, solitary wave profiles of the new equations are numerically generated. The comparison with the wave profiles of other approximations (some with exact formulas) and of the Euler equations (by using Tanaka's algorithm and Fenton's asymptotic solution) is established in terms of the amplitude-speed and amplitudeshape relations. The interactions of solitary waves for the new models, with head-on and overtaking collisions, are studied in Section 4. Finally, the main conclusions of this study are outlined in Section 5.

Mathematical models

Consider an ideal fluid of constant density along with a cartesian coordinate system in two space dimensions (x, y). The y-axis is taken vertically upwards and the x-axis is horizontal and coincides traditionally with the still water level. The fluid is bounded below by an impermeable horizontal bottom at y = -d and above by an impermeable free surface at y = η(x, t). We assume that the total depth h(x, t) ≡ d + η(x, t) remains positive h(x, t) h 0 > 0 at all times t.

Assuming that the flow is incompressible and irrotational, the governing equations of the classical water wave problem, [START_REF] Stoker | Water waves, the mathematical theory with applications[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF], are the following:

φ 2 xx + φ 2 yy = 0 -d y η(x, t), (2.1)
η t + φ x η x -φ y = 0 y = η(x, t), (2.2)
φ t + 1 2 (φ x ) 2 + 1 2 (φ y ) 2 + g η = 0 y = η(x, t), (2.3)
φ y = 0 y = -d, (2.4)
with φ being the velocity potential (by definition, the irrotational velocity field u = φ x and g the acceleration due to the gravity force. The water wave problem possesses Hamiltonian, [START_REF] Petrov | Variational statement of the problem of liquid motion in a container of finite dimensions[END_REF][START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF][START_REF] Broer | On the Hamiltonian theory of surface waves[END_REF][START_REF] Salmon | Hamiltonian fluid mechanics[END_REF], and Lagrangian, [START_REF] Luke | A variational principle for a fluid with a free surface[END_REF][START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF], variational structures. The symmetry group of the complete water wave problem (2.1) -(2.4) was described by [START_REF] Benjamin | Hamiltonian structure, symmetries and conservation laws for water waves[END_REF] in [START_REF] Benjamin | Hamiltonian structure, symmetries and conservation laws for water waves[END_REF]. In particular, the full formulation of the water wave equations admits the Galilean boost symmetry and the invariance under the vertical translations (the latter issue will be addressed by the authors in a future work). However, the water wave theory has been developed from the beginning by constructing various approximate models which may conserve or break some of the symmetries, [START_REF] Craik | The origins of water wave theory[END_REF]. Below we consider several classical models and discuss their Galilean invariance property.

The KdV equation

The unidirectional propagation of long waves in the so-called Boussinesq regime (where the nonlinearity and the dispersion are of the same order of magnitude), [START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF]7], can be described by the celebrated Korteweg-de Vries (KdV) equation [42,[START_REF] Johnson | A Modern Introduction to the Mathematical Theory of Water Waves[END_REF], which in dimensional variables can be written in the form:

u t + gd u x + 3 2 uu x + d 2 6 gd u xxx = 0, (2.5)
where u(x, t) is the horizontal velocity variable which is usually defined as the depthaveraged velocity of the fluid, [60], or the fluid velocity measured at some specific water depth, [START_REF] Bona | A model for the two-way propagation of water waves in a channel[END_REF]. Some well-known properties of (2.5) are reviewed (see e.g. [34,45,51]). First, the KdV equation is an integrable model with the following two-parameter family of solitary wave solutions:

u(x, t) = u 0 sech 2 1 2 κ(x -c s t -x 0 ) , c s = gd + u 0 2 , (κd) 2 = 3u 0 √ gd , u 0 > 0, x 0 ∈ R.
The initial value problem of the Korteweg-de Vries (KdV) equation possesses a Hamiltonian structure

u t = J δH δu ,
(where δ denotes the variational derivative) in a suitable phase space of functions vanishing, along with some of their derivatives, at infinity. The skew-symmetric operator J and the Hamiltonian functional H are

J = -∂ x , H = 1 2 ˆR gdu 2 + 1 2 u 3 -gd d 2 6 u 2 x dx.
The Hamiltonian H is the third conserved quantity of the well-known hierarchy of invariants for (2.5), [45].

The central question in our study is the Galilean invariance of model equations. We recall that the KdV equation (2.5) possesses this property. For a systematic derivation of symmetries of (2.5) we refer to [51] and [54]. To verify the invariant under Galilean transformations we choose another frame of reference which moves uniformly, for example, rightwards with constant celerity c. This symmetry is expressed by the following transformation of variables:

x

→ x -3 2 ct, t → t, u → u + c. (2.6)
In this moving frame of reference (2.5) becomes:

u t - 3 2 cu x + gd u x + 3 2 (u + c)u x + d 2 6 gd u xxx = 0.
After some simplifications one can recover the original KdV equation, which completes the proof of the invariance.

In order to assess the relative magnitude of various terms in equation (2.5) scaled variables are introduced. The classical long wave scaling is the following:

x ′ := x ℓ , y ′ := y d , t ′ := g d t, η ′ := η a , u ′ := u √ gd , (2.7)
where h 0 , a, ℓ are the characteristic water depth, wave amplitude and wave length respectively. Using these three characteristic lengths we can form three following important dimensionless numbers:

ε := a d , µ 2 := d ℓ 2 , S := ε µ 2 .
(2.8)

Parameters ε ≪ 1 and µ 2 ≪ 1 characterize the wave nonlinearity and dispersion, while the so-called Stokes number S measures the analogy between these two effects. In the Boussinesq regime the Stokes number is of order of one, S = O(1), which establishes that the dispersion and the nonlinear effects are comparable. We note that in [START_REF] Bona | Solitary waves and other phenomena associated with model equations for long waves[END_REF] the range of validity of the KdV had been studied and appeared that for the values of the Stokes number 0.5 ≤ S ≤ 10 appears to have excellent performance, while even for larger values the results are acceptable too. Using these dimensionless and scaled variables the KdV equation (2.5) can be written in the form:

u t + u x + 3 2 εuu x + µ 2 6 u xxx = 0. (2.9)
where the primes have been dropped. Formulas (2.7) and (2.8) will also be used in some of the developments below.

The BBM equation

Benjamin, Bona & Mahony (1970), [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF], (see also [START_REF] Peregrine | Calculations of the development of an undual bore[END_REF]) proposed the following modification of the KdV equation, known as the BBM equation:

u t + gd u x + 3 2 uu x - d 2 6 u xxt = 0. (2.10)
The main idea for the derivation of this model is to use the lower order relation between time and space derivatives from (2.9) u t = -u x + O(ε, µ 2 ) to modify the higher-order dispersive term u xxx = -u xxt + O(ε, µ 2 ), cf. [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF]. One of the main practical motivations for this modification is to improve the dispersion relation properties of the KdV equation. Specifically, unlike the KdV equation, the phase and group velocities of the BBM equation have a lower bound. It is also noted that the BBM equation has the following solitary wave solutions:

u(x, t) = u 0 sech 2 1 2 κ(x -c s t -x 0 ) , c s = gd + u 0 2 , (κd) 2 = 3u 0 √ gd + 1 2 u 0 , x 0 ∈ R.
The BBM equation is not integrable, but it can also be written as an infinite-dimensional Hamiltonian system

u t = J δH δu ,
where the operator J and the Hamiltonian functional H are defined as:

J = (1 - d 2 6 ∂ xx ) -1 • (-∂ x ), H = 1 2 ˆR gdu 2 + 1 2 u 3 dx. (2.11)
and the structure is defined on a phase space similar to that of the KdV equation.

As far as the Galilean invariance is concerned, the change of variables (2.6) applied to (2.10) leads to

u t - 3 2 cu x + gd u x + 3 2 (u + c)u x - d 2 6 u xxt + d 2 4 cu xxx = 0,
and after some algebraic simplifications we obtain:

u t + gd u x + 3 2 uu x - d 2 6 u xxt + d 2 4 cu xxx = 0.
Since there is at least one new term ( d 2 4 cu xxx ) appeared in the previous moving frame of reference, the BBM equation is not Galilean invariant. (The relevance of this drawback always puzzled the researchers, cf. [18].).

The iBBM equation

We now present a modification of the classical BBM equation which allows us to recover the Galilean invariance property. Furthermore, the idea behind the arguments below can be applied to other models. The strategy is to add a new ad-hoc term which will vanish the non-invariant contribution of the BBM dispersion u xxt under the transformation (2.6). The resulting equation, which will be called invariant Benjamin-Bona-Mahony (iBBM) equation, takes the form:

u t + gd u x + 3 2 uu x - d 2 6 u xxt - d 2 4 uu xxx = 0. (2.12)
It is straightforward to see that (2.12) is invariant under the Galilean transformation (2.6).

The modification proposed above becomes more transparent in scaled variables. The application of the long wave limit (2.7) to (2.12) leads to

u t + u x + 3 2 εuu x - µ 2 6 u xxt - εµ 2 4 uu xxx = 0.
One can observe that the last term on the left hand side, responsible for the Galilean invariance of (2.12), is a nonlinear term of order O(εµ 2 ) and consequently, it is asymptotically negligible in the BBM formulation. Since this additional term is nonlinear, the linear dispersion relation of (2.10) is not modified. Its effect will be studied thoroughly in the following sections. We stress out that the original BBM and iBBM equations are coincide up to order O(ε, µ 2 ), cf. [START_REF] Dullin | On asymptotically equivalent shallow water wave equations[END_REF].

Remark 1. Unlike the BBM equation (2.10), this invariant version (2.12) does not possess, to our knowledge, a Hamiltonian structure. However, it is possible to propose an invariantization which preserves this variational formulation along with the Galilean invariance.

The alternative given by the equation

u t + gd u x + 3 2 uu x - d 2 6 u xxt - d 2 4 2u x u xx + uu xxx = 0. (2.13)
has an additional higher-order nonlinear term which allows for a non-canonical Hamiltonian structure. In this case, the operator

J = (1 -d 2 6 ∂ xx ) -1 • (-∂ x )
is the same as for the BBM equation and the Hamiltonian H is

H = 1 2 ˆR gdu 2 + 1 2 u 3 + d 2 4 uu 2 x dx.
We underline that the new terms in (2.13) can be also found in several models such as the Camassa-Holm [13], Burgers-Poisson [START_REF] Fetecau | Aproximate model equations for water waves[END_REF] and Degasperis-Procesi [24] equations. Both equations (2.12), (2.13) are Galilean invariant versions of the BBM equation. For the comparisons performed in Section 3, (2.12) has been considered, with the purpose of focusing exclusively on the effects of the Galilean invariance property, which is the main goal of the paper. A similar study could be done with the alternative (2.13).

We now look for travelling wave solutions of (2.12) of the form:

u(x, t) = u(ξ), ξ := x -c s t, (2.14)
where c s is the solitary wave speed. We also assume that u(ξ) decays to zero along with all derivatives when |ξ| → ∞. Substituting (2.14) into the iBBM equation (2.12) leads to the ordinary differential equation:

( gd -c s )u ′ + 3 4 (u 2 ) ′ + c s d 2 6 u ′′′ - d 2 4 uu ′′′ = 0, (2.15)
where the prime denotes differentiation with respect to ξ. Using the boundary conditions at infinity, the identity uu ′′′ = (uu ′′ -1 2 (u ′ ) 2 ) ′ , and after an integration, equation (2.15) becomes: that can be written as a system

( gd -c s )u + 3 4 u 2 + c s d 2 6 u ′′ - d 2 6 1 2 (u ′ ) 2 -uu ′′ = 0, (2.16)
u ′ = v,
(2.17)

v ′ = 2 d 2 cs 3 -u 2 (c s -gd)u - 3 4 u 2 - d 2 8 v 2 .
(2.18)

Now it can be verified that when c s > √ gd, the origin u = v = 0 is a saddle point, as depicted in Figure 1(a), which shows the corresponding phase plane. The homoclinic trajectory O → A → B → O represents a solitary wave. (The MATLAB code for this figure can be found in [55].). In Section 3 we compute solitary wave solutions of (2.16) by numerical means.

The Peregrine system

Under the assumptions described above, D.H. Peregrine in 1967 derived the following system of equations governing the two-way propagation of long waves of small amplitude in the Boussinesq regime [60]:

η t + (d + η)u x = 0, (2.19)
u t + uu x + gη x - d 2 3 u xxt = 0, (2.20)
where u(x, t) is now defined as the depth averaged fluid velocity and η(x, t) is the deviation of the free surface of the water from its rest position. This system is also known as the classical Boussinesq system, cf.

[7] and will be denoted by cPer. In [START_REF] Bona | Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF][START_REF] Chen | Solitary-wave and multi pulsed traveling-wave solution of Boussinesq systems[END_REF], the existence and some properties of solitary wave solutions of (2.19)-(2.20) are obtained, without explicit formulas. On the other hand, to our knowledge, a Hamiltonian structure has not been found, [7].

We now study the Galilean invariance of (2. [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF])- (2.20). In this case, the Galilean transformation takes the following form:

x → x -ct, t → t, η → η, u → u + c.
(2.21)

The mass conservation equation (2.19) in new variables reads:

η t -cη x + (d + η)(u + c) x = 0.
After simplifications one can see that this equation remains invariant under the transformation (2.21). Now let us consider the momentum balance equation (2.20). In the moving frame of reference this equation becomes:

u t -cu x + (u + c)u x + gη x - d 2 3 u xxt + c d 2 3 u xxx = 0,
and after some simplifications

u t + uu x + gη x - d 2 3 u xxt + c d 2 3 u xxx = 0.
As in the BBM case, a new dispersive term c d 2 3 u xxx appears, showing that the system (2. [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF])-(2.20) is not Galilean invariant.

The iPeregrine system

Following the same technique as in the case of the BBM equation, it is possible to derive a modification of the classical Peregrine system (2. [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF])-(2.20) which will allow us to recover the Galilean invariance property. It can be done in a similar way leading to the iBBM equation (2.12). The corresponding system reads:

η t + (d + η)u x = 0, (2.22)
u t + uu x + gη x - d 2 3 u xxt - d 2 3 uu xxx = 0. (2.23)
Note that since the mass conservation equation is invariant, it is not modified in the new version. Now it is straightforward to check the invariance of equation (2.23). Therefore, system (2.22)-(2.23), which will be called invariant Peregrine system or iPer for the sake of conciseness, is Galilean invariant. In dimensionless and scaled variables, the system reads:

η t + (1 + εη)u x = 0, u t + εuu x + η x - µ 2 3 u xxt - εµ 2 3 uu xxx = 0,
and one can see that the new ad-hoc term is of higher-order and, asymptotically speaking, negligible. As in the case of the

η(x, t) = η(ξ), u(x, t) = u(ξ), ξ := x -c s t,
where η and u decay to zero, along with their derivatives, as |ξ| → ∞. After substituting this representation into the governing equations (2.22)-(2.23) they become:

-c s η ′ + (d + η)u ′ = 0, -c s u ′ + 1 2 (u 2 ) ′ + gη ′ + c s d 2 3 u ′′′ - d 2 3 uu ′′′ = 0.
An integration of the mass conservation equation and the decay at infinity lead to

u = c s η d + η . (2.24)
Then the momentum balance equation can be integrated once and we have

-c s u - d 2 3 u ′′ + 1 2 u 2 + gd • u c s -u - d 2 3 1 2 (u ′ ) 2 -uu ′′ = 0. (2.25)
Similarly to the case of the iBBM, one can see that, (2.25) written as a first order system and when c 2 s -gd > 0, the origin is a saddle point; the phase plane sketched in Figure 1(b) also shows a solitary wave, in the form of a trajectory O → A → B → O.

The Serre equations

In order to complete the presentation of our model equations, we consider the fullynonlinear system of equations referred to as the Serre equations, [64, 3, 29]:

h t + [hu] x = 0, (2.26)
u t + uu x + gh x = 1 3 h -1 h 3 (u xt + uu xx -u 2 x ) x , (2.27)
where h(x, t) := d + η(x, t) is the total water depth variable. Solitary wave solutions of (2.26)-(2.27) are explicitly known. They are given by the formulas:

η(x, t) = a 0 sech 2 1 2 κ(x -c s t -x 0 ) , u = c s η d + η , (2.28)
where c s = g(d + a 0 ) and (κd) 2 = 3a 0 d+a 0 . As pointed out in [START_REF] Johnson | Camassa-Holm, Korteweg-de Vries and related models for water waves[END_REF], [START_REF] Li | Linear stability of solitary waves of the Green-Naghdi equations[END_REF] and [START_REF] Li | Hamiltonian structure and linear stability of solitary waves of the Green-Naghdi equations[END_REF], the Serre equations possess a Hamiltonian structure:

qt h t = - ∂ x [q + q∂ x [ h∂ x [ ∂ x [h 0 • δ H / δ q] δ H / δh] ,
where the Hamiltonian functional H is given by

H = 1 2 ˆR h u 2 + 1 3 h 3 u 2 x + g η 2 dx,
The variable q is sometimes referred to as the potential vorticity flux and is defined by

q := h u -1 3 [ h 3 u x ]
x . The Serre equations (2.26)-(2.27) can be shown to have the Galilean invariance property. For the mass conservation equation (2.26) we refer to Section 2.5. Thus it remains to check this property for the momentum conservation equation (2.27). If we make the change of variables t → t, x → x -ct, h → h and u → u + c as before, equation (2.27) becomes:

u t -cu x + (u + c)u x + gh x = 1 3 h -1 h 3 (u xt -cu xx + (u + c)u xx -u 2 x ) x ,
and after simple algebraic simplifications one can recover the original equation (2.27).

Numerical computation of travelling waves

In the previous section we presented several models arising in water wave theory. Moreover, we proposed two novel equations, the iBBM equation and the iPer system, with the aim of incorporating the property of invariance under the Galilean transformation, lost by the original BBM equation and the cPer model. The purpose of this and the next sections is to compare these models through the computation of their respective solitary wave solutions and, whenever possible, with the solitary waves of full Euler equations. In some cases, approximate solitary waves must be generated numerically. In Section 3.1 the great lines of the numerical procedure to this end is given, while the methods used to construct approximate solitary waves of the Euler equations are mentioned in Section 3.1.1, where a comparative study is carried out. In the sequel we consider all the models in nondimensional but unscaled variables (i.e. when ε = µ 2 = 1).

Computation of travelling-wave profiles. The Petviashvili method

The investigation for travelling wave solutions in one-dimensional systems typically leads to a set of differential equations of the form

LU = N(U), (3.1)
for some differential operators L (linear) and N (nonlinear). The numerical resolution of the preceding system can be done in many different ways (see [START_REF] Yang | Nonlinear Waves in Integrable and Nonintegrable Systems[END_REF] and the references therein as a modest representation of the related literature). Among all the possibilities, the socalled Petviashvili method will be used in our computations. This method stems from the pioneering work of V.I. [START_REF] Petviashvili | Equation of an extraordinary soliton[END_REF], [62]. It is based on a modification of the classical fixed point iteration (which in these cases is usually divergent) and it is formulated as follows. Given an initial profile U 0 , the Petviashvili iteration generates approximations U n of the original solution of (3.1) following the formulas

M n = LU n , U n N(U n ), U n , (3.2)
LU n+1 = M γ n N(U n ), (3.3)
where •, • denotes the usual L 2 inner product and γ is a free parameter that controls the convergence of the method. The quantity M n in (3.2) is called the stabilizing factor. See [58,43] for details, generalizations and local convergence results for some model equations.

In this study, the iteration (3.2)-(3.3) is applied to compute solitary wave profiles in the following cases with the corresponding operators, namely: where c s is the solitary wave speed. In order to reconstruct the free surface elevation profile from the horizontal velocity distribution u(ξ) we use the following exact formula which can be derived for the travelling wave solutions to the full Euler equations provided that u(ξ) is defined as the depth-averaged horizontal velocity:

• (iBBM): Lu = ( √ gd -c s )u + c s d 2 6 u ′′ , N(u) = d 2 4 (u ′ ) 2 2 -uu ′′ -3 4 u 2 . • (cPer): Lu = c s u -d 2 3 u ′′ , N(u) = u 2 2 + gdu cs-u . • (iPer): Lu = c s u -d 2 3 u ′′ , N(u) = u 2 2 + gdu cs-u -d 2 3 (u ′ ) 2 2 -uu ′′ ,
u(ξ) = c s η(ξ) d + η(ξ) .

Generation and numerical evolution of the profiles

In many cases, the method (3.2)-(3.3) can be efficiently implemented by using Fourier techniques, [58,43]. Specifically, our implementation for the three systems has been performed by considering the corresponding periodic problem and using a pseudospectral representation for the approximations to the profiles. As an initial approximation, a solitary wave solution (2.28) of the Serre equations or the third-order asymptotic solution of Grimshaw,[37], can be considered. The iterative procedure is continued until the difference between two consecutive iterations in the L ∞ -norm, or the L ∞ -norm of the residual is less than a prescribed small tolerance which, in our case, is of order O(10 -15 ). The convergence is reached within 10-20 iterations. In order to illustrate better the transformations that solitary waves undergo while we gradually increase the propagation speed parameter c s , we superpose several profiles on the same In order to assess the accuracy of the computations, the three models considered in this paper have been numerically integrated in time using the computed solitary wave profiles as initial conditions. Some error indicators measuring the accuracy of the numerical approximation of the solitary waves have been computed. The numerical method for the corresponding initial-periodic boundary value problem consists of a pseudospectral method for the semi-discretization in space and the classical, explicit fourth-order Runge-Kutta scheme for the time integration, [START_REF] Dougalis | Theory and numerical analysis of Boussinesq systems: A review[END_REF].

In the experiments below, we study the propagation of a solitary-wave profile, generated by (3.2)-(3.3) with speed c s = 1.1 in the interval [-128, 128], with N = 2048 nodes for the pseudospectral approximation and spatial and time step sizes ∆x = 1.25 × 10 -2 , and ∆t = 1.25×10 -3 respectively. The results correspond to cPer and iPer systems. The iBBM equation has also been implemented, with similar results.

We study the evolution of two parameters: the normalized amplitude error and the shape error. The first one is computed by comparing, at each timestep, the initial amplitude of the profile (generated by the Petviashvili method) with the corresponding amplitude of the numerical solution computed using Newton's method, [START_REF] Dougalis | A numerical study of the stability of solitary waves of Bona-Smith family of Boussinesq systems[END_REF]. The L 2 based, normalized shape error is also defined at each timestep for, lets say, the solitary wave u as SE n = inf τ U n -u(•, τ ) / u(•, 0) . For the computation of the shape error we compute the time τ * neat the timestep t n such that d dτ ξ 2 (τ * ) = 0, where ξ(τ ) := U n -u(•, τ ) / u(•, 0) using Newton's method and an initial approximation τ 0 = t n -∆t. Then the shape error is the quantity SE n = ξ(τ * ). Their computation is implemented as in e.g. [START_REF] Dougalis | A numerical study of the stability of solitary waves of Bona-Smith family of Boussinesq systems[END_REF]. Figure 3(a) shows the temporal evolution of this amplitude error up to a final time T = 100, for both cPer and iPer. We observe that for the specific values of ∆x and ∆t the amplitude is conserved up to 10 decimal digits in both cases.

On the other hand, the shape error is computed by comparing the numerical solution with time translations of the initial profile with the prescribed speed and minimizing the differences (see [START_REF] Dougalis | Theory and numerical analysis of Boussinesq systems: A review[END_REF] for the details). The results displayed in Figure 3(b) show a virtually constant evolution of this error, which in both cases is of order O(10 -10 ). These results confirm the accuracy of the technique used to generate the solitary wave profiles and of the numerical code for the time evolution. The latter will be used for the experiments below.

Numerical results

In this section we compare the solitary waves of the following models:

• The Korteweg-de Vries (KdV) equation (2.5).

• The Benjamin-Bona-Mahony (BBM) equation (2.10).

• The invariant Benjamin-Bona-Mahony (iBBM) equation (2.12).

• The Peregrine classical Peregrine (cPer) system (2. [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]), (2.20).

• The invariant Peregrine invariant Peregrine (iPer) system (2.22), (2.23).

• The Serre equations (2.26), (2.27), from either the analytical formula (when possible) or the computations with the Petviashvili method. The comparison is established between them and with those of the Euler equations. The construction of the approximations to travelling wave solutions for the 2D Euler equations with free surface (to be considered as reference solutions for the approximate models) will be computed with two techniques. One by using the Tanaka's algorithm, [67], and the other by using asymptotic solutions given in, [START_REF] Fenton | A ninth-order solution for the solitary wave[END_REF][START_REF] Longuet-Higgins | On the Mass, Momentum, Energy and Circulation of a Solitary Wave[END_REF]. The study is focused on the amplitudes and shapes of the computed free surface elevation η(ξ), provided by the models for the same prescribed value of the propagation speed parameter c s .

Solitary wave speed-amplitude relation

Figure 4(a) shows an amplitude-wave speed diagram for the models considered in this paper. First, an approximate relation between the amplitude and speed for the solitary waves of the full Euler system is computed, by using the Tanaka's algorithm, [67], and Fenton's 9th order, [START_REF] Fenton | A ninth-order solution for the solitary wave[END_REF][START_REF] Longuet-Higgins | On the Mass, Momentum, Energy and Circulation of a Solitary Wave[END_REF]. They virtually give the same results and these are compared with the relation obtained by each of the models. We observe that the Serre equations (and, consequently, the other systems, that have weaker nonlinearities) are known to provide a relatively good approximation to the solitary wave solutions of the full Euler equations in a range of amplitudes not greater than 0.5, cf. e.g. [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF]14]. Therefore, our attention is focused on solitary waves with these amplitudes, as it is observed in Figure 4(a). Figure 4(b) shows a magnification for the largest amplitudes.

Corresponding to the comparison between the solitary waves of the models at hand, we note that the non Galilean invariant models, the Benjamin-Bona-Mahony (BBM) equation and the classical Peregrine (cPer) system, tend to underestimate the solution speed for a given amplitude. On the other hand, the curves corresponding to their invariant counterparts along with the fully-nonlinear Serre equations lie above the reference solution. In particular, the results for the invariant Benjamin-Bona-Mahony (iBBM) equation are very close to those of Serre equations. A surprising fact is that the amplitude-speed relation given by the invariant Peregrine (iPer) system is superposed with that of the Serre equations up to the graphical resolution.

Solitary wave shape

A second comparison between the shape of the computed solitary waves, is presented in Figures 5(a)-5(f). They illustrate the cases of small (a/d ≈ 0.1, see Figure 5(a) and a magnification on the wave crest in Figure 5(b)), moderate (a/d ≈ 0.22, see Figures 5(c wave amplitudes of the models (within the range mentioned above). According to these results, it is observed that the invariant Benjamin-Bona-Mahony (iBBM) equation and the invariant Peregrine (iPer) system approximate much better the amplitude of the reference solution (represented, in this case, by the Tanaka's solution) than the non-Galilean invariant counterparts, and they stay very close to the results of the Serre system near the crest.

As a measure of the level of approximation to solitary wave solutions of the Euler system, these and the previous results show the benefits of taking into account the invariantization process in the approximate models.

Solitary waves dynamics

We complete the numerical experiments by studying the effects of the Galilean invariance property in the evolution of solitary waves. Specifically, we first compare, by numerical means, head-on collisions of two solitary waves of the invariant BBM and Peregrine equations with those of their corresponding not invariant models. Then, we will study the numerical solution of the four systems when a solitary wave profile of the Euler equations is used as initial condition.

It is known that the tails produced by the interaction of two solitary waves are sensitive to both the linear terms (characterizing the linear dispersion relation) and the nonlinearities. For example, two solitary waves of the KdV equation interact in an elastic way without producing dispersive tails at all, [START_REF] Zabusky | Interaction of solitons in a collisionless plasma and the recurrence of initial states[END_REF], while the collision of solitary waves of the BBM equation will produce dispersive tails and probably small-amplitude nonlinear pulses as the main indication of an inelastic interaction, [START_REF] Bona | Solitary-wave interaction[END_REF][START_REF] Bona | An Evaluation of a Model Equation for Water Waves[END_REF]. In the new Galilean invariant models, the new nonlinear terms are of order εµ 2 and their effects on the interaction will be studied here. The same code introduced in Section 3 is used for the numerical computations, as well as the Petviashvili method (3.2)-(3.3) to generate solitary wave profiles when necessary.

Head-on collisions of solitary waves

A first group of experiments concerns head-on collisions. The classical Peregrine (cPer) and the invariant Peregrine (iPer) systems have been considered, by constructing, in both cases, two solitary waves on the interval [-256, 256] with speeds c s,1 = 1.15 and c s,2 = 1.05 (translated appropriately such that their maximum values are achieved on x = -50 and x = 50 respectively) and travelling in opposite directions. These solitary waves are of small amplitude and their shapes are almost the same for both models. The code uses N = 4096 nodes, a spatial meshlength of ∆x = 1.25 × 10 -1 and ∆t = 5 × 10 -3 as the time step.

Figure 7 shows the η-profile of the head-on collision for both models and at several times. A tail behind each solitary pulse after the collision is observed. This is larger in the case of the invariant Peregrine system (of the order of 10 -4 ) than in the case of the classical Peregrine system (approx. 10 -4 ), see Figure 7(d) and (f). During the collisions, a similar phase shift also takes place, see Figure 8.

A final comparison is established in terms of the degree of inelasticity of the interaction. This can be measured by using several parameters, [START_REF] Abdulloev | One more example of inelastic soliton interaction[END_REF][START_REF] Lewis | Resonant production of solitons in the RLW equation[END_REF]. In each case, a symmetric head-on collision has been implemented; that is, two solitary waves with the same speed c s travelling in opposite directions. After the interaction, both solitary waves emerge with similar amplitudes A af ter , but below the initial one A init , as can be observed in Table 1.

Then a ratio K of the amplitude of the waves after the collision to their amplitude before the collision has been computed. Figure 6 shows the behavior of this value, as a function of the speed parameter c s and for both systems. The results reveal a higher inelastic collision in the case of the invariant system, in accordance to what is observed in the full water wave model [START_REF] Craig | Solitary water wave interactions[END_REF].

Overtaking interactions of solitary waves

In order to study the overtaking collision of solitary waves, we first consider the classical and the invariant Peregrine systems and compute two solitary waves on the interval [-1024, 1024] with speeds c s,1 = 1.15 and c s,2 = 1.05 (translated appropriately such as their maximum values are achieved at x = -50 and x = 50 respectively) travelling in the same direction, with N = 16384 and timestep of ∆t = 5 × 10 -3 . In Figure 9 we observe that the basic characteristics of the interaction are similar for both systems, (cf. also [START_REF] Antonopoulos | Numerical solution of Boussinesq systems of the Bona-Smith family[END_REF]). Specifically, the interaction is again inelastic; after the collision a tail, of apparent dispersive nature, behind the smallest wave (moving to the right) and a small N-shape wavelet (moving to the left) are generated. Moreover, a small phase shift and a change in shape can be observed in the solitary pulses (Figure 8). The wavelet generated in the case of the classical Peregrine system has an inverse N-shape while the invariant version has an N-shape. With the same input data, an overtaking collision for the case of the BBM and iBBM equations is shown in Figure 10 at several times. In this case, on the contrary, no N-shape wavelet is observed and only a tail behind the waves appears. The tail following the solitary waves of the iBBM appears to be larger again as an effect of the high-order nonlinear terms. Also, the phase shift observed in the case of the iBBM is larger, indicating again a larger amount of nonlinearity and inelasticity in the specific interaction. We can conclude that these two groups of results show that, as expected, the O(εµ 2 ) nonlinear terms of the invariant models do not change significantly, in a qualitative sense, the behaviour of the interactions of solitary waves provided by the corresponding not Galilean invariant system.

Comparison with laboratory data

After observing the inelastic collisions in the previous section it rises the question of which model is more realistic. To answer this question we compare the numerical solution of the cPer and the iPer systems during a head-on collision with some experimental data from [START_REF] Craig | Solitary water wave interactions[END_REF]. Specifically, we consider two solitary waves of speeds c s,1 = 0.7721 m/s and c s,2 = 0.7796 m/s, respectively and translated such as fit with the laboratory data at t = 18.3. These solitary waves are almost identical for both models and the head-on collision appeared to be be very similar compared to the laboratory data. A closer look of the dispersive tails shows that, again, after their interaction a larder dispersive tail is generated in the case of the iPer system. Comparing with the experimental data (see Figure 11) we observe that the collision described by the iPer system is very close to the one described by the cPer and to the laboratory data as well. Moreover, the iPer system demonstrates a better performance in describing the dispersive tail showing that the iPer system could be realistic too.

Comparison with Euler equations

Finally, we study the evolution of a solitary wave of the Euler equations when we use it as initial condition to the approximate models. Specifically, we consider an approximate solitary wave Φ h (x) of amplitude A = 0.2 (and speed c s ≈ 1.095490471188718) obtained by using Fenton's ninth order asymptotic solution. In the case of the Peregrine and the iPeregrine systems we use for initial velocity u 0 (x) = csη 0 (x) /d+η 0 (x). We remind that due to the mass conservation property this formula is exact (see equation (2.24) for the velocity given the surface elevation η for both Peregrine and iPeregrine systems). In Figure 12 we present the solution at T = 100. We observe that the initial condition is resolved into a new solitary wave followed by a dispersive tail. In the case of the classical models, the dispersive tails appear to be smaller. Figure 13 (a) shows the shape error of the solution (i.e. how much different is the solution from being the exact solitary wave of the Euler equations) while Figure 13 (b) presents the amplitude of the solution as a function of time t. From these two figures we observe that in the case of classical models the emerging solitary waves are closer in shape and amplitude to the original solitary wave solution of the Euler equations than the respective solitary waves of the invariant models.

Conclusions

In the present work the influence of Galilean invariance in several equations arising in water wave modelling is studied. We modify not Galilean invariant models in order to include this fundamental property. The technique introduced here consists of adding higher-order terms from the approximation of the governing equations. These terms are asymptotically negligible and consequently, the modified models are still valid in the appropriate regime. As a case study, corresponding modifications of two not Galilean invariant models, the Benjamin-Bona-Mahony (BBM) equation and the classical Peregrine system, are presented. The comparison with reference solutions to the full Euler equations shows that this extra-term is beneficial for the description of the travelling wave solutions in several ways. First, this modification improves the solitary waves amplitude-speed relation which lies closer to the Tanaka's and Fenton's solutions. In this regard, we obtain a surprising performance of the iPeregrine system (2.22), (2.23) with the amplitude-speed relation undistinguishable from the fully-nonlinear Serre equations (2.26), (2.27). Moreover, the amplitudes of solitary wave solutions to the invariant models are closer to the corresponding full Euler solutions than the classical counterparts. The comparison is finished off with a numerical study of head-on and overtaking collisions. Compared to the behavior observed in the not Galilean invariant equations, the higher order nonlinear terms incorporated in the new models do not affect qualitatively the inelastic character of the interactions. However, a relevant difference in the degree of inelasticity is observed, being higher in the case of the Galilean invariant versions. This behavior is closer to what has been observed in the case of the full Euler equations and in laboratory experimentals.

Figure 1 .

 1 Figure 1. Phase plane of the invariant models with c s > √ gd: (a) iBBM. (b) iPer. In both cases, the trajectory O → A → B → O represents a solitary wave.

  invariant Peregrine (iPer) system

Figure 2 .

 2 Figure 2. Approximate solitary wave profiles for different speeds c s , generated by the Peviashvili method (3.2)-(3.3).

  Figure (Figure 2(a) corresponds to the invariant Benjamin-Bona-Mahony (iBBM) system and Figure 2(b) to the invariant Peregrine (iPer) system).

Figure 3 .

 3 Figure 3. Some error indicators for the propagation of a solitary for the Peregrine and the invariant Peregrine systems.

Figure 4 .

 4 Figure 4. The curves corresponding to the iPeregrine system and the Serre system are superposed up to the graphical resolution. The Tanaka solution is represented with red circles, while the Fenton solution is depicted with the red solid line.

Figure 5 .

 5 Figure 5. Solitary waves of different speeds c s .

Figure 6 .

 6 Figure 6. Symmetric head-on collision. Ratio K of amplitudes before and after the collision, for Peregrine and invariant Peregrine systems.

Figure 7 .

 7 Figure 7. Head on collision of two solitary waves for the Peregrine and invariant Peregrine systems.

Figure 8 .

 8 Figure 8. Phase diagram of the overtaking collision: (a) Peregrine system. (b) Invariant Peregrine system.

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Overtaking collision of two solitary waves traveling to the right for classical Peregrine (cPer) and invariant Peregrine (iPer) systems.

  (a) Solution of the BBM equation at T = 100 (b) Solution of the iBBM equation at T = 100 Solution of the iPeregrine system at T = 100

Figure 12 .

 12 Figure 12. Evolution of a solitary wave of the Euler equations, with BBM, iBBM, cPer and iPer models. (A = 0.2)

Figure 13 .

 13 Figure 13. Shape error and amplitude of the solution as a function of time

  )

							Amplitude -Wave speed diagram						Amplitude -Wave speed diagram
		1.2												1.19
														1.18
	1.15
														1.17
	c s /(gd) 1/2	1.1												c s 1.16 /(gd) 1/2
														1.15
	1.05									Tanaka solution Fenton 9th order Classical Peregrine		1.14		Tanaka solution Fenton 9th order Classical Peregrine
											BBM-KdV					BBM-KdV
											iBBM iPeregrine			1.13		iBBM iPeregrine
											Serre					Serre
		1	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.3	0.35	0.4	0.45
								a 0 /d							a 0 /d
		(a) Solitary wave amplitude-speed diagram			(b) Magnification of (a)

Table 1 .

 1 Symmetric head-on collision.

	1.014	Classical Peregrine
		iPeregrine
	1.012
	1.01
	1.008
	K
	1.006
	1.004
	1.002
	1.1	1.15	1.2	1.25	1.3	1.35	1.4
				Speed c s

Acknowledgments

D. Dutykh acknowledges the support from French "Agence Nationale de la Recherche", project "MathOcéan" (Grant ANR-08-BLAN-0301-01) along with the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. The support from the University of Valladolid during the stay of D. Dutykh in July, 2011 is equally acknowledged. A. Duran has been supported by MICINN project MTM2010-19510/MTM. The authors would like to thank Professors Didier Clamond, Vassilios Dougalis and Cesar Palencia for very helpful discussions.


Sun Wuliang

Benoit Combemale

Steven Derrien

Robert B France

Wuliang Sun

Using Model Types to Support Contract-Aware Model Substitutability

Keywords: SLE, Modeling Languages, Model Typing, Contract Matching, Model Substitutability

Model typing brings the benefit associated with well-defined type systems to model-driven development (MDD) through the assignment of specific types to models. In particular, model type systems enable reuse of model manipulation operations (e.g., model transformations), where manipulations defined for models typed by a supertype can be used to manipulate models typed by subtypes. Existing model typing approaches are limited to structural typing defined in terms of object-oriented metamodels (e.g., MOF), in which the only structural (well-formedness) constraints are those that can be expressed directly in metamodeling notations (e.g., multiplicity and element containment constraints). In this paper we describe an extension to model typing that takes into consideration structural invariants, other than those that can be expressed directly in a metamodeling notation, and specifications of behaviors associated with model types. The approach supports contract-aware substitutability, where contracts are defined in terms of invariants and pre-/post-conditions expressed using OCL. Support for behavioral typing paves the way for behavioral substitutability. We also describe a technique to rigorously reason about model type substitutability as supported by contracts, and apply the technique in a usage scenario from the optimizing compiler community.

Introduction

In Model Driven Engineering (MDE), developers of complex software systems create and transform models using model authoring and transformation technologies. The rise in the number of new modeling languages, however, presents a challenge because it requires software engineers to create complex transformations that manipulate models expressed in the new languages. Building these transformations from scratch requires significant effort. To address this problem, various approaches [START_REF] Steel | On model typing[END_REF][2][3][4] [START_REF] De Lara | Generic meta-modelling with concepts, templates and mixin layers[END_REF] have recently been proposed to facilitate the reuse of model transformation across different languages.

Model substitutability rules that are based on model typing [START_REF] Steel | On model typing[END_REF] can be used to support model transformation reuse. For example, a subtyping relation that supports model substitutability allows a model typed by A to be safely used where a model typed by B is expected, where B is the supertype of A. The transformation used for models typed by B can thus be reused on models typed by A.

Current approaches to model type definition and implementation, however, only consider MOF-based metamodels as model types. In MOF, contracts (e.g., pre-conditions, post-conditions and invariants) are externally defined, that is, they are defined in another language, for example, the Object Constraint Language (OCL) [START_REF]OMG: UML Object Constraint Language (OCL) 2.0 Specification[END_REF]. Neither the original paper on model typing [START_REF] Steel | On model typing[END_REF] nor the follow-up paper [START_REF] Guy | On Model Subtyping[END_REF] considers externally defined contracts in subtyping relations. This limits the utility of model subtyping in model-based software development approaches that are contract based (e.g., design by contract [START_REF] Meyer | Applying'design by contract[END_REF]). There is thus a need for model typing that provides support for typing models with contracts.

In this paper we propose a form of model typing that supports contract-aware substitutability, where contracts are defined in terms of invariants and pre-/post-conditions expressed using OCL. We add invariants to model types that specify additional structural properties, and use operation pre-/post-conditions to specify the transformation rules on model types. We also describe a technique for rigorously reasoning about the substitutability of models with contracts.

The rest of the paper is organized as follows. Section 2 illustrates the need for contract-aware substitutability using motivating examples from the high-performance embedded system design domain. Section 3 presents background material needed to understand the work described in this paper. Section 4 presents a formal definition of the subtyping relation between two model types that include contracts, and describes tool support for reasoning about substitutability on model types. Section 5 describes limitations of the approach. Section 6 discusses related work, and Section 7 concludes the paper with a discussion of planned future work.

Motivating Examples

In this section we describe two motivating examples from the high-performance embedded system design domain. Modern heterogeneous embedded hardware platforms are notoriously difficult to design and to program. In this context, tool-supported model based approaches (e.g., Simulink, Ptolemy) are now widely acknowledged as some of the most effective approaches to designing embedded systems.

Typically, these model-based approaches use tool chains that manipulate many different types of models. For example, structural platform description models range from system level models that abstract over processing and storage resource with their interconnections, to very low level Register-to-Logic level circuit models that are used to describe the structure of hardware accelerators within the platforms.

Similarly, behavioral description models range from application level modeling of the application using Models of Computation such as Synchronous Data Flow Graphs or Kahn Process Networks, to fine grain scalar operation level representations such as the basic-block level instruction dependence graph used in an optimizing compiler back-end.

Most of these tool chains share a common goal: They aim to produce highly optimized implementations. This requires the use of advanced algorithms that implement very complex model manipulations. It is also the case that these manipulations often have similar algorithmic patterns. These patterns can be used as the basis for developing reusable model transformations.

Example 1: Using Model Types to Support Structural Substitutability

In the optimizing compiler domain, a variety of models describing different aspects of languages are manipulated (i.e., analyzed and transformed) at different stages of the compilation process. While the analyses and transformations are different, they also share many common characteristics. For example, consider algorithms for schedule optimization. Obtaining an optimized implementation of an application on a target platform involves performing several static scheduling optimizations. Many of these algorithms have common characteristics, for example, they are often expressed as an acyclic graph resource constrained scheduling problem, for which many techniques (heuristics or MILP-Mixed Integer Linear Programming-solver based) have been proposed. Because these scheduling algorithms involve very sophisticated algorithms, reusable algorithms that can be tailored to the different types of representations (models) are highly desirable. For example, it would be useful to have a reusable scheduling algorithm that can be used to derive a schedule for an Application level Synchronous Data-flow graph on a multi-processor based implementation, as well as for generating efficient code for a customized VLIW (Very Long Instruction Word) embedded processor.

However in this case, structural substitutability based only on constraints that can be expressed directly in a metamodel (e.g., multiplicity or element containment constraint) is not sufficient; other structural constraints need to be specified. For example, a classical static scheduling toolset can only operate on acyclic dependence graphs and the acyclicity property cannot be expressed directly in a metamodel. A language such as the Object Constraint Language (OCL) is needed to specify properties of acyclic graphs. In this case, model substitutability requires that a substitute model enforces the acyclicity constraint expressed in OCL. Model typing based on metamodels with OCL constraints can be used to enable such structural substitutability.

Example 2: Using Model Types to Support Contract-based Behavioral Substitutability

Behavioral substitutability for model transformations: A consistent scheduling transformation must ensure that every node in the dependence graph is scheduled at least once. This property can be expressed as a post-condition on the scheduling transformation and thus any scheduler implementation should enforce this post-condition. The effective post-condition could even be stricter; in our case we could consider a postcondition that restricts a node to be scheduled exactly once. The same holds for the pre-condition. For example, most schedulers operate on acyclic graphs and this can be translated as a pre-condition for the transformation. However, there also exists a class of pipelined schedulers that operate on cyclic graphs, in which cycles implement delays to preserve causality. For such pipelined schedulers, the pre-condition would not forbid cycles in the dependence graph. That would, however, prevent a pipelined scheduler from being used to schedule acyclic graphs in a design flow. Contract based tool chain validation: An optimizing compiler custom tool chain consists of a sequence of analyses and transformations (called compiler passes) executed in a very carefully chosen order. They can hence be seen as a model transformation chain. Compiler passes cannot be combined arbitrarily, as each pass usually assumes that the program representation at hand has very specific properties.

For example, consider a compiler tool chain for generating software code from Synchronous Data Flow Graph (SDF) model specifications on an embedded platform. Such a tool chain is illustrated in Figure 1. Before any code can be produced, the SDF first needs to be scheduled on this platform. Depending on whether the target system consists of a single or several processors, it is likely that different scheduling algorithms will be used. Similarly, different code generators (i.e. pretty printers) will have to be used depending on whether we target a mono-processor or multiprocessor. Two backend code generators are shown in Figure 1. The mono-processor code generator can only be used after a mono-processor scheduling stage, whereas the second back-end is more general and can be used for both types of scheduling.

These constraints -that is targeting one or several processing resources -apply to the result of the scheduling stage and to the input on the code generation stage. They can hence be modeled as pre-conditions (resp. post-conditions) expressed using OCL. When chaining a given scheduling and code generation pass, we can ensure the consistency of the flow by checking if the pre-/post-conditions of two chained transformation are satisfiable.

In this section, we describe the concepts underlying our use of model types to support model substitutability. We first present the MOF (Meta-Object Facility) meta-language, the basis for metamodels, and thus model manipulation operators. We then give an overview of model types as currently defined and implemented [START_REF] Steel | On model typing[END_REF][START_REF] Guy | On Model Subtyping[END_REF], and describe the limitations addressed by the approach presented in this paper.

Metamodeling

The Meta-Object Facility (MOF) [START_REF]OMG: Meta Object Facility (MOF) 2.0 Core Specification[END_REF] is the OMG's standardized meta-language, i.e., a language to define metamodels. As such, it is a common basis for a vast majority of modeling languages and tools. A metamodel defines a set of models on which it is possible to apply common operators. The model substitutability approach presented in this paper is applicable to models expressed in languages with MOF metamodels.

MOF supports the definition of metamodels using Classes and Properties. Classes can be abstract (i.e., they cannot be instantiated) and have Properties and Operations, which respectively declare attributes and references, and the signatures of methods available to the modeled concept. A Property can be composite (an object can only be referenced through one composite Property at a given instant), derived (i.e., calculated from other Properties) and read-only (i.e., cannot be modified). A Property can also have an opposite Property with which it forms a bidirectional association.

Metamodels can be viewed as class diagrams in which each metamodel element can be instantiated to obtain objects representing model elements. However, metamodel elements are themeslves instances of MOF elements and thus a metamodel can be drawn as an object diagram where each concept is an instance of one of the MOF elements (e.g., Class or Property classes).

Model Typing

Model Types were introduced by Steel et al. [START_REF] Steel | On model typing[END_REF], as an extension of object typing to provide abstractions about the object type level and enable the reuse of model manipulation operators. Informally, a model type is a substructure (referred to as a type group) of the metamodel's class diagram. It is important to distinguish the usage of the term metamodel from model type. We use the term metamodel to refer to the class diagram used to define a language, and when the same class diagram is used to define the type of a model it is called an exact type. It is also important to note that a model has one and only one metamodel to which it must conform, but the same model can have several model types, where each model type is a substructure of the metamodel. Because model types and metamodels share the same structure, it is possible to extract the exact type of a model from its metamodel. Figure 2 represents a model m 1 that conforms to a metamodel MM 1 and is typed by model types MT A and MT B , where MT B is the exact type of m 1 that is extracted from MM 1 . Both metamodels and model types conform to MOF. Given the above, a model type can be defined as follows: Substitutability is the ability to safely use an object of type A where an object of type B is expected. Substitutability is supported through subtyping in object-oriented languages. However, object subtyping does not handle specializations of model substructures (or type groups) 1 . One way to safely reuse a model manipulation operation created for a model typed by MT A on a model typed by MT B is to ensure that MT A contains elements that can be substituted by elements defined by MT B . However, it is not possible to achieve model type substitutability through object subtyping. Thus, model typing uses an extended definition of object type matching introduced by Bruce et al. [START_REF] Bruce | Polytoil: A type-safe polymorphic objectoriented language[END_REF], namely MOF Class Matching.

(MT B , MT A ) ∈ (also written MT B MT A ) iff ∀ T A ∈ MT A , ∃ T B ∈ MT B such that T B <# T A .
We recently introduced four extended subtyping relations between model types that take into account two additional criteria: The presence of heterogeneities between two model types (using adaptation) and the considered subset of the model types (using model type pruning) [START_REF] Guy | On Model Subtyping[END_REF].

The subtyping relation as currently defined has shortcomings. In particular, the current model typing definition and implementation only considers MOF-based metamodels as model types (through the MOF class matching relation). Unfortunately, MOF delegates the definitions of contracts (e.g., pre and post-conditions or invariants) to other languages (e.g., OCL, the Object Constraint Language [START_REF]OMG: UML Object Constraint Language (OCL) 2.0 Specification[END_REF]). This limits the applicability of model typing for safely reusing model manipulations where OCL contracts are needed to precisely specify the applicability of the model transformation or the structure on which the model transformation can be applied (see motivating examples in Section 2). The approach described in this paper addresses this limitation.

Contribution

In this paper we extend the subtyping relation described in [START_REF] Guy | On Model Subtyping[END_REF] by taking into account OCL contracts for a safe substitutability of models conforming to metamodels including contracts. This provides a safe reuse of model transformations expressed on metamodels that include contracts. Specifically, we extend the MOF class matching (cf. Def. 2 of Section 3) by considering contracts matching (Section 4.1) and provide a technique for analyzing the matching of OCL contracts associated with two classes with different model types (Section 4.2). In this section we describe the contract matching technique we developed to support contract-aware model substitutability. We also describe an Alloy-based prototype tool that supports contract matching (Section 4.3), and illustrate the use of contract-aware substitutability using the motivating examples (Section 4.4).

Contract-aware MOF Class Matching

We consider the use of OCL invariants added to MOF classes to specify additional structural properties, and OCL pre-/post-conditions defined in the context of MOF class operations to specify the model manipulation rules (e.g., transformation) associated with model types. The MOF class matching relation is thus determined by two aspects: the structural features specified using MOF (e.g., classes, properties, operation signatures, etc.) and the contracts expressed using OCL (e.g., invariants and pre-/post-conditions).

The substitutability through model subtyping is a specialization of the Liskov Substitution Principle [START_REF] Liskov | A behavioral notion of subtyping[END_REF] on the model type system. Specifically the contract matching relation that enables contract-aware model substitutability must abide by the following rules: (1) invariants of the supermodel type cannot be weakened in a sub model type, (2) pre-conditions cannot be strengthened in a sub model type, and (3) post-conditions cannot be weakened in a sub model type. The extended MOF class matching relation is formalized as follows:

Definition 4 (Contract-aware MOF Class Matching). Class T matches T (written T <# T) iff their structures match (cf. Def. 3 of [START_REF] Guy | On Model Subtyping[END_REF]), their invariants match and their operation pre-/post-conditions match, where The invariants of T and T match if Models(invs) ⊇ Models(invs ), where Models(invs) returns all models that satisfy invs and Models(invs ) returns all models that satisfy invs .

2 Pre-/post-conditions Match is defined as follows:

∀ op ∈ T.ownedOperation, ∃ S ∈ SuperClasses(T ) such that ∃ op ∈ S .ownedOperation and:

2.1 let op.ownedPrecondition = {pre 1 , pre 2 , ..., pre k } be the pre-conditions defined for op; let pres = pre The approach uses the Alloy Analyzer at the back-end to analyze the relation between two contracts, and it thus requires a translation from OCL expressions to Alloy specifications. The OCL to Alloy translation used in the prototype tool we developed is based on translation rules described in work by Bordbar et al. [START_REF] Anastasakis | On challenges of model transformation from uml to alloy[END_REF].

Contract Matching Checking Tool

The contract matching approach described in the previous subsection has been implemented in a prototype tool. Figure 3 shows an overview of the prototype tool. It consists of an OCL parser, an Ecore2 /OCL transformer and the use of the Alloy Analyzer. The Ecore/OCL transformer is developed using Kermeta [START_REF] Muller | Weaving executability into object-oriented metalanguages[END_REF], an aspect-oriented metamodeling tool. The inputs of the prototype are (1) an Ecore file that specifies two model types, and (2) a textual OCL file that specifies the contracts from each model type. The model types and contracts are automatically transformed to an Alloy model consisting of signatures and predicates.

The prototype provides several interfaces to check contract matching. For example, matchInv(inv1: Constraint, inv2: Constraint) is used to check whether inv1 matches inv2. In addition, matchInvs(cls1: Class, cls2: Class) can be used to check whether the invariants defined in cls1 and the invariants defined in cls2 match.

Case Study

In this section we illustrate how to use our approach to define model types and subtyping relations between them to ensure a safe reuse of model transformations.

A Simple Case Study of Structural Substitutability

Let us reconsider the scheduling example described in Section 2.1. A model transformation performs a static scheduling on an acyclic dependence graph. The model transformation needs a metamodel for "Acyclic Graph" (due to space limitation, the metamodel is not shown in the paper). The model type AcyclicGraph (see Figure 4) shows a simple example of model type definition for the dependency graph used in the example. Its definition consists of metaclasses that specify a graph structure, an invariant that specifies the acyclicity property, and a model transformation that takes as input an acyclic graph. A Simple Case Study of Behavioral Substitutability In the optimizing compiler community, the daily task for software engineers is to design compilation chains in the right partial order, that is, scheduling the various passes (i.e., optimization, translation, code generation, analysis, etc.). Designing compilation chains would benefit from behavioral substitutability by opening the way to describe "abstract" compilation chains, capitalizing a given knowledge in terms of constraints (pre-/post-conditions) to schedule a set of passes for a given purpose, where each pass would be then implemented in various ways, but conforming to the pre-/post-conditions defined in the abstract compilation chain.

Figure 5 shows a simple example of model types used for the compilation chain. Suppose that the abstract model transformation transfo defined in MT is used for optimization purpose and define a post condition stating that the model must conform to the Static Single Assignment (SSA) form. MT also contains transfo2 as the next pass of the compilation chain and states as precondition that the model must conform to the SSA form. The two model types subMT 1 and subMT 2 implement the model transformation transfo but only subMT 1 ensures as postcondition the SSA form. While subMT 2 is not Fig. 5. A Simple Example of Behavioral Substitutability in Kermeta in this case a sub model type to MT , a compilation error for m2.transfo().transfo2() will be returned. This shows that the model returned by transfo in subMT 2 (typed by subMT 2) is not of type MT , and can not reuse transfo2.

Discussion

In this section we discuss limitations of our work, and its scope of application. We first discuss the supported contracts in the subtyping relation of model typing (Section 5.1), and the corresponding model substitutability provided by our approach (Section 5.2).

On the Support of Contracts in Model Typing

In this paper we consider contracts in addition to the object oriented structure described in a metamodel. The object-oriented structure is usually defined using Ecore, an implementation aligned with OMG MOF. Contracts can then be invariants expressed in the context of the concepts (i.e., classes) defined in the MOF metamodel, and pre-/postconditions expressed in the context of operations specified in concepts. While invariants restrict the structure of conforming models and their possible structural substitutability, pre-and post-conditions specify the behavior of the conforming models (i.e. manipulation by model operations) and their possible behavioral substitutability.

In our approach, we assume that the first order logic is used to express contracts in metamodels, and we have chosen OCL to express them. We rely on the provided binding between MOF and OCL as defined by OMG to link OCL expressions to a given MOF metamodel.

To test the feasibility of our approach, we implemented a prototype tool that is integrated into the Kermeta workbench. The tool checks OCL-based contract-aware subtyping relations between model types. While the substitutability related to the MOF structure is computed directly using Kermeta, the one related to the contracts is computed using Alloy through a translation from OCL expressions to Alloy specifications, and then an analysis of the output provided by the Alloy Analyzer. The tool only provides support for translating a subset of OCL to Alloy.

Most of the OCL operators have corresponding Alloy constructs. For example, OCL operator f orAll corresponds to Alloy construct all, exists corresponds to some, includes corresponds to in, excludes corresponds to !in, sum corresponds to sum, and closure corresponds to * . OCL contracts that involves such operators can be directly transformed into Alloy specifications.

However, as pointed out by Anastasakis et al. [START_REF] Anastasakis | On challenges of model transformation from uml to alloy[END_REF], the translation from OCL to Alloy is not seamless. There are some OCL operators that do not have corresponding Alloy constructs, and thus OCL contracts including such operators cannot be easily transformed into Alloy specifications. Some of them can be partially supported by the tool using the Alloy libraries. For instance, OCL operators like select and collect are translated by the tool described in the paper using Alloy functions that implement their semantics. Consequently, the operator iterate is partially supported by the transformation tool. The tool provides support for OCL contracts including iterate expressions that can be rewritten as f orall with select/collect operators. However, the tool cannot be used to deal with iterate expressions that involve arithmetic accumulation since Alloy is a purely declarative language that does not provide support for imperative accumulators. Finally, the translation cannot deal with OCL casting operators like oclAsType since Alloy has a very simple type system that has little support for type casting.

On the Support of Modeling Language Substitutability

The research work described in the paper builds upon our previous work in [START_REF] Guy | On Model Subtyping[END_REF], and paves the way for reasoning about the subtyping relation between two model types that include contracts. Specifically it can be used to reason about the contract-aware subtyping relation that involves structural subtyping (including not only MOF-based Object-Oriented structure but also OCL-based first order invariants) and behavioral subtyping (including a behavioral semantics in an axiomatic way using pre-/post-conditions on operations).

We implement our approach in a (Kermeta-based) tool included in the Kermeta language workbench to check advanced (i.e., including contracts) subtyping relations between modeling languages based on Ecore and OCL. These two meta-languages are supported by the Kermeta language workbench and are used for describing the abstract syntax and the static semantics respectively. This approach and its corresponding implementation addresses the need illustrated by the motivating examples from the high-performance embedded system community used throughout the paper. The scope of the structural substitutability we offer is bounded by OCL and its translation to Alloy, and its applicability is well founded, e.g., in model transformation reuse in model-driven development [START_REF] Guy | On Model Subtyping[END_REF].

The actual scope of behavioral substitutability is more difficult to define. The difficulty is twofold: while we support the motivating examples described in the paper, complex situations of OCL based typing could be considered, such as type propagation in model transformation chains. Such challenges will be addressed in future work. Moreover, the scope itself of behavioral substitutability is more difficult to delimit.

Related Work

The technical contribution of this paper is the integration of contract matching in the subtyping relation of model typing to enhance the substitutability supported between modeling languages. As discuss in the previous section, we rely for that on the most established translation to Alloy. Then, the work related to our contribution discussed in this paper is the applicability of the substitutability as illustrated in the motivating examples, namely on model transformation reuse.

Substitutability is supported through subtyping in object-oriented languages, including the support of contracts (e.g., Eiffel [START_REF] Meyer | Design by contract. the eiffel method[END_REF]). However, object subtyping does not handle type group specialization (i.e., the possibility to specialize relations between several objects and thus groups of types) 3 . Such type group specialization have been explored by Kühne in the context of MDE [START_REF] Kühne | On model compatibility with referees and contexts[END_REF]. Kühne defines three model specialization relations (specification import, conceptual containment and subtyping) implying different level of compatibility. We are only interested here in the third one, subtyping, which requires an uncompromised mutator forward-compatibility, e.g., substitutability, between instances of model types.

Several approaches have been proposed during the last decade for model transformation reuse. Strict substitutability relation, such as the first version of model type matching presented in [START_REF] Steel | On model typing[END_REF], offers the possibility to reuse model transformation through isomorphic metamodels, i.e., metamodels with MOF-based equivalent structures. Such possibility was first proposed in [START_REF] Varró | Generic and meta-transformations for model transformation engineering[END_REF] where the authors introduce variable entities in patterns for declarative transformation rules. These entities express only the needed concepts (e.g., types, attributes, etc.) to apply the rule, allowing any tokens with these concepts to match the pattern and thus to be processed by the rule. Latter, Cuccuru et al. introduced the notion of semantic variation points in metamodels [START_REF] Cuccuru | Templatable metamodels for semantic variation points[END_REF]. Variation points are specified through abstract classes, defining a template, and metamodels can fix these variation points by binding them to classes extending the abstract classes. Patterns containing variable entities and templates can be seen as kinds of model types where the variability has to be explicitly expressed and thus anticipated. Sanchez Cuadrado et al.

propose in [START_REF] Sánchez Cuadrado | Approaches for model transformation reuse: Factorization and composition[END_REF] a notion of substitutability based on model typing and model type matching, but rather to use an automatic algorithm to check the matching between two model types, they propose a DSL that allows users to declare the matching by hand. Finally, De Lara et al. present in [START_REF] De Lara | Generic meta-modelling with concepts, templates and mixin layers[END_REF] the concept mechanism, along with model templates and mixin layers leveraged from generic programming to MDE. Concepts are really close to model types as they define the requirements that a metamodel must fulfill for its models to be processed by a transformation, under the form of a set of classes. The authors also propose a DSL to bind a metamodel to a concept and a mechanism to generate a specific transformation from the binding and the generic transformation defined on the concept.

In the context of model transformation chains, existing approaches deal with explicit relationships between model transformations. Vanhooff et al. [START_REF] Vanhooff | Uniti: A unified transformation infrastructure[END_REF] proposed a domain specific language to model and execute a transformation chain. Aranega et al. [START_REF] Aranega | Using feature model to build model transformation chains[END_REF] used feature models to classify model transformations involved a transformation chain and specified the constraints between them. The user thus can design a new transformation chain by reusing the classified transformations. Yie et al. [START_REF] Yie | Realizing model transformation chain interoperability[END_REF] advocated the use of several independently developed model transformation chains to convert a high-level model into a low-level model. The interoperability among model transformation chains is achieved by deriving correspondence relationships between the final models generated by each model transformation chain.

Unlike the above approaches, contract-aware model subtyping offers a unified and formal type theory to facilitate the safe reuse of model transformations involved in a transformation chain. It follows a declarative fashion to specify a model transformation chain in an abstract way using pre-/post-conditions on abstract model types. This promotes then the reuse of various implementations that match the conditions for a safe execution of the model transformation chain.

Conclusion and Perspective

We propose in this paper a model typing theory where model types include contracts. This includes a formally defined subtyping relation between model types, and a toolsupported approach supporting a safe contract-aware substitutability of models conforming to metamodels including contracts. This ensures a safe reuse of model transformations expressed on metamodels including contracts.

Contracts are defined in terms of invariants and pre-/post-conditions expressed using OCL on MOF-based metamodels. The invariants are added on the classes of a metamodel to specify additional structural properties of the metamodel, and pre-/postconditions are added on the operations of classes to specify model transformations. Consequently, the support of invariants in the subtyping relation ensures a safe reuse of model transformations where OCL contracts are needed to precisely specify the structure on which the model transformation can be applied. The support of pre-/postconditions paves the way for behavioral substitutability to safely reuse model transformations where OCL contracts are needed to precisely specify the applicability of the model transformation.

The subtyping relation is based on a matching relation between two MOF classes that include OCL contracts, and is checked thanks to a technique based on Alloy. The actual scope of the provided contract-aware substitutability is mainly determined by the OCL-to-Alloy translation.

We are currently extending the prototype by providing support for model types and contracts expressed using the Kermeta language workbench. We also explore how we can extend the approach by using SMT solvers at back-end to analyze the OCL contracts that include more complex arithmetic calculation.

Fig. 1 .

 1 Fig. 1. A Model based compiler tool chain for embedded multiprocessors.

Fig. 2 .

 2 Fig. 2. Conformance, model typing and model subtyping relations

Definition 2 .

 2 (MOF class matching) MOF class T matches T (written T <# T ) iff their names are equal, and for each property (respectively method) in T there is a corresponding property (respectively method) in T .The MOF class matching relation can be seen as a kind of object type matching relation that is tailored to MOF concepts. Based on the MOF class matching relation, we can achieve model type substitutability by defining a subtyping relation as follows: Definition 3. (Subtyping relationship for model types) The model type subtyping relation is a binary relation on ModelType, the set of all model types, such that

Fig. 3 .

 3 Fig. 3. Contract Matching Checking Tool Overview

Fig. 4 .

 4 Fig. 4. A Simple Example of Structural Substitutability in Kermeta

1

  Invariants Match is defined as follows: let T.ownedInvariant = {inv T 1 , inv T 2 , ..., inv T k } be the invariants defined for T ; let result T = inv T 1 ∧ inv T 2 ∧ ... ∧ inv T k ; let SuperClass(T) = {cls 1 , cls 2 , ..., cls n } where cls i is a superclass of T; let cls i .ownedInvariant = { inv i1 , inv i2 , ..., inv ik } be the invariants defined for cls i , for i = 1, .., n; let result i = inv i1 ∧ inv i2 ∧ ... ∧ inv ik , for i = 1, .., n; let invs = result 1 ∧ result 2 ∧ ... ∧ result n ∧ result T ; let T .ownedInvariant = {inv T 1 , inv T 2 , ..., inv T k } be the invariants defined for T ; let result T = inv T 1 ∧ inv T 2 ∧ ... ∧ inv T k ; let SuperClass(T ) = {cls 1 , cls 2 , ..., cls n } where cls i is a superclass of T ; let cls i .ownedInvariant = { inv i1 , inv i2 , ... , inv ik } be the invariants defined for cls i , for i = 1, .., n; let result i = inv i1 ∧ inv i2 ∧ ... ∧ inv ik , for i = 1, .., n; let invs = result 1 ∧ result 2 ∧ ... ∧ result n ∧ result T ;

  one to formally analyze the relation between contracts (e.g., to check if the models satisfying one contract includes the models satisfying the other). To do this, a query function Models(MT , C) is used to compute all models that both conform to a model type MT and satisfy an OCL contract C defined in MT . Thus given contract C 1 in a candidate supermodel type MT 1 and contract C 2 in a candidate sub model type MT 2 , C 1 matches C 2 iff (1) C 1 , C 2 are invariants, and Models(MT 1 , C 1 ) ⊇ Models(MT 2 , C 2 ), (2) C 1 , C 2 are pre-conditions, and Models(MT 2 , C 2 ) ⊇ Models(MT 1 , C 1 ), and (3) C 1 , C 2 are post-conditions, and Models(MT 1 , C 1 ) ⊇ Models(MT 2 , C 2 ).Checking the contract matching requires a tool to implement the functionality of the query function Models(MT , C). We use the Alloy Analyzer[START_REF] Jackson | Alloy: a lightweight object modelling notation[END_REF] for this purpose. The Alloy Analyzer is used to analyze Alloy specifications. It is supported by a SAT-based model finder. The Alloy Analyzer can generate models that conform to a model type expressed in Alloy in terms of signatures and fields that specify the model type structure and a predicate that expresses the contracts. In this paper we use the Alloy Analyzer at the back-end to check whether two contracts match. For example, given a candidate supermodel type MT 1 and a candidate sub model type MT 2 , with two OCL invariants respectively, C 1 and C 2 , the procedure below can be used to check if C 1 matches C 2 . 1. (preprocess) Since model subtyping requires each element in the supermodel type to be matched by an element in the sub model type (see Definition 4), the contract defined in the supermodel type refers to elements that also exist in the sub model type. Thus we can move C 1 to MT 2 , and use only the sub model type (i.e., MT 2 ) to check whether C 1 and C 2 match. 2. Transform MT 2 to an Alloy model using the technique described in [14]. Convert C 1 and C 2 into two Alloy predicates, P 1 and P 2 , respectively. 3. Run an empty predicate in the Alloy Analyzer to search for a model conforming to the model type MT 2 . If the Analyzer returns no model satisfying the empty predi-Run P 1 and P 2 respectively. If the Alloy Analyzer returns no model for each predicate (i.e., Models(MT 2 , C 1 ) = / 0 and Models(MT 2 , C 2 ) = / 0), then C 1 matches C 2 ; if the Alloy Analyzer returns a model (or models) for only P 1 , then C 1 matches C 2 ; if the Alloy Analyzer returns a model (or models) for only P 2 , then C 1 does not match C 2 ; otherwise, continue to the next step. 5. Run a predicate to search for a model satisfying both P 1 and P 2 . If the Alloy Analyzer returns a model satisfying the predicate, continue to the next step; otherwise, C 1 does not match C 2 . 6. Run a predicate P 3 to search for a model satisfying both P 1 and ¬P 2 (i.e., the negation of P 2 ), and another predicate P 4 to search for a model satisfying both P 2 and ¬P 1 . If the Alloy Analyzer returns no model for both P 3 and P 4 (i.e., Models(MT 2 , C 1 ) = Models(MT 2 , C 2 )), C 1 matches C 2 ; if the Alloy Analyzer returns a model (or models) satisfying only P 3 , Models(MT 2 , C 1 ) ⊃ Models(MT 2 , C 2 ) and C 1 matches C 2 ; otherwise, C 1 does not match C 2 .

	such relations requires cate (i.e., Models(MT 2 , / 0) = / 0), Models(MT 2 , C 1 ) = / 0 and Models(MT 2 , C 2 ) = / 0. In this case C 1 matches C 2 since / 0 is a subset of / 0; otherwise, continue to the next step. 4. Preprocess Model Types with OCL Contracts, and subtyping relation Sub Model Types and OCL Contracts Result of the Eclipse Ecore + OCL Alloy subtyping OCL Parser Transformer relation Analyser checking
	The operation specifications of T and T match if Models(pres ) ⊇ Models(pres)
	and Models(posts) ⊇ Models(posts )
	4.2 Analyzing the Matching of Contracts
	Definition 4 can be used to formally reason about the matching relation between two
	MOF classes with contracts. The MOF class matching relation in Definition 4 includes
	the matching of the contracts from classes of two model types. Consequently, analyzing

1 ∧ pre 2 ∧ ... ∧ pre k ; let op .ownedPrecondition = {pre 1 , pre 2 , ..., pre k } be the pre-conditions defined for op ; let pres = pre 1 ∧ pre 2 ∧ ... ∧ pre k ; 2.2 let op.ownedPostcondition = {post 1 , post 2 , ..., post k } be the post-conditions defined for op; let posts = post 1 ∧ post 2 ∧ ... ∧ post k ; let op .ownedPostcondition = {post 1 , post 2 , ..., post k } be the post-conditions defined for op ; let posts = post 1 ∧ post 2 ∧ ... ∧ post k ;

For further information on type groups see Ernst's paper[START_REF] Ernst | Family polymorphism[END_REF].

Ecore is an implementation aligned with MOF included in the Eclipse Modeling Framework[START_REF] Steinberg | EMF: Eclipse Modeling Framework[END_REF].

We refer the reader interested in the type group specialization problem to the Ernst's paper[START_REF] Ernst | Family polymorphism[END_REF].

Acknowledgment This work was supported by the National Science Foundation grant (CCF-1018711), the ANR INS Project GEMOC (ANR-12-INSE-0011), and the CNRS PICS Project MBSAR.