Naweiluo Zhou
email: naweiluo.zhou@inria.fr

Gwenaël Delaval
email: gwenael.delaval@inria.fr

Bogdan Robu
email: bogdan.robu@gipsa-lab.grenoble-inp.fr

Éric Rutten
email: eric.rutten@inria.fr

Jean-Franc ¸ois Méhaut
email: jean-francois.mehaut@imag.fr

Control of Autonomic Parallelism Adaptation on Software Transactional Memory

Keywords: autonomic, transactional memory, feedback control, synchronization, parallelism adaptation

Parallel programs need to manage the trade-off between the time spent in synchronization and computation. A high parallelism may decrease computing time while increase synchronization cost among threads. A way to improve program performance is to adjust parallelism to balance conflicts among threads. However, there is no universal rule to decide the best parallelism for a program from an offline view. Furthermore, an offline tuning is error-prone. Hence, it becomes necessary to adopt a dynamic tuning-configuration strategy to better manage a STM system. Software Transactional Memory (STM) has emerged as a promising technique, which bypasses locks, to address synchronization issues through transactions. Autonomic computing offers designers a framework of methods and techniques to build automated systems with well-mastered behaviours. Its key idea is to implement feedback control loops to design safe, efficient and predictable controllers, which enable monitoring and adjusting controlled systems dynamically while keeping overhead low. We propose to design feedback control loops to automate the choice of parallelism level at runtime to diminish program execution time.

I. INTRODUCTION

Multicore processors are ubiquitous, which enhance program performance through thread parallelism (number of simultaneous active threads). Although a high parallelism degree shortens execution time, it may also potentially increase synchronization time. Therefore, it is crucial to find the trade-off between synchronization and computation. The conventional way to address synchronization issues is via locks. However, locks are notorious for various issues such as the likelihood of deadlock and the vulnerability to failure and faults. Additionally it is not straightforward to analyse interactions among concurrent operations.

Transactional memory (TM) emerges as an alternative parallel programming technique, which addresses synchronization issues through transactions. Accesses to shared data are enclosed in transactions which are executed speculatively without being blocked by locks. Various TM schemes have been developed [START_REF] Herlihy | Transactional memory: architectural support for lock-free data structures[END_REF], [START_REF] Felber | Dynamic performance tuning of word-based software transactional memory[END_REF], [START_REF] Damron | Hybrid transactional memory[END_REF] including Hardware Transactional Memory (HTM), Software Transactional Memory (STM) and Hybrid Transactional Memory (HyTM). In this paper, we present the work on runtime program parallelism adaptation under STM systems where the synchronization time originates in transaction aborts. There are different ways to reduce aborts, such as conflict resolutions, ways to detect conflicts, designs of version management and the level of thread parallelism.

Online parallelism adaptation began to receive attention recently in TM systems. The level of a suitable parallelism in a program can significantly affect system performance. However it is onerous to set a suitable parallelism degree for a program offline especially for the one with online behaviour variation. When apropos of the program with online behaviour fluctuation, there is no unique thread number can enable its optimum performance. Therefore, the natural solution is to monitor a program at runtime and alter its parallelism when necessary.

We introduce autonomic computing [START_REF] Kephart | The vision of autonomic computing[END_REF] to STM systems to automatically regulate online program parallelism. In this paper we argue that online adaptation is necessary and feasible for parallelism management in STM systems. We demonstrate that the program performance is sensitive to the parallelism. We present two effective profiling frameworks for the parallelism adaptation on TinySTM [START_REF] Felber | Dynamic performance tuning of word-based software transactional memory[END_REF]. The main contributions of our paper are as follows:

1) We present two adaptive profiling frameworks that detect the suitable parallelism degree in order to optimize system performance at runtime. 2) We dynamically resolve a phase detection method.

The rest of the paper is organized as follows. Section II summaries the background and related work. Section III details the profiling procedures and the online parallelism adaptation methods. Section IV presents the implementation details. Section V shows the results. Section VI discusses the pros and cons of our adaptive framework and Section VII concludes the paper and gives future work.

II. BACKGROUND AND RELATED WORK

A. Background on Software Transactional Memory

Transactional memory (TM) is an alternative synchronization technique. Its accesses to shared data are enclosed in transactions which are executed speculatively without being blocked by locks. Each transaction makes a sequence of tentative changes to shared memory. When a transaction completes, it can either commit making the changes permanent to memory or abort hence discarding the previous changes [START_REF] Herlihy | Transactional memory: architectural support for lock-free data structures[END_REF]. Two parameters are often used in TM to indicate system performance, namely commit ratio and throughput. Commit ratio (CR) equals the number of commits divided by the sum of number of commits and number of aborts; it measures the level of conflict or contention among the current transactions. Throughput is the number of commits in one unit of time; it directly indicates program performance. TM can be implemented in software, hardware or hybrid. Different mechanisms explore the trade-off that impact on performance, programmability and flexibility. In this paper, we focus on STM systems and utilise TinySTM [START_REF] Felber | Dynamic performance tuning of word-based software transactional memory[END_REF] as our experimental platform. TinySTM is a lightweight STM system which adopts a shared array of locks to control the concurrent accesses to memory and applies a shared counter as its clock to manage transaction conflicts.

The performance of STM systems has been continuously improved. Studies to improve STM systems mainly focus on the design of conflict detection, version management and conflict resolution. Conflict detection decides when to check read/write conflicts. Version management determines whether logging old data and writing new data to memory or vice versa. Conflict resolution, which is also known as contention management policy, handles the actions to be taken when a read/write conflict happens. The goal of the above designs is to reduce wasted work. The amount of wasted work resides in the number of aborts and the size of aborts. The higher contention in a program, the larger amount of wasted work. The time spent in wasted work is the synchronization time in a STM view. Apart from diminishing wasted work, one way to improve STM system performance is to trim computation time. A high parallelism may accelerate computation but resulting in high contention thus high synchronization time. Hence parallelism can significantly affect system performance.

B. Background on Autonomic Computing

Autonomic computing [START_REF] Huebscher | A survey of autonomic computing -degrees, models, and applications[END_REF] is a concept that brings together many fields of computing with the purpose of creating computing systems that self-manage. A system is regarded as an autonomic system if it supports one of the following features [START_REF] Kephart | The vision of autonomic computing[END_REF]: (1) self-optimization, the system seeks to improve its performance and efficiency on its own; (2) self-configuration, when a new component is introduced into a system, the component is able to learn the system configuration. (3) selfhealing, the system is able to recover from failures; (4) selfprotection, the system can defend against attacks.

In this paper, we concentrate on the first feature: selfoptimization. We introduce feedback control loops to achieve autonomic parallelism management. A classic feedback control loop is illustrated in Fig. 1 in the shape of a MAPE-K loop [START_REF] Kephart | The vision of autonomic computing[END_REF].

In general, a feedback control loop is composed of (1) an autonomic manager, (2) sensors (collect information), (3) effectors (carry out changes), (4) managed elements (any software or hardware resource). An autonomic manager is composed of five elements: a monitor (used for sampling), an analyser (analyse data obtained from the monitor), knowledge (knowledge of the system), plan (utilise the knowledge of the system to carry out computation) and execute (perform changes). It is worth noting that an autonomic manager can only incorporate a part of the five elements.

Managed Elements sensors Effectors

Fig. 1. A MAPE-K control loop. It incorporates an autonomic manager, a sensor, an effector as well as a managed element among which the autonomic manager plays the main role.

C. Related Work

It has been been addressed in previous work [START_REF] Ansari | Advanced concurrency control for transactional memory using transaction commit rate[END_REF], [START_REF] Ravichandran | F2C2-STM: Flux-based feedbackdriven concurrency control for STMs[END_REF], [START_REF] Didona | Identifying the optimal level of parallelism in transactional memory applications[END_REF], [START_REF] Rughetti | Machine learningbased self-adjusting concurrency in software transactional memory systems[END_REF] to dynamically adapt parallelism degrees via control techniques to reduce wasted work in TM systems. Ansari et. al. [START_REF] Ansari | Advanced concurrency control for transactional memory using transaction commit rate[END_REF] proposed to adapt parallelism online by detecting the changes of the application's CR. The regulation action is made to the parallelism if the CR falls out of the pre-set CR range or is not equal to a single preset CR value. This is based on the fact that CR falls during highly contended phases rises with low contended phases. Ansari et. al. gave five different algorithms which decide the profile length and the level of parallelism. Ravichandran et. al. [START_REF] Ravichandran | F2C2-STM: Flux-based feedbackdriven concurrency control for STMs[END_REF] presented a model which adapts the thread number in two phases: exponential and linear with a feedback control loop. Rughetti et. al. [START_REF] Rughetti | Machine learningbased self-adjusting concurrency in software transactional memory systems[END_REF] utilise a neural network to estimate performance of STM applications. The neural network is trained to predict execution time of wasted transactions which in turn is utilised by a control algorithm to regulate the parallelism. Didona et. al. [START_REF] Didona | Identifying the optimal level of parallelism in transactional memory applications[END_REF] introduces an approach to dynamically predict the parallelism based on the workload (duration and relative frequency, of read-only and update transactions, abort rate, average number of writes per transaction) and throughput, through one feedback control loop its prediction can be continuously corrected.

III. AUTONOMIC PARALLELISM ADAPTATION

In this section, we detail the design of two feedback control loops. We present two approaches to dynamically determine optimum parallelism. We firstly introduce a simple model which searches optimum parallelism, then we present a more sophisticated model (probabilistic model) based on probability theory which predicts the optimum parallelism. We measure three parameters from the STM system, namely the number of commits, the number of aborts and physical time. The number of commits and the number of aborts are addressed as commits and aborts subsequently. We choose CR and throughput to denote program performance, as CR and throughput are both sensitive to parallelism variation. Either is by itself not sufficient enough to represent program performance as:

• A high throughput shows fast program execution whereas a low throughput indicates slow program progress. Nevertheless a low throughput may be caused by low parallelism or simply just a low number of transactions taking places.

• CR indicates the conflicts among threads. A high CR means low synchronization time whereas a low CR mean high synchronization cost. But a low CR can bring a high throughput when a large number of transactions are executing concurrently, whereas a high CR may give low throughput due to a small number of transactions executing concurrently.

The controller observes the CR to detect contention fluctuation and enable corresponding control actions. The correctness of the control actions are verified by checking if the throughput is improved after the actions for parallelism adaptation.

A. Overview of Profiling Algorithm on simple model

To describe the simple model, we firstly give an overview of the profiling procedure and later describe the model through the prism of control theory.

To achieve autonomic parallelism adaptation which provides a program with its optimum parallelism, we propose to periodically profile the parallelism and select the value that achieves the highest throughput. By observing CR, we can obtain the contention information of an application. CR usually fluctuates in a certain range within the same phase. When a program enters a new phase, the current parallelism produces a different CR which falls out of the current CR range. The CR fluctuation triggers a new parallelism adaptation action. Initially the two CR thresholds (upper CR and lower CR thresholds) are both set to be 0 and are trained in the later profile stage. The detail of the profiling procedure is illustrated in Fig. 2.

The parallelism profiling procedure starts once the program starts. Initially the program creates a pool of threads, among which only 2 threads are awaken and the rest are suspended. At each decision point, which corresponds to one state of the automaton in Fig. 3(b), the control loop (see section III-B) is activated to adapt the parallelism or suspend the parallelism adaptation. A profile length is a fixed period (marked by logic time: the number of commits) for information gathering, such as commits, aborts and time. A parallelism profile interval is composed of a continuous sequence of profile lengths within which the parallelism is adjusted and the CR range is computed. The non-action interval consists of one or a continuous sequence of profile lengths, within which the parallelism adaptation is suspended. The duration of parallelism profile interval and non-action interval are not fixed values as shown in Fig. 2. The above procedure continues until the program terminates. Fig. 3(a) gives the structure of the complete platform that forms a MAPE-K feedback control loop. The autonomic element is composed of the STM system, benchmarks, inputs, outputs and the autonomic manager. The autonomic manager, which can be also regarded as the controller, is described as an automaton as shown in Fig. 3(b). The automaton consists of four states, and the program can only reside on one state at each decision point.

B. Feedback Control Loop of the Simple Model

1) Control Objective: Under control theory terminology, the control objective of the feedback control loop is to maximize the throughput and diminish the global execution time.

2) Inputs and Outputs:

As shown in Fig. 3(a), the inputs are commits, aborts, physical time and current active thread number. The outputs are the optimum parallelism, the parallelism profile flag. 3) Decision Functions: Three decision functions cooperate to make decisions: a parallelism decision function (adjusts parallelism), a profile decision function (enable the parallelism profile actions) and a CR range decision function (dictates the phase variation). We describe design of the automaton as it elucidates the relation among the decision functions, as well as how the parallelism decision function and profile decision function are designed. The CR range decision function is presently lastly.

The automaton commences at the state increase tn to increase the thread number, since the parallelism is set to be the minimum at the starting point. In each parallelism profile interval, the parallelism can either continuously increase or decrease. The direction of parallelism regulation (increase or decrease) is determined by the profile decision function. If the current throughput is greater than the previous throughput, one thread is awaken or suspended and the current throughput is recorded as the maximum value. The state transfers to stop profile when the current throughput is less than the maximum throughput. At the final decision point of a parallelism profile interval (the state stop profile), the parallelism is set to be the value which yields the maximum throughput. A new CR range may be determined at the end of a parallelism profile interval, as we will detail later in this section. The automaton then enters from the state stop profile to the state no tn control which corresponds to the non-action interval in Fig. 2. At each decision point of a non-action interval, the profile decision function checks if it should start parallelism adaptation. More specifically, if the CR falls into the CR range, the program stays in the no tn control state. Otherwise a boolean value is set indicating the direction of the parallelism The automaton jumps into increase tn if CR is higher than upper CR threshold or decrease tn if CR is lower than the lower CR threshold. It is worth noting that, in case the value of the upper threshold is 100%, and the program CR is also 100% (it happens when only read operations or no conflicts across transactions), higher parallelism is assigned to the program.

It is worth noting that the throughput often fluctuates before reaching the optimum value (shown in Fig. 4). To prevent a parallelism profiling procedure from terminating at a local maximum throughput, parallelism profiling procedure continues until the throughput decreases over 10% of the maximum value (10% is an empirical value which can be tuned).

The parallelism decision function resides in the state increase tn and decrease tn. Parallelism adaptation is activated when a program enters a new phase. A new phase is denoted by when CR fluctuates out of a certain range. It is onerous to determine such a CR range offline, especially for some programs with online performance variation. Additionally, a constant CR range impedes programs to search its optimum parallelism. Therefore it becomes necessary to dynamically resolve a CR range. We add a CR range decision function. Therefore at the end of a parallelism interval, a new CR range is prescribed. The function is activated at the state stop profile. The two thresholds of the CR range are the CR values produced by running with one more or one less parallelism degree than the optimum one.

C. Probabilistic Model

The approach to manipulate one thread number at each decision point engages long profiling time. This section presents a probabilistic model which predicts favourable parallelism after one profile length based on the CR and current active thread number. The profiling procedure is similar to the simple model as shown in Fig. 2.

The probabilistic model shares the same inputs, outputs and control objective as the simple model. It also incorporates three decision functions. As the CR range decision function is equivalent to that in Section III-B3, we only describe the parallelism decision function and the profile decision function in this section.

1) Decision Functions: We firstly present the parallelism decision function. It is based on two assumptions:

• the same amount of transactions are executed in each active thread during a fixed period, as every thread shows similar behaviour in our TM benchmarks

• the probability of one commit (see Section II-A for definition) approaches a constant, as there is a large amount of transactions executed during the fixed period making the probability of conflicts between two transactions approaches a constant.

The detail of the derivation of the decision function is presented in [START_REF] Zhou | Autonomic Parallelism Adaptation on Software Transactional Memory[END_REF], we only provide the final equation here. The optimum thread number is calculated by (1):

nopt = - n -1 ln CR (1)
Where n opt represents the optimum parallelism, n represents the number of running threads and CR is the commit ratio.

In this paragraph, we describe the automaton to elucidate the relation of the three decision functions as well as the design of the profile decision function. As illustrated in Fig. 5, the automaton commences (with the maximum parallelism) from the predict tn state which yields an estimation of optimum parallelism degree. The predicted parallelism is applied for one subsequent profile length and the automaton unconditionally enters the verify state to verify the correctness of the predicted parallelism. The new parallelism is only applied subsequently when the current throughput is larger than the throughput recorded before the new parallelism is applied. This leads to the state transfer to the CR range state where a new CR range is prescribed. The stop profile state only disables the parallelism profile action when the parallelism does not alter. Otherwise it recovers the previous value (CR range remains unchanged in this case). Contrary to the simple model, the probabilistic model requires an individual state to obtain the CR range.

IV. IMPLEMENTATION

There are two methods of collecting profile information in a parallel program. A master thread can be utilised to the interesting information of itself. An alternative way is to collect the information via all threads. The first method requires little synchronization cost to gather information but the obtained information may not represent the global view. Also the master thread must be active during the whole program execution possibly making it terminate earlier than the other threads, meaning that the fair execution time among threads can not be guaranteed. The second method may suffer from synchronization cost but the profile information gathered represent the global view. More importantly, a strategy to enforce fair time slides can be employed among threads. We choose the second method. The synchronization cost of information gathering is negligible for most of our applications.

We implement a monitor to collect the profile information, adjust the parallelism and the race condition. The monitor is a cross-thread lock which consists of concurrent-access variables by threads. The main variables of the monitor are commits, aborts, two FIFO queues recording the suspended and active threads, the current active thread number, the optimum thread number and the throughput. There are three entry points of the monitor. The first entry point is upon threads initialization, where some initial values (e.g., thread id) are set for the threads and some threads are suspended. The second entry point is upon transaction committing, where commits are accumulated and where the control functions take actions. The third entry point is upon a thread exiting, where one suspended thread is awaken when one thread exits.

We use logic time (number of commits) to mark the profiling length rather than the physical time. As the size of transaction varies in various applications leading to the significant variation of execution time. The choice of the profile length mainly depends on the total amount of transactions of an application. The applications with the same magnitude of transactions share the same profile length. For instance, genome and vacation (two benchmarks from STAMP) share the same profile length as the total number of the transactions in the two applications are on the same magnitude (10 6 ). Time overhead is added to each transaction when calling and releasing a monitor. The overhead is insignificant for the transactions with medium or long length, however it gives nontrivial influence to the short-length transactions. This overhead can be reduced through diminishing the frequency of calling the monitor, i.e. the monitor is called every 100 commits rather than every commit.

V. PERFORMANCE EVALUATION

In this section, we present the results from 6 different STAMP [START_REF] Cao Minh | STAMP: Stanford transactional applications for multi-processing[END_REF] benchmarks and two applications from Eigen-Bench [START_REF] Hong | EigenBench: A simple exploration tool for orthogonal TM characteristics[END_REF]. EigenBench and STAMP are widely used for performance evaluation on TM systems. The data sets cover a wide range from short-length to long-length transaction, short to long program execution time, from low to high program contention. Table I presents the qualitative summary of each application's runtime transactional characteristics: TX length (the number of instructions per transaction), execution time, and contention (the global contention). The classification is based on the application with its static optimal parallelism. A transaction with execution time between 10 us and 1000 us is classified as medium-length. The contention between 30% and 60% is classified as medium. The execution time between 10 seconds and 30 seconds is classified as medium.

A. Platform

We evaluate the performance on a SMP machine with 4 processors of 6 cores each. Every pair of cores share a L2 cache (3072KB) and every 6 cores share a L3 cache (16MB). This machine holds 2.66GHz frequency and 63GB RAM. We utilise TinySTM as our STM platform.

B. Benchmark Settings

We show two different data sets of EigenBench. One with stable behaviour and one with diverse phases. As EigenBench does not shows phase variation, we modify its source code to enable diverse phases. EigenBench include 3 different arrays which provide the shared transactional accesses (Array1), private transactional accesses (Array2) and non-transactional accesses (Array3). We vary the size of Array1 at runtime making the conflict rate vary at runtime. More specifically, within the first 40% amount of the transactions, the size of Array1 keeps the value given by the input file. From 40% to 70%, the array size is shrunk to be 16% of the original value and afterwards the size is set to be 33% of the given value. The main inputs of both data sets are given in Fig. 6.

We have evaluated [START_REF] Ansari | Advanced concurrency control for transactional memory using transaction commit rate[END_REF] Fig. 6. Two data sets of EigenBench inputs for 24 threads labyrinth. Two applications namely bayes and kmeans from STAMP are not taken into account in the paper. As bayes exhibits non-determinism [START_REF] Ruan | STAMP need not be considered harmful[END_REF]: the ordering of commits among threads at the beginning of an execution can dramatically affect the execution time. The aforementioned applications have represented a wide range of characteristics of TM applications, therefore we do not present the results from kmeans due to the page limit. The inputs of the six selected applications are detailed in Fig. 7.

intruder -a8 -l176 -n109187 ssca2 -s20 -i1.0 -u1.0 -l3 -p3 genome -s32 -g32768 -n8388608 vacation -n4 -q60 -u90 -r1048576 -t4194304 yada -a15 -i inputs/ttimeu1000000.2 labyrinth -i random-x1024-y1024-z7-n512.txt

C. Results

We firstly present the results of the execution time comparison of the two autonomic parallelism adaptation approaches with static parallelism. The maximum parallelism is 24 which is the number of the available cores. The minimum parallelism is restricted to be 2, as we are only concerned with parallel applications. All the applications are executed 10 times and results are the average execution time. Fig. 10 elucidates the runtime parallelism variation with simple and probabilistic model of genome. genome experiences three phases at runtime. The first phase is short (two or three profile lengths) which contains both read and write operations. During the second phase, the transactions only thread number execution time (second) include read operations resulting in 100% of CR, hence the maximum parallelism is applied. The third phase generates high contention, hence low parallelism is given. As shown in Fig. 10, the simple model spends some time before reaching the optimum parallelism, thus some staircases reflect in the figure . The probabilistic model reacts fast to respond the CR and phase change, thus some abrupt parallelism changes are shown. As less time spent to reach the optimum parallelism, the performance of probabilistic model outperforms that of the simple model.

• • • • • • • • • • • 0
• • • • • • • • • • • • 0 2 4
• • • • • • • • • • • 0
The autonomic parallelism adaptation aims to regulate the parallelism which retains throughput at the optimum level at each phase. Ideally the throughput from the adaptive models should rival the one with the static parallelism which achieves the maximum throughput. Due to the page limit, we only present the online throughput change of one application genome. Fig. 11 elucidates its online throughput change with static and adaptive parallelism.

VI. DISCUSSION

The overhead of our approaches mainly originate in three aspects. (1) Thread migration. This can introduce a large overhead especially when the parallelism is adapted at runtime.

(2) The choice of the thread number to manipulate at each decision point in the simple model. We simply choose to manipulate one thread number each point which delays the procedure to reach the optimum parallelism. (3) The choice of throughput variation rate in the simple model. We keep profiling even if the current throughput is slightly lower than the recorded maximum value in order to avoid the parallelism profiling to be terminated at a regional maximum point.

Performance penalty can occur when parallelism varies, yet is trivial on shared memory. In addition, to reduce the penalty, a thread is only suspended when its current transaction commits. The CR range decision function shows its limitation on phase detection for vacation. As its CR tend to fluctuate frequently over the CR range, yet the program remains in the same phase. Therefore both model perform insufficiently vacation.

The simple model only manipulates one thread number at each decision point thus long time spent in searching the optimum parallelism. However this eludes the possibility of skipping the optimum parallelism. The application starts with two threads activated rather than the maximum value to avoid excessive contention which may prevent the program from progressing. However this setting brings high parallelism profiling time to the applications which require a high parallelism, this is especially true for labyrinth that requires maximum parallelism to achieve its maximum throughput. Such an overhead is difficult to be compensated by the performance improvement brought by the optimum parallelism. Therefore, the simple model gives heavy overhead against the best case to the application with long-length transactions. genome requires the maximum parallelism at the second phase but minimal parallelism at the third phase leading to a high overhead originating in parallelism descending. Such overhead is especially large when the application generates high contention.

The probabilistic model predicts the parallelism in one step which diminishes the profiling time and often gives better performance than the simple model. However it has the potential risk of overreacting to the phase variation as intruder. Lastly, the probabilistic model relies on two assumptions which are based on the ideal situations, thus some errors are imposed inevitably in reality, meaning that the predicted parallelism may be slightly different as the optimum one. Thereby when phase changes are trivial, the simple model can outperform the probabilistic model (e.g., yada).

VII. CONCLUSION AND FUTURE WORK

In this paper, we investigate two autonomic parallelism adaptation approaches on a STM system. We examine the performance of different static parallelism and conclude that runtime parallelism adaptation is crucial to performance of TM applications. We introduce feedback control loops to manipulate the parallelism. Followed the description of our models, we compare their performance with static parallelism. We then analyse the implementation overhead and discussed the advantages as well as limitation of our work.

Apart from inappropriate parallelism, thread migration impacts on system performance and cause performance degradation. We plan to investigate the issue and design additional control loops that cooperate with the current loops to control thread affinity and further enhance system performance.

Fig. 2 .

 2 Fig. 2. Periodical profiling procedure. At each decision point (marked by dashed red arrow), the actions are taken. Each decision points corresponds to one state of an automaton as in Fig. 3(b) and Fig. 5.

Fig. 3 .

 3 Fig. 3. The feedback control loop of the simple model. th stands for throughput and tn means the number of thread. The boolean value true means unconditional state transfer.

Fig. 4 .

 4 Fig. 4. Throughput fluctuation. The throughput may continuously rise and descend before reaching its maximum point.

Fig. 5 .

 5 Fig. 5. The controller structure of the probabilistic model described as an automaton.

Fig. 7 .

 7 Fig. 7. The inputs of STAMP

Fig. 8 and

 8 Fig. 8 and Fig. 9 illustrate the execution time comparison with different static parallelism and adaptive parallelism of EigenBench and STAMP. The dots represent the execution time with static parallelism. The solid black line represents execution time with the simple model and the dashed red line gives the execution time with the probabilistic model.According to Fig.8and Fig.9, our adaptive models outperform the performance of the majority of the static parallelism. The probabilistic model shows better performance on applications: genome, vacation, labyrinth against the simple model, but it indicates performance degradation on yada and intruder against the simple model. Both models present similar performance on EigenBench and ssca2. TableIIand TableIIIdetail the performance comparison. The results of both models are compared against the static parallelism which presents best, average and worst performance indicating that our models can outperform the performance of static parallelism if an unknown application is given. The digits in the

Fig. 8 .

 8 Fig. 8. Time comparison of EigenBench on static and adaptive parallelisms. The dots represent the execution time with static parallelism

Fig. 9 .Fig. 10 .

 910 Fig. 9. Time comparison for STAMP on static and adaptive parallelism. The dots represent the execution time with static parallelisms.

Fig. 11 .

 11 Fig. 11. Runtime throughput variation of genome. The red line with dots is the simple model and the black line with crosses is the probabilistic model.

TABLE I .

 I Qualitative summary of each application's runtime transactional characteristics. The classification is based on the application with its optimal parallelism applied.

	Application	TX length	Execution time	Contention
	EigenBench stable	medium	long	medium
	EigenBench online	medium	long	medium
	intruder	short	medium	high
	genome	medium	short	high
	vacation	medium	medium	low
	ssca2	short	short	low
	yada	medium	medium	high
	labyrinth	long	long	low

  different applications from STAMP, namely intruder, ssca2, genome, vacation, yada and

	loops	16667	loops	33333
	A1	35536	A1	145530
	A2	1048576	A2	1048576
	A3	8192	A3	8192
	R1	30	R1	30
	W1	30	W1	30
	R2	20	R2	20
	W2	200	W2	200
	R3i	10	R3i	10
	W3i	30	W3i	30
	R3o	10	R3o	10
	W3o	10	W3o	10
	NOPi 0	NOPi 0
	NOPo 0	NOPo 0
	Ki	1	Ki	1
	Ko	1	Ko	1
	LCT	0	LCT	0
	(a) stable	(b) online variation

  Table II and TableIIIdetail the performance comparison. The results of both models are compared against the static parallelism which presents best, average and worst performance indicating that our models can outperform the performance of static parallelism if an unknown application is given. The digits in the

	execution time (second)
	thread number

TABLE II .

 II Performance comparison of simple model against static parallelism on applications. The higher value, the better performance.

	benchmarks	best case	average	worse case
	EigenBench (stable)	-7% (12)	+10%	+50% (2)
	EigenBench (online variation)	+1% (12)	+17%	+58% (2)
	genome	-57% (4)	+95%	+99% (20)
	vacation	-45% (8)	+79%	+92% (24)
	labyrinth	-52% (24)	+5%	+67% (2)
	yada	-3% (8)	+66%	+91% (22)
	ssca2	-14% (24)	+11%	+62% (2)
	intruder	-6% (6)	+62%	+71% (24)

TABLE III .

 III Performance comparison of probabilistic model against static parallelism on applications. The higher value, the better performance

	benchmarks	best case:	average	worse case
	EigenBench (stable)	-5% (12)	+11%	+51% (2)
	EigenBench (online variation)	-1% (12)	+18%	+57% (2)
	genome	+3% (4)	+97%	+99% (20)
	vacation	-18% (8)	+83%	+93% (24)
	labyrinth	+8% (24)	+42%	+80% (2)
	yada	-17% (8)	+61%	+90% (22)
	ssca2	-16% (24)	+10%	+61% (2)
	intruder	-31% (6)	+53%	+64% (24)

ACKNOWLEDGMENT

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir.


Virginie Daru
email: virginie.daru@limsi.fr

Christian Tenaud
email: christian.tenaud@limsi.fr

Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme

Keywords: High-resolution shock capturing scheme, compressible flow, shock tube problem, shock-wave boundary-layer interaction

. With this scheme, converged results are obtained for Reynolds numbers in the range 200 to 1000. The interaction mechanisms are carefully analyzed as well as the flow dynamics.

Introduction

The interaction of a moving shock wave with a boundary layer is encountered in many flow configurations (buffet, air intakes, over-expanded nozzles, shock tubes). These configurations exhibit unsteady boundary layer detachment over large spatial extents which are responsible for unsteady aerodynamic loads and fluid mixing spoiling the system and dramatically reducing the performance. A fine understanding of the physical mechanisms responsible for unsteadiness, requires the development of very accurate numerical methods for the simulation of these flows, with a large number of grid-points. Therefore, low cost numerical schemes are called to obtain grid-converged solutions at an acceptable CPU time.

A model configuration that exhibits these phenomena corresponds to the shock tube problem. It can be found in various real apparatus, such as in high speed aerodynamic flow facilities used to study the chemical relaxation process in high temperature gas mixtures or to generate high enthalpy reservoir conditions. In flow facility configurations, the free piston reflected shock tunnel is ended by a nozzle communicating with the test section. However, the premature contamination of the flow in the test section, through the production of intense vortices by interaction of the reflected shock-wave and the end wall, limits the usefulness of such devices.

In fact, the reflected shock wave interacts with the incident boundary layer with modified gas properties in this end wall region as a result: a complex unsteady interaction develops between the reflected shock wave and the boundary layer. When the stagnation pressure within the boundary layer is lower than that downstream of the normal reflected shock wave, the boundary layer separates and a lambda-like shock wave pattern forms. This flow configuration has been first analyzed by Mark [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF] using a crude model. Experiments were also carried out allowing the verification of the theoretical model in [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF]. Other detailed experimental analyses of the shock-wave / boundary layer interaction in a shock-tube were presented in [START_REF] Davies | Influence of reflected shock and boundary-layer interaction on shock-tube flows[END_REF][START_REF] Matsuo | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF], leading to enhancements of Mark's model.

The numerical simulation of shock-wave boundary layer interaction in a shocktube was also studied numerically by several authors over the last decade. The first computations were performed by Weber et al. [START_REF] Weber | The numerical simulation of shock bifurcation near the end wall of a shock tube[END_REF] who analyzed the interaction of an isolated reflected shock-wave with a laminar boundary layer; in this work, the contact discontinuity and rarefaction waves present in a shock-tube were disregarded. Other numerical studies focused on the flow contamination at the end-wall of the shock-tube [START_REF] Burtschell | Numerical analysis of reducing driver gas contamination in impulse shock tunnels[END_REF][START_REF] Goozée | Simulation of a complete reflected shock tunnel showing a vortex mechanism for flow contamination[END_REF]. The shock-tube configuration at moderate Reynolds number also served as a test-case for numerical method validation [START_REF] Daru | Evaluation of tvd high resolution schemes for unsteady viscous shocked flows[END_REF][START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF][START_REF] Kim | Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part i: Spatial discretization[END_REF][START_REF] Kim | Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part ii: Multi-dimensional limiting process[END_REF][START_REF] Sjogreen | Grid convergence of high order methods for multiscale complex unsteady viscous compressible flows[END_REF]. The numerical simulation of this type of flow requires numerical schemes which are both robust and highly accurate, to simultaneously represent with high precision the smooth regions of the flow and capture the discontinuities. In a previous study [START_REF] Daru | Evaluation of tvd high resolution schemes for unsteady viscous shocked flows[END_REF], we exhibited the convergence issues encountered when using classical shock capturing TVD approaches on such a test-case. Similar conclusions were also drawn by Sjögreen and Yee [START_REF] Sjogreen | Grid convergence of high order methods for multiscale complex unsteady viscous compressible flows[END_REF] for ACM, TVD or MUSCL approaches. We recently developed an original high order coupled time-space scheme (named OSMP7, accurate up to the 7th-order, in the scalar case) which guaranties Monotonicity Preserving (MP) properties [START_REF] Daru | High resolution monotonicity preserving schemes for unsteady compressible flows[END_REF][START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]. The OSMP7 scheme gives very accurate results on numerous classical test cases of 2-D free supersonic flows [START_REF] Daru | High resolution monotonicity preserving schemes for unsteady compressible flows[END_REF], [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF] and on aeroacoustic problems [START_REF] Daru | Aeroacoustic computations using a high-order shock-capturing scheme[END_REF]. Compared to more classical schemes (Runge-Kutta time integration and ENO or WENO space discretization), the OSMP7 scheme better captures the discontinuities, for comparable order of accuracy in the smooth regions. However, the OSMP7 scheme is six time less expensive than the RK3-WENO schemes in terms of CPU times.

The objective of this paper is two fold. First, we apply the higher-order accurate OSMP7 scheme to the shock tube problem for several Reynolds numbers (Re = 200, 500, 750 and 1000). We have restricted ourselves to this range of Reynolds numbers to limit the grid size requirements. The OSMP7 scheme allows for highly converged results, essentially grid independent, that could not be obtained previously with other schemes [START_REF] Daru | Evaluation of tvd high resolution schemes for unsteady viscous shocked flows[END_REF]. The results are confronted with those of Sjögreen and Yee [START_REF] Sjogreen | Grid convergence of high order methods for multiscale complex unsteady viscous compressible flows[END_REF] for comparable computational meshes. Secondly, we provide a detailed analysis of the shock-wave / boundary-layer interaction and flow dynamics inside the shock tube.

The paper is organized as follows. After a brief review of the governing equations, we present in section (2) the different numerical approaches used in this work: the method of lines (based on Runge-Kutta time integration coupled with a WENO space discretization) and the so-called OSMP (based on coupled time and space integration with MP constraint). The flow field configuration of the shock tube problem is presented in section (3). Section ( 4) is devoted to a study of the influence of the numerical parameters, grid convergence assessment and the influence of the different numerical discretization schemes. The dynamics of the shock wave / boundary layer interaction is analyzed in section [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]. The Reynolds number is the key parameter of the study and its influence on the dynamics is detailed in section [START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF]. Conclusions and comments are summarized in section [START_REF] Davies | Influence of reflected shock and boundary-layer interaction on shock-tube flows[END_REF].

2 Numerical procedure

The set of equations

The governing equations are the compressible Navier-Stokes equations. In conservative form and Cartesian coordinates, they write as:

w t + ∇ • (f (w) -f v (w, ∇w)) = 0, (1)
where w = (ρ, ρu, ρE) t is the vector of the conservative variables, using the classical notations, and f (w) and f v are the Euler and the viscous fluxes re-spectively:

f =         ρ u ρ u ⊗ u + P γM 2 0 I ρ u E + u P γM 2 0         , ( 2
)
f v =        0 σ u • σ + Ψ        , ( 3
)
with the strain rate tensor

σ = µ Re ∇u + ∇ t u - 2 3 ∇ • u I ,
and the heat flux

Ψ = µ (γ -1) Re Pr M 2 0 ∇T.
In addition, a perfect gas law is used:

P γ M 2 0 = (γ -1) ρ E - 1 2 ρ u • u , ( 4
)
T = P ρ , (5)
with ρ the fluid density, u the velocity vector, P the static pressure, T the static temperature, E the total energy per unit of mass and µ the dimensionless dynamic viscosity.

These equations are written in dimensionless form by using the reference values of the density (ρ 0 ), the velocity (v 0 ), and the length scale (L 0 ). The Reynolds number is based on the reference values: Re = ρ 0 v 0 L 0 /µ(T 0 ). The Mach number is M 0 = v 0 / γRT 0 with R = 287 J.Kg -1 .K -1 (for air) the gas constant.

For simplicity, this study considers an ideal gas (air) with constant specific heat ratio γ = 1.4, constant dynamic viscosity µ = 1 and constant Prandtl number Pr = 0.73.

Space and time integrations

The resolution of the Navier-Stokes equations (1, 2, 3) rely on a finite volume approach. The viscous fluxes (3) are discretized by means of a second-order accurate centered scheme, whatever coupled or decoupled time and space integrations are. For the coupled time and space approach, details on this viscousflux discretization can be found in Daru and Tenaud [START_REF] Daru | Evaluation of tvd high resolution schemes for unsteady viscous shocked flows[END_REF]. To check the influence of the viscous flux discretization order on the solution, we also performed few simulations using a fourth-order accurate centered scheme. It was found that the viscous scheme has negligible influence on the resulting solutions (not shown). Therefore, all presented results correspond to the centered 2nd-order viscous flux discretization.

Regarding the convective fluxes, two numerical approximations have been compared in this study. The first one is the high-order coupled time and space scheme proposed in [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF]. The second one, using splitted time and space integrations, is a "Method of lines" approach, based on a multi-stage (Runge-Kutta) time integration. In each sub-stage of the time integration, a high order spatial (WENO) dicretization is applied.

This study is conducted on cartesian meshes. We denote by w n j the discrete quantity w(x, t) estimated at a grid point x j = j δx and at time t = n δt, δt and δx being, respectively, the time and space steps.

Coupled time and space approach

The unlimited One-Step scheme developed in [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF] is of Lax-Wendroff type [START_REF] Lax | Systems of conservation laws[END_REF]. These schemes are constructed in the scalar case by correcting the error terms of the successive modified equations, to yield one additional order of accuracy at each level. For stability reasons, the error term, involving highest degree derivatives at each level, is discretized using upwind formulae for odd derivatives and centred formulae for even derivatives. In this way, one obtains a recurrence relation to construct schemes with arbitrary order of accuracy in time and space. In the non-linear system case, the accuracy order cannot be maintained easily as a Cauchy-Kowalewski procedure is required along the derivation of the successive modified equations. Here, a simpler procedure is used, consisting in a local linearization on the eigenvector basis at the cells' interfaces.

To illustrate the procedure, let us consider a 1-D system of equations w t + ∇ • f (w) = 0, and define λ k and r k the eigenvalues and right eigenvectors of the Roe-averaged The one-step scheme reads:

w n+1 j = w n j - δt δx (F p j+1/2 -F p j-1/2 ) ( 6
)
where F p j+1/2 is the pth-order accurate numerical flux of the scheme at the cell interface (j + 1/2), which can be written:

F p j+1/2 = F Roe j+1/2 + 1 2 k (ψ p r) k,j+1/2 (7)
For clarity, the superscript n has been omitted in the expression of the fluxes. F Roe j+1/2 is the first order Roe flux defined as follows:

F Roe j+1/2 = 1 2 (f j + f j+1 ) - 1 2 k (δ|f |r) k,j+1/2 , (8)
with δ|f | k,j+1/2 = |λ| k,j+1/2 δα k,j+1/2
The function ψ can be decomposed in odd, corresponding to the odd derivatives in the modified equations, and even parts:

ψ p k,j+1/2 = m n=1 ψ 2n k,j+1/2 -js m1 n=1 ψ 2n+1 k,j+1-js/2 (9)
where m = p/2 , m1 = (p -1)/2 ( is the integer division symbol), and jssign(λ k,j+1/2 ). The odd and even ψ functions are given by the recurrence formulae (valid for n ≥ 1):

ψ 2n k,j+1/2 = 2n-2 l=0 (-1) l C l 2n-2 • (c 2n δα) k,j+1/2+n-1-l ( 10
)
ψ 2n+1 k,j+1/2 = 2n-1 l=0 (-1) l C l 2n-1 • (c 2n+1 δα) k,j+1/2+(n-1-l)•js , ( 11
)
where

C s r = r!/[(r-s)!s!].
The coefficients c q depend on the local CFL number, ν k,j+1/2 , and are given by:

(c 2 ) k,j+1/2 = |λ| k,j+1/2 (1 -|ν| k,j+1/2 ) (c q+1 ) k,j+1/2 |ν| k,j+1/2 + (-1) q (q + 1)/2 q + 1 • (c q ) k,j+1/2 , q > 2 . ( 12
)
This scheme has a stencil of only nine points to get seventh order accuracy in both time and space; it is much less than for multi-stage time integration schemes yielding comparable order of accuracy. The unlimited One-Step highorder scheme (called OSp for the p-th order accurate scheme) has a classical CFL stability condition 0 ≤ |ν k | ≤ 1. Note that this type of scheme yields the exact solution for CFL number equal to 1, in the scalar case.

TVD-MP constraints: An ad hoc discontinuity-capturing feature must be employed to limit the spurious oscillations in the vicinity of the strong gradient regions. Total Variation Diminishing (TVD) schemes are generally considered to be well suited for the capture of discontinuities. Nevertheless, TVD constraints are known to clip the extrema which appears as a serious drawback in a shock-capturing feature. To avoid this loss of accuracy near extrema, it is necessary to satisfy the Monotonicity-Preserving (MP) criteria introduced by Suresh et al. [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF], that enlarge the TVD intervals to provide room for the numerical flux to maintain an accurate value. In Daru et al. [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF], we derived a Monotonicity Preserving version of our one-step p-th order scheme (named OSMPp) that preserves accuracy near extrema. In particular, we recast the original MP constraints of Suresh et al. [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF] in the TVD framework which allows us to generalize the MP conditions in terms of flux limiting, and avoid any CFL restriction. The MP conditions that preserves accuracy can be expressed directly as constraints on the ψ p function [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF].

Let us define the function φ p such that :

φ p k = ψ p k (1 -|ν| k )δ|f | k ( 13
)
for each characteristic field. The numerical flux of the OSp scheme can be written, using φ p , as:

F p j+1/2 = F Roe j+1/2 + 1 2 k (φ p (1 -|ν|)δ|f |r) k,j+1/2 (14)
If we define : [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF] the TVD conditions derived by Harten can be expressed as a constraint on φ p :

φ 1 k,j+1/2 = 2r k,j+1/2 |ν| k,j+1/2 • 1 -|ν| k,j+1/2-js 1 -|ν| k,j+1/2 φ 2 k,j+1/2 = 2 1 -|ν| k,j+1/2
φ p,T V D k,j+1/2 = max 0, min φ 1 k,j+1/2 , φ p k,j+1/2 , φ 2 k,j+1/2 (16) with r k,j+1/2 = δα k,j+1/2-js δα k,j+1/2 ( 17
)
The MP constraints are based on the principle that the TVD conditions should be activated only in monotone regions. The following flux variations can be defined :

                   δF p k,j+1/2 = 1 2 ψ p k,j+1/2 δf ul k,j+1/2 = 1 2 φ 1 k,j+1/2 • (1 -|ν| k,j+1/2 )δ|f | k,j+1/2 δf md k,j+1/2 = 1 2 δ|f | k,j+1/2 -1 2 d k,j+1/2 δf lc k,j+1/2 = 1 2 δf ul k,j+1/2 + 1 2 1-|ν| k,j+1/2-js |ν| k,j+1/2 d k,j+1/2-js (18) with d k,j+1/2 = minmod(d k,j , d k,j+1 ), where d k,j = (λδα) k,j+1/2 -(λδα) k,j-1/2 is a curvature term.
The TVD conditions constrain δF p k,j+1/2 to be inside the interval [0, δ|f | k,j+1/2 ] [0, δf ul k,j+1/2 ] (this is equivalent to the TVD constraint ( 16)). The MP conditions provide an enlargement of these intervals : the TVD constraint ( 16) is activated only when

δF p k,j+1/2 is outside the intervals [0, δ|f | k,j+1/2 , δf md k,j+1/2 ] [0, δf ul k,j+1/2 , δf lc k,j+1/2 ], where we noted by [f 1 , f 2 , ..., f k ] the interval [min(f 1 , f 2 , ..., f k ), max(f 1 , f 2 , ..., f k )].
In practice, as was done in [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF], the criteria is strengthened by choosing the following definition for the curvature measurement d k,j+1/2 :

d k,j+1/2 = d M 4 k,j+1/2 = minmod (4d k,j -d k,j+1 , 4d k,j+1 -d k,j , d k,j , d k,j+1 )(19)
The resulting scheme (OSMPp) will be pth-order accurate almost everywhere, except around discontinuities where it becomes first order accurate as is the case of all TVD schemes.

The multidimensional extension is even more delicate than in the case of separate time-space discretization (which is straightforward when treated dimension by dimension [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF]), since we have to consider cross derivative terms that appear in the second and higher order terms, which are left uncontrolled if one applies a direction by direction MP correction to a Lax-Wendroff unsplit scheme. A Locally Extremum Diminishing scheme can be obtained if the mixed terms are discretized by using upwind formulae, but is very difficult to implement. The simplest way to avoid this problem of cross derivatives and to recover good properties of the one-dimensional scheme is to use a Strang directional splitting strategy [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF], which is only second order accurate. In two dimensions, the system of equations is written:

∂w ∂t + ∂f (w) ∂x + ∂g(w) ∂y = 0 ( 20
)
where f and g are the fluxes in each direction. We implement the splitting as follows:

w n+2 j = L δx L δy L δy L δx w n j (21) L δx (resp. L δy ) being a discrete approximation of L x (w) = -f (w) x (resp. L y (w) = -g(w) y ).
In such a way, the second order accuracy is recovered every two time steps.

While reducing the order of accuracy of the original scheme, this simple extension to multidimensional non-linear systems, proposed here, provides a scheme that can compete very well with high order WENO schemes at low cost, as we will see hereafter.

Method of lines

Besides the high-order coupled approach, a separate time and space integrations has also been used to solve the governing equations (1, 2).

Temporal integration

The time integration is performed by means of a third order Runge-Kutta method [START_REF] Shu | High order eno schemes applied to two and three dimensional compressible flows[END_REF]:

w 0 = w n w 1 = w 0 + δtL (w 0 ) w 2 = 3 4 w 0 + 1 4 w 1 + 1 4 δtL (w 1 ) (22)
w (n+1) = 1 3 w 0 + 2 3 w 2 + 2 3 δtL (w 2 )
where

L(w) is a discrete approximation of L(w) = -∇ • (f (w) -f v (w, ∇w)) .
This time integration has been chosen since it does not increase the Total Variation of L(w).

Space discretization

The integration of the viscous fluxes is performed by means of a second order centered approximation.

The convective terms are treated following a WENO procedure based on Roe's approximate Riemann solver [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF]. The WENO scheme is based on a weighted (convex) combination of the p stencil candidates of the original ENO scheme [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF], [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes, II[END_REF], to increase the order of accuracy of the latter. The weights [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF] depend on the degree of regularity of the solution; i.e. in regular regions, the weight coefficients allow one to achieve a (2p -1) th -order of accuracy while in regions where discontinuities occur, they are set to zero leading to a standard ENO scheme. The WENO fluxes are estimated by:

F W EN O j+1/2 = k   p-1 l=0 ω (p) l q p l (f j+l-p+1 , ...., f j+l ) r k   ( 23
)
where q p l is the ENO reconstruction on the l th stencil, using the coefficients (ζ p l,k ) given in Table 1:

q p l (f j+l-p+1 , ..., f j+l ) = p-1 k=0 ζ p l,k f k , ( 24
)
and ω

(p) l are the weights defined by:

p-1 l=0 ω (p) l = 1. ( 25
)
ω l = β l p-1 k=0 β k (26) β l = C p l (ε + IS l ) 2 . ( 27
)
ε is a small positive number to avoid denominator to be zero (hereafter we set ε = 10 -6 ) and IS l is a regularity measurement of the flux on the l th stencil candidate. The evaluation of the smoothness measurement (IS l ) is based on the undivided-differences [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF].

Hereafter, the calculations were performed by using p = 5 that gives a theoretical 9 th -order of accuracy in space. The corresponding ENO coefficients are

p l k = 0 k = 1 k = 2 k = 3 k = 4

Boundary conditions

On the solid wall boundaries, a no-slip condition is prescribed on the velocity components. The density is calculated through the use of the continuity equation written on the wall points. An adiabatic condition is prescribed on the wall temperature. A standard fourth order upwind formula is used to estimate the normal to the wall derivatives.

Flow field configuration

We consider a square shock tube with a unit side length and insulated walls. The diaphragm is initially located at the middle of the tube (x = 0.5). The initial states, in terms of dimensionless quantities, are on the left-side of the diaphragm: ρ L = 120, P L = ρ L /γ, u L = v L = 0, and on the right-side

: ρ R = 1.2, P R = ρ R /γ, u R = v R = 0.
The reference velocity is based on the initial speed of sound, corresponding to a reference Mach number M 0 = 1. The Prandtl number is assumed constant (Pr = 0.73, for air). Four Reynolds numbers have been considered: Re = 200, 500, 750, 1000. At the initial time, the diaphragm is broken. The 1-D inviscid solution is presented in Fig. 1, showing the evolution of the density (left) and the velocity (right) in the x -t plane. A shock wave, followed by a contact discontinuity, moves toward the low-pressure region (i.e. to the right, in the present case), while a sonic rarefaction wave moves toward the high-pressure region at the left. The shock Mach number is equal to 2.37. The rarefaction wave reflects on the left wall at a dimensionless time t = 0.5. The incident shock wave is weak, and reflects at the right wall at time t ≈ 0.21. After reflection, it interacts with the contact discontinuity.

Complex interactions then occur leading to transmitted and reflected waves as one can see in figure 1. For later times, the contact discontinuity stays stationary, close to the right wall. On the contrary, the reflected shock wave moves on and begins to interact with the rarefaction wave at time t = 0.4, and later with the reflected rarefaction wave. During the interaction the reflected shock has an approximately constant velocity, but the fluid velocity decreases continuously after its interaction with the incident rarefaction wave (Fig. 1 right). In the 2D viscous case, the propagating incident shock wave and contact discontinuity interact with the horizontal wall, creating a thin boundary layer. After its reflection on the right wall, the shock wave interacts with this boundary layer. As the stagnation pressure in the boundary layer is lower than the one within the outflow region, a separation region (named "bubble") appears over a large extent within the boundary layer, resulting in a major modification of the flow pattern and the formation of "a lambda-shape like shock pattern", as one can see in figure 6. The triple point emerging from the lambda-shape like shock pattern generates a slip line that rolls up in the right end corner.

Grid convergence study and influence of the schemes

The simulations have been carried out for four Reynolds numbers: Re = 200, 500, 750 and 1000. Simulations with increasing grid resolutions have been performed to verify grid-independency of the solutions. All meshes use equally spaced points with cell aspect ratio equal to unity. The grid convergence is shown for the density distribution along the bottom wall at t = 1 (Fig. 2). To make sure that other quantities are well converged, we also provide the distributions of the skin friction coefficient along the bottom-wall for the four Reynolds numbers (Fig. 3). Furthermore, for Re = 1000, we plot the den- sity (Fig. 4-a) and the vorticity (Fig. 4-b) along the horizontal line crossing the separation bubble and the vortices, located at y = 0.05. Simulations for RK3-WENO5 scheme on the finest grid (4000 × 2000) at Re = 1000 were not performed due to prohibitive computational cost.

As expected, greater Reynolds numbers require finer grids to achieve convergence. From this numerical experiment, we found that sufficient grid resolutions for the different Reynolds numbers are: 1000×500 for Re = 200 (Fig. 2 In a previous study [START_REF] Daru | Evaluation of tvd high resolution schemes for unsteady viscous shocked flows[END_REF], we have shown convergence issues when using classical shock capturing TVD approaches for this test-case, especially at Re = 1000. Subsequently, Sjögreen and Yee [START_REF] Sjogreen | Grid convergence of high order methods for multiscale complex unsteady viscous compressible flows[END_REF] have unsuccessfully attempted to obtain converged solutions by using several separate time and space schemes (MUSCL-RK2, WENO5-RK4 and high order filtered difference methods). For instance, the solution for high order filtered difference methods is notably different than the solution for MUSCL-RK2 method. Also, for Re = 1000, they obtained solutions significantly different than ours. They mentioned that one cannot conclude on the respective accuracy of the various schemes and that the stability of the flow is questionable. In our simulations, converge solutions are clearly obtained as seen in previous figures and also in figure 5 where plotted are the density fields at t = 1 for several grids. Again, results for RK3-WENO5 scheme on the finest grid (4000 × 2000) are lacking in figure 5 because of the prohibitively large CPU time required, but it is obvious that the two numerical methods tend to converge to the same solution.

Computations were performed on a NEC-SX8, at the CNRS's National Computing Center (IDRIS). The OSMP7 scheme takes about 3.10 -7 s.CPU per time step and grid point, at a flow-rate of 7 Giga Flops. The CPU time required by the RK3-WENO5, at a flow-rate of 6 Giga Flops, is 2.10 -6 s.CPU per time step and grid point ie six to seven times greater than the OSMP7 scheme. We can observe that the OSMP7 scheme exhibits a faster grid convergence than the RK3-WENO5 scheme on the vortical structure in the bottom-right corner due to the roll-up of the contact discontinuity. This can mainly be attributed to the poor behavior of the WENO scheme for contact discontinuities. However, the convergence of the main vortical structure (downstream of the lambda-like shock pattern) seems to be rather identical for both schemes. Though very small differences are seen between the 3000×1500 and 4000×2000 mesh solutions, we can consider that the solutions on the 3000 × 1500 grid are almost converged for the two schemes. This support the claim of a converged solution for the 4000 × 2000 mesh and OSMP7 scheme, for the highest considered Reynolds number (Re = 1000). We suggest that our converged solutions could be used as a reference in future work.

5 Dynamic flow analysis of a shock wave / boundary layer interaction.

This section is devoted to the analysis of the flow dynamics during the interaction of the reflected shock-wave with the laminar boundary layer. The interaction begins at the bottom-right corner right after the incident shock wave has reflected on the right wall. When the reflected shock interacts with the bottom wall boundary layer, the stagnation pressure within the boundary layer is lower than downstream of the normal shock. After the shock reflection, the hotter fluid in contact with the right wall moves suddenly down to the bottom wall, creating a hot countercurrent jet close in the corner. A recirculating flow is then generated that initiates the boundary layer detachment to form a "bubble". The deviation of the cold flow over the slip line (supersonic shear layer) delimiting the hot fluid "bubble" then induces a lambda-shape like shock pattern. This early stage of the reflected shock/boundary interaction can be seen at a dimensionless time t = 0.4, on figure [START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF] where the static temperature contours are plotted for the solution at Re = 200 obtained on the finest grid. A short time later, the shock interacts with the incident contact discontinuity. The incident flow in between the lambda-shock and the wall jet, rolls up around the bubble and forms a high-speed cold jet. After this jet has impinged on the bottom wall (around t = 0.5 as we can see in figure ( 6)), no more hot shocked fluid enters the bubble, which is isolated from the hot right wall region. The countercurrent part of the jet then alternatively hits and reflects between the bottom wall and the supersonic shear layer originating from the separation point (Fig. 6), thus creating a series of counter-rotating vortices. Each vortex is made up of three "layers" of fluid : at the center the hottest shocked fluid, next a layer of colder fluid originating from the boundary layer, and finally the coldest layer of fluid from the external flow at the periphery (see Fig. 6). The jet flow is very energetic with supersonic velocities at several places. Each impact of the jet, either on the solid wall or on the oblique slip line, creates a relatively high pressure region. This leads to the deformation of the slip line and to the generation of secondary shock waves above it.

6 Influence of Reynolds number on the flow dynamics.

The analysis of the flow dynamics within the separation region has pointed out that the speed of the counter flow jet increases with the Reynolds number. Consequently, during a fixed period of time, the number of jet alternated reflections between the bottom wall and the slip line increases with Re. Also, the number of counter rotating vortices produced in this interaction region increases too. This can be seen from the density gradient and temperature The Reynolds number also affects the location of the triple point of the lambda-shape like shock pattern. One can observe, in figures 7 and 8, that the height of the triple point increases with Re. To better appreciate the motion of the triple point, its coordinates are plotted as a function of time in Fig. 9 for the four Reynolds numbers. Its horizontal velocity is constant and nearly independent of Re (Fig. 9-left). This velocity is slightly greater than that of the normal shock along the tube centerline (y = 0.5). On the contrary, the triple point motion along the vertical direction depends on the Reynolds number (Fig. 9-right). Time evolutions of the triple point distance to the bottom wall, are in agreement with theoretical predictions (e.g. see [START_REF] Davies | Influence of reflected shock and boundary-layer interaction on shock-tube flows[END_REF], [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF] and [START_REF] Matsuo | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF]). Similarly to Matsuo et al. [START_REF] Matsuo | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF], we found two periods before and after t = 0.53 where the Re dependency changes. Surprisingly, at t = 0.53 all the curves cross each other. In the first period (before t = 0.53), the triple point height is lower for higher Re due to the thinner boundary layer. We found that t = 0.53 corresponds to the time when the bubble is "closed", i.e. when the heated vertical jet, coming from the fluid flowing in between the lambda-shock and rolling up around the bubble, impinges on the bottom wall and separates the right wall hot region from the bubble. We have no explanation to propose for the Re independency of the closing time of the bubble. Anyway, before t = 0.53, the separated bubble is fed by both the incoming boundary layer and the heated vertical jet. On the contrary, after that time, the detachment is only fed by a portion of the boundary layer, the bubble being "closed". In accordance with [START_REF] Davies | Influence of reflected shock and boundary-layer interaction on shock-tube flows[END_REF], [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF] and [START_REF] Matsuo | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF], we verified in our results that the pressure downstream the leading oblique shock wave (under the lambda-shock) is the same than the pressure at the beginning of the detachment which is equal to the stagnation pressure in the incoming boundary layer. In a crude theory, Mark [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF] showed that the stagnation pressure in the boundary layer is only related to the Mach number of the reflected shock wave and is independent of the Reynolds number. In the flow condition of our test-case, Marks's theory was corroborated by more sophisticated developments undertaken by Davies et al. [START_REF] Davies | Influence of reflected shock and boundary-layer interaction on shock-tube flows[END_REF] and Matsuo et al. [START_REF] Matsuo | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF]. The ratio between the boundary layer stagnation pressure and the static pressure (above the boundary layer) upstream of the leading reflected shock can be determined using the Mach number of the reflected shock wave [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF] (in the present case M s = 2.37). This Mach number being independent of Re, it explains the constant angle of the leading oblique shock wave for all Re found in the simulations. However, the velocity of the jet increases with Re leading to increasing bubble size, as predicted by Mark [START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF].

Since the angle of the leading oblique shock is constant, it is shifted to the left with increasing distance from the bottom wall of the triple point for increasing Re (Fig. 9-right).

The Reynolds number effects also impacts the vortices created within the lambda shock pattern, i.e. at t = 0.6 (Fig. 7): three vortices are created at Re = 200 and six at Re = 1000. When the number of vortices increases, additional shocks are created inside the lambda shock pattern. To further analyze the dynamical process, we have also plotted, in a coordinate reference frame attached to the reflected normal shock, the vorticity field with some streamlines superimposed. These plots are presented at several dimensionless times and for the four Re studied (Fig. [START_REF] Kim | Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part i: Spatial discretization[END_REF][START_REF] Kim | Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part ii: Multi-dimensional limiting process[END_REF][START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Mark | The interaction of a reflected shock wave with the boundary layer in a shock tube[END_REF]. For the lowest Reynolds numbers, Re = 200 and Re = 500, Fig. 10 and 11 show a regular arrangement of the vortices inside the bubble. This regular arrangement disappears for Re ≥ 750 (Fig. 12 and13). Detailed examination of the dynamics of the flow reveals the following process : eddy structures are created close to the bottom wall, with an increasing enstrophy for increasing Re (Fig. 14). Similar to the vortex/boundary layer interaction, when the enstrophy of the eddy structures is large enough, they are ejected out of the boundary layer. Hence, the structures (mainly the second one from the right as indicated on figures 12-b and 13-b). Due to the shear stress of the outflow, the filament is elongated to the point where it breaks up. After that, the enstrophy intensity decreases close to the wall, resulting in the retraction of the structure with a "spring back effect". This in turn affects the vortex arrangement which becomes irregular, with flattened structures. Neighboring co-rotating structures eventually merge. This study demonstrates that the calculation of such a complex flow requires numerical schemes which are both robust and very accurate. The OSMP7 scheme allows to reach fully converged solutions for Reynolds numbers as high as 1000. Simulations with sufficiently refined grids to obtain converged solutions could not be performed using the classical WENO scheme because of the excessive CPU cost. Using converged solutions, the flow dynamics and its evolution with the Reynolds number were analyzed in detail. The present analysis is only valid for 2-D laminar configurations. For larger Re, 3-D instabilities will change the present flow dynamics. Nevertheless, by considering the Re range we studied, this 2D shock tube problem is an appropriate test-case for validation of high-resolution schemes, and the presented grid-converged results may serve as reference solutions for that purpose.

  l = 1 l = 2 l = 3 l = 4 p = 5 1/126 10/63 10/21 20/63 5/126

Fig. 1 .

 1 Fig. 1. x -t diagram of the density (left), velocity (right).

Fig. 2 .

 2 Fig. 2. Density distributions along the bottom wall of the shock tube, obtained at t = 1 on several fine grids by using the OSMP7 scheme (solid lines) and the RK3-WENO5 scheme (dashed lines), for four Reynolds numbers: (a) Re = 200, (b) Re = 500, (c) Re = 750 and (d) Re = 1000.

Fig. 3 .

 3 Fig. 3. Distributions of the skin friction coefficient along the bottom wall of the shock tube, obtained at t = 1 on several fine grids by using the OSMP7 scheme (solid lines) and the RK3-WENO5 scheme (dashed lines), for four Reynolds numbers: (a) Re = 200, (b) Re = 500, (c) Re = 750 and (d) = 1000.

Fig. 4 .

 4 Fig. 4. Distributions of density (a) and the vorticity (b) along the line y = 0.05, obtained at t = 1 on several fine grids by using the OSMP7 scheme (solid lines) and the RK3-WENO5 scheme (dashed lines), for Re = 1000. and 3-a), 1500×750 for Re = 500 (Fig. 2-b and 3-b), 2000×1000 for Re = 750 (Fig. 2-c and 3-c) and 4000 × 2000 for Re = 1000 (Fig. 2-d, 3-d and 4).

Fig. 5 .

 5 Fig. 5. Contours of ρ obtained on the finest grid at t = 1 for Re = 1000: (a) OSMP7 (4000 × 2000); (b) RK3-WENO5 (3000 × 1500); (c) OSMP7 (3000 × 1500); (d) RK3-WENO5 (2000 × 1000); (e) OSMP7 (2000 × 1000). (21 contour levels between 20 and 120).

Fig. 6 .

 6 Fig. 6. Contours of T obtained on the finest grid for Re = 200 at (a) t = 0.4, (b) t = 0.6, (c) t = 0.7, (d) t = 0.8. Hot fluid in dark, cold fluid in light. 41 niveaux entre 0.4 et 1.2fields in Fig.7 and 8, where the solutions for the four Reynolds numbers are presented at t = 0.6 and t = 1 respectively. Shocklets in the vicinity of the vortices as well as Kelvin-Helmholtz instabilities in the supersonic shear layer originating at the triple point are visible on the contours of density gradients (Fig.7-left and 8-left). The temperature fields (on Fig.7-right and 8-right) give additional information on the flow, notably on the location of hot and cold streams and the flow coming from the initial boundary layer.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Contours of |∇ρ| (left) and T (right, 41 levels entre 0.4 et 1.2) obtained at t = 0.6 on the finest grid for several Reynolds numbers: (a) Re = 200; (b) Re = 500; (c) Re750; (d) Re = 1000.13-b) are stretched and distorted into an intense vortex filament (shown on figures 12-b and 13-b). This filament is then ejected out of the boundary layer and trapped, through a pairing process, with the largest vortical structure (especially, the first structure from the right, as indicated on figures 12-b and

Fig. 9 .

 9 Fig. 9. x (left) and y (right) coordinates of the triple point versus time 7 Conclusion

Fig. 10 .

 10 Fig. 10. Contours of the vorticity field with streamlines superimposed obtained, in a coordinate referential attached to the reflected normal shock, for Re = 200 at time (a) t = 0.6, (b) t = 0.9 and (c) t = 1.

Fig. 11 .

 11 Fig. 11. Contours of the vorticity field with streamlines superimposed obtained, in a coordinate frame attached to the reflected normal for Re = 500 at time (a) t = 0.6, (b) t = 0.9, (c) t = 1.

Fig. 12 .

 12 Fig. 12. Contours of the vorticity field with streamlines superimposed obtained, in a coordinate frame attached to the reflected normal shock, for Re = 750 at time (a) t = 0.6, (b) t = 0.9, (c) t = 1.

Fig. 13 .

 13 Fig. 13. Contours of the vorticity field with streamlines superimposed obtained, in a coordinate frame attached to the reflected normal shock, for Re = 1000 at time (a) t = 0.6, (b) t = 0.9, (c) t = 1.

Fig. 14 .

 14 Fig. 14. Contours of enstrophy density ( 1 2 |w 2 |) obtained at time = 0.9 for (a) Re = 200, (b) Re = 500, (c) Re = 750 and (d) Re = 1000.

Table 2

 2

	Weights (ω	(p)

l ) for the WENO (p=5) procedure reported in Table

[START_REF] Burtschell | Numerical analysis of reducing driver gas contamination in impulse shock tunnels[END_REF]

while the weights are given in Table

[START_REF] Daru | Aeroacoustic computations using a high-order shock-capturing scheme[END_REF]

.

Acknowledgements:

Computations have been carried out on NEC-SX8 of the CNRS's national supercomputing center (IDRIS). The authors greatly acknowledge the support of this institution.