Pedro Valero-Lara
email: pedro.valero.lara@gmail.com

Francisco D Igual

Manuel Prieto-Matías

Alfredo Pinelli

Julien Favier

Accelerating fluid-solid simulations (Lattice-Boltzmann & Immersed-Boundary) on heterogeneous architectures

Keywords:

We propose a numerical approach based on the Lattice-Boltzmann (LBM) and Immersed Boundary (IB) methods to tackle the problem of the interaction of solids with an incompressible fluid flow, and its implementation on heterogeneous platforms based on data-parallel accelerators such as NVIDIA GPUs and the Intel Xeon Phi. We explain in detail the parallelization of these methods and describe a number of optimizations, mainly focusing on improving memory management and reducing the cost of hostaccelerator communication. As previous research has consistently shown, pure LBM simulations are able to achieve good performance results on heterogeneous systems thanks to the high parallel efficiency of this method. Unfortunately, when coupling LBM and IB methods, the overheads of IB degrade the overall performance. As an alternative, we have explored different hybrid implementations that effectively hide such overheads and allow us to exploit both the multi-core and the hardware accelerator in a cooperative way, with excellent performance results.

Introduction

The dynamics of a solid in a flow field is a research topic with a growing interest in many scientific communities. It is intrinsically interdisciplinary (structural mechanics, fluid mechanics, applied mathematics, . . .) and covers a broad range of applications (e.g. aeronautics, civil engineering, biological flows, etc.). The number of works in this field is rapidly increasing, which reflects the growing importance of studying the dynamics in the solid-fluid interaction [START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with movies boundaries[END_REF][START_REF] Calhoun | A cartesian grid method for solving the two-dimensional stream function-vorticity equations in irregular regions[END_REF][START_REF] Russell | A Cartesian grid method for modelling multiple moving objects in 2d incompressible viscous flows[END_REF][START_REF] Lima | Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method[END_REF]. Most of these simulations are compute-intensive and benefit from high performance computing systems. However, they also exhibit an irregular and dynamic behaviour, which often leads to poor performance when using emerging heterogeneous systems equipped with many-core accelerators.

Many computational fluid dynamic (CFD) applications and software packages have already been ported and redesigned to exploit heterogeneous systems. These developments have often involved major algorithm changes since some classical solvers may turned out to be inefficient or difficult to tune [START_REF] Valero-Lara | Block tridiagonal solvers on heterogeneous architectures[END_REF][START_REF] Valero-Lara | Fast finite difference Poisson solvers on heterogeneous architectures[END_REF]. Fortunately, other solvers are particularly well suited for GPU acceleration and are able to achieve significant performance improvements. The Lattice Boltzmann method (LBM) is one of those examples thanks to its inherently data-parallel nature. Certainly, the computing stages of LBM are amenable to fine grain parallelization in an almost straightforward way. This fundamental advantage of LBM has been consistently confirmed by many authors [START_REF] Bernaschi | A flexible high-performance lattice Boltzmann GPU code for the simulations of fluid flows in complex geometries[END_REF][START_REF] Rinaldi | A lattice-boltzmann solver for 3d fluid simulation on {GPU}[END_REF][START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF]11], for a large variety of problems and computing platforms.

In this paper, we explore the benefits of LBM solvers on heterogeneous systems. Our target application is an integrated framework that uses the Immersed Boundary (IB) method to simulate the influence of a solid immersed in a incompressible flow [11]. Some recent works that cover subjects closely related with our contribution are [START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF] and [START_REF] Layton | cuIBM -a GPU-accelerated immersed boundary method[END_REF]. In [START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF], authors presented an efficient 2D implementation of the LBM, which deals with geometries, by using curved boundaries-based methodologies, that is able to achieve a high performance on GPUs. Curved boundaries are taken into account via a non equilibrium extrapolation scheme developed in [START_REF] Guo | An extrapolation method for boundary conditions in lattice Boltzmann method[END_REF]. In [START_REF] Layton | cuIBM -a GPU-accelerated immersed boundary method[END_REF], S. K. Layton et al. studied the solution of two-dimensional incompressible viscous flows with immersed boundaries using the IB projection method introduced in [START_REF] Taira | The immersed boundary method: a projection approach[END_REF]. Their numerical framework is based on a Navier-Stokes solver and uses the Cusp library developed by Nvidia [START_REF] Dalton | Cusp: Generic parallel algorithms for sparse matrix and graph computations[END_REF] for GPU acceleration. Our framework uses a different Immersed Boundary formulation based on the one introduced by Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], which is able to deal with complex, moving or deformable boundaries [11,[START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Wu | Simulating 3d deformable particle suspensions using lattice Boltzmann method with discrete external boundary force[END_REF][START_REF] Huang | Simulation of flexible filaments in a uniform flow by the immersed boundary method[END_REF][START_REF] Zhu | Interaction of two flapping filament in a flow soap film[END_REF][START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-differences and finite-volume Navier-Stokes solvers[END_REF]. Some previous performance results were presented in [START_REF] Valero-Lara | Accelerating solid-fluid interaction using lattice-Boltzmann and immersed boundary coupled simulations on heterogeneous platforms[END_REF]. In this contribution we include a more elaborated discussion about performance results, and as a novelty, we explore alternative heterogeneous platforms based on the recently introduced Intel Xeon Phi device. Special emphasis is given to the implementation techniques adopted to mitigate the overhead of the immersed boundary correction and keep the solver highly efficient.

The remainder of this paper is organized as follows. In Section 2 we introduce the physical problem at hand and the general numerical framework that has been selected to cope with it. In Section 3 we review some optimizations of the baseline global LBM solver. After that, we describe the different optimizations and parallel strategies envisaged to achieve high-performance when introducing the IB correction. In Section 4 we focus on optimising this correction as it was a standalone kernel. Additional optimisations that take into account the interaction with the global LBM solver are studied afterwards in Section 5. Finally, we discuss the performance results of the proposed techniques in Section 6. We conclude in Section 7 with a summary of the main contributions of this work.

Numerical framework

As mentioned above, we have explored in this work a numerical framework based on the Lattice Boltzmann method coupled to the Immersed Boundary method. This combination is highly attractive when dealing with immersed bodies for two main reasons: (1) the shape of the boundary, tracked by a set of Lagrangian nodes is a sufficient information to impose the boundary values; and (2) the force of the fluid on the immersed boundary is readily available and thus easily incorporated in the set of equations that govern the dynamics of the immersed object. In addition, it is also particularly well suited for parallel architectures, as the time advancement is explicit and most computations are local [11]. In what follows, we briefly recall the basic formulation of LBM and then we describe the investigated LBM-IB framework.

The LBM method

Lattice Boltzmann has been extensively used in the past decades (see [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Numerical Mathematics and Scientific Computation)[END_REF] for a complete overview) and today is widely accepted in both academia and industry as a powerful and efficient alternative to classical Navier Stokes solvers for simulating (time-dependent) incompressible flows [11].

LBM is based on an equation that governs the evolution of a discrete distribution function f i (x, t) describing the probability of finding a particle at Lattice site x at time t with velocity v = e i . In this work, we consider the BGK formulation [START_REF] Bhatnagar | A model for collision processes in gases. I: Small amplitude processes in charged and neutral one-component system[END_REF] that relies upon an unique relaxation time toward the equilibrium distribution f (eq) i :

f i (x + e i t, t + t) -f i (x, t) = - t f (x, t) -f (eq) i (x, t) + tF i (1)
The particles can move only along the links of a regular lattice defined by the discrete speeds (e 0 = c(0, 0) ; e i = c(± 1, 0), c(0, ± 1), i = 1, . . ., 4 ; e i = c(± 1, ± 1), c(± 1, ± 1), i = 5, . . ., 8 with c = x/ t) so that the synchronous particle displacements x i = e i t never take the fluid particles away from the Lattice. For the present study, the standard two-dimensional 9-speed Lattice D2Q9 (Fig. 1) is used [START_REF] Aidun | Lattice-Boltzmann method for complex flows[END_REF], but all the techniques that are presented in this work can be easily applied to three dimensional lattices. The equilibrium function f (eq) (x, t) can be obtained by Taylor series expansion of the Maxwell-Boltzmann equilibrium distribution [START_REF] Qian | Lattice BGK models for Navier-Stokes equation[END_REF]:

f (eq) i = ω i 1 + e i • u c 2 s + (e i • u) 2 2c 4 s - u 2 2c 2 s (2)
In Eq. ( 2), c s is the speed of sound (c s = 1/ √ 3), is the macroscopic density, and the weight coefficients ω i are ω 0 = 4/9, ω i = 1/9, i = 1, . . ., 4 and ω 5 = 1/36, i = 5, . . ., 8 according to the current normalization. The macroscopic velocity u in Eq. ( 2) must satisfy a Mach number requirement |u |/c s ≈ M ≪ 1. This stands as the equivalent of the CFL number for classical Navier Stokes solvers.

F i in Eq. ( 1) represents the contribution of external volume forces at lattice level that in our case include the effect of the immersed boundary. Given any external volume force f ( ib), the contributions on the lattice are computed according to the formulation proposed in [START_REF] Guo | An extrapolation method for boundary conditions in lattice Boltzmann method[END_REF] as:

F i = 1 - 1 2 ω i e i -u c 2 s + e i • u c 4 s e i • f ib (3)
The multi-scale Chapman Enskog expansion of Eq. ( 1), neglecting terms of O( M 2 ) and using Eq. ( 3), returns the Navier-Stokes equations with body forces and the kinematic viscosity related to lattice scaling as = c 2 s ( -1/2) t. Without the contribution of the external volume forces stemming from the immersed boundary treatment, Eq. ( 1) is advanced forward in time in two main stages, known as collision and streaming, as follows:

1. Calculation of the local macroscopic flow quantities and u from the distribution functions f i (x, t):

= f i (x, t) and u = e i f i (x, t)
2. Collision:

f * i (x, t + t) = f i (x, t) - t f (x, t) -f (eq) i (x, t)
3. Streaming:

f i (x + e i t, t + t) = f * i (x, t + t)
A large number of researchers have investigated the performance aspects of this update process [START_REF] Rinaldi | A lattice-boltzmann solver for 3d fluid simulation on {GPU}[END_REF]11,[START_REF] Wellein | On the single processor performance of simple lattice boltzmann kernels[END_REF][START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF][START_REF] Habich | Performance engineering for the lattice Boltzmann method on GPGPUs: architectural requirements and performance results[END_REF][START_REF] Wittmann | Comparison of different propagation steps for lattice Boltzmann methods[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF] over the past decade. In Section 3 we review these previous works and describe the implementation techniques that we have opted to explore in our framework.

The LBM-IB framework

Next, we briefly introduce the Immersed Boundary method that we use both to enforce boundary values and to recover the fluid force exerted on immersed objects [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-differences and finite-volume Navier-Stokes solvers[END_REF]. In the selected Immersed Boundary approach (as in several others), the fluid is discretized on a regular Cartesian lattice while the immersed objects are discretized and tracked in a Lagrangian fashion by a set of markers distributed along their boundaries. The general setup of the investigated Lattice Boltzmann-Immersed Boundary method can be recast in the following algorithmic sketch.

1. Given f i (x, t) compute:

= f i (x, t) and u = e i f i (x, t) + t 2 f ib 2. Collision stage: f * i (x, t + t) = f i (x, t) - t f (x, t) -f (eq) i (x, t)
f * i (x, t + t) = f i (x, t) - t f (x, t) -f (eq) i (x, t) + tF i and f i (x + e i t, t + t) = f * i (x, t + t)
As outlined above, the basic idea is to perform each time step twice. The first one, performed without body forces, predicts the velocity values at the immersed boundary markers and the force distribution that restores the desired velocity boundary values at their locations. The second one applies the regularized set of singular forces and repeats the procedure (using Eq. ( 3)) to determine the final values of the distribution function at the next time step. A key aspect of this algorithm is the way by which the interpolation I and the S operators (termed as spread from now on) are applied. Here, following [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-differences and finite-volume Navier-Stokes solvers[END_REF], we perform both operations (interpolation and spread) through a convolution with a compact support mollifier meant to mimic the action of a Dirac's delta. Combining the two operators we can write in a compact form:

f (ib) (x, t) = 1 t U d (s, t + t) - û (y) ĩ (y -s) dy ĩ (x -s) ds (4)
where ĩ is the mollifier, is the immersed boundary, is the computational domain, and U d is the desired value on the boundary at the next time step. The discrete equivalent of Eq. ( 4) is simply obtained by any standard composite quadrature rule applied on the union of the supports associated to each Lagrangian marker. As an example, the quadrature needed to obtain the force distribution on the lattice nodes is given by:

f l ib x i , y j = Ne n=1 F l ib (Xn) ĩ x i -X n , y j -Y n n (5)
where the superscript l refers to the lth component of the immersed boundary force, (x i , y j ) are the lattice nodes (Cartesian points) falling within the union of all the supports, N e is the number of Lagrangian markers and n is a value to be determined to enforce consistency between interpolation and the convolution (Eq. ( 5)). More details about the method in general and the determination of the n values in particular can be found in [START_REF] Pinelli | Immersed-boundary methods for general finite-differences and finite-volume Navier-Stokes solvers[END_REF].

In what follows we will give more details on the construction of the support cages surrounding each Lagrangian marker since it plays a key role in the parallel implementation of the IB algorithm. Fig. 2 illustrates an example of the portion of the lattice units that falls within the union of all supports. As already mentioned, the embedded boundary curve is discretized into a number of markers X I , I = 1, . . ., N e . Around each marker X I we define a rectangular cage I with the following properties: (i) it must contain at least three nodes of the underlying Eulerian lattice for each direction; (ii) the number of nodes of the lattice contained in the cage must be minimized. The modified kernel, obtained as a Cartesian product of the one dimensional function [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF]

ĩ(r) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 6 5 -3|r| --3(1 -|r|) 2 + 1 0.5 ≤ |r| ≤ 1.5 1 3 1 + -3r 2 + 1 0.5|r| ≤ 0.5 0 otherwise (6)
will be identically zero outside the square I . We take the edges of the square to measure slightly more than three lattice spacings (i.e., the edge size is 3 + Á = 3 + Á in the actual LBM normalization). With such choice, at least three nodes of the lattice in each direction fall within the cage. Moreover a value of Á ł 1 ensures that the mollifier evaluated at all the nine (in two dimensions) lattice nodes takes on a non-zero value. The interpolation stage is performed locally on each nine points support: the values of velocity at the nodes within the support cage centered about each Lagrangian marker deliver approximate values (i.e., second order) of velocity at the marker location. The force spreading step requires information from all the markers, typically spaced = 1 apart along the immersed boundary. The collected values are then distributed on the union of the supports meaning that each support may receive information from supports centered about neighboring markers, as in Eq. ( 5). The outlined method has been validated for several test cases including moving rigid immersed objects and flexible ones [11].

Numerical validation of the method

Finally, we close this section presenting several test cases in order to validate the implementation of the code by comparing the numerical results obtained with other studies. One of the classical problems in CFD is the determination of the twodimensional incompressible flow field around a circular cylinder, which is a fundamental problem in engineering applications. Several Reynolds numbers (20, 40 and 100) have been tested with the same configuration. The cylinder diameter D is equal to 40. The flow space is composed by a mesh equal to 40D (1600) × 15D (600). The boundary conditions are set as: Inlet: u = U, v = 0, Outlet: ∂u ∂x = ∂v ∂x = 0, Upper and lower boundaries: ∂u ∂y = 0, v = 0, Cylinder surface: u = 0, v = 0, Volume fraction: 0.5236%. When Reynolds number is 20 and 40, there is no vortex structure formed during the evolution. The flow field is laminar and steady. In contrast, for the Reynolds number of 100, the symmetrical rectangular zones disappear and an asymmetric pattern is formed. The vorticity is shed behind the circular cylinder, and vortex structures are formed downstream. This phenomenon is graphically illustrated in Fig. 3.

Two important dimensionless numbers are studied, the drag

(C D = F D 0.5 U 2 D ) and lift (C L = F L 0.5 U 2 D ) coefficients. F D corresponds
to the resistance force of the cylinder to the fluid in the streamwise direction and F L is the lifting force of the circular cylinder, is the density of the fluid, and U is the velocity of inflow. In order to verify the numerical results, the coefficients were calculated and compared with the results of previous studies (Table 1). The

Implementation of the global LBM update

Parallelization strategies

Parallelism is abundant in the LBM update and can be exploited in different ways. On our GPU implementation, the lattice nodes are distributed across GPU cores using a fine-grained distribution. As shown in Fig. 4 (bottom), we used a 1D Grid of 1D CUDA Block. Each CUDA-thread performs the LBM update of a single lattice node. On multi-core processors, cache locality is a major performance issue and it is better to distribute the lattice nodes across cores using a 1D coarse-grained distribution (Fig. 4, top-left). The cache coherence protocol keeps the boundaries between subdomains updated. On the Intel Xeon Phi, we also distribute the lattice nodes across cores using a 1D coarse-grained distribution, but using a smaller block size (Fig. 4, top-right).

Another important issue is how to implement a single LBM update. Conceptually, a discrete time step consists of a local collision operation followed by a streaming operation that propagates the new information to the neighbour lattice nodes. However, most implementations do not apply those operations as separate steps. Instead, they usually fuse in a single loop nest (that iterates over the entire domain), the application of both operations to improve temporal locality [START_REF] Wellein | On the single processor performance of simple lattice boltzmann kernels[END_REF][START_REF] Habich | Performance engineering for the lattice Boltzmann method on GPGPUs: architectural requirements and performance results[END_REF].

This fused loop can be implemented in different ways. We have opted to use the pull scheme introduced in [START_REF] Wellein | On the single processor performance of simple lattice boltzmann kernels[END_REF]. In this case, the body of the loop performs the streaming operation before collision, i.e. the distribution functions are gathered (pulled) from the neighbours before computing the collision. Algorithm 1 shows a sketch of our implementation. Algorithm 1. LBM pull.

1:

Pull (f1,f2,ϖ,cx,cy) 2:

x, y 3:

xstream, ystream 4: localu x ,localu y ,local 5:

local f [9], feq 6:

for i = 1 →9 do 7: xstream = x -cx[i] 8: ystream = y -cy[i] 9: local f [i] = f1[xstream][ystream][i] 10:
end for 11:

for i = 1 →9 do 12: local + = local f [i] 13: localu x + = cx[i] × local f [i] 14: localu y + = cy[i] × local f [i] 15:
end for 16:

localu x = localu x /local 17:
localu y = localu y /local 18:

for i = 1 →9 do 19: feq = ω[i] • • (1 + 3 • (cx[i] • localu x + cy[i] • localu y ) + (cx[i] • localu x + cy[i] • localu y ) 2 -1.5 × ((localu x ) 2 + (localu y ) 2 )) 20: f2[x][y][i] = local f [i] • (1 -1 ) + feq • 1 21:
end for

Other implementations [START_REF] Bernaschi | A flexible high-performance lattice Boltzmann GPU code for the simulations of fluid flows in complex geometries[END_REF][START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF][START_REF] Schnherr | Multi-thread implementations of the lattice Boltzmann method on non-uniform grids for CPUs and GPUs[END_REF] have used the traditional implementation sketched in Algorithm 2, which performs the collision operation before the streaming. It is known as the "push" scheme [START_REF] Wellein | On the single processor performance of simple lattice boltzmann kernels[END_REF] since it loads the distribution function from the current lattice point and then it "pushes" (scatters) the updated values to its neighbours. Algorithm 2. LBM push.

1:

Collide Stream (ux,uy, ,f1,f2,ϖ,cx,cy) 2:

x, y 3:

xstream, ystream 4: feq 5:

for i = 1 →9 do 6:

feq = ϖ[i] • • (1 + 3 • (cx[i] • ux[x][y] + cy[i] • uy[x][y]) + (cx[i] • ux[x][y] + cy[i] • uy[x][y]) 2 -1.5 × ((ux[x][y]) 2 + (uy[x][y]) 2 )) 7: f1[x][y][i] = f1[x][y][i] • (1 -1 ) + feq • 1 8: xstream = x + cx[i] 9: ystream = y + cy[i] 10: f2[xstream][ystream][i] = f1[x][y][i] 11:
end for 1:

Macroscopic (ux,uy, ,f2,cx,cy) 2:

x, y 3: localu x ,localu y ,local 4:

for i = 1 →9 do 5: local + = f2[x][y][i] 6: localu x + = cx[i] × f2[x][y][i] 7: localu y + = cy[i] × f2[x][y][i] 8:
end for 9:

[

Data layout and memory management

Since LBM is a memory-bound algorithm, another important optimization problem is to maximize data locality. Many groups have considered this issue and have introduced several data layouts and code transformations that are able to improve locality on different architectures [11,[START_REF] Wellein | On the single processor performance of simple lattice boltzmann kernels[END_REF][START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF][START_REF] Habich | Performance engineering for the lattice Boltzmann method on GPGPUs: architectural requirements and performance results[END_REF][START_REF] Wittmann | Comparison of different propagation steps for lattice Boltzmann methods[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF][START_REF] Shet | Data structure and movement for lattice-based simulations[END_REF][START_REF] Shet | On vectorization for lattice based simulations[END_REF]. Here we briefly describe the strategies used by our framework.

Different data structures have been proposed to store the discrete distribution functions f i in memory:

• AoS. This data structure stores all the discrete distribution functions f i of a given lattice point in adjacent memory positions (see Fig. 5(a)). This way, it optimizes locality when computing the collision operation [START_REF] Shet | On vectorization for lattice based simulations[END_REF]. However, it does not provide good performance on GPU architectures since it leads to poor bandwidth utilization [START_REF] Bernaschi | A flexible high-performance lattice Boltzmann GPU code for the simulations of fluid flows in complex geometries[END_REF][START_REF] Rinaldi | A lattice-boltzmann solver for 3d fluid simulation on {GPU}[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF]]. • SoA. In this alternative data structure, the discrete distribution functions f i for a particular velocity direction are stored sequentially in the same array (see Fig. 5(b)). Since each GPU thread handles the update of a single lattice node, consecutive GPU threads access adjacent memory locations with the SoA layout [START_REF] Bernaschi | A flexible high-performance lattice Boltzmann GPU code for the simulations of fluid flows in complex geometries[END_REF][START_REF] Rinaldi | A lattice-boltzmann solver for 3d fluid simulation on {GPU}[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF]. This way, they can be combined (coalesced) into a single memory transaction, which is not possible with the AoS counterpart. • SoAoS. We have also explored a hybrid data structure, denoted as SoAoS in [START_REF] Shet | On vectorization for lattice based simulations[END_REF]. As SoA, it also allows coalesced memory transactions on GPUs. However, instead of storing the discrete distribution functions f i for a particular velocity direction in a single sequential array, it distributes them across different blocks of a certain block size (see Fig. 5(c)). This size is a tunable parameter that trades off between spatial and temporal locality.

Apart from the data layout, the memory management of the different implementations may also differ in the number of lattices that are used internally. We have used a two-lattice implementation, which is denoted as the AB scheme in [START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF]. Essentially, AB holds the data of two successive time steps (A and B) and the simulation alternates between reading from A and writing to B, and vice versa. Other proposals, such as the AA data layout in [START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF], are able to use a single copy of the distributions arrays and reduce the memory footprint. Some previous works have shown that those single lattice schemes outperform the AB scheme on multi-core processors (AA achieved the best results in [START_REF] Wittmann | Comparison of different propagation steps for lattice Boltzmann methods[END_REF]). However, on GPUs the performance benefits of these schemes are less clear. In fact, a recent work has shown that both schemes get similar performance on GPUs [START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF] or the AB scheme is typically a little faster on the latest GPUs. On the other hand, AB simplifies the integration with IB correction, which is the main focus of this research, and is therefore clearly preferred in our framework.

Implementation of the IB correction

In this section we discuss the different strategies that we have explored to optimize and parallelize the IB correction, i.e. the computations related with the Lagrangian markers distributed on the solid(s) surface(s).

Parallelization strategies

The parallelization of this correction step is only effective as long as there is sufficient work, which depends on the number of Lagrangian points in the solid surface and the overlapping among their supports.

On multi-core processors, we have been able to achieve satisfactory efficiencies (even for moderate number of points) using a coarse-grained distribution of (adjacent) Lagrangian points across all cores. This implementation is relatively straightforward using a few OpenMP pragmas that annotate the loops that iterate over the Lagrangian points.

On GPUs, it is better to exploit fine-grained parallelism. Our implementation consists of two main kernels denoted as IF (Immersed Forces) and BF (Body Forces) respectively. Both of them use the 1D distribution of threads across Lagrangian points shown in Fig. 6.

The IF kernel performs three major steps (see Algorithm 3). First, it assembles the velocity field on the supports. Then, it undertakes the interpolation at the Lagrangian markers and finally it determines the Eulerian volume force field on each node of the union of the supports. Note that forces are updated with atomic operations. Despite the overhead of such operations, they are necessary Fig. 6. 1D fine-grained distribution of Lagrangian points across CUDA threads.

to prevent race conditions since the supports of different Lagrangian points can share the same Eulerian points, as graphically shown in Fig. 2.

Once the IF kernel has been completed, the BF kernel (see Algorithm 4) computes the lattice forces and apply a local LBM update on the union of the supports, including now the contribution of the immersed boundary forces. Again, it avoids race conditions using atomic operations. Algorithm 3. Pseudo-code of the IF kernel. F body (Body Force),x, y,velx, vely 3:

for i = 1 → numSupport do 4: x = s.Xsupport[i] 5: y = s.Ysupport[i] 6: velx = s.VelXsupport[i] 7: vely = s.VelYsupport[i] 8:
for j = 1 →9 do 9:

F body = (1 -0.5 • 1 ) • w[j] • (3 • ((cx[j] -velx) • (fx[x][y]) + (cy[j] - vely) • fy[x][y])) 10: AddAtom(f n+1 [j][x][y], F body ) 11:
end for 12:

end for

Data layout

Once again, another key aspect affecting performance in this correction step is the data structure (data layout) that is used to store the information about the Lagrangian points and their supports (coordinates, velocities and forces). Using the locality principle, it is natural to use an array of structures (AoS) to place all this information at a given Lagrangian point in nearby memory locations. This is the data structure that we have used on multi-core processors, since it optimizes cache performance. In contrast, on GPUs it is more natural to use a SoA data structure that distributes the information of all Lagrangian points in a set of one-dimensional arrays. With SoA, consecutive threads access to contiguous memory locations and memory accesses are combined (coalesced) in a single memory transaction. Fig. 7 illustrates the difference between those layouts in a simplified example with only two Lagrangian points.

Coupled LBM-IB on heterogeneous platforms

The complete LBM-IB framework can be fully implemented on many core processors combining the strategies described earlier for the global LBM updated and the IB correction. In this approach, the host processor stays idle most of the time and it is used exclusively for:

• Pre-processing. The host processor sets up the initial configuration and uploads those initial data to the main memory of the accelerator. • Monitoring. The host processor is also in charge of a monitoring stage that downloads the information of each lattice node (i.e., velocity components and density) back to the host main memory when required. On GPUs, this approach consists of three main kernels, denoted in Fig. 8 as LBM, IF and BF respectively, which are launched consecutively for every time step. The first kernel implements the LBM update while the other two perform the IB correction. The overhead of the pre-processing stage performed on the multi-core processor has been experimentally shown to be negligible and the data transfer of the monitoring stage are mostly overlapped with the execution of the LBM kernel.

Although this approach achieves satisfactory results on GPUs, its speedups are substantially lower than those achieved by pure LBM solvers [START_REF] Rinaldi | A lattice-boltzmann solver for 3d fluid simulation on {GPU}[END_REF][START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF][START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF]11]. The obvious reason behind this behavior is the ratio between the characteristic volume fraction and the fluid field, which is typically very small. Therefore, the amount of data parallelism in the LBM kernel is substantially higher than in the other two kernels. In fact, for the target problems investigated, millions of threads compute the LBM kernel, while the IF and BF kernels only need thousands of them. But in addition, those kernels also require atomic functions due to the intrinsic characteristics of the IB method, and those operations usually degrade performance.

LBM-IB on hybrid multicore-GPU platforms

As an alternative to mitigate the overheads caused by the IB correction, we have explored hybrid implementations that take advantage of the host Xeon processor to hide them.

In Fig. 9 we show the hybrid implementation that we proposed on a previous paper [START_REF] Valero-Lara | Accelerating solid-fluid interaction using lattice-Boltzmann and immersed boundary coupled simulations on heterogeneous platforms[END_REF] for heterogeneous multicore-GPU platforms. The LBM update is performed on the GPU as in the previous approach. However, IB and an additional local correction to LBM on the supports of the Lagrangian points are performed on the host processor. Both steps are coordinated using a pipeline. This way, we are able to overlap the prediction of the fluid field for the t + 1 iteration with the correction of the IB method on the previous iteration t. This is possible at the expense of a local LBM computation of the"t + 1 ′′ iteration on the multi-core and additional data transfers of the supports between the GPU and the host processor at each simulation step. Fig. 9. Hybrid multicore-GPU implementation. The LBM update is performed on the GPU, whereas the IB correction and an additional step to update the supports of the Lagrangian points are performed on the multi-core processor.

LBM-IB on hybrid multicore-Xeon Phi platforms

The homogeneous Xeon Phi implementation suffers from the same performance problems that the homogeneous GPU counterpart. But again, multicore-Xeon Phi collaboration allows us to achieve higher performance. It is possible to use the same hybrid strategy introduced above for the multicore-GPU platform. However, we have opted to use a slightly different implementation that simplifies code development. Fig. 10 graphically illustrates the new partitioning. Essentially, it also consists on splitting the computational domain into two subdomains, so that one of them (the solid subdomain) fully includes the immersed solid. With such distribution, it is possible to perform the IB correction exclusively on the host Intel Xeon processor. However, in this case the host also performs the global LBM update on that subdomain, which simplifies the implementation.

As shown in Fig. 10, every time step, it is necessary to exchange the boundaries between both subdomains. To reduce the penalty of such data transfers, we update the boundaries between subdomains at the beginning of each time step. With this transformation, it is possible to exchange those boundaries with asynchronous operations that are overlapped with the update of the rest of the subdomain. In our target simulations, the solid subdomain is much smaller that the Xeon Phi counterpart. To improve performance, the size of both subdomains is adjusted to balance the loads between both processors.

Performance evaluation

Experimental setup

To critically evaluate the performance of the developed LBM-IB solvers, we have considered next a number of tests executed on two different heterogeneous platforms, whose main characteristics are summarized in Table 2.

On the GPU, the on-chip memory hierarchy has been configured as 16 KB shared memory and 48 KB L1, since our codes do not benefit from a higher amount of shared memory on the investigated tests.

Simulations have been performed using double precision, and we have used the conventional MFLUPS metric (millions of fluid lattice updates per second) to assess the performance.

Standalone Lattice-Boltzmann update

Before discussing the performance of our LBM-IB framework, it is important to determine the maximum performance that we can attain. We can estimate an upper limit omitting the IB correction. Fig. 10. Hybrid multicore-Phi implementation. The Lx × Ly lattice is split into two 2D sub-domains so that the IB correction is only needed on one of the domains. The multi-core processor updates the LxIB × Ly subdomain, which fully includes the immersed solid (marked as a circle). The Xeon Phi updates the rest of the lattice nodes (the LxLBM × Ly subdomain). The grey area highlights the ghost lattice nodes at the boundary between both sub-domains. In our target simulations, LxLBM ≫ LxIB. Fig. 11 shows the performance of this benchmark on a Kepler GPU.

Recall from Section 3 that our implementation is based on the pull scheme and uses two lattice with the SoA data layout. As a reference we also show the performance of the Sailfish software package [START_REF] Januszewski | A flexible multi-GPU implementation of the lattice Boltzmann method[END_REF], which includes a LBM solver based on the push scheme, and the following estimation of the ideal MFLUPS [START_REF] Shet | Data structure and movement for lattice-based simulations[END_REF]:

MFLUPS ideal = B × 10 9 10 6 × n × 6 × 8 (7)
where B × 10 9 is the memory bandwidth (GB/s), n depends on LBM model (DxQn), for our framework n = 9, D2Q9. The factor 6 is for the memory accesses, three read and write operations in the spreading step and three read and write operations in the collision step, and the factor 8 is for double precision (8 bytes). Fig. 12(a) and (b) shows the performance on an Intel Xeon server. Although on multi-core processor it is natural to use the AoS data layout, SoA and SoAoS (with a block size of 32 elements) turn to be the most efficient data layouts. The main reason behind these unexpected results lies on the vector capabilities of modern processors. The compiler has been able to vectorize the main loops of the LBM update and both the AoS and SoAoS layouts allow a better exploitation of vectorization. Fig. 13 shows the scalability on this solver. The observed speedup factors over the sequential implementation almost peak with 16 threads, but the application scales relatively well since its performance is limited by the memory bandwidth.

Fig. 14(a) and (b) shows the performance on the Intel Xeon Phi. In this platform, SoAoS is able to outperform the other data layouts. The optimal SoAoS block size in this case is 128 elements, 4× the block size of the Intel Xeon, which coincides with the ratio between the vector widths of both architectures. Fig. 15 analyses the impact of the thread-core affinity on performance using three different pre-defined strategies (see [START_REF] Green | OpenMP thread affinity control[END_REF] for details on the rationale of each strategy). The compact affinity provides the best performance although the differences between them are not significant.

Finally, Fig. 16 shows the scalability on the Xeon Phi. For large problems, performance always improves as the number of threads increases. However, the gap with the ideal MFLUPS estimate is much larger in this case than in the other platforms. Given that the codes used for the Intel Xeon Phi and those used for the Intel Fig. 12. Performance of the LBM update on the Intel Xeon multi-core processor. Xeon are essentially the same (only optimal block sizes and minor optimization parameters change between both implementations) it is expected that either memory bandwidth is not limiting performance in this platform, or better performance values are expected for larger problem sizes.

Overall, the best performance is achieved on the GPU, which is able to outperform both the Intel Xeon and the Intel Xeon Phi. For the largest grid tested, the speed factor is 4.62× and 1.22×, respectively.

Standalone IB correction

Our second test uses a synthetic benchmark that focus exclusively on the IB correction, i.e., it omits the global LBM update and only applies the local IB correction. Fig. 17 shows the speedups of the parallel implementations over a sequential counterpart for increasing number of Lagrangian nodes. We are able to achieve substantial speedups, even for a moderate number of nodes. Notably, the best performance is achieved on the GPU, despite the overheads caused by the atomic updates needed on that implementation. From 2500 Lagrangian markers, the GPU outperforms the 8-core multicore counterpart with a 2.5× speedup. Fig. 17. Speedups of the IB method on multi-core and GPU for increasing number of Lagrangian nodes.

6.4. Coupled LBM-IB on heterogeneous platforms 6.4.1. Multicore-GPU Fig. 18(a) shows the performance of the hybrid LBM-IB solver on a multicore-GPU heterogeneous platform for an increasing number of lattice nodes. As a reference, Fig. 18(b) shows the performance of the homogeneous GPU implementation. We have used the same physical setting studied in Section 2 and we have investigated two realistic scenarios with characteristic volume fractions of 0.5% and 1% respectively (i.e. the amount of embedded Lagrangian markers also grows with the number of lattice nodes).

The performance (MFLUPS) of the homogeneous implementation (Fig. 18(b)) drops substantially over the pure LBM implementation (Fig. 11). The slowdown is around 15% for a solid volume fraction of 0.5%, growing to 25% for the 1% case. In contrast, for these fractions the hybrid GPU-multicore version is able to hide the overheads of the IB method, having similar performance to the pure LBM implementation.

Fig. 19 shows a breakdown of the execution time for the simulation with a solid volume fraction of 1%. On the hybrid implementation, the cost of the IB correction is higher than in the homogeneous counterpart, reaching around 65% of the total execution time. This is expected since the in this case, the IB correction includes additional data transfers and local LBM updates. However, these costs are hidden with the global LBM update.

Multicore-Xeon Phi

Similarly, the performance of the proposed strategy over the multicore-Xeon Phi heterogeneous platform is analyzed in Fig. 18. Performance of the complete LBM-IB solver for an increasing number of lattice nodes. Fig. 20(a). We have focused on the same numerical scenario described earlier with a solid volume fraction of 1%. As in the multicore-GPU platform, the observed performance in MFLUPS roughly matches the performance of the pure-LBM implementation on the same platform (Fig. 14(b)), with a peak performance of 450 MFLUPS for the largest problem size. L x IB, which defines the size of the subdomain that is simulated on the multi-core Fig. 19. Execution time breakdown for a solid volume fraction of 1% of the LBM and IB kernels on the homogeneous GPU implementation (left) and multicore-GPU heterogeneous (right) platforms. Fig. 20. Performance of the complete LBM-IB solver on the multicore-Phi platform as the size of subdomain that is simulated on the multi-core procesor increases.

processor, is a critical parameter in this implementation. As expected, there is an optimal value, which depends on the size of both the immersed solid and the grid, that balances the load in both processors. For example, the peak performance for the smallest grid tested (L x = 2560, L y = 800) is attained when the Xeon subdomain is roughly a 10% of the complete grid (L x IB = 250), growing to 12% for the largest grid tested (L x = 5120, L y = 1400; (L x IB = 620)).

Overall, as we have aimed, the overhead of the solid interaction is mainly hidden thanks to the effective cooperation with the multicore processor. The penalty of the data transfers between both processors is negligible since boundary exchanges are conveniently orchestrated to allow their overlapping with useful computations.

Finally, as a reference, Fig. 20(b) shows the performance attained by an equivalent LBM-IB implementation, but without overlapping the execution on the Intel Xeon and the Intel Xeon Phi, i.e. the Xeon and the Phi subdomains are updated sequentially, one after the other. In this case, the peak performance drops to 300 MFLUPS, which highlights the benefit of overlapping the execution of both processors. Obviously, the lower L x IB, the higher performance we achieve with this synthetic code since it increases the amount of work delivered to the Intel Xeon Phi.

Conclusions

In this paper, we have investigated the performance of a coupled Lattice-Boltzmann and Immersed Boundary method that simulates the contribution of solid behavior within an incompressible fluid. While LBM has been widely studied on heterogeneous platforms, the coupling of LBM with Immersed-Boundary methods has received less attention.

We have reviewed different strategies to enhance the performance on three state-of-the-art parallel architectures: an Intel Xeon server, a NVIDIA KeplerGPU and an Intel Xeon Phi accelerator. Our baseline LBM solver uses most of the state-of-the-art code transformations that have been described in previous work. Notably, performance results (MFLUPS) on the GPU and the multicore processor are close to ideal MFLUPS estimations. On the Xeon Phi, the gap with these ideal estimations is higher, but it also achieves competitive performance. Overall, the best results are achieved on the GPU, which peaks at 550 MFLUPs, whereas the Intel Xeon Phi peaks at 450 MFLUPs.

Our main contribution is the design and analysis of a hybrid implementation that takes advantage of both the accelerator (GPU/Xeon Phi) and multi-core in a co-operative way to solve the coupled LBM-IB problem. For interesting physical scenarios with realistic solid volume fractions, these hybrid solvers are able to hide the overheads caused by the IB correction and match the performance (in terms of MFLUPS) of state-of-the-art pure LBM solvers. This has been possible thanks to the effective cooperation between the accelerator and the multi-core processor and the overlapping of data transfers between their memory spaces with useful computations.

Based on the presented LBM-IB framework, we envision two main research lines as future work. First, we will investigate an adaptation of the framework to problems which deal with deformable and moving bodies. In this type of scenarios, the work partitioning and assignment among computational resources will need to follow a dynamic approach instead of a static strategy as that implemented in this paper, to adapt the workload to the specific configuration at each timestep. Second, we plan to extend our work to 3D simulations; for these problems, memory consumption arises as one of the main problems, naturally driving to distributed-memory architectures, possibly equipped with hardware accelerators. In this case, we believe many of the techniques presented in this work will be of direct application at each node level; the scalability and parallelization strategies across nodes is still a topic under study.

Fig. 1 .

 1 Fig. 1. Standard two-dimensional 9-speed lattice (D2Q9) used in our work.

Fig. 2 .

 2 Fig. 2. An immersed curve discretized with Lagrangian points (•). Three consecutive points are considered with the respective supports.

Fig. 3 .

 3 Fig. 3. Vorticity with Re = 100.

Fig. 4 .

 4 Fig. 4. Fine-grained and coarse-grained distributions of the lattice notes.

  x][y] = local 10: ux[x][y] = localu x /local 11: uy[x][y] = localu y /local

Fig. 5 .

 5 Fig. 5. Different data layouts to store the discrete distribution functions f i in memory.

Fig. 7 .

 7 Fig.7. Different data layouts used to store the information about the coordinates, velocities and forces of the Lagrangian points and their supports. On multi-core processors we use an AoS data structure (top) whereas on GPUs (bottom) we use a SoA data structure.

Fig. 8 .

 8 Fig. 8. Homogeneous GPU Implementation. Both the LBM update and the IB correction are performed on the GPU. The host processors stays idle most of the time.

Fig. 11 .

 11 Fig. 11. Performance of the LBM update on the NVIDIA Kepler GPU.

Fig. 13 .

 13 Fig.[START_REF] Guo | An extrapolation method for boundary conditions in lattice Boltzmann method[END_REF]. Scalability of the LBM update on an Intel Xeon multi-core processor.

Fig. 14 .

 14 Fig. 14. Performance of the LBM update on the Intel Xeon Phi.

Fig. 15 .

 15 Fig. 15. Performance of the LBM update with different thread-core affinity strategies on the Intel Xeon Phi.

Fig. 16 .

 16 Fig. 16. Scalability of the LBM update on an Intel Xeon Phi.

Table 1

 1 Comparison between the numerical results yield by our method and previous studies. Reynolds number of 20 and 40 is equal to the results presented by Zhou et al.[START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF]. The drag coefficient obtained for Reynolds number of 100 is identical to the results obtained by Silva et al.[START_REF] Lima | Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method[END_REF], and the lift coefficient is close to the presented by Zhou et al.[START_REF] Zhou | {GPU} implementation of lattice Boltzmann method for flows with curved boundaries[END_REF].

	Author	Re = 20	Re = 40	Re = 100
		CD	CD	CD	CL
	Calhoun [2]	2.19	1.62	1.33	0.298
	Rusell [3]	2.22	1.63	1.34	-
	Silva et al. [4]	2.04	1.54	1.39	-
	Xu [1]	2.23	1.66	1.423	0.34
	Zhou et al. [9]	2.3	1.7	1.428	0.315
	This work	2.3	1.7	1.39	0.318

Table 2

 2 Summary of the main features of the platforms used in our experimental evaluation.

	Platform	Intel Xeon	NVIDIA GPU	Intel Xeon Phi
	Model	E5520/E5-2670	Geforce GTX780 (Kepler)	5510P
	Frequency	2.26/2.6 GHz	0.863	1.053 Ghz
	Cores	8/16	2304	60
	On-chip Mem.	L1 32 KB (per core)	SM 16/48 KB (per MP)	L1 32 KB (per core)
		L2 512 KB (unified)	L1 48/16 KB (per MP)	L2 256 KB (per core)
		L3 20 MB (unified)	L2 768 KB (unified)	L2 30 MB (coherent)
	Memory	64/32GB DDR3	6GB GDDR5	8 GB GDDR5
	Bandwidth	51.2 GB/s	288 GB/s	320 GB/s
	Compiler	Intel Compiler 14.0.3	nvcc 6.5.12	Intel Compiler 14.0.3
	Compiler Flags	-O3 -fomit-frame-pointer -fopenmp	-O3 -arch = sm 35	1
	1 -O3 openmp -mcmodel=medium fno-alias -mP2OPT hlo use const pref dist=64 -mP2OPT hlo use const second pref dist=32 ′′

Acknowledgments

This research is supported by the Spanish Government Research Contracts TIN2012-32180, Ingenio 2010 Consolider ESP00C-07-20811, MTM2013-40824 and SEV-2013-0323, by the EU-FET grant EUNISON 308874 and by the Basque government BERC 2014-2017 contract. We also thank the support of the computing facilities of Extremadura Research Centre for Advanced Technologies (CETA-CIEMAT) and NVIDIA GPU Research Center program for the provided resources.


K Saleh

A Fernandez

O Llopis

Investigations on an All-Tunable Fiber Ring Resonator

Keywords: fiber ring resonator, tunable optical filters, frequency reference, microwave photonics, CAD of microwave optical systems

The architecture of an all-tunable optical fiber ring resonator is described in detail in this paper. This architecture has been firstly modeled using an original CAD approach. The simulation results demonstrate a total control of both the absolute frequency and the free spectral range of the final optical resonance comb generated by the optical resonator. The different experimental setups used to characterize the tunable resonator are described and the obtained results proving the concept are also provided.

Introduction

Tunability is a desired feature in any optical filter [START_REF] Sadot | Tunable optical filters for dense WDM networks[END_REF]. This also applies to the fiber ring resonator (FRR) case due to the high optical quality factor (QOpt) that can be obtained and to the simplicity in its fabrication [START_REF] Saleh | Fiber ring resonator based opto-electronic oscillator: phase noise optimisation and thermal stability study[END_REF]. Indeed, the FRR is a very useful element for many applications as it can be used as a high selectivity add-drop optical filter, a high resolution gyroscope, a fiber laser etc. We can particularly mention the phase noise performance that can be achieved for microwave signals generated by an optoelectronic oscillator (OEO) when a high QOpt FRR is used as the frequency stability element in the oscillator setup [START_REF] Saleh | Optical scattering induced noise in fiber ring resonators and optoelectronic oscillators[END_REF]. Besides the importance of achieving high spectral purity for the generated microwaves, the tunability is also greatly desired in OEOs. For that reason, numerous studies have been performed to tune the oscillation frequency in different OEO's architectures [START_REF] Ozdur | Low noise optically tunable opto-electronic oscillator with fabry-perot etalon[END_REF][START_REF] Savchenkov | Voltage-controlled photonic oscillator[END_REF][START_REF] Maxin | Widely tunable opto-electronic oscillator based ona dual-frequency laser[END_REF].

Certainly, tuning the FRR's frequency response is possible by changing the resonator's temperature or by using the effect of polarization and birefringence in the resonance loop [START_REF] Seraji | Tunability of optical filter based on fibre-optic ring resonator using polarisation and birefringence effects in the resonator loop[END_REF]. On the one hand, the FRR's response to the temperature is too weak and it is also very slow. In a previous work, [START_REF] Saleh | Fiber ring resonator based opto-electronic oscillator: phase noise optimisation and thermal stability study[END_REF], we have performed a detailed study on this matter. It has been found, for example, that the FSR response to the temperature in a 20m-long FRR is equal to |70𝐻𝑧|/𝐾. This means that a temperature variation of ∆𝑇 = 1.4 × 10 5 𝐾 is theoretically required to tune the FRR's frequency response by a full FSR (FSR ≈ 10 MHz in a 20m-long FRR), which is clearly impractical to achieve. On the other hand, the effect of polarization and birefringence is limited and could be sensitive to external perturbations.

In this paper, we first give a brief description of the FRR and later on we present our theoretical and experimental investigations on a particular architecture of the FRR. This architecture is principally based on the enclosure of a phase shifter and a tunable optical coupler (forming a Sagnac mirror) at the same time in a double coupler add-drop FRR. While some research works studying separately similar architectures were reported in the literature [START_REF] Sadot | Ultra-fast tunable fiber-loop optical filters for dense WDM applications[END_REF][START_REF] Vazquez | Tunable optical filters using compound ring resonators for DWDM[END_REF][START_REF] Vazquez | Sagnac loop in ring resonators for tunable optical filters[END_REF][START_REF] Seraji | Tunable optical filter based on Sagnac phase-shift using single optical ring resonator[END_REF], to our knowledge, no detailed study has been performed to show the effect of the combination of both the phase shifter and the tunable optical coupler on the double coupler adddrop FRR's response. In the present work, our investigations show that such combination allows a complete tuning of the absolute frequency of the generated final optical resonance comb (ORC) and of its free spectral range (FSR) as well. In addition, the results of a study performed on the possibility of tuning independently both of these characteristics of an ORC are also addressed.

Fiber ring resonator

The architecture of the double-direct-coupled FRR, considered in our studies, is shown in Fig. 1(a). This FRR is fabricated using two low loss two-by-two directional optical fiber couplers (C1 and C2) linked together with a single-mode fiber loop. The resonator generates a transverse single ORC with microwave spacing called FSR. The FSR is given by FSR = c/nL, where (L) the fiber loop length, (n) is its refractive index and (c ~ 3.10 8 m/s) is the light speed in free space. The QOpt of this resonator particularly depends on the couplers' coupling coefficients (κ1 for C1 and κ2 for C2) and on the overall losses in the resonator architecture: couplers' excess loss (γ1,2), splices' loss (αs) and fiber linear loss (αf).

Absolute frequency tuning using an optical phase-shifter

As already reported in [START_REF] Sadot | Ultra-fast tunable fiber-loop optical filters for dense WDM applications[END_REF], it is possible to tune the absolute frequency of an ORC generated by a FRR by inserting an optical phase-shifter (PS) within the FRR's loop [see Fig. 1(a)]. The PS enables the control of the phase delay inside the ring loop and, therefore, the control of the phase matching condition. An additional phase-shift added inside the FRR's loop (ΔφPS) will lead to a proportional absolute frequency shift (ΔfOpt; fOpt being the optical resonant frequency) in the ORC as follows, [START_REF] Sadot | Ultra-fast tunable fiber-loop optical filters for dense WDM applications[END_REF]: [START_REF] Saleh | Fiber ring resonator based opto-electronic oscillator: phase noise optimisation and thermal stability study[END_REF] 2

PS PS Opt c f FSR nL         
From the above relation, we can clearly notice that when an additional 2π rad phase-shift is added to the optical signal inside the FRR's loop, the resonant frequency will shift by one complete FSR. In order to simulate the response of such FRR architecture, we have used an upgraded version of a model of the FRR developed in a previous work [START_REF] Saleh | Study of the noise processes in microwave oscillators based on passive optical resonators[END_REF]. This FRR model is based on an original approach using electronic design software (ADS from Agilent), typically used for microwave systems' design. With this FRR's model [depicted in Fig. 1(b)], we have been able to simulate the transmission response of a 20m-long FRR including a phase-shifter. Consequently, the abovementioned frequency tuning technique has been confirmed by changing ΔφPS from 0 to 2π while monitoring the FRR's responses (see Fig. 2; for clarity's sake, the presented results are for a change in ΔφPS's value by π). Furthermore, we have found that such a tuning technique does not modify the FRR's transmission characteristics (insertion loss (IL), QOpt, etc…). The abovementioned theory on the optical absolute frequency shift has been experimentally confirmed using a narrow linewidth laser and the wavelength scanning technique [START_REF] Llopis | High-Q optical resonators: characterization and application to stabilization of lasers and high spectral purity microwave oscillators[END_REF]. This wavelength scanning technique consists of a frequency fine-tuning of a narrow linewidth (~1 kHz) Koheras fiber laser to explore the FRR's ORC by recording the FRR's response on a photodiode followed by an oscilloscope. In our experimental setup, we have used a 23m-long FRR including a PS (abbreviated by FRR-PS) and a second 20m-long FRR as a spectator to witness the absolute frequency shift in the FRR-PS's ORC (see Fig. 3). The results are depicted in Fig. 4 and they confirm that when an additional phase-shift ΔφPS = π is added inside the FRR-PS's loop, its ORC shifts by FSR/2.

Free spectral range tuning using a variable optical coupler

A fiber loop mirror (FLM), or a Sagnac reflector [START_REF] Birks | Jones calculus analysis of single-mode fiber Sagnac reflector[END_REF], is a device made by connecting the ports on one side of a two-by-two fibred directional optical coupler. This FLM is usually used in its simplest form as a reflecting device [START_REF] Feng | Reflectivity characteristics of the fiber loop mirror with a polarization controller[END_REF]. In that case, the input optical signal is split so that two counter-propagating waves are formed in the fiber loop, which then return to recombine (interfere) at the coupler. Since the optical paths for both counter-propagating waves in the loop are the same, the reflectivity of the FLM will be only determined by the coupling ratio of the optical coupler (if the thermal effects and the polarization rotation inside the fiber are disregarded).

In our case, the same effect has been used to tune the ORC's FSR. This has been achieved in a configuration based on a classical FRR architecture including a FLM inside its loop (abbreviated by FRR-FLM). The FLM, of a length LFLM, has been made by connecting the ports on one side of a tunable optical coupler (TC) [see Fig. 5(a)]. A similar architecture to this FRR-FLM has been studied in [START_REF] Vazquez | Tunable optical filters using compound ring resonators for DWDM[END_REF][START_REF] Vazquez | Sagnac loop in ring resonators for tunable optical filters[END_REF], where analytical expressions have been derived for a tunable single-direct-coupled ring resonator. Also, in [START_REF] Seraji | Tunable optical filter based on Sagnac phase-shift using single optical ring resonator[END_REF], the effect of the Sagnac phase-shift by the rotation of the Sagnac loop on the single-direct-coupled ring resonator's response has been studied. In our studies, however, the case of an all-tunable double-directcoupled FRR is considered and investigated in detail. As we will see later this section, this has allowed us to point out many important differences in the characteristics of the second and third outputs of the FRR-FLM. Moreover, the use of a tunable coupler seems to be a robust way to benefit from the Sagnac effect in the FRR-FLM.

In the FRR-FLM that we have designed and studied, we have found that the optical resonance can experience a split leading to the creation of a new ORC. This split value (Δs), or what we can call the new ORC's absolute frequency shift from the main ORC, was depending on the TC's coupling coefficient (κTC). By tuning κTC from 0 to 100 %, we were able to adjust Δs from zero to one complete FSR in the FRR-FLM.

Here it is noteworthy that the FSR for both main and new ORCs in the FRR-FLM is the same. This FSR is given by FSRFRR-FLM = c/n(LFRR+LFLM), where LFRR is the FRR's length excluding the FLM.

Based on our FRR's ADS model, we have developed a model of the FRR-FLM, with an 18m-long FRR and a 2m-long FLM [see Fig. 5(b)]. We have considered a 0.5 dB/km fiber loss (Af) because polarization maintaining fibers have been used in order to get a stable polarization behavior of the optical carrier inside the FRR-FLM [START_REF] Saleh | Optical scattering induced noise in fiber ring resonators and optoelectronic oscillators[END_REF]. The excess losses for the three couplers used were considered as the following: γ1 = γ2 = γTC = 0.03 dB (all these values can be technically achieved). Moreover, an optimal optical coupling has been found for couplers C1 and C2 (κ1= κ1= 8%), allowing to obtain the best achieved transmission characteristics for the FRR-FLM (i.e. low IL, high QOpt, etc…) [START_REF] Saleh | Optical scattering induced noise in fiber ring resonators and optoelectronic oscillators[END_REF]. The optical transmission spectrum at both second and third outputs of the FRR-FLM has been simulated while changing κTC's value from 1 to 100 %. For clarity's sake, in Fig. 6 we only show the simulation results obtained for κTC = 1% and for κTC = 50%. From these results we can notice that the main ORC will remain fixed while the new ORC is uniformly shifted when changing κTC's value. This particularly interesting FRR-FLM's behavior has not been reported in [START_REF] Vazquez | Tunable optical filters using compound ring resonators for DWDM[END_REF][START_REF] Vazquez | Sagnac loop in ring resonators for tunable optical filters[END_REF]. In [START_REF] Vazquez | Tunable optical filters using compound ring resonators for DWDM[END_REF][START_REF] Vazquez | Sagnac loop in ring resonators for tunable optical filters[END_REF], the authors predict that each ORC will experience an equal but opposite frequency shift which also depends on the coupling coefficient κTC but spans a value from 0 to ±FSR/2 when κTC's value varies from 0 to 100 %. Our simulation results show that a frequency shift by FSRFRR-FLM/2 is observed for the new ORC when κTC's value is equal to 50 %. This observed behavior is crucial for many applications if one needs to keep one ORC fixed while tuning the other. On the other hand, if one takes into account the amount of the frequency split reported in [START_REF] Vazquez | Sagnac loop in ring resonators for tunable optical filters[END_REF] in function of κTC's value, the reported results agree with the results obtained using our model: a frequency split that is equal to a complete FSR is obtained when κTC is equal to 100 %.

Aside from that, the results in Fig. 6 show that the combs at the FRR-FLM's second and third outputs seem to have different characteristics. For that reason, we have compared both outputs' characteristics versus κTC's value, particularly regarding the main and new resonances' full width at half maximum (FWHM), IL and worst-case out-of-band rejection (see Fig. 7). The simulation results presented in Fig. 7(a) show that, at both outputs, the FWHM is the same for the main and the new ORC. However, they show that a narrower FWHM can be obtained at the FRR-FLM's third output on a wider range of κTC's variations. This FWHM is in the range of 345 kHz, resulting in a QOpt of 5.6x10 8 at a laser wavelength λl = 1.55 µm. Likewise, the results in Fig. 7(b) show that, for both outputs, the transmission loss is the same for the main and the new ORC. On the other hand, a constant transmission loss can be obtained at the FRR-FLM's second output versus κTC (transmission loss = 8.6 dB), which is not the case at the FRR-FLM's third output. Indeed, the transmission loss at the third output varies from 8 dB to 9.1 dB versus κTc, which is still acceptable if we take into account the behavior of the FWHM and the worst-case out-of-band rejection of the resonances at the third output. Fig. 7(c) results show that, for both outputs, the worst-case out-of-band rejection is the same for the main and the new ORC. However, a higher out-of-band rejection can be obtained at the FRR-FLM's third output on a wider range of κTC's variations.

Experimental results for a FRR-FLM

In order to experimentally confirm the aforementioned theoretical results, we have fabricated a FRR-FLM with some optical devices available at our laboratory. To this purpose, we have firstly characterized a polarization-maintaining variable-ratio optical coupler. This TC's coupling coefficient was adjustable by two different means: a fast tuning, by means of a micrometric screw, and a fine tuning that is remotely controlled using a piezoelectric actuator. The main drawback of this device was its high excess loss [~20% (~1 dB)]. Moreover, the TC's outputs were not well equilibrated and the TC was not hysteresis-free.

Here, it is noteworthy that this FRR-FLM was made for a single test. Therefore, fiber connectors were used to connect the different optical devices in the FRR-FLM instead of splicing them. Moreover, the input and output couplers of the FRR part in the FRR-FLM had the following characteristics: κ1 = κ2 = 30 % (no 8% couplers were available at the laboratory to achieve an optimal coupling in the FRR-FLM) and γ1 = γ2 = 0.1 dB. As a result, the fabricated FRR-FLM had very high transmission loss [see Fig. 5(a)]. This FRR-FLM was 7m-long (LFRR = 5 m and LFLM = 2 m), resulting in a theoretical FSRFRR-FLM of about 29 MHz.

Once fabricated, the 7m-long FRR-FLM has been characterized using a microwave frequency domain characterization technique (see the illustration in Fig. 8). To be able to use this technique, the laser carrier must be first stabilized onto the center of one of the FRR-FLM's optical modes. This is achieved using the Pound-Drever-Hall laser stabilization technique [START_REF] Black | An introduction to Pound-Drever-Hall laser frequency stabilization[END_REF]. Later on, the laser carrier is modulated by means of a frequency-shifter modulator (FSM) driven by a sweeping radio frequency (RF) signal delivered by a microwave vector network analyzer (VNA). As a result, the obtained sweeping single modulation sideband (MSB) will explore the different optical modes of the FRR-FLM. Later on, the beat note of the laser carrier and this optical MSB on a fast photodiode will transcribe the FRR-FLM's transmission characteristics in the RF domain on a given frequency bandwidth. This is done by measuring the S21 transmission parameter on the VNA. More details on this microwave frequency domain characterization technique can be found in [START_REF] Merrer | Characterization technique of optical whispering gallery mode resonators in the microwave frequency domain for optoelectronic oscillators[END_REF]. Fig. 8. Illustration of the microwave frequency domain characterization technique. This illustration shows the importance of using a FSM to modulate the laser carrier in order to generate a single sweeping MSB (rather than two MSBs if a classical linear amplitude modulation of the optical carrier is used) and therefore accurately transcribe the transmission characteristics of a TFRR (or a FRR-FLM) in the RF domain on a given frequency bandwidth In order to correctly characterize the FRR-FLM in the microwave frequency domain, it is essential to use the FSM. In point of fact, the optical resonances belonging to the main and the new ORC in the FRR-FLM are not equidistant regarding each other (see the illustration in Fig. 8). Therefore, if a classical linear amplitude modulation of the optical carrier is used, this will result in a deformation of the RF response. This is because, in that case, two MSBs are generated, that will not overlap at the same time with two optical resonances belonging each to a different ORC of the two ORCs generated by the FRR-FLM. On the contrary, the FSM produces only one optical MSB besides the optical carrier. Consequently, it is possible to stabilize the laser carrier onto one of the optical resonances of the main ORC (ideally) where the single optical MSB will therefore properly explore the other resonances of the two ORCs generated by the FRR-FLM and accurately transcribe them in the microwave frequency domain. Aside from that, the transmission loss was always high, thus preventing a correct description of the FRR-FLM's transmission characteristics. Yet, we were still able to describe the FRR-FLM's absorption characteristics using our microwave characterization bench because we were able to detect more optical power at the absorption port of the FRR-FLM [the FRR-FLM's first output; see Fig. 5(a)]. The S21 transmission parameter measurements are presented in Fig. 9. These results are in very good agreement with the theory as they confirm a total FSRFRR-FLM of 31 MHz and, particularly, a frequency shift in the new ORC by half of the FSRFRR-FLM when κTC's value is around 50%. The observed difference in the absorption resonances' depths is likely due to the non-equilibrated outputs of the variable optical coupler we were using.

Likewise, the 7m-long FRR-FLM has been characterized in the optical domain using a narrow linewidth laser and the wavelength scanning technique [START_REF] Llopis | High-Q optical resonators: characterization and application to stabilization of lasers and high spectral purity microwave oscillators[END_REF]. Because of the high transmission loss issue, an optical amplifier had been added at the FRR-FLM's third output and we were therefore able to measure the transmitted signal while changing κTC's value (κTC = 1%, 30%, 70 % and 87 %). The measurements are depicted in Fig. 10 and they show a qualitatively good agreement with our theoretical predictions on the new ORC's tunability versus κTC in a FRR-FLM. The amplitude deformation in the transmitted signal is most probably due to the non-equilibrated outputs of the variable coupler we were using.

Finally, we can establish that these different theoretical and experimental studies strongly confirm the feasibility of an all-tunable FRR. Obviously, the 7m-long FRR-FLM that we have experimentally characterized was a non-optimized test-resonator used only to validate what our ADS models predicted on the FRR's tunability. Such an all-tunable FRR is certainly very advantageous for many applications in the microwave photonics field (e.g. narrow bandwidth tunable filter, tunable OEO, etc…), particularly if the resonator characteristics can be improved in terms of losses.

All-tunable FRR

The simulation results of section 4 and the experimental results of section 5 are very interesting as they predict, together with section 3's simulation and experimental results, that an all-tunable FRR (TFRR) can be designed if we combine the tunable coupler effect with the phase-shifter effect in the same FRR. This TFRR architecture therefore includes a FLM and a PS as illustrated in Fig. 5(a). Here, it is noteworthy that our model predicts exactly the same TFRR's transmission characteristics either when the PS is inserted inside the FLM's loop or inside the FRR's loop.

The optical transmission spectrum at the third output of a 20m-long TFRR has been therefore simulated using an ADS model. The results are presented in Fig. 11 and they show that a tuning of the new ORC's absolute frequency is obtained by changing the coupling coefficient of the tunable coupler. Also, independently, a uniform tuning of both the new and main ORCs' absolute frequencies is obtained by changing the phase-shift added by the phase-shifter.

Unfortunately such TFRR structure has not been experimentally tested. This was due to the high insertion losses introduced by the PS inside the TFRR. Still, the theoretical and experimental results in this article proves that if one is able to design a TFRR with low loss devices or technology (e.g. a phase shifter and a tunable optical coupler in integrated ring resonators), a complete tunability of the resonator can be obtained as expected in Fig. 11.

Conclusion

In this paper, the optical fiber ring resonator (FRR) tunability has been theoretically studied in detail and we have been able to demonstrate a full tunability of the FRR response. This has been achieved in a particular resonator architecture based on the addition of an optical phase-shifter and a fiber loop mirror inside the FRR's loop. In this architecture, the absolute frequency and the free spectral range of the optical resonance comb generated by the resonator could be both changed independently. This all-tunable FRR has been modeled and a large part of the theoretical results were later on qualitatively validated through different experimental tests. Such an all-tunable FRR can be certainly advantageous for many applications in the microwave photonics field. Our ongoing researches are currently focused on the full experimental characterization of such all-tunable FRR.

Fig. 1 .

 1 Fig. 1. (a) Architecture of the double-direct-coupled FRR (a phase-shifter is inserted to control the phase delay inside the FRR's loop. (b) ADS model of the FRR

Fig. 2 .

 2 Fig. 2. Transmission spectrum at the third output of a 20m-long FRR simulated, using an ADS model, while changing ΔφPS's value by π rad. An absolute frequency shift of the ORC by one half of the FSR is obtained

Fig. 3 .Fig. 4 .

 34 Fig. 3. Experimental setup used to detect the absolute frequency shift in a FRR-PS's ORC. RUT: resonator under test; C3dB: 3dB directional optical coupler; PD: photodiode

Fig. 5 .

 5 Fig. 5. (a) Architecture and (b) ADS model of the FRR including a fiber loop mirror (FRR-FLM). A phase shifter can be inserted into the fiber loop mirror in order to get an all-tunable FRR (TFRR; see section 6)

Fig. 6 .Fig. 7 .

 67 Fig. 6. Optical transmission spectrum at both (a) second and (b) third outputs of the FRR-FLM simulated, using an ADS model, while changing κTC's value. ↑ and↓ notations indicate an increasing or a decreasing κTC's value

Δs

  Photodetector

Fig. 9 .

 9 Fig. 9. Absorption output of a 7m-long FRR-FLM characterized using the microwave frequency domain characterization technique

Fig. 10 .

 10 Fig. 10. Transmitted optical signal measured at the third output of the 7m-long FRR-FLM, using the laser wavelength scanning technique, for different κTC's values: 1%, 30%, 70 % and 87 %. The 2.6 s timescale (or visualization window) on this graphic represents a 96 MHz laser frequency scanning range

Fig. 11 .

 11 Fig.11. Optical transmission spectrum at the third output of a 20m-long TFRR, simulated using an ADS model. A tuning of the new ORC's absolute frequency is obtained by changing the coupling coefficient of the tunable coupler. Independently, a uniform tuning of both the new and main ORCs' absolute frequencies is obtained by changing the phase-shift added by the phase-shifter

Acknowledgements

This work has been supported by the French national research agency (ANR) and the "Direction Générale de l'Armement" (DGA), French MoD procurement agency, under the project ANR-11-ASTR-0029 MINOTOR.