# A preconditioning for the spectral solution of incompressible variable-density flows

L. Reynier, Bastien Di Pierro, Frédéric Alizard, Anne Cadiou, Lionel Le Penven, Marc Buffat

► To cite this version:

L. Reynier, Bastien Di Pierro, Frédéric Alizard, Anne Cadiou, Lionel Le Penven, et al. A preconditioning for the spectral solution of incompressible variable-density flows. Computers and Fluids, In press, 266, pp.106024. 10.1016/j.compfluid.2023.106024. hal-04183956v2

HAL Id: hal-04183956

https://hal.science/hal-04183956v2

Submitted on 8 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A preconditioning for the spectral solution of incompressible variable-density flows

L. Reynier$^{a,*}$, B. Di Pierro$^a$, F. Alizard$^a$, A. Cadiou$^b$, L. Le Penven$^b$, M. Buffat$^a$

## Abstract

In the present study, the efficiency of preconditioners for solving linear systems associated with the discretized variable-density incompressible Navier-Stokes equations with semi-implicit second-order accuracy in time and spectral accuracy in space is investigated. The method, in which the inverse operator for the constant-density flow system acts as preconditioner, is implemented for three iterative solvers: the General Minimal Residual, the Conjugate Gradient and the Richardson Minimal Residual. We discuss the method, first, in the context of the one-dimensional flow case where a top-hat like profile for the density is used. Numerical evidence shows that the convergence is significantly improved due to the notable decrease in the condition number of the operators. Most importantly, we then validate the robustness and convergence properties of the method on two more realistic problems: the two-dimensional Rayleigh-Taylor instability problem and the three-dimensional variable-density swirling jet.

**Keywords:** Variable-density flows, incompressible Navier-Stokes equations, direct numerical simulation, preconditioning, elliptic solver

## 1. Introduction

Flows with large spatial density variations play an important role in widespread industrial and environmental applications such as fluid mixer found in pharmaceuticals processes or pollutant dispersion phenomenon.

Especially, the generation of baroclinic vorticity due to the interaction of non-parallel pressure and density gradients, and mass diffusion effects yield to a large variety of scale motions that presents a significant numerical challenge. In this context, accurate and fast numerical simulations of these types of flows require the development of efficient numerical solvers.

Since the prior work of Bell and Marcus [1], several authors proposed a second-order fractional time-step technique for solving incompressible flows with large density variations where Boussinesq approximation is no longer verified. For instance, Almgren

<footnote>$^a$Université Claude Bernard Lyon 1, LMFA, UMR5509, 69622 Villeurbanne, France</footnote>

<footnote>$^b$École Centrale de Lyon, CNRS, LMFA, UMR5509, 69130 Écully, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email address: loic@loicreynier.fr (L. Reynier)</footnote>
et al. [2] improve the convergence by using adaptive mesh refinement and Calgaro et al. [3] propose a hybrid method where the mass conservation is solved by a finite volume method and the velocity field is computed using a finite element method, just to name a few. The common feature of all the projection-like methods cited above is that at each time-step the pressure is determined by solving an elliptic equation which can be written in its most general form as

$$L_{\rho}(\phi) = \nabla \cdot (a\nabla\phi) = f \qquad (1)$$

due to the divergence-free constraint. Especially, $\phi$ is a scalar quantity related to the pressure, $a = 1/\rho$ with $\rho$ the approximation of the density for a given time and $f$ some right hand-side term that changes with time. Concus and Golub [4] first discuss of how to tackle the numerical solution of equation (1) with an iterative algorithm. For that purpose, the authors suggest solving iteratively the discrete counterpart of equation (1) where the system is rewritten as a Helmholtz-like equation using change of variable. The technique introduces a suitable choice of parameters to accelerate the convergence. In particular, the authors show that for a second-order finite difference scheme, the rate of convergence is independent of the mesh size for a smooth $a(x, t)$ function. However, the choice of an optimal set of parameters for a rapid convergence of equation (1) is directly connected to the condition number of the resulting operators. The high convergence rate is therefore not guaranteed for high-order schemes. Within the context of variable-density flows, Duffy et al. [5], Cook et al. [6], and more recently El Ouafa et al. [7], show that the discrete system associated with equation (1) is indeed ill-conditioned, which has for consequence of significantly slow down the performance of widely used iterative solvers. To overcome this difficulty, Guermond and Salgado [8] develop a numerical scheme where equation (1) is replaced by a Poisson equation using a penalty function to verify the incompressibility constraint. However, the proposed method requires the introduction of an additional term proportional to the divergence of the velocity fields onto the mass equation and its amplitude has to be fixed. More recently, Cook et al. [9] develop a hybrid Fourier spectral high-order compact finite-difference scheme to investigate variable-density flows which achieves tera-scalable computations on massively parallel machines [9]. Using a conservative variable formulation, they solve a Poisson equation for the pressure by introducing an estimation of the update velocity. As a consequence, their projection onto the divergence-free vector space is no longer exact. Finally, several authors suggest the use of preconditioning. As underlined by El Ouafa et al. [7], the incomplete LU (ILU) preconditioning is generally used. Nevertheless, this preconditioner does not scale well in parallel implementation with distributed memory. In an effort to enhance convergence properties, Duffy et al. [5] propose to combine a multigrid technique with a preconditioning. While the method improves the classical multigrid projection used for example by Ravier et al. [10], the efficiency of the preconditioning is not discussed within the context of spectral methods.

The problem of the variable-coefficient discrete Poisson equation is also considered by Knikker [11] for low-Mach-number flows for which an extension of the projection-type method proposed by Bell and Marcus [1] is used. Focusing on high-order finite difference schemes, extensive numerical experiments are then carried out by Knikker [11] to illustrate the performance of various algorithms such as Conjugate Gradient (CG)-like methods with different preconditioning techniques. Besides, the problem of ill-conditioned matrices associated with spectral discretization of the Helmholtz equa-
tion is discussed by Haldenwang et al. [12]. The authors present an iterative algorithm well-suited or Chebyshev polynomial approximation which uses a preconditioning built on a second-order difference discretization of the Laplacian operator. The more general case corresponding to equation (1) is not considered in the work of Haldenwang et al. [12].

From the above discussion, it appears that the evaluation of an efficient preconditioner combined with widely used iterative linear solvers for the resolution of the variable-density incompressible Navier-Stokes equations (VDINSE) with spectral spatial accuracy and second-order time-accurate scheme has not yet been fully addressed. Motivated by this question, it is the objective of the present work to present a preconditioning technique and its performance along with some iterative solvers for solving the VDINSE system.

This paper presents a numerical method for solving VDINSE that includes the motion induced by Fick's mass diffusion law, with spectral spatial accuracy, using a semi-implicit method for the viscous and diffusive term and a second-order fractional time step to hold incompressibility. It extends the previous work of Di Pierro and Abid [13] where the influence of preconditioning is not addressed. In a first section, after having briefly presented the system of equations and numerical schemes, the preconditioning technique based on the constant-density operator is introduced for both the velocity and pressure equations. In the second section, the cost of the preconditioner is investigated through numerical test cases carried out on the implicit systems associated with both the velocity and pressure equations for various iterative solvers. The third section highlights the robustness of the methods by time-marching the VDINSE system for representative numerical flow cases.

## 2. Governing equations and numerical schemes

### 2.1. Mathematical model

We consider hereafter the motion of a viscous fluid in an inhomogeneous medium that takes place in a bounded rectangular domain $\Omega$ (with boundaries noted $\partial\Omega$) where the Cartesian coordinate system is defined by the $x, y, z$ axes and in a time interval $t \in [0, T]$. The mathematical model for the variable-density incompressible Navier-Stokes equations (VDINSE) used in the present contribution is detailed by Frank-Kamenetskii [14], Kazhikhov and Smagulov [15], Antontsev et al. [16] and Guillén-González et al. [17]. As in previous cited studies, we introduce the mean density $\rho(x, t)$ and the mean-volume velocity $\mathbf{u} = (u, v, w)^T(x, t)$, then the dimensionless equations of motion read

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + \frac{\nabla p}{\rho} = \zeta(\rho, \mathbf{u}) + \mathbf{f}, \qquad (2a) $$

$$ \frac{d\rho}{dt} + \mathbf{u} \cdot \nabla \rho = \frac{1}{Re Sc} \nabla^2 \rho, \qquad (2b) $$

$$ \nabla \cdot \mathbf{u} = 0, \qquad (2c) $$

$$ \zeta(\rho, \mathbf{u}) = \frac{1}{\rho Re} \nabla^2 \mathbf{u} + \frac{1}{\rho Re Sc} (\mathbf{u} \cdot \nabla \nabla \rho + (\nabla \rho \cdot \nabla) \mathbf{u}). \qquad (2d) $$

Here, the mass diffusion is modeled according to the Fick's diffusion law. In the momentum equation (2d), $p$ is a potential function analogous to the pressure and $f$ represents
an external body force. In equation (2a) and equation (2b), Re and Sc are the Reynolds and Schmidt numbers, respectively, defined as

$$Re = \frac{\bar{\rho}UL}{\mu}, \quad Sc = \frac{\mu}{\bar{\rho}\lambda}, \qquad (3)$$

with $\mu$ and $\lambda$ the dynamic viscosity and mass diffusivity of the fluid respectively, and $\bar{\rho}$, $U$, $L$ being characteristic density, velocity and length scales.

The system ((2)) is then closed with boundary conditions. Two different types of boundary conditions are used in the present study: no-slip and periodic boundary conditions. For wall-bounded flows, the same boundary conditions as in Guillén-González et al. [17] are imposed:

$$u_{|\partial\Omega} = 0, \quad \frac{\partial\rho}{\partial n_{|\partial\Omega}} = 0. \qquad (4)$$

Finally, one may note that while many fundamental results have been found with periodic boundary conditions [18, 19, 20], such boundary conditions could also be used for spatially developing flows by introducing a fringe region method [21].

The system (2) is solved using a fractional-step scheme with second-order accuracy. The viscous and diffusive terms are discretized by using a semi-implicit Crank–Nicholson scheme (by extension of the prior study of Bell and Marcus [1]):

$$\frac{\boldsymbol{u}^{n+1} - \boldsymbol{u}^n}{\delta t} = -[\boldsymbol{u} \cdot \nabla \boldsymbol{u}]^{n+1/2} - \left[\frac{\nabla p}{\rho}\right]^{n+1/2} + \frac{1}{2} \left(\zeta(\rho^{n+1/2}, \boldsymbol{u}^{n+1}) + \zeta(\rho^{n+1/2}, \boldsymbol{u}^n)\right), \qquad (5a)$$

$$\frac{\rho^{n+1} - \rho^n}{\delta t} = -[\boldsymbol{u} \cdot \nabla \rho]^{n+1/2} + \frac{1}{2} \frac{1}{\mathrm{Re} \, \mathrm{Sc}} (\nabla^2(\rho^{n+1}) + \nabla^2(\rho^n)), \qquad (5b)$$

$$\nabla \cdot \boldsymbol{u}^{n+1} = 0, \qquad (5c)$$

where $\delta t$ is the time step, the $n$ superscript denotes the solution at time $t^n = n\delta t$ and $n+1/2$ represents a second-order approximation at time $t^{n+1/2}$. Here, $[\boldsymbol{u} \cdot \nabla \boldsymbol{u}]^{n+1/2}$ and $[\boldsymbol{u} \cdot \nabla \rho]^{n+1/2}$ are estimated through an Adams–Bashforth scheme whereas the viscous/mass diffusion term $\zeta$ is computed with $\rho^{n+1/2} = (\rho^{n+1} + \rho^n)/2$ for stability considerations [22].

Following Bell and Marcus [1], equations (5a) and (5b) are time-integrated by introducing an intermediate velocity $\boldsymbol{u}^*$:

$$\frac{\boldsymbol{u}^* - \boldsymbol{u}^n}{\delta t} = -[\boldsymbol{u} \cdot \nabla \boldsymbol{u}]^{n+1/2} - \frac{\nabla p^{n-1/2}}{\rho^{n+1/2}} + \frac{1}{2} \left( \zeta(\rho^{n+1/2}, \boldsymbol{u}^*) + \zeta(\rho^{n+1/2}, \boldsymbol{u}^n) \right), \qquad (6a)$$

$$\frac{\rho^{n+1} - \rho^n}{\delta t} = -[\boldsymbol{u} \cdot \nabla \rho]^{n+1/2} + \frac{1}{2} \frac{1}{\mathrm{Re} \, \mathrm{Sc}} (\nabla^2(\rho^{n+1}) + \nabla^2(\rho^n)). \qquad (6b)$$

The velocity $\boldsymbol{u}^{n+1}$ is then updated by performing the projection of $\boldsymbol{u}^*$ onto a divergence-free subspace

$$\boldsymbol{u}^{n+1} = \boldsymbol{u}^* - \delta t \frac{\nabla \phi}{\rho^{n+1/2}} = P_{\rho}(\boldsymbol{u}^*), \qquad (7)$$
where the scalar $\phi$ is solution of the equation

$$L_{\rho}(\phi) = \nabla \cdot \left( \frac{1}{\rho^{n+1/2}} \nabla \phi \right) = \frac{1}{\delta t} \nabla \cdot \mathbf{u}^* \quad (8)$$

with the same boundary condition used for the density field (periodic or $(\partial\phi/\partial\mathbf{n})|_{\partial\Omega} = 0$) to preserve precision order.

At this point, the choice of the update pressure $p^{n+1/2}$ is crucial to truly verify a second-order accurate time scheme as shown by Brown et al. [23] in the case of constant-density flows. By injecting $\mathbf{u}^*$ from equation (7) into equation (6a) and by comparing with the requested scheme from equation (5a) one can get the update pressure $p^{n+1/2}$ in gradient form

$$\nabla p^{n+1/2} = \nabla p^{n-1/2} + \nabla \phi - \frac{\delta t}{2} \rho^{n+1/2} \zeta(\rho^{n+1/2}, \nabla \phi / \rho^{n+1/2}) \quad (9)$$

or in (equivalent) scalar form

$$p^{n+1/2} = p^{n-1/2} + \phi - (\nabla^2)^{-1} \left( \nabla \cdot \left( \frac{\delta t}{2} \rho^{n+1/2} \zeta(\rho^{n+1/2}, \nabla \phi / \rho^{n+1/2}) \right) \right), \quad (10)$$

with $(\nabla^2)^{-1}$ the inverse of Laplacian operator. As mentioned by Di Pierro and Abid [13], the projection operator $P_{\rho}$ is an exact projection operator. Indeed,

$$P_{\rho}(\mathbf{u}) = \mathbf{u} - \frac{1}{\rho} \nabla L_{\rho}^{-1} \nabla \cdot \mathbf{u} \quad (11)$$

and then, for any divergence-free vector field $v$ and any vector field $w$ one gets:

$$P_{\rho}(v) = v, \quad (P_{\rho} \circ P_{\rho})(w) = P_{\rho}(w). \quad (12)$$

### 3. Elliptic solver

The solution of an elliptic equation with constant coefficients using spectral approximation can be computed with fast and accurate algorithms which exploit properties of the spectral basis [24, 12, 25]. The main difficulty for solving VDINSE is that the pressure equation (8) is far more complicated than just a standard Poisson equation. Especially, such an elliptic problem discretized with spectral methods yields to very ill-conditioned full matrices for which direct solvers are inefficient. As underlined by Canuto and Quarteroni [26] and Peyret [25], systems such as equation (8) have to be solved with an iterative process. Then, the choice of an effective preconditioner is crucial for increasing convergence speed and improving the accuracy of iterative solvers [27]. Since the density variations considered here are sufficiently smooth to be accurately projected onto a spectral basis, this paper proposes to use the inverse of the operator associated with constant-density flow as preconditioner for three-dimensional Direct Numerical Simulation (DNS). Within a Galerkin formulation, the inverse matrix can be easily computed. Indeed, constant-coefficients elliptic problems lead to diagonal systems using Fourier decomposition and to quasi-tridiagonal operators with Chebyshev polynomials as basis functions. Hence, for three-dimensional DNSs, these elliptic problems
are reduced to diagonal systems with Fourier-Fourier-Fourier (F-F-F) decomposition or $N^2$ quasi-tridiagonal systems for Chebyshev-Fourier-Fourier (C-F-F) decomposition. As presented by Di Pierro and Abid [28] and Alizard et al. [29], such decompositions are sufficient to study many fundamental configurations. The case of Chebyshev-Chebyshev-Fourier and fully Chebyshev decompositions are not treated here because they lead to very large pentadiagonal or heptadiagonal matrices. For the C-F-F case, the solution of the preconditioned system is determined by using the algorithm of Thual [30]. More specifically, the solution of the tridiagonal system if computed using the recurrence relation detailed in Peyret [25, appendix B]. The algorithm requires $O(N)$ operations, without storing the operator used for the preconditioning. As mentioned by [25], the algorithm is stable is the solution remains nicely bounded. For that purpose, sufficient conditions to ensure such a property are given by Peyret [25], which are satisfied for the following cases. The efficiency of such preconditioning is discussed in the next section.

### 3.1. Implicit viscous solver

For illustration purposes, we focus in this section on the one-dimensional case where quantities are expressed using Chebyshev collocation discretization. For that purpose, the viscous linear operator associated with the discrete equation (6a) can be rewritten in the form

$$V_{\rho}(\mathbf{u}) = \left( \mathbf{I} - \frac{a}{\rho} \left( \mathbf{D}^2 + \frac{\partial^2 \rho}{\partial x^2} + \frac{\partial \rho}{\partial x} \mathbf{D}^1 \right) \right) \mathbf{u} \quad (13)$$

where $\mathbf{D}^k$ is the Chebyshev differentiation matrix of order $k$, $\mathbf{I}$ the identity matrix, and $a = \frac{1}{2}\delta t / Re$ by assuming here $Sc = 1$ for simplicity's sake. We recall that an efficient resolution of $V_{\rho}(\mathbf{u}) = \mathbf{b}$ is strongly correlated with the condition number of the operator $V_{\rho}$. We introduce the constant-density counterpart of $V_{\rho}$

$$V = (\mathbf{I} - a\mathbf{D}^2) \quad (14)$$

which is used for preconditioning $V_{\rho}$. For the numerical tests presented below, a classical top-hat like profile

$$\rho(x) = 1 + \frac{s-1}{2} \left( \tanh \left( \frac{x+x_0}{d} \right) - \tanh \left( \frac{x-x_0}{d} \right) \right) \quad (15)$$

is used for the density field where $d = 0.1$ is the gradient length scale, $s$ is the density ratio, and $x_0$ is chosen as one fifth of the domain length.

Figures 1 to 3 show the condition number $\kappa$ — computed by a singular values decomposition — of $V_{\rho}$ and $V$ when increasing the number of collocation points $N$, the density ratio $s$, and the viscous coefficient $a$. The changes of $\kappa$ with respect to the variations in $N$, $s$ and $a$ for the left- and right-preconditioning ($V^{-1}V_{\rho}$ and $V_{\rho}V^{-1}$, respectively) are also illustrated. The same behavior is obtained with a Fourier expansion, not detailed here for the sake of conciseness. It is clear from the figures that even if $\kappa(V_{\rho})$ and $\kappa(V)$ reach high values ($O(10^4) - O(10^8)$), the proposed preconditioning technique greatly reduces the condition number. Moreover, those last two are quasi-independent of the collocation point number $N$ and remains small ($< O(10^2)$) even for larger values of $a$ (i.e. for viscous dominated flows) or $s$ (i.e. strongly stratified flows).

Let us now focus on the effect of such a preconditioning technique onto the convergence speed of the algorithms for the resolution of $V_{\rho}(\mathbf{u}) = \mathbf{b}$. To this end, three iterative algorithms — based on the minimization of the residual $\|\mathbf{b} - V_{\rho}(\mathbf{u})\|_2^2$ — are investigated:
<figure><img src="image_1.png" /><figcaption>Figure 1: Condition number of the viscous operators versus the Chebyshev collocation points number $N$ with $s = 2$ and $a = 0.01$.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: Condition number of the viscous operators versus the density ratio $s$ with $a = 0.01$ and $N = 256$ Chebyshev collocation points.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Condition number of the viscous operators versus $a$ with $s = 2$ and $N = 256$ Chebyshev collocation points.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: Convergence history for the five preconditioned iterative solvers to solve $V_p(u) = b$ with $a = 10^{-3}$ and $s = 2$ with Fourier differentiation (left side) and Chebyshev differentiation (right side) and $N = 256$ collocation points. Convergence history for different parameters can be found in Appendix A.</figcaption></figure>

• Richardson Minimal Residual (RMR) [26],

• Conjugate Gradient (CG) [27],

• General Minimal Residual (GMRES) [31].

One may note that the Biconjugate Gradient (BiCG) and stabilized Biconjugate Gradient (BiCGStab) methods are not proposed here, since for DNSs, matrices are never formed explicitly. As a consequence, computing the operator transposition used in BiCG method yields some difficulties. The convergence history of such algorithms — assuming a zero function as initial guess for any iterative process — is shown in figure 4 for $s = 2$ and $a = 10^{-2}$, as an example. Similar investigations over a wide range of $(a, s)$ values are reported in Appendix A. For these cases, the convergence history of BiCG and BiCGStab are presented as well for readers interested in other applications. The independence of the convergence rate with the number of collocation points $N$ has also been verified and is not detailed here for the sake of conciseness.

Results are summarized below. The CG method provides the slower convergence and do not reach a residual smaller than $10^{-6}$ (for all numerical tests performed here): it turns out not to be competitive with other methods in terms of accuracy and computational time. This behavior is probably due to the non-symmetric properties of the operator associated with equation (13) For $s = 2$, RMR and GMRES methods converge both towards a residual value of $\approx 10^{-12}$ in about 20 iterations independently of $a$ (see figure B.14 in Appendix A). For $a = 10^{-2}$, it can be seen in figure B.15 (see Appendix A) that the RMR does not converge when $s > 10$. For this value of $a$, the GMRES algorithm is the only method that reaches a residual value smaller than $10^{-9}$ for extremely stratified flows ($s > 100$). One may note that the convergence of the GMRES algorithm
<table><caption>Table 1: Computational time in seconds of the preconditioned RMR and GMRES to solve $V_ρ(\mathbf{u}) = \mathbf{b}$, with $s = 1000$ and $N = 256$ Chebyshev basis collocation points. Computations are performed with an Intel Xeon E5-2670 CPU mounted on a Dell PowerEdge C8220 Compute Node.</caption><thead><tr><td>Number of iterations</td><td>RMR</td><td>GMRES</td></tr></thead><tbody><tr><td>5</td><td>5.6 10<sup>-3</sup></td><td>1.0 10<sup>-2</sup></td></tr><tr><td>10</td><td>1.1 10<sup>-2</sup></td><td>3.4 10<sup>-2</sup></td></tr><tr><td>20</td><td>2.0 10<sup>-2</sup></td><td>9.7 10<sup>-2</sup></td></tr><tr><td>40</td><td>4.2 10<sup>-2</sup></td><td>2.2 10<sup>-1</sup></td></tr></tbody></table>

is here sufficiently fast to not require a restart process. Let $M$ the number of iterations necessary to reach a convergence criterion, the complexity of RMR method is $O(N \times M)$ whereas the GMRES complexity is $O(N \times M^2)$, due to the additional orthonormalization process. Hence, the latter is expected to be much slower even for a few iterations (especially if the numerical method is parallelized with a domain decomposition strategy). This can be seen in table 1 which shows the computational time of RMR and GMRES methods without an explicit construction of differentiation matrices.

Finally, to accelerate the convergence, the previously converged solution $\mathbf{u}^n$ is used as an initial guess for solving $\mathbf{u}^{n+1}$. This leads to a fast convergence for both the RMR and GMRES methods: only 3 to 5 iterations are required to obtain a residual value varying from $O(10^{-12})$ to $O(10^{-14})$. Hence, the RMR method should be adopted for solving equation (6a) when the density ratio is lower than $s = 10$ and the GMRES method should be adopted for steeper cases.

### 3.2. Pressure equation

As previously mentioned, equation (8) for the scalar $\phi$ is very ill-conditioned and has to be solved accurately for ensuring mass conservation. The scalar $\Phi = \phi/\sqrt{\rho}$ is introduced in order to rewrite equation (8) as a modified Helmholtz equation

$$K_{\rho}(\Phi) = \nabla^2 \Phi + \left( \frac{\nabla^2 \rho}{2\rho} - \frac{3}{4} \frac{\|\nabla \rho\|^2}{\rho^2} \right) \Phi = \sqrt{\rho} f \quad (16)$$

with $f = (\nabla \cdot \mathbf{u}^*) / \delta t$. One may recall that Concus and Golub [4] use a similar technique for the general case of elliptic equations with variable coefficients. It's worth noting that when using fully periodic boundary conditions, $\Phi$ remains periodic. When wall boundary conditions

$$\frac{\partial \rho}{\partial \mathbf{n}_{|\partial \Omega}} = 0, \quad \frac{\partial \phi}{\partial \mathbf{n}_{|\partial \Omega}} = 0, \qquad (17)$$

are considered, a Neumann condition is also obtained for $\Phi$ and reads

$$\frac{\partial \Phi}{\partial \mathbf{n}_{|\partial \Omega}} = 0. \qquad (18)$$

As for the viscous operator, the operator $K_\rho$ is solved iteratively and is preconditioned by the inverse of the discrete Laplacian operator $(\nabla^2)^{-1}$ which is nothing else than the inverse of $K_\rho$ with constant density. Following the same line of thought as in the previous section, the behavior of the preconditioned operator is firstly characterized through a one-dimensional test case discretized with Chebyshev collocation matrices [25].
Using the density profile described by equation (15), the condition numbers of these implicit operators (preconditioned or not) are shown in figures 5 and 6 versus either the collocation point number or the density ratio. One may notice that $K_ρ$ and $L_ρ$ are both ill-conditioned even if the use of $K_ρ$ instead of $L_ρ$ leads to some improvements. The figures show that both left- and right-preconditioning provide a condition number $κ = O(10^0-10^2)$. Especially, the condition number is seen independent of the collocation point number $N$. In figure 6, one may also observe that $κ((∇²)^{-1} K_ρ) < κ((∇²)^{-1} L_ρ)$ for $0.1 ≤ s ≤ 100$.

The same algorithms as in the previous section (CG, RMR and GMRES) are studied and compared with the fixed point (FP) algorithm proposed by Di Pierro and Abid [13]. The convergence history of these iterative solvers is shown in figure 7 for $s = 2$ and $a = 10^{-2}$ and in Appendix B for different values of the density ratio $s$. The independence of the convergence rate with the number of collocation points $N$ has been verified and is not detailed here for the sake of conciseness. Results are now summarized. As for the viscous operator, the CG method is inaccurate in comparison to other methods, especially, the residual exhibits a plateau around 0.5 for all test cases. The GMRES has the fastest convergence, leading to a decay of the residual by 12 decades in about 10 iterations for all test cases. Besides, the RMR needs 4 times more iterations in order to reach the same level of convergence as the one obtained using the GMRES for $s = 2$ and does not converge for higher density ratios. Finally, the FP method [13] converges faster than the RMR method for $s = 2$ and has the same convergence rate as the GMRES for more stratified cases while being less accurate (its residual saturates around $O(10^{-9})$). In that regard, the RMR method should be adopted for solving equation (16) when the density ratio is around $s = 2$ while the GMRES should be used for higher density ratios.

# 4. Results and validation

## 4.1. Time accuracy and mass conservation

The derivation of an exact solution of the VDINSE system is not straightforward. To the best of the authors' knowledge, none has been found until now. Hence, to estimate the accuracy of the presented method, mass source and body-force terms are added to the VDINSE equations such that

$$u_e(x, y, z, t) = \cos(t) \sin(x) \cos(y) \cos(z) \qquad (19a)$$

$$v_e(x, y, z, t) = \cos(t) \cos(x) \sin(y) \cos(z) \qquad (19b)$$

$$w_e(x, y, z, t) = -2 \cos(t) \cos(x) \cos(y) \sin(z) \qquad (19c)$$

$$p_e(x, y, z, t) = \sin(t) \sin(2x) \sin(2y) \sin(2z) \qquad (19d)$$

$$\rho_e(x, y, z, t) = 1 + \frac{s-1}{2} (1 + \sin(t) \sin(x) \sin(y) \sin(z)) \qquad (19e)$$

defines an exact solution of the modified VDINSE, where $s$ is (an arbitrary positive) density ratio. One may note that this corresponds to a Taylor-Green vortex with an imposed pressure and density field.

Let $(u(T), p(T), \rho(T))$ be the numerical solution extracted at a time $t = T$ computed by marching forward in time the equations (6a), (6b), (7) and (10) with a time step $\delta t$ and an initial condition matching the system (19) for $t = 0$. Then, the temporal errors at $t = T$ are defined with respect to the exact solution as $\epsilon_u = ||u(T) - u_e(T)||_2$,
<figure><img src="image_5.png" /><figcaption>Figure 5: Condition number of the pressure operators versus Chebyshev collocation point number *N* with *s* = 2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Condition number of the pressure operators versus density ratio *s* with *N* = 256 Chebyshev collocation points.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Convergence history of the different preconditioned iterative solvers to solve the modified pressure equation (16) with *s* = 2, and *N* = 256 Chebyshev collocation points. Convergence history for different parameters can be found in Appendix B.</figcaption></figure>
$\epsilon_p = ||p(T) - p_e(T)||_2$ and $\epsilon_\rho = ||\rho(T) - \rho_e(T)||_2$. The corresponding errors are shown in figure 8 for $128 \times 128 \times 128$ grid discretization points and with $Re = 1$, $Sc = 1$, $s = 10$ and $T = \pi/4$. Figure 8 shows that the numerical solution is, as expected, globally second-order accurate in time. Additionally, the spectral accuracy has been verified by running the same simulation with different grid sizes. Machine precision ($10^{-14}$) is reached from a $32^3$ mesh for the divergence of the velocity field.

From a more physical point of view, the total mass of the considered system $M = \int_{\Omega} \rho \, d\Omega$ has to be strictly conserved throughout the simulation. From Fick's law and equation (2b), the mass evolution over time is directly given by the divergence of velocity field

$$ \frac{dM}{dt} = \int_{\Omega} \nabla \cdot \mathbf{u} \, d\Omega. \qquad (20) $$

Figure 9 shows the time evolution of the $L^2$ norm of $\nabla \cdot \mathbf{u}$ for the previous test case. One can observe that this quantity is nearly equal to the machine precision all along the simulation.

## 4.2. 3D solver performance

The performance of the presented method is now discussed by running three-dimensional DNSs associated with the Taylor-Green vortex discussed above (see system (19)). The computational times for the pressure iterative solver are shown in table 2 for the different solvers investigated in the previous section (FP, RMR and GMRES) and with $128 \times 128 \times 128$ discretization points. Unlike the 1D case detailed in the previous section, a Fourier (periodic) decomposition is used. For the preconditioning, a singular value decomposition is performed to compute a pseudo-inverse of the discrete Laplacian. It is equivalent to remove the first spectral coefficient associated with the polynomial of degree one. It imposes that $\phi$ and $\Phi$ are zero-mean.

One can see from table 2 that RMR and GMRES solvers are faster than the FP method. This is due to the problem formulation: $K_\rho$ needs only two physical/spectral transforms whereas four are needed to evaluate $L_\rho$. To accelerate the convergence, the solution $\phi^{n-1/2}$ is used as an initial condition for evaluating $\phi$ at the time-step $n + 1/2$. In agreement with the condition number of the preconditioned operators, the performance of each solver is also independent of the density ratio $s$. Hence, it validates the performance of the method even for large density ratios. Moreover, the computational time per iteration per thread and per cell is in the order of magnitude of the CPUs clock time. Hereafter, the GMRES algorithm is adopted for the DNS test cases for illustration purposes.

## 4.3. Rayleigh-Taylor instability

First, the effectiveness and robustness of the present method is illustrated through a classical configuration of a 2D Rayleigh-Taylor instability (RTI): a heaver fluid is maintained above a lighter one and gravity acceleration is in the opposite direction of the density gradient. No-slip boundary conditions are used on the top and bottom walls, while periodic boundary conditions are used on the left and right boundaries. The initial density profile is a smoothed Heavyside function in the vertical direction disturbed horizontally with a small amplitude perturbation $\eta$

$$ \rho(x, z, t = 0) = 1 + \frac{s-1}{2} \left( \tanh \left( \frac{z}{\delta} + \eta(x, L_x) \right) + 1 \right). \qquad (21) $$
<figure><img src="image_8.png" /><figcaption>Figure 8: Temporal errors at $t = \pi/4$ with respect to the exact solution of the modified VDINSE (19) versus the time-step $\delta t$ with $Re = 1$, $Sc = 1$, $s = 10$ and $N = 128^3$ discretization grid points.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: $L_2$ norm of the divergence of the velocity field versus time during a simulation of the modified VDINSE equation (19) with $Re = 1$, $Sc = 1$, $s = 10$, $\delta t = 0.01$ and $N = 128^3$ discretization grid points.</figcaption></figure>
<table><caption>Table 2: Computational time in seconds for iterative solvers for the modified pressure equation (16) and for the test case flow (system (19)) Different values of the density ratio *s* are investigated, with 128 × 128 × 128 discretization points. Computations are performed on the 8 CPU cores of an Intel Xeon E5-2670 mounted on a Dell PowerEdge C8220 Compute Node. CPU Times are averaged over 10 time steps (10δt ≈ π/2).</caption><thead><tr><th>Iterative solver</th><th>FP ([13])</th><th>RMR</th><th>GMRES</th></tr></thead><tbody><tr><td colspan="4" style="text-align:center;">s = 2</td></tr><tr><th scope="row">Iterations</th><td>7</td><td>5</td><td>6</td></tr><tr><th scope="row">Time</th><td>0.5973</td><td>0.2124</td><td>0.2480</td></tr><tr><th scope="row">Time/iteration</th><td>0.0853</td><td>0.0425</td><td>0.0413</td></tr><tr><th scope="row">Time/iteration/thread</th><td>0.0107</td><td>0.0053</td><td>0.0052</td></tr><tr><th scope="row">Time/iteration/thread/cells</th><td>5.1e-09</td><td>2.5e-09</td><td>2.5e-09</td></tr><tr><td colspan="4" style="text-align:center;">s = 6</td></tr><tr><th scope="row">Iterations</th><td>9</td><td>7</td><td>8</td></tr><tr><th scope="row">Time</th><td>0.7737</td><td>0.3045</td><td>0.3573</td></tr><tr><th scope="row">Time/iteration</th><td>0.0860</td><td>0.0435</td><td>0.0447</td></tr><tr><th scope="row">Time/iteration/thread</th><td>0.0107</td><td>0.0054</td><td>0.0056</td></tr><tr><th scope="row">Time/iteration/thread/cell</th><td>5.1e-09</td><td>2.6e-09</td><td>2.7e-09</td></tr><tr><td colspan="4" style="text-align:center;">s = 10</td></tr><tr><th scope="row">Iterations</th><td>10</td><td>8</td><td>8</td></tr><tr><th scope="row">Time</th><td>0.8232</td><td>0.3477</td><td>0.3547</td></tr><tr><th scope="row">Time/iteration</th><td>0.0823</td><td>0.0435</td><td>0.0443</td></tr><tr><th scope="row">Time/iteration/thread</th><td>0.0103</td><td>0.0054</td><td>0.0055</td></tr><tr><th scope="row">Time/iteration/thread/cells</th><td>4.9e-09</td><td>2.6e-09</td><td>2.6e-09</td></tr></tbody></table>

The computational box size is fixed to [2L<sub>x</sub>, 2L<sub>z</sub>] = [2, 2] (discretized with N<sup>2</sup> grid points varying from 256<sup>2</sup> to 512<sup>2</sup>). The initial thickness is δ = 0.2 and *s* is the density ratio. The initial velocity field is set to zero. Gravity is added through a volumetric force

$$f = -Fr^{-1}e_z \qquad (22)$$

with *Fr* the Froude number. Within this configuration, computations without preconditioning do not converge, and numerical oscillations take over in a few iterations. Figure 10 shows density contours at different times for both a single-mode ((a1), (a2), (a3))

$$\eta = 0.1 \cos \left( \pi \frac{x}{L_x} \right) \qquad (23)$$

and a multi-mode initial perturbation ((b1), (b2), (b3))

$$\eta = 0.1 \left( \sin \left( 13\pi \frac{x}{L_x} \right) + \cos \left( 17\pi \frac{x}{L_x} \right) + \cos \left( 19\pi \frac{x}{L_x} \right) \right) \qquad (24)$$

with dimensionless numbers *Re* = 1000, *Fr* = 1, *Sc* = 1 and *s* = 10. In the single-mode case, one can see the characteristic mushroom of the Rayleigh-Taylor instability growing downward in the middle of the computational domain. As the heavier fluid penetrates the lighter fluid, the interface rolls up into two vortices as presented in Bell
<figure><img src="image_10.png" /><figcaption>Figure 10: Contour plot of the density field at different times from 2D RTI simulations using the preconditioned RMR for both the viscous and pressure operators with $s = 10$, $Re = 1000$, $Fr = 1$, $Sc = 1$ and $256^2$ grid points. In (a1), (a2), and (a3) a single-mode perturbation is injected in the initial density field. In (b1), (b2), and (b3) a multi-mode perturbation is injected in the initial density field. Figure 11 presents the convergence history for the pressure operator corresponding to (b3).</figcaption></figure>

and Marcus [1]. In the multi-mode case, the waves begin to grow independently of one another before they start to interact strongly with each other reproducing patterns similar to those presented in Bell and Marcus [1]. In both cases, it can be seen that the interface is clearly defined without numerical oscillations, which validates the robustness of the method. Figure 11 presents the convergence history of the three studied methods (FP, RMR and GMRES) when dealing with the modified pressure equation (16) for the same range of parameters and a multiple-mode initial perturbation at time $t = 3.0$. It can be seen that the convergence rate is independent of the mesh size and both the RMR and GMRES algorithms converge towards a residual value lower than $10^{-7}$ which demonstrate the effectiveness of the method in dealing with steep problems. For this case, the FP method does not achieve a good convergence.

The reliability of the proposed method is now illustrated through a representative three-dimensional variable-density flow case. In that respect, the variable-density swirling jet studied by Di Pierro and Abid [13] is selected. The Reynolds number — based on
<figure><img src="image_11.png" /><figcaption>Figure 11: Convergence history of the different preconditioned iterative solvers to solve the modified pressure equation equation (16) in the case (b3) of figure 10.</figcaption></figure>

the base flow stream-wise velocity — is fixed to $Re = 100$ and the others dimensionless control parameters are set to $Sc = 1$, $s = 2$ and $q = 0.5$. The domain size is $2\pi \times 2\pi \times 4\pi$. We consider the following equilibrium state:

$$ \rho_0(r) = 1 + (s - 1) \exp -r^2, \qquad (25a) $$

$$ \omega_x(r) = 2q \exp -r^2, \qquad (25b) $$

$$ V_x(r) = \exp -r^2, \qquad (25c) $$

where $V_0 = (0, V_\theta, V_x)^T$ is the base velocity field (with $V_\theta e_\theta = -(\nabla^2)^{-1} (\nabla \times \omega_x e_x)$), $\omega_x$ is the axial vorticity, and $q$ is the swirl number. We then capture the dynamics of a perturbation $\phi(r) \exp(i(kx + m\theta))$ superimposed to the base flow with a very small amplitude. In that purpose, $256^3$ discretization points are used. It is found that the dynamics is driven for long times by the azimuthal mode $m = 4$ and the fundamental stream-wise wave number $k = 1$ in agreement with the linear stability theory (LST, see [13] for details about the linear stability analysis). Figure 12 shows the evolution of the axial velocity perturbation. We then compare the exponential growth rate extracted from the DNS database with the linear stability results. In that case, the linear stability analysis gives a temporal growth rate $\omega_i = 0.274$. For long times, we found a temporal growth rate from the DNS

$$ \omega_{i, \text{DNS}} = \frac{1}{2(t_1 - t_0)} \log \left( \frac{\|w(t_1)\|_2^2}{\|w(t_0)\|_2^2} \right) = 0.2746 \approx 0.275 \qquad (26) $$

which nearly matches the one computed from the LST. The latter is calculated from the norm of the axial velocity perturbation $w$ with $t_0 = 2.0$ and $t_1 = 5.0$. Figure 13 shows
the isosurfaces of the density field in the non-linear saturated regime associated with such flow case. The interface is well-defined and no oscillations appear. It further illustrates the robustness of the present method to deal with representative variable-density flow case.

## 5. Conclusions

Performances of preconditioners based on the constant-density operators for the numerical resolution of the variable-density and incompressible Navier-Stokes equations with spectral accuracy and second-order time accuracy are investigated. We first observed that test cases presented here did not converge without any preconditioning, which motivated this study. It is shown that these preconditioners highly reduce the condition number of the variable-density elliptic operators for both the pressure and velocity fields equations. Coupled with widely used iterative solvers (Conjugate Gradient, Richardson Minimal Residual, and General Minimal Residual), we give strong evidence that the proposed numerical methods highly enhance the convergence of the implicit systems. Finally, the precision and the robustness of the method is further illustrated on some representative variable-density flow cases.

We thus believe that the present study can serve as a guide for the development of faster and more accurate DNS incompressible solvers that take into account large density variations. For the extension of our method to high-performance simulation on distributed-memory computers, one may recall that the major drawback of spectral methods lies in their global approximation. However, some solutions are suggested in the literature and successfully numerically tested. For instance, one may cite the work of some members of the same team dealing with high-scalability spectral code on high-performance distributed-memory [32]. Especially, the authors show that the parallelization strategy for spectral schemes based on a domain decomposition method — where the computational domain is subdivided along spatial directions into subdomains — exhibits a good scalability and a very fast wall-clock time per iteration on HPC platforms. In the present work, the three-dimensional computations use the same strategy.

## 6. Acknowledgements

The authors thank the "Fédération Lyonnaise de Modélisation et Sciences Numérique" for providing computational facilities on the computer center P2CHPD at "Université Claude Bernard Lyon 1".

The authors also thank the group "LyonCalcul" for fruitful discussions.

# Appendix A. Viscous operator iterative solvers convergence history

In this appendix, we show the convergence history of the CG, BiCG, BiCGStab, RMR and GMRES algorithms for the resolution of the viscous equation $V_p(\mathbf{u}) = \mathbf{b}$, where the viscous operator $V_p$ is defined equation (13). A wide range of $(a, s)$ is investigated.
<figure><img src="image_12.png" /><figcaption>Figure 12: Evolution of the norm of the axial velocity perturbation during a DNS of a dense swirling jet with $Re = 100$, $Sc = 1$, $s = 2$, $q = 0.5$, $m = -4$ and $k = 1$. The growth rate obtained is $\omega_{i,\text{DNS}} = 0.275$, the growth rate predicted by the linear instability analysis is $\omega_i = 0.274$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: Isosurfaces of the density field perturbation in the non-linear saturated regime at $t = 40$ from a DNS of a dense swirling jet with $Re = 100$, $Sc = 1$, $s = 2$, $q = 0.5$, $m = -4$ and $k = 1$.</figcaption></figure>
# Appendix B. Pressure operator iterative solvers convergence history

In this appendix, we show the convergence history of the CG, BiCG, BiCGStab, RMR, GMRES and FP algorithms for the resolution of the modified pressure equation (16) A wide range of density ratios $s$ is investigated.

## References

[1] J. B. Bell, D. L. Marcus, A second-order projection method for variable-density flows, Journal of Computational Physics 101 (1992) 334–348. URL: https://doi.org/10.1016/0021-9991(92)90011-m. doi:10.1016/0021-9991(92)90011-m.

[2] A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, M. L. Welcome, A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations, Journal of Computational Physics 142 (1998) 1–46. URL: https://doi.org/10.1006/jcph.1998.5890. doi:10.1006/jcph.1998.5890.

[3] C. Calgaro, E. Creusé, T. Goudon, An hybrid finite volume–finite element method for variable density incompressible flows, Journal of Computational Physics 227 (2008) 4671–4696. URL: https://doi.org/10.1016/j.jcp.2008.01.017. doi:10.1016/j.jcp.2008.01.017.

[4] P. Concus, G. H. Golub, Use of fast direct methods for the efficient numerical solution of non-separable elliptic equations, SIAM Journal on Numerical Analysis 10 (1973) 1103–1120. URL: https://doi.org/10.1137/2F0710092. doi:10.1137/0710092.

[5] A. Duffy, A. Kuhnle, M. Sussman, An improved variable density pressure projection solver for adaptive meshes, 2002. Unpublished, see https://www.math.fsu.edu/sussman/MGAMR.pdf.

[6] A. W. Cook, W. Cabot, P. L. Miller, The mixing transition in Rayleigh–Taylor instability, Journal of Fluid Mechanics 511 (2004) 333–362. URL: https://doi.org/10.1017/s0022112004009681. doi:10.1017/s0022112004009681.

[7] M. El Ouafa, S. Vincent, V. Le Chenadec, Monolithic solvers for incompressible two-phase flows at large density and viscosity ratios, Fluids 6 (2021) 23. URL: https://doi.org/10.3390/fluids6010023. doi:10.3390/fluids6010023.

[8] J.-L. Guermond, A. Salgado, A splitting method for incompressible flows with variable density based on a pressure Poisson equation, Journal of Computational Physics 228 (2009) 2834–2846. URL: https://doi.org/10.1016/j.jcp.2008.12.036. doi:10.1016/j.jcp.2008.12.036.

[9] A. W. Cook, W. Cabot, P. L. Williams, B. J. Miller, B. R. de Supinski, R. K. Yates, M. L. Welcome, Tera-scalable algorithms for variable-density elliptic hydrodynamics with spectral accuracy, in: ACM/IEEE SC 2005 Conference (SC’05), IEEE, IEEE, 2005, pp. 60–60. URL: https://doi.org/10.1109/sc.2005.70. doi:10.1109/sc.2005.70.

[10] S. Ravier, M. Abid, M. Amielh, F. Anselmet, Direct numerical simulations of variable-density plane jets, Journal of Fluid Mechanics 546 (2005) 153. URL: https://doi.org/10.1017/s0022112005006993. doi:10.1017/s0022112005006993.

[11] R. Knikker, A comparative study of high-order variable-property segregated algorithms for unsteady low Mach number flows, International Journal for Numerical Methods in Fluids 66 (2011) 403–427. URL: https://doi.org/10.1002/2Ffld.2261. doi:10.1002/fld.2261.

[12] P. Haldenwang, G. Labrosse, S. Abboudi, M. Deville, Chebyshev 3-D spectral and 2-D pseudospectral solvers for the Helmholtz equation, Journal of Computational Physics 55 (1984) 115–128. URL: https://doi.org/10.1016/2F0021-9991/2884/2990018-4. doi:10.1016/0021-9991(84)90018-4.

[13] B. Di Pierro, M. Abid, A projection method for the spectral solution of non-homogeneous and incompressible Navier-Stokes equations, International Journal for Numerical Methods in Fluids 71 (2012) 1029–1054. URL: https://doi.org/10.1002/fld.3700. doi:10.1002/fld.3700.

[14] D. A. Frank-Kamenetskii, Diffusion and Heat Exchange in Chemical Kinetics, volume 2171, Princeton University Press, 1955. URL: https://doi.org/10.1515/9781400877195. doi:10.1515/9781400877195.

[15] A. V. Kazhikhov, S. S. Smagulov, The correctness of boundary-value problems in a diffusion model of an inhomogeneous liquid, in: Soviet Physics Doklady, volume 22, 1977, p. 249.

[16] S. N. Antontsev, A. V. Kazhiktov, V. N. Monakhov, Boundary Value Problems in Mechanics of Nonhomogeneous Fluids, Elsevier, North Holland Amsterdam, 1990. URL: https://doi.org/10.1016/s0168-2024(08)x7006-7. doi:10.1016/s0168-2024(08)x7006-7.
<figure><img src="image_14.png" /><figcaption>Figure B.14: Convergence history for the five preconditioned iterative solvers to solve $V_\rho(u) = b$ with $s = 2$ and different values of $a$, with Fourier differentiation (left side) and Chebyshev differentiation (right side) and $N = 256$ collocation points.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure B.15: Convergence history for the five preconditioned iterative solvers to solve $V_\rho(u) = b$ with $a = 10^{-2}$ and different values of $s$, with Fourier differentiation (left side) and Chebyshev differentiation (right side) and $N = 256$ collocation points.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure B.16: Convergence history of the different preconditioned iterative solvers to solve the modified pressure equation (16) with different values of $s$, and $N = 256$ Chebyshev collocation points.</figcaption></figure>
[17] F. Guillén-González, P. D. Damázio, M. A. Rojas-Medar, Approximation by an iterative method for regular solutions for incompressible fluids with mass diffusion, Journal of Mathematical Analysis and Applications 326 (2007) 468-487. URL: https://doi.org/10.1016%2Fj.jmaa.2006.03.009. doi:10.1016/j.jmaa.2006.03.009.

[18] M. E. Brachet, D. I. Meiron, S. A. Orszag, B. G. Nickel, R. H. Morf, U. Frisch, Small-scale structure of the Taylor-Green vortex, Journal of Fluid Mechanics 130 (1983) 411. URL: https://doi.org/10.1017/s0022112083001159. doi:10.1017/s0022112083001159.

[19] R. M. Kerr, Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence, Journal of Fluid Mechanics 153 (1985) 31. URL: https://doi.org/10.1017/s0022112085001136. doi:10.1017/s0022112085001136.

[20] A. Vincent, M. Meneguzzi, The spatial structure and statistical properties of homogeneous turbulence, Journal of Fluid Mechanics 225 (1991) 1-20. URL: https://doi.org/10.1017/s0022112091001957. doi:10.1017/s0022112091001957.

[21] P. Schlatter, N. A. Adams, L. Kleiser, A windowing method for periodic inflow/outflow boundary treatment of non-periodic flows, Journal of Computational Physics 206 (2005) 505-535. URL: https://doi.org/10.1016/j.jcp.2004.12.015. doi:10.1016/j.jcp.2004.12.015.

[22] C. Tadjeran, Stability analysis of the Crank-Nicholson method for variable coefficient diffusion equation, Communications in Numerical Methods in Engineering 23 (2006) 29-34. URL: https://doi.org/10.1002/cnm.879. doi:10.1002/cnm.879.

[23] D. L. Brown, R. Cortez, M. L. Minion, Accurate projection methods for the incompressible Navier-Stokes equations, Journal of Computational Physics 168 (2001) 464-499. URL: https://doi.org/10.1006/jcph.2001.6715. doi:10.1006/jcph.2001.6715.

[24] S. A. Orszag, Spectral methods for problems in complex geometries, Journal of Computational Physics 37 (1980) 70-92. URL: https://doi.org/10.1016/0021-9991(80)90005-4. doi:10.1016/0021-9991(80)90005-4.

[25] R. Peyret, Spectral Methods for Incompressible Viscous Flow, number 148 in Applied mathematical sciences, Springer New York, New York, 2002. URL: https://doi.org/10.1007/978-1-4757-6557-1. doi:10.1007/978-1-4757-6557-1.

[26] C. Canuto, A. Quarteroni, Preconditioned minimal residual methods for Chebyshev spectral calculations, Journal of Computational Physics 60 (1985) 315-337. URL: https://doi.org/10.1016/0021-9991(85)90010-5. doi:10.1016/0021-9991(85)90010-5.

[27] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in C, Cambridge university press Cambridge, 1992.

[28] B. Di Pierro, M. Abid, Spatiotemporal instability of a variable-density Batchelor vortex, Journal of Fluid Mechanics 703 (2012) 49-59. URL: https://doi.org/10.1017%2Fjfm.2012.185. doi:10.1017/jfm.2012.185.

[29] F. Alizard, A. Cadiou, L. Le Penven, B. Di Pierro, M. Buffat, Space-time dynamics of optimal wavepackets for streaks in a channel entrance flow, Journal of Fluid Mechanics 844 (2018) 669-706. URL: https://doi.org/10.1017%2Fjfm.2018.191. doi:10.1017/jfm.2018.191.

[30] O. Thual, Transition vers la turbulence dans les systèmes dynamiques apparentés la convection, Ph.D. thesis, Université de Nice-Sophia Antipolis, 1986.

[31] Y. Saad, M. H. Schultz, Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing 7 (1986) 856-869. URL: https://doi.org/10.1137/0907058. doi:10.1137/0907058.

[32] J. Montagnier, A. Cadiou, M. Buffat, L. Le Penven, Towards petascale spectral simulations for transition analysis in wall bounded flow, International Journal for Numerical Methods in Fluids 72 (2012) 709-723. URL: https://doi.org/10.1002%2Ffld.3758. doi:10.1002/fld.3758.