{'fraction_non_alphanumeric': 0.02145689976610962, 'fraction_numerical': 0.008270906070980645, 'mean_word_length': 4.957592891760905, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A parallel and matrix free framework for global stability analysis of compressible flows

An numerical iterative framework for global modal stability analysis of compressible flows using a parallel environment is presented. The framework uses a matrix-free implementation to allow computations of large scale problems. Various methods are tested with regard to convergence acceleration of the framework. The methods consist of a spectral Cayley transformation used to select desired Eigenvalues from a large spectrum, an improved linear solver and a parallel block-Jacobi preconditioning scheme.


Introduction
In fluid dynamics, it is crucial for many technical applications to obtain stability information of the flow. While the stability of certain flows like the inviscid hyperbolic-tangent velocity profile is well studied [Blu70,Mic64], the global stability analysis of general flows is still a daunting task. Traditionally a stability analysis by e.g. ARPACK [LSY98] is performed by obtaining a linearized matrix representation of the discretized operator and using it to perform the analysis.
This kind of computation is usually limited by the available memory and computational resources. But certain physical phenomena require a highly resolved discretization and high order discretization schemes, for example the computation of acoustic modes generated by a compressible flow. This leads to an immense memory requirement during a global stability analysis, since it requires information about the describing operator in the whole computational domain. This encourages the use of a matrix free method to perform the stability analysis, which lessens the required memory significantly.
This article presents a framework which allows a matrix free computation of a global linear temporal stability analysis of compressible flows using a high order spatial discretization of the compressible Navier-Stokes equations. Compared to previous applications of this method, for example by [Mac09], the numerical stability analysis can be performed in a parallel manner.
We focus on the the temporal linear stability analysis of compressible flows. Being able to compute the stability analysis without the need of an operator in matrix formulation leads to a considerable reduction in required memory during a computation. In turn, certain challenges which arise because of the matrix free method have to be considered.
The presented framework employs iterative Krylov methods in order to obtain the desired Eigeninformation. Krylov methods are a natural choice for this kind of problem since they can be parallelized well and can operate without explicit knowledge of the operator. The base method of choice is the implicitly restarted Arnoldi method (IRAM) [Sor02]. It is combined with a spectral Cayley transformation to facilitate in the selection of the correct Eigenvalues and to provide a speedup of the computation [MSR94]. The use of this transform introduces a new linear system which has to be solved. An improvement of a common iterative method (GMRES) is tested as well as a preconditioning scheme in order to improve convergence.
Validation of the method is performed on a compressible mixing layer setup using a tangent hyperbolic profile as investigated by Blumen [Blu70]. The growth rate of the most unstable mode of this profile is computed using the proposed framework and compared to literature.
This article is structured as follows: Section 2 of this paper will discuss the physical model used in the stability analysis. Section 3 will explain in detail the concepts used in the stability analysis framework and the parallelization of it. To close, a model problem and results based on it are presented and evaluated in sections 4 and 5.

Governing equations
For the problems considered herein, the compressible Navier-Stokes (1) equations are used as the governing equations. They are given in a characteristic type formulation [Ses00] and are used in a matrix free fashion. Note that the use of this special formulation is not required for the framework, other formulations can be applied as well. The formulation introduces the primitive variables p as pressure, u i as the velocity vector and s as entropy. Additionally, τ ij describes the stress tensor, T the temperature, Φ the dissipation and ρ the density, with µ and µ v describing the viscosity and bulk viscosity of the fluid, respectively.
This set of equations is closed with Sutherland's law and the ideal gas equation and uses the summation convention with i = 1,2,3. Going forward, the short hand operator N (q) will be used to refer to the Navier Stokes equations, with q containing all primitive flow variables q = [p, u i , s].
Within the used scheme, the wave propagation terms are discretized with a sixth order finite difference scheme using compact finite difference schemes by [Lel92] on a staggered grid, which allows an accurate computation required to extract the stability information. The dissipation terms were discretized using a similar scheme of fifth order.

Ansatz and overview
In order to perform a global stability analysis, we need to formulate an appropriate ansatz which allows the computation of the temporal stability properties of the flow.
Therein J(q 0 ) describes the temporal evolution of the operator N lin = J(q 0 ), linearized around a base state q o . To this end we investigate the temporal stability. The goal of the stability analysis is to obtain the Eigenvalues λ of the Operator. For this, an ansatz fitting to the problem has to be chosen. The form of the ansatz depends on the structure of the investigated flow, especially on its degree of homogeneity.
For a general three dimensional flow without any degree of spatial homogeneity, the general approach (3) can be used for the stability analysis, as given by [Mac09].
q(x,y,z,t) = φ(x,y,z) exp(−iαct) Herein, α denotes the wave number of the temporal disturbance, c represents the wave propagation speed, while φ(x,y,z) describes the Eigenvector corresponding to α. Combining the model equation (2) with the general approach (3) results in an Eigenvalue problem (4) which has to be solved.
An Eigenvalue problem of this type can be solved directly if the Jacobian matrix J(q 0 ) is known. Since the memory requirements for storing this matrix are too high for global stability analysis problems, we have to employ matrix free methods to solve the it.

Matrix free framework
In order to compute the Eigenvalues needed for the stability analysis, the Jacobian of the operator solving the defining equations (1) is required. Since the defining equations are nonlinear and not available in a matrix formulation, an approximation of the Jacobian has to be computed. Therefore a Frechèt-type derivative (5) is employed for the linearization of the nonlinear operator during the stability analysis.
The correct choice of the parameter required for this derivative is critical to its accuracy. We used the findings on the influence of given by [SSS09] and set it accordingly to = q √ m , with m representing the machine accuracy of the system. The use of this linearization requires the evaluation of the defining equations N (q) every time the Jacobian is required during the computation. This application is usually one of the most numerically expensive operations and one should therefore try to minimize the number of required calls.
The use of a matrix free method to obtain information about J without explicitly storing it leads naturally to the use of iterative Krylov subspace methods, which do not require the Jacobian in explicit form but rather just matrix vector products of the type J(q 0 )φ(x,y). These matrix vector products can be computed in a matrix free fashion by using the previously described Frechèt derivative (5).

Krylov methods
Since the framework is designed to work in a matrix free environment, a method which can operate in such an environment is required for the stability analysis. Krylov methods suite this requirement well and will therefore be used as method of choice. The idea of Krylov methods is to iteratively create a sequence of subspaces K m which meets (6).
As stated before, N lin = J(q o ) is a (linearized) matrix representation of the Jacobian of N (q). Note that due to the iterative nature of Krylov methods only the product given by the Frechèt derivative (5) is required.
The basis V m of these subspaces can then be used to compute the Arnoldi relation (7), which projects a portion of N lin onto the Krylov subspace. This results in a Hessenberg matrix H m that is usually much smaller than N lin and contains an approximation of the eigenvalues and eigenvectors of the operator.

Implicitly restarted Arnoldi method
In order to solve the eigenvalue problem (4) iteratively, the implicitly restarted Arnoldi method (IRAM) developed by [Sor02] is used. It offers a robust method for obtaining Ritz pairs -approximations to the Eigenvalues and Eigenvectors of the operator N . Its main benefits are twofold. On one hand, it is a Krylov-type method and therefore does not need the explicit Jacobian of the operator. On the other hand, the IRAM allows a certain degree of control concerning which part of the spectrum is to be computed. This feature is especially important in order to obtain the desired results in a manageable time frame. In essence, it is a restarted Arnoldi method which allows the removal of unwanted Eigenvalue information during the restart of the iteration by using a QR-shift. A schematic overview of the algorithm is presented in algorithm 1. An important aspect of the IRAM with respect to the stability analysis is to define a suitable convergence criterion to determine when the desired Ritz pairs have been computed. As pointed out by [LS96], choosing an appropriate criterion is difficult for cases involving non-Hermitian operators which occur often when analyzing compressible flows since the residual norm f m of the Arnoldi decomposition does not necessarily imply convergence as opposed to the Hermitian case. To ensure a small backward error, a Ritz estimate (8) proposed by [LS96] is used. This criterion is applied to each computed Ritz pair (x, θ).

Cayley transformation
While the selective properties of the IRAM are very useful for this task, they are sometimes not sufficient when considering stability problems of compressible flows, which lead to problems due to clustering of Eigenvalues close to the origin or when dealing with badly conditioned flows. In order to enhance the selection properties of the framework even further, a Cayley transformation (9), as proposed by [MSR94], is applied to the Eigenvalue problem (4).
This leads to a generalized Eigenvalue problem of the type which has to be solved. The solution of this linear system requires an iterative solver. Due to the nature of the iterative solvers, the Cayley transformation is therefore only applied inexactly.
Algorithm 1 Implicitly restarted Arnoldi factorization [Sor02]. Parallel parts are represented by a "parallel" comment. Start with a m = k + p step Arnoldi factorization (algorithm 7) compute spectrum(H m ), select p shifts µ 1 , µ 2 , . . . ,µ p perform implicit QR-Shifts: convergence check using (8) parallel use obtained k-step factorization as starting point: To illustrate how the Cayley transformation affects the spectrum of the transformed operator, Fig. 2 shows a comparison of a sample spectrum with and without a Cayley transformation. Depending on the parameters µ and σ, the transformation has two effects on the shape of the spectrum. It introduces a rotation as well as a stretching effect. Therefore, it can be used to separate clusters of Eigenvalues and to rotate desired Eigenvalues into the part of the spectrum with positive Eigenvalues, which are favored by Arnoldi based methods. Another important feature of this transform is to map inner Eigenvalues to the outer part of the spectrum. This is important since Arnoldi based method usually converge first onto Eigenvalues distant from the origin.

Iterative solvers
The basic solver for the solution of (10) is a typical restarted GMRES based on [SS86], which was chosen for its robust numerical behaviour. Other methods are also applicable, for example BiCGstab.
Since the solution of (10) usually incurs a large numerical effort, it is beneficial to reduce the iterations required to solve it as much as possible. We investigate the possible speedup of a modified GMRES version and compared them to the standard formulation. The analyzed LGMRES uses error approximations to augment the Krylov subspace between restart cycles [BJM05]. The variant aims at improving the speed of convergence of the traditional restarted GMRES by lessening the impact of the necessary restart.

Parallelization
In order to compute large scale stability analysis, parallelization of the framework is required. Since the framework is based on IRAM, it lends itself well to parallelization since the numerically most expensive steps, the construction of the new Krylov subspaces, can be distributed among available CPUs such that each process does not require information about the entire computational domain. This reduces the required communication between the processes significantly. The steps that require parallelization are marked in the respective algorithm scheme such as Listing 1. The implementation used in the scope of this work achieves parallelization by parallelizing the derivatives using MPI, but other means of parallelization are applicable as well.

Preconditioning
Preconditioning techniques are often necessary to improve the speed of convergence when solving large linear systems arising during the stability analysis. While many methods for preconditioning are well researched, their application to a matrix free and parallel case is not straightforward. Since a global preconditioner is difficult to assemble due to the lack of the operator matrix, a local approach for the preconditioner was chosen for this work to serve as a proof of concept for matrix free preconditioning in parallel stability analysis. The framework uses a local block-Jacobi scheme similar to [HS91] which assembles a local preconditioner matrix on each process and uses it to accelerate the solution of the linear system. While this approach is rather simple, it serves as a starting point for future refinement.

Numerical framework overview
The proposed framework can be roughly divided into two parts, an outer iteration loop consisting of the IRAM-algorithm [Sor02], which obtains the stability information, and an inner iteration loop for the solution of the linear system required for the application of the inexact Cayley transformation. The relationship between these is depicted in Fig. 3.
An application of the operator N lin is required for each step of the inner iteration. Each of these calls of the operator translates to one evaluation of the defining equations (1). This is usually the step which requires the most computational resources. It is also the "innermost" step of the framework and therefore the most called.

Model problem
To illustrate the proposed framework, a model problem is used in the scope of this paper. We investigate the temporal stability of two-dimensional compressible mixing layer using a hyperbolic tangent velocity profile similar to the works of Blumen [Blu70]. The profile used is described by equation (11), which imposes a velocity profile in the x 1 -component of the velocity onto the stream. The thickness of the shear layer is controlled by the parameter δ.
with u ∞ as reference velocity and L 2 as domain length in x 2 direction. The model problem is set up in a computational domain with 128 × 512 points as depicted in Fig. 4. It employs a periodic boundary condition on the left and right boundaries as well as a non-reflective boundary condition on the lower and upper boundaries. Additionally, grid stretching was used to improve the resolution in the critical parts of the domain where the mixing effect occurs. This leads to a minimum cell distance in x 2 -direction of ≈ 0.08δ centered towards to highest velocity gradients present in the profile. The tangent hyperbolic velocity profile is seeded into the domain as a starting condition. In order to test the parallel aspect of the framework, the computation is distributed between at least 4 processors.
This model problem is investigated with regard to its linear temporal stability. The goal is to identify the least stable mode by using the previously outlined framework. The domain length in x 1 -direction L 1 should fit at least one entire wavelength of the desired mode, and was chosen to be exactly one wavelength of the least stable mode scaled by mixing layer thickness δ. Assuming a wavenumber α of the least stable mode, the length of the domain in x 1 -direction would result as L 1 2πδ α . The domain length in y-direction was chosen as 8 times L 1 .
Compared to the traditional matrix based approach on stability analysis, the presented matrix free framework requires significantly less memory. Consider the above model problem with a 128 × 512 points resolution. Assuming double precision, storing the operator alone would require (128·512) 2 ·64Bit ≈ 32GB of memory. When performing a parallel computation, this requirement is shared between all available nodes. The memory requirements of the matrix-free approach are tested in the following Sec. 6.

Validation
The proposed framework is validated using the model problem described in the previous section. The least stable modes for a various settings are computed and their growth rates compared to literature, see [Blu70,Mac09]. All computations are performed in parallel and cover the different types of iterative solvers described in the previous section as well as the preconditioner introduced previously.

Results -Validation
The proposed framework is able to compute growth rates of the least stable mode in good accordance to the results given by Blumen [Blu70] and Mack [Mac09]. The detailed results are compared in Tbl. 2. The different configurations are designed to cover a range of Mach numbers as well as different resolutions at a high Reynolds number. The different configuration parameters are listed in Tbl. 1. In addition to the computed growth rates, an exemplary mode computed during the validation is presented in Fig. 5. From this results we find the growth to be computed in good accordance with literature.

Numerical Experiments
As mentioned previously, the computational costs and required memory of the stability problem increases rapidly when considering problems which require high resolutions. While the memory requirements are mitigated to a great extent by the use of a matrix-free framework, the computational effort required is still significant. It is therefore important to reduce this computational effort, for example by improving the convergence behavior of the stability framework to avoid unnecessary iterations.  This section therefore investigates the effectiveness of the preconditioning method and iterative solvers with regard to speed up the iteration based on the results of numerical experiments. Since one of the main objectives of the presented framework is its matrix-free and parallel implementation, we will discuss the scalability and especially memory aspects based on the numerical experiments.
Lastly, while the choice of Cayley parameters has an influence on the performance of the framework, the presented numerical experiments do not cover a study on these parameters, since the focus of this work is to provide an overview of the framework. For a more detailed description of the effects of the Cayley transform, the reader is referred to [MSR94] or [Mac09]. Still, a few general remarks on the choice of Cayley parameters for this case are summarized at the end of this section.
The numerical experiments are based on the model problem discussed previously and differ in resolution, type of preconditioning and iterative solver used to cover a variety of possible combinations.
All further experiments share the following basic parameters: A Mach number M = 0.5, Reynolds number Re =→ ∞, domain length in x 1 -direction L 1 = 15.83δ, tolerance of the IRAM tol iram=1e−5 and tolerance for the GMRES-type method tol gmres = 1e − 6. The subspace size m of the IRAM was varied from 100 − 250, depending on the resolution of the case, with lower resolutions using smaller subspace sizes.

Preconditioning Scheme
Aside from the growth rate as a measure of the correctness of the computation, the number of required matrix-vector product operations is recorded. Since these correlate directly with the amount of required operator evaluations -the by far most numerically expensive part in these simulations. They serve as a good measure to compare different methods regarding speedup. Fig. 6 compares different computations, using the Block Jacobi preconditioner as well as running without preconditioning. A reduction in needed iterations for the preconditioned cases is clearly visible in the comparison. The reduction differs from case to case, but usually ranges 40-60%.
Aside from the speedup aspect another equally important aspect of the preconditioning became apparent during the numerical experiments. Since the stability problem arising from compressible flows is usually conditioned very poorly, computations without preconditioning often fail to produce accurate results. The use of a preconditioner alleviated this problem by enabling the computation of viable results in these cases and making the computation more robust. The resolutions used for the results above was therefore chosen small to allow a comparison to the unpreconditioned case.

Linear Solvers
Aside from the speedup achieved by the preconditioning, the speedup of the LGMRES solver is also investigated by numerical experiments. Table 3 has exemplary results for the effectiveness of the LGRMESmethod compared to standard GMRES.  In comparison to the effectiveness of the preconditioner, the speed up achieved by this method is rather small. Given the very low additional numerical effort incurred by this method compared to standard GMRES, the use of LGMRES is still advisable.

Memory and computation time considerations
In order to highlight the effects of the different methods with respect to speedup, a comparison between selected setups is displayed in Tbl. 3. The table compares the different configurations in terms of required iterations as well as required memory during the computation. Note that the computation using the higher resolution was not possible without the use of a preconditioner. The memory requirement displayed therein is estimated analytically. The  noted is the memory requirement of the method is heavily influenced by the method of preconditioning. In its current iteration, the use of the Block Jacobi method requires a significant amount of memory when investigating higher resolutions compared to the unpreconditioned case. This is a trade-off for the drastic reduction in required iterations compared to the unpreconditioned computations. For the examples provided in the given table, the preconditioned computation required roughly 20 times less iterations, a significant speedup even when considering the memory requirements of this method. However, the memory requirement is shared over the parallelized nodes.
As an additional metric, the scaling of the computational time and the memory requirement is collected in Fig. 7 to give an example for the scaling behavior of the framework when using the Block-Jacobi preconditioner. The scaling observed in this case massively deteriorates from 64 processors onward. This can be explained by the implementation of the preconditioner which requires a huge amount of calls of N (q) (ergo the right hand side of the considered equation) to compose the preconditioner prior to its application. With increasing number of processors used, the numerical costs of the composition surpass the costs of the solution of the preconditioned linear system, since the calls of N (q) require MPI calls. While this process is parallelized, it still acts as a bottleneck when using a relatively large number of processors -for the parameters given, this number is at 64 processors. Furthermore, the preconditioner cannot be omitted despite deteriorating the scaling of the framework since the preconditioner is a requirement to enable the correct computation in the first place.
Due to this behavior, the improvement of the preconditioner seems to be the most promising route to consider for improving the performance of the framework. As with the framework in general, the goal is to reduce the number of required right hand side calls. This could be achieved by a using non compact, lower order discretization of the operator N (q). This would significantly reduce the MPI communication during operator calls while also having the possible benefit of reducing the memory required during the composition of the preconditioner, [Mac09]. Fig. 8 shows the estimated memory requirement per node when using preconditioning. The actual memory requirements during numerical experiments the scale nearly linearly with the number of used processors, except for a small overhead which results in a sublinear scaling. This is likely dependant on the method of parallelization used in the framework, although this was not investigated in the scope of this work.

Cayley parameters and further remarks
The numerical experiments performed for this work generated further observations which are not covered in the above sections, yet can still be useful to get a better understanding of the problem at hand and to serve as guidelines for common problems.
Firstly, the condition of the linear systems deteriorates with increasing resolution. This is to be expected since an increase in resolution leads to a higher density of computable modes in the mixing layer. This results in an increased clustering of Eigenvalues which can affect the condition of the linear system. Additionally, a higher resolution can also result in an increase of unstable modes introduced by the discretization scheme. There are several ways to lessen the impact of these phenomena.
A first measure to be considered when dealing with the increased resolution is to increase the sizes of the subspaces available to the IRAM and GMRES. This step is usually required to account for the increased complexity of the computed spectrum and comes with an increase in memory and computational resources. For the resolutions investigated in this works, subspace sizes varying from 100 to 250 were found adequate.
In addition to the increase in subspace size, an adjustment of the Cayley parameters can also lead to an improvement of the computation for higher resolutions. Unfortunately, the method for determining the correct parameters is not as straightforward compared to the subspace size since the Cayley parameters are seemingly not correlating with the resolution in a simple fashion. For the resolutions considered herein, the following parameter sets in Tbl. 4 have been used successfully. While the Cayley parameters can influence the number of iterations required by the linear solver, the effect of fine tuning the parameters was not found significant. resolution recommended parameters 32 × 128 µ = 230.5, σ = 800.32 32 × 256 µ = 444.5, σ = 1200.32 64 × 512 µ = 444.5, σ = 1200.32 128 × 1024 µ = 444.5, σ = 1200.32

Conclusions
A framework for the matrix free and parallel stability analysis of compressible flows has been presented and validation results and speedup experiments have been showcased. To enhance the performance of the framework, a straightforward parallel preconditioning method was presented and tested successfully. The methods presented herein can be seen as a proof of concept for the matrix free stability analysis in parallel environments. Furthermore, the numerical experiments provide several ways that can be explored for a further increase in efficiency of the computation. As mentioned previously and judging from the results of the numerical experiments, the preconditioning scheme appears to be the most promising route to explore for further improvements. Two aspects can be considered as an angle for improvements to the preconditioning scheme.
Firstly, the construction of the preconditioner itself could be handled in a more sophisticated manner regarding memory requirements and method of construction. For example, rather than using the whole local operator as base for the preconditioning matrix, one could revert to using a reduced operator which still covers the physics but uses a much lower order finite difference scheme, as employed by [Mac09] in the non parallel case. This would reduce the memory required by the preconditioner while also facilitating its construction.
Secondly, the inversion and application of the preconditioning matrix offers room for improvement, e.g. by using sparse techniques for saving the matrix if the structure of the matrix is suitable.
While two different linear solvers were compared, the impact of the chosen method was found to be small compared to the impact of the preconditioning method. Therefore, while other methods could still be investigated with respect to speed up, the study of the preconditioning method seems to be more efficient, especially.