Anne Benoit

Aurélien Cavelan

Yves Robert

Hongyang Sun

Two-level checkpointing and verifications for linear task graphs

Keywords: resilience, fail-stop errors, silent errors, multilevel checkpoint, verification, dynamic programming

Fail-stop and silent errors are omnipresent on largescale platforms. Efficient resilience techniques must accommodate both error sources. To cope with the double challenge, a twolevel checkpointing and rollback recovery approach can be used, with additional verifications for silent error detection. A failstop error leads to the loss of the whole memory content, hence the obligation to checkpoint on a stable storage (e.g., an external disk). On the contrary, it is possible to use in-memory checkpoints for silent errors, which provide a much smaller checkpointing and recovery overhead. Furthermore, recent detectors offer partial verification mechanisms that are less costly than the guaranteed ones but do not detect all silent errors. In this paper, we show how to combine all of these techniques for HPC applications whose dependency graph forms a linear chain. We present a sophisticated dynamic programming algorithm that returns the optimal solution in polynomial time. Simulation results demonstrate that the combined use of multi-level checkpointing and verifications leads to improved performance compared to the standard single-level checkpointing algorithm.

I. INTRODUCTION

Resilience is one of the major challenges for extreme-scale computing [START_REF] Cappello | Toward Exascale Resilience[END_REF], [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF]. In particular, several types of errors should be considered. In addition to the classical fail-stop errors (such as hardware failures), silent errors, also known as silent data corruptions, constitute another threat that cannot be ignored any longer [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF], [START_REF] O'gorman | The effect of cosmic rays on the soft error rate of a DRAM at ground level[END_REF], [START_REF] Ziegler | Cosmic ray soft error rates of 16-Mb DRAM memory chips[END_REF], [START_REF] Ziegler | IBM experiments in soft fails in computer electronics[END_REF]. In order to deal with both types of errors, a traditional checkpointing and rollback recovery strategy can be used [START_REF] Chandy | Distributed snapshots: Determining global states of distributed systems[END_REF], coupled with a verification mechanism to detect silent errors [START_REF] Benoit | Efficient checkpoint/verification patterns[END_REF], [START_REF] Chen | Online-ABFT: An online algorithm based fault tolerance scheme for soft error detection in iterative methods[END_REF], [START_REF] Sao | Self-stabilizing iterative solvers[END_REF]. Such a verification mechanism can be either general-purpose (e.g., based on replication [START_REF] Fiala | Detection and correction of silent data corruption for large-scale high-performance computing[END_REF] or even triplication [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF]) or application-specific (e.g., based on algorithm-based fault tolerance (ABFT) [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF], on approximate re-execution for ODE and PDE solvers [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF], or on orthogonality checks for Krylovbased sparse solvers [START_REF] Chen | Online-ABFT: An online algorithm based fault tolerance scheme for soft error detection in iterative methods[END_REF], [START_REF] Sao | Self-stabilizing iterative solvers[END_REF]).

Because verification mechanisms can be costly, alternative techniques capable of rapidly detecting silent errors, with the risk of missing some errors, have been recently developed and studied [START_REF] Bautista-Gomez | Which verification for soft error detection?[END_REF], [START_REF] Gomez | Detecting silent data corruption through data dynamic monitoring for scientific applications[END_REF], [START_REF] Berrocal | Lightweight silent data corruption detection based on runtime data analysis for HPC applications[END_REF], [START_REF] Cavelan | Assessing the impact of partial verifications against silent data corruptions[END_REF]. We call these verifications partial verifications, while perfect verifications (with no error missed) are referred to as guaranteed verifications. Furthermore, rather than checkpointing only on a stable storage (e.g., an external disk), a lightweight mechanism of in-memory checkpoints can be provided: one keeps a local copy of the data that has not been corrupted when a silent error strikes, and it can be used to perform a recovery rapidly. However, such local copies are lost if a fail-stop error occurs, and hence copies on stable storage must also be provided.

Designing resilience algorithms by combining all of these techniques is challenging. In this paper, we deal with a simplified, yet realistic, application framework, where a set of application workflows exchange data at the end of their execution. Such a framework can be modeled as a task graph whose dependencies form a linear chain. This scenario corresponds to an HPC application whose workflow is partitioned into a succession of (typically large) tightly-coupled computational kernels, each of which is identified as a task. At the end of each task, we can perform either a partial verification or a guaranteed verification of the task output; or, probably less frequently, we can perform a guaranteed verification followed by a memory checkpoint (we do not take the risk of storing a corrupted checkpoint, hence the need for a guaranteed verification); or again, probably even less frequently, we can perform a guaranteed verification, a memory checkpoint and a disk checkpoint in a row.

The main contribution of this paper is a sophisticated dynamic programming algorithm that returns the optimal solution, i.e., the solution that minimizes the expected execution time, in polynomial time. The originality is that we combine both types of verifications and both types of checkpoints. Furthermore, we present extensive simulations that demonstrate the usefulness of mixing these techniques, and, in particular, we demonstrate the gain obtained thanks to multilevel checkpointing.

To the best of our knowledge, the interplay of verification mechanisms with two types of checkpoints, in-memory and disk-based, has never been investigated for task graphs. Our previous work [START_REF] Benoit | Assessing general-purpose algorithms to cope with fail-stop and silent errors[END_REF] considers linear chains with a single checkpoint type and guaranteed verification (for the record, the pioneering work [START_REF] Toueg | On the optimum checkpoint selection problem[END_REF] for linear chains only deals with a single checkpoint type and no verification). The closest work to this paper is our recent work [START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF] for divisible load applications, where we consider the same combined framework (with two error sources, two checkpoint types and two verification mechanisms); however, in [START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF], we target longlasting executions, which are partitioned into periodic patterns that repeat over time, and we compute the best pattern up to first-order approximations. Here, we do not have the flexibility of divisible load applications, since we can insert resilience mechanisms only at the end of the execution of a task. We may well have a limited number of tasks, which prevents the use of any periodic strategy. Instead, we take a completely different approach and design (quite involved) dynamic programming algorithms that provide the optimal solution for any linear task graph.

The rest of this paper is organized as follows. We detail the model in Section II, before giving the dynamic programming algorithm in Section III. The simulation results are presented in Section IV. Finally, we conclude the paper in Section V with hints for future direction.

II. MODEL

We consider a linear chain of tasks

T 1 → T 2 → • • • → T n ,
where each task T i has a weight w i corresponding to the computational load. For notational convenience, we define W i,j = j k=i+1 w k to be the time to execute tasks T i+1 to T j for any i < j. Furthermore, we assume that hardware faults (fail-stop errors) and silent data corruptions (silent errors) coexist, as motivated in Section I. Since these two types of errors are caused by different sources, we assume that they are independent and that both occurrences follow Poisson process with arrival rates λ f and λ s , respectively. The probability of having at least a fail-stop error during the execution of tasks T i+1 to T j is given by p f i,j = 1 -e -λ f Wi,j and that of having at least a silent error during the same execution is p s i,j = 1 -e -λsWi,j .

To deal with both fail-stop and silent errors, resilience is provided through the use of a two-level checkpointing scheme coupled with an error detection (or verification) mechanism. When a fail-stop error strikes, the computation is interrupted immediately due to a hardware fault, so all the memory content is destroyed: we then recover from the last disk checkpoint or start again at the beginning of the application. On the contrary, when a silent error is detected, either by a partial verification or by a guaranteed one, we roll back to the nearest memory checkpoint, and recover from the memory copy there, which is much cheaper than recovering from the disk checkpoint.

We enforce that a memory checkpoint is always taken immediately before each disk checkpoint. This can be done with little overhead and it has been enforced in some practical multi-level checkpointing systems [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF]. Also, a guaranteed verification is always taken immediately before each memory checkpoint, so that all checkpoints are valid (both memory and disk checkpoints), and hence only one memory checkpoint and one disk checkpoint need to be maintained at any time during the execution of the application. Furthermore, we assume that errors only strike the computations, while verifications, memory copies, and I/O transfers are protected from failures.

Let C D denote the cost of disk checkpointing, C M the cost of memory checkpointing, R D the cost of disk recovery, and R M the cost of memory recovery. Recall that when a disk recovery is done, we also need to restore the memory state. For simplicity, we assume that the cost R M is included in the cost R D . Also, let V * denote the cost of guaranteed verification and V the cost of a partial verification. The partial verification is also characterized by its recall, which is denoted by r and represents the proportion of detected errors over all silent errors that have occurred during the execution. For notational convenience, we define g = 1 -r to be the proportion of undetected errors. Note that the guaranteed verification can be considered as one with recall r * = 1. Since a partial verification usually incurs a much smaller cost and yet has a reasonable recall [START_REF] Gomez | Detecting silent data corruption through data dynamic monitoring for scientific applications[END_REF], [START_REF] Berrocal | Lightweight silent data corruption detection based on runtime data analysis for HPC applications[END_REF], it is highly attractive for detecting silent errors, and we make use of them between guaranteed verifications.

The objective is to determine where to place disk checkpoints, memory checkpoints, guaranteed verifications and partial verifications, in order to minimize the expected execution time (or makespan) of the application.

III. DYNAMIC PROGRAMMING

In this section, we present a sophisticated multi-level dynamic programming algorithm to decide which tasks to verify, which tasks to checkpoint, as well as which type of verification or checkpoint to perform. Recall that we assumed a memory checkpoint always comes with a guaranteed verification to ensure the correctness of the results, and that a disk checkpoint always comes with a memory checkpoint, as motivated in Section II.

For convenience, we add a virtual task T 0 , which is checkpointed on disk (and hence in memory), and whose recovery cost is zero. This accounts for the fact that it is always possible to restart the application from scratch at no extra cost.

We first describe in Section III-A a dynamic programming algorithm for the case where we use only guaranteed verifications, memory checkpoints and disk checkpoints. We then show how to extend this dynamic programming algorithm to include partial verifications in Section III-B.

A. Without partial verifications

Figures 1,2, and 3 illustrate the idea of the algorithm without using partial verifications. The algorithm contains three dynamic programming levels, which are responsible for placing disk checkpoints (Figure 1), memory checkpoints (Figure 2), and guaranteed verifications (Figure 3), respectively. An additional step follows to compute the expected execution time between any two verifications. The following describes each step of the algorithm in detail.

Placing disk checkpoints. The first level focuses on placing disk checkpoints (see Figure 1). Let the function E disk (d 2 ) denote the expected time needed to successfully execute all the tasks from T 1 to T d2 , where task T d2 is verified and checkpointed both on disk and in memory.

In this function, we try all possible locations for the last checkpoint before T d2 . For each possible location d 1 , we call the function recursively on d 1 (to place disk checkpoints before T d1 ), and we add the expected time needed to execute the tasks from T d1+1 to T d2 . This is done through the E mem (d 1 , d 2 ) function, which also decides where to place memory checkpoints, and accounts for the cost of memory checkpoints. The cost of the disk checkpoint C D is finally added after T d2 . Note that a location d 1 = 0 means that no . . . T d1 V * CM CD T d1+1 . . . T m1 V * CM T m1+1 . . .

T 0 V * CM CD T 1 . . . T d1 V * CM CD T d1+1 . . . T d2 V * CM CD . . . E disk (d1) Emem(d1, d2) E disk (d2)
T m2 V * CM . . . Emem(d1, m1) E verif (d1, m1, m2)
Emem(d1, m2) further disk checkpoints are added. In this case, we simply let E disk (0) = 0, which initializes the dynamic program. We can express E disk (d 2 ) as follows:

E disk (d 2 ) = min 0≤d1<d2 E disk (d 1 ) + E mem (d 1 , d 2 ) + C D .
The total expected time needed to execute all the tasks T 1 to T n is given by E disk (n).

Placing memory checkpoints. The second level aims at placing additional memory checkpoints between two disk checkpoints (see Figure 2). The function is first called from the first level between two disk checkpoints, each of which also comes with a memory checkpoint. We define E mem (d 1 , m 2 ) as the expected time needed for successfully executing all the tasks from T d1+1 to T m2 , where there is a disk checkpoint at the end of task T d1 , a memory checkpoint at the end of task T m2 , and no other disk checkpoints. Note that there might be a disk checkpoint after T m2 , for instance when we first call this function, but we do not account for the cost of this disk checkpoint in E mem , only for the cost of the memory checkpoint (the cost of the disk checkpoint is already accounted for in E disk ).

As before, we try all possible locations for the last memory checkpoint between tasks T d1 and T m2 . For each possible location m 1 , we call the function recursively on tasks T d1 to T m1 , and then call the function for the next level, E verif (d 1 , m 1 , m 2 ), which computes the expected time needed to execute the tasks from T m1+1 to T m2 (and decides where to place verifications). Finally, we add the cost of the memory checkpoint C M following T m2 . We can express E mem (d 1 , m 2 ) as follows:

E mem (d 1 , m 2 ) = min d1≤m1<m2 E mem (d 1 , m 1 ) + E verif (d 1 , m 1 , m 2 ) + C M . If m 1 = d 1 ,
there is no extra memory checkpoint between d 1 and m 2 , and therefore we initialize the dynamic program with

E mem (d 1 , d 1 ) = 0.
Placing additional verifications. The third level looks for where to insert additional verifications between two tasks with memory checkpoints (see Figure 3). The function is first called

. . . T d1 V * CM CD . . . T m1 V * CMT m1+1 . . . T v1 V * T v1+1 . . . T v2 V * . . . E verif (d1, m1, v1) E(d1, m1, v1, v2) E verif (d1, m1, v2) Figure 3. Placing verifications.
from the second level between two memory checkpoints, each of which also comes with a verification. Therefore, we define

E verif (d 1 , m 1 , v 2 )
as the expected time needed for successfully executing all the tasks from T m1+1 to T v2 , knowing that the last memory checkpoint is after T m1 , the last disk checkpoint is after T d1 , and there are no checkpoints between T m1+1 and T v2 . Note that E verif (d 1 , m 1 , v 2 ) accounts only for the time required to execute and verify these tasks.

As before, we try all possible locations for the last verification between T m1 and T v2 , and for each possible location v 1 , we call the function recursively on tasks T m1 to T v1 . Furthermore, we add the expected time needed to successfully execute the tasks T v1+1 to T v2 , denoted by

E(d 1 , m 1 , v 1 , v 2 ),
knowing the position of the last disk checkpoint d 1 and the position of the last memory checkpoint m 1 . We express

E verif (d 1 , m 1 , v 2 ) as follows: E verif (d 1 , m 1 , v 2 ) = min m1≤v1<v2 E verif (d 1 , m 1 , v 1 ) + E(d 1 , m 1 , v 1 , v 2 ) . (1)
Again, the case v 1 = m 1 means that no further verifications are added, so we initialize the dynamic program with E verif (d 1 , m 1 , m 1 ) = 0. The verification cost at the end of T v2 is accounted for in the function

E(d 1 , m 1 , v 1 , v 2 ).
Computing the expected execution time between two verifications. Finally, to compute the expected time needed for successfully executing several tasks between two verifications, we need the position of the last disk checkpoint d 1 , the position of the last memory checkpoint m 1 , and the positions of the two verifications v 1 and v 2 .

On the one hand, if a fail-stop error occurs with probability p f v1,v2 , then the execution stops and we must recover from the last disk checkpoint. In this case, we lose T lost v1,v2 time, pay the cost of recovery R D (set to 0 if d 1 = 0), and re-execute the tasks starting from T d1 . The re-execution is done in three steps. First, we call E mem (d 1 , m 1 ) to compute the expected time needed to re-execute the tasks from the last disk checkpoint after T d1 to the last memory checkpoint after T m1 . Then, we call the function E verif (d 1 , m 1 , v 1 ) to account for the time needed to re-execute the tasks between the last memory checkpoint after T m1 to the next verification after T v1 . Finally, we re-execute tasks

T v1+1 to T v2 with E(d 1 , m 1 , v 1 , v 2 ).
On the other hand, with probability 1 -p f v1,v2 , there is no fail-stop error. In this case, we pay W v1,v2 by executing all the tasks from T v1+1 to the next verification after T v2 . Then we add the cost of the guaranteed verification V * . After the verification, there is a probability p s v1,v2 of detecting a silent error, in which case we recover from the last memory checkpoint with a cost R M (set to 0 if m 1 = 0) and only re-execute the tasks from there by calling the function

E verif (d 1 , m 1 , v 1 ) followed by E(d 1 , m 1 , v 1 , v 2 ) as before.
Therefore, we have:

E(d 1 , m 1 , v 1 , v 2 ) = p f v1,v2 T lost v1,v2 + R D + E mem (d 1 , m 1 ) + E verif (d 1 , m 1 , v 1 ) + E(d 1 , m 1 , v 1 , v 2 ) + 1 -p f v1,v2 W v1,v2 + V * + p s v1,v2 R M + E verif (d 1 , m 1 , v 1 ) + E(d 1 , m 1 , v 1 , v 2 ) . (2)
In order to compute the expected execution time, we need to compute T lost v1,v2 , which is the expected time loss due to a fail-stop error occurring during the execution of tasks T v1+1 to T v2 . We can derive:

T lost v1,v2 = ∞ 0 xP(X = x|X < W v1,v2 )dx = 1 P(X < W v1,v2 ) Wv 1 ,v 2 0 xP(X = x)dx ,
where P(X = x) denotes the probability that a fail-stop error strikes at time x. By definition, we have

P(X = x) = λ f e -λ f x and P(X < W v1,v2 ) = 1 -e -λ f Wv 1 ,v 2 .
Integrating by parts, we get:

T lost v1,v2 = 1 λ f - W v1,v2 e λ f Wv 1 ,v 2 -1 . (3)
Now, substituting T lost v1,v2 into Equation ( 2) and simplifying, we obtain:

E(d 1 , m 1 , v 1 , v 2 ) = e λsWv 1 ,v 2 e λ f Wv 1 ,v 2 -1 λ f + V * + e λsWv 1 ,v 2 e λ f Wv 1 ,v 2 -1 (R D + E mem (d 1 , m 1 )) + e (λs+λ f )Wv 1 ,v 2 -1 E verif (d 1 , m 1 , v 1 ) + e λsWv 1 ,v 2 -1 R M . (4)
Complexity. The complexity is dominated by the computation of the table

E verif (d 1 , m 1 , v 2 ), which contains O(n 3
) entries, and each entry depends on at most n other entries that are already computed. All tables are computed in a bottom-up fashion, from the left to the right of the intervals. Hence, the overall complexity of the algorithm is O(n 4 ).

B. With partial verifications

It may be beneficial to further add partial verifications between two guaranteed verifications. The intuitive idea would be to add yet another level to the dynamic programming algorithm, and to replace

E(d 1 , m 1 , v 1 , v 2 ) in Equation (1) by a call to a function E (intuitive) partial (d 1 , m 1 , v 1 , p 2 , v 2 ), with p 2 = v 2 ,
which would compute the expected time needed to execute all the tasks from T v1+1 to T p2 and add further partial verifications (computed from the left to the right).

However, while the dynamic programming approach was rather intuitive without partial verifications, the problem becomes much harder with partial verifications. The main reason is that when computing an interval between two partial verifications, there is a probability g that the error remains undetected after the partial verification. When this happens, we need to account for the time lost executing the following tasks until the error is eventually detected (by the subsequent partial verifications or in the worst case by the guaranteed verification) or until the execution is interrupted by a failstop error. This is only possible if we know the optimal positions of the partial verifications after the interval up to the next guaranteed verification. This requires the dynamic programming algorithm to first compute the values at the right of the current interval, hence progressing the opposite way as what was done so far. Therefore, the function becomes

E partial (d 1 , m 1 , v 1 , p 1 , v 2 ) (
expected time needed to execute all the tasks from T p1+1 to T v2 ), and it tries all positions p 2 for the next partial verification. But then, it also requires to remove some terms that account for re-executed work from the intervals on the left of the current interval (because we do not have this information yet), and to re-inject them later in the computation. In the end, we have quite a complicated algorithm.

Expected lost time in case of silent error. First, we compute

E right (d 1 , m 1 , v 1 , p 1 , v 2 )
, the expected time lost executing the tasks T p1+1 to T v2 , assuming that there was a silent error in this interval. This computation uses p 2 , the optimal position of the verification immediately following p 1 , which is computed with the dynamic programming. Indeed, T p2 may detect the error or not. If the error is detected by T p2 , we lose W p1,p2 +V + R M work, while we use E right (d 1 , m 1 , v 1 , p 2 , v 2 ) if the error remains undetected. Also, we consider fail-stop errors only between T p1+1 and T p2 , because fail-stop errors between T p2+1 and T v2 will be accounted for in

E right (d 1 , m 1 , v 1 , p 2 , v 2 ).
Note that even if we know that there is a silent error in the interval, we may need to recover from a fail-stop error if it strikes before the silent error is detected. Altogether, we have:

E right (d 1 , m 1 , v 1 , p 1 , v 2 ) = p f p1,p2 T lost p1,p2 + R D + E mem (d 1 , m 1 ) + (1 -p f p1,p2 ) W p1,p2 + V + (1 -g)R M + gE right (d 1 , m 1 , v 1 , p 2 , v 2 ) = (1 -e -λ f Wp 1 ,p 2 ) 1 λ f - W p1,p2 e λ f Wp 1 ,p 2 -1 + R D + E mem (d 1 , m 1 ) + e -λ f Wp 1 ,p 2 W p1,p2 + V + (1 -g)R M + gE right (d 1 , m 1 , v 1 , p 2 , v 2 ) . The initialization is E right (d 1 , m 1 , v 1 , v 2 , v 2 ) = R M .
Indeed, in this case, there is no task to execute, and if there was a silent error, it is therefore immediately detected by v 2 (a guaranteed verification), and we just pay R M . Knowing p 2 , we are therefore able to compute all values of E right . We will see later how we use this knowledge in E right . Note that the time to reexecute the tasks after a recovery is omitted here, since it will be accounted for when computing

. . . T d1 V * CM CD . . . T m1 V * CM T m1+1 . . . T v1 V * T v1+1 . . . T p1 V T p1+1 . . . T p2 V T p2+1 . . . T v2 V * . . . E lef t (v1, p1) E -(d1,
E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ),
the expected time needed to successfully execute all the tasks between two partial verifications (from T p1+1 to T p2 ).

Expected time to compute tasks T p1+1 to T p2 . Figure 4 shows all the tasks involved in the computation of an interval consisting of several tasks between two partial verifications at

p 1 and p 2 . Let E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 )
denote the expected time needed to successfully execute all the tasks from T p1+1 to T p2 , knowing that the last disk checkpoint is after T d1 , the last memory checkpoint is after T m1 , the last guaranteed verification is after T v1 , and the next guaranteed verification is after T v2 .

On the one hand, if a fail-stop error occurs with probability p f p1,p2 , then the task stops and we must recover from the last disk checkpoint. We lose T lost p1,p2 time, we pay the cost for the disk recovery R D , and we need to re-execute the tasks starting from T d1 . This is done in three steps: first we call E mem (d 1 , m 1 ) to compute the expected time needed to re-execute the tasks from the last disk checkpoint after T d1 to the last memory checkpoint after T m1 . Then we call the function E verif (d 1 , m 1 , v 1 ) to account for the time needed to re-execute the tasks between the last memory checkpoint after T m1 to the next guaranteed verification after T v1 , and finally we are left with the remaining tasks between T v1+1 and T p1 . Let E lef t (v 1 , p 1 ) denote the expected time needed to re-execute all the tasks from T v1+1 to T p1 . Finally, we can reexecute tasks T v1+1 to T v2 by calling

E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ).
On the other hand, there is a probability (1 -p f p1,p2 ) of having no fail-stop errors. In that case, we execute all the tasks from T p1+1 to the next verification after T p2 and we pay W p1,p2 . Then we add the cost V for the verification. After the partial verification, there is a probability p s p1,p2 of having a silent error. In this case, we pay a recovery from the last memory checkpoint (R M ) and re-executed the tasks from there: we call E verif (d 1 , m 1 , v 1 ), followed by E lef t (v 1 , p 1 ) and E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ). Furthermore, if the error was not detected (with probability g), we use

E right (d 1 , m 1 , v 1 , p 2 , v 2 )
to compute the expected time lost executing the tasks following T p2 , knowing that there is an undetected silent error (as explained earlier). Therefore, we have:

E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) = p f p1,p2 T lost p1,p2 + R D + E mem (d 1 , m 1 ) + E verif (d 1 , m 1 , v 1 ) + E lef t (v 1 , p 1 ) + E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) + (1 -p f p1,p2 ) W p1,p2 + V + p s p1,p2 E verif (d 1 , m 1 , v 1 ) + E lef t (v 1 , p 1 ) + E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) + (1 -g)R M + gE right (d 1 , m 1 , v 1 , p 2 , v 2 ) .
Substituting T lost p1,p2 into the equation above and simplifying, we obtain:

E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) = e λsWp 1 ,p 2 e λ f Wp 1 ,p 2 -1 λ f + V + e λsWp 1 ,p 2 e λ f Wp 1 ,p 2 -1 (R D + E mem (d 1 , m 1 )) + e (λs+λ f )Wp 1 ,p 2 -1 E verif (d 1 , m 1 , v 1 ) + E lef t (v 1 , p 1 ) + e λsWp 1 ,p 2 -1 (1 -g)R M + gE right (d 1 , m 1 , v 1 , p 2 , v 2 ) .
Finally, because we do not know at this point how to compute E lef t (v 1 , p 1 ), we remove the term

e (λs+λ f )Wp 1 ,p 2 -1 E lef t (v 1 , p 1 ) from E(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 )
. This corresponds to the amount of time needed to re-execute all the tasks from T v1+1 to T p1 when there is an error between T p1+1 and T p2 . This time will be added back when computing E partial , as explained below. Therefore, we introduce the modified expression of E, denoted E -, as follows:

E -(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) = e λsWp 1 ,p 2 e λ f Wp 1 ,p 2 -1 λ f + V + e λsWp 1 ,p 2 e λ f Wp 1 ,p 2 -1 (R D + E mem (d 1 , m 1 )) + e (λs+λ f )Wp 1 ,p 2 -1 E verif (d 1 , m 1 , v 1 ) + e λsWp 1 ,p 2 -1 (1 -g)R M + gE right (d 1 , m 1 , v 1 , p 2 , v 2 )
.

Computing E partial (d 1 , m 1 , v 1 , p 1 , v 2 )
, the expected time needed to execute all the tasks from T p1+1 to T v2 (and placing extra partial verifications). Finally, we need to compute E partial and to decide when to use additional partial verifications on tasks that are not yet verified. The function is first called from the third level between two guaranteed verifications, and p 1 is originally set to v 1 . Therefore,

E partial (d 1 , m 1 , v 1 , p 1 , v 2 )
denotes the expected time needed to execute all the tasks from T p1+1 to T v2 , where T p1 is followed by a partial verification (with the exception of the first call) and T v2 is followed by a guaranteed verification, knowing the position of the last disk checkpoint d 1 , the position of the last memory checkpoint m 1 and the position of the last guaranteed verification v 1 .

Contrarily to the expressions derived in Section III-A, note that partial verifications are placed from the left to the right. We use the expression of E -, trying all possible positions p 2 for the partial verification following p 1 , and we account for the fact that tasks between T p1+1 and T p2 may be re-executed several times (because we removed E lef t from E -). In fact, for any number of partial verifications between p 2 and v 2 , we can show that

E -(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) is re-executed e (λs+λ f )Wp 2 ,v 2
times, and hence we obtain:

E partial (d 1 , m 1 , v 1 , p 1 , v 2 ) = min p1<p2≤v2 E -(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) • e (λs+λ f )Wp 2 ,v 2 + E partial (d 1 , m 1 , v 1 , p 2 , v 2 ) .
The initialization case is for p 2 = v 2 . Then, there is no re-execution of the interval induced by errors to the right of p 2 (within an E lef t ), and therefore we compute only once

E -(d 1 , m 1 , v 1 , p 1 , v 2 , v 2 ).
Furthermore, the interval is ended by a guaranteed verification, and therefore we add the corresponding verification cost:

E partial (d 1 , m 1 , v 1 , v 2 , v 2 ) = E -(d 1 , m 1 , v 1 , p 1 , v 2 , v 2 ) + e (λs+λ f )Wp 1 ,v 2 (V * -V ).
Because partial verifications are placed from the left to the right, when implementing the algorithm, we first compute all values of E partial on the right of the interval, which are needed to progress towards the left. This is why we always have the values of the next p 2 when computing E right , which correspond to the minimum value selected by E partial . However, it was not possible to derive the values of the reexecution for the left part of the interval, hence the trick to compute the number of times each interval is re-executed, due to a failure on the right (the term E lef t that is removed from the initial expression of E).

Accounting for re-executions on the left. Finally, let us show that for any number of partial verifications between p 2 and v 2 , the function E -(d 1 , m 1 , v 1 , p 1 , p 2 , v 2 ) re-executes an e (λs+λ f )Wp 2 ,v 2 amount of work. If there are no partial verifications after p 2 , then it is executed once when progressing within the computation, and we also need to account for the e (λs+λ f )Wp 2 ,v 2 -1 re-executed work due to the term E lef t that was suppressed from

E -(d 1 , m 1 , v 1 , p 2 , v 2 , v 2 ).
With one intermediate partial verification p 3 between p 2 and v 2 , the same reasoning shows that there is an amount of e (λs+λ f )Wp 2 ,p 3 re-executed work coming from the initial execution and the E lef t term suppressed from E -(d 1 , m 1 , v 1 , p 2 , p 3 , v 2 ). Furthermore, there is an amount of e (λs+λ f )Wp 3 ,v 2 -1 re-executed work coming from the

E lef t term of E -(d 1 , m 1 , v 1 , p 3 , v 2 , v 2 ).
In turn, this reexecuted work incurs e (λs+λ f )Wp 2 ,p 3 re-executed work (initial execution and re-executions due to the E lef t term coming from E -(d 1 , m 1 , v 1 , p 2 , p 3 , v 2 )). Overall, the number of reexecuted work is finally

e (λs+λ f )Wp 2 ,p 3 + e (λs+λ f )Wp 2 ,p 3 e (λs+λ f )Wp 3 ,v 2 -1 ,
and therefore there is a total amount of e (λs+λ f )Wp

IV. PERFORMANCE EVALUATION

In this section, we conduct a set of simulations to assess the relative efficiency of our approach under realistic scenarios. We instantiate the model with actual parameters from the literature and we compare the performance of three algorithms: (i) a single-level algorithm A DV * with only disk checkpoints (and additional guaranteed verifications), (ii) a two-level algorithm combining memory and disk checkpoints A DM V * (as in Section III-A), and (iii) the complete algorithm using additional partial verifications A DM V (as in Section III-B). The optimal positions of verifications and disk checkpoints can be easily derived for A DV * , using a simplification of the proposed dynamic programming algorithm in Section III-A with no additional memory checkpoints.

Simulation setup. We make several assumptions on the input parameters. First, we assume that recovery costs and checkpoint costs are similar; following [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF], [START_REF] Quaglia | A cost model for selecting checkpoint positions in time warp parallel simulation[END_REF], we set recovery costs to be the same as checkpoint costs, i.e., R D = C D and R M = C M . Then, we assume that a guaranteed verification must check all the data in memory, making its cost in the same order as that of a memory checkpoint, i.e., V * = C M . Furthermore, we assume partial verifications similar to those proposed in [START_REF] Gomez | Detecting silent data corruption through data dynamic monitoring for scientific applications[END_REF], [START_REF] Gomez | Detecting and correcting data corruption in stencil applications through multivariate interpolation[END_REF], [START_REF] Berrocal | Lightweight silent data corruption detection based on runtime data analysis for HPC applications[END_REF], with very low cost while offering good recalls. In the following, we set V = V * /100 and r = 0.8. The total computational weight is set to be 25000 seconds and it is distributed uniformly among up to n = 50 tasks in three different patterns shown as follows.

(1) Uniform: all tasks share the same weight W/n, as in matrix multiplication or in some iterative stencil kernels.

(2) Decrease: task T i has weight α(n + 1 -i) 2 , where α ≈ 3W/n 3 ; this quadratically decreasing function resembles some dense matrix solvers, e.g., by using LU or QR factorization.

(3) HighLow: a set of tasks with large weight is followed by a set of tasks with small weight. In the simulations, we set 10% of the tasks to be large and let them contain 60% of the total computational weight.

We point out that all these choices are somewhat arbitrary and can easily be modified in the evaluations; we believe they represent reasonable values for current and next-generation HPC applications. The code is publicly available at http: //graal.ens-lyon.fr/ ∼ yrobert/chain2levels.zip for the interested readers to experiment with their own parameters.

Platform settings. Table I presents the four platforms used in the simulations and their main parameters. These platforms have been used to evaluate the Scalable Checkpoint/Restart (SCR) library by Moody et al. [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF], who provide accurate measurements for λ f , λ s , C D and C M using real applications. Note that the Hera platform has the worst error rates, with a platform MTBF of 12.2 days for fail-stop errors and 3.4 days for silent errors. In comparison, and despite its higher number of nodes, the Coastal platform features a platform MTBF of 28.8 days for fail-stop errors and 5.8 days for silent errors. In addition, the last platform uses SSD technology for memory checkpointing, which provides more data space, at the cost of higher checkpointing costs. Impact of number of tasks. The first column of Figure 5 presents, for each platform, the normalized makespan with respect to the error-free execution time for different numbers of tasks with the Uniform weight distribution. Note that varying the number of tasks has an impact on both the size of the tasks and the maximum number of checkpoints and verifications an algorithm can place. When the number of tasks is small (e.g., less than 5), the probability of having an error during the execution (either a fail-stop or a silent) increases quickly (more than 10% on Hera) for a single task. As a result, the application experiences more recoveries and re-executions with larger tasks, which increases the execution overhead. However, when the number of tasks is large enough, the size of the tasks becomes small and the probability of having an error during the execution of one task drops significantly, reducing the recovery and re-execution costs at the same time.

Single-level algorithm A DV * . The second column of Figure 5 shows the numbers of disk checkpoints (with associated memory checkpoints) and guaranteed verifications used by the A DV * algorithm on the four platforms and for different numbers of tasks. We observe that a large number of guaranteed verifications is placed by the algorithm while the number of checkpoints remains relatively small (i.e., less than 5 for all the platforms). This is because checkpoints are costly, and verifications help to reduce the amount of time lost due to silent errors. Since verifications are cheaper, the algorithm tends to place as many of them as possible, except when their relative costs also become high (e.g., on Coastal SSD).

Two-level algorithm A DM V * . The third column of Figure 5 presents the numbers of disk checkpoints, memory checkpoints and guaranteed verifications used by the A DM V * algorithm on the four platforms and for different number of tasks. We observe that the number of guaranteed verifications remains similar to that placed by the A DV * algorithm. However, the two-level algorithm uses additional memory checkpoints, which drastically reduces the amount of time lost in reexecution when a silent error is detected. In particular, we observe that the algorithm A DM V * always leads to a better makespan compared to the single-level algorithm A DV * , with an improvement of 2% on Hera and 5% on Atlas, as shown in the first column of Figure 5. This demonstrates the usefulness of the multi-level checkpointing approach.

Two-level algorithm with partial verifications A DM V . The last column of Figure 5 presents the numbers of disk checkpoints, memory checkpoints, guaranteed verifications and additional partial verifications used by the A DM V algorithm on the four platforms and for different numbers of tasks.

Although partial verifications are always more cost-effective than guaranteed ones, due to the imperfect recall, they are only worth it if one can use a lot of them, which is only possible when the number of tasks is large enough. Therefore, the algorithm only starts to use partial verifications when the number of tasks is greater than 30 on Hera, 40 on Coastal and 50 on Atlas, where the silent error rate is the highest among the four platforms. In our setting, adding partial verifications has a limited impact on the makespan, with the exception of the Coastal SSD platform, where the cost of checkpoints and verifications are much higher than on the other platforms. Partial verifications, being 100 times cheaper than guaranteed verifications, remain the only affordable resilience tool on this platform. In this case, we observe an improved makespan (around 1% with 50 tasks) compared to the A DM V * algorithm, as shown in the first column of Figure 5.

Distribution of checkpoints and verifications. Figure 6 shows the positions of the disk checkpoints, memory checkpoints, verifications and partial verifications obtained by running the A DM V algorithm on each of the four platforms and for 50 tasks with the uniform distribution. For all platforms, the algorithm does not perform any additional disk checkpoints. These being costly, the algorithm rather uses more memory checkpoints and verifications. On most platforms, the optimal solution is a combination of equi-spaced memory checkpoints and guaranteed verifications, with additional partial verifications in-between. However, on the Coastal SSD platform, the cost of checkpoints and verifications is substantially higher, which leads the algorithm to choose partial verifications over guaranteed ones.

Decrease pattern. In the following, we focus on the platforms Hera and Coastal SSD, which represent both extremes in terms of size (number of processors) and hardware used for memory checkpointing (RAM and SSD, respectively). The first column of Figure 7 presents the performance of the three algorithms for different number of tasks and for the Decrease pattern. The second column shows the number of disk checkpoints, memory checkpoints, guaranteed and partial verifications given by the A DM V algorithm. The third column is a visual representation of the corresponding solution obtained for 50 tasks and with the same configuration. We observe that the makespan obtained is very similar for all three algorithms (with a slight advantage for A DM V ). Since the large tasks at the beginning of the chain are more likely to fail, they will be checkpointed more often, as opposed to the small tasks at the end, which the algorithm does not even consider worth verifying.

HighLow pattern. Once again, we focus on platforms Hera and Coastal SSD. Similarly to Figure 7, Figure 8 assesses the impact of the HighLow pattern on the performance of the three algorithms as well as on the numbers and the positions of checkpoints and verifications. Recall that we set the first 10% of the tasks to contain 60% of the total computational weight, while the rest of the tasks contain the remaining 40%. With 50 tasks and a total computational weight of 25000s, the first 5 tasks have a weight of 3000s each, while the remaining tasks have a weight of around 222s each. Under this configuration, an error occurring during the execution of a large task would cost T lost ≈ 1500s time loss on average for fail-stop errors (see Equation ( 3)) and 3000s for silent errors, plus an additional 3000s time loss for each preceding task that has not been checkpointed. With the MTBF on Hera, a large task will fail with probability 1.3%, as opposed to the probability of 0.096% for small tasks. As a result, the disk checkpoint, which takes 300s, turns out to be still too expensive, but the memory checkpoint, which takes only 15.4s on Hera, becomes mandatory: on average an error will occur way before the total accumulated cost of our preventive memory checkpoints even adds up to the cost of one task. On Coastal SSD, however, the memory checkpoint is still quite expensive, so that only one of the first 5 tasks is marked for verification and memory checkpointing. On both platforms, since the rest of the tasks are small, the solution is similar to the one we observed for the Uniform pattern, except that memory checkpoints and verifications are less frequent.

Summary of results. Overall, we observe that the combined use of disk checkpoints and memory checkpoints allows us to decrease the makespan, for the three task patterns and the four platforms. The use of partial verifications further decreases the makespan, especially on the Coastal SSD platform where the checkpointing costs are high. To give some numbers, our approach saves 2% of execution time on Hera and 5% on Atlas. These percentages may seem small, but they correspond to saving half an hour a day on Hera, and more than one hour a day on Atlas, with no further overhead.

V. CONCLUSION AND FUTURE WORK

In this paper, we proposed a two-level checkpointing scheme to cope with both fail-stop errors and silent data corruptions on large-scale platforms. Although numerous studies have dealt with either error source, few studies have dealt with both, while it is mandatory to address both sources simultaneously at scale. By combining standard disk checkpointing technique with in-memory checkpoints and verification mechanisms (partial or guaranteed), we have designed a multi-level dynamic programming algorithm that computes the optimal solution for a linear application workflow in polynomial time. Simulations based on realistic parameters on several platforms show consistent results, and confirm the benefit of the combined approach. While the most general algorithm has a high complexity of O(n 6 ), where n is the number of tasks, it executes within a few seconds for n = 50 tasks, and therefore can be readily used for real-life linear workflows whose sizes rarely exceed tens of tasks.

One interesting future direction is to assess the usefulness of this approach on general application workflows. The problem gets much more challenging, even in the simplified scenario where each task requires the entire platform to execute. In fact, in this simplified scenario, it is already NP-hard to decide which task to checkpoint in a simple join graph (n -1 source tasks and a common sink task), with only fail-stop errors striking (hence a single level of checkpoint and no verification at all) [START_REF] Aupy | Scheduling computational workflows on failure-prone platforms[END_REF]. Still, heuristics are urgently needed to address the same problem as in this paper, with two error sources, two checkpoint types and two verification mechanisms, if we are to deploy general HPC workflows efficiently at scale.

Figure 1 .

 1 Figure 1. Placing disk checkpoints.

Figure 2 .

 2 Figure 2. Placing memory checkpoints.

Figure 4 .

 4 Figure 4. Placing partial verifications.

Figure 5 .Figure 6 .

 56 Figure 5. Performance of the three algorithms on each platform with the Uniform pattern. Each row corresponds to one platform.

Figure 7 .Figure 8 .

 78 Figure 7. Performance of the three algorithms, and distribution of disk checkpoints, memory checkpoints and verifications (for the A DM V algorithm) on platforms Hera and Coastal SSD with the Decrease pattern.

  2 ,v 2 reexecuted work. It is easy to extend this reasoning to any number of intervals by induction, assuming that it is true for i intermediate partial verifications pv i , . . . , pv 1 , and adding a partial verification pv i+1 between p 2 and pv i . The same reasoning holds.Complexity. Clearly, the complexity is now dominated by the computation of the tableE partial (d 1 , m 1 , v 1 , p 1 , v 2 ), which contains O(n 5) entries, and each entry depends on at most n other entries that are already computed. Hence, the overall complexity of the algorithm is O(n 6 ).

ACKNOWLEDGMENTS

This research was funded in part by the European project SCoRPiO, by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR), and by the PIA ELCI project. Yves Robert is with Institut Universitaire de France. A short version of this paper was presented at PMBS'15 (without proceedings).


Paul Blinde
email: p.l.blinde@tudelft.nl

Dirk Michaelis

Bas Van Oudheusden

Pierre-Elie Weiss

Roland De Kat

Angeliki Laskari

Jin Yong Jeon

Laurent David

Daniel Schanz

Florian Huhn

Roeland De Kat

Young Jin Jeon

Sebastian Gesemann

Matteo Novara

Cameron Mcphaden

Nathan Neeteson

David Rival

Jan Fg Schneiders

Ferry Schrijer

Comparative assessment of PIV-based pressure evaluation techniques applied to a transonic base flow

Keywords: tomographic PIV, particle tracking, pressure determination

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

Introduction

Fluid pressure is d irectly relevant to phenom ena like surface load ing and sound generation (aeroacoustics) and as such an im portant quantity in m any engineering problem s. The quantification of m ean and flu ctuating pressure is therefore of relevance for efficient d esigns.

While surface pressu re can experim entally be d eterm ined w ith pressure transd ucers and pressure-sensitive paint (PSP), the novel approach of PIV-based pressure d eterm ination (van Oud heusd en 2013) offers unique new capabilities. Com pared to the m ore established m easurem ent techniqu es, PIV-based pressure d eterm ination does not require instrum entation or surface preparation of the w ind tunnel m od el. This allow s for pressure d eterm ination in configurations w here such m od ifications are not practical, e.g. very thin (or m em brane-like) airfoils. In m ore regular configurations, it avoid s the installation of large num bers of p ressure transd ucers w hich is a com m on practice to obtain reliable surface load d istributions. Another unique feature of PIV-based pressure d eterm ination is its inherent ability to provide sim ultaneous velocity and pressure d ata in the full flow field , w hich enables a better insight in the relation betw een fluid d ynam ics and surface load s or sound field s.

Given these beneficial features, there has been a long and increasing interest in PIV-based pressure d eterm ination. In recent years, m oreover, the technique has becom e increasingly feasible and appealing, especially d ue to the d evelopm ent of (tim e-resolved ) volum etric d iagnostic capabilities, such as tom ographic PIV [START_REF] Elsinga | Tom ographic particle im age velocim etry[END_REF][START_REF] Scarano | Tom ographic PIV: principles and practice[END_REF] and Lagrangian particle tracking [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle im age d ensities[END_REF]. Many stud ies have ad d ressed d ifferent im plem entations of this technique (van Oud heusd en 2013). With the m aturing of the m ethod , its scale-up from sm all-size research environm ents is currently receiving attention in collaborative European fram ew ork program s like 'AFDAR' and 'N IOPLEX'. The particular w ork d escribed in the current paper reflects a joint effort by partners in 'N IOPLEX', and to the authors' knowled ge establishes the first system atic com parative stud y, in w hich m u ltiple groups are involved that w ork on the sam e test case for pressure d eterm ination (sim ilar as in the "PIV challenges").

Basic w orking principle

In PIV-based pressure d eterm ination, velocity d ata obtained w ith PIV is used to com p ute the local pressure grad ient (∇p) from the m om entum equation (equation 1), w ith pressure being subsequently obtained through spatial integration.

∇p = -ρ D𝐮 Dt (1)
H ere ρ is the d ensity and D𝐮 Dt ⁄ is the m aterial acceleration, w hich is the acceleration of a fluid particle in a Lagrangian perspective, and w hich is obtained from the flow field m easurem ent.

For incom pressible flow the d ensity is a constant, w hile und er com pressible flow cond itions, the d ensity can be elim inated as ind epend ent variable by com bining the m om entum and energy equations [START_REF] Van Oud Heusd En | Evaluation of integral forces and pressure field s from planar velocim etry d ata for incom pressible and com pressible flow s[END_REF]. Also note that by tim e-averaging of equation 1, m ean pressure field s can be obtained from uncorrelated velocity field s (e.g. van Oud heusd en et al.

2007).

The experim ental d eterm ination of the m aterial acceleration from PIV m easurem ents has been subject of extensive research, for w hich d ifferent m ethod s have been proposed . Using tw o or m ore velocity field s closely separated in tim e, the m aterial acceleration can be d eterm ined via a Lagrangian form ulation by tracing im aginary fluid particles over subsequent velocity field s (e.g. [START_REF] Liu | Instantaneous pressure and m aterial acceleration m easurem ents u sing a fourexposure PIV system[END_REF] or via an Eulerian form ulation by separately d eterm ining local tem poral and spatial velocity d erivatives from the structured velocity m easurem ent grid (e.g. [START_REF] Baur | PIV w ith high tem poral resolution for the d eterm ination of local pressure red uctions from coherent turbulence phenom ena. 3rd Int. Work[END_REF]. To alleviate the need for tim e-resolved velocity d ata, m ethod s have been proposed to calculate the m aterial acceleration from velocity d ata at a single tim e instance by em ploying physical m od els (e.g. de Kat and Ganapathisubram ani 2013; [START_REF] Laskari | Full-field pressure from snapshot and tim eresolved volum etric PIV[END_REF][START_REF] Schneid Ers | Pressure estim ation from single-snapshot tom ographic PIV in a turbulent bound ary layer[END_REF]). Alternative to taking velocity field s as a starting point, the m aterial acceleration can also be obtained from the record ed particle im ages, by tracing ind ivid ual particles (PTV) or particle patterns (e.g. FTEE, [START_REF] Jeon | Fluid trajectory evaluation based on an ensem ble-averaged cross-correlation in tim e-resolved PIV[END_REF]). The PTV-based strategy includ es conventional tracking approaches that rely on relatively low particle seed ing concentrations (e.g. [START_REF] Malik | Particle tracking velocim etry in three-d im ensional flow s[END_REF]) as w ell as m ore recently proposed ad vanced approaches that avoid this d raw back (e.g. N ovara and Scarano 2013; [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle im age d ensities[END_REF]. Also for the integration of the pressure grad ient, d ifferent proced ures have been suggested and tested , such as d irect spatial integration (e.g. [START_REF] Baur | PIV w ith high tem poral resolution for the d eterm ination of local pressure red uctions from coherent turbulence phenom ena. 3rd Int. Work[END_REF]Köngeter 1999, Tronchin et al. 2015) and Poisson form ulations (e.g., de

Kat and van Oud heusd en 2012).

Depend ing on the follow ed approach, the perform ance of

D escription of the test case

The test case consists of the transonic flow over an axisym m etric step, w hich form s a sim plified representation of the flow over the base of a launcher d uring its ascent phase. The flow over the launcher base gives rise to flu ctuating sid e load s on the nozzle and quantifying these load s is critical for the efficient d esign of nozzles and surround ing structures. Apart from this ind ustrial relevance, the choice for the test case is further m otivated by its richness in term s of spatial and tem poral flow scales and its challenging nature that is inherent to high-speed three-d im ensional separating-reattaching flow s. It should be rem arked here that not all applied techniques p erform optim al for the flow cond itions und er consid eration and this test case alone therefore does not constitute an exhaustive perform ance assessm ent of the d ifferent m ethod s. These lim itations should be kept in m ind w hen interpreting the results of the stud y.

To have reference d ata available against w hich the pressure results can be valid ated , the test case is based on a nu m erical sim ulation. A PIV experim ent is sim ulated by creating synthetic particle im ages of a virtual m easurem ent volum e. Mod el geom etry (see Figure 3 Synthetic particle im ages are processed using tom ographic PIV algorithm s [START_REF] Elsinga | Tom ographic particle im age velocim etry[END_REF]) as w ell as w ith the PTV algorithm 'Shake-The-Box' (STB, [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle im age d ensities[END_REF]. This enables the use of PIV-based pressure d eterm ination techniqu es that require velocity input d ata on a Cartesian grid , as w ell as PTV-based techniques w hich take scattered d ata available at particle positions as starting point. To be able to assess the im pact of m easurem ent noise, processing is perform ed using both id ealized , noise-free particle im ages ('clean') as w ell as particle im ages to w hich a controlled am ount of artificial noise is ad d ed ('noisy').

The test case com prises sequences of four subsequent im ages as w ell as continuous tim e-

C p = p p ∞ -1 1 2 γM ∞ 2
(2)

T. I. = √ 1 3 ((𝑢 𝑅𝑀𝑆 ′ ) 2 + (𝑣 𝑅𝑀𝑆 ′ ) 2 + (𝑤 𝑅𝑀𝑆 ′ ) 2 ) U ∞ ⁄ (3)

Calculation of particle tracks

The PIV m easurem ent d om ain is rand om ly seed ed w ith particles so that the synthetic PIV im ages contain 0.05 p articles per pixel (ppp). The particles are propagated using an explicit, fourth-ord er Runge-Kutta m ethod in com bination w ith spline interpolation in tim e and natural neighbour interpolation in space. Three tim e integration steps are perform ed per sim ulation sam pling tim e step of 2 µs, so that the Courant num ber is sm aller than 1 and the Courant-Fried richs-Lew y (CFL) cond ition is m et in the full d om ain w ith the exception of a thin region in the d irect vicinity of the corner of the step. To allow inflow of new particles, a narrow region surround ing the m easurem ent volum e is rand om ly seed ed w ith particles after each tim e step.

The instantaneous particle velocity is taken equal to the instantaneous fluid velocity at the particle position, w hich is a valid approach for the great m ajority of the flow field , w here the flow tim escales are sm aller than the typical seed ing particle relaxation tim e (τ p ) of approxim ately 2 µs [START_REF] Ragni | Particle tracer response across shocks m easured by PIV[END_REF]. To assess the nu m erical errors from integration and interpolation, particles are integrated one sam ple tim e step forw ard and su bsequently backw ard in tim e.

Generation of synthetic particle images

The For the PTV processing, particle paths are d eterm ined using the 'Shake-The-box' algorithm (STB, [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle im age d ensities[END_REF]). An ad apted version of the STB strategy (N ovara et al. 2015) has been applied

to the m ulti-pulse d ata. All cases w ere processed using DLR in-house cod es.

For the initiation phase of the STB PTV processing, the particle d istribution in the first four im ages is d eterm ined using Iterative Particle Reconstruction (IPR, Wieneke 2013). Four norm al triangulation iterations are perform ed (m = 4), follow ed by tw o triangulation iterations (n = 2) using a red uced set of three cam eras. Seven so-called shake iterations are executed after each triangulation (k = 7). The allow ed triangulation error is 0.7 px. in the clean case and 1.1 px. in the noisy case. The intensity threshold for particle id entification is set to 100 counts and 40 counts, respectively. Tracks are searched in the resulting particle d istributions for the four tim e steps, using a pred ictor field obtained from a tom ographic PIV evaluation. For the tim e-resolved cases, from the fifth tim e-step on this pred ictor is com pletely replaced by pred ictions based on alread y id entified particle tracks and the num ber of triangulation and shake iterations is red uced (m = 3, n = 2, k = 7). Three passes of STB are cond u cted going forw ard s and backw ard s in tim e to extend any non-id entified parts of already d etected tracks and to find any previously overlooked tracks.

In the end , the algorithm find s over 99% of all available particle tracks for the clean case and around 96 % for the noisy case. A B-spline function is fit through the resulting particle tracks to red uce noise and to obtain Lagrangian velocity and acceleration as d erivatives of a continuous function.

Concerning the m ulti-pulse d ata, an iterative STB approach strategy is applied , w here the sequential application of IPR (m = 4, n = 2, k = 5) and tracking allow s to progressively increase the num ber of reconstructed tracks. The allow ed triangulation error is 0.6 px. and 1.1 px. for the clean and noisy case respectively. The intensity threshold for particle id entification is set to 20 counts for the clean case, w hile an ad aptive threshold is ad opted for the noisy case w here the threshold value is progressively red uced for each of the STB iterations. Three STB iterations are applied for the clean case and four for the noisy case. Approxim ately 99% and 88% of the real four-pulse particle tracks are successfully id entified by STB for the clean and noisy case respectively. A 3 rd -ord er B-spline function w ith a fixed cut-off frequency has been fit through the resulting particle tracks. A m ore d etailed d escription of the processing m ethod is provided by N ovara et al. (2016). m agnitud e as the error estim ate for the particle integration proced ure, suggesting an alm ost perfect reconstruction. For all other cases, the error typically stays below 0.05 voxel. Larger errors occur at the in-and outflow w here not all tracks are im m ed iately id entified and the particle tracks can benefit less from the inform ation from ad jacent tim e instances. Sim ilar features as for the PIV results can be observed but error levels are an ord er of a m agnitud e sm aller (note the d ifferent scale). For a large part this can be explained by the absence of the spatial filtering that occurs in the PIV processing, w hich is especially noticeable in the shear layer. Furtherm ore, it should be realized that the com parison is som ew hat biased in the sense that the current PIV processing is typical for experim entally obtained particle im ages and no m ulti-fram e crosscorrelation techniques have been applied that w ould exploit the ad d itional tim e inform ation available, w hereas the application of the STB algorithm has not yet been experimentally d em onstrated for the present flow case (i.e. a high-speed transonic flow ).

D escriptions of applied

∇ln ( p p ∞ ) = - 1 RT • D𝐮 Dt = - γM ∞ 2 V ∞ 2 + γ-1 2 M ∞ 2 (V ∞ 2 -V 2 ) • D𝐮 Dt (4)
To provid e a bound ary cond ition to the pressure integration, the isentropic pressure ratio (equation 5) is u sed as Dirichlet bound ary cond ition at the top of the d om ain or as m eans of norm alisation.

p p ∞ = (1 + γ -1 2 M ∞ 2 (1 - V 2 V ∞ 2 )) γ γ-1 (5)
Figure 9 investigates the r.m.s. errors introd uced by the assum ption of ad iabatic flow on tem perature (top -left) and of isentropic flow on the pressure (top -right). The bottom figures show that the im pact of those assum ptions on the pressure solution is an ord er of m agnitud e sm aller than the m agnitud e of the pressure fluctuations present in the flow (com pare Figure 4).

For both assessm ents, the sim ulation d ata w ere first interpolated to the PIV grid points. Material accelerations w ere calculated by taking the spatial grad ient of the reference pressure d ata.

Equation 4 w as su bsequently solved by prescribing the isentropic pressure ratio at the top bound ary and using N eum ann bound ary cond itions at all other bound aries.

PIV-based techniques

The follow ing four techniques reconstruct pressure field s from velocity d ata on a Cartesian grid , as obtained by PIV. These techniques inclu d e an Eulerian approach (EUL) and an (iterative) Lagrangian approach (ILAG) w hich each require at least tw o subsequent snapshots as inp ut d ata for d eterm ining the m aterial acceleration (multi-snapshot approaches) and the Taylor 's hypothesis approach (TH ) and the instantaneous vortex-in-cell m ethod (IVIC) w hich each only require a single velocity snapshot as input d ata (single-snap shot approaches).

Eulerian approach (EUL)

The m aterial acceleration is evaluated in a stationary reference fram e as being com posed of tem poral and spatial velocity d erivatives (see e.g. [START_REF] Baur | PIV w ith high tem poral resolution for the d eterm ination of local pressure red uctions from coherent turbulence phenom ena. 3rd Int. Work[END_REF]. In the present im plem entation these d erivatives are evaluated using second -ord er finite d ifference d iscretization except at the bound aries w here first-ord er d iscretization is used instead . The pressure grad ient is integrated by first taking the d ivergence of equation ( 4) to obtain a Poisson equation (see e.g. de Kat and van Oud heusd en 2012). N eu m ann bound ary cond itions are im posed on all sid es except for the top surface of the dom ain, w here the pressure as obtained from the isentropic flow relations is prescribed as Dirichlet bound ary cond ition (see equation 5).

The problem is d iscretised using a second -ord er finite d ifference schem e and the resulting linear

Instantaneous Vortex-in-Cell (IV IC)
In the IVIC m ethod , the m aterial acceleration is approxim ated from a single-snap shot by application of the vorticity transport equation to the velocity m easurem ent at a single tim e instant. In sum m ary, the m ethod first evaluates the vorticity tem poral d erivative by solution of the (inviscid and incom pressible) vorticity transport equation. In a next step, the velocity tem poral d erivative is calculated by solution of a Poisson equation and then u sed to construct the m aterial acceleration. A m ore d etailed d escription of the approach can be found in [START_REF] Schneid Ers | Pressure estim ation from single-snapshot tom ographic PIV in a turbulent bound ary layer[END_REF]. The subsequent integration of the pressure grad ient is perform ed sim ilarly as for the Eulerian approach d escribed above.

PTV-based techniques

The follow ing four techniques reconstruct pressure field s from acceleration d ata along particle tracks, as obtained by PTV. These techniques includ e Flow Fit and VIC+ w hich first interpolate acceleration d ata to a Cartesian grid and subsequently solve the m om entum equation for pressure, as w ell as Flow Fit2 and Voronoi-based integration that d irectly d eterm ine pressure field s from the sparse acceleration d ata.

FlowFit in combination with FFT integration (FF1)

The acceleration d ata obtained at the particle positions are interpolated to a Cartesian grid by 'Flow Fit'; an iterative optim ization approach that generates a quadratic B-splines representation of the acceleration field , using spline coefficients from m inim izing a cost function that enforces sm oothness and penalizes the curl of acceleration (Gesem ann 2015). Integration is perform ed in Fourier space. A constant pressure offset is ad d ed to the solution to yield the absolute pressure field (H u hn et al. 2016).

FlowFit2 (FF2)

Flow fit2 represents the scalar pressure field as 3D cu bic B-spline function and estim ates its param eters by d efining a cost function that is m inim ized on the basis of a w eighted sum of:  Squared d istances betw een "fitted " and "m easured " pressure grad ients at all particle positions  Squared d istances betw een "fitted " and isentropic pressure at the top of the d om ain  Squared curvatures of the pressure function as a form of Tikhonov regularization.

The resulting linear least squares problem is solved iteratively. After this reconstruction, the Bspline curve is sam pled at the d esired points. Vortex-in-Cell+ (V IC+) VIC+ reconstructs instantaneous velocity and its m aterial d erivative on a regular grid using both the instantaneous velocity and m aterial acceleration evaluated from Lagrangian particle trajectory m easurem ents. The technique follow s an iterative proced ure that m inim izes a cost function and penalizes d eviations from the instantaneous velocity and m aterial acceleration m easurem ent. A d etailed d escription of the approach can be found in Schneid ers et al. (2015).

The subsequent integration of the pressure grad ient is perform ed via a Poisson equation (see Eulerian approach). Pressure grad ients are prescribed as N eum ann bound ary cond itions u sed on all sid es of the d om ain. In a next step, the solution is norm alised by setting the m ean pressure ratio at the top plane equ al to the m ean isentropic pressure ratio (see equation ( 5)). The m ethod is im plem ented using the STB output for the tim e-resolved case and the output of an in-house tom ographic PTV for the m ulti-pulse case.

Voronoi-based pressure integration (V OR)

A Lagrangian finite-volum e m ethod is utilized to spatially integrate equation ( 4) to obtain pressure values at the particle positions provid ed by STB. First, a Lagrangian network is constructed on the particle field [START_REF] De | Pressure-field extraction on unstructured flow d ata using a Voronoi tessellation-based netw orking algorithm : a proof-of-principle stud[END_REF], so that the pressure m ay be extracted d irectly in the Lagrangian fram e. N ext, equation ( 4) is d iscretized to the Lagrangian netw ork using a finite-volum e m ethod . N eum ann bound ary cond itions are enforced im plicitly at all sid es exclud ing the top w here a Dirichlet bound ary cond ition is enforced by prescribing the isentropic pressure ratio (see equation 5). The resulting problem is expressed in a system of linear equations, w hich is solved for pressure.

Pressure results

Pressure from PIV-based techniques

This section presents the results from applying pressure reconstruction techniques to PIV velocity d ata. Figure 10 The left figures show that all m ethod s are able to reconstruct the general features of the reference pressure field . The centre and right figures show that all reconstructed pressure field s have the highest error in the shear layer and reattachm ent region w here the sm allest tem poral and spatial flow scales occur and the precision of the PIV m easurem ent is poorest (see Figure 7). The single- snapshot approaches (TH and IVIC) can be seen to prod uce less precise pressure reconstructions than the m ulti-snap shots approaches (EUL and ILAG); an expected result as EUL and ILAG use m ore tem poral information. N evertheless, in the light of practical d ifficulties in obtaining m ultipulse d ata in high-speed flow s, it is encouraging that from this assessm ent it appears to be quite feasible to reconstruct the m ain pressure field features from a single velocity snapshot.

Com paring the tw o multi-snap shot m ethod s show s that Eulerian approach (EUL) lead s to less precise pressure reconstructions than ILAG (com pare right figures). Investigation of the tim eresolved sequence of instantaneous results furtherm ore show s highly fluctuating error field s for the Eulerian approach (EUL), w hereas the error field for ILAG evolves m ore grad ually over tim e. This is attributed to the better use of tem poral inform ation by ILAG, w hich in the present im plem entation for the tim e-resolved fits a 2 nd -ord er polynom ial over 9 consecutive velocity Com parison w ith the PIV-based approaches shows significantly low er error levels for the PTVbased approaches. This is in line w ith the low er velocity errors from PTV processing com pared to PIV processing for the present test case (com pare Figure 7 and Figure 8). Another indication that the low er errors are d riven by the better input d ata is that the ILAG approach perform s m uch better w hen using PTV d ata than w hen using PIV d ata because the acceleration is d irectly estim ated from the particle m otion.

Overview of results

To further quantify and com pare the perform ance of d ifferent m ethod s, the global errors (ε) is obtained by taking the r.m.s. error for all gridpoints (𝑁 𝑝 ) at all tim e instances (𝑁 𝑡 ): Table 3 and Table 4 show s the resulting errors for the clean (left colum ns) and noisy (right colum ns) cases for the tim e-resolved and m ulti-pulse input d ata, respectively. The m ethod s have been ord ered in term s of increasing error. Errors for the noisy case (right colum ns) can be seen to be slightly higher than for the clean case (left colum ns). Where this is not the case, this is attributed to a lack of statistical convergence. The results for the clean and noisy case follow a sim ilar trend w hich suggests that none of the m ethod s is critically sensitive to the am ount of noise ad d ed in this test case. N ote that Table 4 also includ es the result of the VIC+ m ethod using tom o-PTV input d ata. The error w hen using this inp ut d ata rem ains below the PIV results, but is slightly higher than w hen STB is used for particle tracking. This show s that less ad vanced particle tracking m ethod s than STB can still yield acceptable pressure reconstruction quality. The com parison of d ifferent particle tracking techniques w as left outsid e the scope of this assessm ent.

  d ifferent PIV-based pressure evaluation techniques is influenced by the characteristics of the flow as w ell as the quality of the m easurem ent d ata. The relative strengths, w eaknesses and d epend encies of d ifferent app roaches have been ad d ressed in the past in a num ber of com parative perform ance assessm ents (see van Oud heusd en (2013) for a review ). The present stud y contributes to these efforts by m aking a controlled assessm ent, through the application of a w id e variety of state-of-the-art velocim etrybased pressure d eterm ination techniques (see Table 1) to the sam e test case and em ploying the sam e und erlying synthetic particle-im age d ata. The techniques are classified in PIV-based techniques (1 to 4), w hich rely on the velocity d ata obtained w ith cross-correlation of particle im ages/ volum es betw een tw o tim e instances, and PTV-based techniques (4 to 8), that require tracking of ind ivid ual particles. A further d istinction is m ad e betw een techniques that w ork w ith d ata at a single tim e-instance only (3 and 4) and techniques that require tim e-resolved inp ut d ata (all others).

  ) and flow cond itions are d esigned to represent typical experim ental arrangem ents in the 27 × 27 cm 2 transonic-supersonic w ind tunnel at the Aerod ynam ics Laboratory at Delft University of Technology. The m od el m ain bod y has a d iam eter (D) of 50 mm and has an afterbod y of d iam eter (d ) of 20 mm (0.4D), resulting in a step height of 15 mm (0.3D). Wind tunnel w alls are om itted in the sim ulation for sim plicity and robustness. The free stream flow has a Mach num ber (M ∞ ) of 0.7, total pressure (p t ) of 200 kPa and total tem perature (T t ) of 285 K. Using the isentropic flow relations, this correspond s to a free stream pressure of 144,186 Pa, a free stream tem perature of 260 K, and a free stream velocity (U ∞ ) of 226 m .s -1 . The Reynold s num ber based on the m od el m ain bod y d iam eter (Re D ) is 1.3 m illion. The m easurem ent d om ain of the sim ulated PIV experim ent has a size of 1.2D × 0.47D × 0.08D (60 mm × 24 mm × 4 mm, L × H × W).

Fig. 1

 1 Fig. 1 Sketch of the geom etry; Diam eter D equals 50 mm. Filled colou red contours d epict the m ean stream w ise velocity in the virtual m easurem ent volum e for the sim ulated PIV experim ent.

  resolved d ata. This d istinction is m otived by consid erations w ith respect to the lim itations encountered und er actual m easurem ent cond itions. The m ulti-pu lse d ata can actually be m ad e available for high-speed flow s by u sing PIV system s w ith m ultiple ind epend ent laser and cam era system s (e.g.[START_REF] Souverein | Application of a d ual-plane particle im age velocim etry (d ual-PIV) technique for the unsteadiness characterization of a 18th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 4 -7, 2016 shock w ave turbulent bound ary layer interaction[END_REF] Lynch and Scarano 2014a). Continuous series of tim eresolved d ata on the other hand cannot realistically be obtained for the present flow case (i.e. a high-speed transonic flow ) w ithin the lim its of current m easurem ent capabilities. For low -speed flow s, m od ern high-speed laser and cam era equipm ent do however allow sufficiently high sam ple rates to perform tim e-resolved PIV m easurem ents. It is therefore d eem ed insightful to consid er the tim e-resolved d atasets as a m eans to assess pressure extraction capabilities for (lowspeed ) flow cond itions.

Figure 2

 2 Figure 2 show s the organisation of the extracted d atasets w ith the num bering of the velocity snapshots and im ages as w ell as the associated tim e separation. N ote that the sm aller tim e separation of 2 µs is equal to the tim e separation for PIV processing. The larger tim e separation of 10 µs used in the m ulti-pulse d ata w as fou nd to be su itable as tim e separation for d eterm ining the m aterial acceleration based on a prelim inary assessm ent by Blind e et al. (2014).

Fig. 2

 2 Fig. 2 Tim e-separation d iagram s for consecutive particle im ages, 4-pulse d ata (top) and continuous tim e-resolved d ata sets (bottom )

  -Velocity, d ensity and pressure d ata from the num erical sim ulation interpolated to the PIV grid positions Reference STB -Velocity, d ensity and pressure d ata from the num erical sim ulation interpolated to the particle positions 4. D etails of simulated experiment 4.1 Zonal detached eddy simulation (ZD ES) A d ed icated ZDES sim ulation (Deck 2005, 2012) is perform ed build ing on am ple experience w ith applying ZDES to axisym m etric base flow s (e.g. Deck and Thorigny 2007; Weiss et al. 2009; Weiss and Deck 2011). ZDES is closely related to the classical d etached eddy sim ulation. The m ain d ifference betw een the tw o approaches is that w ithin ZDES the user has to select the RAN S and DES d om ains. The selection of the appropriate m od e d epend s on the nature of the flow problem . Given that the flow separation is fixed by the geom etry, a m od e I (i.e. DES) has been retained for the separated area and a m od e 0 (i.e. URAN S) upstream from the separation occurring at the trailing ed ge of the forebod y.

Fig. 3

 3 Fig. 3 Detail of m esh w ith sketch of the geom etry; Filled coloured contours d epict the m ean stream w ise velocity in the virtual m easurem ent volum e for the sim ulated PIV experim ent.

Figure 4

 4 Figure 4 show s the m ean stream w ise velocity, m ean pressure as w ell as the turbulence intensity and the root m ean square (r.m .s) of the pressure fluctuations in a rad ial plane of the sim ulated d om ain. H ere, the pressure coefficient (Cp) and the turbulence intensity (T.I.) are d efined as:

  In these equations γ d enotes the ratio of specific heats, and 𝑢 𝑅𝑀𝑆 ′ r.m.s. of the fluctuating velocity com ponents in x-, y-and z-d irection, respectively. Results have been obtained on the basis of 5000 tim e sam ples and consid ering all planes in azim uthal d irection to im prove the statistical convergence. The bottom left figure show s the location of the m easurem ent volum e for the sim ulated PIV experim ent (see d ashed box). The m ean flow field (top left figure) show s a shear layer that em anates from the corner of the step at y/ D = 0.3 and grow s in d ow nstream d irection. The m ean reattachm ent of the shear layer occurs at approxim ately x/ D = 1.25, w hich is just outsid e the virtual PIV m easurem ent volum e. Below the shear layer, a recirculation region occurs, w here the m ean velocity reaches a m inim um of -0.30U ∞ , w hile the m inim u m m ean pressure (C p,m in ) is -0.23. Dow nstream of the point of m ean reattachm ent, a high pressure region occurs, w here the m axim um m ean pressure (C p,m ax ) of 0.18 is attained at x/ D = 1.50. The figures on the right show elevated levels of turbulence intensity and pressure flu ctuations in the shear layer and the reattachm ent region, w ith m axim a of T.I. m ax = 0.30 and Cp' RMS,m ax = 0.08 located at about x/ D = 1.0D; y/ D= 0.2D.

Figure 4 .

 4 Figure 4. Mean stream w ise velocity (top -left), m ean pressure (bottom -left), turbulence intensity (top right) and root m ean square of norm alised pressure flu ctuations (bottom -right); Top left figure: black solid line ind icates zero m ean stream w ise velocity; vectors have been subsam pled for clarity; Bottom -left figure: d ashed box ind icates the m easurem ent volum e of the sim ulated PIV experim ent.

Figure 5 d

 5 Figure 5 d epicts a representative realisation of the instantaneous stream w ise velocity, density, pressure and stream w ise m aterial acceleration. Especially the m aterial acceleration field (bottom right) show s sm all-scale flow structures originating from the corner of the step that break d ow n tow ard s the reattachm ent region w hich is characterised by the presence of interacting flow structures w ith a variety of length scales.

Figure 5 .

 5 Figure 5. Instantaneous velocity (top left), d ensity (top right), pressure (bottom left) and stream w ise m aterial acceleration (bottom right)

  Com paring the resulting position w ith the original position yield s a local error estim ate of below 10 -5 U ∞ , (< 0.0001 voxel d isplacem ent). In a second test, targeted at quantifying the errors d ue to interpolation, particle positions are integrated one sam ple tim e step forw ard in tim e using the original sim ulation d ata and the sim ulation d ata subsam pled at half the spatial and tem poral resolution. Com paring the resulting positions and assum ing that the error scales w ith 2 4 of the spatial and tem poral resolution (valid for 4th-ord er m ethod s), yield s typical local r.m.s. error estim ates of <0.003U ∞ (<0.03 voxel d isplacem ent).

  settings used for generating the synthetic im ages are chosen to represent typical PIV m easurem ents of a transonic flow. Particle im ages are record ed every 2 μs by four virtual cam eras w ith a chip size of 1624 × 800 pixels and a pixel pitch of 4.4 µm. The three-d im ensional particle positions are projected onto the two-d im ensional sensors via a pinhole cam era m od el.The four cam eras are placed in a cross configuration w ith yaw and pitch angles of +/ -30° to correspond to a total system aperture of 60°, an optim al configuration for tom ographic reconstruction[START_REF] Scarano | Tom ographic PIV: principles and practice[END_REF]. They are equipped w ith lenses w ith a focal length of 75 mm . The cam eras are placed at a virtu al d istance of 0.83 m resulting in a m agnification of 0.12 and a d igital resolution of 22.9 voxel.m m -1 . The resulting particle d isplacem ent in the free stream is about 11 voxel.Particle im ages are obtained using stand ard 2-D Gaussian integration. Tw o sets of im ages are generated : a set of id ealized particle im ages (Figure6, left) and another set of noisy particle im ages (Figure6, right). The id ealized im ages contain particles w ith a d iffraction spot size of 2 pixels, resulting in a source d ensity (N s ) of 0.16, and a peak intensity of 512 counts. The intensity of the background is zero. The noisy particle im ages contain particles w ith a larger d iffraction spot size of 2.5 pixels, resulting in a source d ensity (N s ) of 0.25, and a low er nom inal peak intensity of 342 counts. The seed ing is m od elled to consist of TiO 2 particles w ith a p rim ary crystal size of 55 nm w hich are know n to form agglom erates w ith a m ean size of about 400 nm (Schrijer et al. 2006). H ere, we rand om ly generate particle d iam eters from a Gaussian d istribution w ith m ean of 400 nm and a stand ard d eviation of 100 nm . The d istribution is truncated at tw o standard d eviations, to eliminate extrem ely sm all and large particles. The peak intensity is m od elled to vary w ith the 4th pow er of the physical p article d iam eter (Ad rian and Yao 1985), resulting in particle peak intensities range from 21 to 1731 counts. The particle d iam eter and peak intensity is attributed to a particle once and d oes not change through tim e or per cam era. Tw o noise sources are ad d ed : shot noise and therm al noise. Shot noise is d epend ent on the pixel intensity and is im plem ented so that for every pixel the intensity is a rand om realisation from a Poisson d istribution w ith m ean that is equal to the intensity of the noise-free image. Therm al noise is ind epend ent of the intensity and is generated as rand om realisations from a Poisson d istribution w ith a constant m ean across the im age. Mean therm al noise is subtracted before further processing.

Fig 6 .

 6 Fig 6. Extracts of a particle im age; clean (left) and noisy after subtraction of m ean noise intensity (right). In this exam ple, the extracts d epict the sam e region at the sam e tim e instance

Fig 7 .

 7 Fig 7. Centre-plane r.m.s error of velocity com ponents, expressed in voxels d isplacem ent obtained by PIV processing of the clean (left) and noisy (right) particle im ages; stream w ise (top), w all-norm al (mid d le) and out-of-plane (bottom ) d irection

Fig 8 :

 8 Fig 8: Centre-plane r.m.s. error of velocity com ponents, expressed in voxels d isplacem ent obtained by STB processing of the clean (right) and noisy (left) particle im ages; stream w ise (top), w all-norm al (mid d le) and out-of-plane (bottom ) d irection. N ote the d ifference in scale w ith respect to Figure 7.

  s rely on the m om entum equation (equation 1) to d eterm ine the pressure grad ient from the m aterial acceleration, w hich is com puted from the flow kinem atic d ata. To account for com pressibility effects, the m om entum equation is rew ritten to elim inate d ensity as ind epend ent variable, und er the assum ption of inviscid and ad iabatic flow (van Oud heusd en et al. 2007). The high Reynold s num ber of the current flow justifies the inviscid flow approxim ation, w hile the absence of heat transfer from the m od el supports the ad iabatic assu m ption.

Fig 9 .

 9 Fig 9. R.m .s. Errors from assu m ing ad iabatic tem perature (top -left) and isentropic pressure (topright) as w ell as r.ms. im pact on pressure solution (bottom figures)

18th

  International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 4 -7, 2016

Fig 10 .

 10 figure. Each other row correspond s to a d ifferent technique. Left figures show sam ple instantaneous (reconstructed ) pressure field s in the centre-plane, w ith the correspond ing instantaneous errors in the centre figures. Right figures show the r.m .s. errors from all available tim e steps and planes in z-d irection. As the error of the reconstructed pressure field s generally increases tow ard s the outer planes, this app roach slightly increases the r.m.s. error level w ith

Fig 11 .

 11 Fig 11. Results and errors from applying PIV-based pressure reconstruction techniques to m ultipulse inp ut d ata. (Reconstructed ) instantaneous pressure field s in the centre-plane (left colum n), instantaneous error in the centre-plane (centre colum n) and r.m.s. errors for all 11 planes (right colum n).

18th

  International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 4 -7, 2016

  figures) has been ad ap ted w ith respect to section 6.1 to better visualize the error d istribution. The scaling of r.m.s. error field s (right figures) rem ains u nchanged to facilitate d irect com p arison w ith the PIV-based techniques.

Fig 12 .

 12 Fig 12. Results and errors from applying STB-based pressure reconstruction techniques to tim eresolved input d ata. (Reconstructed ) instantaneous pressure field s in the centre-plane (left colum n), instantaneous error in the centre-plane (centre colum n) and r.m.s. errors for all 11 z-planes (right colum n); r.m.s. error for Voronoi-based integration consid ers 5 z-planes only.

Fig 13 .

 13 Fig 13. Results and errors from applying STB-based pressure reconstruction techniques to m ultipulse inp ut d ata. (Reconstructed ) instantaneou s pressure field s in the centre-plane (left colum n), instantaneous error in the centre-plane (centre colum n) and r.m.s. errors for all 11 z-planes (right colum n); r.m.s. error for Voronoi-based integration consid ers 5 z-planes only.

  for d ifferent PIV-based and PTV-based pressure evaluation techniques has been d eveloped by constructing a sim ulated experim ent from ZDES sim ulation d ata. Important experim ental error sources are replicated by sim ulating the entire measurem ent chain. To enable an assessm ent of the im pact of m easurem ent noise, the test case includ es both cases w ith and w ithout ad d ed artificial noise. The test case includ es m u lti-pulse d ata (4 particle im ages; 2 velocity snapshots), w hich is representative of the acceleration m easurem ent proced ure available in high-speed flow s, as w ell as longer series of continuous tim e-resolved d ata w hich are typically achievable only for low -speed flow s. Both tom ographic PIV algorithm s and the STB algorithm have been applied to enable the testing of both PIV-and PTV-based pressure d eterm ination proced ures. A variety of state-of-the-art pressure d eterm ination techniques has been applied to a single test case consisting of a transonic flow over an axisym m etric step. The w orking principles as w ell as the im plem entation of these techniques have been d escribed . Using the test case, the techniques have been assessed in term s of their capability to reconstruct the reference pressure field as w ell as the im pact of noise in the input d ata and the benefit of using tim e-resolved d ata. This led to the follow ing m ain conclusions: 1. A range of suitable m ethod s exists that can reconstruct instantaneous pressure field s from a variety of PIV/ PTV input d ata. The m ain features of the p ressure field s can be

Table . 1

 . Overview of pressure evaluation techniques consid ered in this stud y.

	N o. Abbrev. Velocity	D etermination of material	Integration of pressure
			measurement	acceleration	gradient
	1	EUL	PIV	Eulerian approach	Poisson solver
	4	ILAG	PIV	Iterative least-square pseud o-	Poisson solver
				tracking
	2	TH	PIV	Taylor 's hypothesis approach	Poisson solver
	3	IVIC	PIV	Instantaneous Vortex-in-Cell	Poisson solver
	5	FF1	PTV ('Shake-The-Box') + Flow Fit	FFT
	6	FF2		PTV ('Shake-The-Box')	Flow Fit2
	7	VIC+	PTV ('Shake-The-Box') + Vortex-in-Cell-plus	Poisson solver
	8	VOR		PTV ('Shake-The-Box')	Voronoi-based
					integration

Table . 2

 . Overview

	of datasets. 'Clean' and 'N oisy' refer to the cases w ithout or w ith artificial
	noise ad d ed to the particle im ages. The 'clean' and 'noisy' cases consid er d ifferent tim e instances,
	so that the exact noise levels cannot be obtained by com paring the tw o cases.
	Type	D escription
	Continuous tim e-	-Particles im ages at 42 consecutive tim e instances
	resolved d ata	-41 velocity snapshots obtained through PIV processing.
	('clean' and 'noisy')	-Velocity and acceleration d ata at particle positions at sam e 21 tim e
		instance as PIV velocity field s, obtained through STB processing.
	4-pulse d ata	-21 sets of particle im ages at 4 consecutive tim e instances
	('clean' and 'noisy')	-21 sets of 2 consecutive velocity snapshots obtained through PIV
		processing
		-Velocity and acceleration d ata at particle positions at the m id d le tim e
		instances of the sets of particle im ages/ velocities, obtained through
		STB processing.

Table . 2

 . Overview of datasets (continued )

	Type	D escription
	Velocity statistics	-Mean and r.m.s. velocity field s based on full d ataset (5000 snapshots)
		obtained through PIV processing of 'clean' p article im ages
	Cam era calibration	-Calibration im ages for four virtual cam eras
	Reference PIV

Table 3 .

 3 Global errors for tim e-resolved input d ata. Values are based on 171 × 67 × 11 d ata points w ith the exception of VOR w hich is based on 165 × 61 × 6 d ata points

Pressure reconstruction Source of input data Global error [× 10 -2 ]

			Clean	N oisy
	VIC+	STB	0.68	0.63
	FF1	STB	0.72	0.78
	FF2	STB	0.72	0.78
	ILAG	STB	0.90	0.91
	VOR	STB	1.70	1.66
	ILAG	PIV (9 snap shots)	1.43	1.69
	EUL	PIV (2 snap shots)	1.52	2.09
	TH	PIV (1 snap shot)	2.23	2.54
	IVIC	PIV (1 snap shot)	3.44	2.99

Table 4 .

 4 Global errors for m ulti-pulse input d ata. Values are based on 171 × 67 × 11 d ata points w ith the exception of VOR w hich is based on 165 × 61 × 6 d ata points.

	Pressure	Source of input data	Global error [× 10 -2 ]
	reconstruction		Clean	N oisy
	VIC+	STB	0.70	1.02
	FF1	STB	0.81	1.07
	FF2	STB	0.86	1.10
	ILAG	STB	0.99	1.23
	VIC+	Tom o-PTV	0.95	1.36
	ILAG	PIV	1.67	2.01
	EUL	PIV	1.80	2.11
	VOR	STB	3.57	2.35
	TH	PIV (1 snap shot)	2.66	2.78
	IVIC	PIV (1 snap shot)	3.34	3.39

18th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 4 -7, 2016

system is subsequently solved for pressure. The m ethod is im plem ented using pairs of snapshots w ith a tim e separation of 10 μs for both the multi-pulse and the tim e resolved case.

Iterative least-square pseudo-tracing (ILA G)

Im aginary particles are traced through subsequent velocity field s (see e.g. [START_REF] Liu | Instantaneous pressure and m aterial acceleration m easurem ents u sing a fourexposure PIV system[END_REF]. In the present im plem entation, im aginary particles located at the grid points are integrated forw ard and backw ard in tim e in com bination w ith linear interpolation. The m aterial acceleration is obtained from a first-ord er least-square fit through the velocity values at the positions of the im aginary particles, d isregard ing any particle positions outsid e the d om ain. For tim e-resolved PIV, it could also be estim ated w ith the FTEE algorithm [START_REF] Jeon | Fluid trajectory evaluation based on an ensem ble-averaged cross-correlation in tim e-resolved PIV[END_REF]. The m aterial acceleration is used to obtain a better estim ate of the im aginary particle track in a next iteration.

The tem poral length of each trajectory is 9 tim e steps for the case of tim e-resolved input d ata and 2 tim e steps for the case of m ulti-pulse input d ata. In ord er to com pensate the exclud ed outw ard fluid m otion at the volum e bound ary and to prevent correspond ing d ecrem ent of the d ynam ic range, the tem poral trajectory length is extend ed tow ard the opposite sid e. The spatial integration of the pressure grad ient is cond ucted by solving the Poisson equation form ulated as a least-square error cost function on a staggered grid . N eum ann bound ary cond itions are im plicitly integrated . A Dirichlet bound ary cond ition w ith a value of 0 is prescribed in a single grid point. In a next step the solution is offset using the isentropic pressure ratio (see equation 5) sam pled over a region selected based on the Frobenius norm of the m aterial acceleration tensor.

A d etailed d escription of the present im plem entation can be found in [START_REF] Jeon | Least-square reconstruction of instantaneous pressure field around a bod y based on a d irectly acquired m aterial acceleration in tim eresolved PIV[END_REF]. The m ethod is im plem ented for both PIV and STB input d ata. In case of the latter, the tem poral length is 2 tim e steps for both the cases of tim e-resolved and m ulti-pulse input.

Taylor's hypothesis approach (TH)

Taylor 's hypothesis consists of the assum p tion that ad vection d ue to turbulent circulation is sm all and that 'frozen' turbulence therefore ad vects w ith the flow. Ad opting Taylor 's hypothesis enables the estim ation of the m aterial acceleration from a single velocity snapshot. The choice of convection velocity is crucial for the accuracy of this m ethod and is taken here as the m ean stream w ise velocity, averaged in z-d irection. A d etailed d escription of the im plem entation of this approach for pressure d eterm ination can be found in de Kat and Ganapathisubram ani (2013) and [START_REF] Laskari | Full-field pressure from snapshot and tim eresolved volum etric PIV[END_REF]. The subsequent integration of the pressure grad ient is perform ed sim ilarly as for the Eulerian approach d escribed above.

reconstructed from a single PIV velocity snapshot w hile highly accurate pressure reconstructions can be obtained by using PTV approaches in com bination w ith m ore ad vanced techniques.

2. For the present input d ata, the PTV-based techniques prod uce m ore accurate pressure reconstructions than the PIV-based approaches. This d ifference in perform ance is attributed to a com bination of higher spatial resolution of the input d ata and better use of tim e inform ation in the d ata sets. It should be realized how ever that the com parison is biased in the sense that the PIV processing used here is typical for experim entally obtained particle im ages and no m ulti-fram e cross-correlation techniques have been applied w hereas the application of the STB algorithm has not yet been experim entally d em onstrated for high-speed flow cond itions.

3. The use of longer series of tim e-resolved input d ata allow s m ore accurate reconstructions of the reference pressure field . N evertheless, in light of practical d ifficulties in obtaining m ulti-pulse d ata for high-speed flow s, it is encouraging that it appears to be possible to reconstruct the m ain features of the pressure field alread y from a single velocity field . 4. N oise in the input d ata typically red uces the accuracy of the reconstructed pressure field s, but none of the m ethod s proved to be critically sensitive to the amount of noise that w as ad d ed in the present test Given these conclusions on the availability and perform ance of pressure reconstruction algorithm s and their d epend ence on suitable experim ental d ata, future ad vances are best attained through d em onstration and further developm ent of experim ental capabilities.

Acknow ledgements

This w ork is sup ported by the European FP-7 project "N IOPLEX", grant agreem ent 605151.