Hamzah Bakhti
email: hamzah.bakhti@um5s.net.ma

Lahcen Azrar
email: l.azrar@um5s.net.ma

Mahmoud Hamadiche
email: mahmoud.hammadich@gmail.com

Coupled Numerical Scheme for Vascular Fluid-Tube Interaction using Large Deformation Theory. Application to Aneurysm Dynamics under Pulsatile Flow

A methodological approach is elaborated to study the vascular fluid-tube interaction under pulsatile blood flow within an asymmetrical nonlinear aneurysm and large deformation models. The aneurysm dynamics were modeled using a simplified Lagrangian nonlinear system describing the arterial wall motion with large deformation. The flow is governed by the Navier-Stokes equations in two-dimensional domain. A semi-implicit splitting scheme coupled with the finite difference method is developed to solve the flow equation in an irregular domain using a mesh transformation. On the other hand, the wall equations are solved using the Runge-Kutta method combined with the shooting technique by reducing the resulted boundary value problem to an initial value problem. Rigid and elastic aneurysms are considered and the effect of various geometrical and fluid parameters on the flow are investigated, mainly the Reynolds number, the aneurysm length and height. The study focused on the effect of the asymmetric curvatures of the tube walls and their deformations due to the pulsatile flow. The flow is examined under a steady inlet flow as well as under a pulsatile one for various aneurysm forms. The obtained numerical results validate the fluid and structure solvers and demonstrate significant differences between the rigid and elastic models of the structure as well as the effect of the asymmetric propriety of the arterial aneurysm.

Fluid-structure interaction problem (FSI) occurs in many cases when clear deformations of a structure are caused by fluid flow near or around it [START_REF]Fluid-Structure Interaction[END_REF][START_REF] Sun | An accurate fsi-sph modeling of challenging fluid-structure interaction problems in two and three dimensions[END_REF][START_REF] Khayyer | A coupled incompressible sph-hamiltonian sph solver for hydroelastic fsi corresponding to composite structures[END_REF][START_REF] Zimmermann | Quantitative hemodynamics in aortic dissection: Comparing in vitro mri with fsi simulation in a compliant model[END_REF]. Depending on the fluid flow and the properties of the structure's material, these deformations can be large or small [START_REF] Figueroa | A computational framework for fluid-solid-growth modeling in cardiovascular simulations[END_REF][START_REF] Frei | Long-term simulation of large deformation, mechanochemical fluid-structure interactions in ALE and fully Eulerian coordinates[END_REF]. For the case of large deformations, the flow field and the structure deformations affect each other in a bidirectional way. Thus, the problem should be treated as a coupled multiphysics analysis. However, when large deformations and large displacements are involved, careful analysis is needed to tackle the domain displacement and the associated nonlinear deformation. In general, the mathematical modeling of FSI problems using large deformation theory results in complex and nonlinear systems of equations that remain a challenge to many researchers for years.

To better understand some cardiovascular diseases, many numerical techniques for FSI models were developed to study blood flow behavior through different types of arteries. Particularly, with arterial aneurysms, the abnormal dilation of the arterial wall results in blood vessel wall weakness. This can be life threatening as it leads to a higher risk of an arterial wall rupture and a more common type is abdominal aortic aneurysms (AAA) [START_REF] Mcphee | The impact of gender on presentation, therapy, and mortality of abdominal aortic aneurysm in the United States, 2001-2004[END_REF][START_REF] Aggarwal | Abdominal aortic aneurysm: A comprehensive review[END_REF]. By considering the blood flow to be pulsatile, the FSI models are favorable to understanding both the behavior of the blood and the dynamics of the arterial aneurysm [START_REF] Peattie | Pulsatile Flow in Fusiform Models of Abdominal Aortic Aneurysms: Flow Fields, Velocity Patterns and Flow-Induced Wall Stresses[END_REF]. However, the numerical techniques are usually time-consuming due to the complexity of vascular FSI models because of the coupling between the mechanobiology of the wall and the pulsatile blood velocity-pressure fields [START_REF] Roy | A literature review of the numerical analysis of abdominal aortic aneurysms treated with endovascular stent grafts[END_REF]. Therefore, the study of these models, especially in a large deformation domain with the appropriate inlet/outlet boundary conditions, and using timesaving numerical schemes is very useful to understand the dynamics of arterial aneurysm due to pulsatile blood flow.

The study of FSI models is of considerable importance and has direct implications for many industrial, engineering and medical problems. These models draw the attention of researchers for years and many scientific works were published studying a variety of phenomena with applications in various fields, such as pumps [START_REF] Li | Thermal-fluid-structure coupling analysis for valve plate friction pair of axial piston pump in electrohydrostatic actuator (EHA) of aircraft[END_REF], heart valves [START_REF] Kamensky | An immersogeometric variational framework for fluidstructure interaction: Application to bioprosthetic heart valves[END_REF], aerodynamic flutter [START_REF] Gan | Delayed detached Eddy simulation of wing flutter boundary using high order schemes[END_REF], shock absorbers [START_REF] Graczykowski | Mathematical models and numerical methods for the simulation of adaptive inflatable structures for impact absorption[END_REF], turbines [START_REF] Yan | Computational free-surface fluidstructure interaction with application to floating offshore wind turbines[END_REF] and biomechanics [START_REF] Piatti | Hemodynamic and thrombogenic analysis of a trileaflet polymeric valve using a fluid-structure interaction approach[END_REF][START_REF] Chiastra | On the necessity of modelling fluid-structure interaction for stented coronary arteries[END_REF][START_REF] Sundaram | Aortic valve dynamics using a fluid structure interaction model -The physiology of opening and closing[END_REF]. The analytical solutions are almost impossible to obtain, except for rare cases where many simplifications are assumed to the problem [START_REF] Pan | Analytical solutions for the dynamic analysis of a valveless micropump -A fluid-membrane coupling study[END_REF][START_REF] Turkyilmazoglu | Air blast response of compaction foam having a deformable front face panel incorporating fluid structure interactions[END_REF]. Thus, many numerical procedures were developed for this matter [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluidstructure interaction[END_REF][START_REF] Cottet | A semi-implicit level set method for multiphase flows and fluid-structure interaction problems[END_REF][START_REF] Miyauchi | A numerical method for interaction problems between fluid and membranes with arbitrary permeability for fluid[END_REF]. Meanwhile, various FSI models were proposed using the large deformation theory, worth mentioning the work of Tian et. al. [START_REF] Tian | Fluid-structure interaction involving large deformations: 3D simulations and applications to biological systems[END_REF] in which they studied a three-dimensional FSI model involving a large deformation theory with applications in biological systems. Various other works are elaborated using the same theory, as presented in [START_REF] Hart | A twodimensional fluid-structure interaction model of the aortic value[END_REF][START_REF] Paik | Fluid-structure interaction for an elastic structure interacting with free surface in a rolling tank[END_REF][START_REF] Kumar | Bi-directional fluid-structure interaction for large deformation of layered composite propeller blades[END_REF].

Many research works studying the blood flow in arteries are based on the assumption of rigid deformed walls. Bakhti et. al. [START_REF] Bakhti | Steady Flow of Couple-Stress Fluid in Constricted Tapered Artery: Effects of Transverse Magnetic Field, Moving Catheter, and Slip Velocity[END_REF][START_REF] Bakhti | Pulsatile blood flow in constricted tapered artery using a variable-order fractional Oldroyd-B model[END_REF] studied the steady and pulsatile blood velocity profile and the resulting shear stress on a rigid asymmetric stenosed tapered artery using different rheological blood models. Finol and Amon [START_REF] Finol | Blood Flow in Abdominal Aortic Aneurysms: Pulsatile Flow Hemodynamics[END_REF] studied pulsatile blood flow and hemodynamic stresses in arteries with the presence of two, axisymmetric and rigid wall aneurysms using the spectral element method. Later on, Finol et al. used a three-dimensional model of AAAs with a single, asymmetric and rigid aneurysm using the finite element method [START_REF] Finol | The Effect of Asymmetry in Abdominal Aortic Aneurysms Under Physiologically Realistic Pulsatile Flow Conditions[END_REF]. More recently, more realistic numerical studies on blood dynamics based on rigid wall models are elaborated in [START_REF] Rabby | Pulsatile non-newtonian laminar blood flows through arterial double stenoses[END_REF][START_REF] Reorowicz | Numerical simulations of the blood flow in the patient-specific arterial cerebral circle region[END_REF][START_REF] Rabby | Laminar blood flow through a model of arterial stenosis with oscillating wall[END_REF][START_REF] Apostolidis | Non-Newtonian effects in simulations of coronary arterial blood flow[END_REF][START_REF] Shupti | Rheological behavior of physiological pulsatile flow through a model arterial stenosis with moving wall[END_REF][START_REF] Imai | Numerical methods for simulating blood flow at macro, micro, and multi scales[END_REF][START_REF] Shupti | Pulsatile non-newtonian fluid flows in a model aneurysm with oscillating wall[END_REF]. By considering the blood flow to be pulsatile, the FSI models are favorable for understanding both the behavior of the blood and the dynamics of the arterial aneurysm [START_REF] Peattie | Pulsatile Flow in Fusiform Models of Abdominal Aortic Aneurysms: Flow Fields, Velocity Patterns and Flow-Induced Wall Stresses[END_REF]. Worth mentioning the work of Hamadiche [START_REF] Hamadiche | Numerical simulation and mathematical analysis of flow-wall interaction in the large deformation application to the dynamics of the aneurysms[END_REF], where the dynamics of symmetrical aneurysms in large deformations domain due to pulsatile blood flow are studied.

However, due to the high computation costs (run time, required memory) and the mathematical complexity of the vascular FSI models for blood flow, the existing numerical techniques are difficult to code and need powerful computing machines. For this matter, many numerical methods were used to study FSI models in blood flow problems by applying techniques such as semiimplicit approaches and splitting schemes to reduce the computation costs. Also, reduced models seem to be helpful [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Guidoboni | A kinematically coupled time-splitting scheme for fluid-structure interaction in blood flow[END_REF][START_REF] Bukač | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF][START_REF] Formaggia | On the stability of the coupling of 3D and 1D fluid-structure interaction models for blood flow simulations[END_REF]. Actually, the geometrical and fluid parameters, such as the model parameters of the arterial aneurysm and the flow Reynolds number, have a great effect on the stability and convergence of these numerical schemes and that is due to their nature (explicit, semi-implicite, etc...). There is still a large need for more improved numerical techniques for such FSI problems using large deformation theory that are in general challenging in term of the high complexity and strong nonlinearity of these types of problems.

In this paper a mathematical model describing the wall motion in the Lagrangian formulation and using the large deformations theory is formulated. This model leads to a system of nonlinear equations governing the wall motion coupled with the Navier-Stokes equations for the flow motion. An adaptive numerical scheme is presented to simulate the FSI model using multiple techniques, namely, the Runge-Kutta coupled with the shooting method to solve the structure equations. The finite difference method (FDM) with a mesh transformation technique is elaborated to solve the flow equation in a twodimensional domain combined with a semi-implicit splitting scheme. Multiple numerical experiments and tests have been elaborated in order to validate the proposed numerical scheme. The effect of various fluid and structure variables on the physiological parameters such as fluid velocity, pressure and wall displacement are deeply investigated. Results are obtained from elaborated methodological approaches to investigate blood flow and using the proposed vascular FSI model by considering two non-symmetrical aneurysms.

The paper is organized as follows: In section 2, the governing equations for the fluid (Section 2.2) are introduced, the derivation of the structural model for the wall motion is then presented (Section 2.1) and both models are completed by the appropriate initial and the coupling boundary conditions (Section 2.3). In section 3, a semi-implicit numerical scheme is introduced for the FSI solver which couples the Runge-Kutta and shooting method for the structure solver (Section 3.1) with the FDM for the fluid solver using a mesh transformation technique (Section 3.2). In section 4, a series of numerical tests were proposed to validate the structure and fluid solvers with a given exact solution for the structural model (Section 4.1) and classical blood flow problems in rigid arteries with the presence of different types of arterial aneurysms (section 4.2). The results of an FSI problem is then presented by coupling the wall model with the fluid equations to study the blood behavior in elastic arteries with the presence of asymmetrical aneurysms in both artery walls using the proposed nonlinear structure model (section 4.3). Finally, the main observations and results are given in the conclusion.

Mathematical Formulation

For the purpose of this paper, the flow of an incompressible Newtonian fluid in a two dimensional domain defined by a channel with both rigid and elastic parts of the parallel walls, as shown in Figure 1, is considered. The wall is assumed to be an incompressible continuum and its displacement is observed in the two dimensional reference frame (o, ξ, η).

Let d 1 and d 2 be the Lagrangian variables that relate the deformed configuration with the reference configuration of the wall and

d (ξ, η) = (d 1 , d 2 )
t be the displacements vector of a particle in the continuum body.

On the other hand, u and v donate the velocities of the fluid in the x and y directions, respectively. x denotes the vector position of a point in the original frame (o, x, y), which is considered as an orthonormal basis where a two-dimensional flow of an incompressible and viscous fluid takes place in the (x, y)-plane.

A mathematical modeling of the wall kinematics and the governing equations of motion are elaborated herein, as well as the resulting basic flow equations.

Fluid Flow equations

The flow parameters are assumed to depend on only two space variables, x, y. The fluid flow is governed by the two-dimensional nonlinear Navier-Stokes equations for an incompressible viscous Newtonian fluid. The continuity and momentum equations are given in a compact form as follow -Continuity equation

∇.u = 0 (1)
-Momentum equation

ρ ∂u ∂t + u.∇u = -∇p + µ∆u ( 2
)
where ρ is the fluid density, u is the flow velocity, p is the pressure and µ is the fluid viscosity. The system of equations modeling the flow field is given in the following explicit form

                 ∂u ∂x + ∂v ∂y = 0 ρ ∂u ∂t + u ∂u ∂x + v ∂u ∂y = - ∂p ∂x + µ ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ρ ∂v ∂t + u ∂v ∂x + v ∂v ∂y = - ∂p ∂y + µ ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 (3)
The non-dimensional Navier-Stokes equations can be obtained using the following dimensionless variables

x * = x L , u * = u R 0 , t * = t L/R 0 , p * = p ρU 2 (4)
where U is a reference velocity proportional to the maximum inlet velocity at the inlet of the channel, R 0 is the half-width of the rigid part of the artery. Substituting the scales, the non-dimensionalized equations after dropping the stars are obtained

                 ∂u ∂x + ∂v ∂y = 0 ∂u ∂t + u ∂u ∂x + v ∂u ∂y = - ∂p ∂x + 1 Re ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ∂v ∂t + u ∂v ∂x + v ∂v ∂y = - ∂p ∂y + 1 Re ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 (5)
where Re = ρU R0 µ is the Reynolds number.

Structure equation

Lagrangian description of wall kinematics

Following the Lagrangian description of motion and denoting by d(ξ, t) the displacements vector that depends only on ξ and time t, the resulted deformation gradient tensor F can be calculated using the relation F = I + ∇d, that is

F = 1 + ∂d1 ∂ξ 0 ∂d2 ∂ξ 1 ( 6
)
Following the assumption that the wall is an incompressible continuum, the deformation gradient must satisfy the condition det F = 1, thus ∂d1 ∂ξ = 0. Then

F = 1 0 ∂d2 ∂ξ 1 (7)
The displacement gradient tensor is then e = ∇d = F -I. For a large deformation domain, the Lagrangian finite strain tensor is defined as a function of the displacement gradient tensor

E = 1 2 e + e t + e t e (8)
Thus, using [START_REF] Mcphee | The impact of gender on presentation, therapy, and mortality of abdominal aortic aneurysm in the United States, 2001-2004[END_REF], one gets

E = 1 2 ∂d2 ∂ξ 2 ∂d2 ∂ξ ∂d2 ∂ξ 0 ( 9
)
The continuum is assumed to be a hyperelastic material, which is a special case of a Cauchy elastic material. The following St. Venant-Kirchhoff model is used

S = [-p + λ s Tr(E)] I + 2µ s E ( 10
)
where S is the second Piola-Kirchhoff stress, p the hydrostatic pressure of the fluid, λ s and µ s are the structure Lamé constants. The first Piola-Kirchhoff tensor denoted by P is given by

P = FS = S + eS (11)
which yields to

P = [-p + λ s Tr E] F + 2µ s FE (12)
Using a matrix form, one gets

P = -p + λ s 2 ∂d 2 ∂ξ 2 1 0 ∂d2 ∂ξ 1 + µ s ∂d2 ∂ξ 2 ∂d2 ∂ξ ∂d2 ∂ξ + ∂d2 ∂ξ 3 ∂d2 ∂ξ 2 (13)

Governing equations of wall motion

Let us assume that the continuum is subjected to a body force f, and let d denotes the displacement vector, ρ s the mass density in the reference frame. The governing equation of motion for a deformable continuum can be expressed as

ρ s ∂ 2 d ∂t 2 = ∇ • P + f (14)
Let us denote by P int the internal fluid pressure on both sides of the wall, obtained from the fluid equations, P ext is a fixed external pressure applied

P int P ext h n ξ η O θ
Fig. 2: Sketch of the body force applied on the wall.

on the outer tube wall and h is the wall thickness in the direction η of the reference frame. The body force can then be expressed as

f = P int -P ext h n ( 15
)
where n is the outward normal to the wall as shown in Figure 2. For the wall that is defined implicitly by Γ (ξ, η) = ηd 2 (ξ) = 0 the normal vector is then given by

n = ∇Γ ∂Γ ∂ξ 2 + ∂Γ ∂η 2 =    - ∂d 2 ∂ξ 1+ ∂d 2 ∂ξ 2 1 1+ ∂d 2 ∂ξ 2    (16)
Thus, the explicit form of the body force, Eq. 15, is given by

             f 1 = - P int -P ext h ∂d2 ∂ξ 1 + ∂d2 ∂ξ 2 f 2 = P int -P ext h 1 1 + ∂d2 ∂ξ 2 (17)
The momentum equation in ξ-direction can be explicitly expressed as

∂ -p + λs 2 + µ s ∂d2 ∂ξ 2 ∂ξ + f 1 = 0 (18)
For the sake of condensed notations a new variable is defined as fellow

ϕ = -p + λ s 2 + µ s ∂d 2 ∂ξ 2 (19)
Using Eq.19, Eq. 18 reduces to

∂ϕ ∂ξ - P int -P ext h ∂d2 ∂ξ 1 + ∂d2 ∂ξ 2 = 0 (20)
On the other direction η, the momentum equation is explicitly expressed as

ρ s ∂ 2 d 2 ∂t 2 = ∂ ϕ ∂d2 ∂ξ ∂ξ + µ s ∂ 2 d 2 ∂ξ + P int -P ext h 1 1 + ∂d2 ∂ξ 2 (21)
After further developments and using Eq. ( 20), one obtains

ρ s ∂ 2 d 2 ∂t 2 = (ϕ + µ s ) ∂ 2 d 2 ∂ξ 2 + P int -P ext h 1 + ∂d 2 ∂ξ 2 (22)
The final set of equations describing the elastic wall motion are thus given by

                       ϕ = -p + λ s 2 + µ s ∂d 2 ∂ξ 2 ∂ϕ ∂ξ - P int -P ext h ∂d2 ∂ξ 1 + ∂d2 ∂ξ 2 = 0 ρ s ∂ 2 d 2 ∂t 2 = (ϕ + µ s ) ∂ 2 d 2 ∂ξ 2 + P int -P ext h 1 + ∂d 2 ∂ξ 2 (23)
In addition to the dimensionless fluid parameters given in (4), the following dimensionless structural parameters are used

ξ * = ξ R 0 , h * = h R 0 , P * int = P int ρU 2 , P * ext = P int ρU 2 , d * 2 = d 2 R 0 , ϕ * = ϕ ρU 2 , ρ * s = ρ s ρ , µ * s = µ ρU 2 , λ * s = λ s ρU 2 (24)
The obtained dimensionless structural equations after dropping the stars are given by the same previous system of equations

                       ϕ = -p + λ s 2 + µ s ∂d 2 ∂ξ 2 ∂ϕ ∂ξ - P int -P ext h ∂d2 ∂ξ 1 + ∂d2 ∂ξ 2 = 0 ρ s ∂ 2 d 2 ∂t 2 = (ϕ + µ s ) ∂ 2 d 2 ∂ξ 2 + P int -P ext h 1 + ∂d 2 ∂ξ 2 (25)

Initial and Boundary Conditions

The flow equations are completed by the appropriate inlet/outlet initial and boundary conditions. In this work, a periodic axial velocity profile at the inlet of the channel is considered and given by u = u in (y, t), v = 0 at the inlet [START_REF] Paik | Fluid-structure interaction for an elastic structure interacting with free surface in a rolling tank[END_REF] where u in is a known time-depended function to be introduced later. At the other end of the channel, the free outlet boundary condition is imposed as

∂u ∂x = ∂v ∂x = 0, p = p 0 at the outlet ( 27
)
where p 0 is the outlet pressure. On the channel wall the no-slip condition at the rigid solid wall and the continuity of the velocity at the elastic wall are assumed.

   u = v = 0 at the rigid wall u = 0, v = ∂d 2 ∂t at the wall (28)
For the initial conditions, a fixed velocity profile is assumed

u = u 0 (x, y), v = v 0 (x, y) at initial time t 0 ( 29
)
where u 0 and v 0 are known functions to be defined later.

For the structural equations, both ends of the elastic curvature are assumed to be rigid, thus d 2 = 0 at both ends of the elastic curvature (30)

Numerical Method

Numerical methodological approaches and a numerical code have been elaborated to solve the resulting fully coupled fluid and structure models. The code has been tested and several applications of the presented numerical method have been studied to ensure its effectiveness.

In the following subsections, the elaborated numerical techniques used to solve the resulting structural and fluid equations, as well as the algorithm of the FSI coupling, are presented. For the structural system, the coupling of the Runge-Kutta method with the shooting techniques is elaborated to switch the boundary value problem to an initial value one.

As for the flow equations, the FDM method is used to calculate the velocitypressure fields in an irregular domain, where a mesh transformation technique that maps the complex flow domain into a regular one is used. For the time scheme, the pressure-correction method splitting the pressure from the velocity and reducing the CPU time is used. A forward finite differential scheme is then used to evaluate the time derivative.

Wall equations

The exterior pressure P ext is supposed constant and once the interior pressure P int is known, the structural system of equations can be solved using the Runge-Kutta method combined with the shooting technique. The time derivative is descritized using the following centered finite difference scheme

∂d n+1 2 ∂t ≃ d n+1 2 -2d n 2 + d n-1 2 ∆t 2 (31)
The time-discretized wall equations are then expressed as

                           ϕ n+1 = -p + λ s 2 + µ s ∂d n+1 2 ∂ξ 2 ∂ϕ n+1 ∂ξ - P n+1 int -P ext h ∂d n+1 2 ∂ξ 1 + ∂d n+1 2 ∂ξ 2 = 0 ρ s d n+1 2 -2d n 2 + d n-1 2 ∆t 2 = ϕ n+1 + µ s ∂ 2 d n+1 2 ∂ξ 2 + P n+1 int -P ext h 1 + ∂d n+1 2 ∂ξ 2 (32
) which forms a system of equations to be solved as an initial value problem. Thus, the shooting technique is coupled to the Runge-Kutta method to solve the problem on the space variable ξ.

Shooting method

For the structure, one gets

                   ∂ϕ n+1 ∂ξ = P n+1 int -P ext h ∂d n+1 2 ∂ξ 1 + ∂d n+1 2 ∂ξ 2 = 0 ∂ 2 d n+1 2 ∂ξ 2 = 1 (ϕ n+1 + µ s )   ρ s d n+1 2 -2d n 2 + d n-1 2 ∆t 2 - P n+1 int -P ext h 1 + ∂d n+1 2 ∂ξ 2   (33)
For each time-step, one has the following boundary conditions of the structural displacement d 2

d 2 (0) = 0, d 2 (L 0 ) = 0 ( 34
)
where L 0 is the length of the elastic part of the wall. In order to transform the given boundary value problem to an initial value one, the missing initial value ∂d2 ∂ξ (0) = a is defined. Thus, the initial values are then

d 2 (0) = 0, ∂d 2 ∂ξ (0) = a (35)
The missing initial condition is obtained by finding the root of the function F (a) defined as the value of d 2 at the second boundary using the initial value a

F (a) = d 2 (L 0 ; a) (36)
The Newton iteration is used to solve the root problem and the obtained value of a is then the corresponding initial value problem which then is also a solution of the boundary value problem.

Coupling with the flow equations

The Finite Difference Method (FDM) is used to solve the flow equations in a two-dimensional domain and the main difficulty is to apply the method to a complex geometry and with a moving boundary. To overcome this difficulty, the physical domain is mapped onto a rectangular one using the following coordinate transformation

τ = t, x * = x, y * =      y R 0 + w u for y ≥ 0 y R 0 -w l for y < 0 (37)
where w u and w l are, respectively, the displacement w 2 of the upper and lower elastic walls. The derivative rule is applied to compute all the derivatives involved in the Navier-Stokes equations, namely, for y ≥ 0

               ∂ ∂t = ∂ ∂τ - ∂w u (x * , τ ) ∂τ [R 0 + w u (x * , τ )] -1 y * ∂ ∂y * ∂ ∂x = ∂ ∂x * - ∂w u (x * , τ ) ∂x * [R 0 + w u (x * , τ )] -1 y * ∂ ∂y * ∂ ∂y = [R 0 + w u (x * , τ )] -1 ∂ ∂y * (38)
Similar expressions for the lower wall (y < 0) are used. The derivative rule can also be written in an analogous manner. Figure 3 shows an example of a mesh for the flow domain and for the rectangular one after using the given coordinate transformation.

The expressions of the new derivatives in the mapped domain are approximated using centered difference discretization, which is a second order ORIGINAL DOMAIN MAPPED DOMAIN Fig. 3: Meshgrids of the original and mapped domains approximation, for both first and second order derivatives and with respect to x and y, thus

∂u ∂x ≃ u i+1,j -u i-1,j 2∆x (39)
∂u ∂y ≃ u i,j+1 -u i,j-1 2∆y (40)
∂ 2 u ∂x 2 ≃ u i+1,j -2u i,j + u i-1,j ∆x 2 (41)
∂ 2 u ∂y 2 ≃ u i,j+1 -2u i,j + u i,j-1 ∆y 2 (42)
A semi-implicit pressure-correction scheme (splitting method) is used to solve the two dimensional Navier-Stokes equations for incompressible fluid flow. A forward finite differential scheme is used to evaluate the time derivative. The predictor step is first done by solving the momentum equation without the pressure gradient term for a provisional velocity based on the previous outer loop one, namely,

u * -u n ∆t = -u n .∇u n+1 + Re -1 ∆u n+1 (43)
where u * is the predictor velocity which is not divergence free. The viscosity term is treated implicitly while the convection term is treated semi-implicitly. The linear and nonlinear terms are evaluated using the second order central finite difference scheme. The second step is to correct the predictor velocity u * , such that the finally obtained velocity is divergence free. In order to maintain the Navier-Stokes equations as the basic flow equations, the pressure has to satisfy the following equation,

u n+1 -u * ∆t = -∇p n+1 (44)
The flow velocity is sought to be divergence free, i.e. ∇.u = 0, therefore, the pressure is obtained by solving the following Poisson equation

∇.u * = ∆t∆p n+1 (45)
This Poisson equation is solved with the following boundary conditions The divergence free flow velocity can then be obtained from Eq. ( 44). Finally, the procedures followed to obtain the solution for the considered coupled fluid-structure models are presented in the following algorithm.

       ∂p ∂x = 0

Fluid-Structure Interaction Algorithm

Numerical Results

In this section, a series of numerical tests were proposed to validate the structure and fluid solvers with a given exact solution for the structural model and classical blood flow problems in rigid arteries with the presence of different types of arterial aneurysms. The results of an FSI problem is then presented by coupling the wall model with the fluid equations to study the blood behavior in elastic arteries with the presence of asymmetrical aneurysms in both artery walls using the proposed nonlinear structure model.

Validation of structure solver

In order to validate the numerical scheme used for the structure problem, an exact solution is introduced for the wall motion by adding new source terms f 1 and f 2 in the structure equation, Eq. ( 25), as follows

                         ϕ = -p + λ s 2 + µ s ∂d 2 ∂ξ 2 ∂ϕ ∂ξ = P int -P ext h ∂d2 ∂ξ 1 + ∂d2 ∂ξ 2 + f 1 ∂ 2 d 2 ∂ξ 2 = 1 ϕ + µ s   ρ s ∂ 2 d 2 ∂t 2 - P int -P ext h 1 + ∂d 2 ∂ξ 2   + f 2 (47)
The following structural components are considered

d 2 (t, ξ) = sin(πξ) t 2 2 (48) ϕ(t, ξ) = -p + λ s 2 + µ s ∂d 2 ∂ξ 2 (49)
Substituting the exact solutions in [START_REF] Bakhti | Blood flow simulation through two-dimensional complex stenosed arteries using viscoelastic oldroyd-b fluid[END_REF] and with simple calculations, the following source terms are obtained

f 1 = P int -P ext h t 2 π cos(πξ) t 4 π 2 cos(πξ) 2 + 4 -( t 4 π 3 4 sin(2πξ)(λ s /2 + µ s )) (50)
f 2 = -1 µ s -p + t 4 π 2 4 cos(πξ) 2 ( λs 2 + µ s ) ρ s sin(πξ) + P int -P ext h t 4 π 2 4 cos(πξ) 2 + 1 - t 2 π 2 2 sin(πξ) (51)
The convergence of the proposed numerical scheme is tested for the structural problem by computing the resulting error between the exact and numerically obtained solutions. The structure length is L = 1, the space discritization steps vary between ∆ξ = {0.1, 0.05, 0.025, 0.0125} and the maximum time is T = 1 while the time discitization step is a fixed value ∆t = 0.02. The interior pressure is given by P int = 1 + 0.1 cos(ξ). Table 1 presents the physical parameters used in the validation tests of the structure solver.

Figure 4 represents the comparison between the given exact solutions of d 2 and ϕ and the obtained numerical ones. It is clearly seen that the numerical solution converges to the exact one and is almost the same for the given physical parameters and the given values of the time step. Table 2 presents the resulting errors and rates of convergence of the numerical scheme for the structure model. The rate of convergence is calculated using rate = log 2 ∥u h -u∥2 ∥u h/2 -u∥2 , where u h is the numerical solution using the step h and u is the exact solution. The obtained results show that the numerical scheme is of the rate 1.5 in space for both d 2 and ϕ.

Validation of fluid solver

For the fluid solver validation, multiple computational experiments and tests in order to validate the proposed numerical scheme are presented. This is used Table 2: Error analysis and the convergence rates of the structure solver for various ∆ξ. Errors are calculated using the Euclidean norm and for t = 0.5

d 2 ϕ ∆ξ ∥ d2 -d 2 ∥ 2 rate ∥ φ -ϕ∥ 2 rate 2 -1 3.0376531e -05 - 2.7657640e -01 - 2 -2
1.0065493e -05 1.5935392 8.8986147e -02 1.6360254 2 -3 3.4993314e -06 1.5242667 3.0507016e -02 1.5444396 2 to examine the effect of various fluid variables on the physiological parameters such as fluid velocity and pressure. The obtained results are discussed and presented for two separate cases. In the first case considered below (Case 1), the fluid behavior in a rigid structure is studied with the existence of one curvature in a wall side of the tube (Figure 5) or multiple symmetrical or non-symmetrical curvatures at both sides of the tube (Figure 6), with the assumption of a steady fluid velocity at the inlet of the flow domain. Several numerical runs have been made and the obtained results are presented in three sub-cases (1a, 1b and 1c). The effects of fluid and structural parameters on the fluid flow are investigated.

In the second part (Case 2), the obtained results are presented and discussed for the case of rigid structure walls with asymmetrical curvatures at both sides of the tube (Figure 6b) and with the consideration of a pulsatile velocity at the inlet of the channel. The study focuses on the effect of the inlet pulse on the velocity field and fluid pressure in the flow domain.

Case 1 : Rigid wall with steady inlet

A simple fluid flow problem is considered in a tube with the presence of a curvature in the tube wall. The wall geometry functional for the case of a steady structure can be mathematically expressed, with modified inputs from [START_REF] Bakhti | Pulsatile blood flow in constricted tapered artery using a variable-order fractional Oldroyd-B model[END_REF], as follows

R(x) = R 0 1 + 4ϵ L 2 0 L 0 (x -d) -(x -d) 2 , d ≤ x ≤ d + L 0 R 0 , otherwise (52)
where R 0 is the radius of the normal tube, L 0 is the curvature length, d indicates the location of the curvature and ϵ is the maximum height of the curvature outward the tube. This function is used for both upper and lower sides of the wall with various values of the used parameters. For the sake of clarity, this case is divided into three sub-cases in order to clarify the effects of different geometrical and fluid parameters on the flow field.

-Case 1a : Single wall curvature In order to analyze the effect of the wall curvature on the fluid flow, a series of numerical runs were done with a rigid wall and a steady inlet flow. The wall geometry with the presence of one side deformed wall is shown in Figure 5. Therefore, on the rigid wall, the no-slip condition is assumed and the inlet function is defined without the unsteady part as follows

u in = y R 1 - y R (53)
The geometric and mesh parameters, used in Case 1a, are given in Table 3. Numerical results were obtained using a regular grid, containing 51 × 51 fluid and pressure nodes. The time step used in the numerical tests is ∆t = 0.05, which is well defined to assure the accuracy needed due to the splitting error. It is well-known that classical splitting schemes usually require smaller time step. The splitting method is used by separating the pressure and velocity differential operators and the main benefits of this scheme are the simple implementation and low computational costs.

On the other hand, the structural and fluid parameters are given in Table 4, where the exact curvature location, length and the maximum height that its can reach are given. Those values will be changed from case to case for a better understanding of the used model as well as the structural effect on the flow.

Figure 7 represents the fluid streamlines and the pressure contours in the steady regime for the given values in Table 4. It is observed from this figure 4. It is observed that the fluid streamlines begin as a uniform flow and rolls start to develop over time in the curved part of the tube and keep expanding and taking more space till it fill the full curved zone of the tube due to the increase of the fluid velocity until its maximum and the fluid flow reaches a steady regime. In Figure 12, the pressure evolution is presented over time, the pressure contours obtained in the same time interval as streamlines; t ∈ {5, 10, 14.7, 25}, and for the same fixed values of the structural and fluid parameters given in Table 4. It is observed from these results that a higher pressure occurs near the middle of the curvature and continues to move over time around the curvature heading to the end of the deformed wall. When the flow reaches the steady regime, the pressure takes a higher value near the end of the deformed wall and slightly decreases around the wall curvature and that is due to the decrease of the tube diameter till it reaches the normal one of the non-deformed part of the wall.

To understand more the flow behavior in the given flow domain, the effect of various geometrical and fluid parameters on the flow is studied, namely, the Reynolds number (Re), the curvature length (L 0 ) and the curvature max- 4. It can be clearly observed that for an increase in L 0 , the rolls become larger and take more time to reach the steady regime due to the increase in the deformed zone. Finally, the influence of the curvature maximum height on the flow at t = 10 is analyzed. The resulting streamlines for various values of ϵ, ϵ ∈ {0.5, 1, 1.5}, are presented in Figure 11. It is demonstrated that small values of ϵ result in a uniform flow while an increase in the value of ϵ results on a large number of rolls. This is due to the large zone for the rolls to take place and to expand.

Case 1b : Symmetrical two wall curvatures In this case, the behavior of the flow with the presence of two symmetrical curvatures on the two sides of the channel walls is investigated. The considered wall geometry with the presence of two symmetrical curvatures is shown in Fig 6a . Several numerical tests have been made with a fixed wall and steady inlet flow using the same curve function defined in Eq. ( 52) for both upper and lower sides of the wall and with the same structural values for both sides. The no-slip condition is assumed on the rigid wall and the inlet function, without the unsteady part, is defined as follows

u in = 1 - y R 2 (54)
The geometric and mesh parameters, used in Case 1b, are the same ones given in Table 3. The curvature is supposed to be symmetrical for both sides of the wall and the structural and fluid parameters are given in Table 4. The results were obtained using the same mesh parameters and time step as in Case 1a.

Figure 13 shows the streamlined profiles and pressure contours in the steady regime for the given values of the geometrical, structural and fluid parameters values in Tables 34. It can be observed from Figure 13a, as in the first case, that the flow is uniform in the inlet and middle of the tube, while near the curved part of each side of the wall the flow forms symmetrical rolls on both sides of the channel covering the whole deformed zones. Similarly, Figure 13b shows a higher pressure occurring near the ends of the symmetrical curvatures that slightly decreases around the deformed part of the wall toward the inlet of the tube. Figures 14-15 represent the streamline profiles and pressure contours evolution over time for fixed values of structural and fluid parameters given in Table 4. Figure 14 presents the streamlines for the time values t ∈ {5, 9, 13, 25}. It is observed that the fluid begins to form rolls on both curved sides of the walls in a symmetrical way and starts to develop over time in that zone of the tube and keeps taking more space till filling the full deformed zone of the tube. That is due to the increase of the fluid velocity until the flow reaches the steady regime.

The pressure evolution over time in the case of two symmetrical curvatures on both sides of the channel is presented in Figure 15. The pressure contours are obtained for t ∈ {5, 10, 14.7, 25} and for the same fixed values of the structural and fluid parameters given in Table 4. A higher pressure occurs near the middle of both curvatures and continues to move around them in a symmetrical way heading to both ends of the deformed walls till it reaches a steady regime. The pressure then takes a higher value near the ends of the deformed walls and slightly decreases around the wall curvatures. That is due to the decrease in the tube diameter till reaching the normal one of the nondeformed parts of the walls.

While comparing these results to the one side curvature case, 1a, it is noticed that the fluid finds more space to take place in this case because of the new symmetrical deformed side. Thus, the fluid rolls occur less larger in case 1b than 1a while using the same fluid parameters given in Table 4. Also, the highest pressure value in this sub-case is lower than in the previous one.

Case 1c : Asymmetrical wall curvatures

In this final sub-case, the assumption of a steady inlet velocity is kept and the flow behavior in a channel with two non-symmetrical curvatures of the wall to observe the effect of the asymmetry on the fluid flow is investigated. The considered flow domain with asymmetrical curvatures is sketched in Figure 6b. As previously done, numerical tests with the consideration of a fixed wall and steady inlet flow are elaborated using the same curve function defined in Eq. ( 52), for both upper and lower sides of the wall. Different structural parameters values are used to consider better the non-symmetry of the walls geometries. Again, on the rigid wall, the no-slip condition is taken and the inlet function u in is used as in the previous sub-case, 1b, given in Eq. (54). The geometric and mesh parameters, used in this sub-case, are the same ones given in Table 3. As the curvatures are supposed to be non-symmetric, the structural parameters used for the upper side of the wall and fluid parameters are taken from Table 4 while the lower side structural parameters are defined in Table 5. The numerical results were obtained using the same mesh parameters and time step as in the other sub-cases.

Figure 16 shows the steady solution obtained for the streamlined profiles and pressure contours in the steady regime for the geometrical, structural and fluid parameters values given in Tables 345. It is observed from Figure 16a that the flow is uniform in the inlet and quasi-uniform in the middle of the tube. Near the curved part of each side of the wall, the flow forms non-symmetrical rolls in the two sides of the channel covering the whole deformed zone on the upper side, while on the lower side the rolls take place in the right part of the deformed zone. The same observation goes to the pressure contours plotted in Figure 16b that shows a higher pressure near the end of both non-symmetrical curvatures that slightly decreases around the deformed part of the wall toward the inlet of the tube. The pressure is lower at the end of the lower curvature than in the upper one.

With respect to time, the streamline profiles and pressure contours evolution for fixed values of structural and fluid parameters given in Tables 4-5 are presented in Figures 1718. Figure 17 represents the streamlines for the time values t ∈ {5, 9, 13, 25} showing that the fluid starts to form rolls in the upper curvature of the tube and takes some time for the lower curvature, this is due to the distance between the two curvatures and the asymmetrical pro- prieties until the steady regime is reached where the rolls fill the full deformed zone of the upper side of tube and doesn't fill all the space in the lower side that is due to the asymmetrical positions of the curvatures and to the fact that the fluid found more space at the starting position of the lower curvature. Figure 18 presents the pressure evolution over time in the case of two nonsymmetrical curvatures on the two sides of the channel. The same proprieties as in Figure 17 are used. New observations are made that differ from the ones obtained in other sub-cases. Higher pressure occurs near the middle of the curvatures in a non-symmetrical way and continue to move around the curvatures in a non-symmetrical way heading to both ends of the deformed walls till reaching the steady regime where the pressure takes a higher value near the ends of the deformed walls and slightly decreases around the wall curvatures. That is due to the decrease of the tube diameter, while it is remarked that the pressure value at the end of the upper curvature is higher than in the end of the lower curvature due to the non-symmetrical proprieties of the curvatures.

Comparison with the previous sub-cases, 1a and 1b, it can be noticed that, as in Case 1b, the two sides curvatures offer more space to the fluid. The non-symmetrical effect on the fluid flow is clearly shown in the resulting figures where the rolls are not symmetrical and have less appearance in the lower curvature than in the upper one. The pressure contours also display a non-symmetrical behavior on the non-symmetrical walls.

Case 2 : Rigid wall with pulsatile inlet

In this case, the flow behavior is studied in a fixed geometry with the assumption of a pulsatile inlet velocity profile. The asymmetric wall geometry presented in Figure 6b is considered herein. As in the previous cases, several numerical tests with the consideration of a fixed wall and pulsatile inlet flow are performed on the rigid wall, as a boundary condition, the no-slip condition is assumed. Also, the steady solution, found in Case 1, is used as an initial condition in order to better understand the pulsation effect on the flow. The pulsatile inlet function is defined as follows

u in (y, t) = u 0 (t) 1 - y R 2 (55)
where the function u 0 (t) involved in the inlet flow is a periodic function that presents the velocity amplitude variation in the x-axis direction over time, as shown in Figure 19a. For a proper understanding of the fluid flow in the proposed domain, the set of selected time values where the streamlines and pressure contours are plotted and presented in Figure 19b.

The geometric and mesh parameters used in Case 2 are also given in Table 3. The structural parameters used for the upper side of the wall and fluid parameters are taken from Table 4 while the lower side structural parameters are defined in Table 5. The numerical results were obtained using the same mesh parameters and time step as in Case 1. 345. The main obtained results are: the flow becomes periodic due to the periodic inlet velocity, an increase in the inlet velocity results in the expansion of fluid rolls in both curvatures in a non-symmetrical way, and when it decreases, the rolls become less visible and maybe removed in the lower curvature zone. Also, the flow keeps it uniform behavior in the inlet and quasi-uniform in the middle of the tube as in the previous case. The same observations are made for the pressure contours, as obtained in previous case. Higher pressure occurs near the end of the non-symmetrical curvatures and slightly decreases around the deformed part of the wall toward the inlet of the tube.

In comparison with the steady inlet case, it is clearly noticeable that the pulsation effect on the fluid flow behavior. An interesting result of the pulsatile flow is that there is a time difference between the instant where the inlet amplitude reaches its maximum value and the instant where the maximum appearance of the fluid rolls occurs.

It is important to mention that the previous results are well validated with the existing ones from the literature and other numerical techniques.

FSI model for pulsatile blood flow with elastic asymmetric aneurysm

In this subsection, some interesting results are obtained and analyzed from the simulation of the fluid-tube interaction model using the coupled numerical scheme presented in this paper. In this case, the structure walls are assumed not fully rigid and consist of an elastic part. The fluid and structure behaviors were studied for pulsatile inlet velocity and for non-symmetrical walls of the channel.

As the main application, the pulsatile blood flow in asymmetric elastic arterial aneurysm is considered. A two dimensional geometry defined by a tube with both rigid and elastic parts of the wall, as previously shown in Figure 1 is used. Based on the FSI algorithm, several numerical tests have been made to simulate the fluid-tube interaction with pulsatile inlet velocity and asymmetrical elastic parts of the tube walls. The no-slip condition is considered at the rigid part of the wall and the continuity of the velocity at the elastic part, as defined in [START_REF] Bakhti | Steady Flow of Couple-Stress Fluid in Constricted Tapered Artery: Effects of Transverse Magnetic Field, Moving Catheter, and Slip Velocity[END_REF]. The pulsatile inlet function is defined as follows

u in (y, t) = u 0 (t) 1 - y R 2 (56)
where the function u 0 (t) is the inlet velocity amplitude. For a proper understanding of the fluid flow in the proposed domain, the set of selected time values is presented in Figure 22 where the streamlines and pressure contours were taken.

The geometrical and mesh parameters used in the considered FSI problem are given in Table 6, while the structural parameters used for the upper and lower sides of the wall and fluid parameters are given in Table 7. The results were obtained using a regular grid, with 101×51 fluid and pressure nodes. The time step used in these numerical tests is ∆t = 0.1250. During the computation process, no stability problems is seen, or in the obtained solutions.

Figures 2324show the streamlines profiles and pressure contours evolution with respect to the set of time values defined in Figure 22 with fixed values of the geometrical, structural and fluid parameters values presented in Tables 67. It is clearly demonstrated that the behavior of the fluid flow is not the same as discussed in the previous cases. The flow becomes periodic due to the pulsatile inlet velocity and the clear effect on the rolls in both curvature zones is presented. On the other hand, a clear change is observed in the elastic walls displacement over time due to the fluid-tube interaction in a periodic way. Also, when the velocity increases at the inlet, the rolls become less clear and start to expand over time. As for the pressure contours, it is clearly shown that a higher pressure occurs near the end of both non-symmetrical curvatures and increases as the inlet velocity increases and starts to take its initial values over time. The pressure evolution near the tube wall plays a rule in the displacement variation of the upper and lower elastic parts.

Figure 25 shows the upper and lower maximum and minimum displacements of the elastic parts in the selected time values. It is clearly observed that there is an effect of the inlet pulsation on the wall deformation and there is a time difference between the instant where the inlet velocity reaches its maximum. The maximum wall displacement takes place on both the upper and lower sides of the tube walls. The wall displacements and flow behavior are in good agreement with previous research works in [START_REF] Hamadiche | Aneurysm no-linear unsteady dynamic: Mathematical modeling[END_REF][START_REF] Hamadiche | Numerical simulation and mathematical analysis of flow-wall interaction in the large deformation application to the dynamics of the aneurysms[END_REF][START_REF] Bakhti | Blood flow simulation through two-dimensional complex stenosed arteries using viscoelastic oldroyd-b fluid[END_REF].

With comparison to the rigid case, a different behavior of the streamlinespressure profiles due to the increase and decrease in both deformations of the elastic walls. It is observed that the rolls get infected and become less visible when there is an increase in the aneurysm displacement and that gives more space for the fluid to take place. For the pressure variation, it is clearly observed that the effect of the pressure on the wall displacement as it takes maximum form when the pressure reaches its maximum value near the wall and the minimum displacement appears when the pressure takes its lowest value near the elastic walls.

Conclusion

This paper presents a mathematical modeling and numerical simulation for fluid-tube interaction problems using large deformation theory. The application is mainly in cardiovascular systems and particularly abdominal deformed arteries and arterial aneurysms. An FSI model has been elaborated to describe incompressible fluid flow behavior in deformable geometries. This model is based on a simplified non-linear system of equations that are derived to describe the wall dynamics in a large deformation domain coupled with the non-linear Navier-Stokes equations. The pulsatile fluid flow through an elastic tube with non-symmetrical wall geometries is investigated. Efficient mathematical algorithms and numerical simulations are elaborated based on the splitting finite difference method and shooting procedure. A well-adopted semi-implicit and predictor-correction methods are elaborated. A mesh transformation is used to account for irregular domains. The coupled nonlinear fluid-structure interaction problem is numerically solved and the wall displacement and fluid velocity-pressure field are computed. Numerical experiments in three different cases have been investigated to fully understand flow behavior in rigid curved tubes and tubes with non-symmetrical elastic parts. The effects of several structural and fluid parameters on flow are studied and discussed.

The obtained results for different structural and fluid settings of the proposed model are presented. The main observations are summarized as follow -The structural and fluid solver are validated and adapted to the presented model for the elastic aneurysm motion and blood flow behavior. -For the case of rigid deformed walls and steady inlet velocity, the flow behavior changes when changing the wall curvatures parameters or with the Reynolds number. The streamline profiles, in the case of non-symmetrical walls, behave in a non-symmetrical way, as the rolls are not symmetrical and have less appearance in the lower curvature than in the upper one. The same observations are made for the pressure contours compared to other cases of one or two symmetrical curvatures. -For the case of two rigid non-symmetrical curvatures with a periodic pulsatile inlet velocity, it is clearly noticeable that the fluid flow behavior appears to be periodic. Also , there is a gap or a time difference between the instant where the inlet amplitude reaches its maximum value and the instant where the maximum appearance of the fluid rolls occurs in both curvatures zones. -For the FSI model, the case of two non-symmetric elastic curvature with a pulsatile inlet velocity, the pulsation affects the elastic walls displacements and appears to be periodic with a time gap between when the inlet amplitude reaches its maximum value and when the maximum wall displacement of the elastic walls takes place. -Also, for the case of deformed walls, it is observed that the rolls in the deformed zone become less visible due to the increase in the wall displacement and start to build up in a non-symmetrical way over time with clear appearence in the upper aneurysm than in the lower one till it reaches the initial form. -A large variation of the wall displacement with respect to the pressure values near the elastic walls is clearly demonstrated. -The results are in good agreement with previous research work.

The presented fluid-tube interaction model may be very useful in predicting and diagnosing many cardiovascular diseases such as arterial aneurysms and abnormally deformed arteries. The model can be further investigated for higher dimensions and different fluid rheology, as blood is often considered as a non-Newtonian fluid.

Fig. 1 :

 1 Fig. 1: Schematic representation of the fluid flow domain with defomable walls.

  at the inlet of the channel p = p 0 at the outlet of the channel n.∇p = 0 at the rigide and elastic wall[START_REF] Hamadiche | Aneurysm no-linear unsteady dynamic: Mathematical modeling[END_REF]

Fig. 4 :

 4 Fig. 4: Numerical and exact solution of the structure model for different values of t

Fig. 5 :Fig. 6 :

 56 Fig. 5: 2D view of the geometry structure with a curvature at one side of the tube

Table 4 :Fig. 7 :Figures 8 -

 478 Fig. 7: Streamlines profile (a) and pressure contours (b) obtained in the steady regime for the case 1a

Fig. 8 :Fig. 9 : 5 Fig. 10 :Fig. 11 :Fig. 12 :

 895101112 Fig. 8: Flow evolution corresponding to the Case 1a at different time values; (a) t = 5, (b) t = 10, (c) t = 14.7, (d) t = 25.

Fig. 13 :

 13 Fig. 13: Streamlines profile (a) and pressure contours (b) obtained in the steady regime for the case 1b

Fig. 14 :

 14 Fig. 14: Flow evolution over time for the Case 1b : (a) t = 5, (b) t = 9, (c) t = 13, (d) t = 25.

Fig. 15 :

 15 Fig. 15: Pressure contours over time for the Case 1b : (a) t = 5, (b) t = 9, (c) t = 13, (d) t = 25.

Table 5 :Fig. 16 :

 516 Fig. 16: Streamlines profiles (a) and pressure contours (b) obtained in the steady regime for asymmetrical curvatures

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Flow evolution over time for Case 1c : (a) t = 5, (b) t = 9, (c) t = 13, (d) t = 25.

Fig. 19 :

 19 Fig. 19: Amplitude of the inlet velocity as a function of time (a) with zoom in of one single pulse with the set of time samples (b)

Fig. 20 :

 20 Fig. 20: Flow evolution over time for the Case 2 : (a) t = 18, (b) t = 20, (c) t = 22, (d) t = 24, (e) t = 28, (f) t = 28.

Fig. 21 :

 21 Fig. 21: Pressure contours over time for the Case 2 : (a) t = 18, (b) t = 20, (c) t = 22, (d) t = 24, (e) t = 28, (f) t = 28.

Fig. 22 :

 22 Fig. 22: The set of time samples for the Case 3

Fig. 23 :

 23 Fig. 23: Flow evolution over time for the FSI problem : (a) t = 18, (b) t = 20, (c) t = 22, (d) t = 24, (e) t = 26, (f) t = 31.

Fig. 24 :

 24 Fig. 24: Pressure contours over time for the FSI problem : (a) t = 18, (b) t = 20, (c) t = 22, (d) t = 24, (e) t = 26, (f) t = 31.

Fig. 25 :

 25 Fig. 25: Maximum and minimum wall displacements during a single pulse of the upper (a) and lower (b) aneurysm sides with respect to the initial curvature

Table 1 :

 1 Model parameters used in section 4.1.

	Parameters Values Parameters Values
	L			1		T			1
	ρs			0.1		λs			1000
	µs			2000	h			0.1
	p			10		Pext		0.1
	0.3
										t = 0.25
	0.25									t = 0.5 t = 0.75
										Exact Sol.
	0.2
	0.15
	0.1
	0.05
	0
	-0.05
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	2500
										t = 0.25
										t = 0.5
	2000									t = 0.75
										Exact Sol.
	1500
	1000
	500
	0
	-500
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Table 3 :

 3 Geometric and mesh parameters used in the Cases 1 and 2.

	Parameters	Values Parameters	Values
	Initial time t 0	0	Final time T	50
	Time step ∆t	0.05	Normal Radius R0	0.5
	Length of the channel L 10	Number of nodes (N x × N y) 2601
	x-Space step ∆x	0.2	y-Space step ∆y	0.02

Table 6 :

 6 Geometric parameters used in case 1.

	Parameters	Values Parameters	Values
	Initial time t 0	0	Final time T	50
	Time step ∆t	0.1250	Normal Radius R0	0.5
	Length of the channel L 10	Number of nodes (N x × N y) 5151
	x-Space step ∆x	0.1	y-Space step ∆y	0.02

Table 7 :

 7 Fluid & structural parameters for both aneurysm sides.

	Parameters	Values Parameters	Values
	Reynolds number Re	200	h	0.05
	µs	1000	λs	1000
	p 0	15	Pext	5
	Upper Aneurysm location dup	1	Lower Aneurysm location d down	2.5
	Upper Aneurysm length L 0,up	4	Lower Aneurysm length L 0,down

Acknowledgements The authors would like to acknowledge the financial support of the CNRST and the Moroccan Ministry of Higher Education and Scientific Research with the project PPR2/06/2016. The authors acknowledge with thanks the technical and financial supports from PPR2.


Emmanuel Vazquez
email: emmanuel.vazquez@supelec.fr

Julien Villemonteix

Maryan Sidorkiewicz

Éric Walter

Global optimization based on

Keywords: Global optimization, computer simulations, kriging, Gaussian process, noisy evaluations

published or not. The documents may come


Regis Audran

Mogo Wafo

Thomas Messi Nguelé
email: thomas.messi@facsciences-uy1.cm

Xaviera Youh Djam

Clustering-based Graph Numbering using Execution Traces for Cache Misses Reduction in Graph Analysis Applications

Keywords: Graph analysis, Execution trace, Machine learning

Social graph analysis is generally based on a local exploration of the underlying graph. That is, the analysis of a node of the graph is often done after having analyzed nodes located in its vicinity. However, over the time, networks are bound to grow with the addition of new members, which inevitably leads to the enlargement of the corresponding graphs. At this level we therefore have a problem because more the size of the graph increases, more the execution time of graph analysis applications too. This is due to the very large number of nodes that will need to be treated. Some recent work in-faces this problem by exploiting the properties of social networks such as the community structure to renumber the nodes of the graph in order to reduce cache misses. Reducing cache misses in an application allows to reduce the execution time of this application. In this paper, we argue that combining existing graph ordering with a new numbering that exploit execution traces analysis can allow to improve cache misses reduction and hence execution time reduction. The idea is to build graph numbering using execution traces of graph analysis applications and then combine it with an existing graph numbering (such as cn-order). To build this new ordering, we define a new distance and then used it to analyse execution traces with well known clustering algorithms K-means (for Kmeans-order) and hierarchical clustering ( for cl-hier-order). Experiments on a user machine (dual-core) and four cores of Grid'5000 node (Neowise) show that this combination improves slightly existing graph ordering (cn-order, numbaco, rabbit and gorder) in almost all the cases (the two cores of dual-core, all the four cores of neowise), with PageRank graph application and astro-ph dataset. For example, on neowise with one thread and Astro-ph dataset, the best performance is given with the combination kmeans-order_cn-order which allows to reduce by 42.59% the cache misses (compared to the second numbaco with 40.79%) and therefore by 7.27 % the time of execution (compared to 6.89% for the second numbaco).

I INTRODUCTION

Actors and their interconnections in social networks are modeled with graphs where nodes are actors and links are their interconnections. The advent of computers and communication networks allows to analyze data (see on [START_REF] Newman | The structure and function of complex networks[END_REF]) on a large scale and has lead to a shift from the study of individual properties of nodes in small specific graphs with tens or hundreds of nodes, to the analysis of macroscopic and statistical properties of large graphs also called complex networks, consisting of millions and even billions of nodes [START_REF] Newman | The structure and function of complex networks[END_REF].

Social networks welcome new members every day. This contributes to enlarging the size of the corresponding graphs. It should be noted that the larger the graph, the longer the execution of the analysis applications takes time [START_REF] Arai | Rabbit order: Justin-time parallel reordering for fast graph analysis[END_REF]. This usually results from a large number of [START_REF] Nguélé | Social network ordering based on communities to reduce cache misses[END_REF] cache misses, themselves caused by many memory accesses in search of certain nodes during the execution of the graph analysis application.

Graph application consists in analyzing various nodes; process a node usually involves the analysis of the other ones located in its vicinity. Reducing cache misses in an application allows to reduce the execution time of this application. In order to reduce cache misses, the main idea is therefore to renumber the graph such away that the nodes likely to be processed together become close in the memory. Recent numbering algorithms guided by this idea have been proposed such as Cn-order [START_REF] Nguélé | Using Complex-Network Properties for Efficient Graph Analysis[END_REF][START_REF] Nguélé | DSL for Social Network Analysis On Multicore Architecture[END_REF][START_REF] Nguélé | Applying Data Structure Succinctness to Graph Numbering For Efficient Graph Analysis[END_REF], Gorder [START_REF] Wei | Speedup graph processing by graph ordering[END_REF], Rabbit order [START_REF] Arai | Rabbit order: Justin-time parallel reordering for fast graph analysis[END_REF] and NumBaCo [START_REF] Nguélé | Social network ordering based on communities to reduce cache misses[END_REF]. They stood out for their efficiency compared to other existing algorithms ( see in [START_REF] Nguélé | Using Complex-Network Properties for Efficient Graph Analysis[END_REF][START_REF] Nguélé | DSL for Social Network Analysis On Multicore Architecture[END_REF][START_REF] Nguélé | Applying Data Structure Succinctness to Graph Numbering For Efficient Graph Analysis[END_REF]). In this paper we use the execution traces analysis to build these new numbering, and we combine them with one of the recent numbering algorithms in order to have better results.

The remainder of this paper is structured as follows. Section II presents recent works on graph ordering. Section III introduces our main contribution. Experimental results are shown in section IV. Section V is devoted to the synthesis of our contribution, together with the conclusion and future work.

II RELATED WORK

For several years, graph reordering algorithms have attracted the attention of researchers and this attention is constantly increasing. The most recent graph reordering include NumBaCo Messi Nguélé, Tchuente, and Méhaut [START_REF] Nguélé | Social network ordering based on communities to reduce cache misses[END_REF], Gorder Wei, Yu, Lu, and Lin [START_REF] Wei | Speedup graph processing by graph ordering[END_REF], Rabbit Order Arai, Shiokawa, Yamamuro, Onizuka, and Iwamura [START_REF] Arai | Rabbit order: Justin-time parallel reordering for fast graph analysis[END_REF] and Cn-order Nguélé, Tchuente, and Méhaut [START_REF] Nguélé | Using Complex-Network Properties for Efficient Graph Analysis[END_REF]. Table1, gives a summarize of the existing numbering methods with their advantages and their limits.

The main idea of NumbaCo is to exploit the community structure of graph nodes to improve graph analysis applications performances through cache misses reduction (see [START_REF] Nguélé | Social network ordering based on communities to reduce cache misses[END_REF]).

In the case of G-order, Wei, Yu, Lu, and Lin [START_REF] Wei | Speedup graph processing by graph ordering[END_REF]) offered an order that allows nodes in the direct neighborhood to be close in memory.

Rabbit-order proposed by Arai, Shiokawa, Yamamuro, Onizuka, and Iwamura [START_REF] Arai | Rabbit order: Justin-time parallel reordering for fast graph analysis[END_REF]) proposes a numbering method based on two approaches namely:

• -Scheduling based on the community structure of the graphs in the real world: they try to match the hierarchical communities of the graphs in the real world with the hierarchies located at the level of the caches of the central processor unit; • -An incremental aggregation in parallel: this one aggregates in an incremental way the vertices of the same community in parallel, which has the consequence of reducing the number of vertices to be processed. Cn-order proposed by Nguélé, Tchuente, and Méhaut [START_REF] Nguélé | Using Complex-Network Properties for Efficient Graph Analysis[END_REF] is based on a fusion of the advantages of previous models such as:

• grouping in memory nodes appearing frequently in direct neighborhood (based on Gorder) • grouping in memory nodes belonging to the same community or sub-community (based on Numbaco and Rabbit-order).

None of these methods doesn't exploit former executions of target application. In this paper, we argue that exploiting former execution of the target applications (for example, exploiting execution traces analysis) can allow to improve existing performances.

III NUMBERING BASED ON EXECUTION TRACES ANALYSIS

The idea here is to build a graph numbering by analyzing execution traces of graph applications with clustering algorithms. In this section, we first present execution traces in our case (section 3.1), then we present a distance (see section 3.2) that allows us to build our new numbering algorithm (see section 3.3). Figure 1 present the complete principle of numbering using execution traces analyses.

Execution traces

The execution traces allow to see the behavior of the application during its execution. Since it is a graph analysis application, we are interested to keep the way nodes are accessed during execution. To obtain these traces, we add in the graph analysis program some code that allow to print at each line of a file, a node and the list of other nodes directly encountered during the execution program.

The goal of our work is to analyze this execution traces file in other to build a numbering that will allow "closer nodes" to have "closer numbering". To carry out this work, we choose to analyze these execution traces with clustering algorithms such as Kmeans and Hierarchical clustering. Those algorithms use with a distance that measures the closeness between nodes of the execution traces file.

The next section present the approach used to compute the distance between nodes in execution traces file.

Distance between two nodes

Alex Groce [START_REF] Groce | Error explanation with distance metrics[END_REF] defines the distance between two execution traces x and y of the same program as follows: D(x, y) = n i=0 ∆(i) , where

∆(i) = 0, if val x i = val y i 1, if val x i ̸ = val y i val x
i and val y i are the value of x and y at the time i. Inspired by this, we define a metric between two nodes a and b encountered in the execution traces file.

Proximity D between two nodes. Let Γ(a) and Γ(b) the set of nodes encountered directly after a and b respectively. The proximity D between a and b is defined as follows:

D(a, b) = |Γ(a) ∪ Γ(b)| -|Γ(a) ∩ Γ(b)| = (|Γ(a)| + |Γ(b)| -|Γ(a) ∩ Γ(b)|) -|Γ(a) ∩ Γ(b)| = |Γ(a)| + |Γ(b)| -2|Γ(a) ∩ Γ(b)| Theorem III.1: D is a distance
The proximity D between two nodes as defined is a distance.

The complete proof of theorem is on the Annex (see in A).

Proposed Algorithms

In this section, we present numbering based on execution trace analysis and the combinations between these algorithms and the previous existing ones.

Kmeans-order and Cl-hier-order

Algorithm 1 : Kmeans-order Input: G = (V, E), K, F ile, D, Output: G ′′ 1: G ′′ ← Kmeans(K, G, F ile, D){ K is

IV EXPERIMENTAL EVALUATION

The experiments were done two platforms:

• A user machine with the following characteristics: 2 cores at 1.80GHz, 4 GB Ram, L3 of 2048 KB, L2 of 256 KB and L1 of 32KB. • A node of Grid'5000 platform (Neowise) with these characteristics: 10 nodes, each node has 1 CPU AMD EPYC 7642, each CPU has 48 cores, 512GB RAM, L3 of 256 MB. For this evaluation, we present results got with the well known graph analysis application, Pagerank [START_REF] Page | The PageRank citation ranking: bringing order to the web[END_REF]. We used a Posix thread implementation proposed by Nikos Katirtzis1 . This implementation uses adjacency list representation. We use Astro-ph dataset [START_REF] Leskovec | SNAP Datasets: Stanford Large Network Dataset Collection[END_REF] that has n=16,046 nodes and m=242,502 edges (that is almost 2.46 MB with the formula presented at Messi Nguélé and Méhaut [START_REF] Nguélé | Applying Data Structure Succinctness to Graph Numbering For Efficient Graph Analysis[END_REF], the graph size is G_space = 8m + 40N bytes).

In this section, we compare previously proposed graph numerations (cn-order, rabbit, gorder, numbaco) with the numeration induced by execution traces analysis in three ways:cache reference reduction (section 4.1), cache-misses reduction (section 4.2), execution time reduction (section 4.3). We do this task with two tables:

• Table 2 that compares cache references, cache misses and execution time on the user machine described above; • Table 3 that does the same but on four cores of Grid'5000 node (Neowise). Color signification is as follows: colored numbers on a column express the four best, the two first are bold (blue for the time, red for the cache misses and orange for the cache references).

Cache-References Reduction

Cache references correspond to the number of time the program (Pagerank in this case) looks for data at L3 cache memory. So when the data required by the processor is not present at L1 and L2 cache memory, it is looked at L3 cache memory and this causes a cache reference. This means that, the number of cache references increases with the increasing of L1 and L2 cache misses.

In Table 2, we can see that:

• With one thread, the two combinations (Kmean-order_cn-order) and (cl-hier-order_cnorder) produce the best cache-references reduction 50.55% and 50.68% respectively. Cnorder is the third with 50.28%.

• With two threads, the combination (cl-hier-order_cn-order) produces the best reduction with 46.68 %, Cn-order is the second with 46.33% and the combination (Kmean-order_cnorder) is the third with 46.04%. In Table 3, we can remark that:

• With one thread, the combination (cn-order_Kmean-order) produces the best cache-references reduction 39%, Cn-order is the second with 37.95% and Gorder is the third with 35.82%. • With two threads, the combinations (cn-order_Kmean-order, Kmean-order_cn-order, clhier-order_cn-order) produce the best reduction with 36.05%, 34.82% and 34.43% respectively. • With three threads, the combinations (cl-hier-order_cn-order, Kmean-order_cn-order, cn-order_Kmean-order) produce the best reduction with 37.25%, 37.07% and 36.45% respectively. • With four threads, the combination (cl-hier-order_cn-order) produces the best reduction with 33.56%, Cn-order is the second with 33.34% and the combination (Kmean-order_cnorder) is the third with 32.95%. According to these observations, we can say that the combination of cn-order with the numbering induces by execution traces analysis (kmeans-order and cl-hier-order) allows to reduce cache references more than existing graph ordering heuristics. Therefore, using execution traces analysis allows to reduce cache misses on L1 and L2 cache memory.

Cache-Misses Reduction

Cache misses correspond to the number of time the program looks for data at L3 cache memory and doesn't find it. When the data required by the processor is not present at L3 cache memory, it is looked at central memory and this causes a cache miss.

In Table 2, we can see that:

• With one thread, the combination (Kmean-order_cn-order) produces the best cache-misses reduction 20.62%, Cn-order is the second with 19.85% and the combination (cl-hier-order_cn-order) is the third with 16.77%. • With two threads, cn-order is the first with 20.08%, the combination (cl-hier-order_cnorder) is the second with 18.41% and the combination (Kmean-order_cn-order) is the third with 17.99%. In Table 3, we remark that:

• With one thread, the combinations (cl-hier-order_cn-order) and (Kmean-order_cn-order) produces the best cache-misses reduction with respectively 43.38% and 42.59%. Num-BaCo is the third with 40.79%. • With two threads, the combinations (Kmean-order_cn-order, cl-hier-order_cn-order) produce the best reduction with 39.00% and 38.98% respectively, and cn-order is the third with 38.88% . • With three threads, the combination (cl-hier-order_cn-order) produces the best reduction with 41.27%, cn-order is second with 40.33% and the third is the combination (Kmean-order_cn-order,) with 39.45% . • With four threads, the combinations (Kmean-order_cn-order , cl-hier-order_cn-order) produces the best reduction with 39.80% and 39.73% respectively. NumBaCo is the third with 38.80%.

These observations show that the combination of cn-order with the numbering induces by execution traces analysis (kmeans-order and cl-hier-order) is almost always among the best in terms of cache misses reduction compare to existing graph ordering heuristics. Therefore, using execution traces analysis allows to reduce cache misses on L3 cache memory.

Execution Time Reduction

Execution time correspond to duration taken by the program (Pagerank in this case) to analyse the graph passed to it as parameter.

In Table 2, we can see that:

• With one thread, the combination (Kmean-order_cn-order) produces the best time reduction 6.99%, cn-order is the second with 6.76% and the combination (cl-hier-order_cnorder) is third with 6.23%. • With two threads, Rabbit, cn-order and NumBaCo produces the best reduction with 18.09%, 14.12% and 11.83% respectively. The combination (cn-order_cl-hier-order) is fourth with 9.05%. In Table 3, we can remark that:

• With one thread, the combination (Kmean-order_cn-order) produces the best time reduction with 7.27%, NumBaCo is the second with 6.89% and cn-order is the third with 6.61%. • With two threads, the combination (cn-order_cl-hier-order) produces the best reduction with 24.27%, cl-hier-order is second with 22.18%, Rabbit-order is the third with 19.57%. • With three threads, Rabbit-order produces the best reduction with 19.03%, NumBaCo is second with 11.59% and the combination (cl-hier-order_cn-order) is third with 10.68%. • With four threads, Cn-order produces the best reduction with 12.90%, NumBaCo is the second with 12.14% and the combination (Kmean-order_cn-order) is the third with 9.01%. These observation show that, even with the execution time reduction, the combination of cnorder with the numbering induces by execution traces analysis (kmeans-order and cl-hier-order) is almost always among the best compare to the existing graph ordering heuristics.

V CONCLUSION

In this paper, we proposed a new numbering that exploit execution traces analysis in order to improve cache misses reduction and hence execution time reduction in graph application. To build this ordering, we defined a new distance and then we used it to analyse execution traces with well known clustering algorithms K-means (for Kmeans-order) and hierarchical clustering ( for cl-hier-order). This numbering is intended to be combined with existing graph numbering (such as cn-order). We showed on a user machine (dual-core) and four cores of Grid'5000 node (Neowise) that the gotten combination improves slightly existing graph ordering (cn-order, numbaco, rabbit and gorder) in almost all the cases with PageRank and astro-ph dataset. For example, on neowise with one thread, the best performance is given with the combination kmeans-order_cnorder which allows to reduce by 42.59% the cache misses (compared to the second numbaco with 40.79%) and therefore by 7.27 % the time of execution (compared to 6.89% for the second numbaco).

During the executions, we saw that the new numbering takes long execution time due the distance computation. So, build this numbering could be unfeasible with huge graphs. Therefore a direct perspective could be to see the way to reduce the time taken by the distance computation during clustering.

Publications

2|Γ(b)| -2|Γ(a) ∩ Γ(b)| -2|Γ(b) ∩ Γ(c)| + 2|Γ(a) ∩ Γ(c)| ≥ 0 ⇐⇒ 2[|Γ(b)| -|Γ(a) ∩ Γ(b)| -|Γ(b) ∩ Γ(c)| + |Γ(a) ∩ Γ(c)|] ≥ 0 ⇐⇒ |Γ(b)| -|Γ(a) ∩ Γ(b)| -|Γ(b) ∩ Γ(c)| + |Γ(a) ∩ Γ(c)| ≥ 0
We group those who have the negative sign on one side and we have:

|Γ(b)| + |Γ(a) ∩ Γ(c)| -|Γ(a) ∩ Γ(b)| -|Γ(b) ∩ Γ(c)| ≥ 0 ⇐⇒ |Γ(b)| + |Γ(a) ∩ Γ(c)| -[|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|] ≥ 0 So to show that D(a, b) ≤ D(a, b) + D(b, c), it suffices to show that |Γ(b)| + |Γ(a) ∩ Γ(c)| -[|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|] ≥ 0. Let us show that: |Γ(b)| + |Γ(a) ∩ Γ(c)| -[|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|] ≥ 0 |Γ(b)| + |Γ(a) ∩ Γ(c)| -[|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|] ≥ 0 ⇐⇒ |Γ(b)| + |Γ(a) ∩ Γ(c)| ≥ [|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|]
With this inequality, the factors that can have different values are: i

0 ) |Γ(a) ∩ Γ(c)| ≤ |Γ(a)| If |Γ(a)| < |Γ(c)| Or ≤ |Γ(c)| If |Γ(c)| < |Γ(a)| ii 0 ) |Γ(a) ∩ Γ(b)| ≤ |Γ(a)| If |Γ(a)| < |Γ(b)| Or ≤ |Γ(b)| If |Γ(b)| < |Γ(a)| iii 0 ) |Γ(b) ∩ Γ(c)| ≤ |Γ(b)| If |Γ(b)| < |Γ(c)| Or ≤ |Γ(c)| If |Γ(c)| < |Γ(b)|
For this inequality to no longer hold, the right side would have to be greater than the left side. The maximum cases for this to happen are: The result is therefore valid, which allows us to verify our 3 ieme case.

Figure 1 :

 1 Figure 1: Principle of numbering by execution traces analyses

Algorithm 2 :Algorithm 3 :Algorithm 4 :Algorithm 5 :

 2345 the desired number of cluster, G input graph, File contain the execution traces, D distance defined in section3.2} 2: Return G ′′ Algorithm 2 presents clustering based on hierarchical clustering. Cl-hier-order Input: G = (V, E), F ile, D Output: G ′′ 1: G ′′ ← HierClus(G, F ile, D){ G input graph, File contain the execution traces, D distance defined in section3.2} 2: Return G ′′ 3.3.2 Combination between Cn-order with Kmeans-order and Cl-hier-order Algorithm 3 present the combination kmeans-order_cn-order. Kmeans-order_cn-order Input: G = (V, E), K, F ile, D Output: G ′′ 1: G ′ ← Kmeans-order(K, G, F ile, D){ K is the desired number of cluster, G input graph, File contain the execution traces, D distance defined in section3.2} 2: G ′′ ← existing -numbering(G ′ ){ Renumbering graph G' resulting of kmeans-order } 3: Return G ′′ Algorithm 4 present the combination Cl-hier-order_cn-order. Cl-hier-order_cn-order Input: G = (V, E), F ile, D Output: G ′′ 1: G ′ ← HierClus(G, F ile, D){ G input graph, File contain the execution traces, D distance defined in section3.2} 2: G ′′ ← existing -numbering(G ′ ){use cn-order for Renumbering graph G' resulting of Cl_hier-order } 3: Return G ′′ Algorithm 5 present the combination Cn-order_kmeans-order. Cn-order_kmeans-order Input: G = (V, E), K, F ile, D Output: G ′′ 1: G ′ ← existing -numbering(G){use cn-order for numbering graph G } 2: G ′′ ← Kmeans(K, G, F ile, D){ using Kmeans-order for renumbering G'} 3: Return G ′′ Algorithm 6 present the combination Cn-order_Cl-hier-order. Algorithm 6 : Cn-order_Cl-hier-order Input: G = (V, E), F ile, D Output: G ′′ 1: G ′ ← existing -numbering(G){use cn-order for numbering graph G } 2: G ′′ ← HierClus(G, F ile, D){ using cl_hier-order for renumbering G'} 3: Return G ′′

1 st

 1 Case : |Γ(a) ∩ Γ(b)| = |Γ(a)| et |Γ(b) ∩ Γ(c)| = |Γ(b)| 2 nd Case: |Γ(a) ∩ Γ(b)| = |Γ(a)| et |Γ(b) ∩ Γ(c)| = |Γ(c)| 3 rd Case: |Γ(a) ∩ Γ(b)| = |Γ(b)| et |Γ(b) ∩ Γ(c)| = |Γ(b)| 4 th Case: |Γ(a) ∩ Γ(b)| = |Γ(b)| et |Γ(b) ∩ Γ(c)| = |Γ(c)| From these different cases, we can determine the probable maximum value of |Γ(a) ∩ Γ(c)| contained in the left part of the inequality and see if the inequality will always be respected.1.0.1 1 st case : (|Γ(a) ∩ Γ(b)| = |Γ(a)| et |Γ(b) ∩ Γ(c)| = |Γ(b)| )In this first case we have the following information:(|Γ(a)| < |Γ(b)| et |Γ(b)| < |Γ(c)|) =⇒ |Γ(a)| ≤ |Γ(b)| ≤ |Γ(c)| =⇒ |Γ(a) ∩ Γ(c)| = |Γ(a)| The inequality |Γ(b)| + |Γ(a) ∩ Γ(c)| ≥ [|Γ(a) ∩ Γ(b)| + |Γ(b) ∩ Γ(c)|] thus becomes |Γ(b)| + |Γ(a)| ≥ |Γ(a)| + |Γ(b)|And the inequality is therefore verified. 1.0.2 2 nd case : (|Γ(a) ∩ Γ(b)| = |Γ(a)| et |Γ(b) ∩ Γ(c)| = |Γ(c)|) Here |Γ(a)| < |Γ(b)| and |Γ(c)| < |Γ(b)| =⇒ |Γ(a)| ≤ |Γ(b)| ≥ |Γ(c)| In this case, we cannot directly conclude. We can however say that: ∃ q1, q2 ⊂ N / |Γ(b)| = |Γ(a)| + |q1| and |Γ(b)| = |Γ(c)| + |q2| Thereby |Γ(a)| = |Γ(b)| -|q1|et|Γ(c)| = |Γ(b)| -|q2| =⇒ |Γ(a) ∩ Γ(c)| = |Γ(a)| + |Γ(c)| -|Γ(a) ∪ Γ(b)| or |Γ(a) ∪ Γ(b)| = |Γ(b)| car |Γ(b)| > |Γ(a)| The inequality can be written: |Γ(b)| + |Γ(a)| + |Γ(c)| -|Γ(b)| ≥ |Γ(a)| + |Γ(c)| We note that the |Γ(a)| and |Γ(c)| will cancel on either side of the inequality and we will have: |Γ(b)| -|Γ(b)| ≥ 0; The |Γ(b)| will cancel out and we will therefore have a valid expression.1.0.3 3 rd case : (|Γ(a) ∩ Γ(b)| = |Γ(b)| et |Γ(b) ∩ Γ(c)| = |Γ(b)|) For the 3 ith case, we have |Γ(b)| < |Γ(a)| and |Γ(b)| < |Γ(c)| =⇒ |Γ(a)| ≥ |Γ(b)| ≤ |Γ(c)| In this case, we cannot directly conclude. However We can say that: ∃ q3, q4 ⊂ N / |Γ(a)| = |Γ(b)| + |q3| and |Γ(c)| = |Γ(b)| + |q4| . Thus, Γ(a) ∩ Γ(c) = (Γ(b) ∪ q3) ∩ (Γ(b) ∪ q4) = Γ(b) ∪ (q3 ∩ q4) Hence |Γ(a) ∩ Γ(c)| = |Γ(b) ∪ (q3 ∩ q4)| = |Γ(b)| + |q3 ∩ q4| -|Γ(b) ∩ (q3 ∩ q4)|The inequality therefore becomes: |Γ(b)| + |Γ(b)| + |q3 ∩ q4| -|Γ(b) ∩ (q3 ∩ q4)| ≥ |Γ(b)| + |Γ(b)| The |Γ(b)| will all cancel out and we will have: |q3 ∩ q4| -|Γ(b) ∩ (q3 ∩ q4)| ≥ 0 ⇐⇒ |q3 ∩ q4| ≥ |Γ(b) ∩ (q3 ∩ q4)| But by definition, |q3 ∩ q4| ≥ |Γ(b) ∩ (q3 ∩ q4)

1.0.4 4

 4 th case: (|Γ(a) ∩ Γ(b)| = |Γ(b)| et |Γ(b) ∩ Γ(c)| = |Γ(c)| ) For this last case, |Γ(c)| < |Γ(b)| and |Γ(b)| < |Γ(a)| =⇒ |Γ(a)| ≥ |Γ(b)| ≥ |Γ(c)| =⇒ |Γ(a) ∩ Γ(c)| = |Γ(c)|. The inequality therefore becomes: |Γ(b)| + |Γ(c)| ≥ |Γ(b)| + |Γ(c)| This gives us a valid result. All cases could be tested and verified. So we have the proof that D is a distance. B BIOGRAPHY Regis Audran MOGO WAFO Computer Science Student/Phd in University of yaounde 1 Thomas Messi Nguélé Computer Science Doctor/PhD in University of yaounde 1 Xaviera Youth Kimbi Computer Science Doctor/PhD in University of yaounde 1

Table 1 :

 1 Comparison of existing numbering methods

	Reference Year Techniques	Specifications	Limits
	[7]	2015 NumBaCo	Exploits community	Doesn't take into account
				structure of the graph	nodes degrees heterogeneity
				in order to renumber it	Doesn't exploit former
					executions of target applications
	[6]	2016	Gorder	Consider nodes sibling	Does not take into
				nodes	account communities
					Doesn't exploit former
					executions of target applications
	[5]	2016 Rabbit-order Based on community	Doesn't taking into account
				structure + is parallel	nodes degrees heterogeneity
					Doesn't exploit former
					executions of target applications
	[8]	2017	Cn-order	Combines advantages	Doesn't exploit former
				of the others	executions of target applications

Table 2 :

 2 Graph Ordering Comparison (Pagerank, Astro-ph, dual-core machine)

	Heuristic	Time (ms)	Cache misses	Cache references
	Astro-ph 1th	607.29475	2832481.5	16569884.75
	Gorder 1th	570.486 (6.06%)	2389119.25 (15.65%) 9563592.75 (42.28%)
	NumBaCo 1th	571.9995 (5.81%)	2423346.5 (14.44%)	8487259.75 (48.78%)
	Rabbit 1th	575.08125 (5.30%)	2536938.75 (10.43%) 9233646.25 (44.27%)
	Cn-o 1th	566.20725 (6.76%)	2270193.5 (19.85%)	8237936 (50.28%)
	k-means-o 1th	618.528 (-1.85%)	3044132.75 (-7.47%) 17330640.75 (-4.59%)
	cl-h-o 1th	594.1935 (2.16%)	2522649.5 (10.94%)	15282698.25 (7.77%)
	Cn-o_k-means-o 1th	572.231(5.77%)	2398583.25 (15.32%) 8488329.75 (48.77%)
	Cn-o_cl-h-o 1 th	585.329 (3.62%)	2420409.25 (14.55%) 14831642.25 (10.49%)
	k-means-o_Cn-o 1 th	564.867 (6.99%)	2248500 (20.62%)	8193249.75 (50.55%)
	cl-h-o_Cn-o 1th	569.471 (6.23%)	2357463.75 (16.77%)	8172304.5 (50.68%)
	Astro-ph 2th	475.46975	3455467	18489284.5
	Gorder 2th	536.18175(-12.76%) 2964778.75(14.20%)	11210821 (39.36%)
	NumBaCo 2th	419.18575(11.83%)	2958648.25(14.37%)	10173745.5 (44.97%)
	Rabbit 2th	402.621(18.09%)	3026377 (12.41%)	11077611.25 (40.08%)
	Cn-o 2th	408.30775(14.12%)	2761599(20.08%)	9923092 (46.33%)
	k-means-o 2th	520.92475(-9.56%)	3739659 (-8.22%)	19466981 (-5.28%)
	cl-h-o 2th	452.037(4.92%)	3007137.25(12.97%)	17509240.5 (5.30)
	Cn-o_k-means-o 2th	460.95275(3.05%)	2950742(14.60%)	10542184.25 (42.98%)
	Cn-o_cl-h-o 2 th	432.43575 (9.05%)	3000986(13.15%)	16928641 (8.44%)
	k-means-o_Cn-o 2 th	437.171 (8.05%)	2833740.5 (17.99%)	9975563.25 (46.04%)
	cl-h-o_Cn-o 2th	449.738(5.41%)	2819274.75(18.41%) 9858212.75 (46.68%)

Table 3 :

 3 Graph Ordering Comparison (Pagerank, Astro-ph, 4 thread, Neowise[grid'5000])

	Heuristic	Time (ms)	Cache misses	Cache references
	Astro-ph 1th	266.905	10276724.75	42293956
	Gorder 1th	250.212 (6.25%)	6761279 (34.20%)	27140765.5 (35.82%)
	NumBaCo 1th	248.50975 (6.89%)	6084717.75 (40.79%)	29021769.5(31.38%)
	Rabbit 1th	253.89025(4.87%)	7325485.25(34.57%)	27670585(34.57%)
	Cn-o 1th	249.2455 (6.61%)	6220675 (39.46%)	26241990 (37.95%))
	k-means-o 1th	269.25375(-0.87%)	10216968(0.58%)	39365789.5(6.92%)
	cl-h-o 1th	256.76975 (3.79%)	7417636(27.82%)	54311434(-28.41%)
	Cn-o_k-means-o 1th	250.4(6.18%)	6437387.75(37.35%)	25798704.75 (39.00%)
	Cn-o_cl-h-o 1 th	255.4765(4.28%)	7529694.25(26.73%)	48205558.75(-13.97%)
	k-means-o_Cn-o 1 th	247.4955 (7.27%)	5899736.5 (42.59%)	27898282 (34.03%)
	cl-h-o_Cn-o 1th	249.75525 (6.42%)	5818206.25 (43.38%)	27789259 (34.29%)
	Astro-ph 2th	184.07075	11069265	48044309.75
	Gorder 2th	209.2585(-13.68%)	7236403(34.62%)	32165345 (33.05%)
	NumBaCo 2th	158.03575(14.14%)	6797367.75(38.59%)	33022849.75(31.26%)
	Rabbit 2th	148.01475 (19.57%) 7976793.25(27.93%)	35010926.75 (27.12%)
	Cn-o 2th	164.0755(10.86%)	6765081(38.88%)	31541515 (34.35%)
	k-means-o 2th	184.31575(-0.13%)	11660819(-5.34%)	46860805.5 (2.46%)
	cl-h-o 2th	143.2345 (22.18%)	9089591.5(17.88%)	61258435 (-27.50)
	Cn-o_k-means-o 2th	165.903(9.86%)	7194631.5(35.00%)	30722091.25 (36.05%)
	Cn-o_cl-h-o 2 th	139.3905 (24.27%)	8760113.75(20.86%) 56590441.75 (-17.79%)
	k-means-o_Cn-o 2 th	162.779(11.56%)	6751559(39.00%)	31312716.75 (34.82%)
	cl-h-o_Cn-o 2th	157.9145 (14.20%)	6754098.25(38.98%)	31498630.5 (34.43%)
	Astro-ph 3th	136.47175	12759681.25	54466648.25
	Gorder 3th	173.27525(-26.23%)	7975692.5(37.49%)	35949448.5 (33.99%)
	NumBaCo 3th	120.644(11.59%)	7749122.5(39.26%)	36747591.5 (32.53%)
	Rabbit 3th	110.48925(19.03%)	9079863.5(28.83%)	40484777.75 (25.67%)
	Cn-o 3th	124.16175(9.02%)	7613301.75(40.33%)	34903471.25 (35.91%)
	k-means-o 3th	137.493(-0.74%)	12949269.75(-1.48%)	52550052.75 (3.51%)
	cl-h-o 3th	135.2245(0.91%)	9967014.5(21.88%)	61710753.75 (-13.30%)
	Cn-o_k-means-o 3th	125.599(7.96%)	8101565.5(36.50%)	34615031.75 (36.45%)
	Cn-o_cl-h-o 3 th	133.95425(1.84%)	9917300.75(22.27%)	57836840.25 (6.18%)
	k-means-o_Cn-o 3 th	129.6565(4.99%)	7724722(39.45%)	34273040.25 (37.07%)
	cl-h-o_Cn-o 3th	121.89275(10.68%)	7493254(41.27%)	34177652.75 (37.25%)
	Astro-ph 4th	266.07925	12270609.5	54432887.75
	Gorder 4th	245.762(7.63%)	7763059.75 (36.73%)	38749439 (28.81%)
	NumBaCo 4th	233.7525(12.14%)	7508885.25 (38.80%) 38274430.75 (29.68%)
	Rabbit 4th	247.38(7.02%)	8379379.5 (31.71%)	40160100.5 (26.22%)
	Cn-o 4th	231.7295(12.90%)	7408551.5 (39.62%)	36283533.5 (33.34%)
	k-means-o 4th	289.7065(-8.87%)	11881925.25 (3.16%)	52188152.5 (4.12%)
	cl-h-o 4th	305.36675(-14.76%) 9241197.75 (24.69%) 64502363.25 (18.49%)
	Cn-o_k-means-o 4th	251.41425(5.51%)	7714594.25 (37.13%) 36841011.75 (32.31%)
	Cn-o_cl-h-o 4 th	298.334(-12.12%)	8885607.5 (27.59%)	62234059 (-14.33%)
	k-means-o_Cn-o 4 th	242.103(9.01%)	7386400.25 (39.80%) 36492023.25 (32.95%)
	cl-h-o_Cn-o 4th	242.649(8.80%)	7394821.25 (39.73%) 36160439.75 (33.56%)

https://github.com/nikos912000/parallel-pagerank

A ANNEX:PROVING THAT D IS A DISTANCE Proof:. In other to show that D is a distance, we prove that it respects the four properties of a distance. Let N be the set of nodes: We also pass -2|Γ(a) ∩ Γ(c)| on the other side of the inequality in order to have 0 on one side. Thus, we have: