Cristina Silvano

Gianluca Palermo

Giovanni Agosta

Amir H Ashouri

Davide Gadioli

Stefano Cherubin

Emanuele Vitali

Luca Benini

Andrea Bartolini

Daniele Cesarini

Andrea Bartolini

Joao Cardoso

Joao Bispo

Pedro Pinto

Riccardo Nobre

Erven Rohou

Loïc Besnard

Imane Lasri

Nico Sanna

Carlo Cavazzoni

Radim Cmar

Jan Martinovič

Kateřina Slaninová

Martin Golasowski

Andrea R Beccari

Candida Manelfi

An- Drea Bartolini

Autotuning and Adaptivity in Energy Efficient HPC Systems: The ANTAREX Toolbox

Keywords: High-Performance & Power-aware Computing, Autotuning, DSL

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Designing and implementing HPC applications are difficult and complex tasks, which require mastering several specialized languages and tools for performance tuning. This is incompatible with the current trend to open HPC infrastructures to a wider range of users. The current model where the HPC center staff directly supports the development of applications will become unsustainable in the long term. Thus, the availability of effective standard programming languages and APIs is crucial to provide migration paths towards novel heterogeneous HPC platforms as well as to guarantee the ability of developers to work effectively on these platforms. Moreover, to fulfil the target of energy-efficient supercomputers there is the need to combine the current trend of HPC systems with radically new software stacks to exploit the benefits offered by heterogeneity at several levels (supercomputer, job, node).

ANTAREX project addresses these challenging problems through a holistic approach spanning all the decision layers composing the supercomputer software stack and exploiting effectively the full system capabilities (including heterogeneity and energy management) [START_REF] Silvano | The ANTAREX approach to autotuning and adaptivity for energy efficient HPC systems[END_REF][START_REF] Silvano | AutoTuning and Adaptivity appRoach for Energy efficient eXascale HPC systems: the ANTAREX Approach[END_REF]. The main goal of the ANTAREX project is to express by a DSL the application self-adaptivity and to runtime manage and autotune applications for green heterogeneous HPC systems up to Exascale. One key innovation of the proposed approach consists of introducing a separation of concerns (where self-adaptivity and energy efficient strategies are specified aside to application functionalities) promoted by the definition of a DSL inspired by aspect-oriented programming concepts for heterogeneous systems. The new DSL will be introduced for expressing at compile time the adaptivity/energy/performance strategies and to enforce at runtime application autotuning and power management. The goal is to support the parallelism, scalability and adaptivity of a dynamic workload by exploiting the full system capabilities (including energy management) for emerging large-scale and extreme-scale systems. Moreover, ANTAREX enables the performance/energy control capabilities by introducing software knobs (including application parameters, code transformations and code variants), and designs scalable and hierarchical optimal control-loops capable of dynamically leveraging them together with performance/energy control knobs at different time scale to always operate the supercomputer and each application at the most efficient operating point. The key ANTAREX software innovations are designed and engineered since the beginning to be scaled-up to the Exascale level.

The ANTAREX project is driven by two use cases coming from highly social relevant HPC application scenarios: (i) a biopharmaceutical HPC application for drug discovery and (ii) a self-adaptive navigation system for smart cities.

This paper presents an overview of the toolbox developed in the context of the ANTAREX project. In particular, Section 2 presents an overview of the ANTAREX framework, while Sections 3 and 4 describe respectively the modules to be used at compile-time and

ANTAREX TOOLBOX OVERVIEW

The ANTAREX approach and related tool flow (shown in Figure 1) operates both at compile-time and runtime. The application functionality is expressed through C/C++ code (possibly including legacy code), whereas the extra-functional aspects of the application, including requirements and adaptivity strategies are expressed by the DSL developed in the project (namely LARA). This guarantees a separation of concerns between functional and extrafunctional aspects, thus facilitating the code reuse without weighing down the application developer with optimization details. The target optimization space of the ANTAREX toolbox is represented on the right side of Figure 1 with cloud symbols: compiler flags and custom transformations, application software knobs including variable approximation, OpenMP parameters, custom precision alternatives, and Voltge/Frequency levels. To enhance the code thus enabling the possibility of those target optimizations, the Clava source-to-source compiler (see Section 3.1) driven by LARA aspects has been used. This guarantees to the user the possibility to adopt several tools with their own interface by using a single front-end language. The outcome of this phase is not only an enhanced version of the code including tunable knobs and observable points, but also a complete building system for putting all together the several components. The application autotuning is postponed at the runtime phase, where software knobs (application parameters, code transformations and code variants) are configured according to the runtime information coming from the execution environment. The dynamic autotuning phase is based on profile-time information collected to build the application knowledge. Finally, an infrastructure-level monitoring framework is used at run-time to collect information about resource usage and application behavior. The same monitoring framework is also used to build automatically an infrastructure-level powermodel used by the power manager to control and optimize the power usage.

COMPILE TIME FRAMEWORK

This section provides an overview of the main components developed in the ANTAREX compile-time framework (see Figure 1).

Lara DSL and Clava Compiler

One of the main challenges of ANTAREX is to develop a DSL able to integrate several technologies developed during the project [START_REF] Silvano | ANTAREX: A DSL-based Approach to Adaptively Optimizing and Enforcing Extra-functional Properties in High Performance Computing[END_REF]. Given that the use cases and the technologies are implemented in C/C++, we took a two-pronged approach regarding the DSL: 1) create a DSL that could transform and adapt C/C++ source-code, and 2) develop DSL API libraries, one for each technology.

Leveraging some previous work done around the LARA framework, and using the Clang compiler as a parser, we developed the tool Clava during ANTAREX. Clava is a C/C++ source-to-source compiler based on the LARA framework [START_REF] Cardoso | Performance-driven instrumentation and mapping strategies using the LARA aspect-oriented programming approach[END_REF] that executes scripts written in the LARA language, a JavaScript-compatible language created to describe source-code analyses and transformations. The LARA language adds a few keywords to JavaScript such as select, to select points in the given source code (function, loop), or apply, to perform analyses and transformations over these points. However, the LARA language (and related framework) are languageindependent, and they are decoupled from the target language that will be analyzed and transformed. They depend on a specific implementation for each target language, and Clava provides that implementation for C/C++.

For each module in Figure 1 (e.g., mARGOt, LibVC) we developed an API in LARA2 that can be used to integrate that specific module in a target source-code. The API for each module has two purposes: 1) to provide a coherent interface for using each ANTAREX technology and 2) to avoid the users from learning the details of each technology for using it. Currently, it is possible, from a single LARA script, to transform automatically a "clean" C/C++ application to use any of the technologies developed by ANTAREX project. To better leverage the efforts spent in Clava and improve its usability, it has been developed a dedicated IDE including also a documentation generator, CMake integration, and a testing framework.

Dynamic Compilation

Dynamic compilation, or the ability to (re-)compile at runtime part of the application code, is potentially useful in HPC applications when profiling the application for code specialization opportunities would take too long, or when data characteristics may change over time, making a static optimization less effective. Dynamic compilation is a well-known technique in the general purpose computing world, but it is difficult to apply in a controlled way as required in HPC. To allow application developers to control dynamic compilation of multiple versions of a given code fragment, we introduce libVersioningCompiler 3 (abbreviated libVC), an open-source C++ library. libVC can dynamically compile any C/C++ code without the burden of a full dynamic compiler, and can be used to support split compilation, continuous optimization and code specialization based on the input data and on workload changes [START_REF] Cherubin | libVersioningCompiler: An easy-to-use library for dynamic generation and invocation of multiple code versions[END_REF]. libVC has been integrated in the ANTAREX tool flow through several LARA aspects covering the most common usage of the dynamic compilation framework.

Memoization

Memoization is a well known technique to improve program execution time by exploiting the redundancy in program execution [23] [22]. The main idea is to save in a table the results of execution of a section of program so that future execution of the same section with same input set can benefit from the saved results. In memoization schemes, two conditions are important: (1) memoized code should not cause any side-effect; (2) memoization should always produce the same result for the same inputs. At function-level, by applying these two conditions, memoizable functions are deterministic functions that do not modify the global state of the program. In ANTAREX, we implement the works done in [START_REF] Suresh | Compile-Time Function Memoization[END_REF] with some extensions. Considering a memoizable C function foo, the memoization consists of 1) the insertion of a wrapper function foo_wrapper and an associated table, 2) the substitution of the references to foo by foo_wrapper in the application. foo_wrapper calls the foo orignal function when the result of a call is not in the table.

In the ANTAREX project, this technique has been extented for C++ memoizable methods that require to take into account the mangling technique and the references to the objects. Moreover, several parameters have been added to improve several memoization approaches that can be addressed by using runtime autotuning or statically: management of the memoization table (updating policies, size of the table, approximation on the input arguments, pre-load table), dynamical "stop/run" of the memoization. Following the philosophy of the project, the memoization strategies are provided to the user by means of a set LARA aspects. An automatic detection of the potential functions/methods to be memoized is also provided to facilitate the integration of the proposed technique.

Compiler Autotuning

A major challenge in automatically choosing the right set of compiler optimizations is the fact that these code optimizations are programming language, application, and architecture dependent. Additionally, the word optimization is a not appropriate given that there is no guarantee that the transformed code will perform better than the original one. In fact, aggressive optimizations can also cause degradation of the performance of the code. Understanding the behavior of the optimizations, the perceived effects on the source-code, and the interaction of the optimizations with each other are complex modeling problems. This understanding is particularly difficult because compiler developers must consider hundreds of different optimizations that can be applied during the various compilation phases.

In ANTAREX, we faced the problem by developing a compiler autotuning framework based on machine learning techniques [START_REF] Ashouri | A Survey on Compiler Autotuning using Machine Learning[END_REF][START_REF] Ashouri | Automatic Tuning of Compilers Using Machine Learning[END_REF]. In particular, we addressed the problem for GCC [START_REF] Ashouri | COBAYN: Compiler Autotuning Framework Using Bayesian Networks[END_REF][START_REF] Gadioli | SOCRATES -A Seamless Online Compiler and System Runtime AutoTuning Framework for Energy-Aware Applications[END_REF], limiting the problem to compiler flag selection, and for LLVM [START_REF] Ashouri | MiCOMP: Mitigating the Compiler Phase-Ordering Problem Using Optimization Sub-Sequences and Machine Learning[END_REF][START_REF] Nobre | A Graph-Based Iterative Compiler Pass Selection and Phase Ordering Approach[END_REF], where ordering and replication of optimizations enter into the game. In particular, the developed framework is based on an application characterization in terms of features and by off-line learned models to predict the optimal configurations for the specific case.

RUNTIME FRAMEWORK

This section provides an overview of the main ANTAREX components of the runtime framework (see Figure 1). Even if these tools work at runtime, we can find them even at compile time in 1 because they are integrated in the code by using the DSL and source-to-source compiler.

Exascale Monitoring Framework

Today processing elements can monitor performance efficiency by inspecting the activity of microarchitectural components as well as physical parameters (such as power consumption and temperature). This is possible throughout HW performance counters which in x86 systems can be read by privileged users, thus creating practical problems for user-space libraries to access them. Besides sensors which can be read directly from the software running on the core, supercomputing machines embed sensors external to computing elements, but relevant to the overall energy-efficiency. These sensors are used to monitor the node and rack cooling components as well as environmental parameters such as room and ambient temperature.

In ANTAREX, we developed EXAMON4 [START_REF] Beneventi | Continuous learning of HPC infrastructure models using big data analytics and in-memory processing tools[END_REF] (standing for EXAscale MONitoring) for accessing the performance and power monitoring sensors distributely. EXAMON decouples the sensors' reading from the sensor value usage. Indeed, EXAMON uses a scalable approach where each sensor has associated a sensing agent which periodically collects the metrics and send the measured values with a synchronized time-stamp to an external data broker. The data broker organizes the incoming data in communication channels with an associated topic, while maintaining a list of subscribers to these topics to which broadcast each new incoming message on a specific topic. The subscriber registers a callback function to a given topic which gets called each time a new message is received. EXAMON uses a state-of-the-art big-data technologies to let the other layers of the ANTAREX toolchain to take advantage of the collected data. To combine EXAMON monitoring services with application knowledge, we developed a fine-grain profiling library capable of collecting in real-time architectural information during MPI phases for each MPI process independently, thus enabling application-aware power management solutions.

Dynamic Autotuning: mARGOt

Main goal of mARGOt is to enhance an application with an adaptive layer to tune software knobs to satisfy application requirements at runtime. In ANTAREX, we consider an application as a parametric function that elaborates input data to produce an output (i.e. o = f (i, k 1 , . . . , k n ) ), with associated extra-functional requirements. In this context, the parameters of the function (k 1 , . . . , k n ) are software knobs that alter the behavior of the application (such as the parallelism level or the number of trials in a MonteCarlo simulation). Even if the functional behavior of the application is utterly under the control of the developer, extra-functional properties depend on the execution context. In particular, they usually depend on the underlying architecture, the actual input, the system workload and configuration (e.g. Dynamic Voltage Frequency Scaling) and the software knob configuration. Moreover, the application requirements might change at runtime and they typically address properties that are in contrast with each other, thus leading to tradeoffs. For these reasons, selecting the most suitable configuration is a complex task.

To achieve this goal, the mARGOt5 dynamic autotuning framework is rooted in the MAPE-K feedback loop [START_REF] Jeffrey | The vision of autonomic computing[END_REF]. In particular, it relies on the application knowledge, derived either at deploy time or at runtime, that states the expected behavior of extra-functional properties of interest. On one side, mARGOt to adapt uses runtime observations as feedback information for reacting to the evolution of the execution context. On the other side, it considers the features of the actual input to adapt in a more proactive fashion. Moreover, the mARGOt framework is designed to be flexible, defining the application requirements as a multi-objective constrained optimization problem that might change at runtime.

Power Management of HPC Resources

Currently, in HPC systems, computing elements and nodes are power constrained. For this reason, state-of-the-art processing elements embed the capability of fine-tuning their performance to control dynamically their power consumption. This includes dynamic scaling of voltage and frequency, power gating of the main architectural blocks of the processing elements, but also some feedback control logic to keep the total power consumption of the processing element within a safe power budget. This logic in x86 systems is named RAPL. However, demanding the power control of processing elements entirely to RAPL might not be the best choice. Indeed, it has been recently discovered that RAPL is application-agnostic and thus it tends to waste power in application phases which exhibit IOs or memory slacks. Under these circumstances, there are operating points more energy efficient than the one selected by RAPL while respecting the same power budget [START_REF] Cesarini | Benefits in Relaxing the Power Capping Constraint[END_REF]. However, these operating points are only viable if the power capping logic is aware of the application requirements. Therefore, we developed a new power capping runtime based on a set of user space APIs to define the priority for tasks currently in execution on a given core [START_REF] Bartolini | Self-Aware Thermal Management for High Performance Computing Processors[END_REF][START_REF] Cesarini | Prediction horizon vs. efficiency of optimal dynamic thermal control policies in HPC nodes[END_REF]. In parallel applications running on supercomputing systems, MPI primitives are used for exchanging data and for synchronizing processes. During synchronization primitives, significant power is waste without doing useful computation. In ANTAREX, we extended the power manager to be aware of these synchronization primitives without impacting performance.

USE CASE EVALUATION

This section highlights two application scenarios and the benefits introduced by applying the ANTAREX toolbox. Reported results will be further improved during the last six months of final validation.

Molecular Docking for Drug Discovery

The drug discovery process involves several tasks to be performed in vivo, in vitro and in silico. Molecular docking is a task typically performed in silico. It aims at finding the three-dimensional pose of a given molecule when it interacts with the target protein binding site. This task is often used for virtual screening a huge set of molecules to find the most promising ones, which will be forwarded to the later stages of the drug discovery process. Given the huge complexity of the problem, molecular docking cannot be solved by exploring the entire space of the ligand poses. This opens a lot of possible directions to address the molecular docking problem.

In the last ten years, Dompe' has undertaken an action in partnership with CINECA, to develop a commercial product, named LiGen [START_REF] Beccari | LiGen: a High Performance workflow for chemistry driven de novo design[END_REF], to perform in silico drug design. In Ligen, the most computational intensive module is LiGenDock [START_REF] Beato | Use of experimental design to optimize docking performance: The case of ligendock, the docking module of ligen, a new de novo design program[END_REF], that implements the molecular docking phase. The main ANTAREX challenges for this use case are: [START_REF] Ashouri | MiCOMP: Mitigating the Compiler Phase-Ordering Problem Using Optimization Sub-Sequences and Machine Learning[END_REF] to develop an energy and resource efficient docking algorithm that adapts its configuration to the execution environment and (2) to scale-out with the ligand database size.

The ANTAREX contributions to this use case are related to several directions:

(i) We developed a drug design mini-app based on LiGenDock, named GeoDock-MA. It is licensed by DOMPE' free-of-charge for no profit purposes to serve as a demo and benchmarking tool for the drug design domain. The mini-app attempts to capture key computation kernels of the molecular docking application for drug discovery implemented in LiGenDock. By developing GeoDock-MA in parallel with the new version of LiGenDock, application developers can collaborate with system architects and domain experts to evaluate alternative algorithms that can either better satisfy end-user constraints, or better exploit architectural features.

(ii) We applied application autotuning approaches to the denovo computer accelerated drug design based on LiGenDock. The goal was to demonstrate the possibility to optimize both energyto-solution and time-to-solution by means of an improved task processing strategy and exploiting energy-awareness. We exploited the low-intrusiveness and low-overhead features of EXAMON to keep track continuously of the application behavior and its effect on the target architectures. This has been used to support manual and automatic tuning and to deliver application components to available resources that best fit the computational complexity of each task. Figure 2 shows the power and frequency profile of the use case application running on 4 nodes (8 sockets) with and without finegrain power management optimizations. From the results, we can see how the power managemement library is capable of reducing the frequency and thus the power cost during load unbalanced phases.

(iii) We enhanced LiGenDock and the related mini-app with a set of software knobs, derived from approximate computing approaches, for trading off quality of the result and performance of the computation. The trade off has been used to control time-tosolution in the virtual screening task. Figure 3 shows the accuracyperformance trade-off extracted by changing approximation knobs at the application-level. The overlap is the cost function used to evaluate each molecule, the original version is the one in the bottom-left corner showing no overlap degradation. In particular, we can observe that the introduced parameterization, even when considering limited degradation, enables a large speedup in terms of application throughput. The new version of the docking application can be configured in different ways according to the goal of the virtual screening task: a faster and less precise way for large screening tasks or a slower and more precise way when targeting a more focused search.

(iv) We implemented an application model to estimate dynamically the task complexity and to inform the underlying system (and the user) about expected values in terms of time-to-solution or energy-to-solution metrics. In particular, this model provides the capability to tune the application accuracy, and thus the target knobs, according to a specified budget.

(v) Finally, the deep understanding of the application code offered us the possibility to build, aside of LiGenDock, a new module for molecular docking that is based only on geometric evaluations. This reduces a lot the cost of the docking process by enabling the user to increase the size of the ligand database to be processed during the virtual-screening task.

HPC-based Car Navigation System

Current navigation systems can be classified in two main categories:

(1) those which can provide navigation in an offline mode (maps are stored on a device, route calculation is done locally on a device), and (2) those, which extensively and intensively use a server side, called online navigation (both for map loading and route calculation). The main disadvantage of the offline navigation is that it might suffer from not having the maps up-to-date, as well as from not knowing the online traffic situation. The other disadvantage is that the routing algorithm is limited in terms of performance capacity, thus leading to sub-optimal routes, thus heuristics with limited search space must be used. The main disadvantage of online navigation is that the Internet connectivity is not always stable, especially for moving objects, which often leads to slower responsiveness, and might even lead to a total failure of navigation. The best result can obviously be obtained when both worlds are combined in such a way that the best of the two worlds can be applied at the most appropriate moment. This was the situation when the collaboration between Sygic and IT4I started and the use case was designed. In particular, main goal of the use case was to build a navigation infrastructure for making use of (near)real-time data coming from client devices (navigator appliance or smart-phones) and other data sources in order to create a so-called global view on the current traffic situation. Moreover, the challenge was to maintain the efficient operation of the system in response to the dynamic nature of the traffic and the amount of navigation requests.

The improvements driven by ANTAREX tools in the self-adaptive navigation system have been the following:

(i) We designed and developed an infrastructure to collect traffic data by heterogeneous data sources. The system has been designed to handle simultaneously the processing of large volume of data from various sources. The system has been deployed on an HPC infrastructure to handle a large number of incoming routing requests, while maintaining the Service Level Agreement (SLA). Moreover, the high computational power provided by the infrastructure has been used to elaborate more accurately the traffic data, thus providing the highest quality response to the navigation requests.

(ii) We designed the HPC-side of the navigation system to optimize traffic over a particular area. The scenario of a smart city where all cars are self-driving has been considered. This impacted not only the scalability analysis of the deployed service, but also the algorithm implementing the optimal routing that are thought to be not only for the single user, but more towards the entire set of cars.

(iii) We used the mARGOt autotuner to maintain the computation efficiency according to the dynamic nature of the traffic. A first tangible result as the fine-grain tuning of the application was obtained on the Probabilistic Time-Dependent Routing (PTDR) algorithm based on the Monte Carlo method [START_REF] Golasowski | Performance evaluation of probabilistic time-dependent travel time computation[END_REF]. Given that this algorithm is the most computationally intensive part of our routing pipeline, the mARGOt autotuner was used to determine the minimal amount of samples needed to obtain distribution of the travel time with sufficient precision. In the PTDR case, the number of samples used was reduced by more than 50%. An example of the adaptivity in the number of sample selection can be seen in Figure 4 where we plot the number of samples needed to satisfy the extra-functional requirement when the same requests have been made every 15 minutes during the entire week. The daily behavior of rush hours where more complicated is the arrival time prediction is identified by requests with higher samples.

(iv) We adopted a similar approach -considering a broader perspective -to optimize the amount of compute resources occupied by the system on an HPC cluster according to extra-functional properties constraints (e.g. cost, quality of results, response time). The implementation of such approach is currently in progress: We expect it to dynamically respond to variations in the routing requests load by adjusting the amount of occupied compute nodes by (v) We used the toolbox to integrate efficiently the new code features to existing codebase. This was successfully achieved by the integration of the LARA DSL in the built process to modify the source code in a transparent way, thus enabling a seamless integration of the necessary instrumentation needed by the tools. The DSL has been used to integrate the mARGOt autotuner, evaluate the code scalability analysis [START_REF] Golasowski R. Cmar | DSL and Autotuning Tools for Code Optimization on HPC Inspired by Navigation Use Case[END_REF] (see Figure 5) and to generate parts of the data access API which handles the road network graph stored in a HDF5 format [START_REF] Golasowski | Expressing and Applying C++ Code Transformations for the HDF5 API Through a DSL[END_REF].

CONCLUSIONS

Exascale HPC systems need the definition of new software stacks to fully exploit heterogeneity and meeting power efficiency requirements. ANTAREX provides an holistic programming framework capable to decouple functional and extra-functional aspects of the application for energy efficient supercomputing systems. In this paper, we have presented the ANTAREX toolbox and some preliminary results carried out on two different case studies.

Figure 1 :

 1 Figure 1: Overview of the ANTAREX flow

Figure 2

 2 Figure 2: a) Frequency for core0 and 18 for each socket of the 4 nodes of MARCONI supercomputer running the mini-app; b) Power consumption of each socket; c) and d) report the same results when fine-grain power management is used.

Figure 3 :

 3 Figure 3: LiGenDock accuracy-performance trade-off. The color of each point represents the value of one of the several knobs introduced in the application.

Figure 4 :

 4 Figure 4: Adaptive number of samples for the PTDR kernel of the navigation system. The same navigation requests is repeated over the week every 15 minutes. The original static version was designed for the worst case, 3K samples.

Figure 5 :

 5 Figure 5: Example of LARA-DSL adoption

This work has been funded by the EU grant 671623 FET-HPC-ANTAREX.

http://specs.fe.up.pt/tools/clava/doc/

https://github.com/skeru/libVersioningCompiler

https://github.com/fbeneventi/examon

https://gitlab.com/margot_project/core


Federica Caforio
email: federica.caforio@inria.fr

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

A conservative penalisation strategy for the semi-implicit time discretisation of the incompressible elastodynamics equation

Keywords: Elastodynamic equations, Time discretisation, Incompressibility

The principal aim of this work is to provide an adapted numerical scheme for the approximation of elastic wave propagation in incompressible solids. We rely on high-order conforming finite element with mass lumping for space discretisation and implicit/explicit, second-order, energy-preserving time discretisation. The time step restriction only depends on the shear wave velocity and at each time step a Poisson problem must be solved to account for the incompressibility constraint that is imposed by penalisation techniques.

Introduction

Since large applications in computational mechanics concern nearly and pure incompressible elasticity (living tissues, biomaterials), a great effort has been made in the last forty years to provide accurate Finite Element Method (FEM) formulations for the approximation of elasticity in incompressible solids. However, the majority of the works proposed in the literature only deals with static computations. The main contributions to date can be mainly divided into two categories: pure displacement methods and mixed methods. On the one hand, displacement-based FEM can provide accurate solutions of quasi or pure incompressible elasticity problems; nonetheless, the space resolution necessary to provide an accurate solution is far greater than the one required for a compressible material [1; 2]. Indeed, these methods can suffer of undesirable limitations, such as ill-conditioning of the stiffness matrix, spurious or incorrect pressures and numerical locking (severe stiffening near the incompressible limit) [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF], especially if low-order shape functions are adopted, due to the enforcement of the incompressibility constraint -i.e. the requirement that the displacement field is divergence-free. Locking is due to the fact that, in case of incompressible materials, volumetric strains approach zero, while the pressure field is of the order of the boundary traction, therefore it cannot be computed from strains, but it must be calculated directly from the equilibrium equations [2; 3]. Several methods have been proposed to improve accuracy of displacement-based methods. Among them, we cite the reduced/selective integration method proposed in [START_REF] Malkus | Mixed finite element methods. Reduced and selective integration techniques: a unification of concepts[END_REF], the B-bar method [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and nonlinear media[END_REF] and the F-bar method [START_REF] Neto | F-bar-based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF]. These methods circumvent volumetric locking by reducing the number of discrete incompressibility constraints enforced at the quadrature points [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]. We also cite a method adapted to nearly-incompressible materials with volumetric energy penalty function in the framework of the continuum-based absolute nodal coordinate formulation (ANCF) [START_REF] Orzechowski | Nearly incompressible nonlinear material models in the large deformation analysis of beams using ANCF[END_REF]. On the other hand, mixed finite element methods [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] have proven effective and even necessary to obtain accurate results in the resolution of incompressible fluid flows and incompressible elasticity. In these methods, the constrained problem is rewritten in form of an unconstrained saddle-point problem, due to the introduction of a second variable (namely, the pressure). However, not all mixed methods are stable. In fact, the convergence properties of this formulation are governed by stability considerations, involving ellipticity requirements and the famous Ladyzenskaya -Babǔska -Brezzi (LBB) inf-sup condition [10; 11]. For example, equal-order interpolation both for displacement and pressure field does not satisfy the LBB condition for classical mixed FEM [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. If this stability condition is not satisfied, severe unphysical oscillations in the pressure field can appear, named "spurious pressure modes". Stabilised methods have been proposed to overcome the limitations of classical mixed FE formulations in the field of incompressible fluid dynamics (see [12; 13; 14] and references therein for throughout reviews on the subject). Brezzi and collaborators [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF] proposed to extend the equation accounting for incompressibility in Stokes flows by adding a Laplacian of the pressure field. Other methods are based on the addition of artificial high-order differential terms to the discrete continuity equation, in order to let the formulation satisfy the LBB condition (for more details, see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]). We cite for example the Streamline Upwind Petrov-Galerkin method [17; 18; 19]. Similar stabilised methods have been recently extended to the context of linear elasticity [20; 21; 22; 23; 24]. All the aforementioned methods can be straightforwardly extended to the discretisation of dynamic equations using implicit time discretisation (e.g. Implicit Euler scheme or implicit Newmark schemes). However, at each time step, the resolution of the resulting linear system is required for the computation of the velocity field (or most often the displacement field in elasticity). We highlight that, in linear elasticity, this could be done in practice by performing a factorisation of the matrices to invert, since they are constant in time. However, for large scale problems, it is not possible to store any factorisation of the matrices and it is even difficult to store preconditioners. A popular approach to increase the efficiency of dynamic solvers was first proposed in computational fluid dynamics in the late 1960s [25; 26; 27; 28], and it is called fractional-step projection. This family of methods aims at accurately solving the equation governing viscous incompressible fluids by performing a time-discretisation in which viscosity and incompressibility are treated in two separate steps. In more detail, the first half-step (Burgers step) corresponds to an elliptic Boundary Value Problem (BVP) for an intermediate velocity, accounting for viscosity diffusion and advection. The second half-step (projection step) represents an inviscid problem where the end-of-step, divergence-free velocity is computed, along with pressure distribution. This step essentially consists in solving a Poisson problem. In this way, at each time step two decoupled elliptic equations are solved, and this is very advantageous for large scale simulations [29; 30; 31]. Less effort has been made to develop efficient methods for the treatment of the incompressibility constraint in elastodynamics. Since the underlying physics is wave propagation, fully explicit methods seem to be good candidates to obtain efficient schemes (this is possible for nearly-incompressible media). We cite in this context a recent method for quasi-incompressible elastodynamics [START_REF] Scovazzi | A simple, stable, and accurate linear tetrahedral finite element for transient, nearly, and fully incompressible solid dynamics: a dynamic variational multiscale approach[END_REF] with linear finite elements, extended in [START_REF] Rossi | Implicit finite incompressible elastodynamics with linear finite elements: A stabilized method in rate form[END_REF] to the purely incompressible case. Furthermore, a recent velocity/stress formulation for the simulation of linear elastodynamics in nearly-incompressible solids with weakly enforced boundary conditions was proposed by Scovazzi and collaborators [START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF]. However, the stability constraint (CFL condition) is drastically limited by the enforcement of incompressibility. In this regard, the efficiency of the fractional-step projection algorithm mentioned above, along with the similarities with viscous incompressible fluids, suggests the possibility to adopt the main ideas of this method to design an efficient scheme for incompressible elastodynamics. A first method that integrates a fractional time-step method for Lagrangian formulations of elastodynamics problems has been proposed by Lahiri and collaborators [START_REF] Lahiri | A variationally consistent fractional time-step integration method for incompressible and nearly incompressible Lagrangian dynamics[END_REF]. In this paper, the authors use variational integrators that take advantage of the Hamilton variational principle to construct a discrete approximation of the integral of the Lagrangian over a given interval, and they adopt linear finite element discretisation in space. We have developed a numerical scheme that carefully takes into account the intrinsic properties of the wave equation. Namely, we construct a conservative time discretisation, and treat implicitly only the terms corresponding to "informations" travelling at infinite velocity (i.e. the incompressibility constraint) by solving a Poisson problem. The fully discrete scheme that we propose to efficiently solve the incompressible elastodynamic problem is

       ỹn+1 α,h -2 ỹn α,h + ỹn-1 α,h ∆t 2 + A h ỹn α,h + B T h pn α,h = f n h , B h ỹn α,h = α ∆t 2 C h C T h pn α,h , (1)
where A h denotes a stiffness operator (independent of the compressibility parameters), B T h and B h correspond to discrete gradient and divergence operators, and C h C T h represents a discrete laplacian operator. The coefficient α corresponds to a penalisation parameter and (ỹ n α,h , pn α,h , f n h ) denote respectively the displacement field, the pressure and the source term. Therefore, if effective methods for explicit time-discretisation are used, our algorithm requires at each time step the resolution of a scalar Poisson problem (that can be done by several, efficient algorithms) and few matrix-vector multiplications for the explicit methods. Note that the stability condition for scheme (1) reads

∆t 2 4 A h -1 4α -1 4α , α > 1/4.
for some operator norm • that we define later. Contrarily to the standard results one could expect (i.e. ∆t 2 4 A h -1 , see [START_REF] Joly | The mathematical model for elastic wave propagation[END_REF]), the stability condition is slightly more constraining due to the factor (4α -1)/4α.

Our approach is a strategy to approximate elastic wave propagation in quasiincompressible media. In order to do that, we introduce a first good approximation: the pure incompressible formulation. Then, we construct a penalised formulation to approximate the pure incompressible problem. Our procedure is justified by arguments of stability, computational cost and numerical convergence of the resulting schemes.

Finally, the work we propose is closely related to the Selective Mass Scaling Method (SMS) developed in [START_REF] Ye | A selective mass scaling method for shear wave propagation analyses in nearly incompressible materials[END_REF] for wave propagation in quasi-incompressible materials. In this method, the mass matrix is modified: artificial mass is added, hence obtaining very good stability properties at the cost of computing a mass matrix that cannot be lumped anymore. The problem is then implicit for the displacement field, and thus it not possible to use fast solvers (like FFT solvers) for the inversion of a scalar Poisson problem.

The article is organised as follows. In Section 2 we provide two standard formulations of the continuous elastodynamic problem for quasi-incompressible and pure incompressible media along with a novel formulation for the treatment of incompressibility by penalisation. Furthermore, we derive the variational formulation associated with each problem. Section 3 deals with the abstract framework for space discretisation of the incompressible elastodynamics equation by Spectral Element Method. In Section 4 we provide the time discretisations for each formulation by finite difference. Then, a stability analysis based on energy considerations is performed for each scheme in 4.2 and pros and cons are discussed. Numerical results, including convergence curves to the solution of the incompressible elastodynamics equations for different choices of materials, are shown in Section 5. Finally, a three-dimensional illustration in a more realistic test case for elastography imaging is shown in Section 6.

Continuous framework

The equation of elastodynamics

Given a domain Ω ⊂ R d smooth enough, with d = 2 or d = 3, we introduce the following notations to define Hilbert spaces for the elastic displacements

H := v ∈ L 2 (Ω) d , X := H 1 0 (Ω) d , X = H -1 (Ω) d .
For the sake of simplicity, we consider homogenous Dirichlet condition on the boundary of the propagation domain. We also need to consider divergence-free displacements. Hence, we introduce the following subspace of

X V := {v ∈ X | div v = 0}.
Pressure is a variable of interest, and is sought in the spaces

L := L 2 0 (Ω) = q ∈ L 2 (Ω) | Ω q = 0 , M = {q ∈ H 1 (Ω) | q ∈ L}.
As usual, we identify L and H with their dual in what follows. Moreover, for the sake of conciseness, we define, given a function space A on Ω,

C k (A) := C k ([0, T ]; A), k ∈ {0, 1, 2, ..},
where T > 0 is a given final time of observation. Furthermore, we introduce Ω T := [0, T ] × Ω. Our aim is to analyse the propagation of elastic waves in heterogenous, anisotropic, incompressible solids and we consider as a reference model the following partial differential equation (PDE) system: For f given and sufficiently regular, find

y λ ∈ C 2 (H) ∩ C 1 (X) such that    ρ ∂ 2 t y λ -div Cε(y λ ) -λ ∇div y λ = f in Ω T , y λ (t = 0) = 0, ∂ t y λ (t = 0) = 0 in Ω, (QI)
with λ ∈ R + the bulk modulus, that is assumed to be large, due to nearlyincompressibility, ρ(x) the strictly positive density of the medium and C(x) the elasticity tensor which is symmetric, coercive and bounded, i.e. there exist two positive scalars c, C such that

c ε 2 C(x)ε : ε C ε 2 , ∀ ε. (2)

Non-dimensionalisation

Since we are going to consider a limit process where the bulk modulus tends to infinity, a first step is to non-dimensionalise our system of equations. To do so, we introduce a typical length scale L of the domain Ω, a typical observation time τ and a shear modulus µ, and we define a non-dimensionalised displacement as follows:

ŷλ (t, x) := L -1 y λ (τ t, L x), t ∈ [0, T ], x ∈ Ω,
where Ω is a rescaled domain and T = T /τ is of the order of unity. Note that t and x refer now to non-dimensionalised variables. We also introduce the non-dimensionalised quantities

ρ = L 2 ρ τ 2 µ , C = µ -1 C, λ = λ µ , f = L µ f .
Then, the equation of elastodynamics can be recast as

   ρ ∂ 2 t ŷλ -div Cε(ŷ λ) -λ ∇div ŷλ = f in [0, T ] × Ω, ŷλ (t = 0) = 0, ∂ t ŷλ (t = 0) = 0 in Ω.
We also introduce here the pressure field p λ associated with the displacement field y λ . Its non-dimensionalised counterpart is given by

pλ (t, x) := µ -1 p λ (τ t, L x), t ∈ [0, T ], x ∈ Ω.
For the sake of simplicity, we drop the notation • throughout the rest of the paper.

The mixed and penalised formulations

Existence and uniqueness results for problem (QI) are well-known (see [START_REF] Joly | The mathematical model for elastic wave propagation[END_REF]). An alternative, equivalent formulation to (QI) is obtained by introducing the scalar function

p λ := λ div y λ . The couple (y λ , p λ ) ∈ C 2 (H) ∩ C 1 (X) × C 0 (L) satisfies            ρ ∂ 2 t y λ -div (C(x)ε(y λ )) -∇ p λ = f in Ω T , div y λ = λ -1 p λ in Ω T , y λ (t = 0) = 0, ∂ t y λ (t = 0) = 0 in Ω. (QIM)
Since λ is large, it is natural to approximate the solution of (QIM) by the solution obtained at the limit as λ goes to infinity. More precisely, if one defines (y, p) and some corrector functions (y 1 , p 1 ) and (y 2 , p 2 ) such that

y λ = y + λ -1 y 1 + λ -2 y 2 + • • • , p λ = p + λ -1 p 1 + λ -2 p 2 + • • • , ( 3
)
then a standard asymptotic analysis procedure shows that (y, p) satisfies a pure incompressible problem. This formulation reads Find (y, p)

∈ C 2 (H) ∩ C 1 (V) × C 0 (L) such that            ρ ∂ 2 t y -div (C(x)ε(y)) -∇ p = f in Ω T , div y = 0 in Ω T , y(t = 0) = 0, ∂ t y(t = 0) = 0 in Ω. (IM)
Observe that the function p in (IM) acts as a Lagrange multiplier to enforce incompressibility. If we assume that (IM) is the standard problem to solve, then Eq. (QIM) can be seen as an approximate penalised formulation of (IM).

Existence, uniqueness and regularity results of problem (IM) can also be induced from (QI) by a limit process as λ → ∞.

We now introduce another formulation by penalisation of the problem (IM), inspired from existing formulations for the Stokes problem [START_REF] Boffi | Mixed finite element methods and applications[END_REF]. It reads

Find (ỹ α , pα ) ∈ C 2 (H) ∩ C 1 (X) × C 0 (M) such that            ρ ∂ 2 t ỹα -div (C(x)ε(ỹ α )) -∇ pα = f in Ω T , div ỹα = -α ∆p α in Ω T , ỹα (t = 0) = 0, ∂ t ỹα (t = 0) = 0 in Ω. (QIP)
First, note that (QIP) is different from (QIM) due to the introduction of the laplacian operator in the second equation. Furthermore, one can observe that (QIP) differs from the other formulations by several important but subtle aspects. First, the pressure pα is sought in a more regular space, namely C 0 (M), that is mandatory to give an appropriate meaning to the introduced Laplace operator. Second, system (QIP) is not a closed set of equations. Indeed, a boundary condition is required for the second equation for the pressure, that now has a trace (we recall that we use homogeneous Dirichlet conditions on the displacement). A standard choice is to use homogenous Neumann boundary conditions

∇ pα • n = 0 on [0, T ] × ∂Ω. ( 4
)
Then, using the same arguments as before, i.e., writing

ỹα = y + α ỹ1 + α 2 ỹ2 + • • • , pα = p + α p1 + α 2 p2 + • • • , (5)
one can see that (ỹ α , pα ) can be approximated by (y, p), solution of the pure incompressible mixed formulation (IM). Reciprocally, (ỹ α , pα ) represents another approximation of the pure incompressible mixed formulation (IM). Moreover, eliminating the unknown pα , it is possible to rewrite (QIP) as

ρ ∂ 2 t ỹα -div Cε(ỹ α ) -α -1 ∇(-∆) -1 divỹ α = f , ( 6
)
where ∆ -1 : L → M stands for the inverse Laplace operator equipped with a homogeneous Neumann boundary condition. It is possible to prove that the operator -div Cε(•) -α -1 ∇(-∆) -1 div(•) defines a self-adjoint coercive bilinear form. Consequently, existence and uniqueness of the solution can be retrieved from ( 6) by standard theory.

Remark 1 Note that the choice ( 4) is arbitrary. In the context of the Stokes equations it was observed (see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] and reference therein) that it results in a boundary layer that deteriorates the approximation of the gradient of the solution. Correcting terms can be introduced in specific cases (see again [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]) but their analysis is more difficult.

Weak formulation of the continuous PDE

Let us first consider the weak formulation associated with (IM). Given f sufficiently regular, it reads:

Find (y, p) ∈ C 2 (H) ∩ C 1 (X) × C 0 (L) such that for all w ∈ X, q ∈ L          m(∂ 2 t y, w) + a(y, w) + b(p, w) = (f , w), b(q, y) = 0, y(t = 0) = 0, ∂ t y(t = 0) = 0, (7)
where we have defined the bilinear forms m :

X × X → R, a : X × X → R and b : L × X → R such that m(y, w) := Ω ρ y • w dΩ, a(y, w) := Ω C ε(y) : ε(w)dΩ, b(p, w) := Ω p div w dΩ.
Note that the bilinear form m(•, •) is symmetric and positive. Furthermore, due to Eq. ( 2) and the Korn inequality, the bilinear form a(•, •) is symmetric and coercive in X. We can write [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF] as a set of equations written in X × L for all t ∈ [0, T ]. To do so, we introduce the linear continuous operators M :

X → X , A : X → X , such that ∀ (y, w) ∈ X × X M y, w = m(y, w), Ay, w = a(y, w), (8)
and we define the divergence operator B : X → L and its transpose B

T : L → X such that, ∀ (y, q) ∈ X × L, (By, q) L = b(q, y) = B T q, y , (9)
where (•, •) L denotes the standard scalar product of L 2 (Ω). Note that the operator M can be continuously extended to an operator from H to H and, if ρ ≡ 1, then it is the identity operator. Finally, the variational formulation can be equivalently rewritten Find (y, p)

∈ C 2 (H) ∩ C 1 (X) × C 0 (L) such that for all t ∈ [0, T ]    M ∂ 2 t y + A y + B T p = f in X , B y = 0 in L, (IMv)
and

y(t = 0) = 0, ∂ t y(t = 0) = 0.
Analogously, the variational formulation associated with system (QIM) reads, equivalently,

Find (y λ , p λ ) ∈ C 2 (H) ∩ C 1 (X) × C 0 (L) such that for all t ∈ [0, T ]    M ∂ 2 t y λ + A y λ + B T p λ = f in X , B y λ = λ -1 I p λ in L, (QIMv)
and

y λ (t = 0) = 0, ∂ t y λ (t = 0) = 0,
where I is the identity operator from L to L. Concerning the problem (QIP) with the boundary condition ( 4), we propose a variational formulation that reads

Find (ỹ α , pα ) ∈ C 2 (H) ∩ C 1 (X) × C 0 (M) such that for all (w, q) ∈ X × M      m(∂ 2 t ỹα , w) + a(ỹ α , w) + b(p α , w) = (f , w), b(q, ỹα ) = α (∇p α , ∇q) L , (10)
and ỹα (t = 0) = 0, ∂ t ỹα (t = 0) = 0.

It is not straightforward to write equations in dual spaces from the variational formulation [START_REF] Brezzi | A discourse on the stability conditions for mixed finite element formulations[END_REF]. This is due to the fact that we have changed the functional space in which the pressure is sought. To do so, we introduce the divergence operator C : H → M and its transpose C T : M → H defined by

c(q, w) = - Ω ∇ q • w dΩ Cw, q = (C T q, w) H = c(q, w), (11)
where (•, •) H denotes the standard scalar product on L 2 (Ω) d . Then, by identification of the operator B as an operator from X to M (instead of an operator from X to L), one can show that ( 10) is equivalent to

Find (ỹ α , pα ) ∈ C 2 (H) ∩ C 1 (X) × C 0 (M) such that, for all t ∈ [0, T ],    M ∂ 2 t ỹα + A ỹα + B T pα = f in X , B ỹα = αC C T pα in M , (QIPv)
and ỹα (t = 0) = 0, ∂ t ỹα (t = 0) = 0.

Note that the operator C corresponds to an extension of the operator B. Inversely, the gradient operator B T represents the extension of C T in a larger space. This distinction in the notation is not relevant in the continuous framework, but it will be fundamental at the discrete level.

Space discretisation

Let us now consider a regular finite-dimensional space X h ⊂ X for the discretisation of the displacement field and M h ⊂ M for the discretisation of the pressure field. These spaces are obtained by finite element approximation of X and M, respectively. Furthermore, inspired by [START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF], in order to discretise accordingly the variational formulation ( 10) that we propose, we introduce a third finite-dimensional space denoted Y h ⊂ H that should satisfy, for the sake of simplicity,

X h ⊂ Y h and ∇M h ⊂ Y h .
We define the embedding i h :

X h → Y h and its transpose i T h : Y h → X h such that for all (y h , w h ) ∈ X h × Y h (y h , i T h w h ) H := (i h y h , w h ) Y h := (y h , w h ) Y h , (12)
where (•, •) Y h stands for the approximation of the scalar product in H that is defined using quadrature formulae (it is a symmetric coercive and continuous bilinear form in Y h for the norm in H). In a more general way, we assume that quadrature formulae do not corrupt symmetry and positivity properties of the bilinear forms. We introduce the discrete divergence operator C h : Y h → M h and the discrete gradient operator

C T h : M h → Y h such that for all (p h , w h ) ∈ M h × Y h (p h , C h w h ) L := (C T h p h , w h ) Y h := -(∇p h , w h ) Y h .
Observe that C T h corresponds to the operator C T ≡ -∇ applied to functions in M h , considering the resulting functions in the larger space Y h . Then, we define another discrete divergence operator B h : X h → M h and another discrete gradient operator B T h : M h → X h as

B h := C h i h and B T h := i T h C T h . ( 13
)
One can observe that, for all (p h , w h ) ∈ M h × X h , we have

(B h w h , p h ) L = (C h i h w h , p h ) L = (i h w h , C T h p h ) Y h = (w h , i T h C T h p h ) H = (w h , B T h p h ) H . ( 14
)
Note also that the following commutative diagrams (taken from [START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF]) hold

X h M h Y h B h i h C h X h M h Y h B T h C T h i T h
We introduce the linear continuous operators

M h : X h → X h and A h : X h → X h such that for all (v h , w h ) ∈ X h × X h (M h v h , w h ) H = m h (v h , w h ), (A h v h , w h ) H = a h (v h , w h ). ( 15
)
where the subscript h in the notation of the bilinear forms stands for the use of quadrature rule in the computation of integrals. The finite element approximation of (IMv) reads

Find (y h , p h ) ∈ C 2 (X h ) × C 0 (M h ) such that for all t ∈ [0, T ]    M h ∂ 2 t y h + A h y h + B T h p h = f h in X h , B h y h = 0 in M h , (IMh)
and

y h (t = 0) = 0, ∂ t y h (t = 0) = 0,
where f h (t) ∈ X h denotes some approximation of f . Analogously to the pure incompressible mixed formulation, we can retrieve the space discretisation associated with (QIMv). It reads

Find (y λ,h , p λ,h ) ∈ C 2 (X h ) × C 0 (M h ) such that for all t ∈ [0, T ]    M h ∂ 2 t y λ,h + A h y λ,h + B T h p λ,h = f h in X h , B h y λ,h = λ -1 I h p λ,h in M h , (QIMh)
and y λ,h (t = 0) = 0, ∂ t y λ,h (t = 0) = 0.

The operator

I h : M h → M h is such that for all (p h , q h ) ∈ M h (I h p h , q h ) L = (p h , q h ) L h ,
where (•, •) L h stands for the approximation of the scalar product in L by quadrature formulae. Finally, we are able to give the space discretisation associated with the novel formulation we propose (QIPv). It reads

Find (ỹ α,h , pα,h ) ∈ C 2 (X h ) × C 0 (M h ) such that for all t ∈ [0, T ]    M h ∂ 2 t ỹα,h + A h ỹα,h + B T h pα,h = f h in X h , B h ỹα,h = α C h C T h pα,h in M h , (QIPh)
and ỹα,h (t = 0) = 0, ∂ t ỹα,h (t = 0) = 0.

Note that the operator C h C T h corresponds to a discrete laplacian operator on

L h . Indeed, for all (p h , q h ) ∈ M h × M h (C h C T h p h , q h ) L = (C T h p h , C T h q h ) Y h = (∇p h , ∇q h ) Y h .
Therefore, C h C T h : M h → M h is invertible for any reasonable choice of finite element spaces and quadrature rule. This is obviously also true for I h and M h . Let us insist on the importance of the introduction of the space Y h . First, it is related to the definition of the quadrature formulae in H in the definition of [START_REF] Franca | Stabilized finite element methods: II. The incompressible Navier-Stokes equations[END_REF]. Second, even if exact integration is performed, the introduction of the space Y h allows us to take into account the fact that, in general, the gradient of functions in M h does not belong to X h . Indeed, if this was the case, then C h and B h would be the same operator and the penalisation strategy would be useless in terms of computational efficiency (see the discussion of Section 4.1).

In the numerical results we present in this work we use high-order Spectral Finite Element, as in [START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF] and [START_REF] Cohen | Mixed spectral finite elements for the linear elasticity system in unbounded domains[END_REF]. Since we consider a simple geometry for our purposes, we construct a quasi-uniform triangulation of Ω composed of quadrangles or hexahedra

Ω = N i=1 K i , K i ∩ K j = ∅, F i ( K) = K i ,
where K is the unit square or the unit cube and ∀ i ∈ {1, 2, .., N }, F i denotes the invertible transformation of the reference element K to the deformed element K i . Then, we define

X h = φ ∈ C 0 (Ω) d | φ| Ki • F i ∈ Q r ( K) d ,
where Q r is the set of polynomials with degree r 1 in each variable of space. To obtain mass-lumping (meaning that M h can be inverted trivially) one must choose the quadrature points in the computation of m h (•, •) at the same location as the interpolation points (see [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]). Sufficient accuracy is obtained if the interpolation and quadrature points correspond in the reference elements to the Gauss-Lobatto points of same order (r in our case, that gives (r + 1) d quadrature/interpolation points in the reference elements). For the computation of a h (•, •) we use the same quadrature rule as for the computation of m h (•, •) (which gives sufficient accuracy on non-distorted mesh, see [START_REF] Durufle | Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain[END_REF]). Then, we choose

M h = ϕ ∈ C 0 (Ω) | ϕ| Ki • F i ∈ Q r-1 ( K) .
Note that the major difference here is that we choose a lower-order finite element space, but M h is still constructed using continuous finite element.

It is mentioned in [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] that this choice is compatible with X h in the sense that a discrete inf-sup condition is satisfied (the importance of this condition is further detailed below). In addition, we assume that (•, •) L h is computed using Gauss-Lobatto quadrature points of order r -1. Hence, mass-lumping is achieved and the operator I h is easily invertible. Finally, we define

Y h = X h + ψ ∈ L 2 (Ω) d | ψ| Ki • F i ∈ Q r-1 ( K) d
and (•, •) Y h is computed using the Gauss-Lobatto quadrature points of order r -1. Note that we have the inclusion ∇M h ⊂ Y h only if all the F i are affine. We believe that this is only a technical limitation. However, numerical results are provided only in that case.

Time discretisation

This section deals with the time discretisation of the semi-discrete formulations obtained by FE approximation in space. We consider only finite difference schemes that are centred, in order to preserve energy conservation at the discrete level. In what follows, the fully discrete schemes for the standard formulations (IM) and (QIM) are provided. Moreover, we propose the fully discrete scheme for the novel formulation (QIP).

Fully discrete schemes

Let us consider a time interval [0, T ], with T > 0, and define the partition t n = n ∆t, with n ∈ {0, 1, .., N }, and ∆ t = T /N . The fully discrete scheme corresponding to (IM) for n ∈ {0, 1, .., N } is constructed based on a simple second-order finite difference scheme, namely a leapfrog scheme. We shall consider two sequences of approximate displacement fields {y n h ∈ X h } and pressures

{p n h ∈ L h } such that (y 0 h , y 1 h ) is given and for n ∈ {1, , N }        M h y n+1 h -2 y n h + y n-1 h ∆t 2 + A h y n h + B T h p n h = f n h in X h , B h y n h = 0 in M h , (IMnh) with f n h = f h (t n
). The implementation of scheme (IMnh) can be done using Schur complement techniques. It reads as follows: given (y

n-1 h , y n h ), first compute p n h from (B h M -1 h B T h ) p n h = B h M -1 h (f n h -A h y n h ); (16)
then, y n+1 h is given by

y n+1 h = ∆t 2 M -1 h f n h + 2 y n h -y n-1 h -∆t 2 M -1 h A h y n h -∆t 2 M -1 h B T h p n h .
One can see that the pressure field is an intermediate unknown that acts as a Lagrange multiplier to enforce the constraint B h y n+1 h = 0. Note that we use the notation p n h by analogy with the quasi-incompressible schemes that we present in what follows. Moreover, the system ( 16) is well-posed if B T h is injective and B h surjective. This corresponds to verify that the LBB condition is satisfied, i.e. there exists a constant c > 0 such that inf

q h ∈M h sup w h ∈X h (B h w h , q h ) w h X q h L c, (17)
where, if c does not depend on the discretisation parameter, one recovers an optimal convergence behaviour. Note that for the incompressible, linear Stokes problem it is standard to show [42; 43; 44; 45] that standard second-order discrete space-time discretisation can be achieved as soon as the LBB condition holds. We assume that similar results hold for the elastodynamic problem (IMh). Finally, note that Eq. ( 16) can be solved at each time step by iterative algorithms and, since the underlying problem is symmetric, the Conjugate Gradient method is a good candidate. However, it is important to highlight the fact that the operator

B h M -1 h B T h = C h i h M -1 h i T h C T h is much more com- plicated to invert than C h C T
h due to the presence of the term i h M -1 h i T h that has the effect to widen the bandwidth of the corresponding finite element matrix. Our aim is to specifically tackle this problem by avoiding -at each time step -the inversion of B h M -1 h B T h , and inverting instead C h C T h , for which we have efficient solvers.

Remark 2 If non-zero initial displacement y λ (t = 0) = y 0 and/or initial velocity ∂ t y λ (t = 0) = v 0 are considered, or if the source term f h is not 0 at time t = 0, then, to preserve the expected second-order consistency, the computation of the first two iterates is performed as follows:

y 0 h = y h,0 , y 1 h = y h,0 + ∆t v h,0 + ∆t 2 2 M h f 0 h -A h y h,0 -B T h p 0 h ) (18)
where p 0 h is computed from [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF] with n = 0, and (y h,0 , v h,0 ) belong to X h × X h and correspond to an approximation of (y 0 , v 0 ) .

In order to write the fully discrete scheme corresponding to (QIM) for n ∈ {1, .., N }, we define two sequences of approximate displacement fields {y n λ,h ∈ X h } and pressures

{p n λ,h ∈ M h } such that (y 0 λ,h , y 1 λ,h ) is given and for n ∈ {1, .., N }        M h y n+1 λ,h -2 y n λ,h + y n-1 λ,h ∆t 2 + A h y n λ,h + B T h p n λ,h = f n h in X h , B h y n λ,h = λ -1 I h p n λ,h in M h . (QIMnh)
Observe that the scheme (QIMnh) is fully explicit, due to the use of the masslumping technique (we recall that M h and I h are easily invertible). However, we show in the following section that, due to stability considerations, the maximum time step allowed is considerably reduced by the fact that the pressure term is treated explicitly. Observe that the first two iterates can be computed using [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF] with p 0 h replaced by p 0 λ,h = λ I -1 h B h y 0 .

Finally, we provide the fully discrete scheme corresponding to (QIP) for n ∈ {1, .., N }. We define two sequences of approximate displacement fields {ỹ n α,h ∈ X h } and approximate pressures

{p n α,h ∈ M h } such that (ỹ 0 α,h , ỹ1 α,h ) is given and for n ∈ {1, , N }        M h ỹn+1 α,h -2 ỹn α,h + ỹn-1 α,h ∆t 2 + A h ỹn α,h + B T h pn α,h = f n h in X h , B h ỹn α,h = α ∆t 2 C h C T h pn α,h in M h . (QIPnh)
Note that here, for consistency reasons, we have rescaled the penalisation parameter by ∆t 2 and assume α is independent of ∆t. This choice should guarantee the second-order consistency in time that is expected from the leapfrog time discretisation. Note that we can directly rewrite the second equation in (QIPnh) as

pn α,h = (C h C T h ) -1 α ∆t 2 B h ỹn α,h . (19)
Consequently, this step is equivalent to solving a discrete Poisson problem for the pressure at each time step, with homogeneous Neumann boundary conditions on the boundary ∂Ω. One of the main advantages of this formulation is that the Poisson problem is very standard and its solution can be retrieved by fast solvers. In addition, observe that the first two iterates can be computed using [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF] with p 0 h computed with Eq. ( 19) for n = 0 and ỹ0 α,h ≡ y h,0 . Finally, observe that by injecting [START_REF] Franca | A new family of stable elements for nearly incompressible elasticity based on a mixed Petrov-Galerkin finite element formulation[END_REF] in (QIPnh), we retrieve

M h ỹn+1 α,h -2 ỹn α,h + ỹn-1 α,h ∆t 2 + A h ỹn α,h + B T h (C h C T h ) -1 B h α ∆t 2 ỹn α,h = f n h . ( 20
)
Remark 3 The SMS method, first introduced in [START_REF] Olovsson | Selective mass scaling for explicit finite element analyses[END_REF], consists in adding an inertial term. In the specific case of quasi-incompressible materials, it is suggested in [START_REF] Ye | A selective mass scaling method for shear wave propagation analyses in nearly incompressible materials[END_REF] to add the volumetric contribution of the stiffness operator. In our framework, starting from (QIMnh), after eliminating the pressure term p n λ,h , by SMS strategy one could obtain the scheme

M h + βλ B T h I -1 h B h y n+1 β,h -2 y n β,h + y n-1 β,h ∆t 2 + A h y n β,h + λ B T h I -1 h B h y n β,h = f n h , ( 21
)
where β is some well-chosen parameter. Note that the operators I h and B h embed the definition of adequate quadrature formulae to avoid numerical locking effects (typically, reduced integration must be used). Furthermore, observe that

β y n+1 β,h -2 y n β,h + y n-1 β,h ∆t 2 = -y n β,h + β ∆t 2 y n+1 β,h + 1 -2 β ∆t 2 y n β,h + β ∆t 2 y n-1 β,h ,
and therefore the scheme ( 21) can be rewritten, with θ := β ∆t 2 , as

M h y n+1 β,h -2 y n β,h + y n-1 β,h ∆t 2 + A h y n β,h + λ B T h I -1 h B h θy n+1 β,h + (1 -2 θ)y n β,h + θy n-1 β,h = f n h .
The scheme above is an explicit-implicit hybrid scheme (see [START_REF] Rylander | Stability of explicit-implicit hybrid timestepping schemes for Maxwell's equations[END_REF]) which has a stability condition depending only on (M h , A h ) (i.e. the operators related to the shear wave propagation only) as soon as θ 1/4. Although very good stability properties can be proved, at each time step one has to find the solution of a linear system for the displacement field and, thus, it not possible to use fast solvers for the inversion of a scalar Poisson problem. This last point motivated the introduction of the scheme (QIPnh).

Stability analysis

The aim of this section is to find uniform estimates of the discrete energy of the different schemes, i.e.

sup n∈{0,1,..,N } E n+ 1 2 h C, ( 22
)
where the constant C depends on the final time step T = N ∆t and on the data of the continuous problem, but is independent of ∆t and h. However, if explicit schemes are employed for time discretisation, the time step is limited by a stability condition depending on h. We refer to [START_REF] Joly | The mathematical model for elastic wave propagation[END_REF] for further reading.

For the sake of simplicity, we consider a time t n such that the source term has vanished (i.e. f h (t) = 0, ∀t t n ). Then, the energy of the continuous problem is constant in time. In order to retrieve a discrete counterpart of this energy, we consider for every formulation, as a test function, the centred discrete approximation of the time derivative of the displacement at time t n .

Stability of scheme (IMnh).

Let us first consider formulation (IMnh). By scalar product in the first equation in (IMnh) with

v h = y n+1 h -y n-1 h 2 ∆t ( 23
)
as well as by discrete differentiation of the second equation in (IMnh) and by scalar product with p n h we obtain

         M h y n+1 h -2 y n h + y n-1 h ∆t 2 + A h y n h + B T h p n h , y n+1 h -y n-1 h 2 ∆t H = 0, B h y n+1 h -y n-1 h 2 ∆t , p n h L = 0. ( 24
)
Then, due to symmetry of B h and Eq. ( 14), we can simplify the first equation in [START_REF] Zienkiewicz | The finite element method: solid mechanics[END_REF]. We have

M h y n+1 h -2 y n h + y n-1 h ∆t 2 + A h y n h , y n+1 h -y n-1 h 2 ∆t H = 0. ( 25
)
We define the discrete energy at time n + 1 2 as

E n+ 1 2 h = E n+ 1 2 k - ∆t 2 4 E n+ 1 2 kp + E n+ 1 2 p , ( 26
)
where the kinetic energy reads

E n+ 1 2 k := 1 2 M h y n+1 h -y n h ∆ t , y n+1 h -y n h ∆ t H ,
the potential energy reads

E n+ 1 2 p := 1 2 A h y n+1 h + y n h 2 , y n+1 h + y n h 2 H ,
and the mixed energy term is defined as

E n+ 1 2 kp := 1 2 A h y n+1 h -y n h ∆ t , y n+1 h -y n h ∆ t H .
Then, after some computations, using the symmetry properties of the operators M h and A h , we obtain from ( 25) the discrete conservation property

E n+ 1 2 h -E n-1 2 h ∆t = 0. ( 27
)
It now remains to prove that E

n+ 1 2 h is positive.
Proposition 1 A sufficient condition for the stability of scheme (IMnh) is

∆t 2 4 M -1 h A h -1 , ( 28
)
with

M -1 h A h = sup 0 = y h ∈X h a h (y h , y h ) m h (y h , y h ) . ( 29
)
Proof The proof is very standard (see [START_REF] Joly | The mathematical model for elastic wave propagation[END_REF]). Consequently, we only prove the positivity of the energy. We provide here some details for the sake of completeness. Let us first consider the definition of E

n+ 1 2 h . Since E n+ 1 2 s
is positive by definition, we can easily retrieve the lower bound for the energy

E n+ 1 2 h 1 2 M h - ∆t 2 4 A h y n+1 h -y n h ∆ t , y n+1 h -y n h ∆ t H . ( 30
)
Hence, the energy is positive if

m h (y h , y h ) - ∆t 2 4 a h (y h , y h ) 0, ∀ y h ∈ X h . ( 31
)
Finally, Eq. ( 31) can be rewritten as

∆t 2 4 sup 0 = y h ∈X h a h (y h , y h ) m h (y h , y h ) 1, (32)
concluding the proof.

Remark 4 Note that Proposition 1 introduces an abstract CFL condition. Moreover, we expect the following estimation

M -1 h A h ∼ c 2 s h 2 , ( 33
)
with c s a positive constant, depending on the elasticity tensor driving the shear wave propagation. Consequently, we obtain the sufficient stability condition

∆t h c s . ( 34
)
Therefore, the time step is not affected by the pressure wave propagation, that is travelling at an "infinite" velocity at the incompressible limit.

Stability of the scheme (QIMnh).

By similar reasoning, we can retrieve en energy estimation for the formulation (QIMnh). We retrieve, after some computations,

M h y n+1 λ,h -2y n λ,h + y n-1 λ,h ∆t 2 +A h y n λ,h +λB T h B h y n λ,h , y n+1 λ,h -y n-1 λ,h 2 ∆t H = 0. ( 35
)
We introduce A h,λ := A h + λB T h B h . Then, if we define the discrete energy E n+ 1 2 h as in Eq. ( 26), with A h,λ instead of A h , we can assert the following proposition.

Proposition 2 A sufficient condition for the stability of scheme (QIMnh) is

∆t 2 4 M -1 h A h,λ -1 , ( 36
)
with

M -1 h A h,λ = sup 0 = y h ∈X h a h (y h , y h ) + λ B h y h 2 L m h (y h , y h ) . ( 37
)
Remark 5 Note that, by definition,

M -1 h A h,λ M -1 h A h + λ M -1 h B T h B h , with M -1 h B T h B h = sup 0 = y h ∈X h B h y h 2 L m h (y h , y h ) .
By similar reasoning to Eq. ( 33), we can now introduce a constant c p , related to the maximum generalised eigenvalue of the operator (

B T h B h , M h ) such that λ M -1 h B T h B h ∼ c 2 p h 2 . ( 38
)
Therefore, we can assert

M -1 h A h,λ ∼ c 2 s + c 2 p h 2 . ( 39
)
Consequently, the stability condition (36) imposes a significant restriction on the time step.

∆t 2 h 2 c 2 s + c 2 p
Note that, because of the non-dimensionalisation, we expect c s to be close to the unity, whereas, c p is given by the velocity ratio between pressure and shear waves. For soft tissues this ratio is around 10 3 . This makes in practice the scheme (QIMnh) not efficient and justifies our need to formulate more adequate methods for the limit -incompressible -problem.

Stability of the scheme (QIPnh).

Finally, let us analyse the stability estimates related to the novel formulation (QIPnh). By analogous reasoning, we get from ( 20)

M h ỹn+1 α,h -2ỹ n α,h + ỹn-1 α,h ∆t 2 + A h + 1 α∆t 2 Q h ỹn α,h , ỹn+1 α,h -ỹn-1 α,h 2 ∆t = 0, (40)
where we have defined

Q h := B T h (C h C T h ) -1 B h . Then, the discrete energy at time n + 1 2 reads E n+ 1 2 h = E n+ 1 2 k - ∆t 2 4 E n+ 1 2 kp + E n+ 1 2 p , (41)
with kinetic energy

E n+ 1 2 k := M h ỹn+1 α,h -ỹn α,h ∆ t , ỹn+1 α,h -ỹn α,h ∆ t H , potential energy E n+ 1 2 p := A h + 1 α∆t 2 Q h , ỹn+1 α,h + ỹn α,h 2 , ỹn+1 α,h + ỹn α,h 2 H ,
and mixed energy term

E n+ 1 2 kp := A h + 1 α∆t 2 Q h ỹn+1 α,h -ỹn α,h ∆ t , ỹn+1 α,h -ỹn α,h ∆ t H .
Before providing a stability estimate for (40), we introduce the following lemma.

Lemma 1 Let P

h := C T h (C h C T h ) -1 C h : Y h → Y h . Then P h is a projection and P h L(Y h ) 1.
Proof The proof consists in demonstrating that P 2 h = P h . That follows easily by the definition of P h .

We are now able to assert a stability condition for scheme (QIPnh). Proposition 3 A sufficient condition for the stability of scheme (QIPnh) is ∆t 2 4 sup

0 = y h ∈X h 4 α a h (y h , y h ) 4 α m h (y h , y h ) -y h 2 Y h -1 ( 42
)
and

α > 1 4 sup 0 = y h ∈X h y h 2 Y h m h (y h , y h ) . ( 43
)
Proof Again we prove the positivity of the energy E

n+ 1 2 h
. By definition of E

n+ 1 2 h , since E n+ 1 2 p
is positive, we can obtain the estimation

E n+ 1 2 h M h - ∆t 2 4 A h - 1 4α Q h ỹn+1 α,h -ỹn α,h ∆ t , ỹn+1 α,h -ỹn α,h ∆ t H . ( 44
)
Hence, we need to satisfy

m h (y h , y h ) - ∆t 2 4 a h (y h , y h ) - 1 4α (Q h y h , y h ) H 0, ∀ y h ∈ X h . ( 45
)
However, by definition of P h and Lemma 1, we find Q h = i T h P h i h and, since we easily deduce the inequality i h y h Y h y h Y h from Eq. ( 12), we have

(Q h y h , y h ) H = (i T h P h i h y h , y h ) H = (P h i h y h , i h y h ) Y h i h y h 2 Y h y h 2 Y h .
Finally, if (42) holds, then Eq. ( 45) can be rewritten as

∆t 2 4 sup 0 = y h ∈X h a h (y h , y h ) m h (y h , y h ) -y h 2 Y h /(4α) 1, (46)
concluding the proof.

Corollary 1 If exact integration is used and the density ρ is constant, then a sufficient condition for the stability of scheme (QIPnh) is

∆t 2 4 ρ A h -1 4 α ρ -1 4 α ρ , α > 1 4 ρ , ( 47
)
with

A h = sup y h ∈X h y h =0 a(y h , y h ) y h 2 H . Proof If exact integration is used and ρ = const, then m h (y h , y h ) = ρ y h 2 H
and y h Y h = y h H . Consequently, we can simplify Eq. ( 42) to obtain the result of the corollary.

Remark 6

We observe that the CFL condition ( 47) is well-defined for α > 1/(4 ρ). Moreover, it is slightly worse than condition [START_REF] Temam | Navier-Stokes equations: theory and numerical analysis[END_REF], since

4 α ρ -1 4 α ρ 1
in the allowed range for α. Nevertheless, this condition depends only on the space discretisation and on the tensor C and the density ρ. Therefore, it is still very advantageous with respect to condition [START_REF] Joly | The mathematical model for elastic wave propagation[END_REF].

As a final comment, note that if the scheme (QIPnh) is stable, i.e. ( 22) holds and ∆t is sufficiently small, then

sup n∈{0,1,..,N } E n+ 1 2 p C, (48)
with C independent of ∆t. Therefore, denoting ỹn+

1 2 α,h = (ỹ n+1 α,h + ỹn α,h )/2, we have sup n∈{0,1,..,N } Q h ỹn+ 1 2 α,h , ỹn+ 1 2 α,h H α C ∆t 2 .
As a consequence, there exists another constant C independent of ∆t such that, by definition of Q h , we have

sup n∈{0,1,..,N } (C h C T h ) -1 2 B h ỹn+ 1 2 α,h H C ∆t.
This shows that B h ỹn+ 1 2 α,h goes to 0 with ∆t, but for a weak norm that involves the inverse Laplace operator.

Two-dimensional numerical convergence results

In order to perform the numerical validation of the properties of scheme (QIPnh), we consider as a model problem the elastic wave propagation in a 2D medium, and we take into account different constitutive laws. In the interest of clarity, we provide the physical and numerical parameters used for the simulations with their original units of measure. For the sake of simplicity, we assume constant density ρ = 1050 kg m -1 and homogeneous Dirichlet boundary conditions on all the boundaries of the domain. Space discretisation is performed by high-order Spectral Finite Elements (of order 7 for the displacement, 6 for the pressure). The computational grid is a 1 m 2 square composed of N uniform elements of size h = 1/N in each direction. Concerning time discretisation, we adopt the discretisation introduced in scheme (QIPnh) with penalisation coefficient α = 1/3ρ and we choose the time step as

∆t 2 = (1 -ε) M -1 h A h -1 , ( 49
)
where we set ε = 0.2 to account for the fact that the expression above is an approximation of the CFL condition given by Eq. ( 42) when considering that ρ y h 2 Y h = m h (y h , y h ) and accounting for the approximation of the norm in (49) by a power iteration algorithm. Note that the value of this parameter is not tight: ε = 0 also gives stable results. For each iteration of this scheme, the scalar pressure field [START_REF] Franca | A new family of stable elements for nearly incompressible elasticity based on a mixed Petrov-Galerkin finite element formulation[END_REF] is computed by means of an in-house fast solver for the Poisson problem, based on Higher-Order Fourier Transform. Note that even when the medium is heterogeneous, the problem for the pressure wave remains homogeneous, due to incompressibility.

Henceforth, we present a space/time convergence analysis. To do so, we test scheme (QIPnh) against the pure incompressible scheme (IMnh) for different values of N in [START_REF] Brezzi | A discourse on the stability conditions for mixed finite element formulations[END_REF][START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]. For each time step of each simulation, the pressure field is evaluated in Eq. ( 16) by means of the Conjugate Gradient method, with euclidian norm of the residue lower than 1e -14 and maximum number of iterations N iter = 2000. The source term considered is a standard Gaussian profile in space multiplied by a profile in time corresponding to the first derivative of the Gaussian function. In greater detail, it reads

f (x, y, t) := - 2 t σ 2 t exp - (t-t pulse ) 2 σ 2 t •100 exp - (x-x F ) 2 σ 2 s + (y-y F ) 2 σ 2 s , ( 50
)
with centre (x F , y F ) = (0.5, 0.5) m, covariance σ 2 s = 0.0005, σ 2 t = 0.0005 and mean t pulse = 0.6 s. We consider several constitutive laws, and we present four convergence curves for each example, that illustrate the error in L 2 (Ω T ), C 0 (L 2 (Ω)), L 2 (H 1 (Ω)) and C 0 (H 1 (Ω)) norms, respectively.

Homogeneous isotropic material

As a first example, we study a homogenous isotropic medium. A standard constitutive law for this type of medium is

C ε(y) := µ ε(y). ( 51
)
with shear modulus µ = 80 kPa. Note that from Eq. ( 51) we retrieve in (QI) the standard elastodynamic problem for isotropic law. The convergence curves in Figure 1 and Figure 2 confirm that second-order convergence is preserved in L 2 (Ω T ). Note that it is slightly degraded in and L ∞ in time and L 2 norm in space. Furthermore, the error does not vary (or it varies marginally in L ∞ (L 2 (Ω))) if we consider a larger t end such that the wave has reached the boundaries of the computational domain. If we consider H 1 norm in space, the error is slightly degraded with respect to the previous norms. Moreover, when we consider a larger time interval (so that we take into account the effects at the boundaries), the error increases. Figures 3 illustrates three snapshots at t = 1.0 s, t = 2. s and t = 2.65 s, of the displacement. The mesh is composed by N = 40 elements in each direction. Furthermore, Figure 4 depicts the time evolution of the displacement in three locations, indicated in Figure 3. Figures 5 and6 are related to the acceleration. We do not plot the results concerning the velocity, since they do not provide relevant information. Note that the time evolution profiles in x and y direction of the location denoted P 1 correspond to the profiles in y and x direction of the location denoted P 2 , that is the symmetric to P 1 , due to isotropy. Concerning the pressure, Figure 7 indicates that this term corresponds to a correction mostly at the boundaries of the medium. In addition, symmetric points correspond to the same pressure time evolution, as depicted in Figure 8.

Heterogeneous transversely isotropic material

In the perspective of an application to elastography imaging of biological tissues, we need to analyse more complex constitutive laws. For example, striated muscle tissue can be modelled as a transversely isotropic medium, i.e. there exists, at every point, a privileged direction represented by the unit vector τ 1 , [48] is

C ε(y) := µ ε(y) + η 1 9 tr ε(y) 1 - 1 3 tr ε(y) τ 1 ⊗ τ 1 - 1 3 (τ 1 • ε(y) • τ 1 ) 1 + (τ 1 • ε(y) • τ 1 ) τ 1 ⊗ τ 1 , (52)
with µ = 80 kPa, η = 3400 kPa. In what follows, we present the results for a transversely isotropic medium in which the fibre direction τ 1 varies linearly along the direction y. This configuration is inspired by the structure of myocardial tissue. In fact, it has been experimentally validated that muscle fibres are arranged in laminar structures, denoted sheets, of three to four muscle fibres in the thickness, that are oriented transversely to the heart wall [49; 50]. Moreover, the fibre orientation in human left ventricle myocardium changes linearly throughout the wall thickness, from -60 • close to the epicardium to +60 • near the endocardium [START_REF] Sommer | Biomechanical properties and microstructure of human ventricular myocardium[END_REF]. In order to model such a material, we consider an orientation of τ 1 varying linearly from -60 • at the bottom of the domain (y = 0) to +60 • at the top of the domain (y = 1). The resulting medium is highly heterogeneous. Figures 9 and10 show that the error worsens moderately in transversely isotropic media w.r.t. isotropic media. In addition, the order of convergence remains two for L 2 norm in space, and it is slightly degraded in L 2 (H 1 (Ω)) and in L ∞ (H 1 (Ω)). Figure 11 illustrates three snapshots at t = 0.5 s, t = 0.75 s and t = 1.1 s of the displacement field., respectively. The mesh is composed by N = 40 elements in each direction. We do not provide the results concerning the velocity and the acceleration, since they not provide additional information. Note that the pressure field gives a significant contri-10 -1.5 10 -1 10 -0.5 bution also in the interior of the domain, as it is illustrated in Figure 12, due to anisotropy.

A three-dimensional test case

The main application we have in mind is the propagation of elastic waves in nearly-incompressible biological tissues in the context of ultrasound transient elastography [52; 53]. Therefore, by way of illustration, we propose here a three-dimensional test case, inspired by elastography imaging of the myocardial tissue. Note that, due to dissipation, shear waves are fully attenuated in a few wavelength distance, i.e. some millimetres from the focal point of the ultrasonic beam [START_REF] Sarvazyan | Biomedical applications of radiation force of ultrasound: historical roots and physical basis[END_REF]. As a consequence, we can consider a small region of interest and a simple geometry. However, there is a need for a fine space discretisation, compatible with the wavelength of shear waves.

In this example we consider a transversely isotropic domain (defined by Eq. ( 52)), with fibre direction τ 1 oriented in the xy-plane and varying linearly along the direction z from -60 • w.r.t. the x-axis at the bottom of the domain (z = 0) to +60 • at the top of the domain (z = 0.02 m). The geometry of the domain is a parallelepiped of length m in x and y directions, and 0.02 m in z direction. This configuration represents an approximation of a region of interest in the left myocardium, selected in the anterior wall at the middle ventricular level. Density is set equal to ρ = 1050 kg m -1 and homogeneous Neumann boundary conditions are imposed on all the boundaries of the domain. High-order Spectral Finite Elements (of order 7 for the displacement, 6 for the pressure) are used for the space discretisation. The computational grid is composed of 24 uniform elements of size h = 1/24 in each direction, for a total of 14480427 degrees of freedom (DOF) for the displacement field, and 3048625 DOF for the pressure field. We adopt the time discretisation introduced in scheme (QIPnh) with penalisation coefficient α = 1/3ρ and we choose the time step as in Eq. ( 49), with ε = 0.2, as in Section 5. For each time step of the scheme, the Poisson problem for the scalar pressure field ( 19) is computed by means of our in-house fast solver based on Higher-Order Fourier Transform, as in Section 5. In this way, the time step used in the explicit time discretisation is constrained by the shear wave velocity only. Figure 13 depicts four snapshots of the solution at different time steps.

We have also considered an improved approximation by post-processing of the problem, by combination of the solutions associated with α 0 = 1/3ρ and α 1 = 2/5ρ, respectively, considering N = 24 elements in each direction. Since we have not noticed a qualitative change in the behaviour of the solution, we are confident that the observed solution is a good approximation of the pure incompressible problem. Note also that, due to the resolution of the Poisson problem by our in-house fast solver, we have drastically reduced the computational cost of the scheme. In particular, the resolution of 1600 time steps of the problem on a 12-cores workstation (cores at 2.7 GHz and 64 GB of RAM at 1867 MHz), considering higher-order extension of the scheme (i.e. compu- tation of the two solutions corresponding to α 0 and α 1 ) takes around 4 hours.

Conclusions

In this article we have outlined a new numerical scheme that is suitable for approximating elastic wave propagation in incompressible media, based on a mixed explicit/implicit time discretisation. We have demonstrated the time step of the resulting algorithm is only constrained by the shear wave velocity, and it is not limited by the enforcement of incompressibility (as it is the case for fully-explicit time discretisation). Furthermore, if effective methods are adopted for explicit time-discretisation, our algorithm only entails at each time step one resolution of a scalar Poisson problem -that can be performed by various, efficient algorithms -and few matrix-vector multiplications for the explicit part of the scheme. Moreover, we have presented a two-dimensional numerical test case, to demonstrate the favourable convergence properties of the scheme, and one three-dimensional example, to illustrate a realistic application to elastography imaging. A theoretical demonstration of second-order convergence in time of the scheme is out of the scope of this article and will be provided in future work.

List of abbreviations

FEM: Finite Element Method; LBB: Ladyzenskaya -Babuška -Brezzi condition; BVP: Boundary Value Problem; PDE: Partial differential equation.

Declarations

Fig. 1 :Fig. 2 :

 12 Fig.1: Convergence of scheme (QIP) w.r.t. the space step h, for final time t = 2 s, t = 2.65 s for a homogeneous isotropic material. Left: Relative L 2 (Ω T )error on the displacement. Right: Relative L ∞ (L 2 (Ω))-error on the displacement.

Fig. 3 :

 3 Fig. 3: Elastic wave propagation in a homogeneous isotropic medium (absolute value of the displacement field).

Fig. 4 :

 4 Fig. 4: Time evolution of the displacement field in a homogeneous isotropic medium in three locations (See Figure 3 for locations). Left : x direction. Right: y direction.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Elastic wave propagation in a homogeneous isotropic medium (absolute value of the acceleration field).

Fig. 7 :

 7 Fig. 7: Elastic wave propagation in a homogeneous isotropic medium (absolute value of the pressure field).

Fig. 8 :Fig. 9 :

 89 Fig.8: Time evolution of the pressure field in a homogeneous isotropic medium in three locations (See Figure7for locations).

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Convergence of scheme (QIP) w.r.t. the space step h, for final time t = 0.8 s, t = 1.1 s for a heterogeneous transversely isotropic material. Left: L 2 (H 1 (Ω))-error on the displacement. Right: L ∞ (H 1 (Ω))-error on the displacement.

  (a) t = 0.5 s (b) t = 0.75 s (c) t = 1.1 s

Fig. 12 :

 12 Fig. 12: Elastic wave propagation in a heterogeneous transversely isotropic medium (absolute value of the pressure field).

Fig. 13 :

 13 Fig. 13: Elastic wave propagation in a three-dimensional transversely isotropic medium (absolute value of the displacement field).

Acknowledgements

Not applicable.

Availability of data and material

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Funding

Not applicable.

Authors' contributions SI proposed the subject; FC and SI developed the method and implemented the scheme on an in-house finite element code (previously developed by SI); FC performed the finite element simulations and the numerical analysis under the technical supervision of SI; FC is the main author of the paper; critical revision was provided by SI.