Didier Clamond
email: didierc@unice.fr

New exact relations for easy recovery of steady wave profiles from bottom pressure measurements

Keywords:

This paper is about the determination of a free surface wave profile from a given pressure distribution at the bottom. For a two-dimensional irrotational steady surface water wave propagating over a flat horizontal bed, this problem is solved analytically in an implicit form. Explicit solutions can be easily obtained numerically via fixed point iterations, whose convergence is outlined mathematically and numerically. This new surface reconstruction procedure does not involve the resolution of a differential or an integral equation, i.e. the surface is given by algebraic local relations. Thus, this formulation permits the free surface recovery in a simpler and more efficient way than previous methods.

Introduction

Within the framework of incompressible, homogeneous and inviscid flow modelling of water waves, the pressure plays an important role in establishing various qualitative properties of travelling waves. The pressure is essential in quantitative studies because often in practice the elevation of a surface water wave is determined from pressure data obtained at the sea bed (see, e.g., [START_REF] Tsai | On the recovery of surface wave by pressure transfer function[END_REF]. Indeed, measuring a free surface wave motion without intrusive probes is of great importance for practical applications, and one way to achieve this is to use pressure gauges at the sea bed. From these pressure measurements, the question is then how can one reconstruct the free surface. This problem has long been of interest and has recently attracted attention (e.g. [START_REF] Constantin | On the recovery of solitary wave profiles from pressure measurements[END_REF]Deconinck, Oliveras & Vasan 2012;[START_REF] Oliveras | Recovering surface elevation from pressure data[END_REF][START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF].

The first approaches to reconstructing the free surface from pressure measurements consisted of assuming that the hydrostatic approximation is sufficiently accurate. However, for waves of moderate amplitude, prediction errors are significant [START_REF] Bishop | Measuring waves with pressure transducers[END_REF]. These inaccuracies are due to the fact that the hydrostatic approximation ignores even linear wave effects. Taking into account the wave effect to the linear approximation, one can derive a better approximation [START_REF] Escher | On the recovery of the free surface from the pressure within periodic traveling water waves[END_REF], but the benefit is limited for waves of moderate amplitude [START_REF] Tsai | On the recovery of surface wave by pressure transfer function[END_REF].

These considerations motivated the quest for a reconstruction formula that accounts for all nonlinear effects and that is thus applicable to waves of moderate and large amplitude. Nonlinear non-local equations relating the dynamic pressure on the bed and the wave profile were obtained without approximation from the governing equations [START_REF] Constantin | On the recovery of solitary wave profiles from pressure measurements[END_REF][START_REF] Oliveras | Recovering surface elevation from pressure data[END_REF].

Very recently, other exact tractable relations have been obtained [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF], allowing mathematical analysis and straightforward numerical procedures for deriving the free surface from the pressure at the bed. This surface reconstruction method requires the resolution of an ordinary differential equation. Although the numerical resolution of this differential equation is numerically easy, it is actually possible to solve it analytically in an implicit form, thus providing an even simpler reconstruction procedure. The purpose of the present paper is to show how and to provide mathematical insights of the efficiency of this new approach, as well as some numerical evidence.

This paper is organized as follows. In § 2, the hypothesis and notation are introduced for steady irrotational surface waves on constant depth. In order to make the present paper self-contained, the main results of [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF] are given in § 3. A new, purely local, expression for the free surface is then derived in § 4. It requires the introduction of a new holomorphic which can be easily calculated, as demonstrated in § 5. The convergence properties of the new formula for the surface are outlined in § 6. Finally, numerical examples are provided in § 7.

Hypothesis and notation

The physical hypothesis and most mathematical notation are identical to those of [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. For the sake of completeness, they are briefly recalled here.

We consider a steady irrotational wave motion due to a surface gravity wave of permanent form. Let (x, y) be a Cartesian coordinate system moving with the wave, x being the horizontal coordinate and y the upward vertical coordinate. The wave is (2π/k)-periodic in the x direction, solitary waves corresponding to the limiting case k → 0. Here y = -d, y = η(x) and y = 0 denote the equations of the horizontal impermeable bottom, of the impermeable free surface and of the mean water level, respectively. The latter equation expresses the fact that η = 0 for the smooth (2π/k)-periodic wave profile η, where • is the Eulerian average operator over one wavelength, i.e.

η ≡ k 2π π/k -π/k η(x) dx = 0. (2.1)
Let u, v, φ and ψ be the horizontal velocity, the vertical velocity, the velocity potential and the stream function, respectively, such that u = φ x = ψ y and v = φ y = -ψ x . Also let the complex potential f = φ + iψ and the complex velocity w = uiv (with i 2 = -1) be holomorphic functions of the complex abscissa z = x + iy, i.e. f = f (z) and w = df /dz. The mass conservation and the irrotational character of the flow are identically fulfilled with the requirement that f and w are holomorphic functions throughout the fluid domain. The flow being steady, the bottom and surface boundaries being impermeable, the stream function ψ = Im{f } is constant at y = -d and at y = η. Thus, v s = η x u s (with η x = dη/dx) and v b = 0, where subscripts 's' and 'b' denote the restrictions at the free surface and at the bed, respectively.

The momentum equation can be expressed by means of the Bernoulli condition

2p + 2gy + u 2 + v 2 = B, x ∈ R, -d y η, (2.2)
for some Bernoulli constant B, and where p = p(x, y) is the pressure divided by the (constant) density. At the free surface, the surface tensions are neglected and the pressure is zero, i.e. p s = 0. From (2.1) and (2.2), we get relations for the Bernoulli constant We also define the wave phase velocity c such that c = -u b , so that the wave travels with phase speed c in the frame of reference where the mean horizontal velocity is zero at the bed, and where c > 0 if the wave travels toward the increasing x direction. For smooth waves, the inequality u/c < 0 occurs everywhere in the fluid and at the free surface, a stagnation point (where u = v = 0) being possible only at the angular crests of the highest waves. Note that B = c 2 in deep water (d → ∞) and for solitary waves (k → 0), but B = c 2 in general.

B = u 2 s + v 2 s = u 2 b > 0, ( 2

Principle of the surface recovery

Following [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF], let be the holomorphic function P such that

P(z) ≡ 1 2 B + gd -1 2 w 2 (z) = 1 2 B + gd -1 2 (u 2 -v 2 ) + iuv, (3.1)
having thus restriction at the flat bed y = -d with zero imaginary part and real part p, i.e. p(x) = P b = P(xid) = gd + (Bu 2 b )/2, thus P(z) = p(z + id). Here p coincides with the real function Re{P} only on y = -d because p is not a harmonic function [START_REF] Constantin | The trajectories of particles in Stokes waves[END_REF][START_REF] Constantin | Pressure beneath a Stokes wave[END_REF].

Using the surface impermeability and (2.2) on the free surface y = η(x), where p = 0 and v s = u s η x , we have

(u s -iv s ) 2 = (1 -iη x ) 2 u 2 s = (1 + η 2 x )u 2 s (1 -iη x )/(1 + iη x ) = (u 2 s + v 2 s )(1 -iη x )/(1 + iη x ) = (B -2gη)(1 -iη x )/(1 + iη x ). (3.
2) Multiplying this relation by (1 + iη x ) and using (3.1), we obtain at once

gη(1 -iη x ) + iBη x = (P s -gd)(1 + iη x ), (3.3)
with P s = P(z s ) = P(x + iη(x)). The real and imaginary parts of (3.3) give the two equations

gη = Re{P s } -gd -η x Im{P s } = 1 2 (B -u 2 s -v 2 s ), (3.4) (B -gη)η x = (Re{P s } -gd)η x + Im{P s } = 1 2 (B + u 2 s + v 2 s )η x (3.5)
where the second equalities are obtained using (3.1) and v s = η x u s .

The wave amplitude a is obtained by evaluating (3.4) at the wave crest located at x = 0 (where η x = 0 and η = a 0) leading to the implicit equation

a = g -1 Re{P(ia)} -d. (3.6)
Similarly, at the wave trough located at x = π/k (where η x = 0 and η = -b 0, a + b being thus the total wave height), we have the corresponding implicit equation

b = d -g -1 Re{P(π/k -ib)}. (3.7)
The crest height a is obtained as the unique solution to (3.6), as proven by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. The convergence of (3.7) is discussed in § 6 below.

With the wave amplitude a determined by solving (3.6) iteratively, η is subsequently obtained by re-expressing (3.5) as the ordinary differential equation

η x = Im{P s }/ (B -gη -Re{P s } + gd) , (3.8)
with initial data η(0) = a. The right-hand side of (3.8) being smooth, the solution can be obtained by a standard iterative procedure, as proven and illustrated by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. It is however possible to derive a local simpler expression which does not require the resolution of the differential equation. In other words, the differential equation (3.8) can be solved analytically in an implicit form, leading to an algebraic equation for η, as shown in the next section. The surface reconstruction procedure described above is mathematically valid for all Stokes waves, as proven by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. Although the procedure is likely also valid for all regular waves, a rigourous mathematical proof is lacking. Note that, by 'Stokes waves', [START_REF] Constantin | The trajectories of particles in Stokes waves[END_REF], [START_REF] Constantin | Pressure beneath a Stokes wave[END_REF] and [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF] mean all regular periodic symmetric waves with one crest and one trough per wavelength. According to this definition, solitary waves are a limiting case of Stokes waves. This is not the case with the more conventional definition, where 'Stokes waves' are approximations derived from a Poincaré-Lindstedt perturbation scheme that breaks down in shallow water and thus cannot describe solitary waves.

New equations for the surface recovery

Let Q be yet another holomorphic function such that

Q(z) ≡ z z 0 P(z ) -gd dz = z z 0 1 2 B -w(z ) 2 dz , (4.1)
where z 0 is an arbitrary constant. Taking z 0 at the origin of the free surface (i.e. x = 0, y = a, thus z 0 = ia) and choosing the free surface as integration path, the definition (4.1) yields

Q s (x) = x 0 P(x + iη(x )) -gd [1 + iη x (x )] dx , (4.2)
thus, substituting (3.3) and after some elementary algebra,

Q s (x) = x 0 gη(x ) dx + i [η(x) -a] B - 1 2 ga - 1 2 gη(x) . (4.3)
It is worth noting that the choice of z 0 in the definition of Q is obviously unessential, with z 0 = ia being chosen above for convenience. In particular, if z 0 is any point at the free surface (i.e. z 0 = x 0 + iη 0 with η 0 = η(x 0 )), then one obtains a relation similar to (4.3) with

x 0 replaced by

x x 0
and a replaced by η 0 . The imaginary part of (4.3) provides a quadratic algebraic equation in η, relating the latter to Im{Q s }. Solving this quadratic equation for η, one immediately obtains

η = g -1 B -(B -ga) 2 -2gIm{Q s } , (4.4)
where the minus sign in front of the radical is chosen (instead of the plus sign) in order to fulfil the relation η(0) = a. In (4.4), substituting Im{Q s } obtained from (4.2), we have a local algebraic (i.e. neither differential nor integral) relation defining the free surface η because, in practice, the integral in (4.2) can be computed analytically in closed form (see the examples in § 5).

The relation (4.4) is an exact solution of (3.8) written in implicit form, as one can easily check directly from the derivative of (4.4) and using (4.2), i.e.

η x = [Im{Q s }] x (B -ga) 2 -2gIm{Q s } = Im{P s } + (Re{P s } -gd) η x B -gη , (4.5)
then solving (4.5) for η x , yielding (3.8) precisely.

For practical applications, the relation (4.4) is obviously a major improvement compared with (3.8). Indeed, using (3.8) requires the numerical resolution of a nonlinear differential equation of the form η x = Φ(x, η), while with (4.4) we have to solve an equation of the form η = Ψ (x, η). Thus, the latter gives access to the free surface at a single point x = x 0 without computing η at the neighbouring points, while using (3.8) requires the determination of η over the whole period.

Note that analytic (explicit or implicit) solutions are not necessarily useful in practice if they involve hard-to-compute special functions or/and slowly convergent series. This is not the case here since, for practical surface recovery, few elementary functions are involved in (4.4), as shown below.

An equation for the Bernoulli constant is derived averaging the imaginary part of (4.3) and using (2.1)

2aB = g a 2 -η 2 -2 Im{Q s } . (4.6)
This equation for B is simpler than that derived by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. However, the relation (4.6) involves integrations and is therefore non-local. An alternative explicit closed-form expression for B is obtained if both a and b are already known from, respectively, the local relations (3.6) and (3.7). Indeed, from (4.4) written at the trough (x = π/k, η = -b), one obtains

B = 1 2 g(a -b) -(a + b) -1 Im{Q(π/k -ib)}, (4.7)
that is a very simple expression for B once a, b and Q are known.

With the pressure at the bottom being known, it is trivial to obtain P and subsequently Q, as indicated in the next section. Once a, b and Q have been determined, an explicit expression for η can be obtained numerically via fixed-point iterations of (4.4), whose convergence is studied in § 6.

Examples of analytic approximations for the function Q

For the recovery of the surface wave profile η, given the pressure p at the seabed, one can proceed as in the two examples below.

Fourier reconstruction

For periodic waves, the pressure measurements at the bed can be approximated by a N th-order Fourier polynomial and the function P is obtained immediately:

p(x) ≈ gd + N |n|>0 p n exp(inkx) ⇒ P(z) = gd + N |n|>0 p n exp(inkz) exp(nkd) , (5.1)
with p -n = p n since p is real. The amplitude a and the trough height b are then determined from (3.6) and (3.7). The function Q is subsequently easily obtained as

Q(z) = N |n|>0 p n ink exp(inkz) -exp(-nka) exp(nkd) , (5.2)
showing that the determination of Q from p is a trivial matter. Note that the antiderivative, required to determine Q, has a natural smoothing effect (i.e. the Fourier coefficient of Q decay faster than those of P) that is an interesting feature for practical applications. The Bernoulli constant B is then computed from (4.7) and, finally, one obtains the free surface for -b < η < a iterating (4.4).

Cnoidal-like reconstruction

The Fourier polynomial approximation (5.1) is not the only possible approximation and it is not suitable in shallow water, especially for solitary waves. [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF] suggested another simple approximation able to deal with waves in shallow water. It is derived from the shallow-water approximation and involves the Jacobian elliptic functions {sn, cn, dn} and the elliptic integrals K and E of parameter m [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], i.e.

p(x) ≈ gd

+ A dn 2 (κx | m) -E(m)/K(m) , (5.3)
where the (real) parameters d, A, κ and m have to be determined from the experimental data, for instance via a least-squares minimization. Note that the approximation (5.3) assumes a wave crest at abscissa x = 0; if this is not the case, a phase shift parameter has to be introduced and determined from the data. Note also that, with (5.3), the condition (2.4) is fulfilled identically because the squared bracket averages to zero exactly (for all m and κ).

The holomorphic function P is then obtained replacing x by z + id in (5.3), i.e.

P(z) = gd + A dn 2 (κ(z + id) | m) -E(m)/K(m) , (5.4)
and the function Q is subsequently derived at once substituting (5.4) into (4.1):

Q(z) = (A/κ) [Z(κ(z + id) | m) -Z(iκ(a + d) | m)] , (5.5)
where Z is the Jacobi zeta function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. For (2π/k)periodic waves, we have κ = kK(m)/π, and for solitary waves (i.e. k → 0, m → 1) the approximation (5.5) becomes

Q(z) = (A/κ) [tanh(κ(z + id)) -i tan(κ(a + d))] .
(5.6)

This example also shows that the calculation of Q is a simple matter.

Remarks

Other basis of functions can be used to represent the pressure at the bottom p, but the two examples above are a priori the more relevant for practical applications. The convergence of (5.1) as N → ∞ has already been addressed by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. Their numerical example also shows that this Fourier approximation can be used in numerical applications. Further numerical examples are provided in § 7.

The validity of approximations (5.4)-(5.5) is not limited to shallow water. Indeed, this kind of approximation is quite accurate for all depths, even for waves of large amplitude [START_REF] Clamond | Steady finite amplitude waves on a horizontal seabed of arbitrary depth[END_REF][START_REF] Clamond | Cnoidal-type surface waves in deep water[END_REF]. Further evidence is given in § 7 below. If more accuracy is required, then higher-order cnoidal-like approximations can be introduced [START_REF] Fenton | A high-order cnoidal wave theory[END_REF][START_REF] Clamond | Accurate simple approximation for the solitary wave[END_REF].

Convergence outline

A remaining open question is whether or not the iterations of (3.7) and (4.4) converge. Numerical investigations suggest that converge occurs (see § 7), this convergence is supported by the theoretical investigation of this section.

Following [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF], the convergence of (3.7) is investigated introducing the functional iterations

y = F(y) ≡ Re{P(π/k + iy)} g -d = y + p(π/k, y) g , y ∈ [-d, -b]. (6.1)
Iterations of (6.1) converge if F is a contraction, i.e. convergence is ensured if

-1 < F y < 1, divergence occurs if |F y | > 1 and if |F y | = 1 iterations may converge
or not depending on F [START_REF] Devaney | An Introduction to Chaotic Dynamical Systems[END_REF]. For (6.1), convergence is thus guaranteed if

-2g < p y < 0 for y ∈ [-d, -b].
With the results by [START_REF] Constantin | The trajectories of particles in Stokes waves[END_REF] and [START_REF] Constantin | Pressure beneath a Stokes wave[END_REF] yielding a rigourous proof of the inequality p y -g under the troughs of Stokes waves, the upper bound for contractivity is fulfilled. The lower bound p y > -2g is obviously satisfied for solitary waves (since then p y → -g as x → ±∞), but a rigourous proof for general waves is apparently not available in the literature. Numerical experiments (see figure 1 below) suggest that iterations (3.7) converge, even faster than iterations (3.6).

In order to investigate the convergence of (4.4), let the functional iterations be

y = G(y) ≡ B g - (B -ga) 2 -2gIm{Q(x + iy)} g , y ∈ [-d, η(x)], (6.2)
for a known holomorphic function Q at a fixed abscissa x and for given a and B.

Iterations of (6.2) converge if G is a contraction (i.e. if -1 < G y < 1), the definition of G giving

G y = [Im{Q}] y (B -ga) 2 -2gIm{Q} = [Im{Q}] y B -gy . (6.3)
From the definitions (3.1) and (4.1) of, respectively, the functions P and Q, using the Cauchy-Riemann relations and the Bernoulli equation (2.2), we have (6.4) so the relation (6.3) can be rewritten

[Im{Q}] y = Re{P} -gd = 1 2 B -1 2 u 2 + 1 2 v 2 = p + gy + v 2 ,
G y = p + gy + v 2 B -gy = p + gy + v 2 2p + gy + u 2 + v 2 .
(6.5)

The upper bound for contraction condition G y < 1 yields the inequality p + u 2 > 0 that is always satisfied, except at the (angular) crest of the highest wave where p = u = v = 0 (stagnation point at the crest). Conversely, the lower bound G y > -1 yields the inequality p + v 2 + B > 0 that is obviously always satisfied. Therefore, G is a contraction everywhere in the bulk and at the free surface and the iterations (6.2) converge to a unique solution [START_REF] Devaney | An Introduction to Chaotic Dynamical Systems[END_REF], except perhaps at stagnation points of the free surface (angular crests forming a 120 • inner angle).

Note that this analysis is valid for all regular waves, periodic or not, with one or more crests and troughs per period, symmetric or not. In other words, this proof is not limited to Stokes waves, unlike the proof for (3.6) and the incomplete proof for (3.7) (although convergences of these iterations likely also occur for more general solutions). Note also that, with G y being independent of a (i.e. independent of the choice of z 0 in the definition of Q), the convergence proof for the iterations (6.2) is valid wherever z 0 is chosen in the fluid.

Numerical example

We illustrate the new reconstruction procedure with the same numerical example used by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. Thus, we consider an exact surface wave with a wave height/depth ratio equal to 0.25 and wavelength/depth ratio equal to 5 (i.e. a/d ≈ 0.25, ka ≈ 0.314), which is a wave of significant amplitude in rather deep water. This exact solution is computed numerically solving the Euler equations [START_REF] Fenton | The numerical solution of steady water wave problems[END_REF] and the exact pressure at the bed is thus obtained. We take 32 samples of the bottom pressure equally spaced over one wavelength. These are the pressure 'measurements' we use for the reconstruction. We suppose that ρ and g are known parameters, the other parameters having to be computed from the data.

Similarly to the approach of [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF], we do not make any attempt to introduce sophisticated numerical techniques. We use basic numerical implementation of the analytical formulae given above.

Fourier recovery

An analytic approximation of the pressure at the bottom is obtained by fitting, via least-squares minimization, a fifth-order Fourier polynomial among the pressure data at the bed. The accuracy of this reconstruction is illustrated by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF] with various data and graphics. Thus, we focus here on the iterations (3.7) required to compute the trough height b.

Starting with initial guess b (0) = -a, so that it is guaranteed to lie in the fluid domain, the fixed-point iterations of (3.7) converge, even faster than those of (3.6), as illustrated in figure 1. This result provides numerical evidence of convergence, but a rigourous mathematical proof is still lacking. Note that the rate of convergence is linear here because we use simple fixed-point iterations, but more sophisticated algorithms (e.g. Newton, Halley) will easily provide faster rates of convergence.

For the reconstruction of the whole surface, the formula (4.4) provides the same result (within tiny numerical errors) as the more involved resolution of (3.8) used by [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF]. This is to be expected because (4.4) is an exact solution of (3.8). The speed of the two methods is very different, however.

Using M = 1024 equally spaced nodes and fast Fourier transforms, the differential equation (3.8) is solved iteratively via a pseudo-spectral method in O(M log M) operations. This is the fastest possible accurate algorithm. Another resolution of comparable accuracy using, say, a Runge-Kutta method, would be much slower (typically in O(M 2 ) operations). Conversely, the fixed-point iterations of (4.4) give the same result in O(M) operations. For the example considered here, the surface is obtained about 13 times faster using (4.4) than using (3.8). This gain is consistent with the value log 2 (1024) = 10 (characterizing the different number of operations) and the fact that the Bernoulli constant is implicit using (3.8), thus requiring nested loops to compute the solution [START_REF] Clamond | Recovery of steady periodic wave profiles from pressure measurements at the bed[END_REF], while it is explicitly known in the new algorithm.

We also investigated the effect of the number of harmonics N in the reconstruction of the free surface. Figure 2 shows that the accuracy increases with the number of harmonics, and that the trough height b is always a fixed point of (3.7).

Cnoidal-like recovery

We analyse here the relevance of the cnoidal-like ansatz (5.3). It should be emphasized that using this ansatz is a priori difficult because we are dealing with a wave in rather deep water.

The parameters involved in (5.3) being obtained by least-squares minimization, the pressure at the bed (figure 3a) is well-approximated by the ansatz (5.3), i.e. the accuracy is about four digits (figure 3b) and the Fourier analysis shows a spectral content matching the exact one (figure 3c). The reconstructed free surface is reasonably accurate (figure 3d,e), although not as accurate as the fifth-order Fourier reconstruction.

However, the fifth-order Fourier approximation has more free parameters to be adjusted than the cnoidal-like approximation. For a fair contest, we should compare two ansatz having exactly the same number of free parameters. Therefore, in that respect, we should compare (5.3) with a second-order Fourier approximation, that is (5.1) with N = 2. We thus find that these two approximations have comparable accuracies (figure 4).

Conclusion

We have derived some new exact relations for the reconstruction of the free surface wave profile from pressure measurements at the bed. These new relations allow an easy and fast reconstruction, which is interesting for practical applications. This was possible thanks to the introduction of the holomorphic function Q. Indeed, the function Q yields a reformulation of the problem into a much simplified, but still exact, form. Efficient algorithms and easy proofs then follow at once. It remains to determine whether the present reconstruction procedure works for limiting waves with a stagnation point at the angular crest. This is a problem of more theoretical than practical interest. For practical applications, the method efficiency is illustrated with some numerical examples.

The recovery procedure described here is promising. It would be interesting to see whether a similar approach could be achieved to treat other water wave problems, such as unsteady flows and varying bottom. These matters will be the focus of future investigations.

  FIGURE 1. (Colour online) Numerical fixed-point iterations of (3.6) (upward triangles) and(3.7) (downward triangles) for the fifth-order Fourier reconstruction.

  FIGURE 2. (Colour online) Error of the recovered amplitude a (upward triangles) and trough height b (downward triangles) for various-order Fourier reconstructions.

FIGURE 3 .

 3 FIGURE 3. (Colour online) Recovery using the cnoidal-like ansatz (5.3): (a) profile of pressure at the bed; (b) difference between fitted and exact pressure at the bed; (c) Fourier spectrum of pressure at the bed; (d) profile of the recovered free surface; (e) recovered free surface error.

  FIGURE 4. (Colour online) Free surface difference between cnoidal-like and second-order Fourier recoveries.


Xudong Qin

Simon Bliudze

Eric Madelaine

Min Zhang

Using SMT engine to generate Symbolic Automata

Open pNets are used to model the behaviour of open systems, both synchronous or asynchronous, expressed in various calculi or languages. They are endowed with a symbolic operational semantics in terms of so-called "Open Automata". This allows us to check properties of such systems in a compositional manner. We implement an algorithm computing these semantics, building predicates expressing the synchronization conditions between the events of the pNet sub-systems. Checking such predicates requires symbolic reasoning over first order logics, but also over application-specific data. We use the Z3 SMT engine to check satisfiability of the predicates, and prune the open automaton of its unsatisfiable transitions. As an industrial oriented use-case, we use so-called "architectures" for BIP systems, that have been used in the framework of an ESA project and to specify the control software of a nanosatellite at the EPFL Space Engineering Center. We use pNets to encode a BIP architecture extended with explicit data, and compute its open automaton semantics. This automaton may be used to prove behavioural properties; we give 2 examples, a safety and a liveness property.

Introduction

In the nineties, several works extended the basic behavioural models based on labelled transition systems to address value-passing or parameterized systems, using various symbolic encodings of the transitions [START_REF] De Simone | Higher-level synchronising devices in MEIJE-SCCS[END_REF][START_REF] Larsen | A context dependent equivalence between processes[END_REF][START_REF] Hennessy | Symbolic bisimulations[END_REF][START_REF] Lin | Symbolic transition graph with assignment[END_REF]. In [START_REF] Lin | Symbolic transition graph with assignment[END_REF], H.M. Lin addressed value-passing calculi, for which he developed a symbolic behavioural semantics, and proved algebraic properties. Separately J. Rathke [START_REF] Hennessy | Bisimulations for a calculus of broadcasting systems[END_REF] defined another symbolic semantics for a parameterized broadcast calculus, together with strong and weak bisimulation equivalences, and developed a symbolic modelchecker based on a tableau method for these processes. Thirty years later, no practical verification approach and no verification platform are using this kind of approaches to provide proof methods for value-passing processes or open process expressions.

Parameterized Networks of Synchronized Automata (pNets) were proposed to give a behavioural specification formalism for distributed systems, synchronous, asynchronous, or heterogeneous. They are used in VerCors [START_REF] Henrio | Verifying the correct composition of distributed components: Formalisation and Tool[END_REF], a platform for designing and verifying distributed systems, as the intermediate language for various high-level languages. The high-level languages in VerCors formalize each component of the distributed system and their composition. pNets provides the core low-level semantic formalism for VerCors, and is made of a hierarchical composition of (value-passing) automata, called parameterized labelled transition systems (pLTS), where each hierarchical level defines the possible synchronization of the lower levels. Traditionally, pNets have been used to formalize fully defined systems or softwares. But we want also to define and reason about incompletely defined systems, like program skeletons, operators, or open expressions of process calculi. The open pNet model addresses this problem, using "holes" as process parameters, representing unspecified subsystems. The pNet model was developed in a series of papers [START_REF] Henrio | pNets: an Expressive Model for Parameterised Networks of Processes[END_REF][START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF] in which many examples have been introduced showing its ability to encode the operators from some other algebras or program skeletons. The operational semantics of an (open) pNet is defined as an Open Automaton in which Open Transitions contain logical predicates expressing the relations between the behaviour of the holes, and the global behaviour of the system. In the previous publication, only a sketch of a procedure allowing to compute these semantics was presented, together with a proof of finiteness of the open automaton, under reasonable hypotheses on the pNet structure.

Implementing these semantics raised several challenges, in order: to get a tool that could be applied to pNets representing various languages, in particular various actions algebras, with their specific decision theories, to separate clearly the algorithm generating the transitions of the open automaton from combination of all possible (symbolic) behaviours, from the symbolic reasoning part, specifically here using an SMT engine to check the satisfiability of the predicates generated by our algorithm, to build a prototype and validate the approach on our basic case-studies, and understand the efficiency of the interaction with the SMT solver.

In the long term, we want to be able to check the equivalence between open systems encoded as pNets. The equivalence between pNets is "FH-bisimulation" [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF], a dedicated version of symbolic bisimulation taking the predicate of the open transitions into account when matching such open transitions. We foresee that the interplay with the SMT solver that we use here for satisfiability of open transitions will be similar with what we need when proving (symbolic) equivalence between open transitions.

Contribution In the article we show how:

-We define the open automaton generation algorithm. We implemented a full working prototype, within the VerCors platform. In the process, we improved the semantics rules from [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF], and add features in the algorithm to deal with the full model, including management of variables and assignments. -We implement the interaction between our algorithm and the Z3 SMT solver, for checking satisfiability of the transitions generated by the algorithm.

-We show the interest of this approach on an industry-inspired case-study, namely one architectural pattern extracted (and extended) from the BIP specification of a nanosatellite on-board software.

Related work. Very few attempts were made to develop symbolic bisimulation approaches for the value-passing process algebras and languages-our long-term goals-especially, there is no algorithmic treatment of the symbolic systems developed by interacting with automatic theorem provers. The closest work is the one already mentioned from J. Rathke [START_REF] Hennessy | Bisimulations for a calculus of broadcasting systems[END_REF], who developed the symbolic bisimulation for a calculus of broadcasting system (CBS). CBS is similar with classic process calculi such as CCS and CSP, but communicating by broadcasting values, transmitting values without blocking. That makes the definition of the symbolic semantic and bisimulation equivalence different from the classic works.

For other applications, such as the analysis of programming languages, there exist dedicated platforms using external automatic theorem provers (ATP) or automatic tactics from interactive theorem provers (ITP), to perform symbolic reasoning, and for example to discharge some subgoals in the proofs. Tools like Rodin [START_REF] Déharbe | Integration of smt-solvers in b and event-b development environments[END_REF][START_REF] Déharbe | Integrating smt solvers in rodin[END_REF] have already integrated several provers, like Z3, as modules for proving the proof obligations generated from a user model. The prover we use, which also happens to be Z3, is developed by Microsoft Research based on the satisfiability modulo theories framework (SMT), is mainly applied in extended static checking, test case generation, and predicate abstraction. In a similar way, there are several ATPs/ITPs we could consider to use for result pruning and bisimulation checking in our algorithm, as an alternative to Z3, such as CVC4 [START_REF] Barrett | Cvc4. In: Computer aided verification[END_REF], Coq [START_REF] Armand | A modular integration of sat/smt solvers to coq through proof witnesses[END_REF] or others.

BIP (Behaviour-Interaction-Priority) [START_REF] Basu | Rigorous component-based system design using the BIP framework[END_REF] is a framework for the componentbased design of concurrent software and systems. In particular, the BIP tool-set comprises compilers for generating C/C++ code, executable by linking with one of the dedicated engines, which implement the BIP operational semantics [START_REF] Baranov | Offer semantics: Achieving compositionality, flattening and full expressiveness for the glue operators in BIP[END_REF]. This approach ensures that any property, shown to hold on a given BIP model, will also hold by construction on the generated code. BIP Architectures [START_REF] Attie | A general framework for architecture composability[END_REF] formalise design patterns, which enforce global properties characterising the coordination among the components of the system. They provide a compositional approach, ensuring correctness by construction during the design of BIP models. In [START_REF] Attie | A general framework for architecture composability[END_REF], it was shown that application of architectures is compositional w.r.t. safety properties, i.e. when several architectures are applied, each enforcing a safety property, the resulting system satisfies their conjunction.

But the interaction feature in architectures does not handle data-sensitive interaction constraints. Using an encoding of architectures, extended with datadependant interactions, into open pNets was an interesting alternative to a direct extension to the architecture semantics.

Structure. In section 2 we give a description and a formal definition of the pNet model, as found in previous publications. Then in section 3 we present our usecase, based on a BIP architecture from the nano-satellite case-study. Section 4 recalls briefly the operational semantics of pNet. Section 5 explains in details the algorithm used to compute this semantics, including the interaction with Z3, and shows the full result of the semantic computation on the running example. Finally we conclude and discuss perspectives in Section 6.

Background: pNets definition

This section introduces pNets and the notations we will use in this paper. Then it gives the formal definition of pNet structures, together with an operational semantics for open pNets.

pNets are tree-like structures, where the leaves are either parameterized labelled transition systems (pLTSs), expressing the behaviour of basic processes, or holes, used as placeholders for unknown processes, of which we only specify their set of possible actions, named sort. Nodes of the tree (pNet nodes) are synchronizing artifacts, using a set of synchronization vectors that express the possible synchronization between the parameterized actions of a subset of the sub-trees.

Notations. We extensively use indexed structures over some countable indexed sets, which are equivalent to mappings over the countable set. a i∈I i denotes a family of elements a i indexed over the set I. When this is not ambiguous, we shall use notations for sets, and typically write "indexed set over I" when formally we should speak of multisets, and write x ∈ a i∈I i to mean ∃i ∈ I. x = a i . An empty family is denoted ∅. We denote classically a a family when the indexing set is irrelevant. is the disjoint union on indexed sets.

Term algebra. Our models rely on a notion of parameterized actions that are symbolic expressions using data types and variables. As we want to encode the low-level behaviour of possibly very different programming languages, we do not want to impose one specific algebra for denoting actions, nor any specific communication mechanism. So we leave unspecified the constructors of the algebra that will allow building expressions and actions. Moreover, we use a generic action interaction mechanism, based on unification between two or more action expressions. This will be used in the semantics of synchronization vectors to express various kinds of communication or synchronization mechanisms.

Formally, we assume the existence of a term algebra T Σ,V , where Σ is the signature of the data and action constructors, and V a set of variables. Within T Σ,V , we distinguish a set of data expressions E V , including a set of Boolean expressions

B V (B V ⊆ E V ). On top of E V we build the action algebra A V , with A V ⊆ T V , E V ∩ A V = ∅;
naturally action terms will use data expressions as sub-terms. The function vars(t) identifies the set of variables in a term t ∈ T .

pNets can encode naturally the notion of input actions as found, e.g. in value-passing CCS [START_REF] Milner | Communication and Concurrency[END_REF] or of usual point-to-point message passing calculi, but it also allows for more general mechanisms, like gate negotiation in Lotos [START_REF]ISO: Information Processing Systems -Open Systems Interconnection -LOTOS -A Formal Description Technique based on the Temporal Ordering of Observational Behaviour[END_REF], or broadcast communications. Algebra presentations. In practice, the parameterization of the pNet model by some specific action algebra is realized by the definition of a many-sorted "algebra presentation". It will be used to check the well-formedness of a pNet system, and to define the translation of the pNet semantics into the SMT engine input language ( [START_REF] Barrett | The SMT-LIB Standard: Version 2[END_REF]).

Definition 1. An algebra presentation is a triple P = Sorts, Constrs, Ops , where -Sorts is a set of data sorts -Constrs is a set of constructor operators: for each Con ∈ Constrs, arity(Con) = n ∈ N is its arity and 'Con : (sel 1 , sort 1 ), . . . , (sel n , sort n ) → sort' is its signature with the associated selectors. For each argument, the pair (sel i , sort i ) defines an auxiliary operator of name sel i with signature sel i : sort → sort i . -Ops is a set of auxiliary operators, with their arity and signature, of the form: Op : sort 1 , . . . , sort n → sort -Constrs(sortname) and Sels(sortname) are, respectively, the sets of constructors and selectors of the sort sortname Constructors of arity 0 are called constants, and denoted Consts(P).

Sorts Bool and Int are predefined with standard operators. Sort Action also, with a constructor Synchro denoting a synchronized action, i.e. an "internal" action that cannot be further synchronized with the environment. It also comes with an overloaded FUN constructor, used to build actions with arguments, that will be instantiated to the required sorts for a given pNet.

The definition of an Algebra Presentation, and a set of variables V fixes the Term algebra elements T Σ,V , B V , A V .

The (open) pNets Core Model

A pLTS is a labelled transition system with variables, which can be manipulated, defined, or accessed inside states, actions, guards, and assignments. Each state has its set of variables called state variables, which can only be modified by the assignment in transitions targeting this state. A global state variable of a pLTS is a state variable defined in all states. Note that we make no assumptions on finiteness of the set of states, nor on finite branching of the transition relation.

We first define the set of actions a pLTS can use. Let a range over action labels, op are operators, and x i range over variable names. Action terms are: Now, we define pNet nodes as constructors for hierarchical structures. A pNet node has a set of sub-pNets that can be either pNets or pLTSs, and a set of holes, playing the role of process parameters (i.e. unknown in the environment).

A composite pNet consists of a set of sub-pNets, each exposing a set of actions. The relation between actions of a pNet and those of its sub-pNets are given by synchronization vectors, which synchronize one or several internal actions, and expose a single resulting global action.

Definition 3 (pNets).

A pNet is a hierarchical structure where leaves are pLTSs and holes: pNet pLTS | pNet i∈I i , J, SV k∈K k , with I, J, K potentialy infinite, where

• pNet i∈I i is the family of sub-pNets; • J is a set of indexes, called holes. I and J are disjoint:

I ∩J = ∅, I ∪J = ∅ • SV k∈K k is a set of synchronisation vectors (K ∈ I V ). ∀k ∈ K, SV k = α l∈I k J k l → α k | g k , where α k ∈ A V , I k ⊆ I, J k ⊆ J, and vars(α k ) ⊆ l∈I k J k vars(α l ).
The global action of a vector SV k is α k . The Boolean expression g k , such that vars(g k ) ⊆ l∈I k J k vars(α l ), is a guard associated to the vector.

In Fig. 2, we show examples of these constructs, with two pLTSs, one hole and one pNet node encoding our running example.

Running example

As a running example we use the Failure Monitor architecture from the Cu-bETH nanosatellite on-board software case-study [START_REF] Mavridou | Architecture-based design: A satellite on-board software case study[END_REF] realised using BIP. The architecture-based design process in BIP takes as input a set of components providing basic functionality of the system and a set of temporal properties that must be enforced in the final system. For each property, a corresponding architecture is identified and applied to the model, thereby potentially introducing additional coordinator components and modifying the connectors that define synchronisation patterns among ports of components. Figure 1 shows a refined version of the Failure Monitor architecture used in [START_REF] Mavridou | Architecture-based design: A satellite on-board software case study[END_REF]. Contrary to standard BIP models, architectures comprise one or several operand components, whereof only the set ports is given. Here, the operand component is B and its interface consists of the ports finish, resume and fail. The two coordinator components-Control and Timer -are standard BIP components insofar as they also have their behaviour specified by finite automata extended with local data variables. Transitions of these automata are labelled with the ports of the corresponding components, Boolean guards and update functions on local variables. For instance the loop transition t 1 tick,[t>0],t:=t-1 ----------→ t 1 in the Timer component is labeled by the port tick, it can be fired only when the current value of the local variable t is greater than 0. Upon firing, this transition decrements the value of t by 1. When omitted, the default guard (resp. update function) is the constant predicate true (resp. the skip operator). The constant Max, in t 0 start,t:=Max --------→ t 1 , is a parameter of the architecture. Connectors are hierarchical, tree-like structures with component ports at the leaves. They define sets of interactions, based on the attributes of the connected ports [START_REF] Bliudze | The algebra of connectors-Structuring interaction in BIP[END_REF], which may be either trigger (triangles in Fig. 1) or synchron (bullets in Fig. 1). If all sub-connectors of a connector are synchrons, then an interaction may be executed by the connector only if each subconnector can contribute. If at least one of the sub-connectors is a trigger, then any interaction consisting of contributions of any set of sub-connectors, involving at least one of the triggers, can be executed. For instance, the two ports Timer .start and Control .fail are always synchronised, since they belong to the same binary sub-connector, where they are both synchrons. In particular, this means that whenever the transition s 0 fail --→ s 1 is fired, so is the transition t 0 start,t:=Max --------→ t 1 , initialising the timer. The binary connector Timer .start • --•Control .fail is a sub-connector of a hierarchical connector, where the port B.fail is a trigger. Thus, the above interaction can only happen together with B.fail, forming a ternary interaction. On the contrary, being a trigger, the port B.fail can fire alone, forming a singleton interaction. The composition semantics of BIP systems consists in firing exactly one interaction, enabled through at least one of the top-level connectors, at each execution round.

fail t := t-1 [t > 0] t 0 t 1 s 2 s 1 s 0 fail t im e o u t
Finally, priorities-defined by a strict partial order on the set of possible interactions-narrow the choice among the enabled interactions at any given round. The default priority is the so-called maximal progress, whereby among any two interactions a ⊂ b (as sets of ports), b has higher priority than a. Application of the Failure Monitor architecture ensures that, whenever a failure is registered in the operand component, the system will be reset, unless a resumption is registered within Max time units (more details in Sect. 5.4).

SV5 <-, -, f inish, -> → f inish SV4 <reset, -, -> → reset SV3 <timeout, timeout, -> ask SV2 <-, tick, -> → tick SV1 <resume(b1), resume(b2), resume(b0)> → resume [b1 = b2 ∧ b1 ∨ b2 =⇒ b0] SV0 <f ail(b1), start(b2), f ail(b0)> → f ail [b1 = b2 ∧ b1 ∨ b2 =⇒ b0] tick[t > 0] {t M ax} {t := t -1}
Figure 2 shows a pNet encoding of the above Failure Monitor architecture. This encoding is structural: each coordinator component is encoded as a pLTS, the operand component-as a hole; connectors of the BIP model are encoded as synchronisation vectors. Each connector that does not involve triggers is trivially encoded by a synchronisation vector comprising the same ports. In order to encode the semantics of the connectors involving triggers, we 1) in the pLTS encoding the coordinator components, add loop transitions to ensure that all ports involved in such connectors are enabled in all states, 2) associate a Boolean value to each of these ports: the original transitions carry the value true (e.g.

s 0 fail(true)
-----→ s 1 ), the added loops carry the value false (e.g. s 2 fail(false) ------→ s 2 ), 3) add to the corresponding synchronisation vectors the Boolean predicate encoding the connector structure. For example, SV 0 encodes the connector discussed above: the predicate (b1 = b2) ∧ (b1 ∨ b2 ⇒ b0) means that the "true" transitions Control .fail and Timer .start can only fire together (b1 = b2) and whenever one of them fires, B.fail must fire also (b1 ∨ b2 ⇒ b0). This encoding can be systematically obtained for any hierarchical BIP connector [START_REF] Bliudze | Causal semantics for the algebra of connectors[END_REF]. Although, for the sake of brevity, we omit priorities from the encoding, this can be easily achieved, by introducing additional Boolean variables for relevant ports [START_REF] Baranov | Offer semantics: Achieving compositionality, flattening and full expressiveness for the glue operators in BIP[END_REF].

Operational Semantics for Open pNets

The semantics of open pNets will be defined as an open automaton, that is an automaton where each transition composes transitions of several LTSs with the actions of some holes; the transition occurs if some predicates hold, and can involve a set of state modifications. Each state of an open automaton has a set of state variables that can be assigned by incoming transitions. Strictly speaking, the LTSs at the leaves of the open automaton are a restricted form of pLTSs, where labels are parametrised actions, but include no guard nor assignments.

Definition 4 (Open transitions

). An open transition OT over a set S i , s 0i , → i i∈I of LTSs, a set J of holes, and a set of states S is a structure of the form:

•••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• {s i ai -→ i s i } i∈I , { bj -→ j } j∈J , Pred , Post s α -→ s ,
where s, s ∈ S, the a i , b j , α are action expressions, s i ai -→ i s i is a transition of the LTS S i , s 0i , → i , b j is an action of the hole j, and α is the resulting action of OT . Pred is a predicate over the variables of the terms, labels, and states s i , b j , s, α. Post is a set of equations that hold after the open transition, represented as a substitution {x k ← e k } k∈K where x k are variables of s and s i , whereas e k are expressions over the other variables of the open transition. States of open pNets. A state of an open pNet is a tuple of the states of its leaves (in which we denote tuples in structured states as . . . ). For any pNet p, let Leaves = S i , s i0 , → i i∈L be the set of pLTS at its leaves, then

States(p) = { s i∈L i | ∀i ∈ L.s i ∈ S i }.
A pLTS being its own single leave: States( S, s 0 , → ) = { s | s ∈ S}. The initial state is defined as: InitState(p) = s i0 i∈L .

Predicates. Let pNet, J, SV k∈K k be a pNet. Consider a synchronization vector SV k , for k ∈ K. We build a predicate MkPred relating the actions of the involved sub-pNets and the resulting actions. This predicate verifies:

MkPred(SV k , a i∈I i , b j∈J j , v) ⇐⇒ ∃(a i ) i∈I , (b j ) j∈J , v . SV k = (a i ) i∈I , (b j ) j∈J → v ∧ ∀i ∈ I. a i = a i ∧ ∀j ∈ J. b j = b j ∧ v = v
Example 1 (An open transition). This transition is generated by application of the vector SV0, synchronizing the initial actions of pLTSs Control and Timer (see Fig. 2), with an action of the hole B equal to fail. The local variable t of the Timer is assigned to its initial value Max. A full output of the use-case is provided in [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF].

•••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• {s0 fail(true) ------→ s1, t0 start(true) -------→ t1}{ hb -→}, hb = fail (true) ∧ v = fail , {t := Max } s0, t0 v -→ s1, t1
Structural Semantic Rules: The semantics of pNet in term of open automata has been defined in [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF], in the form of 2 structural rules, one for pLTSs, one for pNet nodes. These rules are slightly improved, adding guards in the synchronisation vectors and syntax for universal quantifier in the guards (see [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF]).

In [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF] we also proved that a pNet with finitely many pLTSs, holes, and synchronisation vectors, and with finite pLTSs, has a finite semantics. Remark that all the elements of such pNets and open automata being symbolic, they can represent many classes of unbounded systems.

Generation of Open Automata

In this section we describe the algorithm implementing the pNet semantics, the interaction with the Z3 SMT solver, and we show the result on our example.

Algorithm 1 starts with an open pNet, and builds its set of open transitions. Its main loop is a classical residual algorithm: starting from the initial global state, it picks a state in an unexplored set, and computes all possible OTs, adding their target states in the unexplored set, until this set is empty. As the set of global states is finite (under the conditions above), this terminates.

Algorithm 1 Open Automaton Generation

Input: A pNet P (cannot be a hole) 1: Initialize sets U = {InitState(P )} and E = ∅, for unexplored and explored global states, respectively; L = ∅ for the resulting OTs; 2: while !isEmpty(U ) do 3:

Choose S in U ; remove S from U , add S to E; 4:

OTs = MakeTransitions(P, S); 5:

for each OT ∈ OTs do Check satisfiability of OT using the SMT solver; 6:

if SAT (OT ) then 7:

{Add OT to L; 8:

Let S be the target global state of OT 9:

if (S ∈ U ∪ E) then Add S into U ; } 10: end for 11: end while 12: return OA = (InitState(P ), L);

The inside loop (MakeTransitions method) applies recursively the semantic rules following the structure of the pNet. When applied to a pLTSs at the leaves, we simply take the pLTS transitions of the corresponding local state and use the semantic rule to build the OT 4 . When applied to a pNet node we use two methods, combining and matching, to generate the open transitions in a hierarchical manner, as shown in Alg. 2. This method directly manages the holes of the node, so MakeTransitions is never called on a hole.

At the root of the pNet, the predicate of each OT is translated into SMTlib assertions, and checked for satisfiability. The final open automaton includes all satisfiable OTs, and the set of reachable global states.

Algorithm 2 MakeTransitions() for a pNet node

Input: a pNet node P with subnets sn and holes hole; a global state S.

1: Initialize empty list l and set L for sub-transitions and transitions, respectively; 2: for each Subnet in sn do \\ Recursively apply the semantic rules on the subnets 3: Store MakeTransitions(Subnet, S) in l; 4: end for 5: comb = Combining(l); 6: for each sv ∈ SV and each comb ∈ comb do 7:

ot = Matching(sv, comb, hole); 8:

if (ot is defined) then Store ot in L; \\ if Matching() succeeds 9: end for 10: return L;

Combining. The combining method enumerates all the possible behaviours of the subnets as all the possible combinations of their open transitions. Assume that there is a collection of n subnets. We denote ot i the set of open transitions of the i-th subnet (obtained in line 3 of the algorithm); "-" means that the subnet is not involved. The combination comb, a set of n-tuples, is the cartesian product

: comb = ({-} ∪ ot 1 ) × ({-} ∪ ot 2 ) × • • • × ({-} ∪ ot n ).
Matching. The Matching method builds the OTs of a pNet node from those of its subnets. For each synchronisation vector and each possible combination of behaviours of the subnets, as generated by the Combining method, it builds the corresponding open transition. Here, we only detail the construction of the predicate. From a synchronization vector sv = (a i ) i∈I (b j ) j∈J → v ∈ SV and its guard G k ; a tuple of open transitions C = (ot i ) i∈ [1,n] ∈ comb, such that, for each i ∈ [1, n], either ot i = -, or the result action of ot i is a i ; the hole behaviours Hole = (b j ) j∈J ; and a fresh variable v, representing the result action of the OT under construction, we build the predicate: Filtering. While matching a vector with a combination tuple, Matching tries to filter out some incompatibilities; there may be several reasons why the matching would fail:

MkPred (sv, C, hole) = (∀i ∈ I, a i = a i ) ∧ (∀j ∈ J, b j = b j ) ∧ (v = v ) ∧ G k . ot = •••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• {s0 f ail(true) ------→ s0, t0 resume(f alse) ---------→ t0}, { hb -→}, f ail(true) = f ail(b1) ∧ resume(f alse) = start(b2) ∧ hb = f ail(b0) ∧ v = f ail ∧ b1 = b2 ∧ (b1 ∨ b2) ⇒ b0, {} s0, t0 v -→ s0, t0
if some subnet is marked as inactive in the vector, and the chosen combination has an active behaviour at this position, if some subnet action expression in the vector does not match (by patternmatching) the corresponding action expression in C, or if the whole set of active subnet actions in the vector cannot be matched (by unification) with the corresponding action expressions on the tuple C. Even when unification succeeds, it is still possible that the resulting predicate would be unsatisfiable, because of some incompatibility involving the guards. In our algorithm, we choose to apply only the simplest filter inside the Matching method (before applying the predicate and OT construction). Matching and unification will be checked later, together with the guards collected from synchronisation vectors and from pLTS transtions, using the satisfiability check in the SMT engine.

Management of state variable assignments

In a pLTS, there may be several incoming transitions that assign potentially different values to a state variable. To handle such cases, the algorithm manages, for each pLTS state, a list of expressions collected from the assignments of each state variable. For a global state in the open automaton, the set of state variables (which may be used in a transition) is the disjoint union of sets of state variables of the individual pLTS states constituting this global state.

Pruning the unsatisfiable results

Our matching/filtering strategy builds some transitions where the predicates express incompatible constraints. Even if having an unsatisfiable (symbolic) transition would not be incorrect, we choose to minimize the open automaton (i.e. its number of transitions and states), by checking the predicates for satisfiability. In Fig. 3, we show an unsatisfiable open transition from the open automaton of our running example. It shows the case where the failure controller performs a "fail" action, while the timer executes a "resume". The chosen synchronization vector (SV0 from Fig. 2) does not match with these actions, since it expects Timer .start. This mismatch is materialised by the predicate fragment "resume(false) = start(b 2 )". Checking satisfiability requires some symbolic computation on the action expressions and the predicates, which may depend on the specific theory of the action algebra datatypes. The "Modulo Theory" part of SMT solvers is important here, so that the solver can use specific properties of each action algebra.

Translation to SMTlib

We check the satisfiability of each open transitions using the SMT solver Z3. In this section, we describe the translation of the algebra presentation, of assignments, and of the predicates.

Our implementation submits satisfiability requests to Z3 using its JAVA API. Here, for readability, we show the Z3 code using its SMT-LIB input language. Note that in the previous sections, the OTs were displayed in a simplified, human readable form. The input and output of our tool, and also the generated SMTlib fragments, are slightly more difficult to read, in particular because of structured names for the fresh variables generated by the algorithm, allowing tracability of the result [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF].

Figure 4 shows the translation of the transition of Fig. 3 in the SMTlib syntax. It contains the declaration of the BIP action algebra sorts and constructors, then the declaration of variables, and finally the predicate to be checked, encoded as a set of assertions. Here we display also the diagnosis ("sat" or "unsat") generated by Z3; Production of the SMT-lib code To build the input submitted to Z3 for each OT, we translate the algebra presentation, the predicates and the variable assignments into Z3 (Java-API) calls.

Translation of action algebra presentation. In [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF], we define the translation of an algebra presentation into SMTlib declarations (declare-datatypes and declare-fun). We also formalize a condition of well-formedness of pNets ensures that the generated code is correct, and will not raise runtime errors in the SMT engine. Note that the declare-datatypes command comprises both the action constructors from table 1 and also the constant action names from our example. Additionally, we will include axiomatisation of any required functions and predicates of the presentation data-types.

ask [t = 0] ⊳s2, t0⊲ f inish resume [h B = resume(b0)] ⊳s0, t0⊲ ⊳s1, t1⊲ resume [h B = resume(true)] f inish f inish [h B = resume(b0)] resume f ail [h B = f ail(b0)] tick [t = 0] {t ← t -1} f ail [h B = f ail(true)] {t ← M ax} reset
Translation of open transitions. In [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF] we formally define all steps of the translation of each open transition, including:

collect all variables in the transition, and declare them (using declare-const) check well-formedness and correct typing of expressions, translate the predicate into a conjunction of assertions, if present in the source state, translate the state-variable assignments into a disjunctive assertion. This translation ensures that no runtime error will occur in the SMT engine.

Figure 4 shows the decomposition of the predicate into a set of asserts, each encoding an elementary equality, inequality, or a guard. The result (sat or unsat) of the final check-sat command in the translation () is then decoded. To improve the readability of this figure, we used the following conventions: we omit the transitions of the two pLTS, and the set of "working" holes; and we directly write the resulting action as first element of each OT, rather than including it as an equality inside the predicate.

Result for the running example

Failure Monitor enforces 1) the safety property "the system reset never happens, unless asked for by a timeout following a failure", formalised in CTL by

ϕ ∧ AG(reset → ϕ), where ϕ = A[¬reset W ask] ∧ A[¬ask W fail],
(W being the weak until operator) and 2) the liveness property "a reset will be fired when asked for by a timeout": AG(ask → AF reset). The satisfaction of the safety property could be established by applying symbolic model checking techniques. However, in this example, it is obvious by inspection of the open automaton. The satisfaction of the liveness property relies on the above observation that in the state s2, t0 only the reset transition involves Control . Therefore, under resonable scheduling assumptions, reset will always be fired.

Conclusion and Discussion

The formal definitions and properties of the open pNets model were published in [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF]. In this new work we describe an implementation of the model and its semantics construction, including its interaction with the Z3 SMT engine. The implementation has two parts: the first is a finitary algorithm that builds all possible combinations of synchronisations through the pNet hierarchical structure. The result is a so-called open-automaton, which transitions contains predicates relating the actions of the pNet holes and controllers. Some of the open transitions obtained at this step, may contain predicates which do not represent any possible concrete instantiations. In the second part of the tool we use the SMT solver Z3 for checking the satisfiability of the predicate in each open transition. To this end, we encode into Z3 the representations of the action algebra and the predicates before submitting them to the Z3 solver. In this paper, we used a running example, based on a BIP architecture from an earlier nanosatellite case study [START_REF] Mavridou | Architecture-based design: A satellite on-board software case study[END_REF]. This example shows that open-automata-based semantics can be instrumental in verifying the properties enforced by the architectures through an encoding into open pNets. This encoding-which we intend to formalise and prove correct in a separate paper-also opens the way for an extension of BIP architectures with the transfer of data among variables of different components. Indeed, such data transfer can be easily encoded using the predicates associated to synchronisation vectors in open pNets. The encoding of open transitions into SMTlib and the availability of theories can guide the definition of such an extension. Our case-studies show that our encoding successfully identifies the unsatisfiable open transitions and that the resulting automata correctly reflect the expected movements of the encoded process expressions.

Naturally, our next goals after the generation of the open automata will be to model-check logical properties, and to check equivalence of pNets. While model-checking open automata seems easy to define, equivalence checking is more challenging. In [START_REF] Henrio | A Theory for the Composition of Concurrent Processes[END_REF], we have already found the FH-bisimulation, to be a suitable definition. But weak equivalences, or refinements, will definitely be useful when comparing different pNets with different structure. For bisimulation, we foresee that SMT methods will be the basis for comparison of open transitions.

Scaling up. One important motivation of this work is to attack the complexity of verification of realistic systems by a compositional and parametric approach. Still one may wonder if the price for analysing our symbolic transitions will not make the approach too expensive in term of computing time. We tried a slighly bigger example, assembling 2 Failure controllers. In [START_REF] Qin | Using SMT engine to generate Symbolic Automata -extended version[END_REF], we show that Z3 can check the satisfiablility of a 90K open transitions in a couple of minutes.

α

  ∈ A ::= a(p 1 , . . . , p n ) action terms p i ::= Expr parameters Expr ::= Value | x | op(Expr 1 , . . . , Expr n ) expressions Definition 2 (pLTS). Given a term algebra T Σ,V , a pLTS is a tuple pLT S S, s 0 , → where: • S is a set of states, with s 0 ∈ S the initial state. • →⊆ S × L × S is the transition relation, with L the set of labels of the form α, e b , (x j := e j ) j∈J , where α ∈ A V is a parameterized action, e b ∈ B V is a guard, and expressions E V ∪ A P are assigned to x j . If s α, e b , (xj:=ej ) j∈J -------------→ s then vars(e b ) ⊆ vars(s)∪vars(α), and ∀j ∈ J. vars(e j ) ⊆ vars(s)∧x j ∈ vars(s ).

Fig. 1 :

 1 Fig. 1: The BIP specification of the Failure Monitor architecture

Fig. 2 :

 2 Fig. 2: pNet encoding of the Failure Monitor architectureFor example, the port B.fail will never fire alone in a global state, where both Timer .start and Control .fail are enabled.Application of the Failure Monitor architecture ensures that, whenever a failure is registered in the operand component, the system will be reset, unless a resumption is registered within Max time units (more details in Sect. 5.4).Figure2shows a pNet encoding of the above Failure Monitor architecture. This encoding is structural: each coordinator component is encoded as a pLTS, the operand component-as a hole; connectors of the BIP model are encoded as synchronisation vectors. Each connector that does not involve triggers is trivially encoded by a synchronisation vector comprising the same ports. In order to encode the semantics of the connectors involving triggers, we 1) in the pLTS encoding the coordinator components, add loop transitions to ensure that all ports involved in such connectors are enabled in all states, 2) associate a Boolean value to each of these ports: the original transitions carry the value true (e.g.

Definition 5 (

 5 Open automaton). An open automaton is a structure A = LTS i∈I i , J, S, s 0 , T where: • I and J are sets of indices, • LT S i∈I i is a family of LTSs, • S is a set of states and s 0 ∈ S the initial state, • T is a set of open transitions and, for each t ∈ T , there exist I , J with I ⊆ I, J ⊆ J, such that t is an open transition over LT S i∈I i , J , and S.The states and the shape of predicates in the transitions of an open automaton representing the semantics of a pNet have the following specific structure.

Fig. 3 :

 3 Fig. 3: One of the unsatisfiable open transitions in the Failure Monitor pNet

Fig. 4 :

 4 Fig. 4: The input of the Z3 solver in SMT-LIB language and the output result

Fig. 5 :

 5 Fig. 5: Open Automaton for the Failure Monitor architecture

For

  this example, the tool builds 184 open transitions, whereof 173 are detected unsatisfiable by Z3. The resulting open automaton, with 3 reachable global states (out of the possible 6) and 11 open transitions, is shown in Fig. 5.

Table 1 :

 1 Algebra Presentation: predefined Sorts and Operators Extension for the BIP use-case of Fig. 2 Action FUN_Action_Bool, fail, resume, timeout, reset, start, tick, ask

	Sort	Constructors	Auxiliary Operators
	Bool	true, false	∧, ∨, ¬, =⇒ , =, =
	Action	Synchro, FUN
	Int	0, {i, -i}i∈Nat	-(unary), +, -(binary), ×, ÷ etc.

We omit detailed presentation of this case for the sake of brevity.

This work was partially funded by the Associated Team FM4CPS between INRIA and ECNU, Shanghai