Randomized resolvent analysis

TL;DR

This paper introduces a randomized numerical linear algebra method to efficiently perform resolvent analysis on high-Reynolds-number turbulent flows, reducing computational costs while capturing dominant flow features.

Contribution

It presents a novel randomized approach using sketching and low-rank approximations to accelerate resolvent analysis for complex turbulent flows.

Findings

Significant reduction in computational time and memory usage.

Accurate extraction of dominant flow response and forcing modes.

Application demonstrated on turbulent airfoil flow.

Abstract

Performing global resolvent analysis for high-Reynolds-number turbulent flow calls for the handling of a large discrete operator. Even though such large operator is required in the analysis, most applications of resolvent analysis extracts only a few dominant resolvent response and forcing modes. Here, we consider the use of randomized numerical linear algebra to reduce the dimension of the resolvent operator for achieving computational speed up and memory saving compared to the standard resolvent analysis. To accomplish this goal, we utilize sketching of the linear operator with random test matrices with a Gaussian distribution and with insights from the base flow incorporated to perform singular value decomposition on a low-rank matrix holding dominant characteristics of the full resolvent operator. The strength of the randomized resolvent analysis is demonstrated on a turbulent separated flow over an airfoil. This randomized approach clears the path towards tackling resolvent analysis for higher-Reynolds number bi- and tri-global base flows.

Figures

Tables

Table 1. Table 1: Comparison of the computational time and memory consumption for the implicitly restarted Arnoldi iteration ( svds in MATLAB) with the present randomized approach for different parameter setups with the number of singular values ( n ev subscript 𝑛 ev n_{\text{ev}} ), Krylov subspace dimension ( dim ​ ( 𝒮 ) dim 𝒮 \text{dim}(\mathcal{S}) ), and tolerance.

Iteratively restarted Arnoldi method (MATLAB svds)
$n_{\text{ev}}$	$\text{dim}\left(𝒮\right)$	tolerance	time (sec)	memory
$10$	$30$	1E-14	4185	78.6 GB
$5$	$15$	1E-14	2764	78.6 GB
$2$	$6$	1E-14	1486	78.6 GB
$10$	$30$	1E-05	2245	78.6 GB
$10$	$15$	1E-14	4194	78.6 GB
Randomized resolvent (present)
$n_{\text{ev}}$	$k$		time (sec)	memory
$2$	$2$		354	28.2 GB
$5$	$5$		462	28.4 GB
$10$	$10$		615	28.8 GB

Equations

$$\frac{\partial q^{\prime }}{\partial t}={\mathbf{L}}_{\bar{q}}{q}^{\prime }+f^{\prime },$$

$$-i\omega {\hat{q}}_{\omega }={\mathbf{L}}_{\bar{q}}{\hat{q}}_{\omega }+{\hat{f}}_{\omega },$$

$${\hat{q}}_{\omega }=\mathbf{A}{\hat{f}}_{\omega },$$

$$\mathbf{A}={\left[-i\omega \mathbf{I}-{\mathbf{L}}_{\bar{q}}\right]}^{-1}\in {\mathbb{C}}^{m\times m}$$

$$\mathbf{A}=\mathbf{U}\Sigma {\mathbf{V}}^{\ast },$$

$$\mathbf{Y}=\mathbf{A}\mathbf{\Omega }.$$

$$\text{boldmath}B=\tilde{\mathbf{U}}\mathbf{\Sigma }\text{boldmath}{V}^{\ast }$$

$$\mathbf{A}\approx \mathbf{Q}\tilde{\mathbf{U}}\mathbf{\Sigma }{\mathbf{V}}^{\ast }$$

$$\mathbf{A}{\hat{\mathbf{v}}}_1={\hat{\mathbf{u}}}_1{\sigma }_1.$$

$${\mathbf{\Omega }}_p=\mathsf{diag}(\mathbf{\Phi })\mathbf{\Omega }.$$

Full text

Keywords

resolvent analysis, instabilities, randomized methods

Suggested keywords

I Introduction

One of the central questions in turbulence and flow control is concerned with the evolution of perturbations. Gaining detailed understanding of the perturbation dynamics in turbulent flows is a daunting task, due to the complex nonlinear dynamics that takes place over a broad range of spatial and temporal scales. To modify the global flow characteristics with control, it is necessary that the nonlinear interactions involving the actuation input (perturbation) become appropriately large to alter the base flow. For flow control, we need not track all possible ways in which the actuation input can modify the flow, but instead we can focus on the dominant directions in which the perturbation can be amplified. This notion has led to the modal analysis based approaches [1, 2, 3, 4], including the global stability analysis [5] and the resolvent analysis [6].

In the presence of sustained perturbations or forcing inputs, the linear system response can be described by the transfer function from control theory. This linear analysis is greatly simplified when the input to the system is sinusoidal and leads to the well-known Bode plots that reveal the gain and phase response of the system over a range of forcing frequencies. The transfer function that relates the system input to the output is called the resolvent and its analysis has been extended to fluid flows by Trefethen et al. [6]. The resolvent analysis is based on the pseudospectral analysis and has been used to study the transient energy growth [6] as well as the harmonic response of the system [7]. These initial studies of resolvent analysis were performed about stable laminar flows.

An extension to resolvent analysis to study turbulent flows was presented by McKeon and Sharma [8]. They considered the nonlinear advection term to be the self-sustained input within the natural feedback loop of the fluid flow. This viewpoint has enabled the use of time-averaged base flows to reveal the input-output dynamics of turbulent flows. Moreover, discounting or finite-time horizon based extension of the resolvent analysis has enabled resolvent analysis to study flows with unstable base states [9, 10]. Because resolvent analysis can determine the most amplified forcing and response directions, it serves as a powerful analytical tool to find effective active and passive flow control techniques [10, 11].

The resolvent analysis needs two key ingredients: (i) the base flow and (ii) the linearized Navier–Stokes operator. It is known that the accuracy of the base flow and the spatial discretization of the linear operators is critical for extracting response characteristics correctly [10]. The need for accurate discretization of the linearized Navier–Stokes operators calls for sufficient grid resolution and appropriate computational domain size. As such, the discrete resolvent operator becomes large with size $m\times m$. Here, $m$ is essentially the number of variables times the size of the grid, which can easily be upward of $\mathcal{O}(10^{6})$ for turbulent flows [12]. For resolvent analysis, the singular value decomposition (SVD) needs to be performed on the large resolvent operator with a taxing operation count of $\mathcal{O}(m^{3})$. To enable resolvent analysis at high Reynolds numbers, we must find a relief to perform SVD of the resolvent operator.

Although the resolvent analysis is performed on a very large matrix, only the leading forcing and response modes are generally sought. Based on the amount of necessary matrix data used to perform the analysis, the desired output is only a very small fraction of the input data size. For this reason, it would be natural to consider that all elements of the resolvent matrix are not necessary to determine the leading resolvent modes. Within the resolvent framework, the core of the computations lies with the SVD. In order to handle a large operator for SVD to find the leading modes, we can consider subsampling the matrix of interest and perform the SVD on the low-order representation of this large matrix. Randomized numerical linear algebra has recently emerged as an effective technique to reduce a large matrix to its low-order representation [13, 14, 15], for applications including big data compression and data transfer. The key idea is to pass a randomly generated low-rank test matrix through the large matrix to obtain the so-called sketch of the full matrix. This sketch is low-rank but holds key information about the full matrix and can be used to derive appropriate bases to represent the full matrix in a low-dimensional manner [17, 16, 13, 15]. Randomized techniques can be incorporated into SVD to achieve tremendous computational and memory savings [16, 18, 13].

In recent years, modal analyses are tackling flows over complex geometries and high-Reynolds number flows with increasingly large degrees of freedom [3, 4, 19, 22, 23, 20, 21]. To further aid this endeavor, randomized SVD have been incorporated into data-based modal analysis techniques, including the proper orthogonal decomposition [18] and dynamic mode decomposition [24]. For global operator-based analyses of high-Reynolds-number flows [19, 22, 23], the leading singular values and modes can be determined with significant reduction in computational costs with the aid of randomized techniques. The randomized technique presented in this study can greatly expand the applicability of the resolvent analysis to high-Reynolds number flows with multiple inhomogeneous spatial directions.

In fact, randomized SVD has been adopted in the one-dimensional resolvent analysis of turbulent channel flow by Moarref et al. [25]. Since the details on the use of the randomized techniques for large-scale resolvent analysis is not available, one of our objectives is to provide detailed guidance on the use of randomized resolvent analysis for multi-dimensional turbulent base flows. Moreover, the error bound of randomized SVD is generally not guaranteed and needs to be further characterized. We also note that the sketching strategy can be further improved with insights from the base flow. In the present study, we aim to address these issues on the use of randomized SVD in the context of resolvent analysis for fluid flows.

In what follows, we present the randomized resolvent analysis and demonstrate its use for turbulent flow over a canonical airfoil. In section II, we introduce the randomized resolvent analyses and discuss how randomized SVD can be implemented without the explicit generation of inverse matrices. We also propose an alternative route of recovering the response modes (left singular vectors) with significantly enhanced accuracy. In section III, the present randomized technique is applied to the resolvent analysis of turbulent flow over a NACA 0012 airfoil. We assess the level of error and the convergence behavior with respect to resolvent modes and gains obtained from the use of randomized SVD. To further enhance the accuracy of randomized resolvent analysis, we introduce a physics-informed sampling technique that leverages the insights from the base flow. At last, we provide concluding remarks in section IV.

II Approach

II.1 Full resolvent analysis

Let us consider the flow state $\pmb{q}\in\mathbb{R}^{m}$ as a sum of the time-invariant base state $\bar{\pmb{q}}$ and the statistically stationary fluctuating component $\pmb{q}^{\prime}$. With this Reynolds decomposed flow variable and appropriate discretization, we can express the discrete Navier–Stokes equations as

[TABLE]

where $\bm{L}_{\bar{\pmb{q}}}\in\mathbb{R}^{m\times m}$ is the linearized Navier–Stokes operator about the base state $\bar{\pmb{q}}$ and $\pmb{f}^{\prime}$ collects the nonlinear terms and the external forcing inputs. We gather the nonlinear terms as external forcing in the turbulent mean flow following the perspective of McKeon and Sharma [8], Farrell and Ioannou [26], and Schmid [27]. For the traditional resolvent analysis, $\bar{\pmb{q}}$ is chosen to be the stable laminar equilibrium state such that $\pmb{f}^{\prime}$ can be considered as the forcing input to the system with the nonlinear term neglected [7]. More recently, turbulent mean flows has been used for ${\bar{\pmb{q}}}$ with $\pmb{f}^{\prime}$ representing the nonlinear terms as sustained forcing input within the natural feedback system [8].

We can consider the Fourier transform $\left[\pmb{q}^{\prime}(\pmb{x},t),\pmb{f}^{\prime}(\pmb{x},t)\right]=\int_{-\infty}^{\infty}[\hat{\pmb{q}}_{\omega}(\pmb{x}),\hat{\pmb{f}}_{\omega}(\pmb{x})]e^{-i\omega t}{\rm d}\omega$ and express the relationship between $\pmb{q}^{\prime}$ and $\pmb{f}^{\prime}$ in frequency space as

[TABLE]

where $\omega$ is the frequency. Note that spatial Fourier transform can also be incorporated if directional homogeneity is present. For stable base flows, $\omega$ can be chosen to be real. To extend resolvent analysis to unstable base flows, we can consider the use of finite-time/discounted analysis [9, 10] by choosing a complex frequency $\omega=\omega_{r}+i\beta$, where both $\omega_{r}$ and $\beta$ are real and $\beta$ discounts the modal growth rate of $\bm{L}_{\bar{\pmb{q}}}$. The input-output relationship between $\hat{\pmb{f}}$ and $\hat{\pmb{q}}$ can be found from (2) as

[TABLE]

where

[TABLE]

is referred to as the resolvent operator [6, 7, 8]. It serves as a transfer function that amplifies (or attenuates) the harmonic forcing input $\hat{\pmb{f}}_{\omega}$ and maps it to the response $\hat{\pmb{q}}_{\omega}$. The goal of resolvent analysis is to identify the dominant directions along which $\hat{\pmb{f}}_{\omega}$ can be most amplified through $\bm{A}$ to form the corresponding responses in $\hat{\pmb{q}}_{\omega}$. This question is addressed by the SVD of

[TABLE]

where $\bm{V}^{*}$ denotes the Hermitian of $\bm{V}$. Resolvent analysis interprets left and right singular vectors $\bm{U}=[\hat{\pmb{u}}_{1},\hat{\pmb{u}}_{2},\dots,\hat{\pmb{u}}_{m}]\in\mathbb{C}^{m\times m}$ and $\bm{V}=[\hat{\pmb{v}}_{1},\hat{\pmb{v}}_{2},\dots,\hat{\pmb{v}}_{m}]\in\mathbb{C}^{m\times m}$ respectively as response modes and forcing modes, with the magnitude-ranked singular values $\pmb{\Sigma}={\sf{diag}}(\sigma_{1},\sigma_{2},\dots,\sigma_{m})\in\mathbb{R}^{m\times m}$ being the amplification (gain) for the corresponding forcing-response pair. For unstable base flows, it is important that a finite-time window is chosen with $\beta$ larger than the highest growth rate such that the resolvent analysis reveals the input-output relationship on a shorter time scale than that of the base flow instability [9, 10].

Performing the SVD of $\bm{A}\in\mathbb{C}^{m\times m}$ requires an theoretical operation count of $\mathcal{O}(m^{3})$. In practice, some algorithms can reduce this operation count when only a few singular values are to be recovered, while still being computationally taxing for large $m$ [28]. Such cases are encountered in high-Reynolds number flows and bi/tri-global analysis settings. However, we note that many applications of resolvent analysis call only for the dominant forcing and response modes $[\hat{\pmb{v}}_{1},\hat{\pmb{u}}_{1}]$ associated with the highest gain $\sigma_{1}$. This is appropriate when the first gain $\sigma_{1}$ is much larger than the rest of the gains $\sigma_{j>1}$ and shows a quick roll off. Such condition is related to non-normality of linear operator $\bm{L}_{\bar{\pmb{q}}}$, which is encountered for flows with strong shear and separation [29].

When the contributions from the higher-order modes are neglected, it is referred to as the rank-1 assumption for which the flow response from forcing $\hat{\pmb{f}}$ is approximated as $\hat{\pmb{q}}\approx\hat{\pmb{u}}_{1}\sigma_{1}\langle\hat{\pmb{v}}_{1},\hat{\pmb{f}}\rangle$, provided that $\sigma_{1}\gg\sigma_{2}$ and $\hat{\pmb{f}}$ has reasonable magnitude along $\hat{\pmb{v}}_{1}$. For seeking only the dominant modal insights from resolvent analysis, we discuss a remedy for performing large-scale resolvent analysis [19, 22] in a computationally tractable manner below.

II.2 Randomized resolvent analysis

For a flow that have dominant structures, we consider a low-rank representation of the resolvent operator $\bm{A}$. That is, instead of directly performing SVD of $\bm{A}$ and obtaining the leading-mode representation, we seek a low-rank representation of $\bm{A}$ and perform the SVD of the low-rank version of $\bm{A}$. We can consider finding an appropriate low-dimensional basis to project the large resolvent operator on a suitable subspace to derive the low-rank resolvent approximation.

The action of a full matrix on a vector should reveal some insights on which components are modified in the dominant directions. In the case of flow that can be described with the rank-1 approximation, there should be a low number of dominant directions. This very point can be taken advantage of through what is known as sketching in numerical linear algebra. Sketching refers to a procedure in which a tall and skinny test matrix $\bm{\Omega}\in\mathbb{R}^{m\times k}$ (or $\mathbb{C}^{m\times k}$), where $k\ll m$, is passed through $\bm{A}$

[TABLE]

Here, matrix $\bm{Y}\in\mathbb{C}^{m\times k}$ is called the sketch of the input matrix $\bm{A}$ [16, 13, 15]. The test matrix $\bm{\Omega}$ can be constructed using random values with Gaussian distribution [30] and be weighted by any input matrix insight, as will be discussed in section III.4 for a physics-inspired random test matrix. As the sketch holds the dominant influence of $\bm{A}$, we can consider orthonormalizing $\bm{Y}$ using a QR decomposition to form the orthonormal basis with $\bm{Q}\in\mathbb{C}^{m\times k}$ upon which we can project the full matrix $\bm{A}$ to derive its low-rank approximation. In this way, it is possible to approximate $\bm{A}$ for a rank $k\ll m$ as long as this approximation preserves the features of the leading modes.

Given this $\bm{Q}$, a low-rank approximation of $\bm{A}$ can be found as $\bm{A}\approx\bm{Q}\bm{Q^{*}}\bm{A}$ [13]. We can view this as a low-rank decomposition of $\bm{A}\approx\bm{Q}\bm{B}$, where $\bm{B}=\bm{Q^{*}}\bm{A}\in\mathbb{C}^{k\times m}$. It is this reduced matrix $\bm{B}$ upon which we can perform the SVD

[TABLE]

Hence, as a low-rank approximation, we now have

[TABLE]

where we can consider $\bm{U}\approx\bm{Q}\bm{\tilde{U}}$. This process is the randomized SVD [13], where the sketch $\bm{Y}$ was used to derive $\bm{Q}$. This approximation almost always satisfies $\|\bm{A}-\bm{Q}\bm{Q^{*}}\bm{A}\|\leq\left(1+9\sqrt{k+p}\sqrt{m}\right)\sigma_{k+1}$, where $p$ is the oversampling parameter, which is applied to build an approximation of rank $k$ while projecting the matrix $\bm{A}$ to the low-dimensional subspace with $(k+p)$ vectors. With this overall approach, the computational cost for SVD is reduced to $\mathcal{O}(mk^{2})$ instead of $\mathcal{O}(m^{3})$ for the full SVD. In our implementation, we use $\bm{v}_{1}$ from (8) and retrieve the leading singular value and left singular vector through

[TABLE]

The singular value and vector can be separated by noticing that $\|\hat{\bm{u}}_{1}\|=1$. Applied to resolvent analysis, the last equation provides more accurate leading singular value $\sigma_{1}$ and left singular vector $\hat{\bm{u}}_{1}$ compared to the original randomized SVD algorithm [13]. The same operation can be used to recover the higher-order modes, with better accuracy than using the original algorithm [13]. For applications where high-order modes and orthogonality are desired, we can solve for $\bm{U}\bm{\Sigma}$ and compute its SVD. In the present randomized resolvent analysis, we emphasize that only the discrete linear operator $\bm{L}_{\bar{\pmb{q}}}$ is needed for sketching $\bm{Y}$ and to find the reduced matrix $\bm{B}$. Unlike the original resolvent method, matrix linear solvers are used to avoid calling for the inverse within the resolvent operator. The resulting algorithm constitutes the randomized resolvent analysis summarized in Algorithm 1.

To utilize the randomized SVD for resolvent analysis, we must be aware that the resolvent operator $\bm{A}=\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]^{-1}$ contains an inverse operation in its definition, which need not be numerically performed. We do not intend to perform an inverse operation within $\bm{A}$ in the present work. In the full resolvent analysis, when the matrices become too large and the inverse can not be performed (which is likely the case for 2D and 3D problems), one can focus on modes corresponding to the smallest singular values of $\bm{A}^{-1}$ to find those for the largest singular values of $\bm{A}$. Similar approaches have avoided the inverse computation, including the work by Jeun et al. [19].

Both full resolvent and randomized resolvent analyses are shown schematically in figure 1. Notice that we are not interested in all singular values and vectors of $\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]^{-1}$, but only in a few subset of the largest $\sigma_{j}$ and their corresponding $\hat{\bm{u}}_{j}$ and $\hat{\bm{v}}_{j}$. In the randomized resolvent analysis, we can approximate a low-rank representation of it using $\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]$. Figure 1 shows an adaptation of the procedure from Halko et al. [13] in order to compute the largest singular values of the resolvent without performing its inverse. In the randomized resolvent analysis, we solve a linear system with $\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]$, the columns of the random matrix $\bm{\Omega}$ form the right-hand side and the sketch columns of $\bm{Y}$ are the unknowns. By doing so, we sketch $\left[-i\omega-\bm{L}_{\bar{\pmb{q}}}\right]^{-1}$ without finding the actual inverse matrix. The same procedure is performed to project the matrix to the low-dimensional subspace. The matrices are re-arranged in a way that the projection is performed using $\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]$, but results in the low-dimensional projection of $\left[-i\omega-\bm{L}_{\bar{\pmb{q}}}\right]^{-1}$ instead. Figure 1 also illustrates the procedures for recovering the left singular vectors and the singular values from the original algorithm [13] and from the present implementation.

II.2.1 Oversampling and power iteration schemes

Randomized algorithms can incorporate two additional procedures to improve performance and accuracy. Namely, they are oversampling [31] and power or subspace iterations [18, 13, 32]. Oversampling sketches the input matrix using $(k+p)$ vectors (with $p$ extra vectors) and increases the low-dimensional subspace to accurately recover a smaller quantity of singular values $k$. For the randomized resolvent analysis, oversampling has the same outcome, in practice, of selecting a larger $k$ and the influence of $k$ will be discussed at the end of section III.2. When $k$ becomes large, it should be noticed that the memory consumption increases. Even for large sparse matrices, the sketch matrix and the subsequently reduced matrices that are formed are generally dense, which adds a computational burden.

The second procedure is the power or subspace iterations. These methods are a powerful tool when the singular values of the matrix decay slowly. For example, this type of spectral behavior appeared in the input-output analysis performed by Jeun et al. [19] for jet flows. The method consists of performing additional iterations after the sketch $\bm{Y}$ is evaluated. It should however be realized that such procedure calls for additional linear solvers, which is the most time consuming operation in Algorithm 1. For power iterations, one must compute the adjoint $\bm{A}\bm{A}^{*}$, where $\bm{A}^{*}$ is the Hermitian of $\bm{A}$, $q$ times and solve the linear system $k$ times. For subspace iterations, additional $QR$ decompositions are necessary, and the number of additional linear systems to be solved will be $q$ times $k$. In practice, small values of $q$ improves the accuracy of the results substantially. More information on the general applicability of the subspace and power iterations are discussed by Halko et al. [13] and Erichson et al. [31]. In section III.4, we present another option for improving accuracy of the overall technique in a computationally inexpensive manner by constructing a physics-informed random test matrix $\bm{\Omega}$.

II.2.2 Test matrix $\bm{\Omega}$

The standard choice for the test matrix $\bm{\Omega}$ is the random matrix generated with a normal distributions. Such matrix is known to present excellent performance and accuracy [15]. For some cases, especially for very large matrices or when the singular values present slow decay, larger values of $k$ may be necessary to better approximate the matrix in the low-dimensional subspace. When large values of $k$ are used, orthonormalization of the columns of the test matrix can be considered to improve numerical stability [33, 13]. The test matrix can also be generated using a Rademacher distribution [17]. It is also possible to build an ultrasparse matrix with Rademacher distribution in the non-sparse entries which allows for the control of cost, stability and reliability in the operations [15, 17].

When randomized SVD is applied, there is no a priori knowledge of the structure of the matrix. However, in the present application, we know how the resolvent operator is constructed. This theoretical insight can be used to build a random test matrix that outperforms the standard normal distribution matrix. We later propose a physics-informed test matrix $\bm{\Omega}$ that can focus our sketching operation for regions of physical importance. In our application, the dominant directions are related to regions with the presence of high shear. The results from this approach will be discussed in Section III.4.

III Randomized resolvent analysis of turbulent post-stall flow

We demonstrate the use of randomized resolvent analysis on turbulent flow over a NACA 0012 airfoil. In this example, the randomized resolvent analysis will be applied a resolvent operator of size $m\times m$, where $m\simeq 7\times 10^{5}$, to reveal the dominant gain and modal structures with a thin sketching matrix having as little as $k=10$ columns. The convergence of the gain and resolvent modes will also be reported with respect to the size of the sketching matrix. Influence of the ratio between the first and the second singular values of the resolvent operator will also be examined.

III.1 Problem setup

We consider the spanwise-periodic turbulent flow over a NACA 0012 airfoil at an angle of attack of $9^{\circ}$, a chord-based Reynolds number of $Re_{L_{c}}\equiv v_{\infty}L_{c}/\nu_{\infty}=23,000$ and a free stream Mach number of $M_{\infty}\equiv v_{\infty}/a_{\infty}=0.3$. Here, $v_{\infty}$ is the free-stream velocity, $L_{c}$ is the chord length, $a_{\infty}$ is the free-stream sonic speed, and $\nu_{\infty}$ is the kinematic viscosity. The time- and spanwise-averaged turbulent flow is considered as the base flow for the full and randomized resolvent analyses. For this 2D base flow, we adopt the bi-global setting that decomposes $\pmb{q}^{\prime}$ into spanwise Fourier modes with the wavenumber $k_{z}$.

To obtain the base flow, large-eddy simulation (LES) is performed using a finite-volume compressible flow solver CharLES [34, 35], which is second-order accurate in space and third-order accurate in time. Vremen’s sub-grid scale model [36] is utilized in the LES. The LES is performed on a C-shaped mesh with the domain extent of $x/L_{c}\in[-19,26]$, $y/L_{c}\in[-20,20]$ and $z/L_{c}\in[-0.1,0.1]$ in the streamwise, transverse and spanwise direction, respectively, with the airfoil leading edge at $x/L_{c}=y/L_{c}=0$. Dirichlet boundary condition is specified at the far-field boundary as $(\rho,v_{x},v_{y},v_{z},T)=(\rho_{\infty},v_{\infty},0,0,T_{\infty})$, where $\rho$ is the density, $v_{x}$, $v_{y}$ and $v_{z}$ are respectively the streamwise, transverse and spanwise velocity, and $T$ is the temperature. Over the airfoil, the no-slip adiabatic boundary condition is prescribed. Along the outlet boundary, a sponge layer [37] is applied with a running-averaged state being the target state. The simulation has been validated with respect to the time-averaged pressure, lift and drag over the airfoil. The turbulent separated flow over the airfoil is visualized in figure 2. The visualization of the instantaneous flow shows the laminar separation from the leading edge. We have found that the shear layer physics dominates the pseudospectral behavior of the linearized Navier–Stokes operator, as shear is the main source of nonnormality in the operator. Further details regarding the computational setup, flow physics, and resolvent analysis based flow control of this setup are reported in Yeh and Taira [10].

The full and randomized resolvent analyses are performed on a separate mesh from that used in the LES. This mesh has a 2D rectangular domain with the extent of $x/L_{c}\in[-15,16]$ and $y/L_{c}\in[-12,12]$, comprising approximately $0.15$ million cells. Compared to the LES mesh, the mesh for resolvent analysis is coarser over the airfoil and in the wake, but is much finer in the upstream of the airfoil in order to resolve the forcing mode structures. The time- and spanwise-averaged flow $\bar{\bm{q}}$ obtained from LES is interpolated onto this mesh. At the far-field boundary and over the airfoil, Dirichlet conditions are set for density and velocities and Neumann condition is prescribed for pressure in $\pmb{q}^{\prime}$. At the outlet boundary, Neumann condition is set for all flow variables. With these boundary conditions for $\pmb{q}^{\prime}$ and the base flow $\bar{\bm{q}}$, we construct the linearized Navier–Stokes operator $\bm{L}_{\bar{\bm{q}}}(k_{z})$ for a chosen $k_{z}$. The size of $\bm{L}_{\bar{\bm{q}}}$ and the resolvent operator is approximately $0.75$ million $\times$ $0.75$ million.

For this large operator, we summarize in Table 1 the computational costs of performing resolvent analysis using the Krylov-based Arnoldi-iteration method with a range of parameter setups (i.e., number of singular values ($n_{\text{ev}}$), Krylov subspace dimension ($\text{dim}(\mathcal{S})$), and tolerance) and compare them with those for the present randomized algorithm. The former was conducted by simply calling the svds command in MATLAB. It requires almost $80$ gigabytes of memory and takes approximately $30$ to $70$ minutes (single-core) for each SVD. The high-memory demand necessitates the use of high performance computing resource to conduct the full resolvent analysis. In contrast, the randomized resolvent approach (Algorithm 1) achieves significant reductions in computational time and memory consumption. The present method only requires a third of the memory usage of the Arnoldi-iteration and cuts down the computational time by an order of magnitude. We also note that, in Algorithm 1, the linear systems solvers are the operations with higher computational cost. Since all the three linear systems solvers are conducted for the same operator, the LU decomposition of $\left[-i\omega\bm{I}-\bm{L}_{\bar{\pmb{q}}}\right]$ is performed in the beginning of the algorithm and is passed through the three solvers. This decomposition becomes the main source of the memory consumption.

III.2 Results

We perform the full and randomized resolvent analyses for spanwise wavenumbers of $k_{z}L_{c}=0$ and $20\pi$. Since the base flow is found to be unstable [10], the finite-time approach is adopted with $v_{\infty}/\beta L_{c}=3$ to ensure that the resolvent analysis is performed on a shorter time scale than that associated with the leading growth rate of the instability. Initially, for the randomized analysis, we consider $k=10$ for the width of the test matrix $\bm{\Omega}$ with random Gaussian distribution. Later in Section III.4, we show that this value of $k$ can be further reduced without compromising accuracy.

The leading response and forcing modes obtained from both full and randomized analyses are compared in figure 3 for representative frequencies $St$ and spanwise wavelength $k_{z}L_{c}$. Although $k/m=1.3\times 10^{-5}$, we observe excellent agreement between the modes from the full resolvent analysis and the randomized algorithm. We observe that randomized forcing and response modes are very similar to full resolvent ones. Only at $St=1$ and $k_{z}L_{c}=20\pi$ we observe the appearance of spatially distributed errors in the background, which we refer to as background noise.

The forcing modes are recovered directly from the SVD of the low-dimensional subspace projection. As they are used to recover the response modes using the linear operator, the accuracy of the forcing modes affect the results of the response ones. In the particular case of $St=1$ and $k_{z}L_{c}=20\pi$, when forcing mode is affected by noise, the randomized approach returns some structures emanating from the trailing edge in the response mode $\hat{\pmb{u}}_{1}^{\text{rand}}$, which was not present from the full resolvent analysis. This behavior is related to leakage from high-order modes, as $\sigma_{1}$ and $\sigma_{2}$ are close in the energy spectrum. This remarkable level of agreement over all frequencies and spanwise wavenumbers ensures that the randomized approach presented in Algorithm 1 can help extract insights into the spatial structures to identify regions of sensitivity and guide flow control efforts.

These results were obtained using the present implementation that extends the original randomized SVD algorithm. In figure 4, using the original randomized algorithm [13] within the resolvent analysis to recover the left singular vectors and singular values, the response modes contain background noise. Here, we use Algorithm 1 up to line 5, then SVD is performed as $\bm{B}=\bm{\tilde{U}}\bm{\Sigma}\bm{V}^{*}$ and $\bm{U}$ is recovered a posteriori using the original procedure for the randomized SVD with the resolvent analysis [13, 25], by $\bm{U}=\bm{Q}\bm{\tilde{U}}$. Yeh and Taira [10] showed that the present problem setup presents a peak in the singular values near $St=6$ for both spanwise wavenumbers, influenced by the eigenmodes associated with the shear-layer structure over the separation bubble. These eigenmodes are highly nonnormal and induce high-energy amplification through pseudoresonance [6]. For the frequencies in a narrow region near $St=6$, both implementations present similar results. Far from this band, the response modes obtained by the original procedure [13] are contaminated by random background noise or leakage from higher-order modes, as shown in figure 4 for $St=0.5$ and $15$ and $k_{z}L_{c}=0$. In these critical cases, the original randomized approach may not provide meaningful insights into flow physics. The present implementation shown in Algorithm 1 improves solving for the response modes. Algorithm 1 does not enhance the forcing modes, as they are already accurate. When utilizing this technique to generate reduced-order models, one may perform steps 11 and 12 in Algorithm 1 to orthogonalize the left singular vectors. Considering the results obtained by our implementation, the modes are found very accurately for almost all frequencies and wavenumbers. To provide a concrete assessment, we quantitatively assess the accuracy of the randomized resolvent analysis.

The agreement between the full and randomized analyses with respect to the gain (leading singular value) and modes over a range of frequencies is presented in figure 5. When Algorithm 1 is applied to recover left singular vectors and singular values, the randomized analysis accurately captures the trend of gain distribution over $1\lesssim St\lesssim 15$ in figure 5(a,b) for both wavenumbers. At the low and high-frequency ends, the gain shows deviations. The resemblance of the modal structure is quantified in figure 5(c,d) with the cosine similarities, i.e., the inner products, $\langle\hat{\pmb{u}}_{1}^{\text{full}},\hat{\pmb{u}}_{1}^{\text{rand}}\rangle$ and $\langle\hat{\pmb{v}}_{1}^{\text{full}},\hat{\pmb{v}}_{1}^{\text{rand}}\rangle$. As singular vectors are normalized, the cosine similarity of $1$ suggests that perfect match is attained between the modes from full and randomized resolvent analyses. Since these modes are complex, the cosine similarity removes dependence on the phase difference. For almost the entire range of frequencies the cosine similarities are near unity, which means the agreement between full and randomized modes is excellent. When this value is reduced, the modes may be affected by noise, as seen for $\hat{\pmb{v}}_{1}^{\text{rand}}$ at $St=1$ and $k_{z}L_{c}=0$ in figure 4. For the original approach, using the randomized SVD algorithm [13], the modes have good agreement for a narrow band of frequencies only. By comparing the results from figures 3 and 4 to the values in figure 5(c,d), we observe that the noise affects the modes when cosine similarity is below $0.5$. For frequencies and wavenumbers with cosine similarity up to 0.8 or higher, there is no noise and the results for full and randomized resolvent agree well. For this reason, it is desirable to search for solutions that provide a reliable agreement up to this scale to a broad range of frequencies and both wavenumbers. With the high-gain frequency range well captured, randomized resolvent analysis has demonstrated its capability of predicting the dominant pathway for energy amplification over the spectral space with reduced computational cost.

As stated in section II.2, the use of low-rank approximation in the randomized approach is built upon the assumption of the low-rank nature of the resolvent operator. The randomized resolvent analysis shows its strength when the singular values exhibit fast decay, as evident from figure 6. The accuracy of the modal structured captured by randomized analysis is examined with respect to the ratio of the leading and second singular values, $(\sigma_{1}/\sigma_{2})^{\text{full}}$ from the full resolvent analysis. The error in the modal structures exhibits a decreasing trend as this ratio increases. When this ratio is close to unity, the randomized technique may not accurately separate the first and second modes. In fact, the aforementioned trailing edge structure that appeared in the randomized response mode for $k_{z}L_{c}=20\pi$ and $St=1$ is caused by the leakage of the structures from the second response mode (see figure 3). When the ratio $(\sigma_{1}/\sigma_{2})^{\text{full}}$ is above $30$, the error decreases to $\lesssim 10^{-5}$ for forcing modes and $\lesssim 10^{-8}$ for response modes.

Next, we study the influence of the width of the test matrix $k$ on the error in the leading singular values and modes, as presented in figure 7. When the value of $k$ is varied from $2$ to $500$, the error from the use of randomized analysis decreases. For three representative frequencies, we observe the same rate of convergence $\approx\mathcal{O}(k)$ for both the gain and cosine similarity. As stated in section II.2, increasing $k$ has the same practical effect of oversampling, in the present application. For this flow, we observe that $k=10$ is sufficient to achieve sufficient accuracy with $\lesssim 1\%$ error, which is remarkably low when compared to the high dimensionality of the resolvent operator.

The computational cost of the randomized resolvent technique can be further reduced. For instance, the biconjugate gradient stabilized (BiCGStab) or the generalized minimum residual (GMRES) methods can be utilized to solve linear systems with appropriate preconditioners (e.g., incomplete LU and Jacobian).

III.3 Higher-order modes

Let us discuss the performance of the randomized technique with respect to the high-order modes. For some cases, the second largest singular value may also be spaced apart from the higher-order singular values and be determined accurately. In figure 8, we show for $k_{z}L_{c}=0$ and $St=4$ a case where both the leading and second singular values are spread from the rest of the singular values. In this case, the randomized algorithms accurately capture the second modes. The flow structures at the trailing edge are perceived in the response modes. The forcing modes appear over the pressure side near the trailing edge. For the results obtained from Algorithm 1, the modes are the same for randomized resolvent and for the full resolvent. However, for the resolvent analysis using the randomized SVD algorithm [13] within the resolvent analysis, the secondary singular values are poorly captured and the modes are polluted by background noise and leakage from other modes (not shown here). For this reason, when applying randomized resolvent, the present implementation shown in Algorithm 1 must be considered as they can approximate the detached high-order modes accurately.

III.4 Choice of the test matrix $\bm{\Omega}$

For the randomized resolvent analysis with a test matrix size of $k=10$, Gaussian, orthonormal, Rademacher and ultrasparse Rademacher test matrices provide similar results and no observable difference in computational savings. By using the implementation shown in Algorithm 1, all test matrices present similar accuracy as shown in figures 5 and 6. For very low and very high $St$ numbers in the range of frequencies analyzed in this work, where the singular values decay slowly, one can increase the size $k$ or apply subspace or power iterations. However, it is possible to obtain more accurate results by constructing a random test matrix $\bm{\Omega}$ that incorporates physical insights from the base flow.

While the random test matrix is effective in yielding accurate results, we can consider constructing a test matrix that can generate the entries in a smart manner by incorporating the knowledge of the base flow. We know that regions of strong shear are important in amplifying forcing inputs. Moreover, regions with minimal velocity gradients are not that important. For these reasons, the velocity gradient at each grid point can be used to scale the test matrix. Here, we propose a physics-informed test matrix, $\bm{\Omega}_{p}$, scaled by the 2-norm of the velocity gradient, $||\bm{\nabla}{\bm{v}}||_{2}$, where ${\bm{v}}$ is the velocity vector. We construct a scaling factor $\bm{\Phi}_{j}=||\bm{\nabla}{\bm{v}}_{j}||_{2}$ at each grid point $j$. The scaling vector $\bm{\Phi}$ has to be stacked according to the number of variables to reach the size $m$ of the linear operator. The physics-informed test matrix then becomes

[TABLE]

The results based on the physics-informed test matrix are shown in figure 9. While the results obtained from the use of a normal distribution test matrix have shown excellent accuracy, the use of physics-informed test matrix further improves the accuracy by a few orders of magnitude for the considered frequencies and spanwise wavenumbers. More importantly, we achieve results with comparable or higher level of accuracy using a extremely low width of the test matrix of $k=2$. This results in a considerable reduction in computation time (see Table 1). Using $k=10$, linear systems are solved at least 30 times. Now, using $k=2$, only 6 linear solvers are needed to obtain the same accuracy, which is achieved only with a physics-informed scaling of the test matrix. By combining randomized numerical linear algebra and some physical insights, we are now empowered to perform the input-output analysis for ever more complex 2D and 3D turbulent base flows on a standard computer, or on a high-performance computing cluster to expand the envelope of resolvent analysis.

IV Conclusion

Resolvent analysis has proven to be a powerful technique to reveal the input-output characteristics of fluid flows. However, the computational cost and memory allocation of the resolvent analysis can be taxing for high-Reynolds number flows, making it prohibitive to be applied to complex turbulent base flows. The major computational cost of the analysis is associated with the SVD of the resolvent operator. To remove this bottleneck, the randomized approach has been adopted to reduce the computational cost of SVD by considering the low-rank approximation of the resolvent operator. This was achieved by constructing the low-rank basis based on the insights gained from the sketch of the resolvent, which is obtained from a linear system solver. For flows with fast singular value decays, e.g., flows with strong shear and separation, the randomized resolvent analysis reveals its power to accurately capture the response and forcing modes as well as the gain. Moreover, we consider the use of the velocity gradient to scale the random test matrix. Such scaling enhanced the accuracy of the randomized resolvent analysis. The necessary computation time was significantly reduced as the number of linear systems to be solved is considerably smaller. To demonstrate the capability of the randomized resolvent analysis, we analyzed the turbulent separated flow over a NACA 0012 airfoil. Excellent agreement of the leading forcing and response modes and the gains were shown between the full and randomized resolvent methods. By incorporating the knowledge of the base flow in terms of its velocity gradient into the randomized test matrix, additional speed up and accuracy enhancement were achieved. With the computational cost and memory allocation being relieved with the randomized approach, the application of the resolvent analysis can be significantly extended to higher-Reynolds number 2D and 3D base flows.

Acknowledgments

We gratefully acknowledge the support from the Office of Naval Research (N00014-16-1-2443), Air Force Office of Scientific Research (FA9550-18-1-0040), and Army Research Office (W911NF-17-1-0118). We thank L. Mathelin, S. L. Brunton, and N. B. Erichson for the enlightening discussions. Some of the computations presented here were supported by the Department of Defense High Performance Computing Modernization Program.