George Bosilca

Aurelien Bouteiller

Thomas Herault

Yves Robert
email: yves.robert@ens-lyon.fr

Jack Dongarra

Assessing the impact of ABFT & Checkpoint composite strategies

Algorithm Based Fault Tolerant (ABFT) approaches promise unparalleled scalability and performance in failure-prone environments. With the advances in the theoretical and practical understanding of algorithmic traits enabling such approaches, a growing number of frequently used algorithms (including all widely used factorization) have been proven ABFTcapable. In the context of larger applications, these algorithms provide a temporal section of the execution when the data is protected by it's own intrinsic properties, and can be algorithmically recomputed without the need of checkpoints. However, while typical scientific applications spend a significant fraction of their execution time in library calls that can be ABFTprotected, they interleave sections that are difficult or even impossible to protect with ABFT. As a consequence, the only faulttolerance approach that is currently used for these applications is checkpoint/restart. In this paper we propose a model to investigate the efficiency of a composite protocol, that alternates between ABFT and checkpoint/restart for effective protection of an iterative application composed of ABFT-aware and ABFTunaware sections. We validate this model using a simulator. The model and simulator show that this composite approach drastically increases the performance delivered by an execution platform, especially at scale, by providing means to rarefy the checkpoints while simultaneously decreasing the volume of data needed to be checkpointed.

I. INTRODUCTION

As processor count increases with each new generation of high performance computing systems, the long dreaded reliability wall is materializing, and threatens to derail the efforts and milestones on the road toward Exascale computing. Despite continuous evolutions, such as improvements to the individual processor reliability, the integration of a large number of components leads, by simple probabilistic amplification, to a stern decrease in the overall capability of High Performance Computing (HPC) platforms to execute long running applications spanning a large number of resources. Already today, leadership systems encompassing millions of nodes experience a Mean Time Between Failures (MTBF) of a few hours [START_REF] Schroeder | Understanding failures in petascale computers[END_REF], [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF], [START_REF] Kondo | The failure trace archive: Enabling comparative analysis of failures in diverse distributed systems[END_REF]. Even considering an optimistic scenario with "fat" nodes, featuring heavily many-core systems and/or GPU accelerators, projections of Exascale machine exhibit unprecedented socket counts and will thereby suffer in terms of reliability [START_REF] Dongarra | The International Exascale Software Project: a Call To Cooperative Action By the Global High-Performance Community[END_REF].

The high performance community is not without resources to face this formidable threat. Under the already serious pressure that failures pose to currently deployed systems, checkpointing techniques have seen a large adoption, and many production quality software effectively provide protection against failures with application-level rollback recovery. During the execution, periodic checkpoints are taken that capture the progress of the application. When a failure occurs, the application is terminated, but can be later restarted from the last checkpoint. However, checkpointing techniques inflict severe overhead when failure frequency becomes too high. Checkpoints generate a significant amount of I/O traffic and often block the progression of the application; in addition, they must be taken more and more often as the MTBF decreases in order to enable steady progress of the application. Analytical projections clearly show that sustaining Exascale computing solely with checkpointing will prove challenging [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF], [START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF].

The fault-tolerance community has developed a number of alternative recovery strategies that do not employ checkpoint and rollback recovery as their premise. Strategies such as Algorithm Based Fault Tolerance (ABFT) [START_REF] Huang | Algorithm-based fault tolerance for matrix operations[END_REF], naturally fault tolerant iterative algorithms [START_REF] Shantharam | Characterizing the impact of soft errors on iterative methods in scientific computing[END_REF], resubmission in master slave applications, etc., can deliver more scalable performance under high stress from process failures. As an example, ABFT protection and recovery activities are not only inexpensive (typically less than 3% overhead observed in experimental works [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF], [START_REF] Davies | High Performance Linpack Benchmark: A Fault Tolerant Implementation without Checkpointing[END_REF]), but also have a negligible asymptotic overhead when increasing node count, which makes them extremely scalable. This is in sharp contrast with checkpointing, which suffers from increasing overhead with system size. ABFT is a useful technique for production systems, offering protection to important infrastructure software such as the dense distributed linear algebra library ScaLAPACK [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF]. In the remainder of this paper, without loss of generality, we will use the term ABFT broadly, so as to describe any technique that uses algorithm properties to provide protection and recovery without resorting to rollback recovery. However, typical HPC applications do spend some time where they perform computations and data movements that are incompatible with ABFT protection. The ABFT technique, as the name indicates, allows for tolerating failures only during the execution of the algorithm that features the ABFT properties. Moreover, it then protects only the part of the user dataset that is managed by the ABFT algorithm. In case of a failure outside the ABFT-protected operation, all data is lost; in case of a failure during the ABFT-protected operation, only the data covered by the ABFT scheme is restored. Unfortunately, these ABFT-incompatible phases force users to resort to generalpurpose (presumably checkpoint based) approaches as their sole protection scheme.

Yet, many HPC applications also spend quite a significant part of their total execution time inside a numerical library, and in many cases, these numerical library calls can be effectively protected by ABFT. We believe that the missing link enabling fault tolerance at extreme scale is the ability to effectively compose broad spectrum approaches (such as checkpointing) and algorithm based recovery techniques, as is most appropriate for different phases within a single application.

Possible target applications are based on iterative methods applied across an additional dimension such as time or temperature. Examples of such applications range from heat dissipation to radar cross-section, all of them being extremely time consuming applications, with the usual execution time for real-size problems ranging from several days to weeks. At the core of such applications, a system of linear equations is factorized, and the solution is integrated into a larger context, across the extra dimension. Upon closer inspection of the execution of such application, it becomes obvious that the most costly step is the factorization. Conveniently, factorization algorithms are some of the first algorithms to be extended with ABFT properties, both in the dense and sparse [START_REF] Shantharam | Fault tolerant preconditioned conjugate gradient for sparse linear system solution[END_REF], [START_REF] Agullo | Towards resilient parallel linear Krylov solvers: recover-restart strategies[END_REF], [START_REF] Chen | Algorithm-based recovery for iterative methods without checkpointing[END_REF] linear algebra world.

The main contribution of this paper is the design and evaluation of a new composite algorithm that allows for taking advantage of ABFT techniques in applications featuring phases for which no ABFT algorithm exists. We investigate a composition scheme corresponding to the above mentioned type of applications, where the computation alternates between ABFT protected and checkpoint protected phases. This composite algorithm imposes forced checkpoints when entering (and in some cases leaving) library calls that are protected by ABFT techniques, and uses traditional periodic checkpointing, if needed, between these calls. When inside an ABFT-protected call, the algorithm disables all periodic checkpointing. We describe a fault tolerance protocol that allows for switching between fault tolerance mechanisms, and depicts how different parts of the dataset are treated at each stage. Based on this scheme, we provide a performance model and use it to predict the expected behavior of such a composite approach on platforms beyond what is currently possible through experimentation. We validate the model by comparing its predicted performance to that obtained with a discrete event simulator.

The rest of the paper is organized as follows. We start with a brief overview of related work in Section II. Then we provide a detailed description of the composite approach in Section III, and derive the corresponding analytical performance model in Section IV. Section V is devoted to evaluating the approach, comparing the performance of traditional checkpointing protocols with that of the composite approach with realistic scenarios. This comparison is performed both analytically, instantiating the model with the relevant parameters, and in simulation, through an event-based simulator that we specifically designed for this purpose. We obtain an excellent correspondence between the model and the simulations, and we perform a weak-scalability study that demonstrates the full potential of the composite approach at very large scale. Finally, we provide concluding remarks in Section VI.

II. RELATED WORK

Both hardware and software errors can lead to application failure. The consequence of such errors can take various forms in a distributed system: a definitive crash of some processes, erroneous results or messages, or, at the extreme, corrupted processes exhibiting malignant behavior. In the context of HPC systems, most memory corruptions are captured by ECC memory or similar techniques, leaving process crashes as the most commonly observed type of failures.

The literature is rich in techniques that permit recovering the progress of applications when crash failures strike. The most commonly deployed strategy is checkpointing, in which processes periodically save their state, so that computation can be resumed from that point when some failure disrupts the execution. Checkpointing strategies are numerous, ranging from fully coordinated checkpointing [START_REF] Chandy | Distributed snapshots : Determining global states of distributed systems[END_REF] to uncoordinated checkpoint and recovery with message logging [START_REF] Elnozahy | A survey of rollback-recovery protocols in message-passing systems[END_REF]. Despite a very broad applicability, all of these fault tolerance methods suffer from the intrinsic limitation that both protection and recovery generate an I/O workload that grows with failure probability, and becomes unsustainable at large scale [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF], [START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF] (even when considering optimizations such as diskless or incremental checkpointing [START_REF] Plank | Diskless checkpointing[END_REF].)

In contrast, Algorithm Based Fault Tolerance (ABFT) is based on adapting the algorithm so that the application dataset can be recomputed at any moment, without involving costly checkpoints. ABFT was first introduced to deal with silent error in systolic arrays [START_REF] Huang | Algorithm-based fault tolerance for matrix operations[END_REF]. In recent work, the technique has been employed to recover from process failures [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF], [START_REF] Davies | High Performance Linpack Benchmark: A Fault Tolerant Implementation without Checkpointing[END_REF], [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF] in dense and sparse linear algebra factorizations [START_REF] Shantharam | Fault tolerant preconditioned conjugate gradient for sparse linear system solution[END_REF], [START_REF] Agullo | Towards resilient parallel linear Krylov solvers: recover-restart strategies[END_REF], [START_REF] Chen | Algorithm-based recovery for iterative methods without checkpointing[END_REF], but the idea extends widely to numerous algorithms employed in crucial HPC applications. So called Naturally Fault Tolerant algorithms can simply obtain the correct result despite the loss of portions of the dataset (typical of this are master-slave programs, but also iterative refinement methods, like GMRES or CG [START_REF] Shantharam | Characterizing the impact of soft errors on iterative methods in scientific computing[END_REF], [START_REF] Chen | Fault tolerant high performance computing by a coding approach[END_REF]). Although generally exhibiting excellent performance and resiliency, ABFT requires that the algorithm is innately able to incorporate fault tolerance and therefore stands as a less generalist approach. Another aspect that hinders its wide adoption and production deployment is that it can protect an algorithm and its dataset, but applications assemble many algorithms that operate on different datasets, and which may not all have a readily available ABFT version or employ different ABFT techniques.

To the best of our knowledge, this work is the first to introduce an effective protocol for alternating between generalist (typically checkpoint based) fault tolerance for some parts of the application and custom, tailored techniques (typically ABFT) for crucial, time consuming computational routines.

Many models are available to understand the behavior of checkpoint/restart [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF], [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF], [START_REF] Ozaki | Distribution-free checkpoint placement algorithms based on min-max principle[END_REF], [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF], and thereby to define an optimal checkpoint period. [START_REF] Zheng | Reliability-aware scalability models for high performance computing[END_REF] proposes a scalability model to evaluate the impact of failures on application performance. Compared with these works, we include several new key parameters to refine the model. A significant new contribution is to propose a generalized model for a protocol that alternates between checkpointing and ABFT sections. Although most ABFT methods have a complete complexity analysis (in terms of extra-flops, communications incurred by both protection activity and per-recovery cost [START_REF] Davies | High Performance Linpack Benchmark: A Fault Tolerant Implementation without Checkpointing[END_REF], [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF]), modeling the runtime overhead of ABFT methods under failure conditions has never been proposed. The composite model captures both the behavior of checkpointing and ABFT phases, as well as the cost of switching between the two approaches, and thereby permits investing the prospective gain from employing this mixed recovery strategy on extreme scale platforms.

III. COMPOSITE APPROACH

We consider a typical HPC application whose execution alternates GENERAL phases and LIBRARY phases (see Figure 1). During GENERAL phases, we have no information about the application behavior, and an algorithm-agnostic fault-tolerance technique, like checkpoint and rollback recovery, must be used. On the contrary, during LIBRARY phases, we know much more about the application, and we can apply special-purpose faulttolerance techniques, such as ABFT, to ensure resiliency.

During a GENERAL phase, the application can access the whole memory; during a LIBRARY phase, only the LIBRARY dataset (a subset of the application memory, which is passed as a parameter to the library call) is accessed. We call REMAIN-DER dataset the part of the application memory that does not belong to the LIBRARY dataset. A strong feature of ABFT is that, in case of failure, the ABFT algorithm can recompute the lost ABFT-protected data based only on the LIBRARY dataset of the surviving processors. The main goal of this paper is to compare two fault tolerant approaches: PUREPERIODICCKPT Pure (Coordinated) Periodic Checkpointing refers to the traditional approach based on coordinated checkpoints taken at periodic intervals, and using rollback recovery to recover from failures. ABFT&PERIODICCKPT Algorithm-Based Fault Tolerance & Periodic Checkpointing refers to the proposed algorithm, that combines ABFT techniques in LIBRARY phases with Periodic Checkpointing techniques in GENERAL phases.

It is described below.

Both approaches use PERIODICCKPT techniques, but to a different extent: while PUREPERIODICCKPT uses PERIODIC-CKPT throughout the execution, ABFT&PERIODICCKPT uses it only within GENERAL phases of the application.

A. ABFT&PERIODICCKPT Algorithm

The ABFT&PERIODICCKPT composite approach consists in alternating between periodic checkpointing and rollback recovery on one side, and ABFT on the other side, at different phases of the execution. Every time the application enters a LIBRARY phase (that can thus be protected by ABFT), a partial checkpoint is taken to protect the REMAINDER dataset. The LIBRARY dataset, accessed by the ABFT algorithm, need not be saved in that partial checkpoint, since it will be reconstructed by the ABFT algorithm inside the library call.

When the call returns, a partial checkpoint covering the modified LIBRARY dataset is added to the partial checkpoint taken at the beginning of the call, to complete it and to allow restarting from the end of the terminating library call. Said otherwise, the combination of the partial entry and exit checkpoints forms a split, but complete, coordinated checkpoint covering the entire dataset of the application.

If a failure is detected while processes are inside the library call, the crashed process is recovered using a combination of rollback recovery and ABFT. ABFT recovery is used to restore the LIBRARY dataset before all processes can resume the library call, as would happen with a traditional ABFT algorithm. The partial checkpoint is used to recover the REMAINDER dataset (everything except the data covered by the current ABFT library call) at the time of the call, and the process stack, thus restoring it before quitting the library routine, see Figure 2. The goal of this strategy is that ABFT recovery will spare some of the time spent redoing work, while periodic checkpointing can be completely de-activated during the library calls.

During GENERAL phases, regular periodic coordinated checkpointing is employed to protect against failures. In case of failure, coordinated rollback recovery brings all processes back to the last checkpoint (at most back to the split checkpoint capturing the end of the previous library call).

B. Efficiency Considerations and Application-Specific Improvements

A critical component to the efficiency of the PERIODIC-CKPT algorithm is the duration of the checkpointing interval. A short interval increases the algorithm overheads, by introducing many coordinated checkpoints, during which the application experiences slowdown, but also reduces the amount of time lost when there is a failure: the last checkpoint is never long ago, and little time is spent re-executing part of the application. Conversely, a large interval reduces overhead, but increases the time lost in case of failure. The PERIODICCKPT protocol has been extensively studied, and good approximations of the optimal checkpoint interval exist (known as Young and Daly's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF], [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]). These approximations are based on the machine MTBF, checkpoint duration, and other parameters. We will consider two forms of PERIODICCKPT algorithms: the PUREPERIODICCKPT algorithm, where a single checkpointing interval is used consistently during the whole execution, and the BIPERIODICCKPT algorithm, where the checkpointing interval may change during the execution, to fit different conditions (see Section IV-C, Figures 5 and6).

In the ABFT&PERIODICCKPT algorithm, we interleave PERIODICCKPT protected phases with ABFT protected phases, during which periodic checkpointing is de-activated. Thus, different cases have to be considered:

• When the time spent in a GENERAL phase is larger than the optimal checkpoint interval, periodic checkpointing is used during these phases in the case of ABFT-&PERIODICCKPT (see Figure 3) • When the time spent in a GENERAL phase is smaller than the optimal checkpoint interval, the ABFT&PERIODIC-CKPT algorithm already creates a complete valid checkpoint for this phase (formed by combining the entry and exit partial checkpoints), so the algorithm will not introduce additional checkpoints (see Figure 4).

Moreover, the ABFT&PERIODICCKPT algorithm forces (partial) checkpoints at the entry and exit of library calls; thus if the time spent in a library call is very small, this approach will introduce more checkpoints than a traditional PERIODICCKPT approach. The time complexity of library algorithms usually depends on a few input parameters related to problem size and resource number, and ABFT techniques have deterministic, well known time overhead complexity. Thus, when possible, the ABFT&PERIODICCKPT algorithm features a safeguard mechanism: if the projected duration of a library call with ABFT protection (computed at runtime thanks to the call parameters and the algorithm complexity) is smaller than the optimal periodic checkpointing interval, then ABFT is not activated, and the corresponding LIBRARY phase is protected using the PERIODICCKPT technique only.

Furthermore, since only a subset of the entire dataset is modified during a library call (the LIBRARY dataset), incremental checkpointing techniques can benefit PERIODICCKPT approaches. This consists of saving only the subset of the memory that has been modified since the last checkpoint, when taking a new process checkpoint. This influences the duration of the checkpointing operation, and thus the optimal checkpoint interval. In our models, we will take this parameter in consideration.

IV. MODEL

In this section, we detail the application model and the various parameters used to quantify the cost of checkpointing and ABFT protection. Then we analytically derive the minimal overhead for all scenarios. In Section IV-A, we start by defining some parameters, and then proceed in Section IV-B with determining the cost of the composite approach. We compare this cost to that of classical approaches in Section IV-C.

A. Application and checkpoint parameters

The execution of the application is partitioned into epochs. Within an epoch, there are two phases for the application: the first phase is spent outside the library (it is a GENERAL phase), and only periodic checkpointing can be employed to protect from failures during that phase. Then the second phase (a LIBRARY phase) is spent into a compute intensive library routine that has the potential to be protected by ABFT. Such a scenario is very general, and many scientific applications obey this scheme, alternating phases spent outside and within a library call that can be protected by ABFT techniques. Since each epoch can be analyzed independently, without loss of generality, we focus on a single epoch. Let us introduce some notations. The total duration of the epoch is T 0 = T G + T L , where T G and T L are the durations for the GENERAL and LIBRARY phases, respectively. Let α be the fraction of time spent in a LIBRARY phase: then we have

T L = α × T 0 and T G = (1 -α) × T 0 .
As mentioned earlier, another important parameter is the amount of memory that is accessed during the LIBRARY phase (the LIBRARY dataset). This parameter is important because the cost of checkpointing in each phase is directly related to the amount of memory that needs to be protected. The total memory footprint is M , and the associated checkpointing cost is C. We write M = M L + ML, where M L is the size of the LIBRARY dataset, and ML is the size of the REMAINDER dataset. Similarly, we write C = C L + CL, where C L is the cost of checkpointing M L , and CL the cost of checkpointing ML. We can define the parameter ρ that defines the relative fraction of memory accessed during the LIBRARY phase by M L = ρM , or, equivalently, by C L = ρC.

B. Cost of the composite approach

We now detail the cost of resilience during each phase of the composite approach. We start with the intrinsic cost of the method itself, i.e., assuming a fault-free execution. Then we account for the cost of failures and recovery.

1) Fault-free execution: During the GENERAL phase, we separate two cases. First, if the duration T G of this phase is short (smaller than P G , defined below), then we simply take a partial checkpoint at the end of this phase, before entering the ABFT-protected mode. This checkpoint is of duration CL, because we need to save only the REMAINDER dataset before switching modes. Otherwise, if T G is larger than P G , we rely on periodic checkpointing during the GENERAL phase: more specifically, the regular execution is divided into periods of duration P G = W + C. Here W is the amount of work done per period, and the duration of each periodic checkpoint is C = CL + C L , because the whole application footprint must be saved during a GENERAL phase. The optimal value of P G will be computed below. Without loss of generality, we assume an integer number of periods, and the last periodic checkpoint replaces that of size CL preceding the switch to ABFT-protected mode.

Altogether, the length T ff G of a fault-free execution of the GENERAL phase is the following:

• If T G < P G , then T ff G = T G + CL • Otherwise, we have T G
W periods of length P G , so that

T ff G = T G P G -C × P G (1)
Now consider the LIBRARY phase: we use the ABFTprotection algorithm, whose cost is modeled as an affine function of the time spent: if the computation time of the library routine is t, its execution with the ABFT-protection algorithm becomes φ × t. Here, φ > 1 accounts for the overhead paid per time-unit in ABFT-protected mode. This linear model for the ABFT overhead fits the existing algorithms for linear algebra, but other models could be considered.

In addition, we pay a checkpoint C L when exiting the library call (to save computed data). Therefore, the fault-tree execution time is

T ff L = φ × T L + C L (2)
Finally, the fault-free execution time of the whole epoch is

T ff = T ff G + T ff L ( 3
)
where T ff G and T ff L are computed according to the Equations ( 1) and (2).

2) Cost of failures: Next we have to account for failures. During t time units of execution, the expectation of the number of failures is t µ , where µ is the mean time between failures of the platform. Note that if the platform comprises N identical resources whose individual mean time between failures is µ ind , then µ = µind N . This relation is agnostic of the granularity of the resources, which can be anything from a single CPU to a complex multi-core socket.

For each phase, we have a similar equation: the final execution time is the fault-free execution time, plus the number of failures multiplied by the (average) time lost per failure:

T final G = T ff G + T final G µ × t lost G ( 4
)
T final L = T ff L + T final L µ × t lost L (5)
Equation ( 4) reads as follows: T ff G is the failure-free execution time, to which we add the time lost due to failures; the expected number of failures is

T final G µ ,

and t lost

G is the average time lost per failure. We have a similar reasoning for Equation [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF]. Then,t lost G and t lost L remain to be computed. For t lost G (GENERAL phase), we discuss both cases:

• If T G < P G : since we have no checkpoint until the end of the GENERAL phase, we have to redo the execution from the beginning of the phase. On average, the failure strikes at the middle of the phase, hence the expectation of loss is T ff G 2 time units. We then add the downtime D (time to reboot the resource or set up a spare) and the recovery R.

Here R is the time needed for a complete reload from the checkpoint (and R = C if read/write operations from/to the stable storage have the same speed). We derive that

t lost G = D + R + T ff G 2 (6)
• If T G > P G : in this case, we have periodic checkpoints, and the amount of execution which needs to be re-done after a failure corresponds to half a checkpoint period on average, so that

t lost G = D + R + P G 2 (7)
For t lost L (LIBRARY phase), we derive that t lost L = D + RL + Recons ABFT Here, RL is the time for reloading the checkpoint of the REMAINDER dataset (and in many cases RL = CL). As for the LIBRARY dataset, there is no checkpoint to retrieve, but instead it must be reconstructed from the ABFT checksums, which takes a time Recons ABFT .

3) Optimization: We check from Equations ( 2) and ( 5) that T final L is always a constant. Indeed, we derive that

T final L = 1 1 -D+R L+ReconsABFT µ × (φ × T L + C L ) (8)
As for T final G , it depends on the value of T G : it is constant when T G is small. In that case, we derive that

T final G = 1 1 - D+R+ T G +C L 2 µ × (T G + CL) (9)
The interesting case is when T G is large: in that case, we have to determine the optimal value of the checkpointing period P G which minimizes T final G . From Equations ( 1), ( 4) and ( 7), we derive that

T final G = T G X where X = 1 - C P G 1 - D + R + P G 2 µ (10) We rewrite X = (1 - C 2µ ) - P G 2µ - C(µ -D -R) µP G
The maximum of X gives the optimal period P opt G . Differentiating X as a function of P G , we find that it is obtained for

P opt G = 2C(µ -D -R) (11)
Plugging the value of P opt G back into Equation ( 10) provides the optimal value of T final G when T G is large.

We have successfully computed the final execution time T final of our composite approach in all cases. In the experiments provided in Section V, we report the corresponding The waste is defined as the fraction of time when platform resources do not progress the application's computation (due to the intrinsic overhead of the technique and to failures that strike the application during execution). The waste is given by

WASTE = 1 - T 0 T final (12)
We conclude this section with a word of caution: the optimal value of the waste is only a first-order approximation, not an exact value. Equation ( 11) is a refined version of well known formulas by Young [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF] and Daly [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]. But just as in [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF], [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF], the formula only holds when µ, the value of the MTBF, is large in front of the other parameters. Owing to this hypothesis, we can neglect the probability of several failures occurring during the same checkpointing period. However, to further assess the accuracy of the model, when doing simulations in the experiments, we account for all unlikely scenarios and re-execute the work until each period is successfully completed.

C. Comparison with conservative approaches

A fully conservative approach, agnostic of the ABFT library, would perform periodic checkpoints throughout the execution of the whole epoch. As already mentioned, we call this approach PUREPERIODICCKPT (see Figure 5). Let T final P C be the final execution time with this PUREPERIODICCKPT approach; it can be computed from the results of Section IV-B as follows:

• No ABFT: α = 0 and T final L = 0 • We optimize T final P C = T final G
just as before, with the same optimal period P opt P C = P opt G , employed throughout the epoch.

One can reduce the cost of PUREPERIODICCKPT by noticing that during the LIBRARY epoch, only the LIBRARY dataset is modified (namely M L ). Employing incremental checkpointing would, in this case, yield a checkpoint cost reduction (down to C L ). Obviously, with a different cost of checkpointing, the optimal checkpoint period is different. Therefore, a semi-conservative approach (called BIPERIODICCKPT, see Figure 6) assumes that the checkpoint system can recognize that the program has entered a library routine that modifies only a subset of the dataset, and switches to the optimal checkpoint period according to the application phase. During the GENERAL phase, the overhead of failures and protection remains unchanged, but during the LIBRARY phase, the cost of a checkpoint is reduced to C L (instead of C); however, the cost of reloading from a checkpoint remains R (since the different incremental checkpoints must be combined to recover the entire dataset at rollback time). This leads to two different checkpointing periods, one for each phase. The new optimal checkpoint period can be modeled as follows:

• T final P C = T final G + T final LP C , where T final G is computed as before • T final
LP C is computed similarly as T final G , but with different parameters:

T final LP C = 1 1 - D+R+ P BP C 2 µ × P BP C P BP C -C L × T L (13)
and the optimal period is

P opt BP C,L = 2C L (µ -D -R) (14)
V. EVALUATION

In this section, we evaluate the ABFT&PERIODICCKPT protocol in simulation, and compare its performance to PURE-PERIODICCKPT and BIPERIODICCKPT in different scenarios. We start with a description of the simulator and experiments in Section V-A. Then we detail the results of the comparison of the different protocols in Section V-B. In Section V-B, we compare simulation results and predicted performance results analytically computed from the models presented in Sections IV-B and IV-C, and we do obtain a very good correspondence. Finally, we conduct a weak scalability study in Section V-C, in order to assess the performance of the various protocols at very large scale.

A. Validation

To validate the performance models, we have implemented a simulator, based on discrete event simulation, that reproduces the behavior of the different algorithms, even in cases that the performance models cannot cover. Indeed, as mentioned in Section IV-B3, a few approximations have been made when considering the mathematical models, to make their expressions tractable. For example, the models assume that a single failure may hit the system, until its recovery. The effect of events like overlapping failures, which is uncommon when the MTBF is large enough, is neglected in the proposed performance model. The simulator, however, takes these events into account, accurately reproducing the corresponding costs.

In the simulator, failures are generated following an Exponential distribution law parameterized to fix the MTBF to a given value. Then the application, and the chosen fault tolerance mechanism, are unfolded on that set of failures, triggering rollbacks, and other protocol-specific overheads, to measure the duration of the execution. For each scenario, and each parameter, the average termination time over a thousand executions is returned by the simulator.

We present in [START_REF] Bosilca | Assessing the impact of ABFT and checkpoint composite strategies[END_REF] an exhaustive evaluation of the different parameters independently, comparing the results as predicted by the models, and the simulation. In this paper, we focus the analysis on a smaller subset. We consider an application that executes for a week when there is neither a fault tolerance mechanism nor any failure. The time to take a checkpoint and rollback the whole application is 10 minutes (C, R), a consistent order of magnitude for current applications at large scale [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF]. We consider that the ratio of the memory that is modified by the LIBRARY phase (ρ) is fixed at 0.8 (to vary a single parameter at a time in our simulation), and the overhead due to ABFT is φ = 1.03 (again, typical from production deployments [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF]).

Figure 7 presents 6 evaluations of that scenario. The MTBF of the system varies on the x-axis, and the ratio of time spent in the LIBRARY phase (α) on the y-axis. In Figures 7a to 7f, we present the waste predicted by the model, and validate the model by observing the difference between the model prediction and the waste measured from the simulator for a given combination of parameters and protocol. From the validation perspective, the figures on the right side show an excellent correspondence between predicted (from the model) and actual (obtained from simulation) values. For small MTBF values, the model tends to slightly underestimate the waste. That underestimation does not exceed 12% in the worst case and quickly decreases to below 5%. Qualitatively, this underestimation is expected, because an approximation that must be done to allow a closed formula representation is to assume that failures will not hit processors while they are recovering from a previous failure. In reality, when the MTBF is very small, this event can sometimes happen, forcing the system to start a new recovery, and introducing additional waste.

B. PUREPERIODICCKPT, BIPERIODICCKPT, and ABFT-&PERIODICCKPT

Consider Figures 7a and7b, that represent the waste of PUREPERIODICCKPT as a function of the MTBF (µ) and the amount of time spent in the LIBRARY routine (α): it is obvious that the PUREPERIODICCKPT protocol, which is oblivious of the different phases of the application, presents a waste that is only a function of the MTBF. As already evaluated explained in many other works, when the MTBF increases, the waste decreases, because the overheads due to failure handling tend toward 0, and the optimal checkpointing period can increase significantly, reducing the waste due to resilience in a fault-free execution.

Comparatively, for the protocol BIPERIODICCKPT presented in Figures 7c and7d, the parameter α affects the optimal periods used both in the LIBRARY and general phases. Since the cost of checkpointing for these phases differs by 20% (C L = 0.8C), when the relative time spent in the GENERAL routine increases (α is closer to 0), then the protocol behaves more and more as PUREPERIODICCKPT. When α is almost 1, on the contrary, the behavior is similar to PUREPERIODIC-CKPT, but with a checkpoint cost reduction of 20%. Thus, the waste becomes minimal when α tends toward 1.

In Figures 7e and7f, we present the waste for the ABFT-&PERIODICCKPT protocol. When α tends toward 0, as above, the protocol behaves as PUREPERIODICCKPT, and no benefit is shown. When 50% of the time is spent in the LIBRARY routine, the benefit, compared to PUREPERIODICCKPT, but also compared to BIPERIODICCKPT, is already visible: for 50% of the failures (when the failure hits during a LIBRARY phase), the cost of recovery is reduced to 20% of the rollback cost, plus the constant overhead of ABFT recovery. Moreover, periodic checkpointing is disabled 50% of the time, producing yet another gain compared to BIPERIODICCKPT which still requires saving 80% of the memory periodically. In this case, the gain in checkpoint avoidance compensates for the waste induced by additional computations done during the LIBRARY phase to provide the ABFT protection. When considering the extreme case of 100% of the time spent in the LIBRARY phases, the overhead tends to reach the overhead induced by the slowdown factor of ABFT (φ = 1.03, hence 3% overhead).

C. Weak Scalability

As illustrated above, the ABFT&PERIODICCKPT approach exhibits better performance when a significant time is spent in the LIBRARY phase, and when the failure rate implies a small optimal checkpointing period. If the check- pointing period is large (because failures are rare), or if the duration of the LIBRARY phase is small, then the optimal checkpointing interval becomes larger than the duration of the LIBRARY phase, and the algorithm automatically resorts to the BIPERIODICCKPT protocol. This can also be the case when the epoch itself is smaller than (or of the same order of magnitude as) the optimal checkpointing interval (i.e., when the application does a fast switching between LIBRARY and GENERAL phases).

However, consider such an application that frequently switches between (relatively short) LIBRARY and GENERAL phases. When porting that application to a future larger scale machine, the number of nodes that are involved in the execution will increase, and at the same time, the amount of memory on which the ABFT operation is applied will grow (following Gustafson's law). This has a double impact: the time spent in the ABFT routine increases, while at the same time, the MTBF of the machine decreases. In this section, we evaluate quantitatively how this scaling factor impacts the relative performance of the ABFT&PERIODICCKPT, PURE-PERIODICCKPT and BIPERIODICCKPT algorithms. Owing to the good correspondence between results from the model and results from the simulation, we (confidently) use only the model in this scalability study.

First, we consider the case of an application where the LIBRARY and GENERAL phases scale at the same rate. We take the example of linear algebra kernels operating on 2D-arrays (matrices), that scale in O(n 3 ) of the array order n (in both phases). Following a weak scaling approach, the application uses a fixed amount of memory M ind per node, and when increasing the number x of nodes, the total amount of memory increases linearly as M = xM ind . Thus O(n 2 ) = O(x), and the parallel completion time of the O(n 3 ) operations, assuming perfect parallelism, scales in O( √ x).

To instantiate this case, we take an application that would iterate over a thousand epochs, each epoch consisting of 80% of a LIBRARY phase, and 20% of a GENERAL phase. At 10,000 nodes, the duration of a single epoch is arbitrarily set to 1 minute, and the scaling factor of the corresponding O(n 3 ) operation is applied, when varying the number of nodes that participate in the computation. We set the duration of the complete checkpoint and rollback (C and R, respectively) to 1 minute when 10,000 nodes are involved, and we scale this value linearly with the total amount of memory, when varying the number of nodes. The MTBF at 10,000 nodes is set to 1 failure every day, and this also scales linearly with the number of components. The ABFT overheads, and the downtime, are set to the same values as in the previous section, and 80% of the application memory (M L ) is touched by the LIBRARY phase.

Given these parameters, Figure 8 shows (i) the relative waste of PUREPERIODICCKPT, BIPERIODICCKPT, and ABFT&PERIODICCKPT, as a function of the number of nodes, and (ii) the average number of faults that each execution will have to deal with to complete. The expected number of faults is the ratio of the application duration by the platform MTBF (which decreases when the number of nodes increases, generating more failures). The fault-free execution time increases with the number of nodes (as noted above), and the fault-tolerant execution time is also increased by the waste due to the protocol. Thus, the total execution time of PURE-PERIODICCKPT or BIPERIODICCKPT is larger at 1 million nodes than the total execution time of ABFT&PERIODICCKPT at the same scale, which explains why more failures happen for these protocols.

When comparing BIPERIODICCKPT and PUREPERIODIC-CKPT, one can see the benefit of incremental checkpointing, which spares about 20% of the checkpoint time during 80% of the checkpoints: this benefit shows up by a small linear reduction of the waste for BIPERIODICCKPT. However, both approaches perform similarly with respect to the number of nodes in this weak-scaling experiment.

Up to approximately 100,000 nodes, the fault-free overhead of ABFT negatively impacts the waste of the ABFT-&PERIODICCKPT approach, compared to BIPERIODICCKPT or PUREPERIODICCKPT. Because the MTBF of the platform is very large compared to the application execution time (and hence to the duration of each LIBRARY phase), periodic checkpointing approaches have a very large checkpointing interval, introducing very few checkpoints, thus a small failurefree overhead. Because failures are rare, the cost due to time lost at rollbacks does not overcome the benefits of a small failure-free overhead, while the ABFT technique must pay the linear overhead of maintaining the redundancy information during the whole computation of the LIBRARY phase.

Once the number of nodes reaches 100,000, however, two things happen: failures become more frequent, and the time lost due to failures starts to impact rollback recovery approaches. Thus, the optimal checkpointing interval of periodic checkpointing becomes smaller, introducing more checkpointing overheads. During 80% of the execution, however, the ABFT&PERIODICCKPT approach can avoid these overheads, and when they reach the level of linear overheads due to the ABFT technique, ABFT&PERIODICCKPT starts to scale better than both periodic checkpointing approaches.

All protocols have to resort to checkpointing during the GENERAL phase of the application. Thus, if failures hit during this phase (which happens 20% of the time in this example), they will all have to resort to rollbacking and lose some computation time. Hence, when the number of nodes increases and the MTBF decreases, eventually, the time spent in rollbacking and re-computing, which is linear in the number of faults, will increase the waste of all algorithms. However, one can see that this part is better controlled by the ABFT&PERIODICCKPT algorithm.

Next we consider the case of an unbalanced GENERAL phase: consider an application where the LIBRARY phase has a cost O(n 3 ) (where n is the problem size), as above, but where the GENERAL phase consists of O(n 2 ) operations. This kind of behavior is reflected in many applications where matrix data is updated or modified between consecutive calls to computation kernels. Then, the time spent in the LIBRARY phase will increase faster with the number of nodes than the time spent in the GENERAL phase, varying α. This is what is represented in Figure 9. We took the same scenario as above for Figure 8, but α is a function of the number of nodes chosen such that at 10,000 nodes, α = T final L /T final = 0.8, and everywhere,

T final L = O(n 3 ) = O( √ x), and
T final P C = O(n 2 ) = O(1)
. We give the value of α under the number of nodes, to show how the fraction of time spent in LIBRARY phases increases with the number of nodes.

The PUREPERIODICCKPT protocol is not impacted by this change, and behaves exactly as in Figure 8. Note, however, that T final = T final L + T final P C progresses at a lower rate in this scenario than in the previous scenario, because T final P C does not increase with the number of nodes. Thus, the average number of faults observed for all protocols is much smaller in this scenario. Because more and more time (relative to the duration of the application) is spent in the LIBRARY phase, where 20% of the memory does not need to be saved, the BIPERIODICCKPT algorithm increases its benefit, compared to PUREPERIODICCKPT: less overhead is paid for checkpoints that happen during LIBRARY phases, and the optimal period of checkpointing during these phases are longer. The cost of failures, however, remains the same, since the state of the entire application (LIBRARY memory, and REMAINDER memory) must be restored at rollback time.

The efficiency on ABFT&PERIODICCKPT, however, is more significant. The latter protocol benefits from the increased α ratio in both cases: since more time is spent in the LIBRARY phase, periodic checkpointing is de-activated for relatively longer periods. Moreover, this increases the probability Figure 10: Total waste for ABFT&PERIODICCKPT, BI-PERIODICCKPT and PUREPERIODICCKPT, when considering the weak scaling of an application with variable ratio of time spent in a LIBRARY routine, and constant checkpointing time that a failure will happen during the LIBRARY phase, where the recovery cost is greatly reduced using ABFT techniques. Thus, ABFT&PERIODICCKPT is capable of mitigating failures at a much smaller overhead than simple periodic checkpointing, and more importantly with better scalability.

In both previous evaluations, we have always considered that the checkpointing (and rollback recovery) time is proportional to the global amount of memory that needs to be saved in these checkpoints. This is realistic, if the checkpoint needs to be stored in a remote place, to guarantee its availability after a failure occurs. In this case, the interconnect (or the bandwidth capacity of the disks) eventually becomes a bottleneck, and the saving time becomes proportional to the number of computing resources saving their state simultaneously. To mitigate the negative effect of this bottleneck, system designers are studying a couple of alternative approaches. One consists of featuring each computing node with local storage capability, ensuring through the hardware that this storage will remain available during a failure of the node. Another approach consists of using the memory of the other processors to store the checkpoint, pairing nodes as "buddies," thus taking advantage of the high bandwidth capability of the high speed network to design a scalable checkpoint storage mechanism [START_REF] Zheng | FTC-Charm++: an in-memory checkpoint-based fault tolerant runtime for Charm++ and MPI[END_REF], [START_REF] Ni | Hiding checkpoint overhead in HPC applications with a semi-blocking algorithm[END_REF], [START_REF] Dongarra | Revisiting the double checkpointing algorithm[END_REF], [START_REF] Rajachandrasekar | A 1 PB/s file system to checkpoint three million MPI tasks[END_REF].

Thus, it might be reasonable to consider in the future that the checkpoint storage time will not increase with the number of nodes, but on the contrary will remain constant. This is the scenario that we contemplate in Figure 10. The scenario is identical to the previous scenario of Figure 9, but the checkpoint time (C) and rollback recovery time (R) is independent of the number of nodes that checkpoint, and is fixed at 60s. One can see a noteworthy benefit on both periodic checkpointing protocols: even at 1 million nodes, the waste due to the protocols and the few faults that do occur during the execution (up to 6 failures in average during the whole execution) both add up to below 15%. At the same time, the ABFT technique continues to introduce its constant overhead (due to additional computation) during the whole execution, and appears to present a waste that is almost constant when the number of nodes increases.

Figure 10 shows that PUREPERIODICCKPT and BI-PERIODICCKPT are less efficient than ABFT&PERIODIC-CKPT at 1 million nodes, despite the perfectly scalable checkpointing hypothesis. To reach comparable performance, we must reduce checkpointing overhead by a factor of 10 and use C = R = 6s. Such low figures can only be achieved through new hardware (like NVRAM), and new incremental and hierarchical checkpointing protocols.

VI. CONCLUSION

In this paper, we have formalized and quantified a novel method of composing fault tolerance approaches for applications that alternate between ABFT-aware and ABFT-unaware sections, where each of these sections is protected by its own mechanism, ABFT in one case and checkpoint/restart in the other. A performance model has been derived for such methods and thoughtfully validated using a simulator developed for this scope. We have compared our composite approach with a traditional periodic checkpointing approach using rollback and recovery, under different plausible scenarios. Our model predicts that the cost of a "checkpoint only" approach will maintain a reasonable overhead only under highly optimistic assumptions, where the checkpointing cost stagnates when the number of computational resources increases. Under more realistic assumptions, where the checkpointing cost increases with the number of resources, the composite approach will provide significantly greater benefits compared with checkpoint/restart, by minimizing the waste and thus increasing the platform throughput. Our weak scalability study shows that the gain of the composite approach will continue to grow with the increase in the number of computing resources, making it a more plausible and desirable approach at very large scale.

Figure 1 :

 1 Figure 1: Typical Application

Figure 2 :

 2 Figure 2: ABFT&PERIODICCKPT composite approach

Figure 3 :Figure 4 :

 34 Figure 3: ABFT&PERIODICCKPT composite (time spent in GENERAL phases is large)

Figure 5

 5 Figure 5: PUREPERIODICCKPT

Figure 7 :

 7 Figure 7: Waste as a function of MTBF and fraction of time α spent in LIBRARY phase; difference of the measured waste by simulation WASTE simul minus the predicted waste by the model WASTE model . Here W =1 week, C =R =10 minutes, C L =0.8C, φ =1.03, Recons ABFT =2s.

Figure 8 :

 8 Figure 8: Total waste for ABFT&PERIODICCKPT, BI-PERIODICCKPT and PUREPERIODICCKPT, when considering the weak scaling of an application with a fixed ratio of 80% spent in a LIBRARY routine.

Figure 9 :

 9 Figure 9: Total waste for ABFT&PERIODICCKPT, BI-PERIODICCKPT and PUREPERIODICCKPT, when considering the weak scaling of an application with variable ratio of time spent in a LIBRARY routine.

ACKNOWLEDGEMENTS

This work was supported in part by the National Science Foundation (NSF #0904952 and #1063019), JST Japan, and the French Research Agency (ANR) through the Rescue project. Yves Robert is with Institut Universitaire de France.


P Bouchon

M Makhsiyan

J Jaeck

J.-L Pelouard

Minao

COMMUNICATION A CONGRES Emission thermique à base de nano-antennes

EMISSION THERMIQUE A BASE DE NANOANTENNES

Patrick Bouchon 1 , Mathilde Makhsiyan 1 , Julien Jaeck 1 , Jean-Luc Pelouard 2 , Riad Haïdar

ENCODAGE D'IMAGES POLARISEES ET MULTISPECTRALES

Nanophotonique

OPTIQUE 2015

Les nanoantennes MIM ont fait l'objet de nombreuses études, mais qui se sont limitées à des répétitions périodiques. Ces nanoantennes se comportent comme des résonateurs Fabry-Perot, dont la longueur va fixer la fréquence de résonance pour la polarisation transverse magnétique. La figure 1(a) montre une photographie de la métasurface inhomogène fabriquée, et une image au microscope électronique à balayage (MEB) de neuf cellules élémentaires, chacune contenant plusieurs antennes métalliques différentes. Les spectres calculés entre 3 et 6 µm pour chaque polarisation montre qu'il est possible d'obtenir une émissivité unitaire, mais également de juxtaposer ces émissivités à différentes longueurs d'onde et polarisation.

CONCLUSION

La démonstration de cette métasurface inhomogène ouvre plusieurs perspectives pour la création de mires infrarouge, de dispositifs anti contrefaçons ou pour la détection d'espèces biochimiques. Par ailleurs, cette approche peut être appliqué à d'autres antennes (split-ring, résonateur de Helmholtz optique,…)

  MOTS-CLEFS : Emission thermique ; nanoantennes métalliques1. INTRODUCTIONLes métasurfaces, l'équivalent 2D des métamatériaux, permettent de contrôler la réponse électromagnétique sur des plages spectrales variées. En particulier, il est possible de manipuler l'émissivité thermique, et des donner des propriétés d'émission monochromatiques ou directionnelles à des surfaces nanostructurées[1]. Dans la suite, nous présentons l'utilisation d'une métasurface inhomogène, composée de 100 millions de nano-antennes metal-isolant-métal (MIM), pour obtenir l'émission d'images à différentes polarisation et longueur d'onde[2-3]. Chaque antenne agit comme un émetteur à une polarisation et une longueur d'onde données, indépendant des autres, avec une section géométrique très faible devant sa section efficace d'émission.

Fig. 1 :

 1 Fig. 1 : (a) Image photo de l'échantillon et image MEB de neuf cellules élémentaires contenant des antennes. (b) Spectre d'émissivité calculé, associé à chaque cellule élémentaire périodisé de l'image MEB, pour chacune des polarisations.

Fig. 2 :

 2 Fig. 2 : (a) Image photo de l'échantillon. (b) Emission infrarouge de l'échantillon dans la polarisation horizontale et peinture originale de Nicolas Mignard. (c) Emission infrarouge de lettres N, A, O dans la polarization vertcale à plusieurs longueurs d'onde (toutes les barres d'échelle font 5 mm). La figure 2 montre l'émission thermique de cet échantillon chauffé à 100°C, et observé par une caméra infrarouge à haute résolution dans la bande 3-5 µm. Un polariseur est ajouté au montage pour observer l'image de Molière (Fig. 2(b)), et lorsque le polariseur est tourné à 90°, des filtres spectraux sont rajoutés qui permettent d'observer les lettres N, A et O à différentes longueurs d'onde.

1 1

 1 MinaoONERA, The French Aerospace Lab, F-91761 Palaiseau, France France 2 MiNaO, Laboratoire de Photonique et de Nanostructures (LPN), CNRS, Université Paris-Saclay, route de Nozay, F-91460 Marcoussis, France

	Patrick.bouchon@onera.fr
	RÉSUMÉ Les nanoantennes métalliques permettent de manipuler la lumière, et en particulier de contrôler son absorption. D'après la loi de Kirchhoff, l'émissivité est égale à l'absorptivité, et des nanoantennes peuvent devenir des sources de lumière fonctionnant par émission thermique. L'originalité de ces sources est d'avoir des propriétés qui s'écartent de celles du corps noir idéal décrit par la loi de Planck. Dans cet exposé, je montrerai qu'il est possible de développer une métasurface à base d'ensemble de nanoantennes métal-isolant-métal, où chaque antenne a des dimensions petites devant la longueur d'onde, et va se comporter comme un émetteur à une polarisation et une longueur d'onde donnée

, indépendant des antennes voisines. Il est ainsi possible d'obtenir un contrôle spatial, spectral et en polarisation de la lumière, et d'encoder des images. D'autres concepts de nanoantennes seront présentés, comme le résonateur de Helmholtz optique, qui permettent de contrôler la finesse spectrale de l'émission thermique ou le taux de réjection.


Clara Issanchou
email: issanchou@lam.jussieu.fr

Jean-Loïc Le Carrou

Stefan Bilbao

Cyril Touzé
email: cyril.touze@ensta-paristech.fr

Olivier Doaré

A MODAL APPROACH TO THE NUMERICAL SIMULATION OF A STRING VIBRATING AGAINST AN OBSTACLE: APPLICATIONS TO SOUND SYNTHESIS

A number of musical instruments (electric basses, tanpuras, sitars...) have a particular timbre due to the contact between a vibrating string and an obstacle. In order to simulate the motion of such a string with the purpose of sound synthesis, various technical issues have to be resolved. First, the contact phenomenon, inherently nonlinear and producing high frequency components, must be described in a numerical manner that ensures stability. Second, as a key ingredient for sound perception, a fine-grained frequencydependent description of losses is necessary. In this study, a new conservative scheme based on a modal representation of the displacement is presented, allowing the simulation of a stiff, damped string vibrating against an obstacle with an arbitrary geometry. In this context, damping parameters together with eigenfrequencies of the system can be adjusted individually, allowing for complete control over loss characteristics. Two cases are then numerically investigated: a point obstacle located in the vicinity of the boundary, mimicking the sound of the tanpura, and then a parabolic obstacle for the sound synthesis of the sitar.

INTRODUCTION

In many musical instruments, across various cultures, the interaction between a vibrating structure and an obstacle is a key feature leading to amplitude-dependent timbral modification, and is essential in order to replicate the resulting sound. This contact may arise due to the excitation of the instrument [START_REF] Boutillon | Model for piano hammers: Experimental determination and digital simulation[END_REF][START_REF] Chaigne | Numerical modeling of the timpani[END_REF], in which case the exciting mechanism may be considered as lumped, or during its consequent vibrations [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF]. In the latter case, the contact can be pointwise (e.g. the case of string/fret contact in an electric bass) or distributed (e.g. string/bridge contact in a sitar). Such interactions are strongly nonlinear, which complicates significantly its numerical study.

Following analytical studies of the contact between a string and a rigid obstacle [START_REF] Cabannes | Cordes vibrantes avec obstacles[END_REF][START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF][START_REF] Burridge | The sitar string, a vibrating string with a one-sided inelastic constraint[END_REF]], a number of numerical methods have been developed in order to simulate the interaction between a vibrating string and an obstacle. Waveguides are used in [START_REF] Rank | A waveguide model for slapbass synthesis[END_REF][START_REF] Evangelista | Player-instrument interaction models for digital waveguide synthesis of guitar: Touch and collisions[END_REF][START_REF] Kartofelev | Modeling a vibrating string terminated against a bridge with arbitrary geometry[END_REF][START_REF] Kartofelev | Pitch glide effect induced by a nonlinear string-barrier interaction[END_REF], coupled with finite difference schemes in [START_REF] Krishnaswamy | Methods for simulating string collisions with rigid spatial objects[END_REF], where the string is ideal, and [START_REF] Siddiq | A physical model of the nonlinear sitar string[END_REF], where the string is damped and stiff and the obstacle is located at one end of the string. A modal description is presented in [START_REF] Vyasarayani | Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar[END_REF] for modelling an ideal string vibrating against a parabolic obstacle at one boundary and in [START_REF] Valette | Mécanique de la corde vibrante[END_REF] for a dispersive lossy string where an obstacle consolidated at the bridge of a tanpura is considered. This instrument is also modelled in [START_REF] Chatziioannou | Numerical simulation of tanpura string vibrations[END_REF][START_REF] Chatziioannou | Energy conserving schemes for the simulation of musical instrument contact dynamics[END_REF], where the motion of a stiff damped string against an obstacle is obtained by discretising Hamilton's equations of motion. Finally, a finite difference method is developed in [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF], also allowing the simulation of a stiff, damped string with an obstacle having an arbitrary shape. The case of the interaction between a string and a fretboard is in particular described in [START_REF] Bilbao | Numerical simulation of string/barrier collisions: the fretboard[END_REF]. This obstacle is also considered in [START_REF] Trautmann | Multirate simulations of string vibrations including nonlinear fret-string interactions using the functional transformation method[END_REF] where the Functional Transformation Method is used. However, in these models, eigenfrequencies and damping parameters cannot be arbitrary, but must follow a distribution specified by a small number of tuning parameters.

In this study, we present a conservative numerical scheme to model a stiff damped string vibrating against an arbitrarily shaped obstacle. The key features of the scheme are as follows. A modal expansion is used as a starting point for the linear case (i.e., in the absence of contact). By using an equal number of modes and discretization points, a linear transformation relates the spatial displacement and the modal coordinates, so that the contact force is treated directly with the displacement. Finally a regularized contact force together with an energy-conserving time-stepping scheme are implemented. Eigenfrequencies and damping parameters of the string can be adjusted at ease, and in particular according to experimental measurements so that sound synthesis can be more realistic.

Numerical results of the scheme are illustrated by considering two different obstacles for synthesizing the sounds produced by a string vibrating against a point obstacle and a distributed obstacle, mimicking the bridge of the tanpura and a flat bridge respectively. Sound examples are available at www.lam.jussieu.fr/Membres/Issanchou/Sounds_DAFx16.html.

MODEL SYSTEM

Consider a stiff string of length L (m), tension T (N/m), and with linear mass density µ (kg/m). Stiffness effects in the string are characterised by Young's modulus E (Pa ) of the material, and the moment of inertia I = πr 4 /4, where r is the string radius, in m. The string vibrates in the presence of an obstacle assumed not in contact with the string at rest and described by a fixed profile g(x), x ∈ [0, L]. See Figure 1. Under the assumption of small displacements, the dynamics of the string is described by the following equation of motion:

µutt(x, t) -T uxx(x, t) + EIuxxxx(x, t) = f (x, t), (1)
where u(x, t) is the transverse displacement of the string in a single polarisation perpendicular to the barrier. Partial differentiation with respect to time t and coordinate x are indicated by multiple subscripts. Simply supported boundary conditions at the string endpoints are assumed:

u(0, t) = u(L, t) = uxx(0, t) = uxx(L, t) = 0, ∀t ∈ R + . ( 2
)
No damping is considered as yet, and a detailed model of loss will be introduced once modal analysis has been perfomed. See Section 3.1.

u(x,t) z y x 0 L g(x)

Figure 1: A string of length L vibrating against an obstacle g(x). f (x, t) represents the contact force of the barrier upon the string. In general, it should be positive at times and locations along the string for which the string and barrier are in contact. A power-law expression is selected for the contact force. This choice has already been shown to be successful in the realm of musical acoustics [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF][START_REF] Bilbao | Numerical simulation of string/barrier collisions: the fretboard[END_REF][START_REF] Chatziioannou | Energy conserving schemes for the simulation of musical instrument contact dynamics[END_REF][START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF][START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF]. The interaction force thus reads:

f (x, t) = f (η(x, t)) = K [η(x, t)] α + , (3)
where α ≥ 1 is a constant value, η(x, t) = g(x) -u(x, t) is a measure of interpenetration of the string into the barrier, and

[η]+ = 0.5 (η + |η|) is the positive part of η [3]
. Such a powerlaw model of a repelling force due to a collision may be interpreted as a compression of the medium in the case of flexible objects (such as, e.g., a piano hammer [START_REF] Chaigne | Numerical simulations of piano strings. i. a physical model for a struck string using finite difference methods[END_REF]). In the present case, it is best thought of as a regularized (smooth) force penalising such interpenetration [START_REF] Harmon | Robust, efficient, and accurate contact algorithms[END_REF]. As such, the constant K ideally will take on a very large value. This regularized approach contrasts with nonsmooth methods for which no penetration is allowed [START_REF] Banerjee | Historical origin and recent development on normal directional impact models for rigid body contact simulation: A critical review[END_REF][START_REF] Brogliato | Numerical Methods for Nonsmooth Dynamical Systems[END_REF].

One particular advantage of the present choice relies in the fact that the force f derives from a potential ψ:

f = dψ dη , with ψ(η) = K α + 1 [η] α+1 + . ( 4
)
This property is of special interest in the design of energy-conserving numerical methods. See Section 3.3.

Energy Balance

Energy techniques play an important role in the construction of numerical methods for highly nonlinear systems, as in the present case of the string in contact with a barrier. For instance, it is used in the present case of distributed collisions in [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF], in the case of finite difference schemes. The basic strategy is to associate, with a given numerical method, a numerical conserved or dissipated energy quantity, which is itself a positive semi-definite function of the state. As such, it can be used to bound the dynamics of the system, and to find sufficient numerical stability conditions. The continuous energy expression associated with (1) is obtained by multiplying (1) by ut and then employing integration by parts over the spatial domain. It may be written as:

H = L 0 µ 2 (ut) 2 + T 2 (ux) 2 + EI 2 (uxx) 2 + ψ dx. ( 5
)
It satisfies H ≥ 0 and the following equality:

dH dt = 0, (6)
implying that energy is conserved. The first three terms in the expression correspond to stored energy due to the effects of inertia, tension and stiffness, respectively. The final term denotes the energy stored in the collision mechanism. Note that, as losses have not yet been introduced, the system is Hamiltonian. When losses are introduced, one should expect a balance of the form

dH dt + Q = 0, (7)
for some function Q(t) ≥ 0, with the interpretation of power loss, implying that

0 ≤ H (t) ≤ H (0) (8)
for t ≥ 0. It is not easy to give a simple expression for Q in the case of realistic models of loss in strings, which are usually expressed in the frequency domain [START_REF] Valette | Mécanique de la corde vibrante[END_REF], and not in terms of a spatiotemporal PDE system. Thus our expression for power loss is postponed until a modal analysis has been carried out. See Section 3.1.

NUMERICAL SCHEME FOR A STRING VIBRATING AGAINST AN OBSTACLE

The main characteristics of the numerical scheme are the following:

1. An exact scheme for a lossy linear oscillator without obstacle, as described in in [START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF], is used as a building block.

2. Each mode of the string can be described, in isolation, with such an oscillator. Therefore we apply the exact scheme to each mode, ensuring a fine description of frequencies and losses, adjusted for each mode. Then we add a force term F (corresponding to a modal representation of f in equation ( 1)). At this point, we obtain an equation in terms the modal coefficients of u.

3. Taking as many modes as interior points of the spatial mesh, we can rewrite the equation on u directly, through Fourier transformation. Then the force term is expressed as in [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF], in order to obtain a conservative scheme.

Modal analysis

The modal expansion for the displacement of the string is as follows:

u(x, t) = ∞ j=1 qj(t)φj(x), with φj(x) = 2 L sin jπx L (9)
for simply supported boundary conditions.

Inserting the expansion of u in (1), one obtains:

µ(q + Ω 2 q) = F, ( 10
)
where q is a vector containing modal coefficients, and Ω is a diagonal matrix such that Ωj,j = ωj = 2πνj.

Eigenfrequencies are given by νj = j c Equation ( 10) describes a lossless string, where the linear part corresponds to the description of a lossless oscillator for each mode. Therefore, losses can now be introduced by associating each mode with a lossy oscillator. Then [START_REF] Kartofelev | Pitch glide effect induced by a nonlinear string-barrier interaction[END_REF] becomes [START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF]:

µ(q + Ω 2 q + 2Σ q) = F, ( 11
)
where Σ is a diagonal matrix such that Σj,j = σj. A damping parameter σj is now associated to each modal equation, and can be tuned at ease in order to consider any frequency dependence.

Let us now introduce the theoretical model for losses proposed in [START_REF] Valette | Mécanique de la corde vibrante[END_REF], which will allow us to determine realistic values to damping parameters in [START_REF] Krishnaswamy | Methods for simulating string collisions with rigid spatial objects[END_REF]. This model takes into account the three main dissipation mechanisms in musical strings, namely friction with surrounding air, viscoelastic and thermoelastic behaviour of the material as internal losses. The following expression of the quality factor Qj = πν j σ j has been given as:

Q -1 j = Q -1 j,air + Q -1 j,ve + Q -1 te , (12)
where subscripts ve and te refer to viscoelastic and thermoelastic losses, and:

Q -1 j,air = R 2πµνj , R = 2πηair + 2πd √ πηairρairνj, Q -1 j,ve = 4π 2 µEIδve T 2 ν 2 j .
In these expressions, ηair and ρair are, respectively, the dynamic viscosity coefficient and the air density. In the rest of the paper, they are set to the following values, assumed constant: ηair = × 10 -5 kg m -1 s -1 and ρair = 1.2 kg m -3 . Finally this loss model depends on two parameters that can be fitted from e.g. experimental measurements, as performed for example in [START_REF] Valette | Mécanique de la corde vibrante[END_REF][START_REF] Paté | Predicting the decay time of solid body electric guitar tones[END_REF][START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF]: the viscoelastic losses angle δve, and the constant value Q -1 te characterizing the thermoelastic damping. It results in a frequency-dependent loss model which accurately takes into account the damping mechanism present in musical string vibrations [START_REF] Valette | Mécanique de la corde vibrante[END_REF][START_REF] Paté | Predicting the decay time of solid body electric guitar tones[END_REF][START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF].

Spatial discretization

The spatial discretization is defined as xi = iL N , i ∈ {0, ..., N }. Boundary conditions are u(x0, t) = 0 and u(xN , t) = 0 ∀t ∈ R + . In the following, only the values of u on the grid with i ∈ {1, 2, ..., N -1} are thus needed in any calculation.

Considering N -1 modes, the following relationship is fulfilled, ∀i ∈ {1, 2, ..., N -1}:

ui(t) = N -1 j=1 qj(t)φj(xi) = N -1 j=1 qj(t) 2 L sin jπi N . (13)
This can be written in matrix form as u = Sq, where Si,j = φj(xi), ∀(i, j) ∈ {1, ..., N -1} 2 . The inverse of S can easily be calculated: S -1 = L N S T . This linear relation between u and q will be useful in the following analysis.

Time discretization

Some notations are first introduced:

u n i = u(n∆t, xi) δt-u n = u n -u n-1 ∆t δt+u n = u n+1 -u n ∆t δt.u n = u n+1 -u n-1 2∆t δttu n = u n+1 -2u n + u n-1 ∆ 2 u, v = ∆x j∈{1,...,N -1} ujvj.
Let us now consider the following time discretization of [START_REF] Krishnaswamy | Methods for simulating string collisions with rigid spatial objects[END_REF]:

µ ∆t 2 (q n+1 -Cq n + Cq n-1 ) = F n , (14)
where C, C are diagonal matrices such that :

Ci,i = e -σ i ∆t e √ σ 2 i -ω 2 i ∆t + e - √ σ 2 i -ω 2 i ∆t
Ci,i = e -2σ i ∆t .

When the contact force is not present, the modal approach may be viewed as an assembly of independent linear oscillators. This temporal integration scheme gives an exact solution in this case, as shown for example in [START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF], ensuring that at least the linear part is well-approximated (indeed, perfectly). The contact force, through which modes are coupled, remains to be determined. In order to avoid the difficulty linked to this coupling, the dynamical equation for the modal displacements vector q is rewritten for the physical displacement vector u. Thanks to the linear relationship stated above between u and the modal displacement q, the discrete equation on u can be written as:

µ ∆t 2 (u n+1 -Du n + Du n-1 ) = f n , (15)
where D = SCS -1 and D = S CS -1 . The contact force expression used in the present approach is the same as in [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF]:

f n = - δ t-ψ n+ 1 2 δ t. u n , where ψ n+ 1 2 = 1 2 (ψ n+1 + ψ n ) and ψ n = ψ(η n
). This formulation allows for a conservative scheme when there is no loss, and a dissipative one otherwise. Therefore, at each time step, the following equation must be solved:

r + b + m ψ(r + a) -ψ(a) r = 0, (16)
where

r = u n+1 -u n-1 is the unknown, a = u n-1 , m = ∆t 2 µ and b = -Du n + Du n-1 + u n-1 .
The Newton-Raphson algorithm may be used to this end [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF].

Compared to a finite difference approach, the major advantage is the consideration of damping parameters and frequencies that can be adjusted mode by mode. However, the matrices D, D are full in the modal case, as opposed to sparse in the case of local finite difference approximations, and computation time increases accordingly (see section 4.2).

Numerical energy analysis and stability conditions for the scheme are provided in the next section.

Stability analysis

The numerical scheme [START_REF] Chatziioannou | Numerical simulation of tanpura string vibrations[END_REF] is in a form close to actual implementation. However, in order to derive the discrete energy associated to the scheme, it is more convenient to rewrite it in terms of discrete temporal operators. Following the given in the case of an oscillator in [START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF] to get an equivalent scheme, one obtains:

µ Č1δttq n + Č2q n + Č3δt.q n = F n , (17)
with Č1, Č2 and Č3 diagonal matrices satisfying:

Č1 ii = 1 + (1 -γi) ω 2 i ∆t 2 2 1 + (1 -γi) ω 2 i ∆t 2 2 + σ * i ∆t Č2 ii = ω 2 i 1 + (1 -γi) ω 2 i ∆t 2 2 + σ * i ∆t Č3 ii = 2σ * i 1 + (1 -γi) ω 2 i ∆t 2 2 + σ * i ∆t
.

The coefficients γi and σ * i may be written as

γi = 2 ω 2 i ∆t 2 - Ai 1 + ei -Ai σ * i = 1 ∆t + ω 2 i ∆t 2 -γi ω 2 i ∆t 2 1 -ei 1 + ei where Ai = e -σ i ∆t e √ σ 2 i -ω 2 i ∆t + e - √ σ 2 i -ω 2 i ∆t
and ei = e -2σ i ∆t . (18) The scheme for the displacement u thus may be written as:

µ Ď1δttu n + Ď2u n + Ď3δt.u n = f n , (19)
where Ď1 = S Č1S -1 , Ď2 = S Č2S -1 and Ď3 = S Č3S -1 . The force term is expressed as previously.

A discrete energy balance can be obtained by taking the inner product between equation [START_REF] Bilbao | Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics[END_REF] and δt.u n :

δt-H n+ 1 2 = -µ δt.u n , Ď3δt.u n . ( 20
)
where :

H n+ 1 2 = µ 2 δt+u n , Ď1δt+u n + µ 2 u n+1 , Ď2u n + ψ n+ 1 2 , 1 . (21
) Ď3 is positive semi-definite, so that the scheme is strictly dissipative. Therefore it is stable if the energy is positive.

Since the force potential is non-negative, the stability condition is given by:

δt+u n , ( Ď1 - ∆t 2 4
Ď2)δt+u n ≥ 0. ( 22
)
It is therefore sufficient to have

( Ď1 -∆t 2 4 Ď2) positive semi- definite, which is true if ( Č1 -∆t 2 4 Č2
) is positive semi-definite. Consequently, the condition that must be satisfied may be written as:

1 + ei + Ai 1 + ei -Ai ≥ 0 ∀i. (23)
Equation ( 23) is satisfied if 1 + ei ± Ai > 0.This is always true, and hence the scheme is unconditionally stable. The same conclusion is obtained in the limiting case σi = 0 ∀i, using the same reasoning leading to 22 with a reduced expression of γi.

Contact losses

Nonlinear losses due to contact can be added in the presented framework, following the considerations developed in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF][START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF]. The contact force given by (3) may be modified as:

f = dψ dη - du dt Kβ[η] α + , (24)
with β ≥ 0.

The energy (5) of the system (with no loss inherent to the string) then satisfies [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF]:

dH dt = -Qcontact (25)
where

Qcontact = L 0 (ut) 2 Kβ[η] α + dx. (26)
The additional damping term may be discretised using the following expression [START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF]: δt.u n Kβ[η n ] α + . Instead of ( 16), the equation to be solved at each time step is then:

(1 + c)r + b + m ψ(r + a) -ψ(a) r = 0, (27)
where c = ∆t 2µ Kβ[g -u n ] α + . The discrete energy balance is given by:

δt-H n+ 1 2 = -µ δt.u n , Ď3δt.u n -δt.u n , δt.u n Kβ[η n ] α + . (28
) Since δt.u n , δt.u n Kβ[η n ] α + ≥ 0, the dissipation in the system is then increased.

APPLICATION TO MUSICAL INSTRUMENTS AND SOUND SYNTHESIS

In this section, the numerical scheme presented previously is used to simulate the motion of a string vibrating against an obstacle. The cases of a point obstacle and a distributed obstacle are considered.

The string to be considered here is of length L = 1.002 m, under tension T = 180.5 N and with linear mass density µ = 1.17 × 10 -3 kg/m. Musical strings are known to have a very small stiffness. Consequently, the inharmonicity coefficient B introduced in Section 3.1 is chosen to be B = 1.78×10 -5 , corresponding to measured values according to the protocol described in [START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF]. In order to show the ability of the model to incorporate a complex damping law, the theoretical loss model introduced in section 3.1 is used. The values δve = 0.0045 and Q -1 te = 0.000203 have been chosen, and correspond to experimental data, obtained using the method exposed in [START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF]. This model has a complex frequency dependence and accounts for the main loss mechanisms present in strings.

First, a time step, or sample rate has to be chosen. To this end, a convergence study has been undertaken, as elaborated in the following section.

Convergence study

Contact problems lead to the spontaneous generation of high frequencies in the system, due to the very short timescales resulting from collisions between very stiff objects, as in the present case. This implies, at least in theory, a need for a high sampling rate in order to obtain accurate results, which should decrease as the contact stiffness does. The presence of damping, particularly at high frequencies, is expected to be an ameliorating factor. In order to highlight such accuracy issues, as they relate to the time step, convergence results are presented here.

Consider a point obstacle located at x = 6 mm (see Section 4.3). In Figure 2, time history and spectrograms of simulation output are shown, where output is taken as the velocity at the bridge. For the contact restoring force, K = 10 13 and α = 1.3 in order to obtain a stiff contact, limiting penetration of the string in the obstacle to the range of 1 × 10 -8 m. With K = 10 9 for instance, the penetration attains about 1 × 10 -5 m and the resulting sound richness is significantly altered. The string is initially plucked at x = L/2, with a maximal displacement of 1.8 mm. Figure 2 clearly highlights that for a sample rate Fs = 1/∆t = 128 kHz, the simulated sound is far from convergence (it is substantially different from the result with Fs = 1024 kHz), and thus not reliable. From Fs = 256 kHz, the general shape of the time history, together with the spectral dynamics, seems to be correctly reproduced. However, detailed view inside the signal and auditory comparisons definitely evidenced that the high frequencies com-ponents are not well reproduced, such that this sample rate is still insufficient. Finally, careful inspection crossed with listening tests shows that a sample rate of at least Fs = 1 MHz is necessary for convergence of numerical results. To be confident with the accuracy of the simulations presented below, the sampling rate has been fixed for all the simulations to Fs = 2 MHz. Note also that careful comparisons with experimental results have been presented in [START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF], for which a sample rate of 2 MHz was also necessary to obtain very satisfactory results over long simulation times.

Computation cost

In this part, computation times of the numerical scheme are discussed. The computations have been realized with MATLAB on a single CPU with a clock at 2.4 GHz, and time costs are given in table 1 for the computation of one second of sound, rounded up to the nearest minute when N = 1001. Steps which are controlled are the Newton-Raphson loop, the computation of b, a, r n in ( 16), the computation of the energy given in [START_REF] Chaigne | Numerical simulations of piano strings. i. a physical model for a struck string using finite difference methods[END_REF] and finally the total time is given. It clearly appears that the most costly steps are the computation of b, a, r n and the energy, i.e. products of matrices and vectors, most probably because involved matrices are full. Therefore, the computation time is mostly driven by N , so that a judicious way of reducing computation cost may be to adapt the spatial grid to the obstacle, making the space grid finer around the obstacle and larger elsewhere.

N -1 =

Point obstacle: case of the tanpura

The tanpura is an Indian instrument which is played by plucking open strings. The strings are connected to a curved bridge over which a thread is carefully installed makes the sound very specific to this instrument. This bridge and its thread (fully modelled in [START_REF] Chatziioannou | Numerical simulation of tanpura string vibrations[END_REF]) can mostly be described as a two point bridge [START_REF] Valette | Mécanique de la corde vibrante[END_REF]. Therefore in this part, this two point bridge is considered: assuming that the bridge is at x = 0, a point obstacle is located at x = 6 mm. Two initial conditions are considered, by plucking the string at x = L/5 and x = 4L/5, and the velocity of the string is computed at the bridge. Parameters of the contact force are K = 10 13 and α = 1.3 as in Section 4.1, and the spatial grid is such that N = 1001. At the initial time a smoothed triangle with a maximal amplitude of about 1 mm is imposed as the initial displacement, with no velocity. This small value has been chosen in order to show that even for very small displacements, the contact can introduce strong nonlinear effects which are clearly audible and visible on signals.

When there is no obstacle, one can observe rejection of modes which are not excited by the initial condition, as well as a faster damping of high frequencies compared to low frequencies (see Figure 3). When the obstacle is added, all modes are present and energy is transmitted from modes to others. Spectrum tendencies are similar to those encountered in [START_REF] Chatziioannou | Numerical simulation of tanpura string vibrations[END_REF] (numerical results) and [START_REF] Valette | Mécanique de la corde vibrante[END_REF] (experimental results).

The velocity signal for the string with an obstacle is mostly positive (see Figure 4), this is due to the close position of the measure position to the obstacle. Moreover, the effect of the string stiffness implies dispersion which is clearly visible on temporal signals, and constitutes a precursor. When the plucked point is at x = 4L/5 rather than x = L/5, the precursor needs more time to reach the obstacle, therefore it is much more developed when arriving, which explains the high frequencies richness of the temporal signal in Figure 4.

Moreover, according to [START_REF] Valette | Mécanique de la corde vibrante[END_REF], the frequency sliding depends on the plucking position, which is also observed here.

The temporal evolution of energy H n+ 1 2 is presented in Figure 5. It decreases faster when there is an obstacle, probably because energy is transmitted to higher modes, which are more damped. For a plucking point at x = 4L/5, the energy decreases faster than at x = L/5, since more high frequencies are generated by the bridge.

Distributed obstacle

In this part, a distributed parabolic obstacle is considered at one end of the string, that could mimic the case of a sitar. Its shape is as follows:

g(x) = -ax 2 , (29)
where we set a = 0.0065. The length of the obstacle is 19 mm (see Figure 6). Such a value, small compared to realistic ones (see for instance [START_REF] Vyasarayani | Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar[END_REF], where a measure gives a = 0.5778), is chosen in order to make the contact between the string and the full length of the bridge arise. No particular adjustment is made in terms of inclination, however interesting observations can already be made. Numerical parameters K, α and the initial condition are the same as previously. Temporal and spectral evolutions are pre- sented in Figures 7 and8, as well as the energy of the signal. Because of the obstacle, the velocity signal at the bridge as a minimum value close to 0.

Similarly to the case of the point obstacle, a precursor is visible from the first periods and high frequencies are highly excited by the contact, which results in narrow peaks on temporal signals (Figure 7). As mentioned in [START_REF] Siddiq | A physical model of the nonlinear sitar string[END_REF] for experimental results on a complete instrument, the absence of rejection and a descending formant can clearly be observed in Figure 8. Moreover, as in the point obstacle case, the energy decrease is faster than when there is no obstacle, which may be due to the transfer of energy from lower to higher modes, combined with larger damping at high frequencies.

CONCLUSION

A numerical method combining a modal approach and an energyconserving scheme for the case of the string in contact with a barrier has been introduced. Due to the modal approach, the linear parameters of the strings (eigenfrequencies and damping coefficients) can be set independently for each mode, allowing for flexible control over frequency-dependent loss. In particular, linear characteristics of a measured string can be used in order to obtain very realistic sounds, as proposed in [START_REF] Issanchou | Méthode modale mixte pour le contact unilatéral corde / obstacle : application au chevalet de la tampoura[END_REF]. Such complete control over damping rates for the string constitutes one of the major advantages of a modal approach to synthesis. The scheme itself, though with modal characteristics, operates ultimately in the spatial domain, much like a finite difference scheme, though with the special property of being exact under linear conditions. Such a method can thus be viewed as a type of spectral method [START_REF] Trefethen | Spectral Methods in MATLAB[END_REF], accompanied by a time-stepping method which is tuned according to the modal frequencies (unlike, e.g., more typical methods such as leap-frog, or members of the Runge-Kutta family). As with spectral methods, though the updates are no longer sparse, as in the case of FD schemes, it may be possible to employ fast transforms (such as a variant of the FFT) in order to perform the updates (which in this case, follows directly from the sinusoidal structure of the matrix S).

It has been shown that a high sample rate is necessary in order to obtain reliable results for simulations over a long duration. Computation time, which does not allow a real-time simulation at present on a standard machine, could be improved by considering a variable spatial step, finer in the region surrounding the obstacle, and possibly a variable time step [START_REF] Flores | On the contact detection for contact-impact analysis in multibody systems[END_REF], since high temporal precision is only necessary when the string is in contact with the obstacle. In audio applications, though, such algorithmic refinements must be treated with care-the use of a variable grid spacing will entail a loss of structure in the resulting update matrices, and the use of a variable time step must necessarily be accompanied by some form of on-line sample rate conversion which is perceptually transparent.

The next steps of this research will consider more closely the comparison between the outcomes of these numerical methods and measurements realized on a real string and/or on real instruments. To this end, and also in the interest of higher quality synthesis, a parametric study on the contact force parameters should be carried out, and the incorporation of various additional features of the string/barrier system would be of interest. One such feature is the extension of string vibration to two polarisations, as recently explored in the context of bowed string synthesis [START_REF] Desvages | Two-polarisation finite difference model of bowed strings with nonlinear contact and friction forces[END_REF] and in case of the tanpura [START_REF] Bridges | Investigation of tanpura string vibrations using a two-dimensional time-domain model incorporating coupling and bridge friction[END_REF]. At present, there is not a model of an excitation mechanism, and a refined plucking model could be included, perhaps modeling the dynamics of the player's fingers [START_REF] Chadefaux | A model of harp plucking[END_REF]. Ultimately, a complete instrument will require a model of coupling to the instrument body, acoustic radiation, and possibly sympathetic strings [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF].

2L 1 +

 1 Bj 2 , where B = π 2 EI T L 2 describes the inharmonicity created by taking into account the stiffness of the string. Finally the vector F represents the modal projection of the contact force, with Fj = L 0 f (x, t)φj(x)dx.

Figure 2 :

 2 Figure 2: At left : time-domain output signals. At right : spectrograms. From top to bottom : Fs = 128, 256, 512 and 1024 kHz.

Figure 3 :

 3 Figure 3: Spectral evolution of velocities at the bridge (a) no obstacle, plucking at x = L/5 (b) point obstacle, plucking at x = L/5 (c) point obstacle, plucking at x = 4L/5.

Figure 4 :

 4 Figure 4: Temporal evolution of velocities at the bridge for a point obstacle mimicking the case of the tanpura (a) no obstacle, plucking at x = L/5 (b) point obstacle, plucking at x = L/5 (c) point obstacle, plucking at x = 4L/5

5 Figure 5 :

 55 Figure 5: Energy decay for different configurations.

Figure 6 :

 6 Figure 6: Shape of the distributed obstacle.

Figure 7 :

 7 Figure 7: Temporal evolution of the velocity at the bridge and its energy

Figure 8 :

 8 Figure 8: Spectral evolution of the velocity at the bridge, distributed obstacle

Table 1 :

 1 Computation

	Newton-Raphson	0.7	1.3	2.6	9.8
	b, a, r n in (16)	1.5	3.1	6	31.4
	Energy	0.8	1.7	3.5	16.8
	Total time	3.1	6.4	12.5	59.5
	Newton-Raphson	1	2	4	14
	b, a, r n in (16)	5	7	14	119
	Energy	3	6	12	54
	Total time	10	16	30	189

500 Fs 44.1 kHz 88.2 kHz 176 kHz 1 MHz N -1 = 1001 Fs 44.1 kHz 88.2 kHz 176 kHz 1 MHz times, in minutes, for N -1 = 500 and N -1 = 1001.

ACKNOWLEDGMENTS

This work was supported by the European Research Council, under grant number ERC-2011-StG-279068-NESS.