Ariane Ducellier
email: a.ducellier@brgm.fr

Hideo Aochi

Interactions between Topographic Irregularities and Seismic Ground Motion Investigated Using a Hybrid FD-FE Method

Keywords: Finite difference, Finite element, Hybrid methods, Seismic wave propagation, Topographic effects

A hybrid method combining finite element and 4 th -order finite difference techniques is developed to model SH and P-SV seismic wave propagation in a 2D elastic medium with irregular surface topography. Both the classic staggered grid finite difference scheme and the partially staggered grid scheme are tested. The accuracy of the hybrid method is studied by comparison with a semi-analytical and another numerical method. Subsequently, to study the amplification, numerical simulations of seismic wave propagation in a series of hills are carried out and compared with the single-hill case. Depending on the position of the source in relation to the topography, the ratio between the heights and lengths of the hills or the ratio between the lengths of the hills and the wavelength, the presence of several hills as opposed to a single one can increase the amplification effect due to topography. This study highlights the fact that, when evaluating topographic site effects, surrounding topography must be taken into account in addition to local topography.

Introduction

Simulating ground motion at given sites is an essential task in the framework of seismic hazard studies. In most practical cases the source-site distance is large enough and local conditions are fairly simple from a geometrical standpoint; soil behaviour is therefore simulated in a 1D elastic/anelastic column under seismic excitation at the base of the column. However, some cases are much more complex, precluding a 1D approximation. It is known that a 2D or even 3D structure close to the site considered can significantly affect ground motion (e.g. [START_REF] Griffiths | The effect of Appalachian Mountain topography on seismic waves[END_REF][START_REF] Hartzell | Initial investigation of site and topographic effects at Robinwood Ridge, California[END_REF][START_REF] Spudich | Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge, California, earthquake: Implications for mainshock motions[END_REF]. Furthermore, the characteristics of the incident wave also become important factors when the source-site distance is short. To deal with a 2D (and 3D) context, numerical simulations are needed, since the characteristic ground motion features vary on a case by case basis, and it is accordingly difficult to quantify the effect empirically.

Ground motion can be simulated using volumetric methods such as those involving finite differences or finite elements. One of the most popular methods for simulating the seismic wave propagation in an elastic medium numerically is the finite difference method (FDM) (e.g. [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF][START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF]) due to its simplicity, efficiency and ease of implementation. However, it is not capable of simulating complex boundary conditions like surface topography numerically [START_REF] Ohminato | A free-surface boundary condition for including 3D topography in the finite-difference method[END_REF]. This is quite inconvenient, for taking into account sedimentary basins with strongly contrasted impedances or irregular topography is extremely useful when studying site effects for seismic hazard evaluation. As the finite difference method generally uses a uniform grid, accurate modelling of sedimentary layers near the surface with slow wave velocities involves a dramatic increase in the computational dimensions of the model in a large majority of instances (see [START_REF] Aoi | 3D finite-difference method using discontinuous grids[END_REF] for an exception). An alternative method is the finite element method (FEM). It allows the use of irregular grids with elements of different sizes and geometries and is therefore quite useful to simulate complex geometry and to model strong heterogeneities in the velocity structure.

However, this method requires significantly more memory than the finite difference method. A derivative version of the finite element method, of higher order, called the spectral element method (SEM) [START_REF] Komatitsch | The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF], has been developed and is very accurate and faster than any classic finite element method. Nonetheless, the finite difference method is still more attractive from the standpoint of the required computing resource, and it seems difficult to implant in the SEM the non-linear dynamics required in soil dynamics.

Various combinations of numerical methods for simulating seismic wave propagation in an elastic/anelastic medium were proposed in a large number of studies: hybrid discrete wave-number -finite difference method [START_REF] Fäh | Effects of two-dimensional topographies using the discrete wavenumber -boundary integral equation method in P-SV cases[END_REF][START_REF] Zahradník | Hybrid seismic modeling based on discrete-wavenumber and finitedifference methods[END_REF], hybrid boundary element -finite element method [START_REF] Bielak | Domain reduction method for threedimensional earthquake modelling in localized regions. Part I: Theory[END_REF][START_REF] Gatmiri | Seismic site effects by an optimized 2D BE/FE method I. Theory, numerical optimization and application to topographical irregularities[END_REF], and many others. Such a combination confers advantages not provided by a single method on its own. In particular, the combination of finite difference and finite element methods is expected to retain the efficiency of the finite difference method at large scale and to introduce the geometrical flexibility of the finite element method locally. Therefore, in this paper, we try to adapt the hybrid method previously proposed by [START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF] and [START_REF] Ma | Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach[END_REF], to study ground motion resulting from irregular surface topography. We use various finite element and finite difference methods to simulate wave propagation in a 2D elastic medium. The finite difference method makes it possible to model both a point source and an extended fault: we thus will be able to model the effects of the topography on the wave generated by a complex source. Previously, De Martin et al. (2007) were able to simulate the wave propagation sequentially from the FDM to FEM domains by transferring the wave field information along the boundary. This is a practical and convenient way to simulate ground motion locally for the site considered. However, it should be noted that the local wavefield in the FEM does not impact the regional wavefield in the FDM. This assumption and scheme are accurate enough in many cases; however, in this study, we highlight cases where local disturbance in the wavefield is also important at the regional scale. This is why we need a complete hybrid simulation technique. Many earlier studies have numerically modelled seismic wave scatter around a single canyon or hill (e.g. Gaffet and Bouchon 1989, and many others), mainly in order to compare the results of numerical approaches with an analytical or semi-analytical solution. Also, some authors have applied numerical methods to study wave propagation in more complex real cases. However, few systematic studies have been conducted to determine the effects on ground motion of the interaction of a series of irregular topographies.

To do so, we need to understand how the surface irregularities interact with each other. In the following sections, we begin by verifying our numerical approach and continue with several parametric studies. We will discuss how topographic amplification and resonance phenomena appear in various topographical settings.

The numerical method

In this section, we briefly explain the numerical approach of the hybrid finite difference -finite element method following [START_REF] Ma | Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach[END_REF]. The computational domain is divided into three sub-computational areas: one where only a standard finite difference method (FDM) is used to solve the seismic wave equation, a second where only an explicit finite element method (FEM) is used, and a third interface area where both FDM and FEM are used and the finite element mesh is made up of regular square meshes to facilitate interpolation of the wave field between the FD and the FE schemes. Examples of the domain partition are shown in Figures 1a andb. [START_REF] Ma | Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach[END_REF] used the standard staggered finite difference grid [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF][START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF]). As collocation points (stress and velocity) are different, i.e., shifted, in each scheme, the velocity field is assumed to be a 3 rdorder polynomial function and interpolated. This negatively impacts accuracy, but the use of higher-order interpolations might cause instabilities [START_REF] Ma | Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach[END_REF]. This is a weak point in this hybrid method. Although [START_REF] Badea | Schwarz method for earthquake source dynamics[END_REF] do not specify what their finite difference method entails, their previous work [START_REF] Wolf | Finite element -finite difference coupling for elastodynamic rupture propagation[END_REF]) tested different schemes, including a partially staggered finite difference grid [START_REF] Saenger | Modeling the propagation of elastic waves using a modified finite-difference grid[END_REF]. The advantage of this latter grid scheme is that both horizontal and vertical velocities are computed for the same place, so no interpolation is needed. Here, we have adopted both finite difference schemes. As to the finite element scheme, we have tested higher order meshes, that is, spectral elements as well. However, this did not work in the trials we conducted. Indeed, the use of spectral elements results in high-order interpolation schemes in the interface between FDM and SEM areas. We believe that this does not ensure stability in the computation, as was the case in the study by [START_REF] Ma | Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach[END_REF].

The influence of boundary conditions on the results of the computation is well known. For all the model boundaries excepting the free surface, where only finite differences are used, we introduce a convolutional PML (Komatitch and Martin 2007) to avoid artificial reflections. On the free surface, both finite difference and finite element methods are used. In the FD part, we use a stress-imaging method [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[END_REF] to model the free surface. In the FE part, the free surface condition is obtained by imposing a zero value of stress in the variational formulation.

Numerical verification for simple cases

Each part of the hybrid method was numerically verified for various wave propagation settings. Now we need to verify the hybrid method for some simple cases by comparison with analytical solutions and another numerical method.

Figures 2 and3 show classic examples where the plane wave propagates upward in a homogeneous elastic medium towards a semi-circular canyon of radius a (model geometry shown in Figure 1) for which the solution is (semi-) analytically given by an eigenfunction expansion method [START_REF] Trifunac | Surface motion of a semi-cylindrical alluvial valley for incident plane SH waves[END_REF][START_REF] Trifunac | Scattering of plane SH waves by a semi-cylindrical canyon[END_REF] for a SH problem and [START_REF] Wong | Diffraction of P, SV, and Rayleigh waves by surface topographies[END_REF][START_REF] Wong | Effect of surface topography on the diffraction of P, SV, and Rayleigh waves[END_REF] for a P-SV problem) and by a wavenumber method [START_REF] Kawase | Time-domain response of a semi-circular canyon for incident SV, P, and Rayleigh waves calculated by the discrete wavenumber boundary element method[END_REF]) for both cases. We suppose an incident wave in the FDM domain with a second-order Ricker function:

( ) ( ) ( )                 - -                 - - = 2 2 exp 5 . 0 2 p s p s t t t t t t t u π π (1)
where the characteristic angular frequency is given by

p t π ω 2 0 =
, and for the calculation we take

( ) a πβ ω 2 0 =
, where β is the shear wave velocity. The Poisson ratio is taken as 1/3 for the P-SV problem. The size of the finite difference grid spacing is 0.05 * a. The size of the irregular finite element meshes is approximately equal to the finite difference grid spacing (see Figures 1a andb).

The time step is taken as 0.0005 * (a / β). Figure 2 shows the case of a null wave angle of incidence (top panels) with the fully staggered grid (FSG) finite difference scheme (left panel) and the partially staggered grid (PSG) finite difference scheme (right panel), and the case of a 30° anticlockwise wave angle of incidence (bottom panels) for the FSG (left panel) and PSG (right panel) cases. The line shows the numerical results obtained with the hybrid method, while the points represent the results with the (semi-)

analytical method taken from [START_REF] Kawase | Time-domain response of a semi-circular canyon for incident SV, P, and Rayleigh waves calculated by the discrete wavenumber boundary element method[END_REF]. Figure 3 shows the case of a horizontal incident wave (top panels) and a vertical incident wave (bottom panels) with a null wave angle of incidence. The solid and dotted lines show respectively the horizontal and vertical amplifications, while the points represent the results with the (semi-) analytical method taken from [START_REF] Kawase | Time-domain response of a semi-circular canyon for incident SV, P, and Rayleigh waves calculated by the discrete wavenumber boundary element method[END_REF]. As the edges of the canyon are not smooth (discontinuous in their spatial derivatives), some influence does remain, especially outside the canyon in the P-SV case. This occurs because the reflected waves propagate away from the canyon following the given geometry. The same difficulty is encountered even with the discrete wave-number method [START_REF] Kawase | Time-domain response of a semi-circular canyon for incident SV, P, and Rayleigh waves calculated by the discrete wavenumber boundary element method[END_REF]) with respect to Wong's solution (1979;[START_REF] Wong | Effect of surface topography on the diffraction of P, SV, and Rayleigh waves[END_REF]. The SH problem is always easier to reproduce by any numerical methods. We observed that a slightly better result is obtained if a partially staggered grid is used for the FDM as opposed to a fully staggered grid. This may be because the partially staggered grid does not require additional interpolations. However, this difference remains minor.

Further numerical tests are carried out through the comparison with a spectral element method (code EFISPEC2D, De [START_REF] Martin | Verification of a Spectral-Element Method Code for the Southern California Earthquake Center LOH.3 Viscoelastic Case[END_REF]. The simulation of PSV-wave propagation inside a hill is performed (Figure 4). The height of the hill is normalized to 1 and the ratio of hill-height to length is 0.3. Thus the P-wave velocity is 1.732, the S-wave velocity is 1 and the density is 1. At the lower boundary of the medium, a horizontal displacement is imposed for which the time function is given in Equation ( 2),

( ) ( ) ( )         - -         - - = 2 2 2 2 2 2 exp 2 1 2 1 p s p s t t t t t t t u π π π (2)
with t s = 2 and t p = 1. The size of the finite difference grid spacing is 0.05. The time step is 0.0005. The finite element mesh is composed of 5704 nodes and 9912 elements (950 square elements in the FD-FE interface and 8962 triangular elements in the FE unstructured mesh). It corresponds to about 20 nodes per wavelength for a frequency of 1/t p . The boundaries of the finite difference domain are set far enough away from the canyon to avoid artificial reflections. The wave angle of incidence is α = 0° (horizontal motion only). The ratio of the amplitude of the frequency response on the free surface and the incident wave amplitude at main frequency f 0 = 1/t p = 1 is computed for the hybrid method and the spectral element method. Figure 5 shows the displacement amplification at main frequency f 0 versus the horizontal position. The top and bottom panels correspond respectively to the horizontal and vertical components, normalized by the horizontal amplitude for the flat free-surface case. It should be noted that no vertical component exists in the case of flat free surface. The solid line represents the results with EFISPEC2D, while the dotted line represents the results with the hybrid method. Good agreement is observed between the two different numerical methods, thereby validating the implementation carried out in this study.

Application: the interaction of topographic irregularities

Model setting

As explained in the previous section, a single topographic irregularity (a canyon or a hill) has been studied by many researchers, and many realistic applications are also found. Nevertheless, to the best of our knowledge, few studies deal with the interaction of topographic irregularities. This section aims to ascertain how a series of hills interact as compared to a single irregularity (a hill).

Figure 6 shows the model geometry and Table 1 summarizes the model parameters. Henceforth, for the purposes of our discussion, we are designating the hills No. 1 to 7 from left to right. The central hill, No. 4, is located at the middle of the model, i.e., distance 0 km. Each hill is 400 metres long and 200 metres high (see Figure 6), so that the hill length is initially much smaller than the wavelength of around 2100 metres. The series is composed of seven identical hills, with the same geometry. The total length of the topographic irregularity is thus 2800 metres. Although our focus is on a series of irregularities, for the sake of comparison, we also examine the case of a single hill.

The earthquake source is given as a 1250 m-long finite fault in the FDM domain. The source position varies to allow us to study the effect of the angle of incidence on the amplification factor. The source mechanism is given by dip δ = 45° and rake λ = 90° (reverse faulting) on the 2D plane considered (see Figure 6).

Simulation results

Figure 7 shows snapshots of the wave propagation, magnified around the surface topography, for the case where the source is located 5 km to the left of the centre of the model and at a depth of 5 km. At time step 2.5 s, we can observe the P-wave arrival at the level of the hills and at time step 3.5 s the S-wave arrival.

Then, the S-wave is trapped inside the hills by the reflections along the edges of the hills from time step 4.5 s onwards. The wave field is deformed by the reflections even outside the hills (time step 6.5 s). Figure 8 shows the velocity waveforms along the ground surface for the horizontal (top) and vertical (bottom) components. The abscissa represents time and the ordinate, the distance from the centre of the topographical model. We can identify the P and S-wave arrivals and a reflected S-wave (in the negative direction). The ground motion remains significantly high well after the S-wave arrival between -1400 and 1400 metres, where the topography is irregular.

To quantify the interaction of topographic irregularities, we have selected peak ground velocity (PGV) at the surface for our discussion, one of the most important parameters in earthquake engineering. First Figure 9 shows PGV values for different situations. The first three panels represent the effect of the hill location (a single hill corresponding to the positions of Nos. 1, 4 or 7, cf. Figure 6) with respect to the source, then, for the sake of comparison, the case where there is no topography (a flat free surface), as a reference ground motion, and finally the series of seven hills. First of all, we note that every peak of PGV corresponds to each hill, except for the third example (single hill location at No. 7 = 6.2 km from the source), where the hill is located around the node of the seismic wave radiation pattern.

Moreover, the hills close to the source (towards No. 1) seem to have a higher PGV than the hills far from the source (towards No. 7). However, when we look at the ratio of PGVs in the topographic free surface case and PGVs in the flat free surface case (Figure 10), we note that the amplification due to the topography is lower for the hills located close to the source than for the hills located far from the source. Indeed, the PGV ratio is lower for the hills nearer from the source than for the hills farther from the source in this particular geometrical configuration, contrary to what would naturally be expected. We also note a motion deamplification at the right-hand side of the farthest hill from the source (No. 7).

However, no de-amplification occurs at the trough between the hills. In the case of the series of hills, the PGVs on top the hill farthest to the left (No. 1), the one in the centre (No. 4) and the one farthest right (No. 7) are 0.123, 0.107 and 0.065, respectively, whereas in the cases of the single hill in locations corresponding to the former case, the PGVs are 0.127, 0.085 and 0.035 respectively. The motion amplification due to the neighbouring topography is thus very different depending on the geometrical configuration. When there is no topographic irregularity between the source and the site considered, as in the case of the hill farthest to the left (No. 1), the ground motion amplification is negligible. However, when the topographic irregularity (hills) lies between the pass, as in the case of the hill farthest to the right (No. 7), the motion amplification increases (about twice in this example).

The influence of source location (incident wave direction)

The effects of the distance of the hills from the source, of the ridge steepness and of the ratio between wavelength and the lengths of the hills are investigated. All sets of parameters are given in Table 2. In the first set of parameters, the distance between the source and the centre of the topography varies: the source is located successively at 2, 5 and 10 km to the left of the centre of the topography. The focal depth remains at 5 km.

Figure 11 shows horizontal PGV versus distance from the centre of the topography for the single hill located at the centre of the model (top panels) and the series of hills (bottom panels). The hill closest to the source is not necessarily the one that has the highest PGV on top, due to the directivity of the extended fault source.

Figure 12 shows the ratio of horizontal PGV with and without the topography. We can observe that the ratio for the series of hills is higher than for the single hill. In both topographic configurations, the PGV amplification is highest for the source located 5 km from the centre of the topography. Indeed, as the dip of the finite fault is 45°, the hill is then located in the rupture direction. In the case of the source located 2 km from the centre, for the hills nearer the source, the wave incidence is nearly vertical and the PGV ratio is rather low, whereas for the farther hills, the wave incidence angle is higher, inducing a higher PGV ratio due to reflections on the edges of the hills. This phenomenon does not appear for the case with the farthest source.

Finally, Figure 13 shows ground velocity amplitude at frequencies of 0.5 Hz (top panels) and 1.5 Hz (bottom panels), still for the same calculations. For the low frequency, each peak corresponds to the top of one hill, but the variation is rather small compared to the ground-motion tendency over a large scale. For the higher frequency, we observe a sharper response corresponding to the hills together with a de-amplification for some of the hills. The hills affected change according to the source distance.

The influence of hill geometry

In the second set of parameters, the length (L) of the hills is unchanged (400 m) but the height (h) is set successively at 150, 200 and 300 m, corresponding to height-to-length ratios of 0.375, 0.5 and 0.75. The source is located 5 km to the left of the centre of the topography and at a depth of 5 km.

Figure 14 shows horizontal PGV for the single hill located at the centre of the model (top panels) and for the series of hills (bottom panels) for the three heightto-length ratios. In all cases, the spatial variation in PGV is similar, that is, on top of the hills closer to the source, the PGV is higher. Figure 15 shows the ratio of horizontal PGV with and without the topography. Once again, it is found that the ratio for the series of hills is higher than for the case of the single hill. For both topographies, the ratio is amplified when hill height increases. The PGV ratio is lower for the hills nearer the source than for those farther away for this particular source location, contrary to what would naturally be expected. Finally, Figure 16 shows ground velocity amplitude at frequencies of 0.5 Hz (top) and 1.5 Hz (bottom).

The influence of incident seismic wavelength

It is known that the frequency content of the input motion strongly affects site amplification. While it is considered that topographic effects are negligible for long wavelengths, the effects are highest for wavelengths near the horizontal length of the topography. In this study, the length of a single hill is smaller than the wavelength; however, the size of the total topography comprising the seven hills is comparable to the wavelength. We made the ratio between wavelength and the lengths (L) of the hills vary (third model set in Table 2). The finite fault source mechanism, the medium parameters and the hill shapes remain unchanged, but the length of the hills is set successively to 200, 400 and 800 m while the source depth and distance to the left of the centre of the topography are set to 2.5, 5 and 10 km, in order to preserve the proportional geometry with respect to the given frequency. The characteristic source frequency equals the inverse of the duration of the source time history function, i.e., 1 Hz. The characteristic wavelength is thus the S-wave velocity (c) multiplied by the source period (T), that is, equal to 2100 metres. The ratio of the lengths of the hills over this wavelength is thus successively 0.095, 0.19 and 0.38. Figure 17 shows horizontal PGV for the single hill located at the centre of the model (top panels) and for the series of hills (bottom panels) for the three hill-length/wavelength ratios. In all cases, the PGV on top of the hills closer to the source is higher than the PGV on the other hills.

Figure 18 shows ratio of horizontal PGV with and without the topography. As in both previous sets of parameters, the ratio for the series of hills is higher than for the single hill case. The PGV ratio is once again lower for the hills nearer from the source than for the hills farther from the source. Finally, Figure 19 shows ground velocity amplitude at frequencies of 0.5 Hz (top panels) and 1.5 Hz (bottom panels). For the largest hills (largest L/cT), we observe a de-amplification as compared to the two other cases.

Discussion

In all previous parametric studies, we chose a series of seven identical hills. The height of the hill is now constant and set to 100 m, but the lengths of the hills increase from 250 m to 550 m by an increment of 50 m from left to right and from right to left (fourth model set). The source is located 2.5 km to the left of the centre of the topography and at a depth of 2.5 km. Figure 20 shows horizontal PGV versus distance from the centre of the topography on top. On the left panel, the longest hill (550 m long) is located farthest to the right and the shortest hill (250 m long) farthest to the left, and vice versa on the right panel. On the centre panel, all hills are identical (400 m long). The bottom panels show the ratio of horizontal PGV with and without the topography. In the three configurations, we observe a de-amplification for the hills located opposite from the source. For the hills located near the source, there is little difference between the three geometrical configurations.

It is difficult to establish systematic rules about amplification resulting from topography on the basis of this study, as the wave field depends strongly on the given parameters both at earthquake source and local site. The presence of a series of hills can either amplify or de-amplify the PGV as compared to the case of a single hill, depending on the geometrical configuration. Table 3 summarizes the PGV on top of a single hill for each set of parameters, the PGV on top of the hill at the centre of the series of hills and the ratio between these two PGVs.

Simulations with the first set of parameters show that when the topography is aligned with the rupture direction, the signal amplification due to the interaction between the hills increases.

Moreover, the hill where the PGV ratio is the highest can be the farthest from the source or in the middle of the topography, depending on the geometrical configuration. No set of parameters gave a maximal PGV ratio on the hill closest to the source for this particular source model. Nevertheless, we can conclude that, for both the single hill and the series of hills, the PGV ratio is amplified when the ratio between height and length of the hill increases.

It is useful to mention some computing aspects for further applications. In the example of the last section, the number of finite difference grid nodes is set to 220 x 70 = 15400 for all sets of parameters. The number of finite element nodes varies, depending on the size of the hills. For the 200 m high and 400 m long hills, our mesh is composed of 8077 nodes and 10563 elements. With these parameters, our computations last about two hours on a single CPU, dominated by the finite element calculation. However, we did not try to optimize the code nor the size of the meshes. Indeed, we have around 20 grids per wavelength in the FD part of the domain, whereas 5 or 6 grids would be enough to avoid numerical dispersion [START_REF] Moczo | 3D fourth-order staggered-grid finite-difference schemes: Stability and grid dispersion[END_REF]. In the FE part, the number of nodes per wavelength is much higher (see Figure 6), and it would be possible to reduce the number of nodes while maintaining a good accuracy of the computation. For further applications, a key task should be optimizing the finite element meshing and its solver. As the finite element method is powerful not only for the simulation of a complex geometry but also for various nonlinear problems, this numerical approach would be useful.

Conclusion

A hybrid method combining finite element and 4 th -order finite difference was developed to model SH and P-SV seismic wave propagation in a 2D elastic medium with topographic free surface. Both the classic staggered grid finite difference scheme and the partially staggered grid scheme were tested. The accuracy of the hybrid method was studied by comparison with semi-analytical and numerical methods. To illustrate the method, numerical simulations of seismic wave propagation across a series of hills were carried out and compared to the case of a single hill. Depending on the position of the source in relation to the topography, the ratio between the heights and the lengths of the hills or the ratio between the lengths of the hills and the wavelength, the presence of several hills as opposed to a single hill can increase the amplification effect due to the topography. This study highlights the fact that, when evaluating topographic site effects, surrounding topography must be taken into account in addition to local topography.

Figures

Fig. 1a. Mesh used for the simulation with the fully staggered grid. For the area where both FDM and FEM are computed, the finite-element mesh size is half the finite-difference grid spacing.

  Figures 2 and 3 represent the SH and P-SV cases respectively, showing the displacement amplification at main frequency along the ground surface. Along the abscissa we have plotted the distance from the centre of the canyon, normalized by the canyon radius a. On the ordinate axis is the displacement amplification.

Fig. 1b .

 1b Fig. 1b. Mesh used for the simulation with the partially staggered grid. For the area where both FDM and FEM are computed, the finite-element mesh size is the same as the finite-difference grid spacing.

Fig. 2 .

 2 Fig. 2. Displacement amplification of SH waves versus distance from the canyon centre.

Fig. 3 .

 3 Fig. 3. Displacement amplification of P-SV waves versus distance from the canyon centre.

Fig. 4

 4 Fig. 4 Mesh used for the comparison with the spectral elements method.

Fig. 5 .

 5 Fig. 5. Frequency response versus the non-dimensional horizontal coordinate.

Fig. 6 .

 6 Fig. 6. The geometry of the model composed of a series of seven hills.

Fig. 7 .

 7 Fig. 7. Snapshots of the seismic wave propagation for the horizontal (left) and vertical (right) components of velocity.

Fig. 8 .

 8 Fig. 8. Synthetic seismograms of the simulation.

Fig. 9 .

 9 Fig. 9. Horizontal PGV versus distance from the centre of the topography. The three single hills correspond respectively to the positions of Nos. 1, 4 and 7 of the seven-hill case.

Fig. 10 .

 10 Fig. 10. Ratio of PGV versus distance from the centre of the topography

Fig. 11 .

 11 Fig. 11. Horizontal PGV versus distance from the centre of the topography.

Fig. 12 .

 12 Fig. 12. Ratio of PGV versus distance from the centre of the topography.

Fig. 13 .

 13 Fig. 13. Ground velocity versus distance from the centre of the topography for selected frequencies.

Fig. 14 .

 14 Fig. 14. Horizontal PGV versus distance from the centre of the topography.

Fig. 16 .

 16 Fig. 16. Ground velocity versus distance from the centre of the topography for selected frequencies.

Fig. 18 .

 18 Fig. 18. Ratio of horizontal PGV versus distance from the centre of the topography.

Fig. 19 .

 19 Fig. 19. Ground velocity versus distance from the centre of the topography for selected frequencies.

Fig. 20 .

 20 Fig. 20. Horizontal PGV (top) and ratio of horizontal PGV (bottom) versus distance from the centre of the topography.

Fig. 15. Ratio of horizontal PGV versus distance from the centre of the topography.

Fig. 17. Horizontal PGV versus distance from the centre of the topography.

Acknowledgements

We thank H. Kawase and F. De Martin for communicating to us their calculation results for our verification. The comments from the reviewers and the associate editor were useful for improving our manuscript. This work was supported under a BRGM internal research fund program.

Tables

Model parameters Value

Medium


Ricardo Bedin Franc
email: ricardo.bedin-franca@airbus.com

Sandrine Blazy
email: sandrine.blazy@irisa.fr

Denis Favre-Felix
email: denis.favre-felix@airbus.com

Xavier Leroy
email: xavier.leroy@inria.fr

Marc Pantel
email: marc.pantel@enseeiht.fr

Jean Souyris
email: jean.souyris@airbus.com

Formally verified optimizing compilation in ACG-based flight control software

Keywords: Safety critical systems, Optimized code generation, Toolset performance evaluation

This work presents an evaluation of the CompCert formally specified and verified optimizing compiler for the development of DO-178 level A flight control software. First, some fundamental characteristics of flight control software are presented and the case study program is described. Then, the use of CompCert is justified: its main point is to allow optimized code generation by relying on the formal proof of correctness and additional compilation information instead of the current un-optimized generation required to produce predictable assembly code patterns. The evaluation of its performance (measured using WCET and code size) is presented and the results are compared to those obtained with the currently used compiler.

I. INTRODUCTION

Flight Control Software (FCS) development is a very challenging task: not only it has the typical constraints of other software projects, such as budget, delivery schedule and available hardware, but also those that apply to hard realtime, safety-critical software. There are specific regulations for the avionics domain -in particular, the DO-178/ED12 [START_REF]B: Software Considerations in Airborne Systems and Equipment Certification[END_REF] that enforce rigorous development of embedded software in avionics systems.

In this context, it is not trivial to develop software that can perform optimally while complying with all requirements. This paper focuses on the compilation, which is a very important step in the generation of good performance critical software: compilers are very complex tools that carry out the delicate task of generating low-level code that behaves exactly as expected from the high-level source code. Thus, the compiler has clear influence over software performance and safety, but it is commonly developed by third parties (Commercial Off-The-Shelf "COTS" software) and are not necessarily oriented towards the needs of avionics software.

We present experiments carried out by Airbus with the use of a formally-verified compiler, CompCert, in flight control software development. The goal of these experiments is to evaluate its performance in a realistic environment, with software modules similar to actual FCS ones and the same development and verification tools as real FCS. This paper extends the performance evaluation presented in [START_REF] Franc ¸a | Towards Formally Verified Optimizing Compilation in Flight Control Software[END_REF], using a new version of CompCert that includes an annotation mechanism (for traceability and timing analysis purposes) and using more criteria to compare a CompCert-based development with the currently used approach. This paper compares Worst-Case Execution Times (WCET) taking into account the size and function of the modules.

The paper is structured as follows: Section II presents the fundamentals of flight control software development and the method we use to assess software performance. Section III presents the CompCert compiler and its features that led to its choice for this work. Section IV presents the results of its performance evaluation and Section V draws conclusions and research perspectives.

II. FLIGHT CONTROL SOFTWARE, COMPILERS AND PERFORMANCE

A. An Overview of Flight Control Software While older airplanes had only mechanical, direct links between the pilots' inputs and their actuators, modern aircraft rely on computers and electric connections to transmit these inputs. As Traverse and Brière [START_REF] Brière | AIRBUS A320/A330/A340 Electrical Flight Controls: A Family of Fault-Tolerant Systems[END_REF] describe, digital, electrical flight control systems are used on all aircraft control surfaces since the development of the Airbus A320.

Hardware and software used in flight control systems are subject to very strict requirements, just as any other component. Every kind of avionics software is subject to the DO-178 (currently, version B) regulations, and the DO-178B guidelines become more numerous and more stringent according to the criticality of a given software program. The correct operation of flight control software is essential to a safe flight, hence they belong to "software level A" and their planning, development, verification and project management must comply with very strict regulations.

In addition, aircraft manufacturers usually have their own internal development constraints, ranging from additional safety considerations (e.g. dissymetry and redundancy) to industrial ones, such as delivery delays.

B. The Case Study

In order to carry out a realistic study, we have chosen to use a case study that closely resembles actual flight control software: not only the source code is representative (in functionalities and size) of flight control laws, but the target computer is also representative of true flight control hardware, so as to illustrate typical constraints on hardware usage. The hardware and software used in this work are similar to those described in [START_REF] Souyris | Computing the Worst Case Execution Time of an Avionics Program by Abstract Interpretation[END_REF]: the relevant hardware in the scope of this work comprises the MPC755 microprocessor, and an external RAM memory. The MPC755 is a single-core microprocessor, which is much less complex than modern multi-core ones but does have pipelines, an internal cache and superscalar architecture, three elements that make its behavior less predictable. Naive timing analysis of this microprocessor could lead to the "timing anomalies" described by Lundqvist and Stenström [START_REF] Lundqvist | Timing anomalies in dynamically scheduled microprocessors[END_REF].

It must be noted that the choices of hardware and software are led by a combination of factors -besides performance needs, there are other constraints, such as weight, size, power dissipation, cost and -most importantly in the scope of this paper -verifiability. Thus, choosing the MPC755 for this case study is consistent with flight control systems of modern aircraft. The development process described below also reflects the intent of developing deterministic and verifiable software.

The development described in the next paragraphs follows the basic steps that are recommended by the DO-178B: specification, design, coding/integration and verification. One must take into account, though, that specification and design are treated together, due to the highly detailed software specification.

1) Specification and Design: The largest part of the software program -the "application" subset, which contains the implementation of the flight control laws -is specified as a set of sheets with the graphical formalism SCADE, each sheet being composed of interconnected basic operators (addition, filter, etc). There is no "main sheet": they all belong to the same hierarchical level and they communicate via their input and output parameters. In order to simplify the specification and the code generation, all the symbols used in the sheets are custom-made by the developers. SCADE (V6) state machines are not used, mainly for determinism purposes -conditional statements are kept inside some library symbols.

Figure 1 depicts instances of the custom-made symbols ABS, HWACQDSI, BPO and AFDX FSOUT. Each symbol has inputs and outputs, which are connected to their left and right sides, respectively. At the bottom of the HWACQDSI and AFDX FSOUT symbols, there are some "hidden inputs", which have the same semantics as a normal input but are used to underline symbol parameters (in our example, integer constants) that are not related to the data flow. The program also contains a manually-coded part that goes through distinct specification and design phases, but further details about this part are beyond the scope of this paper, as it makes no use of automatic code generation.

2) Coding: The specification is translated to source code by an automatic code generator (ACG). Automatic code generation, when applicable, has the advantages of being less error-prone and offering much smaller coding times than manual code generation. In order to lighten the burden of source code verification activities, the ACG is qualified as a development tool, according to the DO-178B standards for level A software.

In this study, we use C as a source code language, as it is widely used in critical systems and there are many development and verification tools for C-coded critical systems. Each symbol is represented in C as a macro: a "symbol library" (a set of macros) is manually coded in order to implement each SCADE operator, and the ACG-generated code consists in a sequence of macro instantiations that respect the data flow specified in SCADE. A SCADE sheet is represented by a C function that contains a sequence of macro instantiations -all the data-flow constraints are taken into account by the ACG in order to make sequential C programs that are consistent with the parallel SCADE ones. It must be noted that all the sheets are activated periodically, but their activation period may vary. Thus, the execution cycle is divided into several sequential "tasks" and each sheet may be activated in one or more tasks during a cycle.

The C code is finally compiled with a COTS compiler1 and linked to produce an executable file. It must be noted that the compiler -like the vast majority of compilers industrially used -is seen as a "black box" by the development and verification teams. In this case, the object code must be verified thoroughly, and the safest solution to carry out a complete verification taking into account the use of a COTS compiler and the high reactivity of the ACG process is to forbid compiler optimizations in order to force the generation of constant code patterns for each symbol. As our ACG-generated code is a (potentially long) sequence of a limited number of symbols and the symbol library code changes much less often than the application in actual flight control programs, it is less onerous to carry out thorough verification activities over each possible code pattern for this symbols than verifying all code "sheets" in each compilation.

3) Verification: Every development phase must be verified and this verification must meet the DO-178B requirements. As this paper focuses on the compilation, we shall describe the main activities that verify software coding and integration:

• Source Code Verification: The source code must be traceable and compliant to the design (in our case, the SCADE specification). Also, it must respect the software resource and time constraints.

• Object Code Verification: Object code also must be traceable and compliant to the SCADE specification and the integration of software modules, as well as their integration with the target computer, must be verified. The DO-178B demands requirement-based verification: a program must be verified with respect to its high-level and low-level requirements. In this paper, we suppose that lowlevel verification is carried out at symbol level (e.g. tests and/or formal proofs of symbol outputs), hence the compiler must not optimize away the symbol outputs, even if they are, indeed, intermediate results of a function.

Usually, these verification activities (especially object code verification) involve testing. For level A software, the whole code must be tested with Multiple Condition/Decision Coverage (MC/DC) and traceability between source code and object code is necessary to validate code coverage:

• If coverage is measured over the source code, traceability is necessary to ensure that there is no added, unverified functionality in the object code. Typical cases of "added, unverified" functionalities could be found in compilers that add array bound checks or that have a complex management for switch statements.

• The DO-178B report for clarification [START_REF]Final Report for Clarification of DO-178B "Software Considerations in Airborne Systems and Equipment Certification[END_REF] states that if coverage is measured over the object code, traceability is necessary to ensure that the measured coverage is equivalent to MC/DC, as the object code (such as Assembly language) may not contain the multiple conditions found in the source code. Traceability analysis is much less complicated if the object code presents no optimization and no untraceable code added by the compiler. Once again, it is useful to hinder compiler optimizations in order to simplify the verification activities.

C. Estimating Software Performance

Besides being a DO-178B requirement, Worst-Case Execution Time (WCET) analysis is a safe and reliable timing verification in the avionics software context. Hardware and software complexity make the search for an exact WCET nearly impossible; usually one computes time values which are as close as possible to the actual WCET, but always higher than it. In our case study, the main purpose of WCET analysis is to make sure that no application task oversteps its allowed execution time.

As mentioned by Souyris et al [START_REF] Souyris | Computing the Worst Case Execution Time of an Avionics Program by Abstract Interpretation[END_REF], it was once possible to compute the WCET of avionics software by measurement and analysis, but such method is not feasible in state-ofthe-art programs. The current approach at Airbus relies on AbsInt 2 's automated tool a 3 [START_REF] Heckmann | Worst-case Execution Time Prediction by Static Program Analysis[END_REF] to compute the WCET via static code analysis of the executable file. In order to obtain accurate results, the tool requires a precise model of the microprocessor and other relevant components; this model was designed in close cooperation between Airbus and AbsInt.

Sometimes it is important or even essential to give a 3 extra information about loop counts or register value bounds to refine its analysis. As described in [START_REF] Souyris | Computing the Worst Case Execution Time of an Avionics Program by Abstract Interpretation[END_REF], annotations are necessary when memory access address ranges cannot be computed precisely because of limitations in a 3 value analysis (e.g. floating-point). The imprecisions that arise from such limitations degrade WCET analysis and can go as far as stopping a 3 from completing WCET computation. Such a situation is depicted in Algorithm 1: as a 3 is not yet able to carry out the floating-point comparison, it cannot evaluate the range of addresses that may be accessed in line 6. In this case, a 3 has to continue its computation assuming that the access might occur in any memory address, and the great deal of complexity that is added incurs a strongly overestimated -if not unbounded -WCET. For instance, if it is known that variable i is always within the bounds of the array, this information should be provided to a 3 as an annotation.

In our case, annotations are needed only in a few symbols, so as to compute some addresses more precisely -with a this kind of annotation can be assigned only to microprocessor registers, which are depicted in the Assembly code.

Let us assume that Algorithm 1 is part of the C macro of a symbol and that its corresponding (non-optimized) Assembly code is depicted by Algorithm 2. One can notice that the C variable i is stored in r31, since it is loaded with the resulting value of the floating-point to integer conversion. Thus, if we know that i is always between, say, 0 and 9, the annotation should be:

instruction "Checkpoint" + 0x14 bytes is entered with r31 = from 0 to 9;

In order to keep the fast pace of the ACG-based approach (and avoid potential human mistakes), an automatic annotation generator was devised to avoid manual activities and keep the efficiency of the development process. Each symbol that needs annotations will need them repeatedly for all of its instances, but it is not difficult to annotate automatically the correct Assembly lines with a non-optimized compilation that always generates similar code patterns for all instances of the symbol. Whenever the macro containing Algorithm 1 is instantiated, an annotation would be needed at the same offset 0x14 from the tag Checkpoint. Thus, one has to track the possible code patterns for the symbols that need annotations (to make sure that subtle variations in the code patterns do not change the offset of the instruction that needs an annotation) and find the right offsets to assign those value ranges. This annotation strategy is simple and effective, but would not work if the compiler could optimize the code.

Algorithm 2 Example of a loop that needs annotations

Checkpoint: 00 fctiwz f0,f31 04 stfd f0,8(r1) 08 lwz r31,12(r1) ⊲ a3 cannot infer this value 0c addis r11,r0,lookup table@ha 10 addi r11,r11,lookup table@l 14 rlwinm r10,r31,3,0,28 ⊲ we should help a3 here 18 lfdx f30,r11,r10

Annotations are also used in the manually-coded subsets in order to specify -for instance -the behavior of other hardware components, but those are created manually and are not in the scope of this paper.

III. COMPCERT: TOWARDS A TRUSTED COMPILER

One can figure out that, in extremely critical systems, traditional COTS compilers must be used with great caution with respect to code optimization. However, there are recent advances in the compilation field: in the scope of this work, a most promising development is the CompCert3 compiler. Besides working in a more realistic environment (a large C subset as input language, MPC755 as one of the possible target processors) than other experimental compilers, its development is taking into account the needs of critical systems and its own code is available for study if its end users need to know its internal details in order to devise verification strategies for their software.

As described in [START_REF] Leroy | Formal verification of a realistic compiler[END_REF], CompCert is a multiple-pass, moderately-optimizing compiler that is mostly programmed and proved correct using the Coq proof assistant. Its optimizations are not very aggressive, though: as the compiler's main purpose is to be "trustworthy", it carries out basic optimizations such as constant propagation, common subexpression elimination and register allocation by graph coloring, but no loop optimizations, for instance. As no code optimizations are enabled in the currently used compiler, using a few essential optimization options could already give good performance benefits.

The semantic preservation proof of CompCert guarantees that the generated code behaves as prescribed by the semantics of the source program. The observed behaviors in CompCert include termination, divergence and "going wrong". To strengthen the preservation theorem, behaviors also include a trace of the input-output operations performed during the execution of the program. Input-output operations include system calls (if an operating system is used) as well as memory accesses to global variables declared "volatile" (corresponding in particular to memory-mapped hardware devices). The formal verification of CompCert proves, in effect, that the source program and the generated machine code perform the same input-output operations, in the same order, and with the same arguments and results.

A. CompCert annotation mechanism

To strengthen the guarantees implied by CompCert's formal verification, we have introduced a generic program annotation mechanism enabling programmers to mark source program points and keep track of the values of local variables at these points. Syntactically, annotations are presented as calls to a compiler built-in function, taking a string literal and zero, one or several program variables as arguments: __builtin_annot("x is %1 and y is %2", x, y);

The formal semantics of this statement is that of a pro forma "print" statement: when executed, an observable event is added to the trace of I/O operations; this event records the text of the annotation and the values of the argument variables (here, x and y). In the generated machine code, however, annotations produce no instructions, just an assembler comment or debugging information consisting of the text of the annotation where the escapes %1, %2 are replaced by the actual locations (in registers or memory) where the argument variables x, y were placed by the compiler. For example, we obtain # annotation: x is r7 and y is mem(word,r1+16) if x was allocated to register r7 and y was allocated to a stack location at offset 16 from the stack pointer r1.

Despite executing no instructions, this special comment is still treated, from the standpoint of formal semantics, as a pro forma "print", generating an observable event.

The semantic preservation proof of CompCert therefore guarantees that annotations are executed in the same order and with the same argument values both in the source C program and in the generated assembly code.

A typical use of annotations is to track pieces of code such as library symbols. We can put annotations at the beginning and the end of every symbol, recording the values of the arguments and result variables of the symbol. The semantic preservation proof therefore guarantees that symbols are entered and finished in the same order and with the same arguments and results, both in the source and generated codes. This ensures in particular that the compiler did not reorder or otherwise alter the sequence of symbol invocations present in the source program -a guarantee that cannot be obtained by observing systems calls and volatile memory accesses only.

This possibility of finer-grained semantic preservation is most welcome, since some of our verification activities may be carried out at symbol level and semantic preservation needs to be ensured at this level to be useful in our context. In particular, we consider using per-symbol annotations in order to generalize the results of symbol-based tests: the test results for a given symbol remain valid for all possible code patterns generated when instantiating this symbol. This approach is currently under discussion and such discussions are not in the scope of this paper.

Another use of annotations is to communicate additional information to verification tools that operate at the machine code level, such as the WCET analyzer of the a 3 tool suite. Continuing the example of section II-C, we insert a sourcelevel annotation as shown below.

During compilation, this source-level annotation is turned into a special comment in the generated assembly file, where the placeholder %1 is replaced by the machine register containing variable i. Algorithm 4 below shows the assembly code generated by CompCert for two successive instantiations of the symbol containing Algorithm 3.

The two instantiations generate significantly different assembly code fragments, since the second instantiation reuses From these special comments and their locations in the assembly listing, an automatic tool can easily extract the information that at points 20 and 40 from the beginning of the current function, register r3 (holding the array index) is in the range [0, 9], and communicate this information to the WCET analyzer. Some aspects of this annotation mechanism are still under discussion with the CompCert and a 3 developers, but an experimental annotation generator has already been developed and the ease of its development is a testimony to the usefulness of the CompCert annotation mechanism: readilyavailable, formally-verified variable information simplify the task of automating annotation generation for a 3 . One should remember that, in comparison, the annotation generator for the "default" compiler code must be reconfigured for each symbol library change: a new analysis must be carried out in order to verify which are the possible Assembly patterns for all symbols that need annotations, and which are the offsets that need these annotations.

IV. PERFORMANCE EVALUATION OF COMPCERT

The evaluation environment is essentially the same as in our previous work [START_REF] Franc ¸a | Towards Formally Verified Optimizing Compilation in Flight Control Software[END_REF] and is depicted in Figure 2. CompCert is used only to generate Assembly code from the ACG-coded files, as these files are by far the most voluminous part of the program. Compilation of other software subsets, assembling and linking were done with the compiler, assembler and linker that are used in actual FCS. In order to ensure greater realism in the experiments, about 3600 files that are functionally equivalent to a whole flight control program were compiled with CompCert 1.9. These files represent about 3600 SCADE sheets -when compiled with the default compiler, they correspond to 3.96 MB of Assembly code. The symbol library that was used comprises 145 symbols whose sizes vary from less than 10 to more than 100 lines of code. CompCert's source-level annotation mechanism was used to track symbols' inputs and outputs, and also to generate additional information for some variables that need range annotations. As explained in section III-A, this information is available in the generated assembly files, which are examined by the annotation generator to produce an annotation file in the suitable format for a 3 . a 3 was used to compute WCET at two different levels: the most important benchmark is at task level, as it is the measure used for timing analysis in actual programs. While a traditional WCET analysis consists in verifying that each task is performed within its allocated time, we opted to compare the average WCET of all tasks in order to have a synthesis of the results for every task. In addition, we analyze individually the WCET of all SCADE sheets: we do not seek interprocedural optimizations or a register allocation that goes beyond one single module, hence individual WCET computations are meaningful in this context and are useful to find out which kind of algorithms get the most of CompCert's optimizations. The baseline for the benchmark is the WCET of an executable file generated with the default compiler and the compilation options used in a real flight control program. Some analyzed sheets instantiate symbols that need range annotations; CompCert's annotation mechanism was used together with a simple annotation generation script to assign variable ranges for a 3 when needed. Figure 3 In addition to WCET computations, code size measures were carried out as an auxiliary performance indicatorsmaller code size often means better-performing code.

The results of the WCET analysis are quite encouraging, as the average WCET improvement per task was 10.6%, which is a significant improvement by flight control software standards. As already pointed out in [START_REF] Franc ¸a | Towards Formally Verified Optimizing Compilation in Flight Control Software[END_REF], this is mainly due a better register allocation that saves many loads and stores that had to be performed to keep symbol inputs and outputs on stack.

Figure 4 depicts the WCET computed for all sheets, ordering them according to the WCET obtained when they were compiled with the default compiler. The WCET improvement may change from one region of the graph to another (modules with a very low or very high WCET do not always have a visible improvement, whereas CompCert clearly improved the WCET of those in the middle part of the curves) and even inside a region -the WCET curve for CompCert-compiled modules is not smooth. In order to refine the general results obtained by the analysis of this large number of files, special attention was dedicated to files that had extreme values of WCET and code size. The 10% longest and shortest files (in WCET or code size) had results that differed from the average and had specific statistics in order to underline those differences. In addition, some "unexpected" results (e.g. the default compiler performing better than CompCert) were analyzed individually.

The WCET of all analyzed modules was computed for the executable files generated by both compilers, in order to compare them when compiling modules of various WCET and code size value ranges -using the benchmark WCETs and code sizes to classify the modules into categories. The main conclusions from these experiments are:

• In the analyzed program, even if a module is small, there is usually some possibility of optimization but results may vary according to the symbols that are instantiated in a given module. Some symbols have their code vastly improved by CompCert, whereas -in some very exceptional cases -the WCET of a module rises due to the overhead caused by longer function prologues and epilogues. In fact, modules that present very small code size are not quite a reliable source of WCET analysis because even their address in memory becomes a significant factor in WCET analysis.

• Sheets that are not among the fastest or slowest have a slightly better WCET improvement than the overall results. This shows that the optimizations work best when there is enough code to enable their full use, but the code is still compact enough to avoid register spilling.

• A sheet can have a large WCET for two main reasons: either it may have many instructions to execute or it may contain interactions with hardware devices that are time-consuming. In the former case, CompCert usually performs better, except when dealing with spilled variables -the gains become less significant because spilled variables resemble variables compiled with the default compiler. CompCert optimizations can do little or nothing to improve the WCET of a sheet if its symbols spend most of their computation time doing hardware acquisitions and emissions. In our case study, it is more common to have interactions with hardware than register spilling, hence the WCET gain over "long" sheets (larger code size) is more pronounced than the gain over "slow" ones (higher WCET).

• Even with its optimizations turned off, the default compiler sometimes succeeds in selecting more efficient combinations of PowerPC instructions than CompCert. An example is address computations for volatile memory accesses, which CompCert compiles rather naively. We plan to improve the instruction selection phase of CompCert to reduce these inefficiencies.

A. Verification considerations

Since the main reason to avoid most optimizing compilers is the ensuing difficulty to verify traceability and compliance of the object code, the performance evaluation was followed by a study of possible verification strategies that could use CompCert's semantic preservation in order to meet the DO-178B requirements without losing the performance gains obtained with its optimizations. This study is currently under way but it is already clear that the "traditional" analysis mentioned in [START_REF]Final Report for Clarification of DO-178B "Software Considerations in Airborne Systems and Equipment Certification[END_REF] to verify traceability between source code and object code is still feasible with CompCert, as its optimizations remove computational instructions but do not change significantly the code structure (branches, etc). Also, its semantic preservation theorem could be used as a strong argument for traceability and compliance between source code and object code.

V. CONCLUSIONS AND FUTURE WORK

This paper presented an evaluation of the CompCert compiler, based on the characteristics of Airbus flight control software: ACG-based code, modules with different characteristics. Even if we focused on its WCET analysis to assess its performance, CompCert's formal proofs are already seen as a key to bring more confidence in the compilation process, helping to make safe use of code optimizations. Moreover, CompCert's optimizations apply to the vast majority of the modules that are representative of FCS.

The main ongoing work in our CompCert study is the development of a new verification strategy that must be at least as safe as the current one. It is a complex subject on its own but some conclusions drawn from it (e.g. the need for semantic preservation at symbol level) are already being taken into account -as it is likely that we will need semantic preservation at symbol input level, the performance measures were taken using library symbols endowed with CompCert's annotations to preserve the semantics of their inputs and outputs. An important discussion point is the DO-178 interpretation of a tool like CompCert.

The performance evaluation shall not stop at the current state. As the symbol library was coded bearing in mind the current compilation strategy, an interesting work will be recoding it in order to favor optimizing compilation, with fewer intermediate variables and use of Small Data Areas. It is likely that the obtained WCET will be lower and every percent counts if one intends to improve performance.

Another direction for future work is to further improve WCET by deploying additional optimizations in CompCert and proving that they preserve semantics. The WCC project of Falk et al [START_REF] Falk | A compiler framework for the reduction of worst-case execution times[END_REF] provides many examples of profitable WCET-aware optimizations, often guided by the results of WCET analysis. Proving directly the correctness of these optimizations appears difficult. However, equivalent semantic preservation guarantees can be achieved at lower proof costs by verified translation validation, whereas each run of a non-verified optimization is verified a posteriori by a validator that is proved correct once and for all. For example, Tristan and Leroy [START_REF] Tristan | Formal verification of translation validators: A case study on instruction scheduling optimizations[END_REF] show a verified validator for trace scheduling (instruction scheduling over extended basic blocks) that could probably be adapted to handle WCC's superblock optimizations. Rival has experimented the translation validation approach on a wider scope in [START_REF]Symbolic transfer functions-based approaches to certified compilation[END_REF] but, currently, the qualification and industrialization of such a tool seems more complex.

Figure 1 .

 1 Figure 1. SCADE custom-made operators

Figure 2 .

 2 Figure 2. The development chain of the analyzed program

  depicts the flow of annotation data, from the C macros (where the necessary extra information is specified by the user) to the execution of a 3 .

Figure 3 .

 3 Figure 3. Automatic annotation generation for a 3

Figure 4 .

 4 Figure 4. Overall result of WCET comparison

Table I

 I summarizes the WCET analysis results.

		WCET (CompCert) Size (CompCert)
	All application tasks	-10.6%	-13.8%
	Small code size sheets	-2.0%	-14.6%
	Small WCET	-10.6%	-12.9%
	Average WCET	-12.6%	-14.3%
	Average code size	-10.7%	-13.6%
	Large code size	-7.7%	-14.2%
	Large WCET	-3.8%	-12.4%
		Table I
	CODE SIZE AND WCET COMPARISON

For confidentiality reasons, the currently used compiler, linker and loader names are omitted.

http://compcert.inria.fr