Directeur De

Professeur, DynFluid Jean-Christophe Robinet

Jury M Lutz

Lesshafft

Président Mme Taraneh Sayadi

Professeur M Nicolas Gourdain

Professeur, EPFL M François Gallaire

Examinateur M Grégory Dergham

Safran Docteur

Examinateur M Tech

Content, Docteur, ONERA Examinateur M Cédric

Châtillon Examinateur

M Sébastien Bourasseau

Onera Docteur

Examinateur M Châtillon

Samir Beneddine

Invité M Denis Meudon

Sipp

Scale up of efficient global stability tools in order to characterize specific turbomachineries phenomenons

Keywords:

Scaling up of global stability tools in order tocharacterize phenomena specific to turbomachinery

Valentin Fer

To cite this version:

Valentin Fer. Scaling up of global stability tools in order tocharacterize phenomena specific to turbomachinery.

I.1 Rotating stall phenomenon

The prediction of unsteady aerodynamic phenomena specific to turbomachinery is a key issue for aircraft manufacturers. In axial compressors, under the effect of a three-dimensional disturbance, a pocket of fluid with a lower rotational speed than that of the machine might stall on one or more blades. The stalls may coalesce and modify the angle of incidence seen by the upstream fluid, which is then deflected, giving rise to aerodynamic blockages that result in significant flow rate drops and consequent energy losses. The flow becomes unstable and the phenomenon, known as Understanding these mechanisms is crucial to expand operating conditions and ensure system safety. The causes of these phenomena have been extensively investigated over the last century [START_REF] Day | Stall, Surge, and 75 Years of Research[END_REF]. Greitzer provided in [START_REF] Greitzer | Surge and Rotating Stall in Axial Flow Compressors-Part I: Theoretical Compression System Model[END_REF], a theoretical criterion based on the mass flow coefficient for which when the value exceeds a critical level, the system begins to pulse at frequencies characteristic of the compressor surge. This criterion was validated using experimental data in [START_REF] Greitzer | Surge and Rotating Stall in Axial Flow Compressors-Part II: Experimental Results and Comparison With Theory[END_REF]. Courtiade & Gourdain [START_REF] Courtiade | Experimental investigation of rotating stall in a high speed multi-stage axial compressor[END_REF] conducted experiments on rotating stall on the compressor CREATE (for Compresseur de Recherche pour l'Étude des effets Aérodynamiques et Technologiques) and established a stability limit. They observed that the entrance into surge regime caused the appearance of twelve stall cells successively shrinking and growing, spinning at 85% of compressor speed. When loads increased beyond the critical level, this pseudo-stable state degenerated into a full-span cell, triggering the surge. Other studies of the low-subsonic axial compressor CME2 [START_REF] Margalida | Comparison and sensibility analysis of warning parameters for rotating stall detection in an axial compressor[END_REF] were carried out as part of the European ACONITE project [START_REF] Dazin | The ACONIT project: an innovative design approach of active flow control for surge prevention in gas turbines[END_REF], which aims to improve the stable operating range of axial compressor and reduce the surge margin.

Numerical investigations have also been documented. Hamilton [START_REF] Hamilton | Parallel Simulation of Stall Hysteresis in a Two Dimensional Rotor Blade Row[END_REF] shows that a two-dimensional simulation of unsteady incompressible Reynolds-averaged Navier-Stokes (RANS) equations can effectively model the rotating stall on a portion of an isolated rotor. The code simulated the sudden appearance of a stall cell on a 15 rotor channels configuration. Other quasi-three-dimensional approaches, like that of He [START_REF] He | Computational Study of Rotating-Stall Inception in Axial Compressors[END_REF], explores the rotor-stator interaction with a time-marching Navier-Stokes (NS) method and a mixing-length turbulence model. Various configurations were tested and the study shows that the rotating stall is strongly influenced by the stator. It also demonstrates that the introduction of a perturbation upstream from the compressor plays a role in the instability process. In [START_REF] Gourdain | Simulation numérique des phénomènes de décollement tournant dans les compresseurs axiaux[END_REF], Gourdain presents quasi-3-D simulations of the CME2 compressor. The computations were carried out over 60 revolution periods on a 6.5 × 10 5 points mesh. By establishing a "compressor characteristic", he predicts the limit cycle of the rotating stall and shows that the basic mechanism is intrinsically two-dimensional. However, these methods face difficulties in comparison with experimental data. To improve the accuracy of computational experiments and take into account three-dimensional effects and losses, three-dimensional simulations must be performed.

Analyzing the onset of rotating stall on a rotor/stator configuration, Marty et al. showed in [START_REF] Marty | Numerical Investigations on the Rotating Stall in an Axial Compressor and its Control by Flow Injection at Casing[END_REF] that a 3-D RANS simulation on a single blade channel is sufficient to predict the stability margin. However, to capture the phenomenon, unsteady RANS simulations in a computational domain that matches the engine periodicity (at least 1/10 for the CME2 compressor, resulting in 3 rotor channels and 4 stator channels) must be performed. This agrees with the calculations by Gourdain [START_REF] Gourdain | Simulation numérique des phénomènes de décollement tournant dans les compresseurs axiaux[END_REF] on the same configuration. However, above the stability threshold, such reduced computational domains don't accurately predict the flow behavior in the compressor. To better understand the phenomenon and have a simulation that is more representative above the instability onset, 3-D simulations are mandatory. Gourdain conducted such calculations in [START_REF] Gourdain | Simulation numérique des phénomènes de décollement tournant dans les compresseurs axiaux[END_REF] on a 31 × 10 6 points computational domain and observed the complete coalescence of every stall in about 13 machine rotations. These calculations were carried out on a vectorial machine in approximately 900 hours. Marty et al. [START_REF] Marty | Numerical Investigations on the Rotating Stall in an Axial Compressor and its Control by Flow Injection at Casing[END_REF] performed simulations on a full-annulus rotor-stator configuration to investigate the control of the phenomenon. The computational domain consisted of 110 × 10 6 points, and the complete coalescence was observed in 10 machine rotations, which took approximately 390 hours on ONERA's clusters. At the price of high restitution times for design phases, these studies improved the understanding of the phenomenon, but its precise origin remains unknown. To reduce the level of modeling induced by RANS computations, LES (for Large Eddy Simulation) simulations can also be used [START_REF] Gourdain | Prediction of the unsteady turbulent flow in an axial compressor stage. Part 1: Comparison of unsteady RANS and LES with experiments[END_REF]. However, these methods are about 500 times more expensive than RANS computation, and therefore too costful for industrial use.

Alternatively, the emergence of unsteady phenomena is often linked to the formation of at least one unstable mode in the flow. The exact source of these temporally self-sustained oscillations can be understood using linear stability theory. Applying global stability analysis to the flow in the CME2 compressor may help identify precursors to the rotating stall phenomenon. The various concepts of global stability theory are described below.

I.2 Linear stability theory

Linear stability theory focuses on the stability of solutions of differential equation-governed dynamic systems under the influence of small perturbations. A specific branch of stability theory is the study of the stability of fixed points, which are steady solutions of the governing equations in an equilibrium state. These states are considered stable if they return to equilibrium after being perturbed, as with a damped pendulum that returns to its original position after being moved away from its equilibrium state, or a ball that returns to its starting point in a well (see fig. I.3a). If the perturbation results in a motion with large amplitudes, the system is considered unstable, like a ball on the top of a hill (see fig. I.3b). The equilibrium may also be in a neutral state (see fig. I.3c), in which a perturbation would move the ball away from its starting position with dampened motion. Linear stability theory thus examines the behavior of small perturbations in a dynamic system. In fluid dynamics, the stability of equilibrium states have been studied for over a century. For instance, Kelvin and Helmholtz [START_REF] Thomson | XLVI. Hydrokinetic solutions and observations[END_REF][START_REF] Helmholtz | Über discontinuierliche flüssigkeits-bewegungen [on the discontinuous movements of fluids[END_REF] studied an hydrodynamic instability arising when two fluids are superimposed and one has a velocity above a critical level, giving rise to the so-called Kelvin-Helmholtz instability shown in fig. I.4a. Another example studied by Bénard and von Kàrmàn [START_REF] Bénard | 2e semestre[END_REF][START_REF] Von Karman | The mechanism of the uid and air-resistance[END_REF] is the flow past a body that, when reaching a critical Reynolds number [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF], triggers a vortex shedding, known as the von Kàrmàn vortex street, shown in fig. I.4b. Initially, instability was associated with parallel shear flows or axisymmetric jet for instance with the work of Lin [START_REF] Lin | On the stability of two-dimensional parallel flows[END_REF], Batchelor & Gill [START_REF] Batchelor | Analysis of the stability of axisymmetric jets[END_REF] or Drazin & Reid [START_REF] Drazin | Hydrodynamic stability[END_REF]. In the 1990s, the stability theory about parallel flows has been extended with nonmodal perturbation development by Schmid & Henningson in [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]. All these studies were conducted on the basis of local stability analysis that lies on the hypothesis of parallel or weakly nonparallel flows. For more complex configurations, such assumption is seldom valid and the so-called global stability analysis is required. Huerre & Monkewitz, in [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF], made the distinction between convectively unstable flow, i.e. flows that are spatially unstable and convected out of the domain; and absolutely unstable flows, i.e. flows which grow unsteadily all over the domain and are amplified, without getting convected out (by analogy with the rotating stall and surge presented lately). Absolutely unstable flows are a type of unstable global mode.

Direct global stability analysis. Let us consider a fluid governed by the so-called Navier-Stokes (NS) equations, which is a nonlinear system of PDE (for Partial Differential Equations). Considering a state vector w containing the conservative variables of the flow, its evolution can be described by the following conservative equations:

∂w ∂t = R(w), (I.1)
where R is the NS operator. The linear stability analysis aims to study the temporal evolution of a small perturbation around a base flow, which is a steady solution w b of the system such that R(w b ) = 0. This perturbation is governed by a linear equation, expressed as:

∂w ′ ∂t = J (w b )w ′ , (I.2)
where J (w b ) = ∂R/∂w| w b is the linearized NS operator, otherwise known as the Jacobian of the flow. In global stability analysis, solutions are supposed to grow exponentially in time, under the form w ′ = w(x, y, z, t)e λt , resulting in the so-called direct eigenvalue problem:

J (w b ) w = λ w. (I.
3)

The stability of the k-th direct mode w k is described by its corresponding eigenvalue λ k . The latter is composed of an amplification rate σ k = ℜ(λ k ) and an angular frequency ω k = ℑ(λ k ), while the eigenfunction w k (x, y, z) provides the spatial structure of the mode. If at least one eigenmode has a positive amplification rate, the base flow is said to be asymptotically unstable. Global modes are poorly sensitive to external perturbations and describe the asymptotical behaviour of the flow. They give rise to self-sustained oscillations, making the flow act as an oscillator [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. On the other hand, the base flow is asymptotically stable when only stable modes are present. Additionally, if the Jacobian is normal (i.e. J J † = J † J , where J † is the adjoint of J as defined below), its eigenvectors define an orthogonal modal basis. This means that expanding the perturbation w ′ onto such a basis composed only of stable modes will result in a monotonically decreasing amplitude of the perturbation. However, when the Jacobian is non-normal (i.e. does not commute with its adjoint), the modal basis is no longer orthonormal and energy growths may occur in the flow even in the absence of unstable modes. This phenomenon is known as short-term instability and can be studied using resolvent analysis, as described later.

With the improvement of computational resources, the global stability analysis has been increasingly successful over the past decades [START_REF] Theofilis | Advances in global linear instability analysis of nonparallel and threedimensional flows[END_REF]. The first documented investigation was made by Pierrehumbert and Widnall in [START_REF] Pierrehumbert | The two-and three-dimensional instabilities of a spatially periodic shear layer[END_REF], studying a spatially periodic shear layer. Theofilis reviewed, in [START_REF] Theofilis | Global Linear Instability[END_REF], the various types of global stability analysis and some methods to address the problems. When a 3-D domain is considered with two inhomogeneous directions, the analysis is referred to as "biglobal", while it is referred to as "triglobal" if each spatial direction is inhomogeneous. The first reported study of 3-D global stability analysis was carried out on the flow around a spheroid by Tezuka and Suzuki [START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF]. Paladini [START_REF] Paladini | Transonic buffet instability: From two-dimensional airfoils to three-dimensional swept wings[END_REF] also showed that this technique successfully predicted the appearance of the buffet phenomenon. In the context of turbomachinery, Xie et al. [START_REF] Xie | Stall inception analysis of a transonic compressor based on a global stability analysis method[END_REF] recently used the global stability analysis method to predict the stall inception of the transonic axial compressor NASA Rotor 37.

Adjoint global stability analysis. It is also of interest to study the spectral components of the adjoint operator J † in global stability analysis. It is defined relatively to a given scalar product ⟨•, •⟩, as ∀a, b: ⟨a, J b⟩ = ⟨J † a, b⟩. (I.4)

Eigenpairs of J † are found by solving the adjoint eigenproblem:

J † (w b ) w † = λ † w † . (I.5)
Each direct eigenpair (λ k , w k ) has a corresponding adjoint eigenpair (λ † k , w † k ) with λ † k = λ k and w † k the adjoint mode's spatial structure. In fluid mechanics, J is frequently non-normal, making it difficult to expand the perturbation over the nonorthogonal direct modal basis. The difficulty is overcome by normalizing each adjoint global mode such that ⟨ w k , w † k ⟩ = 1. The direct and adjoint modal basis ( w k , k ≥ 1) and ( w † k , k ≥ 1) form a bi-orthogonal basis verifying ⟨ w k , w † l ⟩ = δ kl 1 . The perturbation w ′ can be uniquely decomposed in the dual basis as:

w ′ = k≥1 ⟨ w † k , w ′ ⟩ w k . (I.6)
The adjoint mode provides the exact field to trigger the corresponding direct mode in the flow. Its spatial structure shows the flow regions most sensitive to external forcing or initial condition changes. The adjoint mode is useful for studying the receptivity of the direct mode and can be used to establish a control strategy to delay unsteadiness [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. On the other hand, the computation of both direct and adjoint global modes gives access to the so-called structural sensitivity introduced by Giannetti & Lucchini in [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. It is defined by the equation:

S(x 0 , y 0 , z 0 ) = ∥ w(x 0 , y 0 , z 0 )∥ 2 • ∥ w † (x 0 , y 0 , z 0 )∥ 2 (I.7)

and represents the spatial overlap of two given direct and adjoint modes. It enables the detection of the so-called wavemaker, which identifies the regions directly responsible for generating self-sustained oscillations. Citro has shown an example of adjoint modal analysis and structural sensitivity calculation on a 3-D case for laminar flow around a spherical body in [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF].

Resolvent analysis.

As mentioned earlier, in some cases, transient phenomena associated with short-term instabilities result in an increase of energy in the flow even in the absence of unstable mode [START_REF] Schmid | Analysis of Fluid Systems: Stability, Receptivity, Sensitivity: Lecture notes from the FLOW-NORDITA Summer School on Advanced Instability Methods for Complex Flows[END_REF]. This occurs in the planar Couette flow, which is globally stable [START_REF] Romanov | Stability of plane-parallel Couette flow[END_REF] but transitions to a turbulent regime around Re = 350 [START_REF] Lundbladh | Direct simulation of turbulent spots in plane couette flow[END_REF]. The understanding of such phenomena was first improved by Reddy et al. [START_REF] Reddy | Pseudospectra of the Orr-Sommerfeld Operator[END_REF] through the study of the pseudo-spectrum of the Orr-Sommerfield operator (see [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF] for further details); or by Trefethen et al. [START_REF] Threfethen | Hydrodynamic stability without eigenvalues[END_REF], who proposed studying the pseudospectra of the linearized problem. As briefly introduced in the paragraph on direct global modes, the explanation lies in the non-normal nature of the Jacobian operator. When J is normal, its eigenmodes are orthogonal and the flow dynamics is entirely governed by the direct global modes. This means that an increase of energy in the system corresponds to the existence of at least one unstable mode. However, when J is non-normal, the direct modes are nonorthogonal and the spectrum of J , which provides the asymptotic behavior of the system, is not sufficient to exhaustively study transient phenomena. Perturbations on a stable base flow may thus result in pseudo-resonances, as in the case of the Couette flow. w ′ = w 1 e λ 1 t + w 2 e λ 2 t , with σ 1 , σ 2 < 0. When the modal basis is orthogonal, the amplitude of w ′ , defined by ∥w ′ ∥ 2 = ∥ w 1 ∥ 2 + ∥ w 2 ∥ 2 , decreases monotonically in time. When the modal basis is nonorthogonal, the amplitude is defined by ∥w ′ ∥ 2 = ∥ w 1 ∥ 2 + ∥ w 2 ∥ 2 + 2⟨ w 1 , w 2 ⟩ cos( w 1 , w 2 ), where w 1 , w 2 is the angle between the two damped modes. In this case, depending on the level of non-normality of the global modes, an increase of energy may be observed before the amplitude of w ′ asymptotically decreases to zero. In this case, the flow is said to be a noise-amplifier.

The resolvent analysis, first introduced by Trefethen [START_REF] Threfethen | Hydrodynamic stability without eigenvalues[END_REF] in 1993, is suitable for studying shortterm instabilities. Let us apply a forcing term ϕ ′ to the linearized equations (I.1):

∂w ′ ∂t = J (w b )w ′ + ϕ ′ . (I.8)
Using the Fourier transforms w ′ = w(x, y, z)e iωt and ϕ ′ = ϕ(x, y, z)e iωt , the relation between the harmonic response w and forcing ϕ becomes:

w = R(ω) ϕ, (I.9)
where R(ω) = (iωI -J (w b )) -1 is the so-called resolvent operator, I the identity matrix and ω a given real frequency. The resolvent operator is a selective noise-amplifier that allows the identification of the optimal forcing of the flow, i.e. the forcing that maximizes the energy gain associated to R. This gain is defined given a certain Hermitian norm ∥ • ∥ (usually related to kinetic or acoustic energy), as η 2 = ∥ w∥ 2 /∥ ϕ∥ 2 . Finding the optimal forcing is thus equivalent to solve the following optimization problem:

η 2 = max ϕ ⟨ w, w⟩ ⟨ ϕ, ϕ⟩ , (I.10)
which in turns, is adressed computing the singular value decomposition of R [START_REF] Beneddine | Characterization of unsteady flow behavior by linear stability analysis[END_REF][START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. To do so, the following eigenvalue problem is solved:

R † R ϕ = η 2 ϕ, (I.11)
where R † is the adjoint of R. The operator R † R is Hermitian, meaning its spectrum is real, and the set of optimal forcings ( ϕ k , k ≥ 1) forms an orthonormal basis over which the forcing ϕ can be expanded ϕ = k≥1 ⟨ ϕ k , ϕ⟩ ϕ k . In a similar way, it may be shown that the set of optimal responses ( w k , k ≥ 1) forms an orthonormal basis and w can be uniquely decomposed as w = k≥1 ⟨ w k , w⟩ w k . For a given eigenpair of R † R composed of the k-th optimal forcing ϕ k and gain η 2 k , the k-th optimal response w k is recovered by:

w k = η -1 k R ϕ k . (I.12)
Using both forcing and response spaces, equation (I.9) may be recast as:

w = k≥1 η k ⟨ ϕ k , ϕ⟩ w k (I.13)
As stated by Bennedine in [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF], when a strong convective instability is present in the flow, the optimal gain η 1 is such that η 1 ≫ η k≥2 , meaning the structure of ϕ is close to ϕ 1 , which yields a response w close to w 1 .

In a nutshell, the resolvent analysis improves the understanding of the genesis of nonmodal and nonlinear phenomenon within a linear framework. Such analysis has been used to describe the pseudo-resonance of Blasius boundary layer [START_REF] Brandt | Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer[END_REF], incompressible jets [START_REF] Garnaud | The preferred mode of incompressible jets: linear frequency response analysis[END_REF], flow around OAT15A airfoil [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], or to make separation control on a NACA0012 airfoil [START_REF] Yeh | Resolvent-analysis-based design of airfoil separation control[END_REF]. More recently, Lugrin [START_REF] Lugrin | Study of the viscous interaction with separation at hypersonic regime[END_REF] exploited the resolvent analysis to explain the amplification mechanism of oblique modes, which are linearly amplified convective instabilities appearing on hypersonic vehicles.

The tools of global stability may hence be useful to characterize the rotating stall phenomenon on axial compressors like the CME2. However, as depicted hereafter, there are few methods available for analyzing such large-scale configurations.

I.3 Methods for high dimensional eigenvalue problems

The global stability have been introduced within a continuous framework. In the context of Computational Fluid Dynamics (CFD), each operator and vectors presented earlier must be discretized. If the computational domain consists of N c cells, and the problem involves d conservative variables, the number of degrees-of-freedom (DOF) of the problem is N = N × d. The discretized state vectors are N -length vectors and the Jacobian operator is a N -by-N matrix. When N c is high, the modal and nonmodal global stability problems are therefore based on the resolution of large eigenvalue problems. As only the least stable eigenvalues of the spectrum are of interest, it would be counterproductive to use direct methods like QR factorization [START_REF] Francis | The QR Transformation A Unitary Analogue to the LR Transformation-Part 1[END_REF] that returns the entire spectrum of the considered operator. Iterative eigensolvers such as Arnoldi [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF] or Krylov-Schur (KS) [START_REF] Stewart | A Krylov-Schur Algorithm for Large Eigenproblems[END_REF] algorithms are preferred instead. Many methods related to the resolution of large eigenvalue problems have been described by Saad in [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF].

The counterpart of Arnoldi-like eigensolvers is that they are best at finding the eigenvalues of the largest magnitude of the considered operators. As discussed by Theofilis in [START_REF] Theofilis | Global Linear Instability[END_REF], the operators of modal analysis are preferentially modified with spectral transformations like shift-and-invert as to solve the eigenvalue problem. The latter, first used in the context of global stability analysis by Christodoulou & Scriven [START_REF] Christodoulou | Finding leading modes of a viscous free surface flow: An asymmetric generalized eigenproblem[END_REF], involves selecting a complex scalar (called shift) s, and solving the transformed problem:

(J -sI) -1 w = λ w, (I. [START_REF] Bénard | 2e semestre[END_REF] where λ is linked to the eigenvalue of the original problem by the relation λ = (λ -s) -1 . This technique allows for the calculation of only a portion of the spectrum of interest, as the largest magnitude eigenvalues of the transformed operator are those closest to s. Thus, whether dealing with modal or resolvent analysis, inversions of large complex-valued operators must be carried out. For 2-D problems and small 3-D ones, the usual strategy is to use a direct method like LU decomposition [START_REF] Bunch | Triangular factorization and inversion by fast matrix multiplication[END_REF]. This can be achieved using the efficient parallel MUMPS solver (see https://mumps-solver.org/) may be used. As stated by L'Excellent in [START_REF] Excellent | Multifrontal Methods: Parallelism, Memory Usage and Numerical Aspects[END_REF], if NE(J ) is the number of entries of the Jacobian, the memory cost of its LU factorization is of the order of (NE(J )/7) 4/3 . When dealing with complex configurations exceeding the 10 × 10 6 DOF, such as the full-annulus compressor CME2, storage of LU factorization thus becomes almost impossible due to memory limitations and inversions must be managed with iterative methods instead. As pointed out by Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF], due to the non-symmetrical, non-normal and not necessarily diagonal dominant nature of the operators encountered in compressible CFD, the most efficient current methods are preconditioned Krylov ones. Saad describes such methods in details in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. For example, Mack & Schmid [START_REF] Mack | A preconditioned Krylov technique for global hydrodynamic stability analysis of large-scale compressible flows[END_REF] used an ILU-based preconditioned Krylov technique, the BICGSTAB, to conduct global stability analysis about a swept parabolic body. Timme & Thormann discussed in [START_REF] Timme | Towards Three-Dimensional Global Stability Analysis of Transonic Shock Buffet[END_REF] the efficiency of various preconditioned Krylov methods for solving global stability problems. Timme also carried out the largest known direct and adjoint global stability analysis in [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF] to characterize the buffet onset of a full airplane configuration, resulting in approximately 50 × 10 6 mesh points for the finest grid studied. To the author's knowledge, no resolvent analysis has been performed on such large cases.

I.4 Positioning of the present work

Based on the previous sections, the understanding of the origin of rotating stall can be improved through the use of global stability analysis. However, due to the size of the problem, traditional direct factorization methods cannot be utilized, and the entire global stability process must be adapted.

At ONERA, work is underway to improve the high performance computing (HPC) layers of the compressible elsA CFD code [START_REF] Mary | Flexible aerodynamic solver technology in an HPC environment[END_REF][START_REF] Maugars | Algorithmic Differentiation for an effcient CFD solver[END_REF]. This includes cache-blocking and vectorization optimizations. Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] has also developed preconditioning tools to efficiently solve the large real-valued sparse linear systems involved in computing fixed points in compressible CFD using the GMRES algorithm [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. This is achieved by computing the Jacobian matrix with Algorithmic Differentiation (AD), which preserves the HPC layer, and using preconditioners based on ILU(k) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] and local-LU factorization. These developments could potentially be adapted for solving linear systems related to global stability problems.

The main obstacle to this is that the linear systems involved in global stability analysis are complex-valued, which is incompatible with vectorization as complex arithmetic is not vectorizable. To take advantage of both Guilbert's work and HPC optimizations, the complex operators need to be reformulated in a real-equivalent version, which should be propagated throughout the entire eigenvalue computation process.

I.5 Objectives and organization of the thesis

The objective of this work is to develop a method in elsA software to scale up the overall global stability process to industrial configurations with more than 10 × 10 6 DOF like the compressor CME2. To achieve this, the manuscript is divided into three chapters. The first chapter starts with a presentation of numerical methods used to discretize governing equations of compressible CFD. Then, efficient fixed point computations, a pre-requisite for global stability analysis, are discussed using Newton-like methods. The concept of AD for Jacobian computation is also covered, followed by an overview of various tools for global stability analysis in a discrete framework. Finally, the efficient iterative eigensolver KS algorithm [START_REF] Stewart | A Krylov-Schur Algorithm for Large Eigenproblems[END_REF], in its complex version, is presented. Each concept is illustrated using the viscous Burgers case [START_REF] Burgers | A mathematical model illustrating the theory of turbulence[END_REF], a 1-D nonlinear advection equation that resembles the Navier-Stokes equations.

In chapter II, after establishing the mathematical background for laminar and RANS computations, the fixed point methods as implemented in elsA are discussed. The focus is on the implicit backward Euler scheme solved with the approximate LU-SGS method and the exact GMRES-based [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] method. A good preconditioner, based on the exact Jacobian matrix, is required to ensure convergence of the latter method. The Jacobian is computed using multicoloring [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF], while the preconditioning is based on ILU(k) and hybrid local-LU methods. A KS algorithm compatible with HPC optimizations is also developed. As mentioned earlier, the operators of global stabil-ity analysis have to be reformulated in a real-equivalent version for effective vectorization. To do so, the K formulation, introduced by Day & Heroux in [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF] is selected, which conserves the sparsity pattern of the original complex operator and thus allows the reuse of the preconditioning routines developed by Guilbert [51]. This formulation is first implemented to be operative with the GMRES algorithm. The spectral properties of the K-formulated operators are then described, and detailed in the case of modal and resolvent analysis. Implementations of the KS algorithm applied to global stability in elsA are provided, giving rise to a new method named K-formulated global stability analysis method (KFGSM). The method is further validated on the canonical case of laminar cylinder flow for direct, adjoint, and resolvent computations.

In the last chapter, the KFGSM is tested on more challenging configurations. First, the RANS case of OAT15A airfoil flow is addressed, and the method is validated for direct, adjoint, and resolvent computations. The KFGSM is then tested on the subsonic laminar sphere flow, a 3-D configuration that reaches 10 × 10 6 DOF, representing industrial configurations for which the KFGSM was designed. The method is used to thoroughly characterize the two bifurcations of this case using direct, adjoint, and resolvent analysis. The code is even tested on a refined mesh with about 125 × 10 6 DOF for demonstration purposes.

II.1 Introduction

The present chapter builds a mathematical background to introduce the various tools of global stability analysis. The purpose is to describe each step of the process in a general way, in order to set these ideas before the K-formulated global stability method (KFGSM), as implemented in elsA [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF] software is detailed in chapter II. To this end, general considerations about conservation laws as well as numerical discretization, in the context of compressible Computational Fluid Dynamics (CFD) are first addressed. Secondly, some fixed point computation methods with the strategy developed to get linearized operators are detailed. Finally, the main tools of global stability analysis are discussed, and the Krylov-Schur (KS) [START_REF] Stewart | A Krylov-Schur Algorithm for Large Eigenproblems[END_REF] algorithm with complex arithmetics is introduced.

II.2 Governing equations in compressible CFD

The general conservation equations arising in fluid mechanics and their discretization within the framework of compressible CFD are presented here. Unlike the incompressible version of Navier-Stokes (NS) equations that can be elliptic or parabolic, their compressible counterpart is hyperbolic [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF]. In the latter description, the convective part must be addressed using Riemann solvers [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. This section presents the general concepts of spatial discretizations, particularly for compressible CFD solvers. These are illustrated on the one-dimensional Burgers equation [START_REF] Burgers | A mathematical model illustrating the theory of turbulence[END_REF], which mimics the NS equations. Mathematical formulations of the CFD problem treated with the elsA solver are presented in III.3.

II.2.1 Conservation laws

Let us consider the state vector w(M , t) ∈ R d containing d conservative variables, where M ∈ R N dim is the N dim -dimensional vector representing the spatial position and t ∈ R + is the time parameter. The transformation p is defined in the following such that the primitive variables are contained in the vector p(w) ∈ R d . For the sake of conciseness, the notation p will be used to denote p(w) in the following. From Newton's second law, the evolution of such vector is governed by the following system of partial differential equations (PDE):

∂w ∂t + ∇f = t s , (II.1)
where Considering the fluxes function to be of class C 1 , the time derivative can be decoupled from the space integral, and the divergence theorem can be applied to get:

f = f (p, ∇p) ∈ R d × R N dim is
∂ ∂t ˆΩ wdΩ + ˛∂Ω f • ndS -ˆΩ t s dΩ = 0, (II.3)
with n the unit external normal defined on the frontier ∂Ω of the domain. Physically, equation (II.3) links the temporal variation of w to its flux crossing the frontier ∂Ω and its production or destruction via the source term. Additionally, the flux term can be seen as the sum of two contributions: the convective part f c , and the diffusive part f d . The system can be recast as:

∂ ∂t ˆΩ wdΩ + ˛∂Ω (f c -f d ) • ndS -ˆΩ t s dΩ = 0. (II.4)
As a reminder, PDE can be classified as hyberbolic, parabolic or elliptic depending on their properties. Hirsch presents the physical meaning of each of these classifications in [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF]. A system of equations is said to be hyperbolic if its homogeneous part admits wave-like solutions. The physics of such system is dominated by convection phenomena. For instance, the so-called inviscid 1-D Burgers equation which describes the convection without diffusion of a scalar quantity u is purely hyperbolic:

∂u ∂t + u ∂u ∂x = 0. (II.5)
When the wave-like solution is damped in time, the equation is said to be parabolic. Such property is associated with time-dependent diffusion. A canonical case is the parabolic heat equation that governs the evolution of the temperature T in a fluid at rest of thermal coefficient κ:

∂T ∂t -κ∆T = 0. (II.6)
If the system does not admit any wave-like solution, it is said to be elliptic which is associated to pure spatial diffusion. Such behaviour describes the spatial diffusion of the pressure p in an incompressible 2-D flow:

∆p + ρ   ∂u x ∂x 2 + 2 ∂u x ∂y ∂u y ∂x + ∂u y ∂y 2   = 0. (II.7)
The nature of the conservation laws varies depending on the physical system of interest. In relation (II.4), the convective flux is an hyperbolic operator while the diffusive flux is parabolic. The system is purely hyperbolic when the latter is zero, as in the stationary Euler equations. But in the case of the NS equations, the presence of f d induces a hybrid mix hyperbolic-parabolic describing a convection-diffusion phenomenon. This behaviour can be reproduced by adding a viscosity parameter ν to the Burgers equation (II.5):

∂u ∂t + ∂ ∂x u 2 2 -ν ∂u ∂x = 0. (II.8)
In this case, f c = u 2 /2 and f d = -ν∂u/∂x. For the sake of better comprehension, the latter equation will be used to illustrate the various concepts presented in this chapter.

II.2.2 Spatial discretization and Riemann solvers

The following presents the spatial discretization by means of the finite volumes method. A solution to the problem of interest is searched on the domain Ω, divided into N c cells. The integral conservation law (II.3) applied in a cell Ω k at a time t writes:

∂ ∂t ˆΩk wdΩ + ˛∂Ω k f • ndS - ˆΩk t s dΩ = 0. (II.9)
If Ω k contains N int interfaces where the j-th one is denoted by I j , this relation becomes:

∂ ∂t ˆΩk wdΩ + N int j=1 ˛Ij f • ndS - ˆΩk t s dΩ = 0. (II.10)
We consider in addition the average value w k of w in cell Ω k :

w k = 1 V (Ω k ) ˆΩk wdΩ, (II.11)
where V (Ω k ) stands for the volume of Ω k . The general formulation of the finite volume method is thus given by:

∂(V (Ω k )w k ) ∂t + N int j=1 h c (w L I j , w R I j ) -h d (w L I j , w R I j ) S j -V (Ω k )T k = 0, (II.12)
where h c and h d are functions that compute convective part and diffusive part of the numerical approximations of the physical fluxes ¸Ij (F c -F d )•ndS crossing each interface I j of surface S j . w L I j and w R I j are the left and right (from I j ) state of conservative variables. In compressible CFD, the equations are dominated by convective fluxes, which makes them mainly hyperbolic. Special care must be taken to build h c to avoid oscillations across a discontinuity in the solution. The diffusive part h d can be seen as a numerical dissipation that damps discontinuities. Hereafter, some types of numerical schemes used in this thesis are described in a general way. Hyperbolic convective and parabolic viscous operators are separately considered to build each flux. These schemes are illustrated for Burgers equation. An extension to 3-D configurations, such as implemented in elsA solver, is presented in III.3.

Convective fluxes.

We focus in the following on the homogeneous hyperbolic part of the conservation law (II.1). The finite volume method uses the average value of w in each cell so the solution is a constant piecewise function, which creates discontinuities at each interface. Godunov [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] introduced a first-order upwind scheme to compute the intercell numerical flux h cI j , requiring the resolution of a nonlinear hyperbolic local Riemann problem [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. In 1-D (i.e. N dim = 1), this problem writes:

           ∂w ∂t + A c ∂w ∂x = 0, w(x, 0) = w L if x < 0 w R if x > 0 (II.13)
where A c = ∂F c /∂w ∈ R d × R d stands for the Jacobian matrix of the convective fluxes. The latter is diagonalizable with real eigenvalues, meaning there exist matrices

X , D ∈ R d × R d verifying A c = X DX -1
and where D is a diagonal matrix containing the eigenvalues of A c in its diagonal.

The expression for the Godunov flux based on the center of the interface I j then writes:

h cI j = 1 2 f c (p L I j ) + f c (p R I j ) • n I j - 1 2 |A c | w R I j -w L I j , (II.14)
where |A c | = X |D|X -1 , and |D| is the diagonal matrix containing the modulus of each eigenvalue. Such a discretization conserves the local information of wave propagation through the entire domain, and physical discontinuities like shock waves can be captured without generating additional spurious oscillations. The counterpart is the need to solve as many Riemann problems as there are interfaces in the domain. This strategy is time-consuming as it requires A c at each interface. Instead of looking for the exact solution of problem (II.13), another possibility is to search for approximate solutions by solving approximated Riemann problems. This idea is behind the Roe scheme [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] widely used in this work. The concept is to replace the Riemann problem with an

x1 2 x 3 2 x k-1 2 x k+ 1 2 x Nc+ 1 2 x Nc-1 2 Ω 1 Ω k Ω Nc cell interface . . . . . . Figure II.1: Segment discretization.
approximated linearized one, which is achieved by replacing A c in (II.13) by a so-called Roe matrix A c . The latter only depends on w L and w R and must satisfy the three following properties:

(i) A c (w L , w R ) is diagonalisable with real eigenvalues; (ii) A c (w L , w R )(w R -w L ) = f c (w R ) -f c (w L ); (iii) A c (w, w) = A(w).
(II.15)

By analogy with the Godunov flux, the Roe intercell flux is defined as:

h cI j = 1 2 f c (p R I j ) + f c (p R I j ) • n I j - 1 2 A c w R I j -w L I j , (II.16)
Let us pinpoint that a lot of schemes based on the formulation of various Riemann problems at the interfaces can be used, like the so-called HLLC scheme, and more details about their derivation can be found in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. Such a scheme ensures a first-order-accurate discretization in space. They can be extended to second-order accuracy using the Monotonic Upstream-centered Scheme for Conservation (MUSCL) reconstruction method introduced by Van Leer in [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection[END_REF]. The latter assumes that the left and right states are no longer constant but may vary linearly through the cells. Second-order accuracy is obtained using a reconstruction at the faces of left and right states. If the left and right cells of interface I j respectively are Ω k and Ω l , then one may write:

   p L I j = p k + ∇p| k • d k I j p R I j = p l + ∇p| l • d l I j , (II.17)
where d k I j is the vector going from the center of cell Ω k to the center of interface I j . The use of MUSCL reconstruction leads to spurious oscillations close to flow discontinuities. This can be countered by using limitor functions for the reconstructed values [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF].

In the case of Burgers equation, the domain Ω is a segment [x 1 , x Nc ] in a cartesian coordinates system with a uniform space grid ∆x, and each cell Ω k has two interfaces in x k-

h c k+ 1 2 = 1 2 f c (u L k+ 1 2 ) + f c (u R k+ 1 2 ) - 1 4 |u k + u k-1 | (u R k+ 1 2 -u L k+ 1 2 ), (II.18) with u L k+ 1 2 and u R k+ 1 2
defined as:

         u L k+ 1 2 = u k - ∆x 2 ∇u| k u R k+ 1 2 = u k+1 + ∆x 2 ∇u| k . (II.19)
Diffusive fluxes. Diffusive fluxes are calculated by a discretization purely centered at interface I j , based on the primitive variables:

h dI j = f d p I j , ∇p| I j • n I j . (II.20)
The primitive variables at the center of I j are estimated by arithmetic mean of their values at the center of each adjacent cells. Gradients are also computed at interfaces. Various methods can be applied, but the one used in elsA is a corrected version of the so-called five-points stencil (5p) method. It is described in III.3.2.

In the case of Burgers equation, the numerical diffusive flux writes:

h d k+ 1 2 = ν ∇u| k+ 1 2 • n k+ 1 2 . (II.21)
The gradient term is obtained by a Taylor expansion of u k+ 1

2

. At second-order, it gives:

u k+ 1 2 = u k + ∆x 2 ∇u| k+ 1 2 . (II.22)
From what precedes, the quantity u at the interface in

x k+ 1 2 is u k+ 1 2 = (u k + u k+1 )/2
, from which one deduces the expression of the gradient at the face:

∇u| k+ 1 2 = 1 2 (u k+1 -u k ) . (II.23)
Hence, the numerical diffusive flux reads:

h d k+ 1 2 = ν ∆x (u k+1 -u k ) . (II.24)
In a similar way, we get at fourth order:

h d k+ 1 2 = ν 24∆x (-u k-1 + 26u k -u k+1 ) . (II.25)
Now that the finite volume method and flux discretization have been introduced in the context of compressible CFD, the next section presents methods to compute fixed points. In the following, the discrete representation of all vectors and scalars will be written in upper case. For example, the discretized state vector W ∈ R N , with N = d × N c , is composed of the N c state vectors w. We also introduce the discrete residual operator R(W ) ∈ R N whose k-th component reads:

t sk - 1 V (Ω k ) N int j=1 h(w L I j , w R I j )S j . (II.26)
The general discretized governing equations can thus be recast under the following form:

∂W ∂t -R(W ) = 0. (II.27)

II.3 Fixed point computation

As introduced in I.2, a global stability analysis studies the linear evolution of global harmonic perturbations about a baseflow W b to determine its asymptotic stability. Therefore, being able to compute W b accurately is a prerequisite. In this work, W b is chosen to be a fixed point, i.e. a steady solution of equation (II.1). Note that it is also physically relevant to carry out a stability analysis on the mean flow instead of a fixed point (see for instance [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]). But this latter approach is out of the scope of the present thesis and will not be discussed.

Finding a steady-state solution of equation (II.27) results in the following equivalence:

∂W ∂t = 0 ⇔ R(W ) = 0. (II.28)
Such a problem can be solved by Newton-like methods or a pseudo time-marching loop. However, from time integration, the use of explicit schemes such as the Runge-Kutta scheme [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF] could be very expansive in terms of CPU usage. This is why only implicit Newton-like methods were used in this work.

Hereafter, the Newton method and the linearized first-order backward Euler method with local time-step are derived for the conservative equations. They are then illustrated on the resolution of the viscous Burgers problem defined by equation (II.8). We focus on its resolution on the segment [-10, 10], with u(-10, t) = 1 and u(10, t) = -1, and where N c = 400. A theoretical fixed point of such a boundary-valued problem is known to be u theo = -tanh(x/2) [START_REF] Guba | Continuous spectra and numerical eigenvalues[END_REF]. Both convective and diffusive fluxes are discretized at second-order accuracy using respectively (II.18) and (II.24).

II.3.1 Newton method

Given an initial guess W 0 , the Newton method iteratively evaluates the gradient direction toward a fixed point. The n-th Newton iteration is based on the following linearization:

R(W n+1 ) = R(W n + ∆W ) ≃ R(W n ) + J (W n )∆W , (II.29)
where J (W n ) = ∂R/∂W | Wn is the Jacobian matrix of the residual. Each iteration consists in evaluating ∆W through the resolution of the following linear system:

J ∆W = -R(W n ), (II.30)
and the solution of the nonlinear problem at the next iteration is given by W n+1 = W n + ∆W . The algorithm stops when R(W n+1 ) becomes lower than some user-defined tolerance. Under certain assumptions, the convergence is expected to be quadratic. This method has been developed on the Burgers case. Given the differentiated right Jacobian-vector product obtained using the Algorithmic Differentiation (AD, presented in II.4), the Jacobian J is computed by applying this product to each canonical vector e i = 0

• • • 0 1 0 • • • 0 T . The system (II.

II.3.2 Implicit backward Euler method

The Newton method might not be robust enough when the physical problem is particularly stiff2 , or when the initial guess is too far from the searched solution. In fluid mechanics, a common practice to overcome this difficulty is solving the nonlinear system (II.27) with an implicit scheme in time. In particular, we focus in this thesis on the linearized first-order backward-Euler method with a local time-step that leads to the following equation:

1 ∆t I N -J (W n ) ∆W = R(W n ), (II.31)
where I N is the N -by-N identity matrix. A linear system has to be solved here as well, and one may see that when ∆t → +∞, this scheme is equivalent to a step of the Newton method (II.30).

To compute a steady-state solution, Crivellini et al. [START_REF] Crivellini | An implicit matrix-free Discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations[END_REF] developed a pseudo-transient continuation technique that uses an adaptative Courant-Friedrichs-Lewy (CFL) number [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] defined as:

           CFL = min CFL min r β , CFL max , r = max ∥R(W n )∥ 2 ∥R(W 0 )∥ 2 , ∥R(W n )∥ ∞ ∥R(W 0 )∥ ∞ , (II.32)
where CFL min , CFL max and β are user-defined input parameters. In this work, β is always 1.

During computation, the CFL number grows while the residual norms decrease. Since the timestep ∆t is proportional to the CFL number, its value grows as the algorithm becomes closer to a steady-state solution. This approach is physically relevant since the goal is to capture asymptotical phenomena that do not vary even during large periods of time. Quadratic convergence is obtained during this late part of the computation. This fixed point method has been implemented on Burgers equation as well. Fig. II.3a shows the convergence curve resulting from a computation initialised with U 0 = 1. During the first 450 iterations, the system is more relaxed than that in Newthod method and converges slower with a lot of oscillations. However, the solution reaches then a quadratical convergence toward the basin of attraction. The resulting field provided in fig. II.3b also presents a good match with U theo , proving the robustness of the method.

Newton-like methods perform well at finding steady solutions with good accuracy. However, to ensure a convergence below 10 -10 , the Jacobian operator must be computed with high reliability. As mentioned earlier, the illustrative results from this section were obtained with AD to evaluate Jacobian-vectors products. The AD theory and its practical implementation are detailed in the following section.

II.4 Linearization procedure

Whether we work on fixed point computation or global stability analysis (see II.5), critical is the importance of the linearization procedure for direct and adjoint operators. One solution is differentiate equation (II.1) analytically, then to discretize the set of linearized equations [START_REF] Peter | Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches[END_REF][START_REF] Jameson | Aerodynamic design via control theory[END_REF]. However, as stated in [START_REF] Peter | Contributions to discrete adjoint method in aerodynamics for shape optimization and goal-oriented mesh-adaptation[END_REF], in the compressible case, the consistency from continuous to discrete operator is not guaranteed. The solution used in this thesis is to work in a fully discretized framework in which equation (II.1) is discretized into equation (II.27), and then linearized. Mettot [START_REF] Mettot | Linear stability, sensitivity, and passive control of turbulent flows using finite differences[END_REF] states such strategy has the other advantage of being less error-prone and simpler to implement since all the complexity, such as boundary conditions or numerical schemes, is already accounted for in the discrete linearized equations without additional treatment.

As introduced in [START_REF] Beneddine | Characterization of unsteady flow behavior by linear stability analysis[END_REF], one way to approximate the product between J and a test vector φ is to use the following finite differences relation:

J φ ≈ 1 ϵ (R(W b + φ) -R(W b )) , (II.33)
where ϵ is a small user-defined parameter. Nevertheless, Knoll & Keyes in [START_REF] Knoll | Jacobian-free Newton-Krylov Methods: A Survey of Approaches and Applications[END_REF] explained that the choice of ϵ can be tedious and case-dependent. If too large, the matrix-vector product is poorly approximated, and if too small, floating-point round-off errors could contaminate the result. In this work, the issue has been avoided by using AD tools instead, which provides exact matrix-vector products.

II.4.1 Description of Algorithmic Differentiation with Tapenade

There are two main families of AD tools: the Source Code Transformation (SCT) which replaces the source code for a function by an automatically generated source code that includes statements for calculating the derivatives interleaved with the original instructions; and the Overloading Operator (OO) that computes the derivatives present in a computer code by appending a dual component to all the dependent variables. For performance issues [START_REF] Maugars | Algorithmic Differentiation for an effcient CFD solver[END_REF], the SCT has been retained in this work. In particular, the choice was made to use the Tapenade software, developed by Hascoët & Pascual [START_REF] Hascoet | The Tapenade Automatic Differentiation tool: principles, model, and specification[END_REF] as it is the most mature and reputable open-source tool in the community nowadays.

Given m, l ∈ N and an algorithm that implements a function G : X ∈ R m → G(X) ∈ R l , Tapenade produces differentiated codes in SCT that computes the exact Jacobian-vector product with a column (forward mode) vector by the right, or a line vector (backward mode) by the left, up to machine precision. Assuming G is differentiable, AD relies on the fact that its implementation is composed of a sequence of p elementary instructions, which are mathematical functions g that can be differentiated one by one. The composition of these functions returns G:

G = g p • g p-1 • • • • • g 1 .
(II.34)

From such a relation, the chain rule applies on the derivatives to get the Jacobian operator J G = ∂G/∂X:

J G (X) = g ′ p (X p-1 ) × g ′ p-1 (X p-2 ) × • • • × g ′ 1 (X 0 ) (II.35) with X 0 = X, X k = g k (X k-1 ) and g ′ k = ∂g k /∂X k-1 .
Forward mode. From a given column vector X F ∈ R m , the quantity J G (X)X F is evaluated with successive right matrix-vector products:

J G (X)X F = g ′ p (X p-1 ) × g ′ p-1 (X p-2 ) × • • • × =X F 1 g ′ 1 (X 0 )X F =X F p-1 , (II.36)
which, as shown in [START_REF] Hascoet | The Tapenade Automatic Differentiation tool: principles, model, and specification[END_REF], is implemented using the following standardized sequence of instructions:

X 0 = X X F 0 = X F X 1 = g 1 (X 0 ) X F 1 = g ′ 1 (X 0 )X F 0 . . . X k = g k (X k-1 ) X F k = g ′ k (X k-1 )X F k-1 . . . X p = g p (X p-1 ) X F p = g ′ p (X p-1 )X F p-1 , (II.37) and J G (X)X F = X F p .
Backward mode. This mode computes, for a given X B ∈ R l , the quantity X B T J G (X) with successive left matrix-vector products:

X B T J G (X) = =X B p-1 X B g ′ p (X p-1 ) ×g ′ p-1 (X p-2 ) × . . . =X B 1 ×g ′ 1 (X 0 ). (II.38)
Note that actually, since X B T J G (X) = (J T G (X)X B ) T , the backward mode provides the action of J T on a given vector, which is particularly useful in the framework of global stability analysis (see II.5.2 and [START_REF] Luchini | Adjoint Equations in Stability Analysis[END_REF]). Unlike the forward mode, the vectors X k must all be computed before calculating the first backward quantity X B g ′ p (X p-1 ):

X 0 = X X 1 = g 1 (X 0 ) . . . X k = g k (X k-1 )
. . .

X p = g p (X p-1 ) X B p = X B g ′ p (X B ) X B p-1 = X B p g ′ p-1 (X p-2 ) . . . X B k = X B k+1 g ′ k (X k-1 ) . . . X B 0 = X B 1 g ′ 1 (X 0 ), (II.39)
and

X B T J G (X) = X B 0 T .
Appendix A presents an example of AD code in forward and backward modes obtained with Tapenade on a piece of FORTRAN code that computes the L 2 scalar product. One may check the accuracy of the produced routines by comparing their outputs to finite differences using equation (II.33) or to the complex-step derivative approximation introduced by Martins [START_REF] Martins | The complex-step derivative approximation[END_REF]. To validate the compatibility between forward and backward differentiated codes, a so-called duality test may be performed [START_REF] Paulin | Automatic verification of algorithmically differentiated code[END_REF]. Given two random vectors X F ∈ R m and X B ∈ R l , it consists in checking the accuracy of the following equality:

X B • J G (X)X F forward = X B J G (X) backward •X F , (II.40)
which is in practice achieved by evaluating the relative error ϵ AD rel :

ϵ AD rel = X B • J G X F -X B J G • X F ∥J G X F ∥ 2 .
(II.41)

II.4.2 Application to the conservation equations

Significant efforts have been made at ONERA, based on the work of Mary [START_REF] Mary | Flexible aerodynamic solver technology in an HPC environment[END_REF], to challenge the two limiting factors of CFD applications: their CPU and memory costs. In elsA CFD solver, this so-called high performance computing (HPC) layer consists of cache-blocking and vectorization optimizations. The computational domain is divided into packs of tiles (cells or interfaces, depending on what is computed) independent of each other for vectorization. Once a pack is loaded in memory, the goal is to perform as many operations as possible before other packs have to be treated. Such a technique optimizes loadings and unloadings of memory cache at the L3 level (see fig. III.1), and allows for the vectorization of calculus loops since each tile is independent of others. In practice, the solver is organised as a chain of operators specified by user choices, which ensures these improvements on direct, forward, and backward computations.

The resulting chain of operators can be represented by a graph. For example, the graph of R : W ∈ R N → R(W ) ∈ R N with second-order schemes for fluxes, without source terms (T s = 0) is:

W F R ∇P P
where P and F are respectively the discretized primitive vector and fluxes defined in II.2.1. Such graph illustrates how each operator must be successively computed and nested one by one to finally give R(W ). This formulation allows the use of chain rules, and the Jacobian operator J = ∂R/∂W can be computed with crossed derivatives:

J = ∂R ∂F ∂F ∂P + ∂F ∂∇P ∂∇P ∂P ∂P ∂W . (II.42)
Tapenade is able to generate differentiated code for each operators in forward or backward mode. The graph providing the right Jacobian-vector product, i.e. the mapping φ ∈ R N → J φ ∈ R N , reads:

φ ∂P ∂W φ ∂F ∂P ∂P ∂W φ + ∂F ∂∇P ∂∇P ∂P ∂P ∂W φ J φ ∂∇P ∂P ∂P ∂W φ
while that of the left Jacobian-vector product, i.e. the mapping ψ ∈ R N ψ T → J ∈ R N , reads

ψ ψ T ∂R ∂F ψ T ∂R ∂F ∂F ∂P + ψ T ∂R ∂F ∂F ∂∇P ∂∇P ∂P ψ T J ψ T ∂R ∂F ∂F ∂∇P
where the terms between parenthesis are line or column vectors following the mode used.

Once coded and assembled, the chain of operators can be integrally differentiated in forward or backward mode to get the action of both J and J T (with the relation J T ψ = (ψ T J ) T ) on a given vector at machine precision. The use of SCT allows the preservation of all HPC optimizations as in the direct solver, unlike the OO differentiation [START_REF] Maugars | Algorithmic Differentiation for an effcient CFD solver[END_REF]. As depicted in II.4.1, the validity of each differentiated operator is verified by comparing with finite differences, and the compatibility between forward and backward operators is checked using duality test (see equation (II.40)). Such a strategy has been implemented for the Burgers equation and the resulting duality error associated with J (U theo ) is ϵ AD rel = 1.6467 × 10 -16 . This is below machine precision, validating the differentiated codes produced by Tapenade.

II.5 Global stability analysis

The stability of a base flow can be studied by analyzing its behavior under the effect of a perturbation. As discussed in I.2, one or more global modes are necessarily involved in the asymptotical behaviour of an unstable base flow. Their spatial support, amplification rate and frequency can be computed by finding the least stable eigenpairs of J [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. When the baseflow is stable, the global modes are exponentially and asymptotically damped, and are not sufficient to give explanation of transient energy growth in the system. The resolvent analysis [START_REF] Threfethen | Hydrodynamic stability without eigenvalues[END_REF][START_REF] Schmid | Stability and Transition in Shear Flows[END_REF] examines such short-term instability providing the most amplified frequency (and the associated forcing and response) of a given sub-critical flow. The direct and adjoint global stability as well as the resolvent analysis in their discrete form are detailed hereafter.

II.5.1 Direct modal analysis

By analogy with I.2, the discrete formulation of the direct eigenvalue problem is detailed hereafter. We assume first the existence of a steady solution W b of (II.27) perturbed by a small disturbance W ′ . The discrete residual can be linearized about W b using:

R(W b + W ′ ) ≃ R(W b ) =0 +J (W b )W ′ .
(II.43)

For the sake of conciseness, the notation J will be used to denote J (W b ) in the following. Injecting this last quantity in equation (II.27) results in the so-called linearized Navier-Stokes (LNS) equations that govern the temporal evolution of W ′ :

∂W ′ ∂t = J W ′ . (II.44)
As stated for instance by Sipp et al. in [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF], these disturbances may be decomposed into modes of the form W ′ (x, t) = W (x)e λt , where W is an eigenfunction of J representing the spatial structure of this direct mode. The quantities σ = ℜ(λ) and ω = ℑ(λ) respectively stand for the growth rate and angular frequency associated with W ′ . Equation (II.44) then reduces to the following eigenproblem:

λ W = J W . (II.45)
W b is said to be unstable (respectively stable) if there exists at least one eigenvalue of J such that σ > 0 (respectively σ < 0). The global stability of W b is thus analyzed by computing the spectrum of J . In this work for the sake of generality, we assume that the flow does not exhibit any spatial homogeneity, but when there is one or more, the study simplifies into bi-global or local stability analysis [START_REF] Juniper | Modal Stability Theory: Lecture notes from the FLOW-NORDITA Summer School on Advanced Instability Methods for Complex Flows[END_REF].

Usually, only the least stable modes of J are of interest, and the eigenproblem (II.45) is solved with an iterative eigensolver3 like the Krylov-Schur algorithm (KS, presented later in section II.6). This algorithm is best at finding the eigenvalues of the largest magnitudes of the considered operator, which are not necessarily the least stable ones. To target the least stable part of the spectrum, the spectral shift-and-invert transformation is applied: a complex-valued shift s is introduced, and the following transformed complex eigenproblem is solved instead:

λ W = (J -sI N ) -1 W , (II.46)
where λ = (λ -s) -1 . This strategy is further explained in II.6.5, which highlights its efficiency to filter the eigenvalues closest to s. The stability of the Burgers equation's fixed point W b computed in II.3.2, has been tested. The computation of eigenpairs was carried out by both the direct QR algorithm (given the relatively small dimensions of J ) with LAPACK [START_REF] Anderson | LAPACK: A Portable Linear Algebra Library for High-Performance Computers[END_REF]; and by SLEPc [START_REF] Vicente | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF] with a KS algorithm with shiftand-invert transformation, where only the least stable eigenvalues are looked for. The inversions of the complex operator are managed in SLEPc by LU factorization, as in section II.3. The resulting spectra are showed in fig. II.4a, where a good match between LAPACK and SLEPc for the least stable eigenvalues. We also notice that W b is a stable base flow. The least stable eigenvector computed with SLEPc is presented in fig. II.4b. The largest part of the global mode is located is located at the center of the domain.

II.5.2 Adjoint modal analysis

To describe the adjoint modal analysis in its discrete form, we first introduce a discrete scalar product ⟨

•, •⟩ Q : C N ×N → C, represented by the mass matrix Q ∈ R N × R N . For two given vectors A, B ∈ C N , it reads: ⟨A, B⟩ Q = A H QB, (II.47)
where the superscript H denotes the conjugate transpose. Assuming that Q is a non-singular symmetric matrix, the discrete adjoint operator J † of J is defined relatively to ⟨•, •⟩ Q such that ⟨A, J B⟩ Q = ⟨J † A, B⟩ Q for all vectors A and B. Since J is real and Q is non-singular, an explicit relation between J † and J T is obtained:

J † = Q -1 J T Q. (II.48)
The matrices J † and J T are thus similar and share the same spectrum. We can therefore focus on computing the eigenvalues of J T formed with the differentiated code provided by Tapenade in backward mode, and pre-multiply the resulting modal basis with Q -1 to recover adjoint eigenfunctions. For a given direct eigenpair (λ k , W k ), a corresponding adjoint eigenpair (λ k , W † k ) solution of the following adjoint eigenproblem exists:

λ k W † k = J † W † k , (II.49)
where W † k is the adjoint eigenfunction. Solving (II.49) is achieved similarly to problem (II.45): the eigenpairs of J T are computed using the shift-and-invert transformation, and the resulting eigenvectors are pre-multiplied by Q -1 to get the adjoint global modes.

In fluid mechanics, J is often non-normal 4 and due to the nonorthogonality of eigenmodes, the expansion of W ′ over the direct modal basis ( W k , k ≥ 1) can't be achieved easily [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. Such a difficulty can be overcome by introducing the dual basis made of the adjoint eigenmodes (

W † k , k ≥ 1), which can be normalized such that ⟨ W † k , W l ⟩ Q = δ kl .
The direct and adjoint modal bases form a bi-orthogonal basis over which W ′ can be uniquely expressed:

W ′ = k≥1 ⟨ W † k , W ′ ⟩ Q W k . (II.50)
The adjoint mode provides the receptivity of the direct mode [START_REF] Schmid | Analysis of Fluid Systems: Stability, Receptivity, Sensitivity: Lecture notes from the FLOW-NORDITA Summer School on Advanced Instability Methods for Complex Flows[END_REF] to a given forcing or variation of an initial condition for instance. Additionally, the adjoint mode is the exact initial perturbation that triggers the direct mode in the flow. In finite volumes, the scalar product is usually chosen to correspond to the L 2 -norm. On a domain Ω, the latter is defined as ⟨A,

B⟩ Q Ω = ´Ω A H B dΩ = A H Q Ω B,
where Q Ω is a real diagonal matrix whose coefficients are the volumes V (Ω k ) of each cell. A numerical experiment on the Burgers equation in the same conditions as for the direct eigenproblem has been carried out, and the resulting spectra are showed in fig. II.5a. Again, a good match between LAPACK and SLEPc for the least stable eigenvalues is noticed. The least stable adjoint eigenfunction computed with SLEPc is presented in fig. II.5b. The largest parts of the adjoint global mode is mainly located at the limits of the domain.

II.5.3 Resolvent analysis

Based on the theoretical framework introduced in I.2, the discrete version of the resolvent analysis, introduced by Trefethen & al. [START_REF] Threfethen | Hydrodynamic stability without eigenvalues[END_REF], is depicted hereafter. The goal is to analyse the transient behaviour of a stable base flow under the effect of an external forcing Φ ′ . Adding such a term in equation (II.44) leads to:

∂W ′ ∂t = J W ′ + Φ ′ . (II.51)
The Fourier transform of equation (II.51) gives:

iω W = J W + Φ. (II.52)
Introducing the harmonic response W and forcing Φ at the frequency ω ∈ R, such that W ′ = W (x)e iωt and Φ ′ = Φ(x)e iωt , equation (II.51) can be recast as:

W = R(ω) Φ, (II.53)
with R(ω) = (iωI -J ) -1 the resolvent operator, a transfert function between W and Φ that is a selective noise-amplifier. For the sake of conciseness, the notation R will be used to denote R(ω) in the following. The forcing responsible for the maximum energetic response at the frequency ω in the flow is computed by maximizing the energetic gain η(ω

) = ∥ W ∥ Q 1 /∥ Φ∥ Q 2 associated to R,
where the norms are respectively defined relatively to the scalar products Q 1 and Q 2 . As for the adjoint discrete problem, Q 1 and Q 2 are non-singular symmetric real matrices. In practice, Q 1 = Q Ω and Q 2 is associated to kinetic or acoustic energy. In the computations carried out with elsA in this work, Q 2 is associated with the Chu energy [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF] (see equation (III.110)). Considering iω not to be an eigenvalue of R, meaning R and R H are defined5 , the quotient η(ω) can be expanded as:

η = ⟨ W , W ⟩ Q 1 ⟨ Φ, Φ⟩ Q 2 = ⟨R Φ, R Φ⟩ Q 1 ⟨ Φ, Φ⟩ Q 2 = Φ H R H Q 1 R Φ Φ H Q 2 Φ , (II.54) with R H = (-iωI N -J T ) -1 . Since R H Q 1 R is Hermitian (i.e. R H Q 1 R = (R H Q 1 R) H
), such a relation corresponds to a Rayleigh quotient and the problem can be recast in the following way:

Q -1 2 R H Q 1 R Φ = η 2 Φ. (II.55)
The solutions of this eigenproblem are the eigenpairs of the Hermitian operator

Q -1 2 R H Q 1 R.
Finding the forcing that maximizes η is thus equivalent to calculate the eigenvector associated to the largest eigenvalue of

Q -1 2 R H Q 1 R.
The gain η can alternatively be expanded with the harmonic response, leading to the following eigenproblem:

RQ -1 2 R H Q 1 W = η 2 W . (II.56) Both operators RQ -1 2 R H Q 1 and Q -1 2 R H Q 1 R
share the same real spectrum, and their eigenvectors respectively form the basis of optimal forcing ( Φ k , k ≥ 1) and response ( W k , k ≥ 1). The response W k associated with the k-th gain η k is linked with the forcing Φ k by the relation:

η k W k = R ϕ k , (II.57)
and the harmonic response can be decomposed in the following way:

W = k≥1 η k ⟨ Φ k , Φ⟩ W k . (II.58)
The vectors composing the forcing and response space respectively are right and left singular vectors of R, and the gains are its singular values. The curve of η(ω) is established by solving one of the eigenproblems (II.55) or (II.56). In this work, we focus on the problem (II.56) as to both benefit from the optimized Jacobian-vector products and avoid inversions when recovering the responses from the forcings. For a given frequency ω, Φ 1 is the optimal forcing, i.e. the forcing providing the most energetic response in the flow. In general, η 1 ≫ η k≥2 and the gains η k≥2 are the sub-optimal ones. The KS algorithm is suitable to solve the eigenproblem as it provides the eigenvalues of largest magnitudes of RQ -1 2 R H Q 1 . One may notice that this operator is composed of the inverse matrices R = (iωI N -J ) -1 and R H = (-iωI N -J T ) -1 , meaning its application to a given vector requires two inversions of complex linear systems. The computation of the KS basis (see II.6.2) is hence much more expansive than for the modal analysis. However, as the optimal gain is much larger than sub-optimal ones, the KS process usually requires less iteration to accurately converge eigenpairs, as shown for instance in chapter IV. A resolvent analysis has been carried out on the Burgers equation with This section was devoted to the derivation of the global stability tools in their discrete forms. As mentioned several times, both modal and non modal analyses require the resolution of complexvalued eigenproblems, usually addressed by Arnoldi-like methods like the KS algorithm. The details of the KS algorithm are exposed hereafter.

Q 1 = Q 2 = Q Ω .

II.6 Krylov-Schur algorithm

The KS algorithm, introduced by Stewart in [START_REF] Stewart | A Krylov-Schur Algorithm for Large Eigenproblems[END_REF], is an Arnoldi-based method that projects the eigenvectors of a given complex operator onto a Krylov space. This iterative process provides approximate eigenpairs, called Ritz pairs, such that the eigenvalues are of largest magnitude. In this section, some concepts about the Krylov decomposition through Arnoldi process, as well as the mechanics of the KS algorithm are exposed.

II.6.1 Krylov subspace and Arnoldi algorithm

The Krylov subspace K m (C, v) of order m, of an operator C ∈ C N × C N associated to a given vector v ∈ C N , is defined as follows:

K m (C, v) = span v, Cv, C 2 v, . . . , C m-1 v . (II.59)
The principle of Arnoldi-based methods, as shown in [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF], is to project C onto K m and extract eigenpairs of this projection, which are approximations of the effective eigenpairs of C. To do this, an orthonormal basis of K m is produced using the Arnoldi process (presented in Algorithm 1), based on the modified Gram-Schmidt algorithm [START_REF] Rice | Experiments on Gram-Schmidt Orthogonalization[END_REF], and the matrix C is reduced to an upper Hessenberg m-by-m matrix H m . The output of Arnoldi process is a particular case of Krylov Algorithm 1 Arnoldi process

Input: C, v 1 unitary (i.e. ∥v∥ 2 = 1) and m ∈ N Output: V m , H m , v m+1 and h m+1,m 1: for k ∈ 1, m do 2: v k+1 ← Cv k 3: for l ∈ 1, k do 4: h kl ← ⟨Cv k , v l ⟩ 2 5: v k+1 ← v k+1 -v l h kl 6:
end for 7:

h k+1,k ← ∥v k+1 ∥ 2 8: if h k+1,k = 0 then 9:
v k+1 ← v k+1 /h k+1,k 12:
end if 13: end for decomposition:

CV m = V m H m + v m+1 ĥH m+1,m , (II.60)
where V m is an N -by-m matrix whose columns are an orthonormal basis of K m , and ĥH

m+1,m = h m+1,m e H m .

II.6.2 General Krylov decomposition

Let us now consider an order-m general Krylov decomposition that reduces the matrix C to a matrix B m with an orthonormal basis U m :

CU m = U m B m + u m+1 bH m . (II.61)
Using the orthonormality of the columns of U m and pre-multiplying with U H m provides:

U H m CU m = B m , (II.62)
where B m is an approximate space whose spectrum is composed of the Ritz values λ k ∈ C. The Ritz vectors define the projection of the eigenspace of C onto B m . The k-th Ritz vector is defined as

U m y k , with y k ∈ C m the k-th eigenvector of B m .
The decomposition (II.61) can be post-multiplied by y k to get:

CU m y k = U m B m y k + u m+1 bH m y k . (II.63)
After recasting, taking the L 2 -norm and due to the unitarity of u m+1 , the exact error commited on the k-th Ritz vector is:

∥CU m y k -λU m y k ∥ 2 = | bH m y k |. (II.64)
The term | bH ŷk | is relatively cheap to compute as m is several orders of magnitude smaller than N . When | bH ŷk | is sufficiently small, the Ritz pair (λ k , U m y k ) is considered to be an approximate eigenpair of C:

CU m y k ≃ λU m y k , (II.65)
It can additionally be shown that an orthogonal transformation of B m produces a new Krylov decomposition equivalent to (II.61). We introduce a m-by-m unitary matrix W6 that defines the unitarily equivalent transformation of

B m = W H T m W, with T m ∈ C m × C m .
Both B m and T m have the same Ritz values and equation (II.62) can be recast as:

W H U H m CU m W = T m , (II.66)
which means that T m is also an approximate space. It can further generate the following new Krylov decomposition:

CU m W = U m WT m + u m+1 bH m W, (II.67)
with U m W the new Krylov space. This illustrates the possibility of transforming a Krylov decomposition without modifying the spectral properties of the initial approximate space. In the KS algorithm, this helpful property is used to sort the targeted portion of the spectrum of C following a provided criterion as described hereafter.

II.6.3 Krylov-Schur process

The KS process is based on the iterative expansion and contraction of a Krylov decomposition, by keeping the spectral informations of interest and enriching them at each iteration. More details about this method, as well as its variations (like spectral KS) and the various possible implementations, can be found in [START_REF] Stewart | A Krylov-Schur Algorithm for Large Eigenproblems[END_REF][START_REF] Kressner | Numerical Methods for General and Structured Eigenvalue Problems[END_REF][START_REF] Hernández | Krylov-Schur Methods in SLEPc[END_REF].

Let us first produce an order-m Krylov decomposition of operator C using the Arnoldi process:

CV m = V m H m + v m+1 ĥH m+1,m . (II.68)
With the QR algorithm [START_REF] Francis | The QR Transformation A Unitary Analogue to the LR Transformation-Part 1[END_REF], a Schur matrix S m containing the Ritz values on its diagonal can be calculated from the approximate space H m , resulting in the relation

H m = W 1 S m W H 1 ,
with W 1 a unitary matrix. A sorting represented by the unitary transformation W 2 can be applied without modifying the spectral content of H m . In practice, the KS algorithm performs best at computing the eigenvalues of largest magnitude. To do so, Q 2 is chosen so that the Ritz values are arranged in descending magnitude from the top left corner to the bottom right one of S m , leading to the sorted Schur matrix T m . From the results of II.6.2, by introducing the KS basis

V m = V m W 1 W 2 and bH m = ĥH m+1,m W 1 W 2 , the new decomposition C V m = V m T m + v m+1 bH m (II.69)
is a Krylov decomposition equivalent to (II.68). This relation can also be writen as:

C V m = V m T p ⋆ 0 T m-p + v m+1 bH p bH m-p , (II.70)
where T p contains the spectral information relative to the p wanted Ritz vectors. The upper triangular shape of T m allows the decomposition to be truncated at order p:

C V p = V p T p + v p+1 bH p , (II.71)
with v p+1 = v m+1 . This step of truncation purges the unwanted part of the spectrum before a re-extension at order m through m -p steps of the Arnoldi process is performed at next iteration.

From the second iteration, the approximate space is not a Hessenberg matrix anymore. For instance, if m = 6 and p = 3, the unsorted new approximate space has the following shape:

         

s s s h h h 0 s s h h h 0 0 s h h h b b b h h h

0 0 0 h h h 0 0 0 0 h h           . (II.72)

II.6.4 Deflation of Ritz values

This process of successive extension and contraction is repeated until the desired number of Ritz values are converged. Given a Ritz value λ k and based on equation (II.73), the error associated to the Ritz vector V m y k is easily calculated:

∥C V m y k -λ k V m y k ∥ 2 = | bH m y k |. (II.73)
In this work, the relative error is rather used to test the eigenpairs:

ϵ KS rel = ∥C V m y k -λ k V m y k ∥ 2 |λ k ∥ V m y k ∥ 2 = | bH m y k | |λ k |∥ V m y k ∥ 2 . (II.74)
Algorithm 2 KS algorithm Input: C, v 1 unitary, m and l ≤ m the number of desired eigenpairs Output:

V j , H j such that CV j = V j H j , with j ≤ l 1: initialization of V m ← v 1 0 , k ← 0, p ← 0 2: while j ≤ l do 3:
m -p Arnoldi steps with deflation:

C V m = V m H m + v m+1 bH m 4:
Schur decomposition with sorting

H m ← W H 2 W H 1 H m W 1 W 2 5: V m ← V m W 1 W 2 6:
for k ∈ j, l do

7: compute (λ k , y k ) verifying H m y k = λ k y k and set ϵ KS rel = | bH m y k |/(|λ k |∥ V m y k ∥ 2 ) 8: if ϵ KS rel < tol then 9:
lock (λ k , y k ) and continue The Ritz value is considered as converged if ϵ KS rel is lower than the user-defined tolerance. When a Ritz value is converged, it should no longer participate in the search for other eigenvalues. This is the goal of deflation. Supposing that λ 1 is the Ritz value with the largest magnitude of T m , the decomposition (II.69) can also be writen as:

C ṽ1 V m-1 = ṽ1 V m-1 λ 1 ⋆ 0 T m-1 + v m+1 b 1 bH m-1 . (II.75)
If λ 1 is converged, the coefficient b 1 can be set to 0, and all subsequent operations can be managed exclusively on the order-(m -1) decomposition

C V m-1 = V m-1 T m-1 + v m+1 bH m-1 . (II.76)
The next Ritz values are then tested and when the test fails, the loop breaks and the decomposition is truncated. Finally, if l Ritz values are converged, we get C V j = V j T j , where the approximate eigenvalues are on the diagonal of T m , and the Ritz pairs of T m are approximate eigenpairs of C. Algorithm 2 gives a synthesized implementation of the KS method.

II.6.5 Shift-and-invert transformation

As stated earlier, the Schur matrix resulting from the KS basis expansion is sorted so that the Ritz values of largest magnitude are located at its top left corner. Other critera such as sorting by lowest magnitudes could have also been chosen, but the KS algorithm performs poorly when computing eigenvalues located in the interior of the spectrum [START_REF] Hernández | Krylov-Schur Methods in SLEPc[END_REF]. In global modal stability, the KS algorithm is generally used in conjunction with the spectral shift-and-invert transformation [START_REF] Theofilis | Global Linear Instability[END_REF][START_REF] Schmid | Stability and Transition in Shear Flows[END_REF] already mentioned in previous sections. A shift s ∈ C is selected and the transformed eigenvalue problem involving the shit-and-invert operator (C -sI N ) -1 is solved. Both operators C and (C -sI N ) -1 have the same eigenvectors but the eigenvalues of largest magnitude of (C -sI N ) -1 are directly linked to the eigenvalues of C closest to s. The concept is schematically illustrated in fig. II.7, where the spectrum of a fictive operator C is compared with that of the shift-and-invert operator. Only one unstable eigenvalue, equal to 1, is identified. The transformation, carried out with s = 0.9, clearly shows that the eigenvalues have been separated and the unstable eigenvalue has become the one with the largest magnitude. The counterpart of such a technique is the many inversions of complex linear systems needed during the formation of the KS basis, because of the application of (C -sI N ) -1 on the successive Krylov vectors.

-20 -10 0 10 ℜ(λ) -5 -2.5 0 2.5 5 ℑ(λ) λ (λ -0.9) -1

II.7 Conclusion

This chapter provides the mathematical background for both the numerical discretization of operators in compressible CFD and the overall global stability analysis process. The first-order implicit backward Euler method is used to compute a fixed point of the Burgers equation with high precision, matching the theoretical solution. The stability analysis is then applied to the fixed point, which is shown to be unconditionally stable by modal analysis. The optimal gain computed with resolvent analysis shows a peak at zero frequency, resulting in a steady optimal forcing. The KS algorithm is also introduced in its complex arithmetic version for the sake of introduction.

The main challenge in scaling up global stability analysis for industrial configurations is dealing with extremely large linear systems that are too big to be factored using LU factorization due to memory constraints. The next chapter presents a strategy for conducting global stability analysis on large-scale systems using elsA. Based on previous work by Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF], inversions in elsA will be managed using the iterative GMRES algorithm [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. To maintain HPC optimizations, the linear systems of interest will be reformulated in the real-equivalent K formulation.

Chapter III

Scaling up of the method to industrial cases

III.1 Introduction

Chapter II was devoted to the presentation of the general concepts of global stability analysis and gave an overview of the proposed method used in this thesis. Burgers equation has been used to illustrate fixed point computation as well as linear, adjoint, and resolvent global stability analysis. This chapter gives details about the proposed method to deal with large CFD problems (formulated with Euler, NS, or RANS equations), as implemented in elsA software. Because of the number of degrees of freedom (DOF) of such configurations, the linear systems involved in the analyses may become too large to solve with direct methods. Iterative methods must be used instead, and we focus in particular on the Generalized Minimal Residual (GMRES) algorithm introduced by Saad in [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

III.1.1 Scale up issues

We consider a 3-D RANS case (d = 6) with N c = 10 × 10 6 cells, corresponding to N = 60 × 10 6 complex-valued DOF for the global stability problems. The computation is run on a computer cluster composed of CPUs like Cascade Lake processors (representative of ONERA's clusters), organized in nodes, as shown in fig. III.1. In this case, each node has 28 CPUs, with 4.57 GB/CPU and a total memory of 128 GB. Global stability problems require the inversion of complex operators, which have the same size as the Jacobian matrix, and the number of entries can be estimated as NE(J ) ≃ stencil × d 2 × N c . Using a second-order finite volume discretization 7 , storage of the Jacobian matrix with complex floating-point numbers in double precision requires about 720 GB for the representative case.

We first focus on the memory footprint of the LU factorization of J , written LU J . Based on the work of L'Excellent [START_REF] Excellent | Multifrontal Methods: Parallelism, Memory Usage and Numerical Aspects[END_REF], the number of entries of LU J may be estimated with the following relation:

NE (LU J ) ≃ NE(J ) 7 4/3 . (III.1)
For the considered case, this results in about 500 TB of memory and would requires more than 100 × 10 3 CPU. Most of the existing calculators do not have this amount of memory available during the computation unless resorting to an overly large number of CPUs in parallel. As an indication, the maximum number of CPUs allowed at ONERA for one job is 4788, representing 171 nodes. Instead of such direct methods, the inversions might be addressed using preconditioned iterative Krylov methods (KM) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. With such approach, preconditioners can be based on ILU(k) or local LU factorization built from the Jacobian matrix partitioned on each CPU. If n KM is the 7 for which the approximate size of the numerical stencil is stencil = 5 3 in elsA.

Cluster

Method

Memory

N CPU LU ∽ 500 TB > 100000 KM-ILU(k) ∽ 2.5 TB ∽ 500 KM-local LU ∽ 30 TB ∽ 6000
Table III.1: Comparison of memory costs between direct and preconditioned iterative methods for the considered case (N = 60 × 10 6 ).

number of Krylov vectors involved in the computation and γ(k) a factor depending on the number of fill-in8 , the number of entry of the KM-ILU(k) method is approximated by:

NE(KM-ILU(k)) ≃ γ(k) × NE(J ) + n KM × N. (III.2)
For the considered case, 2.5 TB of memory is needed and the computation can be launched on 500 CPUs. In this case, additional storages might be required if the convergence has to be improved using Restricted Additive Schwarz (RAS) methods [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. Nevertheless, such a cost is much more affordable than that of direct methods. For particularly stiff cases, another approach based on local LU factorization can be used. The number of entries in KM-local LU method is approximated by:

NE(KM-local LU) ≃ N -1/3 CPU × NE(J ) 7 4/3 + n KM × N. (III.3)
This represents about 30 TB of memory and would requires 6000 processors for the computation, considering the representative case. This is still being too large for ONERA's clusters, which is why hybrid preconditioners like Hybrid Local LU (HLU) factorization have been developed (see III.3.4 for more details). Table III.1 summarizes the costs of each method. Based on these observations, the preconditioned KM suit well to circumvent memory issues resulting from direct factorizations. Another important issue that must be addressed is the complexity of the arithmetic used to construct and invert the operators of global stability, which is poorly vectorizable compared to real arithmetic. This prevents optimal use of HPC optimizations. Moreover, Guilbert's previous work [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] introduced powerful preconditioning tools to efficiently solve real linear systems arising in the fixed point computation (see II.31) in the context of compressible CFD. Taking these considerations into account, we choose to reformulate the non-vectorizable complex linear systems of global stability into a real-equivalent formulation that preserves HPC optimizations and reuses the existing preconditioning tools developed in [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. This is achieved in this thesis by using the K formulation introduced by Day and Heroux [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF] that preserves the sparsity pattern of the original systems. Some of its properties and the K-formulated KS algorithm that has been implemented in elsA are detailed in this chapter.

III.1.2 Numerical set-up of laminar cylinder flow

Each concept will be illustrated using the canonical example of laminar cylinder flow. This example is used to validate the proposed algorithms as it has been widely documented (see [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Marquet | Stabilité globale et contrôle d'écoulements de recirculation[END_REF]). The numerical domain, shown in of the cylinder as the reference length, L, and the free-stream density, velocity, and temperature. The other physical quantities are calculated using the following expressions:

a ∞ = 1 M ∞ , (III.4) p ∞ = 1 γM ∞ , (III.5) e ∞ = 1 γ(γ -1)M 2 ∞ = c v , (III.6) R = 1 γM 2 ∞ , (III.7) µ ∞ = 1 Re ∞ . (III.8)
A second-order Roe scheme with MUSCL reconstruction and no limiter is used for convective fluxes, while diffusive fluxes are computed with a 5p cor scheme (see. III.3.2). The boundary conditions are no-slip at the cylinder wall, non-reflective at the farfield, and symmetry conditions on each side of the domain to produce a 2-D-like computation. Despite being inherently 2-D, this case will be treated as a 3-D computation by the elsA code.

III.2 Equations of compressible CFD

The equations solved in elsA are based on the compressible NS equations, which describe the behavior of a fluid flow in terms of mass, momentum, and energy conservation. This section provides the mathematical formulation of these equations in a general form for both laminar and turbulent flows.

III.2.1 Navier-Stokes equations

A reference frame where the space-time coordinate of each point is identified by (x, y, z, t) is considered. The state vector of the conservative variables of the flow, w = ρ ρu x ρu y ρu z ρE T , involves the density field ρ, the absolute velocity field u = u x u y u z T and the total energy per unit of mass E, which is the sum of the internal energy e and the kinetic energy ∥u∥ 2 2 /2. The set of NS equations governs the evolution in space and time of w and writes:

                 ∂ρ ∂t + ∇ • ρu = 0 ∂ρu ∂t + ∇ (ρu ⊗ u + pI 3 -τ ) = 0 ∂ρE ∂t + ∇ • (ρEu + pu -τ • u + q) = 0.
(III.9) with p, τ and q respectively standing for the static pressure, the viscous stress tensor and the heat flux vector. Using the notation of II.2.1, the convective and diffusive parts of (III.9) are the following:

f c =         ρu x ρu 2 ρu 3 ρu 2 x + p ρu x u y ρu x u z ρu x u y ρu 2 y + p ρu y u z ρu x u z ρu y u z ρu 2 z + p ρu x (E + p) ρu y (E + p) ρu z (E + p)         ,
(III.10)

f d =         0 0 0 τ xx τ xy τ xz τ xy τ yy τ yz τ xz τ yz τ zz u x τ xx + u y τ xy + u z τ xz -q 1 u x τ xy + u y τ yy + u z τ yz -q 2 u x τ xz + u y τ yz + u z τ zz -q 3         .
(III.11) Expressions of τ and q have to be explicited with behavior laws. The fluid is assumed to be Newtonian, meaning that the tensor τ which models the internal frictions in the fluid can be expressed as:

τ = µ 2 ∇u + ∇u T - 2 3 ∇ • uI 3 , (III.12)
with µ is the dynamic viscosity of the fluid. Fourier law is used to describe the heat flux vector:

q = -κ∇T, (III.13)
with κ the thermal conductivity coefficient and T the absolute static temperature.

From now, we have a system of five equations for the seven unknowns ρ, u x , u y , u z , E, p and T . This means that two additional relations must be provided to close the problem. In this work, only perfect gas in thermodynamic equilibrium, not chemically reacting and with constant isochoric and isobaric specific heats c v and c p , are considered. The specific heat ratio is γ = c p /c v = 1.4. Such assumptions lead to the following expressions for e and p:

e = c v T (III.14) p = ρRT, (III.15)
with R the ratio of the perfect gas constant to the molar mass of the fluid, set at 287.053 J kg -1 K -1 for air. For perfect gas, the viscosity coefficient µ is considered to be independent of the pressure but strongly depends on the temperature. To this end, the Sutherland law [START_REF] Sutherland | LII. The viscosity of gases and molecular force[END_REF] is used:

µ = µ S T T ref 1 + T S /T ref 1 + T S /T , (III.16)
with µ S = 1.711 × 10 -5 N m -2 s the dynamic viscosity coefficient of air at temperature T ref = 273 K, and T S = 110.4 K the Sutherland temperature. As we mainly deal with gas, the thermal coefficient κ is proportional to µ according to the law:

κ = c p µ Pr , (III.17)
where Pr is the Prandtl number [START_REF] Prandtl | Applications of modern hydrodynamics to aeronautics[END_REF], assumed to be constant and equal to 0.72.

III.2.2 Reynolds Averaged Navier-Stokes equations

In this work, turbulence is modeled using a statistical approach in which each field f is decomposed into a mean field f , and a fluctuating one f ′ , so that f (x, t) = f + f ′ . When f slowly varies in time, it is assumed to be ergodic [START_REF] Walters | An Introduction to Ergodic Theory[END_REF], and its average value can be obtained from the temporal average on a certain period of time T :

f = 1 T ˆT/2 -T /2 f (x, t + t ′ )dt ′ . (III.18)
When f has a variable density, the use of Favre average [START_REF] Favre | The equations of compressible turbulent gases[END_REF] is prefered, giving the decomposition f (x, t) = f + f ′′ in which the mean field writes:

f = ρf ρ . (III.19)
In NS equations, the temporal average is used for variables ρ, p, τ and q while the Favre average is used for u and E. This results into the following set of RANS equations:

                   ∂ ρ ∂t + ∇ • ρ u = 0 ∂ ρ u ∂t + ∇ ρ u ⊗ u + pI 3 -τ -τ r = 0 ∂ ρ E ∂t + ∇ • ρ E + k u + p u -τ • u + q + q t = 0. (III.20)
Three new terms appear in this last system: the kinetic energy of the turbulent movement k = ρu ′′2 /2 ρ, the Reynolds tensor τ r = -ρu ′′ ⊗ u ′′ and the turbulent heat flux

q t = ρu ′′ h ′′ .
The problem is closed providing a turbulence model, which we choose to follow the Boussinesq hypothesis. It states, by analogy with momentum transfers resulting in molecular viscosity µ, that momentum transfers resulting from turbulent eddies are modeled with an eddy viscosity µ t . The Reynolds stress tensor and the turbulent heat flux vector can thus be recast as:

τ r = - 2 3 ( ρk + µ t ∇ • u) I 3 + 2µ t ∇ u + ∇ u T (III.21) q t = - c p µ t Pr t ∇T, (III.22)
where Pr t is the turbulent Prandtl number, set at 0.9 for air. With such hypothesis, we can therefore introduce the following starred quantities:

                           E ⋆ = E + k p ⋆ = p + 2 3 ρk τ ⋆ = 1 + µ t µ τ q ⋆ = 1 + µ t µ Pr Pr t q (III.23)
and replace in equation (III.20), their corresponding non-starred quantities, leading to the RANS equations in a similar form as the NS equations.

The one-equation Spalart-Allmaras (SA) model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] is commonly used in compressible RANS simulations for aerodynamic. It has been designed for aerodynamic applications involving wallbounded flows. It solves a transport equation for a kinematic-like viscosity coefficient ν that is linked to the turbulent dynamic viscosity via:

µ t = ρ νf v1 , (III.24)
where f v1 is the following damping function:

f v1 = χ 3 χ 3 + c v1 , (III.25)
with χ = ν/ν, ν being the molecular kinematic viscosity (ν = µ/ρ), and c v1 = 7.1. The conservative form of the added transport equation governs the evolution of the quantity ρ ν and writes:

∂ ρ ν ∂t +∇• ρ ν u - 1 σ SA ρ (ν + ν) ∇ ν = c b1 (1 -f t2 ) S ρ ν + c b2 σ SA ∇ ν • ∇ ρ ν production -ρ c w1 f w - c b1 κ 2 SA f t2 ν d w 2 destruction .
(III.26) Unlike the other considered equations, this one has a source term on its right-hand side representing a balance between a production and a destruction terms. The various quantities of the model are then provided by:

f t2 = c t3 e -c t4 χ 2 ,
(III.27)

f w = g 1 + c 6 w3 g 6 + c 6 w3 1/6
, (III.28)

g = r + c w2 r 6 -r , r = min ν Sκ 2 SA d 2 w
, 10 (III.29)

S = S + ν κ 2 SA d 2 w f v2 , S = 2 Ω : Ω, Ω = 1 2 ∇ u -∇ u T , (III.30) f v2 = 1 - χ 1 + χf v1 , (III.31)
where d w stands for the distance from the field point to the nearest wall. The constants of the model are:

σ SA = 2 3 , κ SA = 0.41 (III.32) c b1 = 0.1355, c b2 = 0.622, (III.33) c w1 = c b1 κ 2 SA + 1 + c b2 σ SA , (III.34) c w2 = 0.3, c w3 = 2, (III.35) c t3 = 1.2, c t4 = 0.5. (III.36) (III.37)
This original SA model requires ν to be always positive, which is ensured by taking its absolute value. Since this function cann't be differentiated at zero, we have preferred to use the negative SA model [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF] which is a slight modification of the standard model. When ν ≥ 0, equation (III.26) is solved. However, when it becomes negative, the following equation is solved instead:

∂ ρ ν ∂t +∇• ρ ν u - 1 σ SA ρ (ν + νf n ) ∇ ν = c b1 (1 -c t3 ) S ρ ν+ c b2 σ SA ∇ ν•∇ ρ ν+ ρc w1 f w ν d w 2 , (III.38)
with

f n = c n1 + χ 3 c n+1 -χ 3 , c n1 = 16. (III.39)
Next section first provides a brief overview of the numerical discretization techniques utilized in elsA. It then describes the two fixed point computation methods used and offers insight on preconditioning and multicoloring.

III.3 elsA CFD solver

The numerical methods used in elsA for computing steady solutions of the discretized NS equations are described below. The focus is on implicit methods that solve (approximately or exactly) the implicit backward Euler (see equation (II.31)), which are preferred over explicit methods due to their unconditional stability and ability to handle larger time-steps. However, these methods require solving linear systems, which can be computationally expensive.

In some cases, such as when the initial field is uniform, using an exact inversion method at each step can be counterproductive, especially during the first iterations. To address this issue, less efficient but faster methods, like matrix-free ones [START_REF] Corre | Contribution à la simulation et à l'analyse des écoulements compressibles[END_REF] or the so-called LU-Symmetric Successive Overrelaxation Scheme (LU-SSOR) [START_REF] Yoon | An LU-SSOR scheme for the Euler and Navier-Stokes equations[END_REF] that solves an approximation of the linear system can be used to refine the initial uniform field towards a fixed point, reducing the need for costly matrix factorizations during the early iterations. After this, a more accurate computation is performed using the GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] algorithm to invert the exact linear systems and decrease the residual below 10 × 10 -10 . These methods will be described in detail in the following section.

III.3.1 Unstructured computation

In this thesis, only unstructured computational domain are considered. The mesh is contained in a CGNS file (https://www.cgns.org) and, as discussed by Bourasseau in [START_REF] Bourasseau | Contribution à une méthode de raffinement de maillage basée sur le vecteur adjoint pour le calcul de fonctions aérodynamiques[END_REF], elsA, for an unstructured block, performs the calculations by iterating through all the faces of the mesh. The faces are ordered by starting with the internal faces and ending with the boundary faces, grouping both the connection faces and the physical boundary faces to which a physical boundary condition is imposed. However, since the boundary faces are not ordered by physical or connection conditions, elsA generates two lists grouping, for one, all physical boundary faces and, for the other, all connection faces. Each of these lists contains lists whose length is equal to the number of faces of the boundary condition it represents. It is therefore necessary to create, for each list, an indirection table to go from the local numbering of the face corresponding to its position in the list to the global numbering of that face. The correspondence between the faces and the cells is handled by building two indirection tables to determine for each face its right and left cells, knowing that the definition of right and left is imposed by the normal to the face, itself imposed by the ordering of the nodes that make up this face.

III.3.2 Numerical discretization

To describe the numerical discretization as employed in elsA, we rely on the general form of the governing equations of compressible CFD and their discretization with finite volumes detailed in II.12. We work on a flow domain Ω that is partitioned into an unstructured grid containing N c cells, as illustrated in fig. III.3. The system (III.9) is integrated over each cell Ω k , and the finite volume method is used to compute the solution as in equation (II.12). Approximate Riemann solver like Roe's (see equation (II.13)) or HLLC [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] are used to calculate the function h c at each interface.

As for equation (II.20), the diffusive fluxes are computed with a discretization centered at the faces. The main difficulty results from the computation of the gradients of the primitive variables. The 5p method (for five-point stencil) is a simple arithmetic mean of gradients computed at cell centers of both sides of the face I j . Referring to the scheme of fig. III.3, it writes:

∇p 5p I j = 1 2 (∇p| k + ∇p| l ) . (III.40)
However, this expression is modified to get a so-called 5p cor method to improve the discretization quality and avoid oscillations. In the latter, the gradient is computed as if the control volume was centered on I j and composed of cells Ω k and Ω l , providing the following new expression: For RANS configurations, the turbulent source terms are discretized using a second-order integration method.

∇p 5pcor I j = ∇p 5p I j -∇p 5p I j • x k -x l ∥x k -x l ∥ 2 - p k -p l ∥x k -x l ∥ 2 x k -x l ∥x k -x l ∥ 2 . (III.41) Ω k Ω l I j x c k

III.3.3 LU-SGS method for fixed point computation

Iterative methods such as matrix-free, LU-SGS, or LU-SSOR [START_REF] Corre | Contribution à la simulation et à l'analyse des écoulements compressibles[END_REF][START_REF] Outtier | A new dynamic code architecture for CFD computations : application to the development of an overset-grid compact high-order solver for compressible aerodynamics[END_REF] are derived by approximating the implicit operator A = ∆t -1 I N + J to solve the system in equation (II.31). This system is constructed by discretizing the convective component using the Rusanov flux and finding the spectral radius of the resulting operator, which leads to a simple tridiagonal operator A. During each nonlinear iteration, a sub-iteration process is used to calculate ∆W n and the system is solved using the SSOR method [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Yoon | An LU-SSOR scheme for the Euler and Navier-Stokes equations[END_REF]. As pointed out in [START_REF] Weber | Développement de méthodes implicites pour les équations de Navier-Stokes moyennées et la simulation des grandes échelles : application à l'aérodynamique externe[END_REF], the best results are achieved when the relaxation parameter is set to 1, leading to the symmetric Gauss-Seidel (SGS) method. Considering L l , D l and U l being the strictly lower part, the diagonal part, and the strictly upper part of A(W n,l ), the LU-SGS method is obtained by factorizing A(W n,l ) as:

A(W n,l ) = L l + D l + U l ≃ (L l + D l ) D -1 l (U l + D l ) , (III.42)
where L l D -1 l U l is the error committed relatively to the original factorization. The resolution is carried out by the three following steps:

         (U l + D l ) ∆W l+1/3 = R(W n,l ) D l ∆W n,l+2/3 = ∆W n,l+1/3 (L l + D l ) ∆W n,l+1 = ∆W n,l+2/3 , (III.43)
where R stands for the residual associated with the approximated implicit system.

Inversions of the approximate tridiagonal operator at each nonlinear iteration is much less expansive than that of the exact operator, resulting in a much faster computation. It is especially fast for low CFL numbers up to CFL = 100, when the goal is to filter coarse unsteadiness, i.e. decreasing the residual ∥R(W n )∥ 2 to below 10 -3 to 10 -4 (see [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]). However, when the goal is to reach convergence around machine precision, the CFL number must increase to recover the Newton method that achieves a quadratic convergence, and the exact linear system must be solved.

III.3.4 GMRES algorithm for fixed point computation

To increase the accuracy of a computed a fixed point of the governing equations, the implicit backward Euler scheme (II.31) with adaptative CFL must be solved at each nonlinear iteration. This requires the multiple iversions of A during the computation. To this end, given the sizes of 3-D configurations, the use of iterative methods instead of direct ones is mandatory, due to memory limitations. Additionally, as noted by Guilbert in [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF], due to the non-symmetrical, non-normal, and potentially non-diagonal dominant nature of the operator, the focus is placed on the GMRES algorithm introduced by Saad & Schultz in [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

Derivation of GMRES algorithm. Considering a non-singular matrix A ∈ R N ×R N , and a righthand side vector b ∈ R N , we focus on the resolution of the following real-valued linear system:

Ax = b, (III.44)
where x ∈ R N is the unknown. The GMRES algorithm computes a solution to equation (III.44) by projecting x onto a given Krylov space. Until convergence, the Krylov basis is continually improved. The steps of the GMRES algorithm are described hereafter. Given an initial approximation x 0 of x, the initial residual r 0 is defined as r 0 = b -Ax 0 . The first step is to produce an orthonormal basis of the Krylov space K m (A, r 0 ) associated with r 0 , denoted as K m in the following. As in II.6.1, this is achieved with the Arnoldi algorithm 1 which produces the following Krylov decomposition:

AV m = V m H m + v m+1 h T m,m+1 , (III.45)
which can be written in a more compact form as:

AV m = V m+1 H m . (III.46)
Here, V m+1 and H m are N -by-(m + 1) and (m + 1)-by-m real matrices, respectively, and β = ∥r 0 ∥ 2 . The Arnoldi algorithm is initialized with the unitary vector v 1 = r 0 /β. The matrix-vector products required to build the decomposition are computed using AD, which is optimal to benefit from all HPC optimizations (see II.4). The next step is to look for a vector

x m ∈ x 0 +K m that minimizes ∥r m ∥ 2 . As (v k , 0 ≤ k ≤ m-1) is an orthonormal basis of K m , there exists a vector y m ∈ R m verifying x m = x 0 + V m y m . The residual r m = b -Ax m can thus be recast as: b -Ax m = r 0 -AV m y m = r 0 -V m+1 H m y m = V m+1 V T m+1 βv 1 -H m y m , (III.47)
which becomes, after taking the L 2 -norm and using the orthogonality of the Krylov vectors:

∥b -Ax m ∥ 2 = ∥βe 1 -H m y m ∥ 2 . (III.48)
Minimizing the left-hand side of the latter equality is hence equivalent to solve following size m least-square problem:

x m = argmin y ∥βe 1 -H m y∥ 2 . (III.49)
This can be efficiently achieved using the Givens rotations [START_REF] Givens | Computation of Plain Unitary Rotations Transforming a General Matrix to Triangular Form[END_REF]. This technique has the additional advantage of computing the relative residual ϵ GMRES rel = ∥r m ∥/β implicitly, i.e. without requiring an explicit calculation of the approximate solution x m . The approximation is computed once when the relative residual becomes lower than the user-specified tolerance. for l ∈ 1, k do 7:

h lk ← ⟨Av k , v l ⟩ 8: v k+1 ← w k -v l h lk 9:
end for

y m = argmin y ∥βe 1 -H m y∥ 2 20: x m = x 0 + V m y m
When using the GMRES algorithm to compute ∆W n , it can be counterproductive to set the tolerance too low, as it may prevent the adaptive CFL from increasing as it should. Through trial and error, a tolerance of 10 -4 has proven to be a good compromise. In practice, this algorithm, in its current state, is not expected to converge well, as the operator A is generally ill-conditioned. In fact, the solvability of system (III.44) depends on the condition number of A, which, given a convenient norm ∥ • ∥, is defined as:

κ(A) = ∥A -1 ∥∥A∥.
(III.50)

In our systems of interest, the non-normality of A is responsible for increasing κ(A), and the larger it is, the more complicated the resolution is. This is why we need to slightly modify this system to incorporate operators called preconditioners, which aim at lowering κ(A).

Preconditioned GMRES. Preconditioning is a modification of a linear system that makes it easier to solve using an iterative method. Although not guaranteed, the resulting system may need fewer steps to converge. Concretely, the goal is to find an operator M such that M ≃ A and apply it on the left (equation (III.51)) or right (equation (III.52)) of the system (III.44):

M -1 Ax = M -1 b AM -1 Mx = b. (III.51) (III.52)
As discussed in [START_REF] Aditi | A comparison of preconditioned krylov subspace methods for large-scale nonsymmetric linear systems[END_REF], the asymptotic behavior of an iterative method is the same whether we use left or right-preconditioning. However, in the left-handed case, the preconditioning also applies

v k+1 ← AM -1 v k 6:
for l ∈ 1, k do 7:

h lk ← ⟨Av k , v l ⟩ 8: v j+1 ← w j -v k h kj 9:
end for

y m = argmin y ∥βe 1 -H m y∥ 2 20: x m = x 0 + M -1 V m y m
to the residual making it different from that of the original problem. This might result in significant differences between the original residual and the left-preconditioned one, as the calculation proceeds. Instead, we choose to work with right-preconditioning as it leaves the residual vector unchanged and often shows better convergence:

r = b -Ax = b -AM -1 Mx. (III.53)
With right-preconditioning, algorithm 3 becomes algorithm 4. In this version of GMRES, the preconditioner remains the same at each iteration, and the solution is recovered at the end of the process by pre-multiplying V m y m with M -1 . Another version that is popular is CFD is the so-called Flexible GMRES (FGMRES) algorithm [START_REF] Giraud | A Set of Flexible GMRES Routines for Real and Complex Arithmetics on High Performance Computers[END_REF] in which the preconditioner is updated at each step. To this end, two GMRES are nested one into the other. At each step of the external GMRES, an internal one computes the intermediate vector u k+1 with:

u k+1 = M -1 k v k , (III.
u k+1 ← M -1 j v k 7: v k+1 ← Au k+1 8:
for l ∈ 1, k do 9:

h lk ← ⟨Av k , v l ⟩ 10: v k+1 ← w k -v l h lk 11:
end for

v k+1 ← v k+1 /h k+1,k 19:
end if 20: end for 21: y m = argmin y ∥βe 1 -H m y∥ 2 22: x m = x 0 + U m y m external Krylov bases must thus be carefully chosen to be compliant with memory limits. In each case, the preconditioning is initially built upon the Jacobian matrix. In the next paragraph, the procedure employed in elsA to compute J is described.

Formation and storage of J . A review of existing preconditioners is exposed in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Simple ones arise from classical iterative techniques like Jacobi or Gauss-Seidel. However, they are known to show a poor convergence rate when used in the context of compressible CFD. The so-called matrix-free approaches are really cheap in terms of both memory and CPU usage, but the most efficient among them requires a minima the formation of a good approximation of A. When addressing compressible CFD problems, Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] showed that the most convenient way is to precondition from the exact matrix. We recall that the matrices involved in CFD problems, like J , usually are sparse matrices, i.e. matrices in which most of the elements are zero, as shown in fig. III.4.

The first step of preconditioning is therefore to compute the exact operator A, and in particular the Jacobian matrix J . In chapter I, the columns of J are calculated by evaluating successive matrix-vector products with AD on each vector of the canonical base. However, when dealing with larger 2-D or 3-D configurations, meshes could reach millions of cells, and such a strategy becomes unaffordable. In elsA, this issue is addressed by using a multicoloring technique. Full details about this approach can be found in [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. To summarize, the idea is first to identify structurally orthogonal columns, i.e. columns that have no common element in the same line, and then to produce a probing vector through which a Jacobian-vector product returns the identified columns.

Multicoloring can be illustrated with the following schematic sparse matrix whose structurally orthogonal columns has been colored the same way:

A =                . (III.56)
Introducing the three following probing vectors:

t 1 =           1 0 1 0 0 0           , t 2 =           0 1 0 1 0 0          
, and

t 3 =           0 0 0 0 1 1           , (III.57)
the entire matrix can be recovered evaluating At 1 , At 2 and At 3 . The matrix being built by blocks, the final number of matrix-vector products to compute is N color × d, which also depends on both the discretization order and the mesh topology.

Several storage format are available to store sparse matrices and most of them are presented in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. In elsA, A is stored in the Block compressed Sparse Row (BSR) format, that is the block variation of the Compressed Sparse Row (CSR) format. The CSR storage of a matrix A is composed of three tables AA(A), JA(A) and IA(A). The first one contains every nonzero coefficient, the second their column index, and the third provides, for a given index p, the position in AA(A) of the first nonzero coefficient of line p. An example of CSR storage for a 3-by-3 matrix is provided hereafter: (III.58)

A =   
The BSR format is identical except that the coefficients a qr of A are d-by-d blocks instead of scalar quantities.

Preconditioners. This thesis builds upon the previous developments made by Guilbert in elsA, which give access to a large variety of preconditioning tools based on domain decomposition. Incomplete LU factorizations [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] like ILU(0) and ILU(k), where k is the level of fill-in, coupled with FGMRES, are sufficient for a large variety of CFD problems. The quality of those preconditioners may be enhanced by distributed computing extension via the Restricted Additive Schwarz (RAS) methods. When the problem is particularly stiff, like in RANS computations, more sophisticated preconditioning tools must be used. Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] developed an hybrid method based on the mesh partitioning, that identifies the critical zones through the "LU-ILU distance" criterion, and decides whether ILU(0), ILU(k) or local LU factorization (see https://www.pardiso-project.org/) are performed. In the most critical zones, the preconditioner is a local LU factorization, and userdefined incomplete factorizations in the remaining zones. By allowing each domain to be preconditioned differently, this hybrid approach may yield important gain over a classical local LU decomposition. This type of preconditioning will be denoted HLU factorization in the following. Finding a good preconditioner for large sparse systems in CFD remains a challenging task and is still an important research topic at ONERA with for instance the work of Jang [START_REF] Jang | Randomized flexible GMRES with Singular Vectors Based Deflated Restarting[END_REF].

The different variations of preconditioning developed by Guilbert in elsA have been tested on laminar cylinder case. To test and compare the convergence of FGMRES with the various available preconditioners, the focus was made on the resolution of the system J z = b, where J is the Jacobian matrix of a fixed point calculated at Re = 47, and b is some random vector. This system is a particular case of equation (III.44) when CFL → +∞, i.e. when the fixed point is reached and the time term vanishes. At the beginning of a fixed point computation, this term is large and increases the diagonal dominance of A, which makes the system easier to solve. On the contrary, the closer the fixed point is, the stiffer the system becomes, resulting in a more challenging resolution. In the experiments, both internal and external Krylov basis were chosen to be composed of 60 vectors, and ϵ GMRES rel has been set to 10 -14 while the inner GMRES tolerance is set to 10 -3 . Fig. III.5 presents the convergence results using ILU(1), ILU(2), ILU(3), and HLU and table III.3 shows the required memory and elapsed time for each computation. HLU is the best in term of CPU time but also has the highest memory requirement. This is why in practice, we prefer to use ILU [START_REF] Greitzer | Surge and Rotating Stall in Axial Flow Compressors-Part II: Experimental Results and Comparison With Theory[END_REF] preconditioning which provides good performances in both CPU time and memory requirement. HLU will be used for highly stiff problems.

III.3.5 Fixed point calculation of the laminar cylinder flow

Fixed point computation on the laminar cylinder case have been performed by solving the implicit backward Euler scheme presented in II.3.2 with adaptative CFL [START_REF] Crivellini | An implicit matrix-free Discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations[END_REF] and using FGMRES-ILU(3) for inversions. The calculations were started with a uniform field. Considering ∆l k the characteristic length of a cell Ω k , the local time-step in elsA is defined as follows: This section was dedicated to the presentation of the techniques selected in this thesis to accurately compute fixed points with elsA. The GMRES algorithm as well as the preconditioning techniques, as implemented to address inversions of the real-valued linear systems involved in fixed point computation, were detailed. The next natural step is to wonder how these tools could be adapted to solve the complex linear systems involved in global stability analyses. Next section presents the strategy developed in elsA to do so, and focus especially on the use of the K formulation [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF].

∆t k = CFL min ∆l k ∥U k ∥ 2 + a , ρPr∆l 2 k 2γµ k (III.

III.4 K formulation and some of its properties

The tools of fixed point computation developed by Guilbert operate in real arithmetics, which differs from the global stability problems that are complex-valued. The code must therefore be thoroughly rewritten to be adapted to solve complex linear systems. Moreover, one of the main component of the HPC optimizations implemented for the fixed point computation is the vectorization of loop patterns. Under some conditions like the independence of each operation within a loop (see [START_REF] Lecas | Vectorisation SIMD[END_REF] for more details), this provides a significant speed-up during the computation [START_REF] Maugars | Algorithmic Differentiation for an effcient CFD solver[END_REF]. For instance, if a given loop satisfies the vectorization conditions, and the optimization packs the indices five by five, this means that all subsequent operations are managed for the five indices at the same time. In practice, a speed-up factor of three is approximately observed when using vectorization. Another mandatory condition for the vectorization is the memory contiguity of the floats composing the input vector. However, in a complex vector, each component consist of a tuple of real and imaginary parts. From a memory point of view, the complex vector is stored in an alternance of real and imaginary components, and the contiguity is lost9 . The use of complex arithmetics thus breaks the vectorization and consequently the HPC optimizations. To avoid such issue, we choose in this work to reformulate the complex systems involved in global stability in a real-equivalent version. Among all available formulations, the K formulation introduced by Day & Heroux [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF] has been selected. The justification of this choice, as well as some of its mathematical properties are provided in this section.

III.4.1 Real equivalent formulation of a complex problem

For the sake of clarity, the focus is made on the general complex system under the following form:

(δB + γI N )z = t, (III.60) with z, t ∈ C N , γ ∈ C, δ = ±1, and where B stands either for J or J T . We introduce the sparse matrix C = (δB + γI N ). By decomposing each vector and matrix in their real and imaginary part, and by respectively introducing D and E as the real and imaginary parts of C, the system (III.60) becomes:

(D + iE)(z ℜ + iz ℑ ) = (t ℜ + it ℑ ),
which can be recast in the following form:

D -E E D =C K1 z ℜ z ℑ = t ℜ t ℑ . (III.61)
This so-called K1 formulation provides an equivalent real alternative to the complex system where the operator

C ∈ C N × C N has been reformulated into operator C K1 ∈ C 2N × C 2N .
With such a modification, the vectorization is preserved however, the sparsity pattern of the original system is not. This is problematic since the preconditioners developed by Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] are built from the sparsity pattern of the considered operator.

For this reason, we choose instead to focus on the equivalent K formulation [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF], in which each coefficient c qr = d qr + ie qr is replaced by a 2-by-2 block entry c K qr :

c K qr =
d qr -e qr e qr d qr .

(III.62)

The K-formulated matrix C in the K formulation will be noted C K . This mapping is an homomorphism [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF], i.e. it satisfies the following condition:

(i) I K N = I 2N (III.63) (ii) (X Y) K = X K Y K , ∀X , Y ∈ C N × C N . (III.64)
For example, the K formulation of the matrix is the matrix As detailed in II.6, the KS algorithm operates in complex arithmetics. To benefit from the vectorization during the whole eigenproblem resolution, the real-equivalent formulation needs to be conserved within each step of the KS algorithm, and not only for the inversions of linear systems. To this end, the original eigenproblems are reformulated and the new problematic is to be able to link the eigenpairs of the K-formulated operators to those of the original ones.

C =   
C K =           d 11 -
          -→ BSR AA(C K ) = c K 11 c K 13 c K 22 c K 32 c K 33 JA(C K ) = 1 3 2 2 3 IA(C K ) =

III.4.2 Eigenvalues of the K-formulated operator

Considering the eigenvalues λ 1 , . . . , λ N of C in any prescribed order, the goal is to explicit the spectrum Sp(C K ). Provided that C is a complex squared matrix, there exists a so-called complex Schur decomposition verifying:

C = X SX H , (III.67)
where X ∈ C N × C N is a unitary matrix10 , and S ∈ C N × C N is an upper triangular matrix whose diagonal entries are the eigenvalues of C:

S =     λ 1 ⋆ . . . 0 λ N     .
(III.68)

From property (III.64), the following relation holds:

C K = X K S K X H K , (III.69)
where S K reads:

S K =     Λ 1 ⋆ . . . 0 Λ N     , with Λ q = ℜ(λ q ) -ℑ(λ q ) ℑ(λ q ) ℜ(λ q )
, for all q ∈ 1, N .

(III.70)

Proposition 1 Sp(C K ) = Sp(C) ∪ Sp(C).
Proof Let us first focus on the attributes of X K . From properties (III.63) and (III.64), we have:

X K X H K = (X X H ) K = I K N = I 2N . (III.71)
Moreover, if the general term of X is (x qr ) 1≤q,r≤N , then that of X H writes (x rq ) 1≤r,q≤N and the general terms of X K and X H K are, for a given (q, r) ∈ 1, N × 1, N :

x K qr = ℜ(x qr ) -ℑ(x qr ) ℑ(x qr ) ℜ(x qr ) and x K rq = ℜ(x qr ) ℑ(x qr ) -ℑ(x qr ) ℜ(x qr ) . (III.72)
Equation (III.72) highlights the relation X H K = X K T , which combined with (III.71), yields:

X K X K T = I 2N . (III.73)
We can therefore state that the K formulation X K of the unitary matrix X is an orthogonal matrix, and the decomposition (III.69) can be rewritten under the following form:

C K = X K S K X K T . (III.74)
Considering the shape of S K , we can stated, from [START_REF] Aldhahcri | A real schur form method for modeling singularly perturbed systems[END_REF], that the relation (III.74) is a real Schur decomposition. In this real form, the eigenvalues are contained in 2-by-2 diagonal blocks (or scalar for real eigenvalues) instead of complex scalars. The consequence of this result is that Sp(C K ) = Sp(S K ) and the eigenvalues of S K are those of all of its 2-by-2 diagonal blocks, i.e.:

Sp(S K ) = N q=1 Sp(Λ q ) = N q=1
{λ q , λ q }, (III.75) which demonstrates proposition 1:

Sp(C K ) = Sp(C) ∪ Sp(C) . (III.76)
The operator C K thus contains the spectral informations of both C and C. The next objective is to recover the eigenvectors of C or C from those of C K . This point is addressed in the following.

III.4.3 Eigenvectors of the K-formulated operator

As stated in [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF], C K and C K1 are orthogonally similar matrices linked through the relation:

C K = PC K1 P T , (III. 77
)
where P is a permutation matrix that, for a given (q, r) ∈ 1, 2N × 1, 2N , is of general term:

p qr =       
1 if q is odd and r = (q + 1)/2, 1 if q is even and r = q/2 + N , 0 otherwise. In the case N = 3, this leads to:

Pυ =          
1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 Proof This proposition is an equivalence that is demonstrated showing both implications. We first focus on the left-to-right implication. Considering that υ o and υ e are eigenvectors of C associated with the eigenvalue λ, we have Cυ o = λυ o and Cυ e = λυ e . Expanding these two relations leads to the following set of equations: Pυ being an eigenvector of C K1 , equation (III.79) can be recast as:

                    υ 1 υ 2 υ 3 υ 4 υ 5 υ 6           =           υ 1 υ 3 υ 5 υ 2 υ 4 υ 6           . ( III
             Dυ o ℜ -Eυ o ℑ = αυ o ℜ -βυ o ℑ Eυ o ℜ + Dυ o ℑ = βυ o ℜ + αυ o ℑ Dυ e ℜ -Eυ e ℑ =
D -E -E D υ o ℜ + iυ o ℑ υ e ℜ + iυ e ℑ = (α + iβ) υ o ℜ + iυ o ℑ υ e ℜ + iυ e ℑ .
(III.86

)
This system is equivalent to the following new set of equations: Proof The proof is similar to that of proposition 2, except that the system (III.82) to (III.85) is replaced by: Proof λ being real, the system (III.87) to (III.90) simplifies as: For a given eigenpair of C K , whether the eigenvalue is complex or real, we are hence able to recover an eigenpair of C. We now have all the required informations to adapt the KS algorithm to K formulation. The next section provides the implementation details of the K formulation applied to both GMRES and KS algorithms.

             Dυ o ℜ -Eυ e ℜ = αυ o ℜ -βυ o ℑ Eυ o ℜ + Dυ e ℜ = αυ e ℜ -βυ e ℑ Dυ o ℑ -Eυ e ℑ = βυ o ℜ + αυ o ℑ Eυ o ℑ + Dυ e ℑ =
βυ o ℜ + αυ o ℑ = αυ o ℑ -βυ e ℑ . ( III
             Dυ o ℜ + Eυ o ℑ = αυ o ℜ -βυ o ℑ -Eυ o ℜ + Dυ o ℑ = βυ o ℜ + αυ o ℑ Dυ e ℜ + Eυ e ℑ =
             Dυ o ℜ -Eυ e ℜ = αυ o ℜ Eυ o ℜ + Dυ e ℜ = αυ e ℜ Dυ o ℑ -Eυ e ℑ = αυ o ℑ Eυ o ℑ + Dυ e ℑ =

III.5 K-formulated global stability analysis method

To get a complete working implementation of global stability analysis in K formulation, each component of the eigenproblem resolution must be adapted. This section first presents the adaptation of the Jacobian-vector products of Tapenade, which is one of the main elementary component of the GMRES. The preconditioning tools are also adapted and some convergence results of the GMRES in K formulation are provided for the laminar cylinder case.

III.5.1 K-formulated matrix-vector product

Producing differentiated code with Tapenade in K formulation is not only doubling the size of each input vectors. In fact, as mentioned earlier, the contiguity of the real and imaginary components of the input K-formulated vector must be conserved. To this end, the input vector must be split and Jacobian-vector products must be managed on the real and imaginary parts separately. The description of the process is described below.

For the sake of illustration and without loss of generalities, we consider here that B = J and δ = 1. The product between C = J + sI N and a given complex vector υ ∈ C N writes:

Cυ = [J + (s ℜ + is ℑ ) I N ] (υ ℜ + iυ ℑ ) = J υ ℜ + s ℜ υ ℜ -s ℑ υ ℑ + i (J υ ℑ + s ℜ υ ℑ + s ℑ υ ℜ ) =     (J υ ℜ ) 1 + s ℜ υ ℜ 1 -s ℑ υ ℑ 1 . . . (J υ ℜ ) N + s ℜ υ ℜ N -s ℑ υ ℑ N     =ℜ(Cυ) +i     (J υ ℑ ) 1 + s ℜ υ ℑ 1 -s ℑ υ ℜ 1 . . . (J υ ℑ ) N + s ℜ υ ℑ N -s ℑ υ ℜ N     =ℑ(Cυ)
.

(III.102)

The general term of the K formulation C K = (c K qr ) 1≤q,r≤N of C is expressed as follows:

c K qr =            J qr + s ℜ -s ℑ s ℑ J qr + s ℜ if q = r,
J qr 0 0 J qr otherwise.

(III.103)

Now, reformulating υ in υ K = υ ℜ 1 υ ℑ 1 • • • • • • υ ℜ N υ ℑ N T ∈ R 2N
, the K formulation of Cυ reads:

C K υ K =             (J υ ℜ ) 1 + s ℜ υ ℜ 1 -s ℑ υ ℑ 1 (J υ ℑ ) 1 + s ℜ υ ℑ 1 -s ℑ υ ℜ 1 . . . . . . (J υ ℜ ) N + s ℜ υ ℜ N -s ℑ υ ℑ N (J υ ℑ ) N + s ℜ υ ℑ N -s ℑ υ ℜ N             . (III.104)
To achieve the computation of C K υ K , the vector υ K has first to be split into υ ℜ and υ ℑ . Given the Jacobian-vector products provided by Tapenade (see II.4.1), the two products J υ ℜ and J υ ℑ must then be performed. They are finally assembled in a 2N -dimensional vector with the remaining terms of sυ. Using the notation of II.4.2, the graph providing the K-formulated right matrix-vector product, i.e. the mapping φ K ∈ R 2N → J K φ K ∈ R 2N , reads:

J K φ K φ K ∂P ∂W φ ℜ ∂F ∂P ∂P ∂W φ ℜ + ∂F ∂∇P ∂∇P ∂P ∂P ∂W φ ℜ ∂∇P ∂P ∂P ∂W φ ℜ ∂P ∂W φ ℑ ∂F ∂P ∂P ∂W φ ℑ + ∂F ∂∇P ∂∇P ∂P ∂P ∂W φ ℑ ∂∇P ∂P ∂P ∂W φ ℑ
while that of the left matrix-vector product, i.e. the mapping

ψ K ∈ R 2N → J T K ψ K ∈ R 2N , reads: ψ K T J K ψ K ψ T ℜ ∂R ∂F ψ T ℜ ∂R ∂F ∂F ∂P + ψ T ℜ ∂R ∂F ∂F ∂∇P ∂∇P ∂P ψ T ℜ ∂R ∂F ∂F ∂∇P ψ T ℑ ∂R ∂F ψ T ℑ ∂R ∂F ∂F ∂P + ψ T ℑ ∂R ∂F ∂F ∂∇P ∂∇P ∂P ψ T ℑ ∂R ∂F
∂F ∂∇P

A possible implementation consists in extracting and storing the real and imaginary parts from the input K-formulated vector, successively computing the two Jacobian-vector products, and recombining them into a new K-formulated vector. This strategy requires three copies for each matrix-vector product in K formulation, resulting in a large number of additional copies when thousands of matrix-vector products are needed as in global stability computation. Alternatively, by inlining the product in K formulation, the real and imaginary Jacobian-vector products might be computed simultaneously on a 2N -dimensional vector conveniently arranged. This strategy avoids the three copies and improves the vectorization and cache-blocking optimization. It has been the object of a collaboration with Hascoët and resulted in a new feature of Tapenade in which the operations on the real and imaginary parts are done directly on the 2N -dimensional vector.

Examples of K-formulated differentiated codes are provided in appendix B.

III.5.2 GMRES algorithm in K formulation

By doubling the size of the Krylov vectors and using the matrix-vector product in K formulation, the whole GMRES implementation described in III.3.4 can be used on the K-formulated problem C K z K = t K . Likewise, the preconditioning routines already implemented in elsA can be reused on the K-formulated operators, doubling the block size, as shown in III.4.1. The jacobian matrix is formed at the beginning of the computation, and the K-formulated operator is obtained by transforming AA(J ) to AA(J K ).

The resulting implementation of the GMRES in K formulation has been tested on the Kformulated linear systems associated to the laminar cylinder case. As in III.3.4, both external and internal Krylov basis are composed of 60 vectors, and the solver is the FGMRES-ILU(3). Results are presented in fig. III.10, and table III.4 shows the allocated memory and time of each computation. We first observe that the resolution converges faster with the K-formulated adjoint operator than for the direct one. This agrees with the results of Timme [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF] who also noticed such behaviour. On the other hand, when s = 0, convergence is more difficult than when s ̸ = 0. This is explained by the reinforcement of the diagonal dominance (see III.3.4) with a nonzero shift. Additionally, when the shift is really close to an eigenvalue, the operator is almost singular resulting in a large condition number and a more difficult inversion. A balance is thus to be found between setting a shift with a large real part, which makes the inversions easy and fast but requires much more eigensolver iterations to converge toward the desired eigenvalue; and setting a shift

ϵ rel J J K (J -sI N ) K J T K J T -sI N K Figure III
.10: GMRES convergence curves for various operators on the laminar cylinder case at Re = 47. The dashed line is at 10 -5 which is the selected tolerance for global stability analysis about cylinder flow. the shift is s = 1 + 0.9i. very close to the selected eigenvalue, which makes the eigensolver faster but the systems stiffer ant the inversions slower.

Operator J J K (J -sI N ) K J T K J T -

III.5.3 Adaptating the Krylov-Schur algorithm to the K formulation

The KS algorithm also have to be adapted to K formulation. It has been implemented in elsA in its real version, in which the complex Schur decomposition in equation (III.67) is replaced by a real one [START_REF] Aldhahcri | A real schur form method for modeling singularly perturbed systems[END_REF], where the eigenvalues are contained in 2-by-2 diagonal blocks (or scalar for real eigenvalues). The size of the Krylov vectors that compose the KS basis has been doubled to comply with the K formulation. For the shift-and-invert transformation, the convergence test has been modified to test the errors associated with eigenvectors of J instead of J K . Some details of these implementations are provided below.

Modal analysis. Whether we focus on the direct of adjoint eigenproblem, the operator of interest has the form of a shift-and-invert operator, as shown in II.5. The KS algorithm builds a KS basis from K formulation of operator C -1 = (B -sI N ) -1 . The homomorphism properties (III.63) and (III.64) of the K formulation provide:

C -1 K = C K -1 .
(III.105)

The KS basis is generated by successive evaluations of C K -1 on the Krylov vectors through GMRES algorithm. This basis is successively extended and truncated to find the eigenvalues λ K l of C K -1 , which are linked to the eigenvalues λ K l of C K through the relation λ

K l = λ K -1 l
. Both C K -1 and C K have the same eigenvectors. The eigenvalues λ l of B can then directly be recovered from the λ K l by using the K formulation properties and the fact that C -1 = C -1 :

λ l =    λ K l + s if λ K l ∈ Sp(C), λ K l + s otherwise. (III.106)
To discriminate one of the two possibles values, a test must be carried out. By using the results of III.4.3 and provided the eigenvector υ l ∈ C 2N associated to λ K l , a test that compares the quantities ∥υ o l + iυ e l ∥ 2 and ∥υ o l -iυ e l ∥ 2 may be implemented. However, based on the author's experience, the relations υ o l = -iυ e l and υ o l = iυ e l do not necessarily hold when the eigenvector is not perfectly converged. We choose instead to focus on the residual associated to an eigenvector of B extracted from υ l . To do so, let us consider υ o l the eigenvector to test and the quantities

ϵ 1 = ∥Bυ o l -(λ + s)υ o l ∥ 2 and ϵ 2 = ∥Bυ o l -(λ + s)υ o l ∥ 2 .
The lowest value between ϵ 1 and ϵ 2 indicates the correct value of λ l in equation (III.106) and gives the residual associated with the operator B. Hence, the KS algorithm is entirely implemented in K-formulated shift-and-invert, Algorithm 6 K-formulated KS algorithm for modal analysis m -p Arnoldi steps with deflation using preconditioned GMRES:

Input: B, s ∈ C, v 1 ∈ R 2N
(B -sI N ) K -1 V m = V m H m + v m+1 bT m 4:
real Schur decomposition with sorting

H m ← W T 2 W T 1 H m W 1 W 2 5: V m ← V m W 1 W 2 6:
for k ∈ j, l do

7: compute ( λ K k , y k ) verifying H m y k = λ K k y k and set υ k = V m y k 8: compute ϵ 1 = ∥Bυ o k -( λ K -1 k + s)υ o k ∥ 2 9: compute ϵ 2 = ∥Bυ o k -( λ K -1 k + s)υ o k ∥ 2 10: set ϵ KS rel = min(ϵ 1 , ϵ 2 )/(|λ k |∥υ o k ∥ 2 ) 11: if ϵ KS rel < tol then 12:
lock (λ k , y k ) and continue compute b p and insert it into H p 20: end while but the convergence test is done on the original operator instead of the K-formulated one. The implementation is given in algorithm 6. It is worth noting that such a test is associated to an eigenpair of B instead of C -1 , which might give too optimistic bounds and provide an eigenvalue of B that is not really converged [START_REF] Vicente | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF].

Resolvent analysis. The shift-and-invert transformation is not used here, but the operator of interest still involves inverted operators. As stated in II.5.3, the resolvent analysis can be achieved by solving the eigenproblem associated to one of the operator

RQ -1 2 R H Q 1 or Q -1 2 R H Q 1 R. When solving the eigenproblem associated with Q -1 2 R H Q 1 R
, the resulting eigenvectors are the optimal forcings and, as seen with relation (II.57), each optimal response's computation involve an inversion. In this work, the focus is made on the operator C = RQ -1 2 R H Q 1 , whose eigenvectors provide the optimal responses, from which the optimal forcings can be recovered by using a simple matrixvector product with the operator (iωI N -J ). Based on what we know about the K formulation, the expression of C K reads:

C K = R K Q K 2 -1 R H K Q K 1 = (-J -iωI) K -1 Q K 2 -1 (-J T + iωI) K -1 Q K 1 . (III.107)
From [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF], this operator is also Hermitian and has the same eigenvalues as C, each with doubled multiplicity 11 . Moreover, we know that a real matrix with real eigenvalues have complex eigenvectors. However, we ensure the eigenvectors are real by operating with real arithmetics, Algorithm 7 K-formulated KS algorithm for resolvent analysis

Input: J , J T , Q 1 , Q 2 , ω ∈ R, v 1 ∈ R 2N
unitary, m and l ≤ m the number of desired eigenpairs Output:

j eigenpairs of RQ -1 2 R H Q 1 , with j ≤ l 1: initialization of V m ← v 1 0 , j ← 0, p ← 0 2: while k ≤ l do 3:
m -p Arnoldi steps with deflation using preconditioned GMRES:

(-J -iωI) K -1 Q K 2 -1 (-J T + iωI) K -1 Q K 1 V m = V m H m + v m+1 bT m 4:
real Schur decomposition with sorting

H m ← W T 2 W T 1 H m W 1 W 2 5: V m ← V m W 1 W 2 6:
for k ∈ j, l do

7: compute ( η k , y k ) verifying H m y k = η k y k and set ϵ KS rel = | bT m y k |/(| η k |∥ V m y k ∥ 2 ) 8: if ϵ KS rel < tol then 9:
lock ( η k , y k ) and continue compute b p and insert it into H p 17: end while meaning if the eigenvalue is real, the corresponding eigenvector will also be. From proposition 4, by considering an eigenpair (υ l , ηl ) ∈ R 2N × R of C K , the vector υ o l + iυ e l is selected as an eigenvector of RQ -1 2 R H Q 1 and the resulting optimal forcing is calculated by

√ η l -1 R -1 (υ o l + υ e l ).
No spectral transformation is performed in this resolvent case. The error calculated from the Krylov decomposition (see equation (II.73)) directly provides the residual associated to the K-formulated eigenproblem. The calculation of the residual associated to C is avoided since it requires two inversions for each eigenpair to test. An error is thus committed because the convergence is based on the K-formulated problem. However, based on the author's experience, the error of an eigenpair of C K is usually of the same order of magnitude as the error associated to the extracted eigenpair of C. An implementation of the K-formulated KS algorithm for resolvent analysis is presented in algorithm 7.

These adaptations of the KS algorithm to K formulation complete the description of the method to scale up global stability analysis to industrial configurations, denoted as K-Formulated Global Stability Method (KFGSM). In the next section, the KFGSM is tested to draw an exhaustive picture of the stability of the laminar cylinder flow, and the results are compared to the reference SLEPc eigensolver.

III.6 Global stability analysis of the laminar cylinder flow

The numerical experiments carried out on the laminar cylinder case to validate the KFGSM are described hereafter. The various eigenproblems that arise in the global stability analysis are solved for the base flows discussed in III.3.5, and the results are compared to those obtained from the reference eigensolver SLEPc [START_REF] Vicente | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. The latter uses complex arithmetic, forms the operators explicitly through multicoloring (see III.3.4), and solves the complex linear systems through an LU factorization, while the KFGSM uses real arithmetic and the K formulation, and employs FGMRES-ILU(3) with an external tolerance of 10 -6 . The relative tolerance for the eigenpairs for all calculations has been set to 10 -4 . The KS basis includes 20 vectors for direct and adjoint calculations and 10 vectors for resolvent analysis. We also introduce the direct and adjoint perturbation velocities U and U † .

III.6.1 Direct modal analysis

The eigenpairs of J were calculated by solving the eigenproblem in equation (II.45) using both SLEPc and elsA. The relative tolerance for the two solvers was set to ϵ KS rel = 10 -4 . The amplification rate σ, and Strouhal number St = 2π/ω, obtained for the least stable eigenvalues are shown as a function of Re, in fig. III.11a and III.11b, respectively. These values agree with those obtained by Gianetti and Luchini in [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. The critical value of Re at which the global mode becomes unstable is Re c = 47, which is consistent with values found in the literature [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF][START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. These results validate our numerical setup. The spatial structures of the streamwise and crossstream perturbation velocities, normalized by the L ∞ -norm of the perturbation velocity vector, are shown in fig. III.12. Coherent structures, with an anti-symmetry through x-z plane, can be seen to be temporally and spatially amplified with a spatial wavelength of Λ x ≃ 6.3L, and advected downstream. The von Kármán vortex street appears when the Reynolds number is increased beyond Re c , due to this unstable mode, which denotes a linearly amplified mechanism. The mode having a nonzero frequency, it is oscillating and thus denotes a so-called Hopf bifurcation 12 . These linear computations are confronted to nonlinear unsteady ones in III.7.

Moreover, a good match is observed between the eigenvalues computed with elsA and SLEPc. This result validates the implementation of the K-formulated KS algorithm in elsA using real arithmetic, which returns the same eigenvalues as the complex reference eigensolver SLEPc. Table III.5 also shows that elsA returns the same eigenvalue using 1 or 8 CPUs, demonstrating its ability to operate in a massively parallel environment. However, in this particular case, the SLEPc solver is much faster than elsA. This was expected as our solver, designed for industrial configurations, has a much higher arithmetic intensity than the SLEPc counterpart used here, which relies on LU factorizations. This SLEPc strategy may be used for 2-D cases and small 3-D ones, but is unable to handle 3-D configurations with a high number of DOF, as discussed in III.

Solver

SLEPc -

III.6.2 Adjoint modal analysis

The same numerical experiments as in the direct eigenproblem are shown for the adjoint one (see (II.49)). The mass matrix considered is the 5N -by-5N matrix Q Ω (see II.5.2) of general term:

q Ω =         V (Ω) 0 0 0 0 0 V (Ω) 0 0 0 0 0 V (Ω) 0 0 0 0 0 V (Ω) 0 0 0 0 0 V (Ω)         . (III.108)
The least stable eigenpairs was computed with both elsA and SLEPc eigensolvers. The evolutions of σ and St as a function of Re are presented in fig. III.13a and III.13b, respectively. A good agreement between direct and adjoint eigenvalues is observed, demonstrating the duality between of the direct and adjoint eigensolvers. The critical Reynolds number is also found to be Re c = 47. Fig. III.14 displays the spatial structure for the streamwise and cross-stream perturbation velocities, normalized by the L ∞ -norm for the least stable mode at Re c . Both shapes are consistent with previous results from the literature [START_REF] Marquet | Stabilité globale et contrôle d'écoulements de recirculation[END_REF]. Unlike the direct mode, the adjoint is advected upstream. As stated by Sipp [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF], this is due to the opposite signs of the advection operator in the adjoint and direct perturbation equation. The areas with the brightest colors correspond to the zones where the flow instability is the most sensitive to an external forcing. This information is useful to develop flow control strategy in order to delay the instability onset [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF].

As for the direct eigenproblem, a good match is observed between elsA and SLEPc, which validates the K-formulated KS algorithm for the adjoint operator. Table III.6 presents detailed results after various computations. The eigenvalue returned by elsA using 1 or 8 CPUs is the same, further confirming the solver's ability to operate in a massively parallel environment. Additionally, it can be noted that the adjoint eigensolver implemented in elsA converges faster than the direct one, in agreement with Timme [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF]. Again, SLEPc is faster than elsA for the reasons mentioned above.

Given both direct and adjoint unstable modes, and by analogy with equation (I.7), the structural sensitivity can be computed in its discrete form to describe the critical zones in the domain responsible for the instability (the so-called wavemaker). From [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF], it is defined as:

S(x) = ∥ U (x)∥ 2 ∥ U † (x)∥ 2 , (III.109)
where the L 2 norm is complex. The resulting map is shown in fig. III.15 and agrees with the results from [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. It displays two dark parts around x = 2.3, indicating the precise origin of the instability. This means that modifying the base flow outside these zones will have a marginal effect on the instability mode as long as these sensitive zones remain unaltered. In an industrial context, this information is valuable as it also provides insights into how to distribute meshing points economically while still capturing the instability.

III.6.3 Resolvent analysis

The KFGSM was ultimately tested on the resolvent analysis for the subcritical Reynolds number Re = 40. The energetic gain was optimized with respect to Chu's energy [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF] for which the mass matrix Q Chu is of general term: [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF] shows the two optimal gains obtained at this regime. A peak at St = 0.117 is observed, corresponding to the frequency of the least stable mode at this regime. This most amplified frequency corresponds to that of the vortex shedding observed when the von Kàrmàn instability is triggered. Spatial structures of streamwise and cross-stream velocities for Φ opt and W opt are presented in fig. III.17. As expected, the optimal forcing has the same shape as the adjoint mode of fig. III.14, while the optimal response is similar to the direct mode of fig. III.12.

q Chu = V (ω) ρ 2         ∥u∥ 2 2 -u x -u y -u z 0 -u x 1 0 0 0 -u y 0 1 0 0 -u z 0 0 1 0 0 0 0 0 0         . (III.110) Fig. III.
The eigenvalues calculated with elsA match well with those obtained with SLEPc. Table III.7 shows the results of various numerical experiments at St = 0.1165. Although SLEPc is faster than elsA, it requires more memory due to the need to compute and store LU factorizations of R and R H , while elsA only stores their ILU(3) factorizations, which are less expensive. In addition, SLEPc must store the mass matrix, whereas elsA never forms the matrix explicitly, using only matrix-vector products to calculate the actions of Q -1 Ω and Q Chu on a given vector. The results of the numerical experiments show that elsA can work in a massively parallel environment, with the same eigenvalues being obtained using 1 or 8 CPUs.

To ensure the accuracy of the eigenpairs, it is important to set the external tolerance of GM-RES as low as possible, with a maximal value of 10 -6 . This is to ensure the validity of the product (RQ -1 Ω R H Q Chu ) K with a given vector. The KS algorithm is efficient in finding the largest eigenvalues, which is why a large KS basis and multiple iterations may not be necessary in the resolvent case. Indeed, the largest eigenvalue is generally much higher than the other ones and, based on the author's experience, the algorithm usually converges it in one or two iterations, using a relatively small KS basis. This contrasts with direct and adjoint calculations, which use shift-and-invert without knowledge a priori of the eigenvalue, and may require a larger basis and more iterations 0.04 0.08 0.12 0.16 0.2 St

Solver

SLEPc -1 CPU elsA - This section aimed at showing the validation of the KFGSM on both the reference eigensolver SLEPc, and the literature, for modal and nonmodal analyses. A common practice to compare the amplification mechanism in both linear and nonlinear frameworks is to perform nonlinear unsteady computations. In the next section, the second-order Gear method is described and unsteady computations are performed on the laminar cylinder case.

III.7 Unsteady computation

The nonlinear unsteady computation can be performed prior to a global stability analysis in order to determine the frequency of the unsteadiness. This information can then be utilized to choose the shift in a modal analysis or to compute the optimal forcing and response at this frequency through the resolvent analysis. The results of the global stability analysis can also be compared retrospectively with the unsteady computation to verify if the origin of a particular unsteadiness can be attributed to a linearly amplified mechanism. To carry out the temporal integration of the governing equations, the elsA software employs the Gear method, which is explained in more detail below. Unsteady computations have been performed on the laminar cylinder flow and validate the capability of the KFGSM to predict the onset of the von Kármán vortex street.

III.7.1 Gear method

Each unsteady computations in this work are carried out by solving the Gear scheme. The latter, described in [START_REF] Cinnella | High-order implicit residual smoothing time scheme for direct and large eddy simulations of compressible flows[END_REF], is a backward multistep method that is second-order accurate in time. It is associated with the following nonlinear implicit scheme:

G(W n+1 ) = V (Ω) ∆t 3 2 W n+1 -2W n + 1 2 W n-1 -R(W n+1 ) = 0. (III.111)
Similarly to II.3.1, this scheme is equivalent to a Newton method:

J G (W n,l )∆W n,l = -G(W n,l ) (III.112)
where l is the index of the current sub-iteration of the Newton algorithm. As for the LU-SGS method described in III.3.3, the Jacobian matrix is replaced by an approximated operator J and the following system is solved at each nonlinear iteration:

3 2∆t I N + J (W n,l ) ∆W n,l = R(W n,l ) (III.113)
Note that J is updated at each sub-iteration of the Newton resolution. The time-step ∆t corresponds to the physical user-defined one.

III.7.2 Nonlinear unsteady computation of the laminar cylinder flow

The physical time-step for these computations has been set to ∆t = 10 -3 , so that the CFL number at the wall is approximately 10. The initial fields are the base flows presented in III.3.5, and perturbed by least stable mode computed using the direct eigensolver. [START_REF] Drazin | Hydrodynamic stability[END_REF], which shows isocontours of vorticity magnitude at t = 300. A spatial wavelength of Λ x ≃ 6.5L is also observed, which is of the same order as that identified for the direct unstable mode (see Fig. III.12).

The measured signal is analyzed in fig. III.20a and III.20b to extract the corresponding nonlinear amplification rate σ NL and Strouhal number St NL . Only one peak is identified in the magnitude spectrum, which is associated to the frequency of the vortex shedding. Table III.8 summarizes the results obtained for various Reynolds numbers, where (σ L , St L ) are the values resulting from the linear direct eigensolver. The growth rates from the linear analyses are consistent with those obtained from the linear stage of the nonlinear simulations. Differences in both the amplification rate and Strouhal number can be seen between the linear and nonlinear computations. As stated for instance in [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF], while this discrepancy is normal since it is due to nonlinearities, the actual nonlinear frequency may nonetheless be recovered from the linear stability analysis by linearizing the equations about the mean flow instead of the fixed point. This is further discussed for instance in [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. However, at relatively low Reynolds number, these results support the conclusion that the vortex shedding originates from a linearly amplified mechanism.

The unsteady computation is hence useful in identifying unsteadiness in the flow, but when multiple unsteady phenomena coexist, they interact nonlinearly and it is often difficult to isolate them. Additionally, the required computation time to reach the limit cycle can be quite long. This is why global stability analysis is useful, as it allows for the isolation of different mechanisms in a reasonable amount of computation time.

Re

III.8 Conclusion

This chapter first presents the fixed point methods, including preconditioning computation, as implemented in elsA. Then, a mathematical foundation for the KFGSM is established. The code is validated through testing on the laminar cylinder case. In chapter III, the KFGSM is applied to more complex configurations, including the turbulent OAT15A airfoil flow in transonic conditions, which is a challenging RANS case, and the 3-D laminar subsonic sphere flow, which contains a large number of DOF. Bifurcations of axisymmetric and planar-symmetric base flows are successfully characterized.

Chapter IV

Challenging test cases

IV.1 Introduction

Chapter II aims to introduce general concepts related to compressible CFD and global stability. In chapter III, the KFGSM is described and validated on the canonical case of the laminar cylinder flow to predict the emergence of the von Kármán vortex street. In this chapter, the capability of the KFGSM to handle more challenging cases is demonstrated. The first case studied is the RANS OAT15A airfoil flow, which destabilizes through a Hopf bifurcation resulting in the buffet phenomenon, an interaction between the boundary layer and a shock wave located on the suction side of the airfoil [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF]. The linear systems involved in both fixed point computation and global stability analysis are known to be particularly stiff [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF].

The second case considered is the 3-D subsonic laminar sphere flow. This flow has two bifurcations: a regular one around Re = 210 that breaks one symmetry, and a second one around Re = 280 which is a Hopf bifurcation, responsible for the degeneration towards hairpin-like vortex shedding [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]. The computational grid for this case consists of about 10 × 10 6 cells, resulting in 50 × 10 6 DOF. For the super-critical Hopf bifurcation of the planar-symmetric base flow, a refined mesh consisting of about 125 × 10 6 DOF was used. Such sizes are representative of those of industrial configurations that the code was designed for. This demonstrates the capability of the KFGSM to operate on such configurations where direct factorizations are not feasible due to memory limits. As a reminder, the largest known global stability analysis conducted by Timme to characterize the buffet phenomenon on an entire plane was performed on 50 × 10 6 DOF [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF].

IV.2 OAT15A airfoil case

Under given Mach number and angle of attack, the shock on the suction side of a plane wing may interact with the detached boundary layer and gives rise to the so-called unsteady buffet phenomenon [START_REF] Brunet | A complete Experimental and Numerical Study of the Buffet Phenomenon over the OAT15A Airfoil[END_REF]. This causes the shock to oscillate through a self-sustained mechanism, and may lead to severe damage, mainly due to fatigue. Crouch first showed in [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF] that this mechanism could be predicted on a 2-D airfoil through the use of global stability analysis. Sartor conducted such an analysis on the OAT15A airfoil as described in [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], and Paladini [START_REF] Paladini | Transonic buffet instability: From two-dimensional airfoils to three-dimensional swept wings[END_REF] extended the study to 3-D unswept wing, exhibiting buffet cells on the boundary layer. In this section, we focus on the 2-D supercritical OAT15A airfoil to validate the KFGSM on this particularly stiff RANS case [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. As for the laminar cylinder flow, fixed point computation as well as direct, adjoint, and resolvent analysis are detailed below.

IV.2.1 Numerical set-up

The computational domain presented in fig. IV.1 is an unstructured C-grid composed of 70091 hexaedral cells. The first mesh point in the boundary layer on the profile is always below y + = 0.9, and considering the chord length c as the characteristic dimension, the grid definition in the shock region is ∆x/c = 0.003. The selected physical parameters are given in table IV.1. To improve the resolution of linear systems during fixed point computations or global stability analysis, the problem is undimensioned using stagnation density ρ i,∞ , stagnation temperature T i,∞ , and the associated sound-velocity a i,∞ , which means that these quantities are all unitary. Such a modification of the original problem is a form of preconditioning that improves the condition number of the linear systems. Other physical quantities are obtained with the following relations:

a ∞ = 1 + γ -1 2 M 2 ∞ -1 , (IV.1) p ∞ = 1 + γ -1 2 M 2 ∞ -γ γ-1 , (IV.2) ρ ∞ = 1 + γ -1 2 M 2 ∞ -1 γ-1 , (IV.3) e ∞ = 1 γ(γ -1) 1 + γ -1 2 M 2 ∞ -1 , (IV.4) R = 1 γ , (IV.5) µ ∞ = cρ ∞ M ∞ a ∞ Re ∞ . (IV.6)
The angle of attack α is taken in to account by imposing at the farfield boundaries the streamwise u x,∞ = u ∞ cos(α) and cross-stream u y,∞ = u ∞ sin(α) velocities.

Concerning the numerical methods, the negative SA scheme presented in III.2.2 is used to model turbulent effects. Convective fluxes are discretized using the approximate Riemann solver HLLC scheme as described in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF], with second-order accuracy without a limiter function. The diffusive fluxes, as in the case of the laminar cylinder case, are discretized by computing the gradient at faces using the 5p cor method (see III.3.2). As for laminar cylinder flow (see III.1.2), the boundary conditions are set to no-slip velocity at the wall, non-reflective conditions in the farfield, and symmetry conditions on the transverse faces of the domain.

IV.2.2 Fixed point computation

Under unstable regimes, computing a fixed point in this case is particularly difficult. In fact, using the backward Euler scheme with GMRES, as presented in III.3.4, initialized with a uniform field fails and causes the calculation to diverge. We chose instead to first perform a LU-SGS computation (see III.3.3) initialized with a uniform field to get a non-converged field that can then be used to initialize the backward Euler scheme combined with the GMRES. Both simulations were run on 24 CPUs, and the convergence result for the LU-SGS method is presented in fig. IV.2a. The residual ∥ρ∥ 2 was reduced to 7 × 10 -7 in 10 4 iterations.

For GMRES computations, the FGMRES-HLU method is employed due to the stiffness of the linear systems of interest (see III. IV.2 provides quantitative data to compare the two methods. As expected, the memory requirements for the GMRES method are significantly higher than those of the LU-SGS method, as the GMRES method solves the exact linear systems rather than approximate ones, requiring the formation of an effective preconditioner, which consumes a lot of memory.

Such results demonstrate the complementary nature of both fixed point computation methods and the efficacy of using them successively to compute a fixed point in a stiff case without prior knowledge of its shape. However, in an industrial setting, particularly for turbomachinery computations, the LU-SGS method may not converge at all when initialized with a uniform field. This issue can be addressed by approximating the initial field using the meridian model, as described in [START_REF] Godard | Compresseurs et turbines[END_REF]. Based on these successful base flow computations, the focus is made in the following on the global stability analyses of this case.

IV.2.3 Direct modal analysis

With the fixed points computed with high accuracy, the KFGSM can now be tested on this configuration. The computations are still performed on 24 CPUs. The KS basis consists of 30 vectors and the relative tolerance for the eigenvalue is set to 10 -4 . The inversion method used is the FGMRES-HLU with 2-level RAS [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF], where both the internal and external Krylov bases consist of 60 vectors. The internal and external tolerances are set to 10 -2 and 10 -6 , respectively. For this case, the linear systems are so stiff that the real part of the shift was set to 5 to make the inversions operational. The imaginary part of the shift was set near the expected non-dimensional buffet frequency, which is related to the Strouhal number defined as follows:

St = f c u ∞ . (IV.7)
The evolution of both σ and St with respect to α is presented in fig. IV.4a and IV.4b, respectively. The buffet onset is found to occur at α c = 4 • with the corresponding critical Strouhal number St c = 0.062, resulting in a buffet onset frequency around 64 Hz. The flow restabilizes between α = 5.5 • and α = 6 • , which indicates the buffet offset. These results are in agreement with the results obtained by Sartor in [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF][START_REF] Morisco | Dynamique nonlinéaire et stabilité linéaire d'une tuyère sur-détendue[END_REF]. The contour maps of both velocity and density of the direct unstable mode at α c are shown in fig. IV.5. The most energetic part of the mode is located within the shock wave, and a smaller part is identified in the mixing layer.

As with the laminar cylinder flow, the elsA solver is compared to the SLEPc solver to validate the implementation of the K-formulated KS algorithm. A good match is observed for both σ and St. Table IV.3 provides quantitative figures to compare both solvers at α c . The elsA computation requires much more memory than the SLEPc solver, which is not surprising since the GMRES inversions require many arrays to be allocated for the Krylov basis and the two-level overlapping preconditioning is memory-intensive. However, it should be noted that this case serves as a demonstration and is not representative of industrial cases for which it was designed. The direct modal analysis is hence validated for this RANS case, and the adjoint modal analysis is described in the following.

IV.2.4 Adjoint modal analysis

The adjoint operator is defined relative to the scalar product defined by the 6N -by-6N matrix Q Ω (see equation (III.108)). Similar to the laminar cylinder flow and in agreement with Timme [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF], the adjoint problem is less stiff than the direct one, allowing the shifts to be set closer to the expected eigenvalues, resulting in fewer KS iterations for the solver to converge. As shown in table IV.4, both solvers only required one KS iteration with a tolerance of 10 -4 . In comparison with the direct eigensolver, the memory requirements for the elsA adjoint solver are also 70% higher than those of SLEPc due to the HLU preconditioning with two RAS levels.

The evolution of both σ and St with respect to the angle of attack is shown in fig.

IV.6 and compared with SLEPc. The close proximity of the eigenvalues confirms the validity of the adjoint eigensolver under RANS conditions. Furthermore, the eigenvalues computed using the (a)

Solver

SLEPc - adjoint solver are identical to those obtained using the direct eigensolver within a certain tolerance, demonstrating the duality between the two operators. The least stable adjoint mode calculated for α c is displayed in fig. IV.7. In accordance with the findings of Sartor [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], the most energetic part of the adjoint is the oblique component which follows the characteristic line that impacts the suction side at the start of the separation region. This suggests a strong receptivity of the flow to external forcing in this area. As noted by Sartor, this position of the adjoint could indicate a significant influence of the start of the separation region on the flow's dynamics. Additionally, the placement of the adjoint indicates that a perturbation downstream of the airfoil could impact the instability, as a result of the propagation of upstream acoustic waves. The structural sensitivity is calculated using equation (III.109), and the contour map for α c is shown in fig. The wavemaker is identified as having a lambda shape and is mainly located in the shock wave. As a result, the grid definition in this mesh zone should be carefully selected to capture the instability.

In agreement with the direct eigensolver, the KFGSM in its adjoint form is validated for RANS computations. This implies that the K-formulated resolvent operator can be assembled and made operational. In the following, a resolvent analysis of a sub-critical base flow is demonstrated, illustrating the flow's receptivity to the various frequencies calculated.

IV.2.5 Resolvent analysis

The mass matrices are set to the 6N -by-6N matrices Q Chu and Q Ω (see equations (III.110) and (III.108) respectively). The angle of attack is set to the sub-critical value of α = 3.75 • . The tolerance for GMRES inversions must be set lower than for the direct or adjoint eigensolver to ensure accuracy when performing matrix-vector products with the operator

R K Q K Ω -1 R H Q K Chu .
The tolerance is set to 10 -6 for all computations. The KS basis consists of 10 vectors, requiring 20 inversions for the first iteration and 10 for subsequent iterations, as only half of the basis is recycled through the KS iterations. This relatively small number of vectors compared to the size of the problem is due to the fact that the largest eigenvalue is typically well separated from the others, leading to rapid convergence.

The two highest gains computed for various frequencies with elsA are shown in fig.

IV.9 and compared with SLEPc results. The results of are close to each other, but accuracy could be improved by carrying out the same analysis with more vectors in the KS basis. Table IV.5 provides a comparison of both solver computations. elsA require to store two HLU factorizations with two RAS level, which results in an extra 80% memory usage compared to SLEPc. The linear systems in resolvent analysis are stiff, and the relatively high tolerance set for GMRES convergence leads to a much higher CPU time than SLEPc, which operates with the faster LU resolution. This is consistent with the results from the direct and adjoint modal analysis.

From a physical perspective, the evolution of η 1 with respect to St shows two peaks. The

Solver

SLEPc - 3). This indicates that the sub-critical base flow at α c can become unstable by amplifying such frequency. A secondary instability is identified at St = 2.44, resulting in a dimensional frequency of 3430 Hz. The corresponding optimal forcing and response are plotted in fig. IV.11. The optimal response is related to a Kelvin-Helmholtz instability. As stated by Sartor in [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], due to the higher frequency, both the optimal response and forcing exhibit smallscale structure patterns. Furthermore, because of the receptivity of the flow at these medium frequencies, the optimal forcing shape could be used to stabilize the unstable buffet phenomenon that occurs at lower frequencies.

We have just demonstrated that the developed KFGSM is fully operational for RANS cases, ranging from modal to non-modal analysis. As for the laminar cylinder case, nonlinear unsteady computations have been performed, and the results are shown below.

IV.2.6 URANS computation

The accuracy of the linear solver in predicting the unsteady buffet phenomenon was tested by conducting nonlinear URANS computations using the Gear method (see III.7). The physical non-dimensional time-step was set at ∆t = 2.2 × 10 -5 to maintain a CFL number at the wall of approximately 10. The calculations were initiated with the base flow perturbed by the unstable direct mode. Fig. IV.12 shows the exponential temporal evolution of the lift coefficient at α = 4.5 • . The limit cycle is reached around t = 50, where self-sustained oscillations are observed. The motion at two consecutive instances is depicted by plotting contours of the density gradient in fig. IV.13. A strong lambda shock is visible on the suction side of the profile, a characteristic of the buffet phenomenon [START_REF] Morisco | Dynamique nonlinéaire et stabilité linéaire d'une tuyère sur-détendue[END_REF]. This instability causes the shock to move from x/c = 0.25 at t = 86 to x/c = 0.16 at t = 88.

The temporally amplified signal is analyzed in fig. IV.14a and IV.14b to extract σ NL and St NL , respectively. Table IV.6 compares the values of σ and St for linear and nonlinear computations for various angles of attack. Small discrepancies are observed for the amplification rates between linear and nonlinear computations, except for α = 5.5 • , but in this case, we are far from the criticality and this may be improved by lowering the tolerance of the eigensolver. The frequencies calculated with the linear solver are in good agreement with those obtained with URANS computations, demonstrating the capability of the linear RANS solver in predicting instability both onset and frequency. As discussed by Sartor [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], this result suggests that the second-order harmonics generated by the nonlinearities remain weak, even in very unstable regimes. This outcome is not The second part of this chapter focuses on the 3-D subsonic laminar sphere flow. This case was chosen to evaluate the code's ability to solve large 3-D problems and to identify multiple bifurcations in the converged base flows. The bifurcations are investigated using both modal and non-modal analysis, and the results are presented subsequently.

IV.3 Subsonic laminar three-dimensional sphere flow

The unsteadiness in the subsonic laminar sphere flow is a significant area of research because of its similarity to some industrial problems, as noted by Tezuka et al. [START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF]. Sansica et al. have studied both subsonic and supersonic regimes in [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF][START_REF] Sansica | Laminar supersonic sphere wake unstable bifurcations[END_REF] and describe the frequencies and eigenfunctions associated with the unstable modes that cause various instabilities. A schematic bifurcation diagram for this flow is shown in fig. IV.15. When the Reynolds number is greater than Re 1 c ≃ 210, the base flow is axisymmetric (AS) and has a pair of steady streamwise vortices immediately after the sphere body. After Re 1 c , a regular steady bifurcation occurs, breaking the axisymmetry and transforming the base flow into a planar-symmetric (PS) steady flow. A second bifurcation, a supercritical Hopf bifurcation, occurs around Re 2 c = 280, leading to unsteadiness that takes the form of a large-scale hairpin-like vortex shedding behind the sphere, similar to the laminar cylinder flow (see III.7.2). In this work, we focus on the subsonic case, and both bifurcations are analyzed using direct and adjoint global stability analyses. The receptivity of both types of base flows at various frequencies is also studied through resolvent analysis. Additionally, we demonstrate that the pseudo-transient continuation technique thoroughly used in this thesis (see II.3.2) is robust enough to compute AS base flows that do not have physical existence in regimes above the Hopf bifurcation of PS base flows, revealing another specific Hopf bifurcation.

IV.3.1 Numerical set-up

The computational domain, presented in fig. IV.16, is based on that used by Sansica in [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]. It is an unstructured grid composed of 9.3 × 10 6 hexahedral cells and constitutes a good test case to assess the capability of the KFGSM to perform global stability analysis on a 3-D configuration of size representative of industrial applications. The physical parameters are presented in table IV.7. Like in the cylinder case (see III.1.1), the conservative equations are non-dimensionalized by T ∞ , u ∞ , and ρ ∞ . A second-order Roe scheme with MUSCL reconstruction and no limiter is used for convective fluxes, while diffusive fluxes are computed using the 5p cor presented in III.3.2. A no-slip velocity condition is imposed at the wall, non-reflective boundary conditions based on A different strategy was adopted to compute PS base flows beyond the first regular bifurcation. In fact, initializing the computations with a uniform field always result in an AS base flow, even at Reynolds numbers above the supercritical Hopf bifurcation. To overcome this issue, the uniform field was advanced toward the asymptotic PS base flow using the LU-SGS computation (see III. IV.8 provides comparative figures between both resolution methods. As with the OAT15A case (see IV.2.2), the memory requirement of the FGMRES-ILU(3) method is much higher than that of the LU-SGS method, but the base flow is converged with higher precision, reducing numerical noise in the computed fields.

The in contrast to the PS ones, which grow slightly slower in our case. This difference can be explained by both the different convergence rates achieved in both computations and the different numerical modelling. Despite these observations, both tendencies seem similar enough to validate our base flow computations.

Solver

Due to the distinct characteristics of AS and PS base flows, the rest of this study is organized to analyze the modal and non-modal stability of each type of base flow separately. The following sections present a global stability analysis around each of these distinct base flows.

IV.3.3 Regular bifurcation of axi-symmetric base flow

Direct modal analysis. The selected preconditioner for solving the K-formulated linear systems involved in the global stability analyses performed on the laminar 3-D sphere case was the FGMRES-ILU(3) with one RAS level. The computations were carried out on the same cluster as the fixed point computations, but 43 nodes, consisting of 1200 CPUs, were required for the modal analyses to meet the memory requirements of the preconditioner. Within the GMRES, the internal and external Krylov bases consisted of 100 and 30 vectors, respectively. The internal tolerance was set to 10 -3 and the external one was set to 10 -6 . The shift was set to s = 1, and the KS basis was composed of 22 vectors, with a relative tolerance for convergence of eigenpairs of 10 -4 . The least stable eigenpair converged below the tolerance in about 8 KS iterations.

The eigenvalues are shown in fig. IV.20. The onset of symmetry breaking is determined to occur between Re = 200 and Re = 210, which is consistent with previous studies by Sansica [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF], Tezuka & Suzuki [START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF], and Fabre et al. [START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF]. The least stable direct eigenfunction at Re = 210 is depicted in fig. IV.21, normalized by the L ∞ -norm of the perturbation velocity vector. To better highlight the azimuthal nature of the eigenfunctions, a transformation from cartesian coordinates (x, y, z) to cylindrical coordinates (x, r, θ) has been applied, yielding the transformed perturbation velocity components: The steady mode's spatial structure is consistent with previous findings by Sansica [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF], Fabre [START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF] and Tezuka & Suzuki [START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF]. As noted by Fabre, the mode displays reflectional symmetries about longitudinal planes. The dynamic is dominated by the streamwise velocity component, which expands and grows over the domain at x = 20 with infinitely large scale structures due to its non-oscillatory character.

Adjoint modal analysis. The computation of adjoint eigenpairs of AS base flows was performed using the same parameters as for direct modal analysis, except for the tolerance on the eigenpairs, which was set to 10 -3 . The convergence was reached in approximately 5 KS iterations. The resulting least stable adjoint eigenvalues are shown in Fig. IV.20. A good agreement is observed between the direct and adjoint eigenvalues. The duality between the direct and adjoint operators has already been validated in the laminar case using a cylinder, but this result shows that the duality is also valid in three dimensions. The perturbation velocity components of the corresponding unstable eigenfunction in cylindrical coordinates at Re = 210, normalized by the L ∞ -norm, are shown in fig. IV.23. Unlike the direct mode, the dynamic of the adjoint mode is dominated by radial and streamwise perturbation velocity components. The main part of this mode is localized immediately after the sphere body and extends up to x ≃ 2. It also shares the same symmetry properties as the direct mode. To delay the breaking of symmetry for this AS base flow, it is hence important to consider the strong receptivity of the direct mode to external forcing in these regions.

The shape of the structural sensitivity calculated using the direct and adjoint unstable steady modes is displayed in fig. IV.24. It is seen to be plane-symmetric across the x-y and x-z planes, and it indicates that a symmetry break may occur when one of the two dark lobes visible between x = 0.5 and x = 2 is perturbed. The most significant values are observed in the x-z plane. The streamlines of the corresponding AS base flow have been added to help better understand the localization of the wavemaker in the beginning of the recirculation bubble. The two dark lobes in the recirculation bubble are therefore the zones directly responsible for the instability. As previously mentioned, the computational mesh must be small enough in this region. This quantity is of great interest in an industrial context as it enables the reduction of the mesh size while still capturing the instability. To the author's knowledge, both the adjoint eigenfunction and structural sensitivity for this first bifurcation have not been shown in the literature.

Having characterized the first bifurcation of the AS base flow, we will focus in the following on a second bifurcation that occurs in this flow. This second bifurcation is a Hopf bifurcation and results in oscillatory modes.

IV.3.4 Hopf bifurcation of axi-symmetric base flow

Above Re = 210, the AS base flow has already undergone a bifurcation towards a planesymmetric configuration, which, in turn, becomes unsteady through a supercritical Hopf bifurcation around Re = 280 (see fig. IV.3.6 and [107,[START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF]). As stated in the introduction, an AS base flow computed above Re = 280 is not physically meaningful. However, as the Reynolds number is increased, a pair of complex conjugate oscillatory unstable modes appear, related to a Hopf bifurcation. Thus, two types of instability coexist consisting in one steady mode and two unsteady ones. This situation mimics more complex configurations in which multiple bifurcations can occur simultaneously and demonstrates the robustness of our approach in computing them. Being able to compute each of these modes can help to understand the scenario of transition from steady to unsteady. To this end, direct and adjoint modal oscillating modes are calculated for the cases Re = 240, Re = 260, and Re = 280.

Direct modal analysis. The KS bases used in these computations consisted of 30 vectors and half of the basis was reused from one KS iteration to the next. The shift was set to s = 0.5 + 0.82i in order to target the oscillatory unstable mode, which, based on resolvent analysis (to be presented in IV.3.5), should have a frequency close to St = 0.13. The errors associated to the unstable oscillatory eigenpairs drop below 10 -4 in 4 KS iterations. The computations yield both steady and unsteady modes. between Re = 260 and Re = 280. This is an interesting result since the Hopf bifurcation of the PS base flow also occurs in this range of Reynolds numbers (see fig. IV.3.6). To the author's knowledge, this result has never been reported in the literature.

Isosurfaces and sectional contours of the velocity components, normalized by the L ∞ norm, are shown in fig. IV.26 for Re = 280. Since each mode computed using linear stability analysis is defined by one phase and one amplitude, and due to the axisymmetry of the base flow, we chose the phase so that the symmetry plane aligns with the x-z plane. The mode consists of small-scale structures that are advected downstream from the sphere body until x = 20. A symmetry pattern across the longitudinal direction is also observed. It is worth noting that the mode has a similar spatial structure to the directly unstable mode computed by Sansica in [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF] and Citro in [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF] for the planar-symmetric base flow. The steady mode is not displayed here, but it has a shape that is very similar to that calculated at Re = 210 and shown in fig. IV.21.

Adjoint modal analysis. For these adjoint computations, the same parameters were used as for the direct ones. As previously demonstrated for the direct case, the spectra corresponding to the specified Reynolds number are shown in Fig. IV.27. The adjoint eigenvalues match the direct eigenvalues exactly, further confirming the duality between the direct and adjoint eigensolvers.

Isosurfaces and sectional contours of the perturbation velocity components, normalized by the L ∞ -norm, are shown in fig. IV.28. The most energetic part of the mode is located in the recirculation bubble downstream of the spherical body. This eigenfunction shares similarities with the adjoint mode calculated by Citro in [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF] for the Hopf bifurcation of PS base flows. Similar to the direct mode, a symmetry pattern along the longitudinal direction is observed. However, unlike the direct mode, this adjoint mode is convected upstream in the flow domain, which is analogous to the results from the laminar cylinder case (see III.6.2).

The structural sensitivity calculated from the direct and adjoint oscillatory modes have been computed and the resulting spatial structure is shown in with the spatial distributions of both the direct and adjoint modes that are dominated by the streamwise velocity, which is anti-symmetric across the x-y plane, but symmetric across the x-z plane. Once again, to the best of the author's knowledge, these results have not been previously shown in the literature.

To conclude the global stability analysis of the AS base flows, a resolvent analysis is performed and described in the following on a subcritical AS base flow. As demonstrated below, even far from the Hopf bifurcation, a peak in the energetic gain curve is observed, indicating a strong receptivity of the flow at the frequency corresponding to the oscillatory axisymmetric mode.

IV.3.5 Resolvent analysis of axi-symmetric base flow

This study focuses on the subcritical AS base flow at Reynolds number 200. Like the previous two test cases in this work, the mass matrices used to solve the eigenvalue problem are selected as

Q 1 = Q Chu and Q 2 = Q Ω .
As previously stated, the resolvent analysis requires the computation of two preconditioners for the application of

R K Q K Ω -1 R H Q K Chu .
The number of CPUs used for direct and adjoint computation is not sufficient to meet the memory requirements. Thus, there are two possible strategies: either we store only one preconditioner at a time, meaning that for each inversion, one preconditioner is formed and the old one is erased, resulting in lower memory requirements but more CPU time; or we can double the number of processors to increase the total available memory to meet the memory requirements for storing two preconditioners. Given the imposed usage of ONERA's clusters, we chose the second strategy. This results in computations run on 2800 CPUs, using FGMRES-ILU(3) preconditioning with one RAS level. Inversion of the adjoint operator is stiff in this case, so both Krylov bases have been composed of 100 vectors to ensure inversions fall below the tolerance of 10 -6 . The internal tolerance was set to 5 × 10 -2 . Based on the experience of the OAT15A case (see IV.2.5), the KS basis is composed of only 9 vectors, resulting in 18 inversions at each KS iteration. These parameters are sufficient to converge the most energetic gain with a relative precision of 5 × 10 -4 , which is usually well separated from the others. The total memory used during the computations is about 12 TB, and the CPU time is about 9 hours.

The optimal gain for various St is shown in fig. IV.30. There is no peak near zero-frequency, unlike what was expected based on previous modal analysis results. This might be improved by 0.08 0.12 0.16 0.2 St conducting resolvent analysis closer to St = 0, but this would require solving ill-conditioned linear systems that may not be possible with current tools. Tests were conducted at St = 10 -4 , but the GMRES method did not converge. However, a peak is seen near St = 0.13, indicating strong flow receptivity at this frequency, which corresponds to that of the previously calculated unstable oscillatory modes. This means that this mode can be expressed as a convective instability by exciting this frequency, even in a sub-critical regime where the base flow is still AS.

IV.3.6 Hopf bifurcation of planar-symmetric base flow

In the following, the KFGSM is used to predict the onset of the Hopf bifurcation on the PS base flow presented in IV.3.2. This bifurcation has been well described in the literature, and the results of the direct and adjoint global stability analyses are compared to those obtained by Sansica in [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF] and Citro in [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF].

Direct modal analysis. The inversion parameters used in the KFGSM are the same as those used in direct and adjoint computations around the AS base flows. The KS basis consisted of 22 vectors, and the shift was set at s = 0.5 + 0.82i to target oscillatory unstable modes. For each Reynolds number, the relative error on the least stable eigenpairs decreased below 10 -4 after approximately 8 KS iterations. Fig.

IV.33a and IV.33b respectively show the evolution of the growth rate, and the Strouhal number with respect to the Reynolds number. The values of σ obtained are in good agreement with those calculated by Sansica in [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]. However, some discrepancies are observed in the evolution of St between the elsA results and those of Sansica and Citro. This may be attributed to different numerical schemes and levels of convergence. Nevertheless, the overall trend is similar and validates the curve obtained in elsA, which compares well with the results obtained by both Sansica and Citro [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF]. The instability threshold is found to be Re = 280, which is consistent with the literature [START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]. These results obtained from the linear solver are further validated by nonlinear unsteady computations presented in IV. The structural sensitivity is computed from both the direct and adjoint unstable modes at Re = 280, and the resulting isosurfaces and sectional contours are presented in fig. IV.37. The shape is consistent with the result obtained by Citro [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF]. Unlike the structural sensitivities calculated for the AS base flow, only a symmetry across the x-z plane is observed. The wavemaker is located immediately downstream of the spherical body, at the bottom of the recirculation vortex. It differs from the wavemakers identified in the AS case.

As with the AS base flow, this study can be concluded by analyzing the receptivity of the flow to external forcings. To do so, a resolvent analysis is performed on a sub-critical PS base flow, and the results are discussed below.

IV.3.7 Resolvent analysis of planar-symmetric base flow

As previously identified with modal analysis and further confirmed with unsteady computations, the PS base flow evolves toward an unsteady solution above Re = 280. The mode responsible for this unsteadiness has been calculated, as well as the precise region where the instability is triggered, thanks to the structural sensitivity. The natural extension of this work is to conduct a resolvent analysis.

The analysis is conducted on a stable PS base flow computed for Re = 260. The mass matrices defining the resolvent operator are, again,

Q 1 = Q Chu and Q 2 = Q Ω .
Both the KS parameters and the number of CPUs used are the same as those used for the resolvent analysis of the AS base flow (see IV.3.5). Similar to the latter, the frequencies where the linear systems were most difficult to solve were those close to St = 0. For instance, for the case with St = 0.025, an average of 95 external GMRES iterations were required during inversions of the adjoint part of the resolvent operator to reduce the error below 10 -6 .

The resulting gain curve is presented in Strouhal number of the least stable oscillating mode calculated at this Reynolds number is 0.133. It turns out that a convective instability might be triggered by exciting the stable base flow at this frequency. The shapes of both the optimal forcing and response calculated at St = 0.13 are respectively shown in fig. IV.39 and IV.40. As expected, the optimal forcing has the same spatial characteristics as the adjoint eigenmode calculated at Re = 280 (see IV.36). The optimal forcing is mainly localized immediately after the spherical body, is convected upstream, and is dominated by streamwise velocity, which is symmetric across the x-z plane. It highlights the zones most sensitive to external forcing to trigger or control this convective instability. The optimal response, in turn, is similar to the unstable direct eigenmode calculated at Re = 280 (see IV.34). It is convected downstream and shares the same symmetry pattern as the optimal forcing. The spatial wavelength for this eigenfunction is Λ x ≃ 6.1L, which is of the same order as that measured for the direct eigenmode at Re = 280. This description of the resolvent analysis carried out on a sub-critical PS base flow completes the global stability analysis of the subsonic laminar 3-D sphere flow. The last step was to perform unsteady computations to compare the results with the modal analysis, which is the aim of the following.

IV.3.8 Unsteady computation of planar-symmetric base flow.

Unsteady nonlinear calculations are conducted at Reynolds numbers of Re = 300 and Re = 320. The physical non-dimensional time-step is set to ∆t = 10 -3 , resulting in a CFL number around 10 at the wall. As for the two precedent test cases, the Gear method was used for these computations. The base flow, perturbed by the oscillating direct mode, is used as the initial condition. A probe in the wake of the sphere measures the perturbation velocity signal, as shown in fig. IV.41. The temporal evolution of the signal indicates an exponential growth until reaching the limit cycle around t = 85.

St

(σ L , St L ) (σ NL , St NL ) ∆St ∆σ 300 (0.0617, 0.1344) (0.0615, 0.1371) 1.9% 3.2% 320 (0.1113, 0.1354) (0.1052, 0.1373) 1.4% 5.8% Table IV.9: Amplification rate and Strouhal number of linear and nonlinear computations. sequence of interconnected coherent structures extends and grows far from the sphere boundary at x = 20. Additionally, the spatial wavelength Λ x ≃ 6.5L is very close to that identified for the direct unstable mode at Re = 280.

Both σ NL and St NL are extracted from the perturbation velocity signal and shown in fig. IV.43a and IV.43b, respectively. The amplification rates obtained from linear calculations are very close to those of nonlinear computations. The same applies to the Strouhal number, for which a relative difference of only 1.4% is observed, even far from the instability threshold at Re = 320. As stated by Sansica [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF], this finding contrasts with the case of the laminar cylinder flow, where the prediction of shedding frequency diverges quickly as the Reynolds number increases (see table III.8). While the linear stability theory is less predictive for the amplification rate, it may be enhanced by reducing the tolerance of the eigensolver.

Before concluding this section on the laminar sphere flow, some data concerning the scaling up process for this case will be provided below.

IV.3.9 Data about the scaling up process

As mentioned in the introduction of this chapter, a finer mesh consisting of about 25 × 10 6 cells (representing 125 × 10 6 DOF) was also considered. In particular, direct and adjoint computations were carried out on the PS base flow at Re = 280, with the same parameters as for the regular mesh. For both the direct and adjoint unstable modes, the relative error in the associated eigenvalue with respect to that computed with the regular mesh was 4.4%. This small error validates the global stability computations on this refined mesh.

Table IV two types of mesh, with the theoretical memory footprint of the corresponding LU factorizations, based on the estimations built in III.1.1. A factor of 50 is identified between the memory footprint of the KFGSM and the LU factorization. This result is satisfying as it allows us to conduct global stability computations on such 3-D configurations. As a reminder, and to the author's knowledge, the largest global stability analysis ever performed with the type of method used in this thesis was carried out by Timme [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF] on a 3-D plane configuration composed of about 50 × 10 6 DOF.

IV.4 Conclusion

This chapter aims to demonstrate the capabilities of the KFGSM on challenging test cases. First, the code is validated on the RANS case of the flow around an OAT15A airfoil. Fixed points of this stiff problem are computed with high accuracy, and the method performs well for direct and adjoint modal analyses. The compatibility with the linear solver and nonlinear URANS computations with a negative-SA turbulence model is also verified. Resolvent analysis is carried out and characterizes both the buffet phenomenon for sub-critical base flow and the secondary Kelvin-Helmholtz instability, as tackled by Sartor [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF].

The 3-D case of the laminar flow around a sphere body is then addressed. The computational domain for this case is composed of approximately 10 × 10 6 DOF. AS and PS base flows are studied. The code allows the prediction of regular and Hopf bifurcations for the AS base flow. A resolvent analysis at sub-critical Reynolds numbers illustrates the similarity between optimal forcing and response and adjoint and direct modes for the Hopf bifurcation. New results concerning adjoint steady mode of AS base flow, and the corresponding structural sensitivity, were demonstrated. A complete description of the Hopf bifurcation occurring in AS base flow at relatively high Reynolds number is also provided, which, to the author's knowledge, have not been shown in the literature.

Finally, the data provided to illustrate the scaling up of the KFGSM demonstrate its efficiency in terms of memory storage relative to classical global stability methods involving LU factorization. The global stability analysis was successfully scaled up, which was the main objective of this thesis.

Chapter V Conclusion

As discussed in chapter I, global stability analysis is a powerful tool for improving the understanding of unsteady flows. After a base flow has been computed, this method may be seen as a thorough post-processing that returns modes of the flow that may be linearly amplified and potentially responsible for the unsteadiness. In particular, it could be used to predict the appearance of the rotating stall phenomenon occurring in axial compressors like the CME2, and thus help aircraft manufacturers to adapt their design to delay instability, reduce the surge margin, and improve performance. However, this methodology requires the inversion of many linear systems. When addressing large-scale problems exceeding 10 × 10 6 DOF, inversions cannot be managed using direct methods like LU factorization, as is often the case [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF][START_REF] Theofilis | Global Linear Instability[END_REF].

The purpose of this thesis was to develop a robust method capable of scaling up global stability analysis to industrial cases like the CME2 compressor. To tackle this problem, the main constraint was to find a way to solve extremely large complex-valued linear systems while preserving highperformance computing optimizations, with the knowledge that complex arithmetic cannot be vectorized. To do so, a new method was designed in which complex problems were reformulated in the real-equivalent K formulation [START_REF] Day | Solving Complex-Valued Linear Systems via Equivalent Real Formulations[END_REF], a formulation that preserves the sparsity pattern of the original complex operators, and where inversions are managed by the iterative preconditioned GMRES algorithm [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF][START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. A mathematical background was established to build the method, named the K-Formulated Global Stability Method (KFGSM). Even not presented in this work, the latter was first validated on a toy code that solved the global stability problems for the nonlinear advection-diffusion Burgers equation [START_REF] Burgers | A mathematical model illustrating the theory of turbulence[END_REF].

To implement the KFGSM in the elsA software, several elementary bricks needed to be assembled. The first step was to adapt the existing preconditioned GMRES based on obtaining Jacobian-vector products using Algorithmic Differentiation (AD) through the Tapenade software developed by Hascoët and Pascual [START_REF] Hascoet | The Tapenade Automatic Differentiation tool: principles, model, and specification[END_REF], and the efficient preconditioning tools developed by Guilbert [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF] to solve the real-valued problems arising in the computation of fixed points in compressible CFD. A collaboration with Hascoët led to a new feature of Tapenade that computes the forward and backward K-formulated Jacobian-vector products. The preconditioning tools of Guilbert were also adapted so that the resolution of linear systems with preconditioned GMRES is operative on K-formulated problems. The final step was to implement the Krylov-Schur eigensolver in real arithmetic, where the Schur matrix is in real form, incorporating the K-formulated GMRES. Slight modifications in the convergence test needed to be carried out to ensure the algorithm returned the eigenpairs of the original system and not the K-formulated one. The algorithm have been developed in three versions corresponding to the resolution of eigenproblems arising in direct and adjoint modal analysis and resolvent analysis.

The method was challenged on three test cases of gradual difficulty. The simplest is the laminar flow around a cylinder. Fixed points at various Reynolds numbers were computed with Flexible GMRES preconditioned with ILU(3) factorization. A good convergence near machine precision is obtained, demonstrating the robustness of the fixed point computation method developed by [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. The direct, adjoint, and resolvent analysis conducted with the KFGSM are in good agreement with the literature [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF], validating the method on the laminar case and demonstrating its capability to predict the onset of von K'arm'an vortex street. Moreover, a comparison with unsteady computations showed the compatibility between linear and nonlinear solvers.

The second problem addressed is the RANS case of the flow around the OAT15A airfoil. This case has been selected due to its high stiffness [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. Fixed points were computed in sub and supercritical regimes with high accuracy. The KFGSM showed good results at predicting the buffet onset on this RANS case with direct and adjoint modal analysis. Such results were also successfully compared with URANS computations. From the resolvent analysis, two energy peaks were observed. The largest one corresponds to the buffet frequency, while the second is linked with a secondary instability that is a Kelvin-Helmholtz one [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], illustrating the capability of the resolvent analysis to study in detail the receptivity of the flow within a large range of frequencies.

The code was ultimately tested on the three-dimensional subsonic laminar flow around a sphere body. The computational grid was composed of about 10 × 10 6 degrees of freedom (DOF), representative of the industrial cases the KFGSM has been designed for. After fixed point computation near machine precision were carried out for axi and planar-symmetric base flows, the two bifurcations were studied in detail with direct, adjoint, and resolvent analysis. In particular, our fixed point computation method allowed the calculation of axi-symmetric base flow above a second bifurcation that, to the author's knowledge, had never been documented. The resolution of both direct and adjoint eigenproblems allowed the calculation of the structural sensitivity [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. The second bifurcation occurring on the planar-symmetric base flow was also successfully characterized. On this particular type of base flow, a global stability analysis was even performed on a refined mesh composed of about 125 × 10 6 DOF. Such results validated the scale up of global stability analysis, as the largest reported global stability case is the characterization of the buffet phenomenon on a full airplane configuration composed of about 50 × 10 6 , carried out by [START_REF] Timme | Global instability of wing shock-buffet onset[END_REF].

Chapter VI Perspectives

Various perspectives can be proposed for this work. At ONERA, important work is currently being carried out to improve the GMRES solver. For instance, Jang [START_REF] Jang | Randomized flexible GMRES with Singular Vectors Based Deflated Restarting[END_REF] has improved the formation of the Krylov basis by implementing the randomized Gram-Schmidt process in elsA. This process can be less computationally expensive than the classical Gram-Schmidt process, while being at least as numerically stable as the modified Gram-Schmidt process used in this work. Another perspective for improving the GMRES algorithm is to improve the preconditioning tools. As seen with the laminar flow around the sphere body, the memory cost for the storage of the preconditioner remains high when dealing with large problems. For problems stiffer than the sphere case, ILU(k) factorization may not be sufficient to allow the GMRES to converge, and more advanced techniques like the HLU may need to be used. However, the HLU method is currently unusable on the ONERA cluster for such large cases. To circumvent the potential memory issue related to preconditioning, a thesis is currently being carried out by Dubois at ONERA. The aim is to reduce the memory cost of preconditioning.

Improvements to the Krylov-Schur solver can also be cited. In fact, given the time limit imposed on the computing cluster at ONERA, a job may not be sufficient to carry out the entire global stability analysis. To this end, it could be useful to include restarting so that the Krylov-Schur basis can be extracted at the end of the computation if not all the desired eigenpairs have yet converged. This basis could then be read in a subsequent job to restart the computation where it stopped in the previous job.

Another natural perspective of this work is to perform a global stability analysis on the isolated compressor CME2 to characterize the genesis of the rotating stall phenomenon. The full-annulus case represents a domain of about 30 × 10 6 cells. The usual boundary conditions for this case are an injection condition at the inlet, no-slip at walls, and a radial equilibrium with a valve law at the outlet. However, such a radial condition is not yet linearized in the elsA software. To circumvent this difficulty, computations have been carried out with an extended domain representing about 40 chords, and an outflow condition has been applied at the outlet. To check the validity of this configuration, compressor characteristics have been calculated by conducting RANS computations with the LU-SGS method (which doesn't require the linearized matrix-vector products), where the convective fluxes are discretized with a second-order Roe scheme, the diffusive fluxes with the 5p cor method, and the turbulence model is the negative Spalart-Allmaras. The results are given in fig. VI.1. The pressure distribution is slightly lower in the extended vein with outflow condition compared with the normal vein with radial equilibrium condition. However, until the stability limit occurring around Q = 4 kg s -1 , the extended vein with outflow condition is in good agreement with the extended vein with radial equilibrium. Based on these positive results, a URANS computation performing a stability analysis on only three channels.

The continuation of this work is to perform a fixed point computation to obtain a converged base flow with high accuracy. This may be challenging because, in the author's opinion, the instability may involve the rotation of the stall cells. Thus, the unstable base flow has to be already detached and could result in a stiff problem. A possibility is to converge a detached base flow on a single channel case, to duplicate it and use the result as the starting field for the base flow computation. Another possibility would be to benefit from the periodicity of the compressor CME2 and use the Schmid [START_REF] Schmid | Stability analysis for n-periodic arrays of fluid systems[END_REF] approach. The latter, already exploited by Paladini in [START_REF] Paladini | Transonic buffet instability: From two-dimensional airfoils to three-dimensional swept wings[END_REF], allows the computation of eigenpairs of the full-annulus compressor by carrying a stability analysis on only three channels.

More globally, the tools developed in elsA in this work allow for the analysis of modal and nonmodal phenomena in an industrial context. This will allows for the characterization of oscillating phenomena on a linear basis, which is very interesting for predicting undesirable unsteady phenomena. This work also provides robust tools to characterize noise-amplified phenomena through resolvent analysis. In a flow, multiple unsteadiness may coexist with highly variable time-scales. A computation that captures each phenomenon requires performing with all time-scales, which could result in very long and costly calculations. The resolvent is hence of interest for many applications involving super or hypersonic flows and could improve the knowledge about the laminarto-turbulent transition scenario on such configurations. One can cite, for instance, the studies currently carried out at ONERA on the BOLT configuration, presented in fig. VI.3, which was expressly designed to study the laminar-to-turbulent transition on hypersonic configurations. été vraiment formateur pour moi. Durant cette soutenance, la salle Contensou du site de Châtillon était presque pleine. Lorsque j'ai envoyé des invitations pour cet évènement début janvier, j'étais très surpris et touché de voir que 26 personnes allaient répondre présentes. D'abord, les potes de Lyon, ma deuxième famille. Céd, Mamoule, Jé, Anto, Rapace, Rapacette, Camouille, Laure, vous avez fait le déplacement jusqu'à Châtillon city pour venir me soutenir et c'est quelque chose qui restera à jamais gravé dans mon coeur. Ensuite, tous les potes de la Kill et proches affiliés. Dufal, tu crois qu'on savait dans quoi on se lançait il y a 5 ans, au détour d'une discussion au Bar'Bu, quand on a décidé quasiment ensemble de faire une thèse ? Je n'oublierai pas ces trois premiers mois de thèse, incroyablement productifs, passés dans le même bureau, de même que toutes nos discussions méca flu. Tiens bon, tu n'es plus qu'à quelques mois de finir ta thèse !! Je remercie bien chaleureusement tous les autres qui ont fait le déplacement pour être avec moi ce jour-là : Bib, Louis, S'il vint, Gorce, Ch'Arvor, Larrin, Marblé, Thoyus, Rochette, Nounouck, Lucho. Je remercie aussi ceux qui n'ont pas pu être là mais qui m'ont probablement soutenu à distance comme Étienne, Saman, Dr. Labat ou encore futur Dr. Beuz'. Vous êtes tous de vrais amis sur qui je pense toujours pouvoir compter.

J'en arrive à ma famille que j'aimerais remercier pour son soutien indéfectible. Caro, Ana, vous avez toujours été des pilliers dans ma vie, des sources de soutien, de sagesse, de partage de connaissances mais surtout d'amour inconditionnel. Je vous aime si fort. Je pense également à vos familles respectives : Charlotte, dont j'ai eu la chance d'être le coloc pendant 2 mois, votre petite Sashounette qui, j'en suis sûr, ira loin, mais aussi François, mon grand Malo et ma petite Lou. Ensuite, Manu et Laulau, mes p'tits rarons (sans oublier Roro, Yummi et Car bien sûr), vous avez été un bol d'air incroyable lorsque nous nous voyions régulièrement à la fin de la rédaction de mon manuscrit, et le fait que vous veniez me voir le jour J m'a tellement touché. Je vous remercie de tout mon coeur. J'ai une grande pensée pour tous mes oncles, tantes, cousins et cousines qui, pour la plupart, m'ont suivi à distance. Vous n'étiez pas là physiquement, mais bien présents dans mon esprit. Pap, Mam, j'ai beaucoup de chance de vous avoir. Vous êtes des épaules sur qui toujours me reposer en cas de doutes, des soutiens infaillibles, et avec vous, il y a toujours une solution à tous les problèmes ! Je vous aime très fort. Je pense aussi à Jean-Phi et Barbara, merci d'être vous et d'avoir été là pour moi depuis toutes ces années. J'ai également une pensée ému pour mes quatre grands-parents qui, j'en suis sûr, seraient heureux que leur petit-fils soit devenu docteur.

Pour finir, Anne. La réussite de cette thèse, c'est un peu la tienne tant nous avons été proches depuis tout ce temps. Elle a pris beaucoup de place, et je t'en ai fait baver, mais tu as toujours été là pour faire redescendre la tension et dissiper mes angoisses. La bonne tenue de cette soutenance, je la dois en grande partie à toi et je t'en remercie sincèrement. Tu as été un phare constant dans ma vie et je ne l'oublierai jamais. Sache qu'à ton tour, quoiqu'il arrive, tu pourras toujours compter sur moi.

J'étais sûr en les commençant que ces remerciements ne tiendraient pas sur une page. Ce n'est pas grave, je suis content d'avoir pris le temps d'avoir une pensée pour vous tous. Peu importe où la vie me mènera, je sais que ne suis pas seul, et ça c'est quelque chose de très précieux. Désormais, il est l'heure pour moi de prendre un repos que je pense dûment mérité .

Valentin Fer, le 14 février 2023

Figure I. 1 :

 1 Figure I.1: Illustrative representation of the rotating stall phenomenon occuring in a compressor that rotates at speed Ω.

Figure I. 2 :

 2 Figure I.2: Compressor stall vizualization on a Sukhoi Su-57 ptototype, by Rulexip via Wikipedia Commons.

Figure I. 3 :Figure I. 4 :

 34 Figure I.3: Illustration of stable (a), unstable (b) and neutral (c) equilibrium.

w 1 Figure I. 5 :

 15 Figure I.5: Geometric illustration of the temporal evolution of ∥w ′ ∥ in the normal (left) and non-normal (right) cases.

Figure II. 2 :

 2 Figure II.2: Convergence curve (a) and initial, theoretical and converged field (b) after 16 iterations for the Newton method.

Figure II. 3 :

 3 Figure II.3: Convergence curve (a) and initial, theoretical and converged field (b) for the implicit backward Euler with adaptative CFL method.

Figure II. 4 :

 4 Figure II.4: Spectrum of J (a) computed by SLEPc and LAPACK and real part of eigenmode calculated by SLEPc (b) corresponding to the least stable eigenvalue λ = -4.55 × 10 -6 .

Figure II. 5 :

 5 Figure II.5: Spectrum of J † (a) computed by SLEPc and LAPACK and real part of eigenmode calculated by SLEPc (b) corresponding to the least stable eigenvalue λ = -4.55 × 10 -6 .

Figure II. 6 :

 6 Figure II.6: Optimal gain η 1 (ω) (a) and real parts of corresponding optimal forcing and response (b) at ω = 10 -3 Hz.

Figure II. 7 :

 7 Figure II.7: Illustration of shift-and-invert transformation.

Contents III. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.1.1 Scale up issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.1.2 Numerical set-up of laminar cylinder flow . . . . . . . . . . . . . . . . . . . III.2 Equations of compressible CFD . . . . . . . . . . . . . . . . . . . . . . . . . . III.2.1 Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.2.2 Reynolds Averaged Navier-Stokes equations . . . . . . . . . . . . . . . . . III.3 elsA CFD solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.3.1 Unstructured computation . . . . . . . . . . . . . . . . . . . . . . . . . . . III.3.2 Numerical discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.3.3 LU-SGS method for fixed point computation . . . . . . . . . . . . . . . . . III.3.4 GMRES algorithm for fixed point computation . . . . . . . . . . . . . . . . III.3.5 Fixed point calculation of the laminar cylinder flow . . . . . . . . . . . . . III.4 K formulation and some of its properties . . . . . . . . . . . . . . . . . . . . III.4.1 Real equivalent formulation of a complex problem . . . . . . . . . . . . . . III.4.2 Eigenvalues of the K-formulated operator . . . . . . . . . . . . . . . . . . . III.4.3 Eigenvectors of the K-formulated operator . . . . . . . . . . . . . . . . . . III.5 K-formulated global stability analysis method . . . . . . . . . . . . . . . . . . III.5.1 K-formulated matrix-vector product . . . . . . . . . . . . . . . . . . . . . . III.5.2 GMRES algorithm in K formulation . . . . . . . . . . . . . . . . . . . . . . III.5.3 Adaptating the Krylov-Schur algorithm to the K formulation . . . . . . . . III.6 Global stability analysis of the laminar cylinder flow . . . . . . . . . . . . . . III.6.1 Direct modal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.6.2 Adjoint modal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.6.3 Resolvent analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.7 Unsteady computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.7.1 Gear method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III.7.2 Nonlinear unsteady computation of the laminar cylinder flow . . . . . . . . III.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure III. 1 :

 1 Figure III.1: Schematic representation of computing components.

Figure III. 2 :

 2 Figure III.2: Schematic view of the cylinder mesh (left) and computational grid (right).

Figure III. 3 :

 3 Figure III.3: Example of unstructured mesh. Grey zone highlights the cell Ω k of center x c k . Interface I j is the fronter between Ω k and Ω l .

Algorithm 3

 3 Basic GMRES Input: A, b, x 0 , tol and m Output: x m minimizing ∥r m ∥ 2 with actualized m 1: r 0 = b -Aw 0 2: β = ∥r 0 ∥ 2 3: v 1 = r 0 /β 4: for k ∈ 1, m do 5: v k+1 ← Av k 6:

Algorithm 4 3 : v 1 = r 0 /β 4 :

 4314 Right-preconditioned GMRES Input: A, M, b, x 0 , tol and m Output: x m minimizing ∥r m ∥ 2 with actualized m 1: r 0 = b -Aw 0 2: β = ∥r 0 ∥ 2 for k ∈ 1, m do 5:

54 ) 3 : v 1 = r 0 /β 4 : u 1 = M -1 v 1 5 :

 5431415 and the new participant to the external Krylov basis writes v j+1 = Au j+1 . The set of internal Krylov vectors are stored in the basis U m to recover the solution at the end of the main loop:x m = x 0 + U m y m .(III.55)Algorithm 5 presents an implementation of FGMRES. This variation usually shows significant improvement in terms of convergence speed compared to the basic right-preconditioned GMRES. However, it requires the storage of a second Krylov basis. The dimensions of both internal and Algorithm 5 FGMRES Input: A, M, b, x 0 , tol and m Output: x m minimizing ∥r m ∥ 2 with actualized m 1: r 0 = b -Aw 0 2: β = ∥r 0 ∥ 2 for k ∈ 1, m do 6:

Figure III. 4 :

 4 Figure III.4: Jacobian pattern of the laminar cylinder (0.9% of nonzero elements).

a

  11 0 a 13 0 a 22 0 0 a 32 a 33    -→ CSR AA(A) = a 11 a 13 a 22 a 32 a 33 JA(A) = 1 3 2 2 3 IA(A) = 1 3 4 .

Figure III. 5 :

 5 Figure III.5: Comparison of preconditioned GMRES

59 )Figure III. 6 :

 596 Figure III.6: Convergence curve (a) and streamlines (b) for the Re = 47 case. The solid black line indicates the separated region.

Figure III. 7 :

 7 Figure III.7: Streamwise (a) and cross-stream (b) velocities of the baseflow computed at Re = 47. The solid line indicates the recirculation bubble.

Figure III. 8 :

 8 Figure III.8: Evolution of recirculation bubble length with Reynolds number.

c 11 0 c 13 0 c 22 0 0 c 32 c 33 = c 11 c 13 c 22 c 32 c 33 JA

 3333

1 3 4 (

 4 III.66)The block sparsity pattern of C K exactly matches the point sparsity pattern of C, and the CSR storage of C has the exact same structure as the BSR storage of C K . Consequently, applying any pattern-based block preconditioning technique like those implemented in the elsA software (ILU(k), local LU, etc.) to C K corresponds to applying the equivalent point preconditioning to C. This statement stays true if c qr are d-by-d block entries, as in elsA, and the block size of the K-formulated operator becomes 2d-by-2d. This is shown in the case of the cylinder in fig. III.9, on which we note that the patterns of J and J K are identical.

Figure III. 9 :

 9 Figure III.9: Patterns of J (left) and J K (right) for the laminar cylinder.

(III. 78 )

 78 The main consequence of relation (III.[START_REF] Luchini | Adjoint Equations in Stability Analysis[END_REF]) is that Sp(C K ) = Sp(C K1 ). Considering a complex eigenvector υ = (υ) 1≤q≤2N ∈ C 2N of C K associated with the eigenvalue λ = α + iβ, we hence know that Pυ is an eigenvector of C K1 verifying:C K1 Pυ =λPυ, (III.79) and that λ is also an eigenvalue of C K1 . Introducing the two vectors υ o = (υ 2q ) 1≤q<N ∈ C N and υ e = (υ 2q+1 ) 1≤q<N ∈ C N (respectively composed of the odd and even components of υ), and based on the definition (III.78) of P: Pυ = υ o υ e . (III.80)

. 81 )

 81 These relations are introduced in order to simplify the next demonstration. The following demonstrated propositions clarify the shape of υ depending on whether λ belongs to Sp(C), to Sp(C), or to Sp(C) ∩ Sp(C) (which is included in R), and show how an eigenvector of C can be recovered from υ. Proposition 2 υ o and υ e are eigenvectors of C associated with λ ∈ C ⇔ υ o = -iυ e .

. 93 )

 93 Assuming λ / ∈ Sp(C) ∩ Sp(C), and thus λ / ∈ R, then β ̸ = 0, and finally υ o ℜ = -υ e ℑ , which completes the demonstration of υ o = -iυ e . We now focus on the right-to-left implication by considering that υ o = -iυ e , i.e. υ o ℜ = -iυ e ℑ and υ o ℑ = υ e ℜ . Pυ is an eigenvector of C K1 so the set of equations (III.87) to (III.90) holds. Replacing υ e ℜ by υ o ℑ and υ e ℑ by -υ o ℜ in equation (III.87) and (III.88) leads to equation (III.82) and (III.83). And replacing υ o ℜ by -υ e ℑ and υ o ℑ by υ e ℜ in equation (III.89) and (III.90) leads to equation (III.84) and (III.85). Finally, the whole system (III.82) to (III.85) is satisfied, which means υ o and υ e are eigenvectors of C . Proposition 3 υ o and υ e are eigenvectors of C associated with λ ∈ C ⇔ υ o = iυ e .

2 :

 2 unitary, m and l ≤ m the number of desired eigenpairs Output: j eigenpairs of B, with j ≤ l 1: initialization of V m ← v 1 0 , j ← 0, p ← 0 while j ≤ l until convergence do 3:

Figure III. 11 :

 11 Figure III.11: Evolutions of amplification rate (a) and Strouhal number (b) according to Re for direct modal analysis. The critical Reynolds number is Re c = 47.

Figure III. 12 :

 12 Figure III.12: Real part of streamwise (a) and cross-stream (b) velocity perturbation of the direct unstable mode at Re = 47. Solid lines indicate iso-contours at ℜ( u x ) ± 0.1 and ℜ( u y ) ± 0.1.

Figure III. 13 : 5 Table III. 6 :Figure III. 14 :

 135614 Figure III.13: Evolution of amplification rate (a) and Strouhal number (b) according to Re for adjoint modal analysis. The critical Reynolds number is Re c = 47.

Figure III. 15 :

 15 Figure III.15: Structural sensitivity of the laminar cylinder flow at Re = 40. The dark zones indicates the wavemaker.

Figure III. 16 :

 16 Figure III.16: Optimal gain obtained for the resolvent analysis at Re = 40.

23 Table III. 7 :Figure III. 17 :

 23717 Figure III.17: Real part of streamwise and cross-stream perturbation velocity of the optimal forcing (a-b) and response (c-d) adjoint unstable mode at Re = 40 and St = 0.117.

  Fig. III.18 displays the time evolution of the lift coefficient C L at Re = 55. It grows exponentially until the nonlinear saturation leads to a limit cycle around t = 190, which represents the destabilization of the unstable steady base flow toward the von Kármán vortex street. This can be seen in fig. III.

Figure III. 18 :Figure III. 19 :Figure III. 20 :

 181920 Figure III.18: Time evolution of lift coefficient at Re = 55.

Contents IV. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 IV.2 OAT15A airfoil case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 IV.2.1 Numerical set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 IV.2.2 Fixed point computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 IV.2.3 Direct modal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 IV.2.4 Adjoint modal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 IV.2.5 Resolvent analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 IV.2.6 URANS computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 IV.3 Subsonic laminar three-dimensional sphere flow . . . . . . . . . . . . . . . . 85 IV.3.1 Numerical set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 IV.3.2 Fixed point computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 IV.3.3 Regular bifurcation of axi-symmetric base flow . . . . . . . . . . . . . . . . 89 IV.3.4 Hopf bifurcation of axi-symmetric base flow . . . . . . . . . . . . . . . . . 92 IV.3.5 Resolvent analysis of axi-symmetric base flow . . . . . . . . . . . . . . . . 96 IV.3.6 Hopf bifurcation of planar-symmetric base flow . . . . . . . . . . . . . . . . 98 IV.3.7 Resolvent analysis of planar-symmetric base flow . . . . . . . . . . . . . . . 101 IV.3.8 Unsteady computation of planar-symmetric base flow. . . . . . . . . . . . . 104 IV.3.9 Data about the scaling up process . . . . . . . . . . . . . . . . . . . . . . 105 IV.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Free- stream 1 :Figure IV. 1 :

 stream11 Figure IV.1: Computational domain for the two-dimensional OAT15A airfoil. It expands over 100c.

Figure IV. 2 :Table IV. 2 :Figure IV. 3 :

 223 Figure IV.2: Convergence curve of LU-SGS method initialized with a uniform field (a) and FGMRES-ILU(3) method (b) at α = 4 • .

Figure IV. 4 : 5 Table IV. 3 :

 453 Figure IV.4: Evolutions of amplification rate (a) and Strouhal number (b) according to α for direct modal analysis. Critical angle of attack is found to be α c = 4 • .

Figure IV. 5 :

 5 Figure IV.5: Real part of streamwise velocity perturbation (a) and density (b) for the direct unstable mode at α = 4 • .

Figure IV. 6 :

 6 Figure IV.6: Evolutions of the growth rate (a) and Strouhal number (b) according to α for adjoint modal analysis. Critical angle of attack is found to be α c = 4 • .

Figure IV. 7 :

 7 Figure IV.7: Real part of streamwise velocity perturbation (a) and density (b) for the adjoint unstable mode at α = 4 • .

Figure IV. 8 :

 8 Figure IV.8: Structural sensitivity at α = 4 • .

Figure IV. 9 :

 9 Figure IV.9: Optimal gain evolution according to St at α = 3.75 • . Two peaks are identified at St = 0.0613 (buffet) and St = 2.44 (Kelvin-Helmholtz).

Figure IV. 10 :Figure IV. 11 :

 1011 Figure IV.10: Resolvent analysis at α = 3.75 • and St = 0.0613. (top) ℜ(û x ) (a) and ℜ(ρ) (b) for optimal forcing. (bottom) ℜ(û x ) (a) and ℜ(ρ) (b) for the corresponding optimal response.

7 Figure IV. 12 : 7 Figure IV. 13 : 7 Figure IV. 14 :

 712713714 Figure IV.12: Time evolution of the exponentially growing lift coefficient at α = 4.5 • .

α

  (σ L , St L ) (σ NL , St NL ) ∆St ∆σ 4.5°(0.1243, 0.0638) (0.1248, 0.0657) 2.9% 4% 5°(0.1185, 0.0669) (0.113, 0.0695) 3.7% 4.8% 5.5°(0.0717, 0.0728) (0.062, 0.0733) 0.7% 15% Table IV.6: Comparison between SLEPc and elsA computations for resolvent analysis at Re = 40 and St = 0.1165.

Figure IV. 15 :

 15 Figure IV.15: Schematic bifurcation diagram of the laminar 3-D sphere flow.

  3.3) over 10 4 iterations, performed on 200 CPUs. The resulting field was used to initialize the implicit backward Euler solution with FGMRES-ILU(3) preconditioning. The tolerance on ∥ρ∥ 2 was set to 10 -11 and the convergence, shown in fig. IV.17b, was achieved in 41 iterations. Table

Figure IV. 17 :

 17 Figure IV.17: Convergence of steady base flow at Re = 210 (a) and Re = 280 (b).

Figure IV. 18 :Figure IV. 19 :

 1819 Figure IV.18: Contour map of velocity magnitude for the AS base flow at Re = 210, in the x-y (a) and x-z (c) planes, and the PS base flow at Re = 280 in the x-y (b) and x-z (d) planes. Streamlines are plotted and the dashed-lines represent the recirculation bubble.

Figure IV. 20 : 1 Figure

 201 Figure IV.20: Evolution of amplification rate of the least stable mode according to Re for direct global stability analysis around AS base flow in the regular bifurcation. Symmetry breaking occurs at Re 1 c = 210.

Figure IV. 22 :

 22 Figure IV.22: Evolution of amplification rate of the least stable mode according to Re for adjoint global stability analysis around AS base flow in the regular bifurcation.

1 Figure

 1 Figure IV.23: 3-D views of the adjoint unstable mode at Re = 210, with (σ, St) = (0.02028, 0) for AS base flow. (top) Isosurfaces of streamwise (a), radial (b) and azimuthal (c) perturbation velocities for the values u x = ±0.1, u r = ±0.1, u θ ± 0.1. (bottom) Corresponding sectional contour plot at (d) z = 0, (e) z = 0 and (f) y = 0, with solid-lines indicating iso-contours at u x = ±0.1, u r = ±0.1, u θ ± 0.1.

Figure IV. 24 :

 24 Figure IV.24: (a) Isosurface of structural sensitivity for the AS base flow at Re = 210 for the value S = 0.5, and corresponding sectional contour plots at (b) z = 0 and (c) y = 0.

Figure IV. 25 : 1 Figure

 251 Figure IV.25: Spectra consisting of the least stable direct modes at various Reynolds number for AS base flows. The oscillatory mode destabilizes between Re = 260 and Re = 280.

Figure IV. 27 :

 27 Figure IV.27: Spectra consisting of the least stable adjoint eigenvalues at various Reynolds number for AS base flows. The corresponding direct eigenvalues are also displayed with circles.

1 FigureFigure IV. 29 :

 129 Figure IV.28: 3-D views of oscillatory unstable direct mode at Re = 280, with (σ, St) = (0.0145, 0.1235) for AS base flow. (top) Isosurfaces of streamwise (a), radial (b) and azimuthal (c) perturbation velocities for the values u x = ±0.15, u r = ±0.1, u θ ± 0.1. (bottom) Corresponding sectional contour plot at (d) y = 0, (e) y = 0 and (f) z = 0, with solid-lines indicating iso-contours at u x = ±0.15, u r = ±0.1, u θ ± 0.1. (a)

2 × 10 3 3 × 10 3 4 × 10 3 6 × 10 3 η 1 Figure IV. 30 :

 333130 Figure IV.30: Evolution of optimal gain for the resolvent analysis of the sub-critical AS base flow at Re = 200. A peak is identified at St = 0.13.

1 Figure 1 Figure

 11 Figure IV.31: 3-D views of optimal forcing at Re = 200 calculated at St = 0.13. (top) Isosurfaces of streamwise (a), radial (b) and azimuthal (c) perturbation velocities for the values u x = ±5×10 -5 , u r = ±3.5 × 10 -5 , u θ ± 3.5 × 10 -5 . (bottom) Corresponding sectional contour plot at (d) y = 0, (e) y = 0 and (f) z = 0, with solid-lines indicating iso-contours at u x = ±5 × 10 -5 , u r = ±3.5 × 10 -5 , u θ ± 3.5 × 10 -5 .

Figure IV. 33 : 1 FigureFigure IV. 35 : 1 Figure

 331351 Figure IV.33: (a) Growth rate and (b) Strouhal number evolutions for direct modal analysis of PS base flow. Critical Reynolds number is Re 2 c = 280.(a)

Figure IV. 37 :

 37 Figure IV.37: (a) Isosurface of structural sensitivity for the PS base flow at Re = 280 for the value S = 0.4, and corresponding sectional contour plots at (b) z = 0 and (c) y = 0.

Fig. IV. 38 . 1 Figure IV. 38 : 1 Figure 1 Figure

 3813811 Figure IV.38: Evolution of optimal gain for the resolvent analysis of the sub-critical PS base flow at Re = 260. A peak is identified at St = 0.13.

Figure IV. 41 :

 41 Figure IV.41: Temporal evolution of lift coefficient at Re = 300. Calculation was initialized with PS base flow perturbed by direct unstable mode.

Figure IV. 43 :

 43 Figure IV.43: (a) Lift coefficient trend over the time where the dotted line represents σ NL = 0.1248. (b) Magnitude spectrum of the signal, where a pike at St NL = 0.063667 is identified.

Figure VI. 1 :

 1 Figure VI.1: Characteristic curve of the single channel flow for various tested cases, with SVRE for Short Vein Radial Equilibrium; EVRE for Extended Vein Radial Equilibrium; EVO for Extended Vein Outflow.

Figure VI. 2 :

 2 Figure VI.2: (a) Time evolution of a pressure signal recorded on three probes during URANS computation of the full-annulus compressor CME2. (b) Instantaneous entropy field after 12 revolutions.

Figure VI. 3 :

 3 Figure VI.3: Hypersonic BOLT configuration.

2 Fortran code A. 1 :r e a l * 8 : : s c a l p r o d r e a l * 8 :Fortran code A. 3 : 2 : 4 : 5 :

 21883245 L 2 scalar product subroutine s c a l P r o d L 2 ( n , vec1 , vec2 , s c a l P r o d ) i m p l i c i t none ! > INPUT v a r i a b l e s integer : : n r e a l * 8 , dimension ( n ) : : vec1 , vec2 r e a l * 8 , dimension ( n ) : : s c a l P r o d ! > LOCAL v a r i a b l e s integer : : i s c a l P r o d = 0 . d0 do i =1,n s c a l P r o d = s c a l P r o d + vec1 ( i ) * vec2 ( i ) end do end subroutine s c a l P r o d L Implementation of L 2 scalar product. subroutine scalprodL2_d ( n , vec1 , vec1d , vec2 , vec2d , s c a l p r o d , s c a l p r o d d ) i m p l i c i t none ! > INPUT v a r i a b l e s integer , : : n r e a l * 8 , dimension ( n ) : : vec1 , vec2 r e a l * 8 , dimension ( n ) : : vec1d , vec2d : s c a l p r o d d ! > LOCAL v a r i a b l e s integer : : i s c a l p r o d = 0 . d0 s c a l p r o d d = 0 . d0 do i =1,n s c a l p r o d d = s c a l p r o d d + vec2 ( i ) * vec1d ( i ) + vec1 ( i ) * vec2d ( i ) s c a l p r o d = s c a l p r o d + vec1 ( i ) * vec2 ( i ) end do end subroutine scalprodL2_d Fortran code A.2: Forward differentiation of L 2 scalar product. 123 subroutine scalprodL2_b ( n , vec1 , vec1b , vec2 , vec2b , s c a l p r o d , s c a l p r o d b ) i m p l i c i t none ! > INPUT v a r i a b l e s integer , : : n r e a l * 8 , dimension ( n ) : : vec1 , vec2 r e a l * 8 , dimension ( n ) : : vec1b , vec2b ) = vec1b ( i ) + vec2 ( i ) * s c a l p r o d b vec2b ( i ) = vec2b ( i ) + vec1 ( i ) * s c a l p r o d b end do s c a l p r o d b = 0 . d0 end subroutine scalprodL2_b Adjoint differentiation of L 2 scalar product. pd ( i c e l l , 2 ) = (wd( i c e l l , 2)-temp * wd( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 2 ) = temp temp = w( i c e l l , 3 ) /w( i c e l l , 1 ) pd ( i c e l l , 3 ) = (wd( i c e l l , 3)-temp * wd( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 3 ) = temp temp = w( i c e l l , 4 ) /w( i c e l l , 1 ) pd ( i c e l l , 4 ) = (wd( i c e l l , 4)-temp * wd( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 4 ) = temp temp = (w( i c e l l , 2 ) * w( i c e l l , 2)+w( i c e l l , 3 ) * w( i c e l l , 3)+w( i c e l l , 4) * & & w( i c e l l , 4 ) ) /w( i c e l l , 1 ) pd ( i c e l l , 5 ) = (gamma-1 . ) * (wd( i c e l l , 5) -0.5 * (2 * w( i c e l l , 2 ) * wd( i c e l l & & , 2)+2 * w( i c e l l , 3 ) * wd( i c e l l , 3)+2 * w( i c e l l , 4 ) * wd( i c e l l , 4)-temp * wd& & ( i c e l l , 1 ) ) /w( i c e l l , 1 ) ) p ( i c e l l , 5 ) = (gamma-1 . ) * (w( i c e l l , 5) -0.5 * temp ) Forward differentiation of primitive routine. subroutine computeprimitive_d_kform ( n c e l l , w, wdr , wdi , p , pdr , pdi , gamma) i m p l i c i t none ! > INPUT v a r i a b l e s integer : : n c e l l r e a l * 8 , dimension ( n c e l l , 5 ) : : w r e a l * 8 , dimension ( n c e l l , 5 ) : : wdr r e a l * 8 , dimension ( n c e l l , 5 ) : : wdi r e a l * 8 , dimension ( n c e l l , 5 ) : : p r e a l * 8 , dimension ( n c e l l , 5 ) : : pdr r e a l * 8 , dimension ( n c e l l , 5 ) : : p d i do i c e l l =1, n c e l l pdr ( i c e l l , 1 ) = wdr ( i c e l l , 1 ) p d i ( i c e l l , 1 ) = wdi ( i c e l l , 1 ) p ( i c e l l , 1 ) = w( i c e l l , 1 ) temp = w( i c e l l , 2 ) /w( i c e l l , 1 ) pdr ( i c e l l , 2 ) = ( wdr ( i c e l l , 2)-temp * wdr ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p d i ( i c e l l , 2 ) = ( wdi ( i c e l l , 2)-temp * wdi ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 2 ) = temp temp = w( i c e l l , 3 ) /w( i c e l l , 1 ) pdr ( i c e l l , 3 ) = ( wdr ( i c e l l , 3)-temp * wdr ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p d i ( i c e l l , 3 ) = ( wdi ( i c e l l , 3)-temp * wdi ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 3 ) = temp temp = w( i c e l l , 4 ) /w( i c e l l , 1 ) pdr ( i c e l l , 4 ) = ( wdr ( i c e l l , 4)-temp * wdr ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p d i ( i c e l l , 4 ) = ( wdi ( i c e l l , 4)-temp * wdi ( i c e l l , 1 ) ) /w( i c e l l , 1 ) p ( i c e l l , 4 ) = temp temp = (w( i c e l l , 2 ) * w( i c e l l , 2)+w( i c e l l , 3 ) * w( i c e l l , 3)+w( i c e l l , 4) * & & w( i c e l l , 4 ) ) /w( i c e l l , 1 ) pdr ( i c e l l , 5 ) = (gamma-1 . ) * ( wdr ( i c e l l , 5) -0.5 * (2 * w( i c e l l , 2 ) * wdr ( i c e l l & & , 2)+2 * w( i c e l l , 3 ) * wdr ( i c e l l , 3)+2 * w( i c e l l , 4 ) * wdr ( i c e l l , 4)-temp * wdr& ( i c e l l , 1 ) ) /w( i c e l l , 1 ) ) p d i ( i c e l l , 5 ) = (gamma-1 . ) * ( wdi ( i c e l l , 5) -0.5 * (2 * w( i c e l l , 2 ) * wdi ( i c e l l & & , 2)+2 * w( i c e l l , 3 ) * wdi ( i c e l l , 3)+2 * w( i c e l l , 4 ) * wdi ( i c e l l , 4)-temp * wdi& & ( i c e l l , 1 ) ) /w( i c e l l , 1 ) ) p ( i c e l l , 5 ) = (gamma-1 . ) * (w( i c e l l , 5) -0.5 * temp ) end do end subroutine computeprimitive_d_kform Fortran code B.3: K-formulated forward differentiation of primitive routine. subroutine computeprimitive_b ( n c e l l , w, wb , p , pb , gamma) i m p l i c i t none ! > INPUT v a r i a b l e s integer : : n c e l l r e a l * 8 , dimension ( n c e l l , 5 ) : : w r e a l * 8 , dimension ( n c e l l , 5 ) : : wb r e a l * 8 , dimension ( n c e l l , 5 ) : : p r e a l * 8 , dimension ( n c e l l , 5 ) : : pb r e a l * 8 : : gamma ! > LOCAL v a r i a b l e s integer : : i c e l l r e a l * 8 : : tempb wb = 0 . d0 do i c e l l =n c e l l ,1 , -1 wb( i c e l l , 5 ) = wb( i c e l l , 5 ) + (gamma-1.) * pb ( i c e l l , 5 ) tempb = -(0.5 * (gamma-1.) * pb ( i c e l l , 5 ) /w( i c e l l , 1 ) ) wb( i c e l l , 2 ) = wb( i c e l l , 2 ) + 2 * w( i c e l l , 2 ) * tempb wb( i c e l l , 3 ) = wb( i c e l l , 3 ) + 2 * w( i c e l l , 3 ) * tempb wb( i c e l l , 4 ) = wb( i c e l l , 4 ) + 2 * w( i c e l l , 4 ) * tempb wb( i c e l l , 1 ) = wb( i c e l l , 1 ) -(w( i c e l l , 2) * * 2+w( i c e l l , 3) * * 2+w( i c e l l & & , 4 ) * * 2 ) * tempb/w( i c e l l , 1 ) tempb = pb ( i c e l l , 4 ) /w( i c e l l , 1 ) wb( i c e l l , 4 ) = wb( i c e l l , 4 ) + tempb wb( i c e l l , 1 ) = wb( i c e l l , 1 ) -w( i c e l l , 4 ) * tempb/w( i c e l l , 1 ) tempb = pb ( i c e l l , 3 ) /w( i c e l l , 1 ) wb( i c e l l , 3 ) = wb( i c e l l , 3 ) + tempb wb( i c e l l , 1 ) = wb( i c e l l , 1 ) -w( i c e l l , 3 ) * tempb/w( i c e l l , 1 ) tempb = pb ( i c e l l , 2 ) /w( i c e l l , 1 ) wb( i c e l l , 2 ) = wb( i c e l l , 2 ) + tempb wb( i c e l l , 1 ) = wb( i c e l l , 1 ) + pb ( i c e l l , 1 ) -w( i c e l l , 2 ) * tempb/w(Adjoint differentiation of primitive routine. subroutine computeprimitive_b_kform ( n c e l l , w, wb , p , pb , gamma) i m p l i c i t none ! > INPUT v a r i a b l e s integer : : n c e l l r e a l * 8 , dimension ( n c e l l , 5 ) : : w r e a l * 8 , dimension ( n c e l l , 5 ) : : wbr r e a l * 8 , dimension ( n c e l l , 5 ) : : wbi r e a l * 8 , dimension ( n c e l l , 5 ) : : p r e a l * 8 , dimension ( n c e l l , 5 ) : : pbr r e a l * 8 , dimension ( n c e l l , 5 ) : : p b i do i c e l l =n c e l l ,1 , -1 wbr ( i c e l l , 5 ) = wbr ( i c e l l , 5 ) + (gamma-1.) * pbr ( i c e l l , 5 ) wbi ( i c e l l , 5 ) = wbi ( i c e l l , 5 ) + (gamma-1.) * p b i ( i c e l l , 5 ) tempbr = -(0.5 * (gamma-1.) * pbr ( i c e l l , 5 ) /w( i c e l l , 1 ) ) tempbi = -(0.5 * (gamma-1.) * p b i ( i c e l l , 5 ) /w( i c e l l , 1 ) ) wbr ( i c e l l , 2 ) = wbr ( i c e l l , 2 ) + 2 * w( i c e l l , 2 ) * tempbr wbi ( i c e l l , 2 ) = wbi ( i c e l l , 2 ) + 2 * w( i c e l l , 2 ) * tempbi wbr ( i c e l l , 3 ) = wbr ( i c e l l , 3 ) + 2 * w( i c e l l , 3 ) * tempbr wbi ( i c e l l , 3 ) = wbi ( i c e l l , 3 ) + 2 * w( i c e l l , 3 ) * tempbi wbr ( i c e l l , 4 ) = wbr ( i c e l l , 4 ) + 2 * w( i c e l l , 4 ) * tempbr wbi ( i c e l l , 4 ) = wbi ( i c e l l , 4 ) + 2 * w( i c e l l , 4 ) * tempbi wbr ( i c e l l , 1 ) = wbr ( i c e l l , 1 ) -(w( i c e l l , 2) * * 2+w( i c e l l , 3) * * 2+w( i c e l l & & , 4 ) * * 2 ) * tempbr /w( i c e l l , 1 ) wbi ( i c e l l , 1 ) = wbi ( i c e l l , 1 ) -(w( i c e l l , 2) * * 2+w( i c e l l , 3) * * 2+w( i c e l l & & , 4 ) * * 2 ) * tempbi /w( i c e l l , 1 ) tempbr = pbr ( i c e l l , 4 ) /w( i c e l l , 1 ) tempbi = p b i ( i c e l l , 4 ) /w( i c e l l , 1 ) wbr ( i c e l l , 4 ) = wbr ( i c e l l , 4 ) + tempbr wbi ( i c e l l , 4 ) = wbi ( i c e l l , 4 ) + tempbi wbr ( i c e l l , 1 ) = wbr ( i c e l l , 1 ) -w( i c e l l , 4 ) * tempbr /w( i c e l l , 1 ) wbi ( i c e l l , 1 ) = wbi ( i c e l l , 1 ) -w( i c e l l , 4 ) * tempbi /w( i c e l l , 1 ) tempbr = pbr ( i c e l l , 3 ) /w( i c e l l , 1 ) tempbi = p b i ( i c e l l , 3 ) /w( i c e l l , 1 ) wbr ( i c e l l , 3 ) = wbr ( i c e l l , 3 ) + tempbr wbi ( i c e l l , 3 ) = wbi ( i c e l l , 3 ) + tempbi wbr ( i c e l l , 1 ) = wbr ( i c e l l , 1 ) -w( i c e l l , 3 ) * tempbr /w( i c e l l , 1 ) wbi ( i c e l l , 1 ) = wbi ( i c e l l , 1 ) -w( i c e l l , 3 ) * tempbi /w( i c e l l , 1 ) tempbr = pbr ( i c e l l , 2 ) /w( i c e l l , 1 ) tempbi = p b i ( i c e l l , 2 ) /w( i c e l l , 1 ) wbr ( i c e l l , 2 ) = wbr ( i c e l l , 2 ) + tempbr wbi ( i c e l l , 2 ) = wbi ( i c e l l , 2 ) + tempbi wbr ( i c e l l , 1 ) = wbr ( i c e l l , 1 ) + pbr ( i c e l l , 1 ) -w( i c e l l , 2 ) * tempbr /w(& & i c e l l , 1 )wbi ( i c e l l , 1 ) = wbi ( i c e l l , 1 ) + p b i ( i c e l l , 1 ) -w( i c e l l , 2 ) * tempbi /w(K-formulated adjoint differentiation of primitive routine.

Table III . 3 :

 III3 Memory consumption and computation duration for various preconditioners on the laminar cylinder case at Re = 47.

	Preconditioner	ILU(1)	ILU(2)	ILU(3)	HLU
	Allocated memory for preconditioning	0.88 GB 0.96 GB 1.04 GB 1.61 GB
	Allocated memory for computation	1.14 GB 1.21 GB 1.29 GB 1.87 GB
	CPU time	46.5 s	36.84 s	37 s	26.2 s

Table III .

 III sI N

						K
	Allocated memory for preconditioning	1.04 GB 1.61 GB	1.72 GB	1.71 GB	1.84 GB
	Allocated memory for computation	1.29 GB 3.01 GB	3.03 GB	3.06 GB	3.05 GB
	CPU time	37 s	86.2 s	34 s	65.7 s	34.2 s

4: Memory consumption and computation duration for various operators on the laminar cylinder case at Re = 47.

Table III .

 III

5:

Comparison between SLEPc an elsA computation for direct global stability analysis, computation at Re = 47.

(σ L , St L ) (σ NL , St NL ) ∆St

					∆σ
	48 (0.0046, 0.1188) (0.0043, 0.1205) 1.3%	6.5%
	55	(0.03, 0.12)	(0.028, 0.13)	7.7%	7.1%
	60 (0.045, 0.1205)	(0.04, 0.1366)	11.8% 12.5%

Table III .

 III

8: Comparison between SLEPc an elsA computation for resolvent analysis computation at Re = 40 and St = 0.1165.

Table IV .

 IV

4: Comparison between SLEPc and elsA computations for adjoint global stability analysis at α = 4 • .

Table IV .

 IV

5: Comparison between SLEPc and elsA computations for resolvent analysis at α = 3.75 • and St = 0.1165.

Table IV . 8 :

 IV8 Indicative quantities for both fixed point computation methods at Re = 280 for the planar-symmetric base flow.

		LU-SGS	GMRES
	nonlinear iterations	10 4	41
	Final ∥ρ∥ 2	2.4 × 10 -8 9 × 10 -12
	Number of CPUs	300	500
	Memory/CPU after computation	0.72 GB	3.24 GB
	CPU time	300 s	10.65 h

Table IV . 10 :

 IV10 Quantitative figures from modal stability analysis on regular and fine mesh.

	Mesh	N DOF	N CPU Memory elsA Memory LU
	regular 47 × 10 6	1200	4.9 TB	∽ 300 TB
	fine	125 × 10 6 4788	20.5 TB	∽ 1000 TB

.10 compares the memory footprint of global stability computations carried out on the

δ kl being the Kroenecker symbol:δ kl = 1 if k = l, and δ kl = 0 if k ̸ = l

meaning being highly sensitive to numerical modeling.

rather than a direct eigensolver that returns the entire spectrum.

verifying J J T = J T J .

when R is not defined, iω is a marginally stable eigenvalue of J .

verifying WW H = I m .

a typical three levels of fill-in, resulting roughly in γ(3) = 3 and n KM = 120 consitute good parameters.

 ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆

 ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆

The contiguity is achieved when for a given v ∈ C N , each component of v ℜ and v ℑ are aligned in memory.

verifying X X H = I N .

we recall that Sp(C) ∈ R.

unlike a regular bifurcation in which the mode is steady.

Remerciements

Chapter II

Overview of global stability analysis tools and their implementation

 [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]and respresentative vizualisation of the mesh.

characteristics are set at the lateral domain boundaries, a uniform inflow velocity u ∞ is imposed at the x = 0 plane, and an outflow condition is imposed at the x = 20 plane.

As with the two previous test cases, the first step of this study is to compute a base flow, a steady solution of the conservative equations. The following section provides details on the strategy used to compute both AS and PS base flows.

IV.3.2 Fixed point computation

The computations were carried out on 18 Intel Xeon Broadwell E5-2680v4 nodes, utilizing 500 CPUs. For AS base flows, the calculations were converged using the backward Euler scheme combined with the GMRES method, initialized with a uniform field. The preconditioning used was FGMRES-ILU(3) with one RAS level [START_REF] Guilbert | Amélioration de l'inversion de grand systèmes creux pour la simulation numérique en mécanique des fluides[END_REF]. The internal and external Krylov bases were composed of 60 vectors, and the internal and external tolerances were set to 10 -3 and 10 -4 respectively.

Abstract

In specific physical conditions, rotating stall is an unsteady phenomenon that may be triggered in axial compressors of turbomachinery and could lead to drastic loads on the structures. Nowadays, its understanding remains partial and could be enhanced by performing a global stability analysis. After finding a fixed point of the governing equations of the flow, it consists of solving an eigenvalue problem where large sparse complex-valued linear systems need to be inverted. A usual strategy involves direct factorizations that, for systems exceeding 10 × 10 6 degrees-of-freedom, lead to memory overflow.

To overcome such a limit, the use of preconditioned iterative algorithms becomes mandatory. Furthermore, as complex arithmetic loops can't be vectorized, the complex systems need to be rewritten to preserve HPC optimizations, for which important work has been carried out recently at ONERA. In this thesis, a method to scale up the global stability problems to industrial cases has been developed: the KFGSM, for K-formulated Global Stability Analysis Method. The linearized matrix-vector products used to solve the eigenvalue problem are obtained by means of Algorithmic Differentiation, which ensures the preservation of loop optimizations. The complex problems have also been reformulated in the real-equivalent K formulation that preserves the sparsity patterns of the considered operators, allowing the reuse of the most advanced preconditioning techniques of compressible CFD.

The reformulated eigenvalue problem is solved by means of a Krylov-Schur algorithm where inversions are managed by the GMRES algorithm. The KFGSM is validated on three test cases with gradually increasing difficulty.

KEYWORDS: