Yushan Wang
email: yushan.wang@lri.fr

Marc Baboulin
email: marc.baboulin@inria.fr

Jack Dongarra
email: dongarra@eecs.utk.edu

Joël Falcou
email: joel.falcou@lri.fr

Yann Fraigneau
email: yann.fraigneau@limsi.fr

Olivier Le Maître

A parallel solver for incompressible fluid flows

Keywords: Navier-Stokes equations, prediction-projection method, ADI method, partial diagonalization, parallel computing, tridiagonal systems, SIMD extension Equations de Navier-Stokes, méthode de prédiction-projection, méthode ADI, diagonalisation partielle, calcul parallèle, systèmes tridiagonaux, extension SIMD

The Navier-Stokes equations describe a large class of fluid flows but are difficult to solve analytically because of their nonlinearity. We present in this paper a parallel solver for the 3-D Navier-Stokes equations of incompressible unsteady flows with constant coefficients, discretized by the finite difference method. We apply a prediction-projection method that transforms the Navier-Stokes equations into three Helmholtz equations and one Poisson equation. For each Helmholtz system, we apply the Alternating Direction Implicit (ADI) method resulting in three tridiagonal systems. The Poisson equation is solved using partial diagonalization which transforms the Laplacian operator into a tridiagonal one. We present an implementation based on MPI where the computations are performed on each subdomain and information is exchanged at the interfaces between subdomains. We describe in particular how the solution of tridiagonal systems can be accelerated using vectorization techniques.

Introduction

The incompressible Navier-Stokes (NS) equations express the Newton's second law applied to fluid motion, with the assumptions that the fluid density is constant and the fluid stress is the sum of a viscous term (proportional to the gradient of the velocity) and an isotropic pressure term. These equations are widely used to model fluid flows in numerous scientific and industrial applications (e.g. in aeronautics, aerospace, meteorology, thermo-hydraulics, and climatology). Efficient numerical methods for the NS equations have then received considerable attention over the last decades, as their simulation remains a challenging problem especially for flows exhibiting complex transient, turbulent dynamics and features at many different scales. Alternative numerical methods for the NS equations include finite difference (FD) [START_REF] Anderson | Computational Fluid Dynamics: The Basics with Applications[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Harlow | Numerical calculation of time-deendent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Joao | Finite difference simulations of the Navier-Stokes equations using parallel distributed computing[END_REF], finite element [START_REF] Pironneau | The Finite Element Method for Fluids[END_REF][START_REF] Rannacher | Finite element methods for the incompressible Navier-Stokes equations[END_REF], finite volume [START_REF] Anderson | Computational Fluid Dynamics: The Basics with Applications[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF], spectral methods [START_REF] Canuto | Spectral Methods: Fundamentals in Single Domains[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF] and mesh-free methods such as the vortex method [START_REF] Cottet | Vortex Methods Theory and Practice[END_REF] or smoothed particle hydrodynamics [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Li | Meshfree and particle methods and their applications[END_REF][START_REF] Martin | An incompressible Navier-Stokes with particles algorithm and parallel implementation[END_REF]. These methods have pros and cons, depending on the problem considered. In addition, the evolution of computer architectures and the emergence of massively parallel computers call for the development of new algorithms to fully exploit hardware advances.

In this paper, we focus on FD discretizations (of second order in time and space) of the threedimensional NS equations which are solved using the classical prediction-projection method [START_REF] Alexandre | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Naver-Stokes equations[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Brown | Accurate projection methods for the incompressible Navier-Stokes equations[END_REF]. This method relies on an explicit treatment of the non-linear terms and an implicit treatment of the viscous and pressure terms to ensure stability. As detailed in Section 2, it leads to the solution of 3-D diffusion (velocity) and Poisson (pressure) equations that constitute the main computational workload of the simulation. Further, we restrict ourself to spatial discretizations on structured cartesian grids resulting from the tensorization of one-dimensional grid. Although not general, the structured character of these grids greatly facilitates the design and efficient implementation of the algorithms on modern computer architectures, and encompasses many situations of interest. Even for these simple meshes, 3-D parallel solvers are needed since complex flow simulations and targeted applications commonly require grid resolution corresponding to 10 7 -10 9 unknowns to properly capture the flow dynamics. These large simulations require the use of parallel computers. Another important advantage brought by the predictionprojection approach on tensored grids is that it leads to Helmholtz-like or Poisson equations that can be consistently reduced to one-dimensional operators using respectively operator splitting and partial diagonalization [START_REF] Abdi | The eigen-decomposition: Eigenvalues and eigenvectors[END_REF]. The second order discretization results in tridiagonal systems that can solved efficiently by taking advantage of current parallel architectures. The parallel strategy is based on the idea of dividing the computational work between several processors which perform their calculations concurrently. We use a domain decomposition method in which the computational domain is divided into subdomains, each associated to one processor. The Poisson and Helmholtz-like equations are then solved over each subdomain in parallel using the Schur complement method to ensure the continuity of the solution at the subdomains interfaces. This involves data communication between neighboring subdomains for the interface solution. Communications are handled by the Message-Passing-Interface (MPI) [START_REF]MPI : A Message-Passing Interface Standard[END_REF]. We also rely on the BLAS [START_REF]Basic Linear Algebra Subprograms Technical Forum Standard[END_REF], LAPACK [START_REF] Anderson | LAPACK Users' Guide[END_REF] and ScaLAPACK [START_REF] Blackford | ScaLA-PACK Users' Guide[END_REF] kernels from the MKL [START_REF]Intel. Math Kernel Library (MKL)[END_REF] library to perform the computations. In addition to the combination of the numerical methods mentioned above, an important contribution of this work is the implementation of a vectorized version of Thomas algorithm for a fast solution of the tridiagonal systems.

This paper is organized as follows. Section 2 summarizes the solution of the Navier-Stokes equations for incompressible fluid flows and the prediction-projection method. Then, the FD discretization in space and time is introduced together with the directional splitting and partial diagonalization to obtain the tridiagonal systems to be solved. Section 3 discusses the parallel implementation and emphasizes on vectorization techniques to improve the Thomas algorithm

RR n°8247

for the solution of the tridiagonal systems. In Section 4, several performance measurements of our Navier-Stokes solver are provided using a current parallel system. Concluding remarks are presented in Section 5.

Navier-Stokes Problem

We consider the unsteady flow of a incompressible Newtonian fluid governed by the Navier-Stokes (NS) equations that can be expressed as

∂V ∂t + ∇ • (V T ⊗ V ) = -∇P + 1 Re ∆V , (1)
∇ • V = 0, (2)
where V = (u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)) T is the velocity vector and P = P (x, y, z, t) is the pressure. Eqs. ( 1) and ( 2) state the conservation of momentum and mass respectively. We have denoted ∇ the nabla and ∆ the Laplacian operators, while ⊗ is the Kronecker product [START_REF] Graham | Kronecker Products and Matrix Calculus with Application[END_REF].

In the following, we use the shorthand notation CT := ∇ • (V T ⊗ V ) for the convective term. In Eq. ( 1), Re := ρU L/µ is the Reynolds number based on a reference length L, a reference velocity U , the fluid density ρ, and the dynamic viscosity µ. The Reynolds number is a good indicator of the fluid state (laminar, transient, turbulent). We see from Eq. ( 1) that the larger Re is, the more dominant the convection is, with respect to viscous effects. This leads to the onset of flow instabilities and a turbulent state of the flow for high Reynolds numbers. In turbulent flows, the order of magnitude for the ratio between the largest and the smallest length scale of eddy structures is Re . So, a finer discretization is required when the Reynolds number increases and the resulting order of magnitude for the 3D mesh size will be Re 9 4 . Solving the NS equations also involves defining the initial and boundary conditions. The details about the implementations of the boundary condition will not be discussed in this paper. We mention that our code can deal with mixed Dirichlet and Neumann (homogeneous or not), and periodic/symmetric boundary conditions. The initial conditions for velocity and pressure will depend on the physical problem.

Prediction-Projection method

In this section, we present the prediction-projection method applied on the NS equations. The spatial discretization uses second-order finite differences approximations of the differential operator. The domain is discretized using N x × N y × N z cells, in the x, y and z directions. The discrete velocity and pressure unknowns are defined along a staggered arrangement, in which the pressure is defined on the center of each cell, the velocity components are defined on the center of the corresponding edge, see Fig. 1, to prevent spurious oscillations [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF].

We denote the time level with upper-indices and ∆t the time-step size. We rely on a second order approximation of the time derivative (BDF2) [START_REF] Ascher | Computational Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF][START_REF] Hundsdorfer | Partially implicit BDF2 blends for convection dominated flows[END_REF], an implicit treatment of the viscous and pressure term, and a second order explicit extrapolation of the convective term.

Then the discrete NS equations can be expressed as where

3V n+1 -4V n + V n-1 2∆t + CT n+1 = -∇P n+1 + 1 Re ∆V n+1 , (3)
∇ • V n+1 = 0, (4)
Inria Ú , •½ Ù , •½ Ú , Ú ¹½, Ú , ¹½ Ú •½, Ù •½, Ù , Ù ¹½, Ù , ¹½ Ó v ÓÒØÖÓÐ ÚÓÐÙÑ Ó u Ñ × Ð Ñ ÒØ ÚÓÐÙÑ ÓÒØÖÓÐ Ü , •½ Ü •½, •½ Ü ¹½, •½ Ü ¹½, Ü , Ü •½, Ü ¹½, ¹½ Ü , ¹½ Ü •½, ¹½ È •½, •½ È , •½ È ¹½, •½ È •½, È •½, ¹½ È , È , ¹½ È ¹½, È ¹½, ¹½
CT n+1 = CT n -CT n-1 .
The integration for one time-step consists of two steps: a prediction followed by a projection. The latter step is also called the correction step.

Prediction step

In the prediction step, an approximate solution for the momentum equation is sought without accounting for the incompressibility condition in Eq. ( 4). More specifically, we compute an auxiliary velocity

V * 3V * -4V n + V n-1 2∆t + CT n+1 = -∇P n + 1 Re ∆V * . (5)
Note that in Eq. ( 5), the pressure P n+1 has been substituted by P n . Defining σ := 2∆t/3, the Eq. ( 5) for V * can be recast as

I - σ Re ∆ V * = S, S := 4V n -V n-1 3 -σ( CT n+1 + ∇P n ). (6)
To solve Eq. ( 6), we apply the Alternating Direction Implicit (ADI) method which approximates the 3-D operator (I -σ Re ∆) as a product of three 1-D operators:

I - σ Re ∆ = I - σ Re ∆ x I - σ Re ∆ y I - σ Re ∆ z + O(∆t 2 ),
where ∆ i is the second order derivative along the i-direction. With this product, we can transform the 3-D Helmholtz-like equation (Eq. ( 6)) into a system of three 1-D Helmholtz-like equations

I - σ Re ∆ x T 1 = S, I - σ Re ∆ y T 2 = T 1 , I - σ Re ∆ z V * = T 2 . (7)
In the remainder of this paper, Eq. ( 7) will be simply referred to as Helmholtz equation.

RR n°8247

Projection step

In general, V * is not divergence-free so it is corrected in the second step. The correction is sought as a potential velocity field as follows. Taking the divergence of the difference of Eqs. ( 3) and ( 5), it comes:

∇ • (V n+1 -V * ) = -σ∆(P n+1 -P n ) + σ Re ∇ • ∆(V n+1 -V * ). (8)
Since ∇ • V n+1 = 0 and since, for any vector U , we have ∇ • ∆U = ∆∇ • U , then Eq. ( 8) becomes

∇ • V * = σ∆ P n+1 -P n - 1 Re ∇ • V * . ( 9
)
Denoting φ := P n+1 -

P n + 1 Re ∇ • V * , then φ satisfies the Poisson equation ∆φ = 1 σ ∇ • V * , (10)
which is complemented with appropriate boundary conditions. After computing φ, we update the velocity and pressure fields using

V n+1 = V * -σφ, P n+1 = P n + φ - 1 Re ∇ • V * .
There exist several methods to solve the Poisson equation given in Eq. [START_REF] Brown | Accurate projection methods for the incompressible Navier-Stokes equations[END_REF]. For example, the fast Fourier transformation (FFT), fast sine/cosine transformations, and multigrid methods [START_REF] Noury | A fast Poisson solver for OpenCL using multigrid methods[END_REF]. In this paper, we consider a direct method based on partial diagonalization.

Let us denote ∆φ = 1 σ ∇ • V * = S, and rewrite ∆ = ∆ x + ∆ y + ∆ z . Then the second order central schemes for these 1-D operators ∆ x , ∆ y and ∆ z give three tridiagonal matrices for the considered direction. To further simplify the computation, we transform tridiagonal matrices into diagonal ones by using their eigendecomposition [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] .

Given a real square (N × N ) matrix A with N linearly independent eigenvectors q i (i = 1, ..., N ), then A can be factorized as A = QΛQ -1 where Λ is the diagonal matrix whose diagonal elements are the corresponding eigenvalues, and Q is the square (N ×N ) matrix whose ith column is the eigenvector q i of A.

For our problem, we decompose only ∆ x , ∆ y and keep ∆ z , then we have

(Q x Λ x Q -1 x + Q y Λ y Q -1 y + ∆ z )φ = S. (11)
Then we multiply Eq. ( 11) with Q -1 x Q -1 y , and observe that

Q -1 x Q -1 y = Q -1 y Q -1 x , Q -1 x Λ -1 y = Λ -1 y Q -1 x , Q -1 y Λ -1 x = Λ -1 x Q -1 y ,
so we have

(Λ x + Λ y + ∆ z ) φ = S, where φ = Q -1 y Q -1 x φ and S = Q -1 y Q -1 x S. (12)
So the original Poisson equation (Eq. 10) is now transformed into the form of Eq. ( 12) which is a 1-D tridiagonal equation for the z-direction. In our parallelization strategy we divide the domain equally into subdomains and we assign each subdomain to one MPI process in such a way that each process has the same amount of computational workload. Once the subdomains and the boundary conditions have been defined, the computations on each subdomain can be performed in parallel. The boundaries here include the real domain boundaries and the interfaces between the subdomains. The classical Schur complement method [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] is used to ensure continuity of the solution at the interface between subdomains, and the information at the interfaces is exchanged via MPI routines. Fig. 2 shows an example of a 2-D domain decomposition with boundary cells. For instance in subdomain IV, the top and right boundaries are real domain boundaries while the left and lower boundaries are actually interfaces with neighboring subdomains. We suppose now that we have obtained the tridiagonal system for one direction (corresponding for instance to the x-direction, as depicted in Fig. 2). The matrix structure is shown by Fig. 3(a). This tridiagonal system is solved by the Schur complement method. By reordering the variables so that the interface variables are listed last, we obtain a new system where the matrix has a structure as depicted in Fig. 3(b).

We recall below how the solution is computed via the Schur complement method. Given a Wang et al.

block system

A k A ki A ik A i X k X i = S k S i , (13)
where indices k and i refer respectively to interior and interface variables.

From the first block equation in Eq. ( 13) the unknown vector X k can be expressed as

X k = A -1 k (S k -A ki X i ), (14)
and we substitute this into the second equation to obtain

(A i -A ik A -1 k A ki )X i = S i -A ik A -1 k S k . ( 15
)
If the Schur complement matrix A i -A ik A -1 k A ki can be formed and the linear system (15) can be solved, all the interface variables X i are then available and the remaining unknowns can be computed via Eq. ( 14). Due to the particular block structure of A k and A i , we have r separate systems that can be solved in parallel, where r is the number of subdomains along the direction in which we solve the original tridiagonal system.

Solving tridiagonal systems

We describe in this section how the block-tridiagonal systems resulting from the discretizations of Eq. ( 7) (using ADI method) and Eq. ( 10) (using partial diagonalization) can be solved efficiently. There have been several efforts in recent years to develop efficient tridiagonal solvers [START_REF] Bini | The cyclic reduction algorithm: From Poisson equation to stochastic processes and beyond[END_REF][START_REF] George | Parallel Scientific Computing in C++ and MPI[END_REF][START_REF] Polizzi | A parallel hybrid banded system solver: the SPIKE algorithm[END_REF]. A well-known procedure is based on the Thomas algorithm [START_REF] Atkinson | Elementary Numerical Analysis[END_REF] which basically corresponds to a Gaussian elimination without pivoting. Given a tridiagonal system Ax = s, where A ∈ R n×n is diagonally dominant (nonsingular), and x, s ∈ R n , the Thomas algorithm consists of two steps. The first step is a forward elimination where we eliminate the lower diagonal coefficients and transform the original matrix into an upper triangular one. The second step is a backward substitution that solves the upper triangular system. This requires O(n) operations.

To better understand the structure of the linear system, let us take a simple 2-D example as depicted in Fig. 4 which uses a 4 × 3 2D mesh and let us consider the operator ∆ x in Eq. [START_REF] Atkinson | Elementary Numerical Analysis[END_REF]. We order the unknowns first by row and then by column as shown in Fig. 4. Then the matrix corresponding to the operator ∆ x can be represented as the matrix given in Fig. 5 (with m = 3 tridiagonal blocks). Similarly, if we consider the operator ∆ y in Eq. ( 7), the unknowns will be ordered first by column then by row and the resulting matrix will contain 4 tridiagonal blocks.

We apply the Thomas algorithm within each tridiagonal block. If the tridiagonal blocks are identical, then we consider a smaller system with multiple right-hand sides (RHS) that can be expressed as AX = S with A ∈ R n×n and X, S ∈ R n×m . In the following we explain how the solution of this system can be accelerated using vectorization techniques.

Since the late 90's, processor manufacturers provide specialized processing units called multimedia extensions or Single Instruction Multiple Data (SIMD) extensions. This feature allows processors to exploit the latent data parallelism available in applications by executing a given instruction simultaneously on multiple data stored in a single special register. However, taking advantage of SIMD extensions in applications remains a complex task. We can find for instance in [START_REF] Esterie | Boost.simd: generic programming for portable simdization[END_REF] a description of Boost.SIMD, a high-level C++ library to program SIMD architectures. Boost.SIMD provides both expressiveness and performance by using generic programming to handle vectorization in a portable way. A first requirement to achieve performance using SIMD is to ensure memory alignment. In our case, this means that the Fortran arrays have to be allocated by a dedicated C++ function that uses the ISO_C_BINDING interface.

                    1
                   
The SIMD opportunity in our algorithm stands from the fact that solving the tridiagonal systems with multiple RHS can be performed in parallel. This means that once we have an SIMD vector containing the tridiagonal values of a set of rows, the backward substitution of these rows can be issued as a single SIMD instruction. The triangularization followed by the backward substitution will be accelerated by the fact that data handled at each step of the algorithm will be contiguous in memory using a shuffling operation. Shuffling is a typical idiom of SIMD programming that replaces some class of complex memory access patterns by computation on simpler patterns. Depending on the SIMD extension, these shuffling operations can either be limited (like in SSE2 where shuffling can only occur piecewise inside a given register) or be arbitrary (e.g. in Altivec where shuffling corresponds to permutations of bytes). Shuffling also enables deinterleaving of scattered data that can be fetched from main memory and bring back into a single, contiguous SIMD register. In our case, we use shuffling to aggregate values from the sparse representation of the system into a set of contiguous values located at the proper place so that the elimination and backward substitution phases can take place inside an SIMD vector. Expected benefit should be a faster code due to less memory access in different part of the memory (thus limiting cache misses) and augmented SIMD speedup due to the layout of the data after the shuffling and better use of pipelining.

We illustrate this method in Fig. 6 by considering 2 RHS (case of SSE where a double precision vector holds 2 values) where we would perform the following operations:

load 3 registers units with coefficients (a i , a i+1 ), (b i , b i+1 ) and (c i , c i+1 ) coming from respectively the lower, main and upper diagonals,

-shuffle the data to have a 2 × 3 local matrix a i b i c i a i+1 b i+1 ,
perform two Gaussian eliminations to eliminate the coefficients (a i , a i+1 ) in the lower diagonal, -after elimination of the a i 's, we perform a backward substitution to solve the upper triangular

RR n°8247

system with the RHS coefficients loaded in one register (e.g. 2 RHS for SSE or 4 RHS for AVX SIMD extension). Our implementation uses this method to handle multiple RHS. Note that the memory alignment mentioned previously enables here the "load" process to be performed efficiently. We point out that this vectorization approach can be also used in other parts of our solver, for instance to compute the source term faster.

= c3 a3 a4 c1 c2 x 1 1 x 2 1 x 1 2 x 2 2 x 1 3 x 2 3 x 1 4 x 2 4 s 1 1 s 1 2 s 2 1 s 2 2 b1 b2 b3 b4 s 1 3 s 1 4 s 2 3 s 2 4 (a) Current iteration. = c3 c1 c2 x 1 1 x 2 1 x 1 2 x 2 2 x 1 3 x 2 3 x 1 4 x 2 4 s 1 1 s 1 2 s 2 1 s 2 2 b1 b2 b3 b4 s 1 3 s 1 4 s 2 3 s 2 4 (b) After elimination. = c3 c1 c2 s 1 1 s 1 2 s 2 1 s 2 2 b1 b2 b3 b4 s 1 3 s 1 4 s 2 3 s 2 4 x 1 4 x 2 4 x 2 3 x 1 3 x 2 2 x 1 2 x 2 1 x 1 1 (c) Backward substitution.

Performance results

The following experiments were carried out using a MagnyCours-48 system from University of Tennessee. This machine has a NUMA architecture and is composed of four AMD Opteron 6172 (with a SSE-4a instruction set) running at 2.1GHz with twelve cores each (48 cores total) and 128GB of memory. Our solver is linked with the LAPACK and ScaLAPACK routines from the 10.3.6 version of the Intel MKL [START_REF]Intel. Math Kernel Library (MKL)[END_REF] library and communications are performed using OpenMPI 1.4.3. In our experiments, we use one MPI process per core and each MPI process is a single thread.

In Fig. 7(a) we consider a 3-D vortex problem with mesh size 240 3 , i.e. about 1.4 × 10 7 unknowns. We represent the performance (in seconds) for one iteration of the NS solver including the Helmholtz and Poisson equations, and miscellaneous tasks (mainly I/O and velocity/pressure updates). The number of threads varies from 1 to 48. We observe in Fig. 7(a) that the CPU time decreases significantly with the number of threads, showing then a good scalability of the solver. We observe that the Helmholtz equation represents about 30% of the global computational time and this percentage remains the same when the number of threads increases. shows the parallel speedup of the solver. We obtain a speedup up of 33 using 48 cores. We also observe in Fig. 7(b) that the solution for the Poisson equation scales slightly better than for the Helmholtz equation. This is because in the Poisson equation we have only one tridiagonal system to solve (vs three systems in Helmholtz equation) and thus there is less information exchanged (as mentionned in Section 3, when we solve a tridiagonal system via the Schur complement method, we have a reordering of the variables that results in additionnal communication).

In a second experiment, we consider a 3-D vortex problem where the size increases with the number of threads and with a fixed mesh size per subdomain (weak scalability). The local size is set to 240 × 240 × 10 and 10 time iterations are performed. We observe in Fig. 8 that the time spent in solving Helmholtz and Poisson equations does not vary (about 60 and 110 milliseconds respectively) because the number of unknowns computed is always 240 2 × 10 = 5.76 × 10 5 . However the global CPU time increases linearly with the number of processes due to a larger amount of I/O after each iteration. Indeed, the numerical solution is stored after each time iteration in order to generate a detailed animation that visualizes the fluid movement. In this numerical test, we only divide the domain along one direction (z-direction here) to keep the shape of the interfaces. Thus, the amount of information to be exchanged (the number of interface elements) is 240 2 × (number of MPI processes -1).

Let us now study more specifically the performance of the tridiagonal solver described in Section 3.2. We compare the performance of a vectorized Thomas algorithm (as illustrated by Fig. 6) with the LAPACK 3.2 routine DGTSV from Netlib, linked with the MKL BLAS library, which solves a general tridiagonal system using Gaussian elimination with partial pivoting (note that, since our matrices are diagonally dominant, the routine DGTSV does not pivot and only the search for pivot is performed). In Fig. 9(a), we plot the number of cycles per value, i.e. the amount of CPU cycles required to compute one element of the result (we operate here on tridiagonal blocks of size 100 in order to exhibit the cache effects). This metric allows us to do a fine grain analysis of both the impact of vectorization and the impact of memory access. For the curve related to DGTSV, we observe a cycle/value amount that jumps when we hit the cache size (L2 and L3). This is due to the fact that every computation in the elimination and backward substitution processes requires more memory access to be completed. However, when we consider the vectorized Thomas algorithm, we have a constant amount of CPU cycle/value up to the largest sizes, due to the shuffling that replaces memory access by computation. Fig. 9(b) shows the time ratio between DGTSV and the vectorized Thomas algorithm. For large problem sizes (e.g. more than 10 6 unknowns), we observe that the vectorized version of Thomas algorithm outperforms DGTSV with a factor 3.

Conclusion

In this paper, we presented a numerical method and a parallel implementation for the solution of the 3-D Navier-Stokes equations. The resulting solver is scalable on problems of realistic size and the solution of the tridiagonal systems can be significantly enhanced by using vectorization. We point out that this solver can also take into account other characteristics of the fluid or of the domain that are not discussed in details in the paper. For example, several boundary conditions are implemented and the domain can be treated either as continuous or discontinuous. Note that other physical factors like temperature and density are also considered in the solver.

Figure 1 :

 1 Figure 1: Staggered mesh structure.

Inria 3

 3 Parallel implementation for solving Navier-Stokes equations 3.1 Domain decomposition and Schur complement method

Figure 2 :

 2 Figure 2: Domain decomposition and boundary cells.

Figure 3 :

 3 Figure 3: Illustration of variable reordering for Schur Complement method for the case of 2 subdomains. ( : interior variables. : interface variables.)

Inria

Figure 4 :

 4 Figure 4: 2D n×m mesh with n = 4, m = 3.

Figure 5 :

 5 Figure 5: An example of block tridiagonal matrix for operator ∆ x .

Figure 6 :

 6 Figure 6: Thomas algorithm with vectorization.

  Time breakdown for one iteration of the NS solver.

  Scalability of NS solver.

Figure 7 :Fig. 7 (

 77 Figure 7: Performance of the Navier-Stokes solver.

Figure 8 :

 8 Figure 8: Scalability of problems with variable size.

  Cycles per value for Thomas algorithm.

  Time ratio between DGTSV and vectorized Thomas.

Figure 9 :

 9 Figure 9: Performance of vectorized Thomas algorithm.

RR n°8247

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

This work was supported by Région Île-de-France and Digitéo (http://www.digiteo.fr), CAL-IFHA project, contract No 2011-038D. We also thank Pierre Esterie (University Paris-Sud) for his support in using Boost SIMD.

Inria


Lina Ye
email: lina.ye@lri.fr

Philippe Dague
email: philippe.dague@lri.fr

New Results for Joint Diagnosability of Self-observed Distributed Discrete Event Systems

Diagnosability is an important property that determines at design stage how accurate any diagnosis algorithm can be on a partially observable system. Most existing approaches assumed that each observable event in the system is globally observed. Considering the cases where there is no global information, a recent work has proposed a new framework to check diagnosability in a system where each component can only observe its own observable events to keep the internal structure private in terms of observations. However, the authors implicitly assume that the local paths in each component can be exhaustively enumerated, which is not true in a general case where there are embedded cycles. In this paper, we get some new results about diagnosability in such a system, i.e., what we call joint diagnosability in a self-observed distributed system. First we prove its undecidability with unobservable communication events by reducing the Post's Correspondence Problem (PCP) to an observation problem, inspired from an existing work. Then we propose an algorithm to check a sufficient but not necessary condition of joint diagnosability. Finally we briefly discuss about the decidable case where communication events are all observable.

INTRODUCTION

Over the latest decades, with more performance requirements imposed on the complex systems, they are subject to more errors. However, it is not realistic to detect faults manually for such complex systems. Automated diagnosis mechanisms are therefore required to monitor large distributed applications. Generally speaking, diagnosis reasoning aims at detecting possible faults that can explain the observations. The possibility to achieve such a diagnosis reasoning depends on the diagnosability of the system. Diagnosability is an important property that determines at design stage how accurate any diagnosis algorithm can be on a partially observable system. The diagnosability analysis problem has received considerable attention in the literature. In this paper, the systems we discuss about are Discrete Event Systems (DES).

Some existing works analyzed diagnosability in a centralized way [START_REF] Sampath | Diagnosability of discrete event systems[END_REF], [START_REF] Jiang | A polynomial time algorithm for diagnosability of discrete event systems[END_REF], etc.), i.e., a monolithic model of a given system is hypothesized, which is a very powerful information for diagnosability analysis. But such a global model is not always available and moreover the centralized approaches lead to combinatorial explosion of the search space. This is why very recently distributed approaches for diagnosability began to be investigated ( (Pencolé, 2004), [START_REF] Schumann | Scalable diagnosability checking of event-driven systems[END_REF], etc.), relying on local objects. However, all these approaches assume that each observable event in the system can be observed by all components, i.e., globally observed. However, there are some cases where it is not possible to assume the presence of global information. For example, networked control systems are characterized by the fact that multiple distributed components possess their own part of available information instead of a global knowledge. Then [START_REF] Ye | Diagnosability analysis of discrete event systems with autonomous components[END_REF] has proposed a new framework to check diagnosability in a system where each component can only observe its own observable events to keep the internal structure private in terms of observations. However, the authors implicitly assume that the local paths in each component can be exhaustively enumerated, which is not true in a general case where there are embedded cycles. In this paper, we generalize this work to get some new results about the diagnosability of what we call self-observed distributed systems, i.e., systems where locally observable events can only be observed by their own component.

We make several contributions in this paper. The first one is to extend diagnosability of globally observed systems to what we call joint diagnosability of self-observed systems and then to prove its undecidability, in the general case where communication events are unobservable, by reducing the Post's Correspondence Problem (PCP) to a joint observability problem, which is inspired from the undecidability result of joint observability (Tripakis, 2001). The second one is to propose an algorithm for testing a sufficient condition of joint diagnosability with unobservability of communication events, where we first obtain pairs of local trajectories in the faulty component, such that for each pair only one trajectory contains the fault but both trajectories have the same local observations, and then check their global consistency into two phases.

We provide the proof that it is a sufficient condition and point out why it is not necessary. The third one is to discuss about the decidable case where communication events are observable.

PRELIMINARIES

In this section, we model self-observed distributed DES and then recall joint diagnosability features [START_REF] Ye | Diagnosability analysis of discrete event systems with autonomous components[END_REF]. We consider a self-observed distributed DES composed of a set of components {G 1 , G 2 ,..., G n } that communicate by communication events. Each component can only observe its own observable events and thus keeps its internal structure private in terms of observations. Such a system is modeled by a set of finite state machines (FSM), each one representing the local model of one component. The local model of a component G i is a FSM, denoted by

G i = (Q i , Σ i , δ i , q 0 i ), where Q i is the set of states; Σ i is the set of events; δ i ⊆ Q i × Σ i × Q i is
the set of transitions; and q 0 i is the initial state. The set of events Σ i is partitioned into four subsets: Σ io , the set of locally observable events that can be observed only by their own component G i ; Σ iu , the set of unobservable normal events; Σ i f , the set of unobservable fault events; and Σ ic , the set of communication events shared by at least one other component, which are the only shared events between components. Figure 1 depicts a self-observed distributed system with two components: G 1 (left) and G 2 (right), where the events Oi denote locally observable events, the event F denotes an unobservable fault event, the events U i denote unobservable normal events and the events Ci denote communication events.

We denote the synchronized FSM of components G 1 , ..., G n by (G 1 , ..., G n ), where the synchronized events are the shared events between components and any one of them always occurs simultaneously in all components that define it. The state space of the synchronized FSM is the Cartesian product of the state spaces of components. The global model of the entire system is implicitly defined as the synchronized FSM of all components based on their shared events, i.e., communication events. However, the global model will not be calculated in this paper since in a selfobserved distributed system, the global occurrence order of observable events is not accessible. In the following, we call subsystem of G the synchronization of a subset of components of G, i.e., (G s1 , ..., G sm ), where {s 1 , ...s m } ⊆ {1, ...n}. Note that one component or the entire system can also be considered as a subsystem.

Given the system model G = (Q, Σ, δ, q 0 ), the set of words produced by the FSM G is a prefix-closed language L(G) that describes the normal and faulty behaviors of the system. Formally, L(G) = {s ∈ Σ * |∃q ∈ Q, (q 0 , s, q) ∈ δ}, where the transition δ has been extended from events to words. In the following, we call a word of L(G) also a trajectory in the system G and a sequence q 0 σ 0 q 1 σ 1 ... a path in G, where σ 0 σ 1 ... is a trajectory and for all i, we have (q i , σ i , q i+1 ) ∈ δ. Given s ∈ L(G), we denote the post-language of L(G) after s by L(G)/s and denote the projection of s to observable events of G (resp. G i ) by P (s) (resp. P i (s)). For example, if s = O1.U 2.O3 * , then we have P (s) = O1.O3 * , where Oi denotes an observable event. We adopt the assumption described in (Pencolé, 2004), i.e., the projection of the global language on each local model is observable live, in particular there is no unobservable cycle in any component. For the sake of simplicity, our approach is illustrated by dealing with only one fault, which can be extended to the case with multiple faults.

Next we rephrase for our context the definition of reconstructibility introduced in (Cori and [START_REF] Dague | Recognizable subsets of some partially abelian monoids[END_REF]. Definition 1 (Reconstructibility). Given a system G that is composed of several subsystems, i.e., G = (G s1 , ..., G sm ), a set of trajectories in these subsystems is said to be reconstructible with respect to G if it is obtained by projection on this set of subsystems of a trajectory ρ in G.

If there is no common communication event between two subsystems, then any trajectory in one subsystem and any one in the other subsystem are reconstructible.

For the sake of consistency, now we rename what is called cooperative diagnosability in [START_REF] Ye | Diagnosability analysis of discrete event systems with autonomous components[END_REF] as joint diagnosability. We denote a trajectory ending with the fault f by s f . Definition 2 (Joint diagnosability). A fault f is jointly diagnosable in a self-observed distributed system G with components {G

1 , ...G n }, iff ∃k ∈ N, ∀s f ∈ L(G), ∀t ∈ L(G)/s f , (∀i ∈ {1, ..., n}, |P i (t)| ≥ k)⇒ (∀p ∈ L(G) (∀i ∈ {1, ..., n}, P i (p) = P i (s f .t)) ⇒ f ∈ p).
Joint diagnosability of a fault f means that for each trajectory s f in G, for each t that is an extension of s f with enough locally observable events in all components, every trajectory p in G that is equivalent to s f .t for local observations in each component should contain in it f . In a self-observed system, we call a pair of trajectories p and p′ satisfying the three following conditions a (global) indeterminate pair: 1) p contains f and p′ does not; 2) p has arbitrarily long local observations in all components after the occurrence of f ; 3) ∀i ∈ {1, ..., n}, P i (p) = P i (p′). Here arbitrarily long local observations can be considered as infinite local observations. Now we have the following theorem [START_REF] Ye | Diagnosability analysis of discrete event systems with autonomous components[END_REF]. Theorem 1 Given a self-observed distributed system G, a fault f is jointly diagnosable in G iff there is no (global) indeterminate pair in G.

2 3 UNDECIDABLE CASE To discuss about joint diagnosability, we consider two cases: 1) communication events are unobservable; 2) communication events are observable. In this and the next sections, we consider the first general case, i.e., unobservability of communication events.

From theorem 1, we know that checking joint diagnosability boils down to check the existence of (global) indeterminate pairs that witness non joint diagnosability. Inspired from (Tripakis, 2001), we discuss first about whether it is decidable or not. To be selfcontained, we rephrase joint observability (Tripakis, 2001), where the global model of a system is assumed. Definition 3 (Joint observability) Given the global model G of a system and the sets of observable events Σ io ⊆ Σ, i = 1, ..., k, then a fault f is jointly observable in L(G) w.r.t. these sets if the following condition holds, where the projection of s to Σ io is denoted by

P Σi o (s): ∀s ∈ L(G), s′ ∈ L(G), f ∈ s, f / ∈ s′, ∃i = 1, ...k, P Σi o (s) = P Σi o (s′)
Joint observability of a fault in a system means that there are not two system behaviors such that only one of them contains the fault but their projections on each set of observable events are the same. Then the undecidability of joint observability with at least two sets of observable events is proved by reducing the Post's Correspondence Problem (PCP) to an observation problem. Now we briefly describe the outline of the proof with two sets (Tripakis, 2001). 1)PCP: given a finite alphabet Σ, two sets of words v 1 , v 2 , ..., v k and z 1 , z 2 , ..., z k over Σ, then a solution to PCP is a sequence of indices (i m ) 1≤m≤n with n ≥ 1 and 1 ≤ i m ≤ k for all m such that v i1 v i2 ...v in = z i1 z i2 ...z in .

2) Now let Σ′ = {a 1 , ...a k } be a set of new letters, not in Σ. Then consider the language L over Σ ∪ Σ′ ∪ {good, bad}, defined by the regular expression:

good(v 1 a 1 + ... + v k a k ) + + bad(z 1 a 1 + ... + z k a k ) + ,
where Σ + denotes the set of all finite words over Σ except ǫ. All words in L that start with good constitute the normal behaviors and all those that start with bad constitute the faulty behaviors.

3) If there is a solution for the above PCP, i.e., there exist indices i

1 , ..i n ∈ {1...k}, n ≥ 1, such that v i1 v i2 ...v in = z i1 z i2 ...z in ,
then the fault is not jointly observable w.r.t. observable events sets Σ and Σ′. The reason is that in this case, we have a pair of words ρ = goodv i1 a i1 v i2 a i2 ...v in a in and ρ′ = badz i1 a i1 z i2 a i2 ...z in a in such that both ρ, ρ′ ∈ L with ρ′ only containing the fault, and ρ, ρ′ have the same observations both for Σ and Σ′, which violates joint observability. 4) On the other side, if the fault is not jointly observable, there is at least one pair of words violating joint observability, denoted by ρ and ρ′. Since only one of them is a normal behavior, say ρ, then ρ must be of the form goodv i1 a i1 v i2 a i2 ...v in a in and ρ′ must be of the form badz j1 a j1 z j2 a j2 ...z j l a j l . Furthermore, we know that ρ and ρ′ have the same observations both for Σ and Σ′. So we have a i1 a i2 ...a in = a j1 a j2 ...a j l , which means that l = n, i 1 = j 1 , i 2 = j 2 , ..., i n = j n since a im is a letter not a word. And then we also get v i1 v i2 ...v in = z i1 z i2 ...z in , which means that there does exist a solution for the above PCP. 5) Thus checking joint observability, i.e., checking the existence of such a pair violating joint observability, is equivalent to check for a solution to PCP. Since PCP is a well known undecidable problem, then checking joint observability is also undecidable. Now, from definition 2, we know that joint diagnosability is violated iff it exists two infinite trajectories, one with the fault and the other without the fault, and both give the same infinite sequence of observations for each component (observer), i.e. it exists an indeterminate pair. So, for just two components, the example of the proof above shows that joint diagnosability checking boils down to check the existence of an infinite sequence i 1 ...i n ... such that v i1 ...v in .... = z i1 ...z in ..., which is actually the infinite PCP. There are two major differences between joint diagnosability in our framework and joint observability in (Tripakis, 2001). One is that the former assumes that local observers are attached to local components that are synchronized by common communication events to get a global model while the latter separates arbitrarily the observable events in the global model into several sets. The other one is that joint diagnosability consists in separating infinite trajectories while joint observability consists in separating finite ones. Thus, if any communication event is assumed to be unobservable, joint diagnosability checking boils down to infinite PCP. But this one has also been proved to be undecidable [START_REF] Halava | Undecidability of infinite post correspondence problem for instances of size 9[END_REF], which gives the result.

For the sake of simplicity, we give now a simpler and self-contained proof of undecidability of joint diagnosability that relies only on finite PCP. Theorem 2 Given a self-observed distributed system where communication events are unobservable, checking joint diagnosability of a given fault is undecidable. Proof : 1) To prove this theorem, without loss of generality, consider the example depicted in figure 2, where the system is composed of two components G 1 and G 2 . In G 1 , each one of V i, i ∈ {1, ..., k}, and each one of Zi, i ∈ {1, ..., k}, denotes a sequence of observable events all different from O1, C1, ..., Ck are unobservable communication events, F denotes a fault event and O1 is an observable event. In G 2 , each one of Ai, i ∈ {1, ..., k}, denotes an observable event different from O2, C1, ..., Ck are unobservable communication events and O2 is an observable event.

Then the observations in G 1 can be described as

V i 1 V i 2 ...V i n O1 * without fault or Zi 1 Zi 2 ...Zi n O1 *
with fault, where ∀i j , j ∈ {1, ..., n}, i j ∈ {1, ..., k}. In G 2 , the observations are Ai 1 Ai 2 ...Ai n O2 * . In this system, the occurrence of the fault can be confirmed by the observation of O1. Next we discuss about the case without the observation of O1.

2) Without the observation of O1, the local observations are wO1 + for G 1 and Ai 1 Ai 2 ...Ai n O2 * for G 2 , where w = V i 1 V i 2 ...V i n when there is no fault or w = Zi 1 Zi 2 ...Zi n when there is a fault. Clearly, if PCP has a solution, i.e., ∃(i m ) 1≤m≤n such that V i 1 V i 2 ...V i n = Zi 1 Zi 2 ...Zi n , we have two trajectories p and p′ such that the observations of p in G 1 are V i 1 V i 2 ...V i n O1 + , which is a trajectory without fault, while the observations of p′ in G 1 are Zi 1 Zi 2 ...Zi n O1 + , which is a trajectory with a fault. And both p and p′ have the same observations for G 2 , i.e., Ai 1 Ai 2 ...Ai n O2 * . Thus we get that p and p′ have the same observations for both G 1 and G 2 , i.e.,

V i 1 V i 2 ...V i n O1 + =Zi 1 Zi 2 ...Zi n O1 + for G 1 and Ai 1 Ai 2 ...Ai n O2 * for G 2 , then the fault is not jointly diagnosable.
3) On the other hand, if the fault is not jointly diagnosable, then we obtain at least one indeterminate pair, denoted by p and p′ such that the projection of

p on G 1 is Ci 1 V i 1 Ci 2 V i 2 ...Ci n V i n O1 * , on G 2 is Ai 1 Ci 1 Ai 2 Ci 2 ...Ai n Ci n O2 * and that of p′ on G 1 is Cj 1 Zj 1 Cj 2 Zj 2 ...Cj m Zj m F O1 * and on G 2 is Aj 1 Cj 1 Aj 2 Cj 2 ...Aj m Cj m O2 * .
From the fact that p and p′ have the same observations for G 2 , we get Ai 1 Ai 2 ...Ai n O2 * = Aj 1 Aj 2 ...Aj m O2 * and thus we have m = n and i 1 = j 1 , ..., i n = j n . And then from the same observations of p and p′ on G 1 , we get

V i 1 V i 2 ...V i n O1 * = Zi 1 Zi 2 ...Zi n O1 * , i.e., V i 1 V i 2 ...V i n = Zi Zi 2 ...Zi n ,
which means that there is a solution for PCP. This proves that the existence of a solution for PCP is equivalent to the fault being not jointly diagnosable. Since PCP is an undecidable problem, then checking joint diagnosability is undecidable.

ALGORITHM TO TEST A SUFFICIENT CONDITION

We have proved that joint diagnosability in selfobserved distributed systems with unobservable communication events is undecidable. We can nevertheless propose an algorithm to test a sufficient condition, which is still quite useful in some circumstances. The first step is to construct the local diagnoser from a given local model, which allows one to get fault information for any local trajectory. Then we show how to build the local twin plant, which allows one to obtain original information about indeterminate pairs (also called local indeterminate pairs in the following), based on the local diagnoser. The next step is to check the global consistency, i.e., to check whether the local indeterminate pairs can be extended into (global) indeterminate pairs, whose existence testifies non joint diagnosability. We give an algorithm to only test a sufficient but not necessary condition for global consistency. The proof is consistent with the undecidability result obtained in the precedent section. Actually our algorithm remains trivially applicable when the assumption of unobservability of communication events is partially relaxed, i.e.., in the most general case where some communication events are observable and others unobservable.

Original diagnosability information

From theorem 1, we know that joint diagnosability verification consists in checking the existence of indeterminate pairs in the system. In the distributed framework, we use the structure called local twin plant defined in [START_REF] Jiang | A polynomial time algorithm for diagnosability of discrete event systems[END_REF] to analyze joint diagnosability. In particular, the considered fault is assumed to only occur in one component, denoted by G f . Then the local twin plant for G f contains original information for indeterminate pairs: actually this twin plant is a FSM that compares every pair of local trajectories to search for the pairs with the same arbitrarily long local observations, but exactly one of the two containing a fault, which precisely we call local indeterminate pairs. First, we define an operation called delay closure with respect to a subset Σ d of Σ to preserve all information about the events in Σ d by abstracting away irrelevant parts. Definition 4 (Delay Closure). Given a FSM G = (Q, Σ, δ, q 0 ), its delay closure with respect to

Σ d ⊆ Σ is ∁ Σ d (G) = (Q, Σ d , δ d , q 0 ) where (q, σ, q ′ ) ∈ δ d iff ∃s ∈ (Σ\Σ d ) * , (q, sσ, q ′ ) ∈ δ.
We now describe how to construct the local diagnoser of a given component, based on which we build the local twin plant. Given a local model, we get a modified one by attaching fault label, denoted by l ∈ {N, F }, where N for normal and F for fault, to each state. In other words, before the occurrence of the fault, each state is labeled with label N and, after its occurrence, with label F . Definition 5 (Local diagnoser). Given a local model G i , its local diagnoser D i is obtained by operating the delay closure with respect to the set of communication events and observable events on the modified model:

D i = ∁ Σi o ∪Σi c (G m i ), where G m i is the modified ver- sion of G i .
Based on the local diagnoser, the corresponding local twin plant is obtained by synchronizing the local diagnoser with itself based on the locally observable events, allowing one to obtain all pairs of local trajectories with the same observations to search for local indeterminate pairs. To simplify this synchronization, the two identical local diagnosers, denoted by D l i for the left instance and D r i for the right instance, can be reduced as follows: D l i is obtained by retaining only paths with at least one fault cycle and D r i is obtained by retaining only paths with at least one non-fault cycle. This reduction keeps all necessary original diagnosability information since what we are interested in here are only local indeterminate pairs, i.e., pairs of local trajectories with only one containing the fault and both of them with the same observations. However, this reduction is only applicable for the construction of the local twin plant of the faulty component G f ; for other components, the local twin plant is obtained by synchronizing the non reduced left instance and the non reduced right instance since there is no fault information. Considering that this synchronization is based on the set of locally observable events Σ io , the nonsynchronized events are distinguished between the two instances by the prefix L or R: in D l i (D r i ), each communication event c ∈ Σ ic from D i is renamed by L : c (R : c). The names of all locally observable events are left unchanged.

Definition 6 (Local twin plant). Given a local diagnoser D i for the component G i , the corresponding local twin plant is a FSM, denoted by

T i = D l i D r i
, where the synchronized events are locally observable events in G i .

Each state of a local twin plant is a pair of local diagnoser states providing two possible diagnoses with the same local observations. Given a twin plant state ((q l , l l )(q r , l r )), where (q l , l l ) and (q r , l r ) are two diagnoser states and each one contains both a system state and a corresponding pattern state, if the considered fault f ∈ l l ∪ l r but f / ∈ l l ∩ l r , which means that the occurrence of f is not certain up to this state, then this state is called an ambiguous state with respect to the fault f . An ambiguous state cycle is a cycle containing only ambiguous states. In a local twin plant, if a path contains an ambiguous state cycle with at least one locally observable event, then it is called a local indeterminate path, which corresponds to a local indeterminate pair. Note that local indeterminate paths contain original diagnosability information and can be obtained only in the local twin plant of the component G f . If a local indeterminate pair can be extended into a global indeterminate pair, then we say that its corresponding local indeterminate path is globally consistent. Figure 3

Global consistency checking

Joint diagnosability verification consists in checking the existence of globally consistent local indeterminate paths, whose existence proves non joint diagnosability.

To check the global consistency of local indeterminate paths in a subsystem, we consider now two important issues: 1) the global consistency of the corresponding left trajectories of the local indeterminate paths in the local twin plant of the subsystem, shortly called left consistency checking in the following; 2) the global consistency of the corresponding right trajectories of the local indeterminate paths in the local twin plant of the subsystem, shortly called right consistency checking. Definition 7 (Left (Right) consistent plant). Given a subsystem G S composed of components G i1 , ..., G im and their corresponding local twin plants T i1 , ..., T im , to obtain a left (right) consistent plant with respect to the subsystem G S , denoted by T l f (T r f ), we perform the following two steps:

• Distinguish right (left) communication events between local twin plants by renaming them with the prefix of component ID. For example, R:C2 (L:C2) in the local twin plant of G 2 is renamed as G 2 :R:C2 (G 2 :L:C2).

• Synchronize the renamed local twin plants with

the synchronized events being the common left (right) communication events, which works because observable events do not intersect between components and non-synchronized right (left) communication events are distinguished by the prefix of component ID. In the left (right) consistent plant with respect to a subsystem G S , each path p corresponds to a set of paths p i1 , ..., p im in the local twin plants of all components in G S such that the set of left (right) trajectories of p i1 , ..., p im are reconstructible with respect to G S . For our example, the bottom part of figure 4 shows a part of the left consistent plant T l f , which is obtained by synchronizing the renamed local twin plant of G 1 and that of G 2 (top part of figure 4) based on the common left communication events.

Algorithm

Algorithm 1 presents the procedure to verify a sufficient condition of joint diagnosability. As shown in the pseudo-code, algorithm 1 performs as follows. Given the input as the set of component models, the fault F that may occur in the component G f , we initialize the parameters as empty, i.e., G l S , the subsystem for the left consistency checking and G r S , the subsystem for the right consistency checking. The procedure of the algorithm can be separated into two parts: left consistency checking (line 3-12) and right consistency checking (line 13-24). Left consistency checking begins with the local twin plant construction of G f , the subsystem G l S being now G f (line 3-4). As long as both the left consistent plant T l f with respect to the current left subsystem G l S and DirectCC(G, G l S ) are not empty (line 5), where DirectCC(G, G l S ) is the set of directly connected components to the subsystem G l S (a directly connected component being one sharing at least one common communication event with the subsystem) Algorithm 1 Algorithm to check a sufficient condition of joint diagnosability in the general case 1: INPUT: the system model G = (G 1 , ..., G n ); the fault F and the faulty component G f 2: Initializations: G l S ← ∅ (subsystem for left consistency checking); G r S ← ∅ (subsystem for right consistency checking) 3:

T l f ← ConstructLT P (G f ) 4: G l S ← G f 5: while T l f = ∅ and DirectCC(G, G l S ) = ∅ do 6: G i ← SelectDirectCC(G, G l S )
7:

T i ← ConstructLT P (G i ) 8: T l f ← T l f T i 9: G l S ← Add(G l S , G i ) 10: T l f ← RetainConsisP aths(T l f ) 11: if T l f = ∅ then 12:
return "F is jointly diagnosable in G" 13: else 14:

T r f ← AbstractRight(G f , T l f ) 15: G r S ← G f 16: while T r f = ∅ and G l S = G r S do 17: G i ← SelectDirectCC(G l S , G r S ) 18: T r f ← T r f AbstractRight(G i , T l f ) 19: G r S ← Add(G r S , G i ) 20: T r f ← RetainConsisP aths(T r f ) 21: if T r f = ∅ then 22:
return "F is jointly diagnosable in G"

23:

else 24:

return "Joint diagnosability cannot be determined"

the algorithm repeatedly performs the following steps to further check left consistency in an extended subsystem: 1. Select one directly connected component G i to the subsystem G l S and construct its local twin plant T i (line 6-7).

Synchronize T l

f with T i to obtain left consistent plant for this extended subsystem based on the set of common left communication events (line 8). To do this, the set of non-synchronized right communication events are distinguished by the prefix of component ID.

Update the subsystem G l

S by adding G i and reduce the newly obtained T l f by retaining only paths with ambiguous state cycles containing observable events for all components in G l S (line 9-10).

If the left consistent plant T l f is empty, then there is no local indeterminate path that corresponds to a set of paths in the local twin plants of all components in the subsystem such that their left trajectories are reconstructible, which implies the non existence of a globally consistent local indeterminate path. In this case joint diagnosability information is returned (line 11-12). Otherwise, if T l f is not empty (line 13), then we proceed to check right consistency of the corresponding paths in T l f that have been already verified to be left consistent in the whole system. Right consistency checking begins with the function AbstractRight(G f , T l f ) (line 14), which performs delay closure with respect to the set of right communication events and observable events of G f . In this way, what we obtain does not contain left communication events. Then the subsystem G r S is assigned as G f (line 15). As long as the right consistent plant T r f for the current right subsystem G r S is not empty and G l S = G r S (line 16), we repeatedly perform the following steps to check right consistency in an extended subsystem (note that since left consistency checking does explore all connected components, during right consistency checking we only need to consider the subsystem G l S instead of the whole system):

1. Select a directly connected component G i to G r S from G l S (line 17). 2. Perform the function AbstractRight(G i , T l f
), which has already been described above, and then synchronize with T r f based on the set of common right communication events (line 18). To do this, we rename the right communication events by removing the prefix of component ID, e.g., G i :R:C2 renamed as R:C2.

Update the subsystem G r

S by adding G i and reduce the newly obtained T r f by retaining only paths with ambiguous state cycles containing observable events for all components in G r S (line 19-20). If the right consistent plant T r f is empty, then there is no local indeterminate path that corresponds to a set of paths in the local twin plants such that their left trajectories and right trajectories are reconstructible respectively, i.e., there is no globally consistent local indeterminate path. In this case, the algorithm returns joint diagnosability information (line 21-22). Otherwise, if T r f is not empty, actually we cannot determine whether the fault is jointly diagnosable or not. Then the algorithm returns the information about the indetermination of joint diagnosability (line 23-24). In other words, empty left consistent plant T l f or empty right consistent plant T r f is a sufficient condition but not a necessary condition of joint diagnosability.

Theorem 3 In algorithm 1, if the left consistent plant T l

f or the right consistent plant T r f is empty, then the fault is jointly diagnosable, but the reverse is not true.

Proof : (⇒) Suppose that T l f or T r f is empty and that the fault is not jointly diagnosable. From non joint diagnosability, it follows that there exists at least one globally consistent local indeterminate path. Since global consistency of a local indeterminate path implies both left consistency and right consistency, from algorithm 1 we know that, after left and right consistency checking, this local indeterminate path must correspond to a path both in T l f and in T r f . Thus neither T l f nor T r f is empty, which contradicts the assumption. ( ) Now we explain why non emptiness of both T l f and T r f does not necessarily imply that the fault is not jointly diagnosable. Suppose that T l f is not empty and that it contains two paths, denoted by ρ1 and ρ2, corresponding to two local indeterminate paths. ρ1 corresponds to a set of paths ρ 1 i , 1 ≤ i ≤ n in the local twin plants of all components and ρ2 corresponds to a set of paths ρ 2 i , 1 ≤ i ≤ n in all local twin plants. Now suppose that the right trajectories of the set of paths ρ 1 i , 1 ≤ i ≤ n are not reconstructible and the same for that of the set of paths ρ 2 i , 1 ≤ i ≤ n. It ollows that the two local indeterminate paths cannot be extended into global indeterminate pairs and thus are not globally consistent. Then we further suppose that the right trajectories of the set of paths ρ 1 1 , ..., ρ 1 n-1 , ρ 2 n are reconstructible or the same for the set of paths ρ 2 1 , ..., ρ 2 n-1 , ρ 1 n . In this case, from algorithm 1, it follows that finally the right consistent plant T r f is not empty. Now both T l f and T r f are not empty but there is no globally consistent local indeterminate paths, i.e., the fault is jointly diagnosable. Now illustrate on our example the fact that the condition is not necessary. The top part of figure 5 shows the results of performing delay closure with respect to the set of right communication events and locally observable events both for G 1 and G 2 on the left consistent plant depicted in the bottom part of figure 4. Then, to check right consistency, we rename again the right communication events by removing the component ID such that they can be synchronized. Finally the bottom part of figure 5 shows a part of the right consistent plant, which is not empty. Now both left and right consistent plants are not empty, but this does not imply the existence of global indeterminate pairs that witness non joint diagnosability. Actually the part of the left consistent plant depicted here corresponds to two local indeterminate pairs in G 1 with their corresponding left consistent pairs in G 2 , i.e., one local indeterminate pair is ((C1.O1.F.O1 * ), (O1.C2.O1 * )) in G 1 with its left consistent pair ((C1.O3.O5 * ), (O3.U 2.O5 * )) in G 2 and the other local indeterminate pair is ((O2.U 1.C2.F.O1 * ), (C1.O2.C2.O1 * )) in G 1 with its left consistent pair ((O4.C2.O5 * ), (O4.C2.O5 * )) in G 2 . While the right consistent plant shown here corresponds to one local indeterminate pair in G 1 , which is ((C1.O1.F.O1 * ), (O1.C2.O1 * )), with its corresponding right consistent pair in G 2 , i.e., ((O4.C2.O5 * ), (O4.C2.O5 * )). Thus we can see that the same local indeterminate pair does not correspond to the same consistent pair in G 2 in the left consistent plant and in the right consistent plant, which means that this local indeterminate pair cannot be extended into a global indeterminate pair. So our algorithm gives indeterminate information for joint diagnosable systems that satisfy the following condition: for any set of paths, i.e., one path in the local twin plant of each component and one in the plant for faulty component being a local indeterminate path, that are left consistent and right consistent respectively, then we have that their corresponding local trajectories in the components cannot constitute an indeterminate pair through synchronization.

DECIDABLE CASE

We have proved the undecidability of joint diagnosability when communication events are unobservable. If we assume their observability, then this problem becomes decidable. The reason is that, when any communication event is observable, then in the local twin plant of each component, we obtain all pairs of local trajectories with the same observations, including the same observable communication events. In other words, each path in the local twin plant corresponds to a pair of local trajectories with the same sequence of communication events. It follows that to check global consistency of local indeterminate paths, what was a separate checking for left and right consistency becomes only one checking. While in algorithm 1, the checking into two separate phases is the reason why it gives only a sufficient but not necessary condition for joint diagnosability. Actually, the observability of communication events makes joint diagnosability equivalent to classical diagnosability since only one checking for global consistency implies the same global occurrence order of observations for global indeterminate pairs.

Algorithm 2 presents the procedure to check joint diagnosability when communication events are assumed to be observable. Taking the system model and the faulty component as input, the parameter, i.e., the current subsystem G S , is initialized as empty. Then the algorithm begins with the construction of the local twin plant of the faulty component G f , the current subsystem becoming G f (line 3-4). Here we emphasize that, due to the observability of communication events, then the local twin plant should be constructed by synchronizing two instances based on the set of observable events and the set of communication events, i.e., here communication events do not need to be distinguished by changing their names. As long as both T f and DirectCC(G, G S ) are not empty (line 5), then the following steps are repeatedly performed:

• Select one component G i directly connected to the current subsystem G S and then construct its local twin plant T i , which is obtained by first operating delay closure with respect to the set of communication events and observable events and then by synchronizing the two instances based on all events, i.e., the set of communication events and observable events. (line 6-7) return "F is jointly diagnosable in G" 13: else 14: return T f

CONCLUSION

In this paper, we consider self-observed distributed systems such that observable events can only be observed by their own components. Clearly, we do not need the monolithic model of the system, thus the distributed and private (w.r.t. observation) nature of real systems is taken into account. Then we prove the undecidability of joint diagnosability checking when communication events are unobservable, before proposing an algorithm to test a sufficient condition of joint diagnosability in this case. To check the non existence of indeterminate pairs in the system, we start from local indeterminate paths in the local twin plant and then we check both in sequence left consistency and right consistency. Note that, due to the observation-privacy, the global occurrence order of observable events between different components is not known, which is taken into account through constructing left and right consistent plants separately. At the opposite, in the approaches for DES with globally observable events, twin plant is constructed by incrementally synchronizing local twin plants via both left and right communication events at the same time, which means that in their case, knowledge of the global occurrence order is required. Similar to the distributed algorithms for classical diagnosability ( (Pencolé, 2004), [START_REF] Schumann | Scalable diagnosability checking of event-driven systems[END_REF], etc.), our algorithm has to construct some part of a global structure, but much less than in the centralized approach (in particular the components not involved in the algorithm have their model completely not disclosed) and this is normally unavoidable for off-line diagnosability analysis. Then we discuss the decidable case where communication events are observable by giving an algorithm that checks joint diagnosability. However, the decidable case in [START_REF] Ye | Diagnosability analysis of discrete event systems with autonomous components[END_REF] is not the same as that in this paper. The former also deals with unobservable communication events, which becomes decidable because of the assumption of exhaustive enumeration about local paths. While here the decidable case is about observability of communication events. We see that there is a gap between these two cases as the unobservable case is undecidable and the observable case is decidable. Next interesting work is to investigate where is the frontier between these two cases, i.e., to study the decidability of joint diagnosability for partial observability of communication events.

Figure 1 :

 1 Figure 1: A system with two components G1 (left) and G2 (right).

Figure 2 :

 2 Figure 2: A system with two components G1 (top) and G2 (bottom) for proving undecidability of joint diagnosability.

Figure 3 :

 3 Figure 3: Two reduced instances of the diagnoser for G1 (top) and part of the corresponding local twin plant (bottom).

  shows the left and right instances of the local diagnoser for the faulty component G 1 of Figure 1 (top) as well as a part of the corresponding local twin plant (bottom). Clearly, in the local twin plant, we have local indeterminate paths since they have ambiguous state cycles with observable events.

Figure 4 :

 4 Figure 4: Part of the renamed local twin plants for G1 and G2 (top) and part of the left consistent plant T l f (bottom).

Figure 5 :

 5 Figure 5: FSM after delay closure on the left consistent plant (figure 4) for G1 and G2 (top) and part of the right consistent plant (bottom).

Algorithm 2

 2 Algorithm for checking joint diagnosability with observability of communication events1: INPUT: the system model G = (G 1 , ..., G n ); the fault F and the faulty component G f 2: Initializations: G S ← ∅ (subsystem considered for current checking) 3: T f ← ConstructLT P (G f ) 4: G S ← G f 5: while T f = ∅ and DirectCC(G, G S ) = ∅ do 6: G i ← SelectDirectCC(G, G S ) 7: T i ← ConstructLT P (G i ) 8: T f ← T f T i 9: G S ← Add(G S , G i ) 10:T f ← RetainConsisP aths(T f ) 11: if T f = ∅ then 12:

•

  Synchronize the current local twin plant T f and T i based on the common communication events of G S and G i . (line 8) • Update the current subsystem by adding this selected component and keep only the paths in the During this procedure, if the local twin plant T f for the current subsystem happens to be empty, which means that there is no path that contains ambiguous state cycle with observations for all concerned components, thus there is no local critical path that is globally consistent and the algorithm returns joint diagnosability information (line 11-12). Otherwise, if at the end the final FSM is not empty, it is returned by the algorithm as non joint diagnosability information (line 13-14) as any path in it corresponds to a globally consistent local indeterminate path. The reason is that if the communication events are observable, then any path in a local twin plant corresponds to a pair of local trajectories with the same observations, including the same communication events. So with the assumption of observability of communication events, joint diagnosability checking becomes decidable, whose verification is provided by the algorithm 2.

	newly obtained FSM that contain ambiguous state
	cycles with observations for all involved compo-
	nents. (line 9-10)
	7