Yann Barsamian
email: ybarsamian@unistra.fr

Sever A Hirstoaga
email: sever.hirstoaga@inria.fr

Éric Violard
email: violard@unistra.fr

Efficient Data Layouts for a Three-Dimensional Electrostatic Particle-in-Cell Code

Keywords: data structures, space-filling curves, SIMD architecture, hybrid parallelism, strong and weak scaling, three-dimensional Particle-in-Cell simulation, plasma physics

The Particle-in-Cell (PIC) method is a widely used tool in plasma physics.

To accurately solve realistic problems, the method requires to use trillions of particles and therefore, there is a strong demand for high performance code on modern architectures. The present work describes performance results of Pic-Vert, a hybrid OpenMP/MPI and vectorized three-dimensional electrostatic PIC code.

The code simulates 3d3v Vlasov-Poisson systems on Cartesian grids with periodic boundary conditions. Overall, it processes 590 million particles/second

on a 24-core Intel Skylake architecture, without hyper-threading (25 million particles per second per core).

The paper presents extensions in 3d of our preliminary 2d results [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF], with highlights on the difficulties and solutions proposed for these extensions. Specifically, our main contributions consist in proposing a new space-filling curve in 3d (called L6D) to improve the cache reuse and an adapted loop transformation (strip-mining) to achieve efficient vectorization. The analysis of these optimization strategies is performed in two-stages, first on a 24-core socket and second on a super-computer, from 1 to 3,072 cores, demonstrating significant performance gains and very satisfactory weak scaling results of the code.

Introduction

Implementing a Particle-in-Cell (PIC) method is an important step of many applications in the field of computational science. Although particle methods are usually characterized by low accuracy, they are able, when using a very large number of particles, to adequately reproduce complicated phenomena, for instance in plasma physics [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF]. However, in such a situation, a naive implementation can quickly exhibit memory bottlenecks since the overall cost is dominated by data motion and not by computation [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF].

Therefore, a lot of research efforts were recently devoted towards more adapted PIC implementations that utilize efficiently modern super-computing resources [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF][START_REF] Decyk | Optimization of particle-in-cell codes on reduced instruction set computer processors[END_REF][START_REF] Vincenti | An efficient and portable SIMD algorithm for charge/current deposition in particle-in-cell codes[END_REF][START_REF] Germaschewski | The plasma simulation code: A modern particle-in-cell code with patch-based load-balancing[END_REF][START_REF] Decyk | Particle-in-cell algorithms for emerging computer architectures[END_REF][START_REF] Jocksch | A bucket sort algorithm for the particle-in-cell method on manycore architectures[END_REF]. Despite these significant advances, targeting realistic applications in plasma physics is still challenging, due to the complex multi-scale six-dimensional problems to be solved in the phase space. In this direction, our objective is to develop a massively parallel and highly scalable PIC code, in view of physically meaningful realistic simulations.

First steps were achieved in [START_REF] Chacon-Golcher | Optimization of particlein-cell simulations for Vlasov-Poisson system with strong magnetic field[END_REF][START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] where we built a hybrid parallel and vectorized PIC code for a Vlasov-Poisson model in a two-dimensional (2d) physical space. Several data structures for particles and for grid quantities were analyzed in order to enhance data locality and to reduce the execution time with respect to a classic code. More precisely, we showed in [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] that a structure of arrays and a L4D space-filling curve lead to an efficient ordering of the particles and of the electric field and the charge density, respectively. In addition, the performance of parallelization of the loops through distributed and shared memory paradigms was assessed in tandem with the memory channels.

In the present contribution we extend to three-dimensions (3d) for the physical space the PIC code in [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] for simulating electrostatic plasma. Even though reduced 2d models are used in the literature to gain insights about the main behavior of the plasma, full 3d simulations clearly improve the realism of the physical description. Moreover, in some situations, reducing the dimensionality is not even possible and thus a 3d simulation is unavoidable. However, with the aim of keeping a satisfactory accuracy of the computed solutions in this case, we need to consider the significant increase in the amount of data to be processed.

A first obvious difficulty of the extension from 2d to 3d simulations is the need for more storage for the grid quantities [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF]. If, in a 2d simulation, keeping the whole grid quantities in the cache memory is still possible, this target is much more difficult or even impossible to achieve for a fine grid in the 3d case. Secondly, a 3d simulation requires more data traffic and computations. Specifically, passing from 4 grid points to 8 points for the interpolation and accumulation steps (when using a linear approach) becomes more difficult to handle, both for the data flow and for the vectorization. In addition, the difficulty increases significantly if higher order approximations within these steps are in use. The same challenges need to be addressed when vectorization is used for the updating positions step.

The paper is organized as follows: in Section 2 we present the basic kinetic model for the plasma, we detail the steps of the PIC implementation, introduce the related work and explain our contributions. In Section 3 we detail the code optimizations and we present the performance results on 24 cores with OpenMP only. In Section 4 we show the scalability of the code on up to 3,072 cores of the supercomputer Marconi. Section 5 summarizes the work and presents some future directions.

PIC overview for the 3d Vlasov-Poisson model

Description of the problem

A PIC method simulates a plasma by integrating self-consistently the trajectories of charged particles with fields that are generated by the particles themselves [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF]. In the case where there is no other external field and the self-consistent magnetic field is neglected, this relies on solving the following Vlasov-Poisson system:

               ∂ t f + v • ∇ x f + q m E • ∇ v f = 0 Vlasov f (x, v, 0) = f 0 -∆φ = ρ 0
Poisson where ρ(x, t) = q f (x, v, t)dv and E(x, t) = -∇φ(x, t).

In the system above, f = f (x, v, t) stands for the distribution of one species of particles (with charge q and mass m) in a six-dimensional phase space (three dimensions for positions and three dimensions for velocities), ρ stands for the charge density, and E for the self-consistent electric field. A PIC method consists in discretizing (sampling) the distribution function by a collection of macroparticles that move in the phase space following the characteristics of the Vlasov equation. Then, a PIC simulation follows four steps: accumulate on the spatial grid the particle charge, solve the Poisson equation to obtain the grid electric field, interpolate this field to the particles, and finally push in time the particle positions and velocities.

In our PIC code, the particle positions and velocities are initialized randomly using the WELL generator [START_REF] Panneton | Improved long-period generators based on linear recurrences modulo 2[END_REF]. The particles all have the same fixed weight. They are advanced in time with a leap-frog time stepping [2, Section 2.4] (second-order in time). The electric field is computed by solving the Poisson equation on a uniform Cartesian grid, by a Fourier method. Then, for the particle and force weighting we use the Cloud-in-Cell model [START_REF] Birdsall | Clouds-in-clouds, clouds-in-cells physics for manybody plasma simulation[END_REF] (first-order in space), meaning that eight neighboring grid points are used in the interpolation/accumulation steps for a particle in that cell. The PIC pseudo-code is detailed in Fig. 1.

Next we describe the two basic types of data on which the sequential implementation of the code is based.

Particle Data Structure

A particle is given by its position and velocity. While the velocity is classi- Then, we consider the integers i x = floor(x), i y = floor(y) and i z = floor(z),

and the normalized offsets (which are real numbers in [0; 1))

dx = x -i x , dy = y -i y and dz = z -i z . (2)
The cell index i cell is a number in {0, 1, . . . , ncx • ncy • ncz -1}, taken as the image of some one-to-one mapping depending on (i x , i y , i z ). For example, the row-major mapping

(i x , i y , i z ) → i cell = (i x • ncy + i y ) • ncz + i z i cell →          i x = i cell ncz•ncy i y = mod i cell ncz , ncy i z = mod(i cell , ncz) (3)
is commonly used in C. However, several bijection mappings will be analyzed in the following sections with the aim of improving the cache reuse.

Finally, a particle is stored in memory with 1 int (i cell ), 3 floats (dx, dy and dz) and 3 doubles (vx, vy and vz). In our code, these 7 attributes are stored in a Structure of Arrays (SoA), to enable efficient vectorization. In this framework, the update particle position step (line 10 in Fig. 1) can be accomplished in the following four sub-steps: Update positions 1 Compute (i x , i y , i z ) from i cell .

2 Update (x, y, z) using formula (2) and line 10 in Fig. 1.

3 Compute the new values of (i x , dx, i y , dy, i z , dz) using formulas (1) and (2).

4 Compute i cell from (i x , i y , i z ).

The reason behind this approach stems from the fact that on modern architectures, it might be faster to make rapid computations than to store numbers. In this case, we clearly need fast algorithms for computing the bijection functions in sub-steps 1 and 4 above. For example, the bijection mapping in (3) can be computed very fast in both directions, which is not the case for all the mappings analyzed in this paper. In Section 3.2 we show for some bijection mappings that computing them is faster than storing the integers (i x , i y , i z ) in addition to i cell .

E and ρ Data Structure

The standard 3d representation of the electric field E and of the charge density ρ stores their values at the grid points:

double Ex[ncx][ncy][ncz], Ey[ncx][ncy][ncz], Ez[ncx][ncy][ncz]; double rho[ncx][ncy][ncz];
In this case, the interpolation of the 3d arrays Ex, Ey and Ez requires access to memory locations that are not contiguous. A solution to partially overcome this problem consists in storing components of the field in only one array [5, Section IV, Case 3]:

double Exyz[ncx][ncy][ncz][3];
Unfortunately, this data structure still leads to non-contiguous accesses. This problem is solved by using a redundant one-dimensional array of coefficients [START_REF] Bowers | Speed optimal implementation of a fully relativistic particle push with charge conserving current accumulation on modern processors[END_REF][START_REF] Chacon-Golcher | Optimization of particlein-cell simulations for Vlasov-Poisson system with strong magnetic field[END_REF]:

double E_1d[ncx*ncy*ncz][24]; double rho_1d[ncx*ncy*ncz][8];
The redundant array E 1d stores for each cell the values of each of the three field arrays at the grid points on the eight corners of the cell, contiguously in memory. The array rho 1d similarly stores for each cell the values of the charge density on the corners of that cell. They take eight times more memory than the standard layout, but it has been shown that, for the charge density, it is more efficient because it opens the possibility to vectorize the accumulation step (line 11 in Fig. 1) [6, Section 4.1.2.]. We show in the sequel how to gain performance through cache hit improvements, both for E and ρ.

Related Work

In order to scan E and ρ as locally as possible, the particles are sorted by cell index periodically in time (the step in lines 4-5 in Fig. 1) [START_REF] Decyk | Optimization of particle-in-cell codes on reduced instruction set computer processors[END_REF]Section VI]. The cost of this step is made linear in N by using the counting sort [14, Section 8.2], since the number of cells is much smaller than the number of particles.

Space-filling curves can be used to enhance cache performances. Precise results were shown on applications with regular memory accesses such as linear algebra [START_REF] Chatterjee | Nonlinear array layouts for hierarchical memory systems[END_REF], or with irregular accesses such as the n-body problem [START_REF] Mellor-Crummey | Improving memory hierarchy performance for irregular applications using data and computation reorderings[END_REF]. Nevertheless, none of these results can directly apply to a PIC code, which has irregular accesses over the arrays E and ρ and not the particle array itself, as in the n-body problem. In the previous contribution [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] we showed how spacefilling curves can enhance the performance of a 2d PIC code. The present study is an extension to 3d of this work.

The space-filling curves are also of interest in particle codes at the interprocess level, to achieve load balancing when using domain decomposition [START_REF] Germaschewski | The plasma simulation code: A modern particle-in-cell code with patch-based load-balancing[END_REF] or to minimize communication between processes [START_REF] Deford | Empirical analysis of space-filling curves for scientific computing applications[END_REF].

A rather popular technique on modern PIC codes is to organize the particles by so-called super-cells, e.g. PIConGPU [START_REF] Bussmann | Radiative signatures of the relativistic Kelvin-Helmholtz instability[END_REF], UPIC [START_REF] Decyk | Particle-in-cell algorithms for emerging computer architectures[END_REF], ORB5 [START_REF] Jocksch | A bucket sort algorithm for the particle-in-cell method on manycore architectures[END_REF], PI-CADOR [START_REF] Surmin | Co-design of a particle-in-cell plasma simulation code for Intel xeon phi: A first look at knights landing[END_REF], PICSAR [START_REF] Vincenti | An efficient and portable SIMD algorithm for charge/current deposition in particle-in-cell codes[END_REF]. The L6D space-filling curve we propose can be seen as a technique leading to a similar organization of the particles: particles are kept together in memory also in a block-like fashion, cf. Fig. 6; nevertheless we do not apply a reordering at each time step, and when we reorder, we perform this operation by cell and not by super-cell.

Contributions

The main contributions of this paper are focused on the efficiency of a 3d PIC code on a multicore architecture. We also present results in the 2d setting that underline the difficulties when passing from 2d to 3d. Since it is well-known that PIC codes are memory-bound, in this article we propose two directions to reach memory efficiency:

Our previous work in 2d [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] showed that space-filling curves can optimize the cache performance. We design a new space-filling curve, the L6D, and we demonstrate that we achieved this aim in 3d by comparison with classic curves. In addition, we illustrate its efficiency by showing snapshots of real simulations in a two-dimensional setting.

We use a loop optimization, namely the strip-mining, in order to be able to vectorize the code and reduce memory transfers in 2d and 3d. By studying the memory bandwidth of our code from 1 to 24 cores, we provide an argument for the necessity of the strip-mining strategy. Moreover, we report the different outcomes of this loop optimization when going from 2d to 3d and thus, how to adapt this approach to the 3d case in order to achieve performance.

Single Socket Optimizations

In this section we describe different strategies to gain performance in 3d simulations on a single socket of 24 cores. The present study stems from [1, Section IV] where we performed a similar analysis on a single core. Our new approach is motivated by the differences between the optimization strategies on single core and on shared-memory multiple cores, due to the different ratios between computational performance and memory bandwidth of the two configurations.

Note that in many cases, 3d simulations benefit from the same optimizations as 2d simulations [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF]. However, we discuss in the sequel an additional optimization that behaves differently in 2d and in 3d.

All the simulations in this paper were conducted on the Marconi supercomputer, on which we were granted the use of 64 nodes with 2 sockets each. Each socket is an Intel Xeon Platinum 8160 @ 2.1 GHz (Skylake), with 96 GB of RAM, 6 memory channels, and 24 cores. Our C code was compiled using Intel C Compiler 17.0.4, using the FFTW3 library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] for the Poisson solver. The update-positions and accumulate loops are vectorized thanks to OpenMP 4.0 #pragma omp simd.

The results shown in this section were all obtained with the test case presented in Table 1. In addition, we also simulated nonlinear Landau damping and two-stream instability test cases. Theoretical results which allow to verify the code are available in [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF]. Thus, we checked the numerical conservation of the total energy and the numerical evolution in time of the electric field. The first optimization we implemented is the loop fission. More precisely, the loop "Foreach particle" in Fig. 1 This transformation leads to the pseudo-code shown in Fig. 2, and speeds up the code by 8.9%. Starting from now, we call "update-velocities loop" the loop which contains both the interpolation and the velocity update.

Data Structure and Data

Layout for E and ρ

Introduction

The aim of this section is two-fold: (i) first, we compare the relative performance of the standard 3d data structure for the electric field E to the redundant one and (ii) second, we assess the impact of different ordering strategies on both the electric field E and the charge density ρ. In the sequel, we suppose that ρ is always stored using the redundant data structure that has been shown to In PIC codes memory accesses are a well-known major bottleneck. Each time that the code accesses a cell of E or ρ, a contiguous portion of the array is loaded into the cache: ideally, all the computations using these data cells should be performed while the information is still there, avoiding to reload them later from the main memory. In this respect, the contribution we describe in this section shows that using different memory layouts for the redundant data structure decreases the number of cache misses.

It should be noted that a periodic sorting of the particles needs to be applied to improve data locality, since particles are moving at each iteration. In this manner, two particles contiguous in memory are in the same grid cell and thus, they access the same E (resp. ρ) cell during the interpolation (resp. accumulation). Nevertheless, sorting at every iteration would be computationally expensive and therefore we have to find a memory layout of the cells such that the cache benefits from the sorting last as long as possible. More precisely, with the previous notations, we aim at constructing a mapping (i x , i y , i z ) → i cell such that, when a particle moves from a cell to another, the probability that its new cell-index i cell is close to the old one should be high.

We remark that the mapping of the row-major ordering in equation (3) has advantageous data locality when a particle moves along the z-axis: if i z increases by one, the new cell-index also increases by one (except for particles on the top face of the grid), becoming exactly the index accessed by the following particles in the particle array. However, when a particle moves along the other axes, this favorable behavior is lost: if i y (resp. i x ) increases by one, the cell-index changes by ncz (resp. ncy • ncz) which implies cache misses for the values of E and ρ.

Four different strategies for ordering the cells have been tested. They are listed below from the least to the most computational-intensive, in terms of the computation of the mapping (i x , i y , i z ) → i cell :

(i) Scan-order (or row-major order): the canonical C memory layout.

(ii) L6D-order, cf. Fig. 6, which is an extension of the L4D-order, cf. [15, Section 2.1.] and Fig. 5. In this work, we present new algorithms to efficiently convert from and to these orderings (both in 2d and in 3d).

(iii) Morton-order ( N-order) or Lebesgue-order (Z-order), cf. [START_REF] Morton | A computer oriented geodetic data base and a new technique in file sequencing[END_REF] and Fig. 3.

Algorithms to convert from and to this ordering can be found in [START_REF] Raman | Converting to and from dilated integers[END_REF].

(iv) Hilbert-order, cf. [START_REF] Hilbert | Über die stetige abbildung einer linie auf ein flächenstück[END_REF] and Fig. 4. Algorithms to convert from and to this ordering can be found in [25].

Bijection Algorithms

The aim of this section is to detail the formulas and their implementation needed for the application of space-filling curves in a PIC code. When using the cell index plus offset representation for the particles, the update-positions loop can be written as in Listing 1.

Explanations of the correctness of this loop, together with a way of removing the * delta_t / delta_{x,y,z} computations (lines 2-4) can be found in [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF].

Below, we only focus on the i_{x,y,z}_from and i_cell_from functions.

In 2d, the mentioned functions are given in Listing 2 and are defined by

(i x ; i y ) → i cell = ncx • SIZE • i y /SIZE + SIZE • i x + mod(i y , SIZE) i cell →    i x = mod(i cell , ncx • SIZE)/SIZE i y = mod(i cell , SIZE) + SIZE • i cell /(ncx • SIZE) . (4)
They can be set as macros or inline functions. In our code, the constant parameters are part of the macro arguments; they are here omitted to present a shorter code, as it will not affect performance. First, the tile size is given, as global parameter (line 1). Then, depending on the size of the grid, we have a constant which should be set at the beginning of the simulation (line 3): the number of cells per column gives the number of cells that are in a full column of width TILE_SIZE, 64 in Fig. 5.

The algorithm can be explained as follows: first, we have to know in which column the index is. This is computed as i_y / TILE_SIZE. The starting index of this column is thus num_cells_per_column * (i_y / TILE_SIZE). Then, in this column of size ncx * TILE_SIZE, we recognize the row-major curve, which means we have to add i_col * TILE_SIZE + j_col. On the x-axis, the index is just i_x, and on the y-axis, the index is mod(i_y, TILE_SIZE).

In 3d, the mentioned functions are given in equation ( 5) and Listing 3; again, they can be set as macros or inline functions. In our code, the constant parameters are part of the macro arguments; they are here omitted. First, the tile size is given, as global parameter (line 1). Then, depending on the size of the grid, we have constants which should be set at the beginning of the simulation (line 3). The number of cells per tower gives the number of cells that are in a full tower of area TILE_SIZE * TILE_SIZE, 256 in Fig. 6. The number of cells per wall gives the number of cells that are in a full row of towers. In Fig. 6, there are 4 towers per wall (there are 4 towers of size TILE_SIZE = 4 in a row of ncx = 16 cells).

The explanation of the algorithm is similar to the 2d case. We have to know in which wall parallel to the (Oxz) plane the index is, then inside this wall in which tower it is, then this is a column-major ordering inside this tower.

(i x ; i y ; i z ) → i cell = ncz • SIZE 2 • i x /SIZE + i y /SIZE • ncx/SIZE + i z • SIZE 2 + mod(i y , SIZE) • SIZE + mod(i x , SIZE) i cell →                      i x = mod(i cell , SIZE) + SIZE • mod(i cell , ncz • SIZE 2 • ncx/SIZE )/(ncz • SIZE 2 ) i y = SIZE • i cell / ncz • SIZE 2 • ncx/SIZE + mod(i cell , SIZE 2 )/SIZE i z = mod(i cell , ncz • SIZE 2 )/(SIZE 2 )
(5)

Results

We present in Tables 2 and3 the performance gains when using the spacefilling curves described in the previous section. Our preceding contribution [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF] only showed results on one core. On modern architectures, there are usually more cores than memory channels: it is thus not straightforward to extrapolate the one core results on the full multicore architecture; we therefore show here results on the full processor. Moreover, additional arrays to store the indices i x and i y were used in [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF]. We show now that additional gains can be obtained with efficient computations of the bijection functions. In the tables, "arrays" means that we use additional arrays to store the indices, otherwise we recompute them. We focus on three meaningful comparisons in Tables 2 and3. The first one is on the data structure: is it beneficial to use the redundant one for the E 1.

arrays? The only point where the code changes is in the update-velocities loop.

We see in the tables that in 2d, it is detrimental to use it if we stick to the row-major curve, but it is already beneficial with the canonical curve in 3d. We recall that we use here many particles per grid cell, and that when using only a few particles per cell, the redundant data structure is not a good choice, see Section 3.4 for a detailed comparison.

The second comparison concerns the data layout. Is it possible to obtain notable gains by changing the ordering of the grid cells? There are two places in the code where the changes might become significant: in the interpolation and in the accumulation. This time, we can answer positively. Yes, by taking another order than the canonical one, we can save time (thank to a reduction in the cache misses [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF]Table 2]). In 2d, the L4D and Morton curves seem to give similar and optimal timings, while in 3d, the L6D curve allows additional gains and seems to be the best one.

The last comparison is on the particle data structure needed for the noncanonical orderings. We can either store the indices in additional arrays (here, arrays of short int), or re-compute them at each time step. When storing them, it requires more memory for the particles, therefore we need more time in the update-positions loop and in the sorting step.

Last but not least, we see a surprising timing of the update-positions loop when using the Hilbert ordering without additional arrays. It is due to the fact that the computation (i x , i y , i z ) from i cell is expensive and not vectorized. This is thus the only curve for which using additional arrays is profitable. Nevertheless, even with additional arrays, this curve does not improve performances compared to the standard layout, and consequently has to be discarded.

Additional Remarks

We have to choose carefully TILE_SIZE depending of the cache sizes. In our tests, TILE_SIZE = 8 led to the best timings. It can be replaced with other values, as long as they are not too large for the cache.

It should be noted that choosing a value of TILE_SIZE that does not divide the grid sizes is possible: then, there will be a few allocated cells that correspond to physical positions outside the boundaries and that will never be accessed.

As a side note, we can remark that if grid sizes are powers of two and if the architecture represents integers with two's complement, we can save some computations on each modulo operation (lines 11-13 in Listing 1). For example for the modulo in the x-axis, we can use the variable int ncx_minus_one = ncx -1

and then compute mod(i_x, ncx) as i_x & ncx_minus_one. This is more efficient than (i_x + ncx) % ncx on one core. However, when using the full 24-cores architecture, this small optimization brings no significant gains. The same goes for all the modulo computations in the L4D and L6D bijection functions, whenever we compute modulo a power of two.

We also note that it is possible to use space-filling curves without using the redundant data structure for E. We do not show here the corresponding results, but the conclusion is that when using a high number of particles per cell (as is our case) the redundant data structure turns out to be the best approach. This is clearly not true when using a low (e.g. less than a hundred) number of particles per cell.

In the previous paper [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF], we showed that in 2d we could save up to 43% cache misses on E and ρ. This was confirmed by measurements. We here give an intuition of why this works, in the particular case of our Landau damping test case, which is roughly isotropic. We show in Figs. 7, 8, 9 and 10 snapshots of the same 2d simulation, one using the row-major curve and the other using the L4D curve. For each ordering, two snapshots are taken: a first snapshot just after a sorting, and another one 3 iterations later. For each ordering, a coloring is applied on the particles to indicate which particles are close in memory. We will here focus on the black particles in the center of the row-major ordering just after the sorting, and on the deep purple particles in the center of the L4D ordering just after the sorting. Just after the sorting, the black particles occupy only the center cells. Three iterations later, they are spread on 7 as many cells (some move up to 3 cells up, some move up to 3 cells down -the ones that move horizontally do not incur the use of new cells). Just after the sorting, the deep purple particles occupy only the center cells. Three iterations later, they are spread on 3 as many cells (they are spread on a 14 x 14 square instead of a 8 x 8 square). This explains why, when fetching the E values and when updating the ρ values, the L4D curve will lead to less cache misses than the row-major one. In 3d, the explanations are similar, but harder to show on a figure.

Strip-mining [21, Section 9.8]

Even though the loop fission gives satisfactory results, the resulting code still needs to scan some particle arrays three times, thus putting a lot of pressure on the memory bus. A natural loop transformation that comes after is thus the strip-mining. Instead of having three loops each scanning all the particles, we split the particle arrays in sub-arrays of size k (where k has to be chosen, depending on the architecture), and have the three loops operate only on k particles. Thus, for each particle, instead of having to fetch its properties in the main memory for each loop, it is now possible to fetch its properties in the cache for the two last loops. This transformation leads to the pseudo-code shown in Fig. 11, and speeds up the code by 22% in 2d.

We can note that if we choose k = 1, we are back to the base code, which was not optimal. If we choose k = N , we are back to the previous code with loop fission. In our 2d experiments, choosing a strip-size k between 64 and 256

gives similar optimal results.

Unfortunately, in 3d this strip-mining does not improve performances. This is explained by the fact that the cache is filled with too many E values, thus the expected gain in performance coming from the cache reuse of the particle arrays is out of reach. Thus, in 3d, to be able to efficiently reuse the particle data, the strip-mining has to be done only on the two last loops. This transformation leads to the pseudo-code shown in Fig. 12, and speeds up the code by 12%.

In our 3d experiments, choosing a strip-size k between 32 and 128 gives similar optimal results. This strip-mining technique was more or less already used in VPIC [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF], which 20 advanced 4 particles at a time for vectorization (k = 4 with our notations). However, in VPIC, an additional assumption was made -none of those 4 particles should cross cell boundaries -and scalar code was generated for the particles that crossed cell boundaries. In our code, the update-positions loop is vectorized without exceptions.

Overall gains and comparisons

The optimizations presented in this section are summarized in Table 4. In this table, the baseline is a version of the code with the standard 3d data structure for E, the redundant one for ρ, and a loop hoisting [START_REF] Severance | High Performance Computing[END_REF]Section 2.3.5.3] already applied that, in particular, remove the * delta_t / delta_{x,y,z} computations (lines 2-4 in Listing 1) [START_REF] Barsamian | Efficient data structures for a hybrid parallel and vectorized particle-in-cell code[END_REF]. The gains (in %) are computed with respect to the previous line of the table and the accumulated gains are computed with respect to the baseline.

Loop fission: better memory management for E and ρ, allows to vectorize the update-positions loop.

Redundant arrays for E together with appropriate use of space-filling curves: less cache misses in the interpolation and in the accumulation.

We note that the redundant arrays are only useful when using at least a hundred particles per grid cell, see Table 5.

Strip-mining: allows to reuse particle data between loops. Those performances are compared in Table 5 to another recent PIC code [START_REF] Jocksch | A bucket sort algorithm for the particle-in-cell method on manycore architectures[END_REF],

solving the same equations with the same precision (both codes use double precision and first order interpolations). We ran simulation with the parameters chosen in this other paper, which differ from our previous ones. Simulations presented in that paper were conducted on the Piz Daint supercomputer consisting of 8 Sandy-Bridge cores @2.6 GHz with a theoretical memory bandwidth of 51.2 GB/s (or 6.4 GB/s/core), and we recall that we used the Marconi supercomputer consisting of 24 Skylake cores @2.1 GHz with a theoretical memory bandwidth of 127.99 GB/s (or 5.3 GB/s/core). PIC codes being memory bound, comparing performances per core on those two architectures makes sense.

In those cases, there are only 7.6 particles per cell (when using 1 million particles), and 122 particles per cell (when using 16 million particles). The redundant data structure for E is a good choice only if there are a lot of particles per cell (in this paper, we use 3, 815 particles per cell for our test case in Table 1).

Of course, when there are roughly as many particles as grid points, multiplying by 8 the data for E almost doubles the memory transfers. With such a low number of particles, using a redundant data structure for E is detrimental, and we would have to make another study for such a configuration. But with 122 particles per cell, the data layouts and code transformations presented in this paper become useful. At the scale of 1 billion particles, we recall that our code needs only 35.47 ns per particle per iteration per core.

4. Thread-level and process-level parallelism

Process-level parallelism

The state-of-the-art approach for parallelizing PIC simulations on distributed memory machines is to decompose the physical domain into smaller sub-domains and to assign the particles inside a sub-domain to a processor (among the wide literature, see e.g. [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF][START_REF] Germaschewski | The plasma simulation code: A modern particle-in-cell code with patch-based load-balancing[END_REF]). In grid-based simulations, this method was established to give good scaling, as long as the work due to the communications through the sub-domain boundaries remains small compared to the computations inside the sub-domains. However, the main drawback of this technique is the difficulty of maintaining the load balance.

In this work, we handle the process-level parallelism with particle decomposition instead of domain decomposition: during the whole simulation, every process holds a fixed amount of particles but it keeps track of the whole grid quantities. Thus, at every iteration, every process accumulates the charge density associated with its particles and an MPI ALLREDUCE gives the total charge density. The Poisson equation is then solved by every process over the whole grid.

The main advantage of this method is its simplicity: the only communication is via MPI ALLREDUCE for the reduction of the charge array and no particle has to move from one process to another during the simulation. Thus, all the computations are automatically work-balanced.

The bottleneck of this approach usually cited in the literature is that the scalability is highly limited by the global reduction step. Of course, when using a large number of processes, the communication becomes too costly. But in practice, if we use the full memory on each MPI process to put particles, this is not the main drawback. The main bottleneck of this approach is that in realistic simulations (using a grid bigger than our 64 x 64 x 64 grid), there are so many cells that the computations on one process will be inefficient, due to the high number of cache misses involved. This behavior is quite common for example in matrix computations, where you can see super-scaling behaviors on large matrices, due to the large reduction of cache misses when computing only on sub-blocks. This is the reason why the technique we propose in this paper should be seen as only a small brick inside a more complex scheme: one should probably use domain decomposition on top of the efficient OpenMP algorithm we provided.

Strong scaling results on 24 threads with OpenMP

We next show a strong scaling up to 24 cores of two versions of our code.

Fig. 14 shows a strong scaling of the code with loop fission (pseudo-code given in Fig. 2) and Fig. 15 show a strong scaling of our code with strip-mining (pseudocode in Fig. 12).

Fig. 13 shows both codes on a single graph, for comparison. On 1 thread, the loop fission code is 6.7% faster, but on 24 threads the strip-mining code is 12% faster.

These figures illustrate the importance of having an efficient code on one core, but also the importance of precise performance analysis to enhance multicore efficiency. On the comparison graph, we see that up to 8 threads, the code with 3 loops performs better, then the code with strip-mining performs better. To understand why, we can look at the two other figures, where we show the memory bandwidth of our code, compared to that of the triad test in the Stream benchmark [START_REF] Mccalpin | Memory bandwidth and machine balance in current high performance computers[END_REF]. On one hand, these histograms underline that the update-velocities, accumulation and sorting steps are far to reach the peak memory bandwidth and thus, they have a good scaling up to 24 threads. On the other hand, the update-positions step reaches the same memory bandwidth as the Stream benchmark (the theoretical peak on 24 threads is 127.99 GB/s).

Accordingly, this step cannot be further accelerated when using 24 threads. The strip-mining idea thus becomes natural when looking at those figures: when the update-positions loop cannot be further sped up by more threads, due to memory bandwidth limits, merging it with another loop becomes a good idea.

Weak scaling on 3,072 cores with MPI

Our parallel results come from simulations executed on the supercomputer Marconi. Each node has 2 sockets of 24 cores each, hence we used one MPI process per socket and 24 threads per process. Fig. 16 shows a weak scaling from 1 core to 3,072 cores (64 nodes). These simulations run with 1 billion particles per socket in order to use the full memory.

We can see that up to 3,072 cores, the overhead due to MPI ALLREDUCE stays acceptable, thanks to the hybrid parallelism OpenMP + MPI.

Conclusion and Outlook

In order to simulate kinetic plasmas, we developed a three-dimensional Particle-in-Cell code and studied its performance. We thus explored several space-filling curves for the data layout of E and ρ, and different loop transformations for the particle loops, with the aim of improving the cache reuse and achieving efficient vectorization. More precisely, we introduced a novel spacefilling curve in 3d, called L6D, which have been demonstrated to be the best strategy in the case of the Landau damping simulations. We compared our results to those of [START_REF] Jocksch | A bucket sort algorithm for the particle-in-cell method on manycore architectures[END_REF] and obtained significant gains when using at least a hundred particles per grid cell. We showed that strip-mining can be used in addition to the vectorization technique from [START_REF] Vincenti | An efficient and portable SIMD algorithm for charge/current deposition in particle-in-cell codes[END_REF] to further improve the code efficiency. We also implemented a hybrid MPI/OpenMP parallelism and we addressed memory bandwidth issues for justifying the scaling results with OpenMP.

In the future, it would be interesting to port our PIC code to Many Integrated Core (MIC) architectures to reinforce the multi-threading. Another challenging question would be to apply these optimization techniques to electromagnetic codes, in which the data structure analysis should be extended when solving the Maxwell equations. The efficient PIC code we developed in this paper thus opens up the possibility to run realistic simulations in plasma physics in a three-dimensional physical space.

ParametersN3

  : number of particles. ncx × ncy × ncz: number of grid cells. ∆t: time step.f 0 : initial distribution function. q and m: particle charge and mass. Variables particles[N ]: set of particles, with position xp and velocity vp. ρ[ncx][ncy][ncz]: charge density. E[ncx][ncy][ncz]: electric field. Algorithm 1 Randomly initialize N particles following f 0 2 Compute initial ρ and E Foreach

  cally represented by three real numbers, the position is identified in the present work with a single cell index i cell and three normalized offsets within this cell. The advantages of this representation are exposed in [4, Section III-E]. In addition to the parameters explained in Fig. 1, we denote the physical space by [x min ; x max ) × [y min ; y max ) × [z min ; z max ), and the grid spacing by ∆x = (x max -x min )/ncx, ∆y = (y max -y min )/ncy and ∆z = (z max -z min )/ncz. Thus, a particle positioned at (x physical , y physical , z physical ) is mapped on the grid at the position (x, y, z) ∈ [0; ncx) × [0; ncy) × [0; ncz), where x = x physical -x min ∆x , y = y physical -y min ∆y and z = z physical -z min ∆z .

  is broken into three parts: one loop to interpolate E and to update velocities, one to update positions, and one loop to accumulate the charge. There are two main reasons to use three loops instead of one: (a) we can efficiently vectorize the update-positions as a stand-alone loop and (b) a separate processing of the arrays of E and ρ in different loops leads to a better overall memory management.

Figure 2 :

 2 Figure 2: Particle-in-Cell pseudo-code with loop fission.

Figure 3 :Figure 4 :

 34 Figure 3: Morton layout of a 16 x 16 matrix.

Figure 5 :Figure 6 :

 56 Figure 5: L4D layout of a 16 x 16 matrix, TILE SIZE=4.

Figure 7 :

 7 Figure 7: Row-major ordering, just after a sorting.

Figure 8 :

 8 Figure 8: Row-major ordering, 3 iterations after a sorting.

Figure 9 :

 9 Figure 9: L4D ordering, just after a sorting.Figure 10: L4D ordering, 3 iterations after a sorting.

Figure 11 :4Figure 12 :

 1112 Figure 11: Particle-in-Cell pseudo-code with strip-mining.

Figure 13 :

 13 Figure 13: Loop fission (3 loops) VS Strip-mining (2 loops).

Figure 16 :

 16 Figure 16: Weak scaling on Marconi.

Table 1 :

 1 Test case for OpenMP optimizations.

	Physical test case	Linear Landau damping [2, 5.15], initial distribution
		1 + 0.01 cos x 2 cos y 2 cos z 2	1 (2π) 3/2 exp -|v| 2 2
	Spatial grid	[0; 4π) 3 decomposed in 64 3 cells, periodic boundaries
	Number of particles	1 billion
	Number of iterations 100 (sorting every 10 iterations)
	Time step	0.05
	Particle crossing: av-	49.4% of the particles move 1 cell away,
	eraged, per iteration	0.00150% of the particles move 2 cells away
	Architecture	24 cores on Intel Skylake

3.1. Loop fission

[START_REF] Wolfe | High Performance Compilers for Parallel Computing[END_REF] Section 9.3]

Table 3

 3 Time spent in the different loops (in seconds). Test case in Table

			: 3d Space-filling Curves Timings
		Up. v Up. x Acc. Sort Total
	3d standard	126.7	55.3	31.5 21.5 235.8
	Row-major	92.6	55.3	31.5 21.4 201.7
	L6D arrays	92.8	79.0	30.4 29.5 232.6
	Morton arrays	96.5	79.0	30.3 27.5 234.2
	Hilbert arrays	95.3	80.4	31.1 26.9 234.8
	L6D	85.5	55.5	29.9 20.9 192.9
	Morton	89.4	56.7	33.5 19.8 200.3
	Hilbert	87.3	244.4 29.2 20.3 382.2

1

  Foreach subset of k particles in particles,

	2	Foreach particle in this subset,
	3	Interpolate E to particle
	4	Update the velocity
	5	Foreach particle in this subset,
	6	Update the position
	7	Foreach particle in this subset,
	8	Accumulate particle charge to ρ

Table 4 :

 4 Gains of Different Optimizations

	Time (s)	Gains	Accum. gains
	Baseline	258.7	0.0%	0.0%
	+ Loop Fission	235.8	8.9%	8.9%
	+ Space-filling curves (L6D)	192.9	18.2%	25.4%
	+ Strip-mining	169.5	12.1%	34.5%
	Total execution time, gains and accumulated gains. Test case in Table 1.
	Overall, these optimizations result in 590 million particles processed per sec-
	ond, on 24 cores on Intel Skylake architecture, without hyper-threading (25 mil-
	lion particles per second per core), or 1.48 ns per particle per time step (35.47 ns
	per particle per time step per core).

Table 5 :

 5 Comparison with Another Code

	X Code X X	X X	X	X X Nb. of particles X X X X X X X	10 6	16 • 10 6
	Jocksch et al.	work [9]	153.92	115.76
	(8 cores, 6.4 GB/s/core)
	This	paper	(24	cores,	854.33	93.07
	5.3 GB/s/core)
	Time spent per particle per iteration per core, in nanoseconds. Test case [9]:
	512 x 256 x 1 grid. Initial particle distribution uniform in space and velocity
	with |v max | = 1.

Table 1 .

 1 To merge two loops, we can either "undo" part of our loop fission, by applying loop fusion [21, Section 9.2]; or we can use strip-mining. To preserve the high efficiency of this vectorized loop, it is better to use strip-mining rather than loop fusion, which would break the vectorization opportunity.

				Marconi (Intel Skylake) -1 billion particles -64 x 64 x 64 Total execution time (s) with respect to the number of threads
	Total execution time (s)	500 1000 1500 2000 2500 3000 3500							Sort Update v Update x + Accumulate
		0	Loop fission	Strip-mined	Loop fission	Strip-mined	Loop fission	Strip-mined	Loop fission	Strip-mined	Loop fission	Strip-mined	Loop fission	Strip-mined
			1 thread	2 threads	4 threads	8 threads	16 threads	24 threads
						Test case in

Table 1 .

 1 Figure 14: Memory Bandwidth of the loop-fission code (3 loops).

	Memory Bandwidth (GB / s)	Stream Update v Update x Accumulate Sort x1,9	x2	x3,7	x3,9	x6,5 x8,1 x6,8 x7,9	x7,6 x15 x8,2 x15 x13 x7,6 x20 x8,2 x22 x16
					x3,9	x7.2
			x2	x2	x2	x4 x3.8
		1 thread	2 threads 4 threads 8 threads 16 threads 24 threads
					Test case in

Table 1 .

 1 Figure 15: Memory Bandwidth of the strip-mined code (2 loops).

					Number of cores
		24	48	96	192	384	768	1536	3072
	Execution time (seconds)	50 100 150 200	Total execution time Communication time
		1 0.04% 0	2 0.16%	4 0.21%	8 0.70%	16 1.26%	32 3.29%	64 4.54%	128 6.62%
					Number of sockets

Test case: 64 x 64 x 64 grid, 1 billion particles per socket, 100 iterations simulation (sorting every 10 iterations). Architecture: Skylake. Communication time is also shown as percentage of the execution time.

Acknowledgment This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom Research and Training Program 2014-2018 under Grant Agreement No. 633053. Simulations were run on the EUROfusion Marconi supercomputer, in the context of the Selavlas project led by K. Kormann. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Charlie Stephan

Moussa Dicko

Paolo Stringari

Christophe Coquelet
email: christophe.coquelet@mines-paristech.fr

Liquid-Liquid Equilibria of Water + Solutes (Acetic Acid/ Acetol/Furfural/Guaiacol/Methanol/Phenol/Propanal) + Solvents (Isopropyl acetate/Toluene) Ternary Systems for Pyrolysis Oil Fractionation

Keywords: pyrolysis oil, liquid-liquid equilibrium, solvents extraction, biorefinery

Liquid-liquid equilibria of water + solutes (acetic acid/ acetol/furfural/guaiacol/methanol/phenol/propanal) + solvents (isopropyl acetate/toluene) ternary systems for pyrolysis oil fractionation

Introduction

The use of renewable and low environmental impact raw materials for energy and chemicals production is one of the major 21st century's concerns. Bio-refineries are particularly adapted to this requirement since their aim is to produce multiple products including bio-sourced fuels and industrial chemicals [START_REF] Cherubini | The biorefinery concept: Using biomass instead of oil for producing energy and chemicals[END_REF]. Biomass can be converted in various chemicals via thermochemical treatments such as gasification, fast pyrolysis and liquefaction. Fast pyrolysis of lignocellulosic biomass leads to solid char, bio-oil and non-condensable gas with different yields depending on the type of biomass treated and on the operating conditions of pyrolysis process [START_REF] Bridgwater | An overview of fast pyrolysis of biomass[END_REF]. Pyrolysis oil, which is the main product of this process, contains hundreds of oxygenated compounds such as acids, ketones, aldehydes, alcohols, furanic and phenolic compounds [START_REF] Chiodo | Pyrolysis of different biomass: Direct comparison among Posidonia Oceanica, Lacustrine Alga and White-Pine[END_REF]- [START_REF] Montoya | Bio-oil production from Colombian bagasse by fast pyrolysis in a fluidized bed: An experimental study[END_REF]. Extraction of organic compounds from pyrolysis oil aqueous fraction is a promising way for the optimization of second generation bio-refineries, as it allows producing valuables chemicals, while the organic fraction can be used to produce bio-fuels [START_REF] Vitasari | Water extraction of pyrolysis oil: The first step for the recovery of renewable chemicals[END_REF], [START_REF] Vitasari | Glycolaldehyde co-extraction during the reactive extraction of acetic acid with tri-n-octylamine/2-ethyl-1-hexanol from a wood-based pyrolysis oil-derived aqueous phase[END_REF].

Concerning the process, liquid-liquid extraction is one of the most suitable operations for pyrolysis oil aqueous phase fractionation. By addition of a solvent in the pyrolysis oil aqueous phase, the organic compounds distribute in certain proportion in the water-rich phase and in the solvent-rich phase. The possibility of recovering the target chemicals depends on the selectivity of the solvent. The choice of the solvent is then one key aspect in the design of the liquid-liquid extraction operation. Many solvents have been investigated based on the methodology developed by Gani et al. [START_REF] Gani | Method for selection of solvents for promotion of organic reactions[END_REF]. The following criterions have been used: price, toxicity, physico-chemical properties (such as viscosity, interfacial tension, boiling temperature…) and efficiency (selectivity, partition coefficients, solvent loss). The price, toxicity and physico-chemical data were retrieved from literature and the efficiency criteria were predicted with the modified UNIFAC (Dortmund) model [START_REF] Gmehling | A Modified UNIFAC Model. 2. Present Parameter Matrix and Results for Different Thermodynamic Properties[END_REF]. The solvents selected were Toluene and Isopropyl acetate. Toluene is used to extract aldehydes, ketones, phenolic and furanic compounds from aqueous phase. Isopropyl acetate is used to extract the other components from aqueous phase (acids, alcohols, cetoalcohols, and glycols).

The process conception and design of separation operations requires accurate thermodynamic models.

In the case of liquid-liquid extraction, the model will give information about the Liquid-Liquid Equilibrium (LLE) between the solutes, the solvents and water (diluent) and will calculate the partition coefficient and the selectivity of the solutes between aqueous and organic phase.

The aim of this work is to provide an accurate thermodynamic model for the representation of liquidliquid equilibria between water, pyrolysis oil organic compounds and extraction solvents. The solutes considered are acetol, acetic acid, furfural, guaiacol, methanol, phenol and propanal, which are representative molecules of the different chemical families composing the pyrolysis oil aqueous phase. Some experimental LLE datasets have already been measured in this way and were available in the literature. This article also provides new sets of LLE data. Table 1 presents the LLE systems data measured in this work and the ones already available in the literature. All these LLE data were used for the modeling. Water was purified by a Millipore system with a 0.22µm membrane.

Experimental apparatus

Experimental work was carried out at CTP laboratory. The apparatus and experimental procedure for the LLE measurements were identical to the ones presented in a previous publication of Auger et al.

[21] and were based on the "static analytical method". The LLE glass cell, with a volume of 100mL, was enclosed with a stainless steel end-cap, and was equipped with three sampling points for the liquid phase I, the liquid phase II and the vapor phase (which was not used). The sampling points were equipped with air-tight septa to prevent any loss of material. The temperature was measured using a Pt100 Ω platinum probe ( ).

All measurements were done at atmospheric pressure. The chemicals were introduced in the cell in adequate amount so that the interface between the liquid phases was above the heavy liquid phase sampling point and under the light liquid phase sampling point. The cell was introduced in a water bath and the temperature set to the desired value and the magnetic stirrer switched on. The two phase mixture was stirred for at least 1 hour and then the stirrer was switched off and the mixture was allowed to reach equilibrium during at least 2 hours. Approximately 0.5mL of each phase was sampled and their compositions were measured by gas chromatography. All samples were injected into the Gas Chromatograph (temperature of injector is equal to 200°C, flow of carrier gas (helium) equal to 20mL per minute). FID temperature was set to 200°C and TCD temperature was set to 180°C. All columns used were 2 meters length and with a diameter of 2mm. Table 3 presents the columns names and the temperature profiles for the systems studied. The GC detectors (TCD and FID) were calibrated by injecting known quantities of each chemical using GC liquid syringes. Uncertainties methods and calculations are provided in appendix.

Othmer and Tobias correlation for tie-lines [START_REF] Othmer | Liquid-Liquid Extraction Data -The Line Correlation[END_REF] was used to assess the smoothed trend of the data measured.

Models

Thermodynamic models

The activity coefficients were calculated using NRTL [START_REF] Renon | Local compositions in thermodynamic excess functions for liquid mixtures[END_REF] and UNIQUAC [START_REF] Abrams | Statistical thermodynamics of liquid mixtures: A new expression for the excess Gibbs energy of partly or completely miscible systems[END_REF] Excess Gibbs Energy models.

For NRTL model, and , where is the non-randomness parameters, and and are respectively the temperature dependent and temperature independent binary interaction parameters. The use of NRTL requires three sets of parameters , and . For the studied systems was set equal to .

For UNIQUAC model, , where and are respectively the temperature dependent and temperature independent binary interaction parameters. The use of UNIQUAC requires two sets of parameters and . The use of UNIQUAC model requires two sets of parameters , , and and parameters that are listed in Table 4.

Data regression

The parameters of NRTL equation were set to the recommended values by Renon and Prausnitz [START_REF] Renon | Local compositions in thermodynamic excess functions for liquid mixtures[END_REF], depending on the type of components (non-polar, polar and self associative, or polar and non-self associative). When the correlation between the data calculated with the NRTL model and the experimental data was not satisfactory the parameters was regressed.

In order to regress the binary interaction parameters for both models, different objective functions were used depending of the available data. For the pairs of components with partial miscibility the binary interaction parameters were regressed on the binary systems LLE data, and for the pairs of miscible components the binary interaction parameters were regressed on the ternary systems LLE data.

Binary systems:

For both models binary interaction parameters between two components with partial miscibility, as it is the case in systems 1 to 6 (see Table 1), were regressed on the binary systems LLE. For these pairs of components both and parameters of each activity coefficient model were regressed as the data available for the corresponding binary systems cover a large range of temperature. The objective function used for the regression of these parameters was the following equation (Eq. 5).

(Eq. 5)

Where is the mole fraction calculated of component for the experimental value;

is the mole fraction measured of component for the experimental value; is the number of points in the data sets; and is the uncertainty on the mole fraction measured of component for the experimental value.

Ternary systems:

The binary interaction parameters of components that are miscible were regressed on ternary systems LLE data sets. For these systems, only temperature independent parameters ( ) were regressed as the data sets are available at one, two, or three temperatures only (see Table 1). It is to be noted that the binary interaction parameters of water-methanol, water-acetic acid and water-acetol systems are involved in two ternary systems, corresponding to the two solvents considered. Thus, these parameters were regressed using both systems datasets at the same time in order to obtain consistent values for the parameters. The objective function (Eq. 6) takes into account the partition coefficient in the calculations of the parameters as it is a key value for the liquid-liquid extraction simulation and conception.

(Eq. 6)

Where Commercial software Simulis™ from PROSIM, France, was used to correlate the experimental data with the UNIQUAC and NRTL models.

Results and Discussion

In the supplementary work, Tables 1 to 8 present the results with the uncertainties obtained for the ternary systems LLE.

The empirical tie-line Othmer-Tobias correlation (Eq. 8, [START_REF] Othmer | Liquid-Liquid Extraction Data -The Line Correlation[END_REF]) was used in order to check the smoothed trend of the experimental LLE data.

(Eq. 8)

Where is the mole fraction of water in the aqueous phase, is the mole fraction of solvent in the organic phase, and are adjusted on experimental data. Two examples are presented in Figures 1 an d 2 for the water-acetol-toluene and water-guaiacol-isopropyl acetate systems respectively. All the charts are available in supplementary work (Figures 1 to 8). We can observe that all charts have a smoothed trend and most of the determination coefficients are close to 1. The only datasets that have low determination coefficients are the ones of the water-guaiacol-isopropyl acetate system. If we take into account the accuracies of the Othmer-Tobias's approach of testing reliability of LLE data, we can see that the datasets reliability is satisfactory (Figure 2). A density inversion phenomenon is observed for systems containing furfural, guaiacol or phenol.

These solutes have higher densities than water, while the solvents used in this work have lower densities than water. The organic phase, mostly containing solvent and solute, has lower density than the aqueous phase when its composition is rich in solvent, and higher density than the aqueous phase when its composition is rich in solute. When organic phase density is equal to aqueous phase density, an emulsion is obtained but decantation is impossible as gravity plays insufficient role as the driving force for phase separation. This phenomenon has not been deeply investigated because it has no impact on the liquid-liquid extraction process of the pyrolysis oil since the solutes concentrations are too low for causing the described phenomenon. However, it has an impact on the measurements.

When organic phase and aqueous phase densities are similar, decantation phenomenon is much slower and may be incomplete when sampling. The density inversion phenomenon can be investigated by analyzing the density variations of each phase in order to find the compositions in which the density of the two phases is the same. For water-acetol-isopropyl acetate system (Figure 3, supplementary work, Table 16, Figure 17) the compositions calculated by the models showed a significant deviation from the experimental data when approaching the liquid-liquid critical point. For acetol lower concentration, the models were accurate. NRTL model was slightly more accurate than UNIQUAC model for this system. For waterguaiacol-isopropyl acetate and water-furfural-isopropyl acetate systems (supplementary work, Tables 17 and18, Figures 18,[START_REF] Tamura | Ternary and Quaternary Liquid-Liquid Equilibria for Fuel Additives of the Water + Methanol + Toluene and Water + Methanol + Toluene + Methyl tert-Butyl Ether or tert-Amyl Methyl Ether Systems at 298[END_REF]20 and 21), an accurate representation of the data with NRTL model was obtained, but UNIQUAC model was unable to represent accurately the experimental LLE data for these systems. For water-furfural-toluene system (supplementary work, Table 24, Figure 32), an accurate representation of the data with UNIQUAC model was obtained, and NRTL model showed significant deviation for high furfural concentrations. For ternary systems containing phenol (supplementary work, Table 20, Figures 25 and26), a significant deviation between calculated and experimental data was obtained for both models.

The average absolute deviation (AAD, Eq.9) for all the systems is presented in Tables 5 and6 (Eq. 9) The AAD on the molar compositions is 0.0089 with the UNIQUAC model and 0.0082 with the NRTL model. We can observe that both models represented the data with a good agreement and that the NRTL model was slightly better for the systems considered. These two activity coefficient models calculate data with a significant deviation with experimental data. Practically, thermodynamic models with parameter adjusted on experimental data should represent the LLE for the considered systems more accurately. This deviation should be contextualized in the frame of extraction process simulations. The use of this model with the regressed parameters still represents an improvement, and it allows estimating the compositions and partition coefficients with an error of less than 20%.

The binary interaction parameters for both models are presented in Table 7.

Figure 1 :Figure 2 :

 12 Figure 1: Data treatment of the ternary system water (1) -acetol (2) -toluene (3) using the Othmer-Tobias correlation (o): 298.15 K, and (•): 323.15 K

For the purpose ofFigure 3 :Figure 4 :

 34 Figure 3: Ternary diagram for water -acetol -isopropyl acetate system LLE at 323.15K and atmospheric pressure, (□): experimental data [13], (◊): experimental data [this work], (-): UNIQUAC model, (---): NRTL model

Figure 5 :

 5 Figure 5: Ternary diagram for water -phenol -isopropyl acetate system LLE at 323.15K and atmospheric pressure. (□): experimental data [13]. (o): experimental data [11]. (◊): experimental data [this work]. (---): NRTL model with recommended parameter [23], (-):NRTL model with regressed parameter.

  , and the parity charts of partition coefficients of solutes (experimental vs. calculated) are presented in supplementary work (Figures 40 to 63).

Table 1 : Systems considered for the modeling of LLE

 1

	N°	Systems	Temperature (K)	Source
	1	Water + Furfural	273 -363	[9]
	2	Water + Guaiacol	298 -323	[10]
	3	Water + Phenol	273 -340	[11]
	4	Water + Propanal	288 -313	[12]
	5	Water + Isopropyl acetate	282 -348	[13]
	6	Water + Toluene	273 -363	[14]
	7	Water + Isopropyl acetate + Acetic acid	298	[15]
	8	Water + Isopropyl acetate + Acetol	298, 323	this work
	9	Water + Isopropyl acetate + Furfural	298, 323	this work
	10	Water + Isopropyl acetate + Guaiacol	298, 323	this work
	11	Water + Isopropyl acetate + Methanol	298, 323	this work
	12	Water + Isopropyl acetate + Phenol	298, 323	[16], this work
	13	Water + Isopropyl acetate + Propanal	298, 313	this work
	14	Water + Toluene + Acetic acid	298, 313	[17]
	15	Water + Toluene + Acetol	298, 323	this work
	16	Water + Toluene + Furfural	303	[18]
	17	Water + Toluene + Methanol	298	[19]
	18	Water + Toluene + Phenol	293, 303, 313	[20]
	19	Water + Toluene + Propanal	298, 313	this work
	2 Experimental
	2.1 Materials

Table 2

 2

	presents the chemicals used in this work. Components purity was assessed by gas
	chromatography and by measuring the refraction index of pure components, and comparing it with the
	reference values. No further purification was done after this analysis showed no significant impurities.

Table 2 : Chemical samples used for experimental work

 2

						Refractive index
	Chemicals	CAS number	Supplier	Purity (area %)	Analysis method a	Measured b at 20°C and	Supplier
						101.25 kPa
	Acetol	116-09-6	Fisher Scientific	> 95.0%	GC	1.4262	-
	Furfural	98-01-1	Sigma Aldrich	> 99.0%	GC	1.5264	1.5261
	Guaiacol	90-05-1	Fisher Scientific	> 99,0%	GC	1.5441	1.5429
	Methanol	67-56-1	Sigma Aldrich	> 99.99%	GC	1.3293	1.3292
	Phenol	108-95-2	Sigma Aldrich	> 99.5%	GC	1.5405 c	1.5408 c
	Propanal	123-38-6	Fisher Scientific	> 99.0%	GC	1.3644	1.3646
	Toluene	108-88-3	Sigma Aldrich	> 99.9%	GC	1.4969	1.4967
	Isopropyl acetate	108-21-4	VWR	> 99.0%	GC	1.3779	1.3773
	a GC: gas chromatography. b Measured with an Anton Paar Abbemat 300 Refractometer (standard uncertainty: refractive index u(nD) = 10-4
	pressure (GE Druck DPI 142 Barometric Indicator) u(p) = 0.029 kPa, and temperature u(T) = 0.03 °C). c Measured at T = 41°C

Table 3 : Chromatograph conditions for the compositions measurements within the different systems

 3

				Temperature profile
	Systems	Column name	Initial temperature (K)	Level time (min)	Increasing rate (°C/min)	Final temperature (K)	Level time (min)
	8, 9, 10, 12	Chromosorb100	90	2	5	110	2
	11, 13	20%OV210	120	2	20	180	5
	15	30%RTX1	85	6	-	-	-
	19	20%OV210	60	2	5	80	-

Table 4 : ri and qi parameters of considered components for UNIQUAC model

 4

	Chemical	r i *	q i *
	Acetol	2.874	2.612
	Acetic acid	2.202	2.072
	Furfural	3.172	2.500
	Guaiacol	4.531	3.512
	Methanol	1.431	1.432
	Phenol	3.552	2.680
	Propanal	2.574	2.336
	Water	0.920	1.400
	Toluene	3.923	2.968
	Isopropyl acetate	4.152	3.652

*Values from Simulis TM software

Table 5 : Average Absolute Deviations (AAD) for the compositions calculated with UNIQUAC and NRTL models with respect to the experimental data for binary systems

 5

				UNIQUAC				NRTL
	Systems*	Organic phase	Aqueous phase		Organic phase	Aqueous phase
		water	solute	water	solute	water	solute	water	solute
	1	0.0055	0.0055	0.00094	0.00094	0.0090	0.0091	0.0013	0.0013
	2	0.0024	0.0024	0.00010	0.00010	0.0038	0.0038	0.000031	0.000031
	3	0.0112	0.0145	0.0059	0.0050	0.0078	0.0097	0.0052	0.0038
	4	0.0139	0.0129	0.0096	0.0089	0.0071	0.0071	0.0078	0.0078
	5	0.0034	0.0026	0.00011	0.00011	0.0024	0.0028	0.00014	0.00015
	6	0.00029	0.00030	0.000023	0.000024	0.00026	0.00026	0.000023	0.000022
	AAD	0.0061	0.0064	0.0028	0.0025	0.0051	0.0055	0.0024	0.0022
	*The system numbers refers to the ones in Table 1

Table 6 : Average Absolute Deviations (AAD) for the compositions calculated with UNIQUAC and NRTL models with respect to the experimental data for ternary systems

 6

	UNIQUAC
	Systems*

*The system numbers refers to the ones in Table

1

Table 7 : Binary interaction parameters for UNIQUAC and NRTL models

 7

	Component i	Component j		UNIQUAC				NRTL
			aij (-)	aji (-)	bij (K)	bji (K)	aij (-)	aji (-)	bij (K)	bji (K)	α (-)
	water	furfural	-1.782	3.109	425.6	-1111	3.471	-3.338	31.44	1264	0.3
	water	guaiacol	-3.988	5.600	1096	-1857	5.884	-5.185	-288.8	1807	0.3
	water	isopropyl acetate	-1.828	1.774	450.7	-1015	7.790	-0.4584	-1015	669.9	0.3
	water	phenol	-0.01308	1.062	-569.5	-55.63	1.228	-5.239	764.7	1662	0.47
	water	propanal	-2.429	1.215	710.5	-710.0	7.604	1.696	-1674	-217.4	0.3
	water	toluene	-1.102	1.766	29.73	-1373	7.119	-2.440	-96.85	1855	0.2
	water	acetol	0	0	446.4	-1997	0	0	-304.7	388.3	0.3
	water	acetic acid	0	0	-20.50	222.3	0	0	-409.4	904.7	0.3
	water	methanol	0	0	497.3	-2023	0	0	-134.6	-674.4	-0.093
	isopropyl acetate	propanal	0	0	17.72	-20.29	0	0	3844	-181.4	0.3
	isopropyl acetate	furfural	0	0	-10.79	10.49	0	0	-99.63	205.7	0.47
	isopropyl acetate	phenol	0	0	-547.8	384.6	0	0	825.0	-275.3	0.3
	isopropyl acetate	guaiacol	0	0	-4.119	45.83	0	0	616.4	-182.2	0.47
	isopropyl acetate	acetol	0	0	91.54	-137.8	0	0	493.4	296.2	0.3
	isopropyl acetate	acetic acid	0	0	283.8	-29.20	0	0	-79.78	14.74	0.3
	isopropyl acetate	methanol	0	0	-445.5	164.4	0	0	614.0	-435.4	0.3
	toluene	propanal	0	0	-23.70	3.452	0	0	-98.03	-276.3	0.3
	toluene	furfural	0	0	-1290	272.1	0	0	-166.1	498.9	0.3
	toluene	phenol	0	0	-291.1	91.56	0	0	-419.8	314.3	0.3
	toluene	guaiacol	0	0	-0.1296	-311.0	0	0	-445.1	265.4	0.3
	toluene	acetol	0	0	-109.5	-74.79	0	0	393.2	439.9	0.47
	toluene	acetic acid	0	0	28.90	-124.7	0	0	73.64	369.3	0.47
	toluene	methanol	0	0	-508.9	43.68	0	0	399.3	384.4	0.47
	5 Conclusion

In this study, new LLE data were measured from 298.15 K to 323.15 K at atmospheric pressure for eight ternary mixtures with static analytical method. These mixtures represent typical systems involved in biorefinery processes. They were combined with other LLE data available in the literature

and were represented accurately with NRTL and UNIQUAC activity coefficient models. The modeling of these systems is valuable for the simulation of liquid-liquid extraction process of oxygenated compounds from pyrolysis oil aqueous phase using two solvents. According to the literature and the measured data, toluene is considered to be a good candidate as extraction solvent for phenol, furfural and propanal. However, toluene is a toxic a petroleum based solvent and alternative solvents should be investigated. Similarly, Isopropyl acetate is considered to be a good candidate as extraction solvent for all the solutes considered.

Appendix

The uncertainties on the molar compositions have been calculated based on the error of measurements on the FID and TCD. The calibration of the detectors was done by injecting known values of products.

Linear correlations between the chromatogram peaks areas and the molar quantities of injected products was done. The relative error of the correlations is less than 6.0%. The distribution of the values of the molar quantities has been considered rectangular. Thus, the uncertainty on the molar quantities is .

The calculation of the uncertainties on the molar composition was calculated as follow:


Deresh Ramjugernath

Alain Valtz

Dominique Richon

Mark D Williams-Wynn

Christophe Coquelet

Isothermal Vapor-Liquid Equilibrium Data

Introduction

There has been little reported in open literature on studies related to the use of fluorinated refrigerants as solvents for hydrocarbon systems. A partial list of authors who have reported the vapour-liquid equilibrium (VLE) data of hexafluoroethane (R-116) with n-alkanes was supplied by Ramjugernath et al. [START_REF] Ramjugernath | Isothermal Vapour-Liquid Equilibrium Data for the Hexafluoroethane (R116) + n-Butane System at Temperatures from[END_REF] Authors who have reported VLE measurements of hexafluoroethane and alkanes include Ramjugernath et al. [START_REF] Ramjugernath | Isothermal Vapour-Liquid Equilibrium Data for the Hexafluoroethane (R116) + n-Butane System at Temperatures from[END_REF][START_REF] Ramjugernath | Isothermal Vapor-Liquid Equilibrium Data for the Hexafluoroethane (R116) + Propane System at Temperatures from[END_REF] and Zhang et al. [START_REF] Zhang | Vapour-Liquid Equilibrium Measurements for Ethane + Hexafluoroethane System over a Temperature Range from[END_REF][START_REF] Zhang | Vapour-Liquid Equilibrium Measurements and Correlations for an Azeotropic System of Ethane + Hexafluoroethane[END_REF] The system of R-116 and ethane was investigated by Zhang et al., whereas VLE data sets for the system of R-116 with both propane and n-butane were investigated by Ramjugernath et al.

The VLE measurements for the R-116 and n-pentane system, as well as R-116 and n-hexane system reported herein are a continuation of a previous study into the use of R-116 as a solvent for the separation of short-chained alkanes, [START_REF] Ramjugernath | Isothermal Vapour-Liquid Equilibrium Data for the Hexafluoroethane (R116) + n-Butane System at Temperatures from[END_REF][START_REF] Ramjugernath | Isothermal Vapor-Liquid Equilibrium Data for the Hexafluoroethane (R116) + Propane System at Temperatures from[END_REF] and have not been previously investigated. The separation of alkanes is of importance for the purification of fuels from natural gas reserves.

These systems would likely have their phase diagrams classified as either type I or type II by the classification of van Konynenburg and Scott. [START_REF] Van Konynenburg | Critical Lines and Phase Equilibria in Binary Van Der Waals Mixtures[END_REF] In this classification scheme, a type I system is the simplest type, with only one critical locus curve between the two pure component critical points. This does not have a three-phase region. The type II system is more complex, and consists of two critical locus curves. One, between the two pure component critical points, and the other, starting at a fictitious 'critical point at infinite pressure', and ending at an upper critical end point (UCEP). This type of system has a region in which there is liquid-liquid immiscibility, and the UCEP is the upper temperature limit of this three-phase region. Low temperature VLE data is often required to distinguish between type I and type II systems. Each isothermal experimental data set reported herein was correlated with the Peng-Robinson Equation of State (PR-EOS). [START_REF] Peng | A New Two-Constant Equation of State[END_REF] The PR-EOS was modified by replacing the original alpha function with the Mathias-Copeman (M-C) alpha function. [START_REF] Mathias | Extension of the Peng-Robinson Equation of State to Complex Mixtures: Evaluation of the Various Forms of the Local Composition Concept[END_REF] This modification allowed a better representation of the vapour pressures of the n-pentane and the n-hexane due to the increase in the number of adjustable parameters that became available. The Wong-Sandler (W-S) mixing rules [START_REF] Wong | Theoretically Correct Mixing Rule for Cubic Equations of State[END_REF] were used in this correlation to represent the mixture properties with greater accuracy. These mixing rules utilised the NRTL excess Gibbs energy equation 9 to account for the liquid phase non-idealities. These measurements will be useful in determining the behaviour of the solutes in supercritical R-116. The interaction parameter values for group contribution methods could also be updated by using the data sets that are reported here.

Experimental

Materials. The suppliers, the CAS numbers, the stated purities and the measured purities for the chemicals used in these measurements are given in Table 1. The n-pentane and n-hexane were degassed in the equilibrium cell, upon loading, but the R-116 was not further purified.

The 'measured purities' of the pure components were estimated from the (gas chromatograph) GC peak area fractions, which were defined as the area under the peak for the given component, divided by the sum of all the peak areas. A VARIAN CP-3800 gas chromatograph fitted with a Porapack Q column with a mesh size of 80/100 was used. A thermal conductivity detector (TCD) was used to analyse these peaks.

The critical properties and acentric factors of R116, n-pentane and n-hexane are given in Table 2.

Experimental Apparatus. The "static-analytic type" apparatus, used for the VLE measurements of R-116 + propane and R-116 + n-butane by Ramjugernath et al., [START_REF] Ramjugernath | Isothermal Vapour-Liquid Equilibrium Data for the Hexafluoroethane (R116) + n-Butane System at Temperatures from[END_REF][START_REF] Ramjugernath | Isothermal Vapor-Liquid Equilibrium Data for the Hexafluoroethane (R116) + Propane System at Temperatures from[END_REF] was again used in this work. This apparatus was commissioned in the CTP laboratory (formerly CEP/TEP laboratory) at MINES ParisTech by Valtz et al., 10,11 based on the original static-analytic design of Laugier and Richon. 12 The uncertainties for the measured data were estimated by the guidelines from NIST, outlined by Taylor and Kuyatt. 13 All of the combined expanded uncertainties were calculated with a coverage factor of k = 2. This value corresponds to a 95 % confidence interval.

The two Pt-100 temperature probes that were used to monitor the temperature of the cell were calibrated according to the 1990 International Temperature Scale Protocol. The reference probes that were used for this calibration were supplied by TINSLEY Precision Instruments and were certified by the Laboratoire National d'Essais (Paris). The combined expanded uncertainty in the measurement of temperature was estimated to be U(T) = 0.03 K. The constant temperature bath, inside of which the VLE cell was submerged, had a thermal stability of 0.01 K.

The pressure transducer employed for the monitoring of the pressure within the cell was a PTX 611 pressure transducer. This transducer had a range of 0 to 20 MPa and was supplied by Druck. The calibration of this transducer was performed against a 5202S dead weight pressure balance which was supplied by Desgranges and Huot in France. The expanded uncertainty in the pressure readings from this transducer was estimated to be U(P) = 0.0006 MPa.

Composition analysis of the vapour and the liquid phases was performed by means of a VARIAN CP-3800 gas chromatograph. The signal from the thermal conductivity detector was integrated with respect to time, to obtain the peak areas for the different components. A 4-m x 1/8'' Porapak Q column with a mesh size of 80/100 (supplied by Restek, France) was used to separate the components from each of the samples. The gas chromatograph calibration was performed by injecting many differing, but known volumes of each pure component into the column with a gas syringe. The relationship between the peak area and the amount of component that was injected to produce this peak could therefore be determined and used as the calibration of the detector. Upon taking all the possible uncertainty sources into account when performing the calibration, the combined expanded uncertainty in the measurement of the composition was estimated to be less than or equal to 0.005. However, as has been shown by Lasala et al., 14 for any data regression that considers the compositions in the fitting procedure, the combined expanded uncertainties for each composition must be estimated individually.

Experimental Procedure. The equilibrium cell and the associated loading lines were all subjected to a vacuum (0.1 Pa, absolute, at room temperature), prior to the measurement of each isotherm. The cell was thereafter loaded with roughly 5 cm 3 of the liquid component (npentane or n-hexane) and the cell was again subjected to a vacuum; for the degassing of the liquid. The pressure was thereafter allowed to equilibrate at the temperature of the particular isotherm being measured. This pressure was recorded as the vapour pressure of the pure component. It was assumed that the equilibrium vapour pressure had been attained when the equilibrium cell temperatures, reported by both Pt-100 temperature probes, were stable and equal for at least ten minutes. The standard uncertainties were considered during this operation.

The R-116 was thereafter added in step-wise amounts, with the equilibrium pressure and phase compositions being measured upon reaching equilibrium after each successive cell loading. Equilibrium conditions were assumed to have been attained as soon as the changes in the cell pressure remained less than 0.1 kPa for a period of ten minutes. The contents of the cell were stirred vigorously throughout.

To ensure repeatability of the measurements, as well as to verify that the samples were representative of the equilibrium phase compositions, each phase was sampled at least five times. ROLSI® pneumatic samplers were used to perform the sampling of both the liquid and the vapour phases.

Correlations

The isothermal VLE data sets reported in this paper were correlated with a PR-EOS based model. [START_REF] Peng | A New Two-Constant Equation of State[END_REF] The PR-EOS, containing the M-C alpha function was fitted to the pure component vapour pressures to obtain the alpha function parameters.

A vapour pressure data set for R-116 was reported by Valtz et al. 15 For the n-pentane, vapour pressure measurements of Rice and El-Nikheli 16 and of Duce et al. 17 were used. For n-hexane, the vapour pressures reported by Rice and El-Nikheli, 16 Duce et al. 17 , Marrufo et al. [START_REF] Marrufo | Isobaric vapor-liquid equilibrium for binary mixtures of 1-hexene + n-hexane and cyclohexane + cyclohexene at 30, 60 and 101.3 kPa[END_REF] and Fernández et al. [START_REF] Fernández | Measurements of the Excess Properties and Vapor-Liquid Equilibria at 101.32 kPa for Mixtures of Ethyl Ethanoate + Alkanes (from C5 to C10)[END_REF] were used. The Mathias-Copeman alpha function coefficients regressed from the reported vapour pressures of the three components are listed in Table 3.

To determine the properties of the mixtures from pure component properties, the W-S mixing rule [START_REF] Wong | Theoretically Correct Mixing Rule for Cubic Equations of State[END_REF] with the NRTL local composition model was used. 9 A value of 0.3 was recommended for the alpha value in the NRTL model, for systems such as those investigated in this study. [START_REF] Fernández | Measurements of the Excess Properties and Vapor-Liquid Equilibria at 101.32 kPa for Mixtures of Ethyl Ethanoate + Alkanes (from C5 to C10)[END_REF] Therefore, the only parameters that were adjusted to fit the model to the binary experimental data were τij and τji, from the NRTL local composition model, and the binary interaction parameter, kij, from the W-S mixing rule. These three parameters were adjusted to obtain a good fit of the model to the experimental data by using a modified Simplex algorithm. The objective function that was used for this algorithm was a flash calculation function. This function is given by Equation 1.

F obj = 100 N [∑ ( x exp -x calc x exp ) 2 + ∑ ( y exp -y cal y exp ) 2 N 1 N 1 ] (1)
where N is the number of data points in the set; xexp is the experimental liquid mole fraction, xcal is the calculated liquid mole fraction, yexp is the experimental vapour mole fraction and ycal is the calculated vapour mole fraction.

The correlation of the experimental VLE data was performed using an 'in-house' data fitting software tool, developed by Armines and Mines ParisTech.

Results and Discussion

The measured vapour pressures for n-hexane [16][17][START_REF] Marrufo | Isobaric vapor-liquid equilibrium for binary mixtures of 1-hexene + n-hexane and cyclohexane + cyclohexene at 30, 60 and 101.3 kPa[END_REF][START_REF] Fernández | Measurements of the Excess Properties and Vapor-Liquid Equilibria at 101.32 kPa for Mixtures of Ethyl Ethanoate + Alkanes (from C5 to C10)[END_REF] were plotted against predicted vapour pressures in Figure 1. The vapour pressures calculated using the PR-EOS with the original alpha function and with the M-C alpha function were both predicted and plotted in Figure 1.

This comparison indicated that the use of the M-C alpha function provided a far better fit of the experimental data to the predicted data than when using the original PR-EOS alpha function. The PR-EOS with the M-C modification was therefore used for the VLE regression.

The VLE data set for the system of R-116 and n-pentane is reported in Table 4, with Figure 2 showing a plot of the experimental vapour-liquid equilibrium data alongside the data calculated from the regressed parameters. There was a good fit of the model to the measured data for this system. This system exhibits zeotropic behaviour. There are also no observable liquid immiscibilities, and thus, either the mixture is a type I system on the van Konynenburg and Scott 5 classification, or the mixture is a type II system that is above the UCEP.

The VLE data set for the system of R-116 and n-hexane is reported in Table 5. Also included in Table 5 are the vapour-liquid-liquid equilibrium (VLLE) and liquid-liquid equilibrium (LLE) data. The VLLE and LLE data points occurred at greater concentrations of R-116. These greater concentrations corresponded to the higher equilibrium pressures. A plot of the experimental VLE data against the data predicted with the model, using the regressed parameters, is given by Figure 3. This system exhibits zeotropic behaviour, with liquid miscibilities occurring at both measured temperatures. This system is therefore most likely to be a type II system. The VLLE and LLE data for the system were not included in the regression.

The regressed parameters that gave the best fit of the PR-EOS with the M-C alpha function and the W-S mixing rules, using the NRTL local composition model, to the experimental VLE data for both systems are given in Table 6.

Two summary statistics were used to indicate the degree of correlation between the experimental data and the models. These were the absolute mean relative deviation (MRD) of the regressed compositions and the bias of the regressed compositions. Equation 2gives the calculation of the absolute mean relative deviation and the calculation of the bias is given by Equation 3.

MRD (u)= 100 N ∑ | (u exp -u cal ) u exp | N i (2) BIAS (u)= 100 N ∑ ( (u 𝑒𝑥𝑝 -u 𝑐𝑎𝑙 ) u exp ) N i (3)
In both Equations 2 and 3, u is either the liquid or the vapour composition, and N is the number of experimental data points. 'exp' and 'cal' denote the experimental and calculated values respectively. The calculated values for the MRD and the bias are given in Table 7.

Conclusions

The vapour-liquid equilibrium measurements for the R-116 and n-pentane system and the R-116 and n-hexane system are reported herein. These measurements were performed using a "static-analytic" type VLE apparatus. Two isotherms were measured for each system, with one isotherm for each system being below the critical temperature of R-116 and the other above it. The measured VLE data was thereafter correlated with the Peng-Robinson EOS coupled with the Mathias-Copeman alpha function and the Wong-Sandler mixing rules. The Wong-Sandler mixing rules used the NRTL local composition model to describe the excess Gibbs energy of the system. Both systems exhibited zeotropic behaviour, with the R-116 + n-pentane system not exhibiting any liquid immiscibilities at the measured temperatures. The R-116 + nhexane system, however, presented a large miscibility gap at the temperatures measured. This miscibility gap was not modelled with thermodynamic models, due to the limited data that was available.

*Corresponding author. E-mail: williams-wynn@ukzn.ac.za.

Figure 2 . 1 Figure 3 .

 213 Figure 2. P-x-y data for the R116 (1) + n-pentane (2) system. ▲, T = 288.25 K; •, T = 296.24

Table 1 .

 1 Chemical Suppliers, CAS numbers, Stated Purities and Measured Purities Further purification of these components was performed by degassing in situ. Estimated from GC area fractions.

	Tel. +27 31 260 3128. Fax:

b c

Table 2 .

 2 Critical Properties (Tc, Pc) and Acentric Factors (ω) a The data values reported in this table were taken from Ref. 20.

		Tc / K	Pc / MPa	ω
	R-116	293.04	3.042	0.229
	n-pentane	469.70	3.369	0.251
	n-hexane	507.40	3.014	0.296

a

Table 3 .

 3 Mathias-Copeman Alpha Function coefficients, ci.

		c1	c2	c3
	R-116 a	0.79629	-0.55124	0.10316
	n-pentane b	0.73224	-0.09516	0.76161
	n-hexane c	0.78670	0.01877	0.50091
	a Obtained from regression of data found in reference 15
	b Obtained from regression of data found in references 16 and 17
	c Obtained from regression of data found in references 16-18

Table 4 .

 4 Vapour-liquid equilibrium (P-x-y) data for the R-116 (1) + n-Pentane (2) system at temperatures, T, of 288.25 K and 296.24 K. a

			T = 288.25 K					T = 296.24 K
	P/MPa	x1	U(x1) b	y1	U(y1) c	P/MPa	x1	U(x1) b	y1	U(y1) c
	0.048	0.000	-	0.000	-	0.066	0.000	-	0.000	-
	0.307	0.021	0.001	0.832	0.003	0.310	0.019	0.001	0.764	0.004
	0.699	0.060	0.001	0.920	0.001	0.705	0.053	0.001	0.896	0.002
	1.105	0.110	0.002	0.946	0.001	1.104	0.096	0.002	0.917	0.002
	1.497	0.183	0.003	0.957	0.001	1.508	0.153	0.003	0.936	0.001
	1.903	0.364	0.005	0.962	0.001	1.905	0.243	0.004	0.944	0.001
	2.152	0.754	0.004	0.965	0.001	2.300	0.477	0.005	0.947	0.001
	2.378	0.896	0.002	0.975	0.001	2.596	0.796	0.003	0.951	0.001
	2.5178 0.952	0.001	0.985	0.001	2.900	0.916	0.002	0.965	0.001
						2.994	0.940	0.001	0.970	0.001
	a Combined expanded uncertainties are U(P) = 0.0006 MPa and U(T) = 0.03 K (with k = 2)
	b Combined expanded uncertainty in liquid phase composition (k = 2)
	c Combined expanded uncertainty in vapour phase composition (k = 2)

Table 6 .

 6 Regressed temperature independent model parameters (𝜏 𝑖𝑗 and kij) for the PR-EOS with the M-C alpha function and the W-S mixing rule using the NRTL local composition model at temperatures, T.

		R-116 (1) + n-pentane (2)	R-116 (1) + n-hexane (2)
	T / K	288.25	296.24	288.24	296.23
	𝜏 12 /𝐽 • 𝑚𝑜𝑙 -1	1323.2	1988.2	4016.5	3583.2
	𝜏 21 /𝐽 • 𝑚𝑜𝑙 -1	3621.8	3213.0	2749.7	2692.9
	𝑘 12	0.356	0.319	0.305	0.332

Table 7 .

 7 Mean residual deviations and bias of compositions obtained in the fitting of the experimental VLE data to the PR-EOS with the M-C alpha function and the W-S mixing rule

	using the NRTL Model
		R-116 (1) + n-pentane (2)	R-116 (1) + n-hexane (2)
	T / K	288.25	296.24	288.24	296.23
	MRD x %	1.3	1.3	1.3	0.7
	bias x %	-0.8	-0.8	-0.8	-0.2
	MRD y %	0.2	0.2	0.2	0.5
	bias y %	0.1	0.1	0.1	-0.4

Figure 1. Vapour pressure of n-hexane. ○, Ref. 16; ∆, Ref. 17; +, Ref. 18; □, Ref. 19; ---, PR-EOS model with the original alpha function; ──, PR-EOS model with the M-C alpha function.

Acknowledgement

The authors would like to thank the NRF International Science Liaison Programme. Prof. D.

Ramjugernath would also like to thank MINES ParisTech for the use of the equipment in the CTP laboratories to enable these measurements to be undertaken.

Funding Sources

This work is based upon research supported by the South African Research Chairs Initiative of the Department of Science and Technology and National Research Foundation.

Author Information

Corresponding Author Table 5. Vapour-liquid equilibrium (P-x-y) data for the R-116 (1) + n-hexane (2) system at temperatures, T, of 288.24 K and 296.23 K. a e Supposed liquid-liquid equilibrium.