Efficient harmonic resolvent analysis via time-stepping

Abstract

We present an extension of the RSVD-$\Delta t$ algorithm initially developed for resolvent analysis of statistically stationary flows to handle harmonic resolvent analysis of time-periodic flows. The harmonic resolvent operator, as proposed by Padovan et al. (2020), characterizes the linearized dynamics of time-periodic flows in the frequency domain, and its singular value decomposition reveals forcing and response modes with optimal energetic gain. However, computing harmonic resolvent modes poses challenges due to $(i)$ the coupling of all $N_{\omega}$ retained frequencies into a single harmonic resolvent operator and $(ii)$ the singularity or near-singularity of the operator, making harmonic resolvent analysis considerably more computationally expensive than a standard resolvent analysis. To overcome these challenges, the RSVD-$\Delta t$ algorithm leverages time stepping of the underlying time-periodic linearized Navier-Stokes operator, which is $N_{\omega}$ times smaller than the harmonic resolvent operator, to compute the action of the harmonic resolvent operator. We develop strategies to minimize the algorithm’s CPU and memory consumption, and our results demonstrate that these costs scale linearly with the problem dimension. We validate the RSVD-$\Delta t$ algorithm by computing modes for a periodically varying Ginzburg-Landau equation and demonstrate its performance using the flow over an airfoil.

1 Introduction

Model reduction plays a critical role in the study of fluid mechanics due to the complex and high-dimensional nature of fluid flows, especially turbulent ones. In particular, modal decompositions, both data-driven and equation-based, have proven to be effective at identifying low-dimensional sets of modes associated with interpretable coherent structures that significantly influence quantities of engineering interest such as kinetic energy, heat and mass transfer, and noise emissions (Taira et al., 2017; Towne et al., 2018). In particular, resolvent analysis has been extensively employed to comprehend, model, and control statistically stationary turbulent flows such as wall-bounded turbulence (Dawson & McKeon, 2019), turbulent boundary layer (Sipp & Marquet, 2013), and jet flows (Pickering et al., 2021). However, real-world fluid systems often exhibit non-stationary behavior, including periodicity, rendering resolvent analysis less effective. Examples of flows exhibiting periodic behavior include internal waves in a stratified fluid (Troy & Koseff, 2005), vortex shedding in the wake of a bluff body (Giannenas et al., 2022), and pulsatile blood flow (Farghadan & Arzani, 2019). Recently, an extension of resolvent analysis, known as harmonic resolvent analysis, has been developed by Padovan et al. (2020). This extension characterizes the perturbation forcing and responses around a periodic base flow, providing a more effective tool for analyzing, modeling, and controlling periodically time-varying flows. The harmonic resolvent operator can be considered a special case of the harmonic transfer function introduced by Wereley (1990), wherein the operator exhibits a dominant frequency.

Investigating the harmonic resolvent operator is one way to identify the optimal spatio-temporal perturbations and their corresponding responses in fluid systems. Unlike the resolvent operator, the harmonic resolvent operator takes into consideration the periodicity of the base flow by representing it in the form of a Fourier series. This allows for characterization of the dynamic behavior of the system, particularly with respect to its spectral components and it provides a means to understand and explore the interactions between input and output modes through the linearized dynamics. Harmonic resolvent analysis determines the triadic interactions between the input frequency $\omega_{1}$, the base flow frequency $\omega_{1}-\omega_{2}$, and the response frequency $\omega_{2}$. The number of Fourier modes of the base flow is determined by its power spectral density (PSD), which in turn determines the number of triads interacting with the forcing. Typically, periodic flows with a dominant periodicity exhibit a multimodal nature with their most energetic modes being harmonics of the base frequency. When the base flow is statistically stationary, linearization around the temporal average is sufficient, and traditional resolvent analysis arises as a special case.

Recent studies have expanded on this research area to model and control the dynamics of turbulent flows with periodic motions. Lin et al. (2023) applied the harmonic resolvent framework for linear time-periodic systems with more than one dominant frequency, offering a tool to study the flow control of a plunging cylinder. In a distinct approach, Franceschini et al. (2022) proposed a novel method aimed at capturing the phase dependency of small-scale turbulent phenomena in flows exhibiting a periodic limit cycle. This method builds upon classical spectral proper orthogonal decomposition (SPOD) and resolvent analyses but introduces a quasi-steady approximation that effectively separates the high-frequency turbulent component from the slow periodic oscillation. Moreover, Heidt & Colonius (2023) have pioneered the theory and algorithm for cyclostationary SPOD (CS-SPOD). This innovative approach extends the conventional SPOD algorithm to effectively capture the characteristics of flows characterized by periodic motions. CS-SPOD is closely tied to harmonic resolvent analysis, akin to the relationship between SPOD and resolvent analysis (Towne et al., 2018).

From a computational standpoint, computing harmonic resolvent modes presents similar challenges to resolvent analysis, including the need to invert sparse matrices in Fourier space. This challenge is further exacerbated in harmonic resolvent analysis due to expanded system size, where all relevant frequencies are merged into a unified operator, in contrast to resolvent analysis, where each frequency has its own independent system. The second shared challenge involves performing the singular value decomposition (SVD). Specifically, turbulent flows characterized by periodic motions require a high-resolution mesh, presenting both challenges as formidable obstacles. A modified randomized SVD (RSVD) can be used to simultaneously overcome both computational challenges (Moarref et al., 2013; Ribeiro et al., 2020; Padovan et al., 2020). Nevertheless, the modified RSVD, referred to as “RSVD-LU” in this study, becomes impractical for tackling large-scale turbulent flows since it requires the lower-upper (LU) decomposition of the harmonic resolvent operator. Recently, RSVD-$\Delta t$, an improvement to RSVD-LU, addresses bottlenecks in resolvent analysis (Farghadan et al., 2023). The approach of time-stepping, initially developed by Monokrousos et al. (2010) and further optimized by Martini et al. (2021), allows for computing the action of an operator on a vector without solving a linear system in Fourier space, making it a key component of RSVD-$\Delta t$.

The main objective of this study is to expand the applicability of the RSVD-$\Delta t$ algorithm, initially crafted for resolvent analysis, to facilitate the computation of harmonic resolvent modes. A streaming approach is adopted for time integration, reducing memory usage while maintaining computational efficiency. Moreover, an efficient transient removal strategy is introduced, cutting the temporal length of the integration and thereby minimizing the CPU cost of RSVD-$\Delta t$ for achieving a desired accuracy. The RSVD-$\Delta t$ algorithm for periodic flows maintains linear scalability in both CPU time and memory consumption. This scalability enables the efficient computation of harmonic resolvent modes for large-scale turbulent flows that were previously out of reach.

The structure of our paper unfolds as follows. We begin by explaining resolvent analysis in $\lx@sectionsign$2, followed by the formulation of harmonic resolvent analysis in $\lx@sectionsign$3, where the computation of harmonic resolvent modes using RSVD-LU is briefly described. The RSVD-$\Delta t$ algorithm is introduced in $\lx@sectionsign$4. A comparison of the computational complexity of algorithms and suggested improvements appear in $\lx@sectionsign$5. Further details on performance and error analysis are presented in $\lx@sectionsign$6. In $\lx@sectionsign$7, we conduct two test cases to demonstrate the accuracy and cost efficiency of RSVD-$\Delta t$ for two periodic systems, and we conclude with final remarks in $\lx@sectionsign$8.

2 Resolvent analysis

Before considering the time-periodic case, we briefly review standard resolvent analysis for stationary flows. The Navier-Stokes equations can be linearized around a base flow via Reynolds decomposition of the state, $\bm{q}=\bar{\bm{q}}+\bm{q}^{\prime}$, where $\bar{\bm{q}}$ represents the temporal average of the flow, and $\bm{q}^{\prime}$ denotes the time-varying fluctuations. The dynamics of perturbations are then described by the linear time-invariant system

$$\frac{d\bm{q}^{\prime}}{dt}=\bm{A}(\bar{\bm{q}})\bm{q}^{\prime}+\bm{f}^{\prime}(\bar{\bm{q}},\bm{q}^{\prime}),$$ (2.1)

where $\bm{A}\in\mathbb{C}^{N\times N}$ is the linearized Navier-Stokes (LNS) operator, $N$ is the state dimension of the discretized system, and $\bm{f}^{\prime}$ represents an external forcing as well as nonlinear terms that are treated as an additional forcing, as suggested by McKeon & Sharma (2010). It should be noted that $\bm{A}$ is time-independent because the base flow is time-invariant. By taking the Fourier transform

$$\mathcal{F}(\cdot)=\hat{(\cdot)}(\omega)=\int_{-\infty}^{+\infty}(\cdot)e^{-i\omega t}\,dt$$ (2.2)

in time, we obtain

$$\hat{\bm{q}}(\omega)=\bm{R}(\omega)\hat{\bm{f}}(\omega),$$ (2.3)

where $\omega$ represents the frequency and $\hat{(\cdot)}$ denotes the frequency-domain representation of the time-domain vector. The resolvent operator

$$\bm{R}=(\text{i}\omega\bm{I}-\bm{A})^{-1}\in\mathbb{C}^{N\times N}$$ (2.4)

is a transfer function that maps the input forcing to the output response. Here, $\text{i}=\sqrt{-1}$ and $\bm{I}$ is the identity matrix of appropriate dimension.

To analyze the perturbation characteristics of a system, the transfer function can be examined across a range of frequencies. Typically, the most amplified responses are identified by computing the left singular vectors $\textbf{U}_{R}$ of the resolvent operator,

$$\bm{R}=\textbf{U}_{R}\bm{\Sigma}_{R}\textbf{V}_{R}^{*},$$ (2.5)

where $(\cdot)^{*}$ denotes complex conjugate transpose. The optimal forcing inputs that lead to the highest amplifications comprise the right singular vectors $\textbf{V}_{R}$, and the degree of amplifications are determined by the singular values $\bm{\Sigma}_{R}$.

The SVD in 2.5 implies that disturbance amplitudes are measured in a Euclidean norm. More generally, $||\bm{x}||^{2}_{f}=\langle\bm{x},\bm{x}\rangle_{f}=\bm{x}^{*}\bm{W}_{f}\bm{x}$ computes the $f$-norm of any vector $\bm{x}$, and input and output norms can be different, i.e., $||\cdot||_{q}=||\cdot||_{f}$ is not necessary. Additionally, one can define an input matrix $\bm{B}$ to restrict the forcing and an output matrix $\bm{C}$ to extract the output of interest from the state. For the sake of notational brevity, we assumed the weight, input, and output matrices to be the identity, i.e., $\bm{W}=\bm{B}=\bm{C}=\bm{I}$, throughout this paper. However, our algorithm, with minor adjustments, can accommodate non-identity weight, input, and output matrices by computing the modes of the modified resolvent operator

$$\tilde{\bm{R}}=\bm{W}^{1/2}_{q}\bm{C}\bm{R}\bm{B}\bm{W}^{-1/2}_{f}.$$ (2.6)

For more details, please refer to Farghadan et al. (2023).

3 Harmonic resolvent analysis

In this section, we provide an overview of harmonic resolvent analysis. After reviewing its basic formulation, we discuss the singularity issue that can arise in some special flows, the truncation of the harmonic resolvent operator, and two variations of harmonic resolvent analysis, namely cross-frequency amplification and subharmonic resolvent. Lastly, a brief overview is provided for the state-of-the-art algorithm used to compute the harmonic resolvent modes and gains.

3.1 Formulation

Analogous to resolvent analysis, by inserting the generalized Reynolds decomposition $\bm{q}(t)=\bar{\bm{q}}(t)+\bm{q}^{\prime}(t)$ into the Navier-Stokes equations, the linearization around a periodic base flow $\bar{\bm{q}}(t)=\bar{\bm{q}}(t+T)$ leads to the linear time-periodic system

$$\frac{d\bm{q}^{\prime}}{dt}=\bm{A}_{p}\bm{q}^{\prime}+\bm{f},$$ (3.1)

where $\bm{A}_{p}(t)=\bm{A}_{p}(t+T)\in\mathbb{C}^{N\times N}$ is the periodic LNS operator, $\omega_{f}$ is the fundamental frequency, and $T=2\pi/\omega_{f}$ is defined as the fundamental cycle length. As both $\bm{A}_{p}$ and $\bm{q}^{\prime}$ are time-periodic, the Fourier transform of (3.1) incorporates interactions between all pairs of frequencies. For simplicity, we drop the prime notation for fluctuations from here on out.

Expanding $\bm{A}_{p}(t)$ and $\bm{q}(t)$ in terms of the Fourier series

$$\begin{gathered}\bm{A}_{p}(t)=\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t},\\ \bm{q}(t)=\sum_{j=-\infty}^{\infty}\hat{\bm{q}}_{j}e^{\text{i}j\omega_{f}t},\end{gathered}$$ (3.2)

where $(\cdot)_{j}$ denotes $j^{th}$ harmonic of $\omega_{f}$, and substituting these expansions into (3.1), we obtain

$$\left[\bm{T}\hat{\bm{q}}\right]_{k}=\text{i}k\omega_{f}\hat{\bm{q}}_{k}-\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{k-j}\hat{\bm{q}}_{j}=\hat{\bm{f}}_{k}.$$ (3.3)

Here, $\bm{T}$ is an infinite dimensional block matrix of the form

$$\bm{T}=\begin{bmatrix}\ddots&\vdots&\vdots&\vdots&\vdots&\udots\\ \dots&\bm{R}_{-1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&\hat{\bm{A}}_{-3}&\dots\\ \dots&\hat{\bm{A}}_{1}&\bm{R}_{0}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&\dots\\ \dots&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{1}^{-1}&\hat{\bm{A}}_{-1}&\dots\\ \dots&\hat{\bm{A}}_{3}&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{2}^{-1}&\dots\\ \udots&\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},$$ (3.4)

where the diagonal entries are in fact the inverse of resolvent operators at various frequencies,

$$\bm{R}_{j}^{-1}=\text{i}j\omega_{f}\bm{I}-\hat{\bm{A}}_{p,0},\forall j\in\mathbb{Z}.$$ (3.5)

The harmonic resolvent operator

$$\bm{H}=\bm{T}^{-1}$$ (3.6)

transfers the harmonic inputs to the outputs for periodic base flows,

$$\hat{\bm{q}}=\bm{H}\hat{\bm{f}}.$$ (3.7)

The SVD of the harmonic resolvent operator

$$\bm{H}=\textbf{U}_{H}\bm{\Sigma}_{H}\textbf{V}_{H}^{*}$$ (3.8)

unveils the most amplified responses $\textbf{U}_{H}$ that correspond to optimal forcing $\textbf{V}_{H}$, each accompanied by associated amplification magnitudes $\bm{\Sigma}_{H}$. The matrix $\bm{H}$ encompasses the base flow, the fundamental frequency, and higher harmonics. Whereas individual frequencies can be analyzed separately in a standard resolvent analysis, all frequencies are coupled within $\bm{H}$ in harmonic resolvent analysis and must be considered simultaneously. Furthermore, the modes of the modified harmonic resolvent operator, defined as

$$\tilde{\bm{H}}=\bm{W}^{1/2}_{q}\bm{C}\bm{H}\bm{B}\bm{W}^{-1/2}_{f},$$ (3.9)

can be computed using our algorithm with the same minor adjustments as detailed in Farghadan et al. (2023). In a special case where the base flow is time-independent, i.e., $\hat{\bm{A}}_{p,j}=0$ for $j\neq 0$, $\bm{H}$ reduces to

$$\bm{H}=\begin{bmatrix}\ddots&\vdots&\vdots&\vdots&\udots\\ \dots&\bm{R}_{-1}&0&0&\dots\\ \dots&0&\bm{R}_{0}&0&\dots\\ \dots&0&0&\bm{R}_{1}&\dots\\ \udots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},$$ (3.10)

so each frequency is decoupled. Therefore, resolvent analysis can be perceived as a special case of harmonic resolvent analysis.

3.2 Singular or near-singular nature of $\bm{T}$

The operator $\bm{T}$ tends to be ill-conditioned for many periodic systems (Padovan et al., 2020). A system is considered ill-conditioned when the ratio of the largest and smallest singular values is large $(\kappa=\sigma_{\text{max}}/\sigma_{\text{min}}\sim O(10^{4})$ or higher). The smallest singular value of $T$ can be shown to reach machine precision zero if the periodic base flow satisfies the Navier-Stokes equations (Padovan et al., 2020; Leclercq & Sipp, 2023), i.e., when

$$\frac{d\bm{\bar{q}}}{dt}=\bm{\mathcal{N}}(\bm{\bar{q}}),$$ (3.11)

where $\bm{\mathcal{N}}$ is the nonlinear Navier-Stokes operator. Taking a time derivative of both sides,

$$\frac{d\left(\frac{d\bm{\bar{q}}}{dt}\right)}{dt}=\frac{d\left(\bm{\mathcal{N}}(\bm{\bar{q}})\right)}{dt}=\frac{\partial\bm{\mathcal{N}}}{\partial\bm{\bar{q}}}\frac{d\bm{\bar{q}}}{dt}=\bm{A}_{p}\frac{d\bm{\bar{q}}}{dt}.$$ (3.12)

Taking a Fourier transform of (3.12) and recalling that $\bm{T}$ is obtained as shown in (3.3), $\widehat{d\bm{\bar{q}}/dt}$ is a non-trivial solution of

$$\bm{T}\bm{q}=0,$$ (3.13)

residing within the null space of $\bm{T}$, thereby proving the singularity of $\bm{T}$. In systems where the base flow approximately satisfies (3.11), $\bm{T}$ becomes nearly singular. Consequently, irrespective of the specific periodic system under consideration, $\bm{T}$ might be poorly conditioned, presenting even more challenges when dealing with large systems.

Padovan et al. (2020) demonstrated that defining the harmonic resolvent operator in a manner that eliminates the phase shift imparted by the Fourier coefficients of $d\bm{\bar{q}}/dt$ is effective, without adversely affecting the other dominant amplification mechanisms. By constraining $\bm{T}$ to a subspace $\bm{U}^{\perp}$ that is orthogonal to the direction of the phase shift $\widehat{d\bm{\bar{q}}/dt}$, the range of the harmonic resolvent operator is limited to $\bm{U}^{\perp}$. The key concept here is to eliminate the singular vectors associated with the smallest singular value of $\bm{T}$ without affecting the other singular values and vectors of $\bm{T}$. Given that the phase shift resides in the null space of $\bm{T}$, the desired right singular vector $\bm{v}=\widehat{\frac{d\bm{\bar{q}}}{dt}}/||\widehat{\frac{d\bm{\bar{q}}}{dt}}||$ is readily available. We can also solve

$$\bm{T}^{*}\bm{\tilde{u}}=\bm{v}$$ (3.14)

and compute the corresponding left singular vector as $\bm{u}=\bm{\tilde{u}}/||\bm{\tilde{u}}||$. Upon obtaining $\bm{v}$ and $\bm{u}$, the process becomes straightforward: projecting out these modes from the forcing and response terms effectively restricts $\bm{T}$ to $\bm{U}^{\perp}$. These steps are elucidated in detail in Padovan et al. (2020).

To compute the action of the harmonic resolvent operator on a given vector (or matrix), our goal is to determine the vector (or matrix) $\bm{q}$ that satisfies $\bm{q}=\bm{H}\bm{f}$. In brief, to compute the action of a singular $\bm{H}$ on a vector $\tilde{\bm{f}}$, the component along $\bm{u}$ must be projected out as

$$\tilde{\bm{f}}_{in}=\tilde{\bm{f}}-\bm{u}(\bm{u}^{*}\tilde{\bm{f}}),$$ (3.15)

and this will provide a response

$$\bm{T}\bm{\tilde{q}}=\tilde{\bm{f}}_{in}\rightarrow\bm{\tilde{q}}=\bm{H}\tilde{\bm{f}}_{in},$$ (3.16)

which is orthogonal to the direction of the phase shift. Nevertheless, to mitigate potential round-off errors, it is advisable to refine the output by projecting out the component along $\bm{v}$ as

$$\tilde{\bm{q}}_{out}=\bm{\tilde{q}}-\bm{v}(\bm{v}^{*}\bm{\tilde{q}}).$$ (3.17)

Similarly, when computing the action of singular $\bm{H}^{*}$ on a vector, the initial step involves projecting out the component along $\bm{v}$. After the response is obtained, to mitigate round-off errors, the component along $\bm{u}$ is projected out.

3.3 Truncating the harmonic resolvent operator

Computing harmonic resolvent modes is contingent upon $\bm{T}$ being of finite size. The size and block sparsity pattern of $\bm{T}$ are determined by two key variables: the number of frequencies within the base flow, $N_{b}$, and within the response, $N_{\omega}$. Here, we elucidate the impact of each variable on the operator.

The rows of $\bm{T}$ consist of $N_{b}$-stencil blocks, dictated by the presence of non-zero elements in $\hat{\bm{A}}_{p,j}$. In the context of periodic flows, the time-dependent base flow can be expressed using a Fourier series expansion

$$\bar{\bm{q}}(t)=\sum_{j\omega_{f}\in\Omega_{\bar{q}}}\hat{\bar{\bm{q}}}_{j}e^{\text{i}j\omega_{f}t}.$$ (3.18)

The set $\Omega_{\bar{q}}$ encompasses the relevant frequencies associated with the dominant flow structures. Typically, $\Omega_{\bar{q}}$ is a subset of $j\omega_{f}$ with $j\in\mathbb{Z}$ and encapsulates the majority of energy within the periodic flows. By retaining the base periodicity and a limited number of higher harmonics, it becomes possible to simplify the base flow representation, resulting in $N_{b}$ frequencies within $\Omega_{\bar{q}}$ and $N_{b}$ non-zero elements within each row of $\bm{T}$.

On the other hand, we aim to resolve perturbations occurring at temporal frequencies $\omega\in\Omega_{q}$, where $\Omega_{q}$ represents a subset of $j\omega_{f}$ with $j\in\mathbb{Z}$, although the specific choice of harmonics for $\Omega_{q}$ may vary, as further elucidated in $\lx@sectionsign$3.5. Typically, $\Omega_{\bar{q}}$ is a subset of $\Omega_{q}$, i.e., the perturbation frequency content spans a wider range of harmonics (Padovan et al., 2020). The number of frequencies $N_{\omega}$ expanding both $\bm{f}$ and $\bm{q}$ determines the dimensions of the block matrices, yielding a size of $(N_{\omega}N)\times(N_{\omega}N)$ for $\bm{T}$.

We provide a simple example to visually observe the harmonic resolvent operator. Suppose we use $\Omega_{\bar{q}}=\{-2,-1,0,1,2\}\omega_{f}$ to expand the LNS operator and $\Omega_{q}=\{-4,-3,\dots,3,4\}\omega_{f}$ to expand the perturbations. In this case, the stencil length would be $N_{b}=5$, and the size of the block matrix would be $(9N)\times(9N)$. The harmonic resolvent operator is then

$$\bm{H}=\begin{bmatrix}\bm{R}_{-4}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0&0&0\\ \hat{\bm{A}}_{1}&\bm{R}_{-3}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0&0\\ \hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{-2}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0\\ 0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{-1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0\\ 0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{0}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0\\ 0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0\\ 0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{2}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}\\ 0&0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{3}^{-1}&\hat{\bm{A}}_{-1}\\ 0&0&0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{4}^{-1}\end{bmatrix}^{-1}.$$ (3.19)

The forcing and response vectors can be represented as

$$\hat{\bm{f}}=\begin{bmatrix}\hat{\bm{f}}_{-4}\\ \hat{\bm{f}}_{-3}\\ \hat{\bm{f}}_{-2}\\ \hat{\bm{f}}_{-1}\\ \hat{\bm{f}}_{0}\\ \hat{\bm{f}}_{1}\\ \hat{\bm{f}}_{2}\\ \hat{\bm{f}}_{3}\\ \hat{\bm{f}}_{4}\\ \end{bmatrix},\hat{\bm{q}}=\begin{bmatrix}\hat{\bm{q}}_{-4}\\ \hat{\bm{q}}_{-3}\\ \hat{\bm{q}}_{-2}\\ \hat{\bm{q}}_{-1}\\ \hat{\bm{q}}_{0}\\ \hat{\bm{q}}_{1}\\ \hat{\bm{q}}_{2}\\ \hat{\bm{q}}_{3}\\ \hat{\bm{q}}_{4}\\ \end{bmatrix}.$$ (3.20)

The accuracy of the obtained resolvent system in Fourier space depends on how well the chosen frequencies capture the relevant flow information, and it is affected by the truncation limit imposed on the frequency range. Although it is feasible to reduce the infinite-dimensional harmonic operator to a finite dimension based on perturbation frequencies, the choice of frequency range is consequential. To demonstrate the influence of truncation on response modes, consider a periodic flow with a dominant first harmonic, such that the base flow can be expressed as

$$\bar{\bm{q}}(t)=\hat{\bar{\bm{q}}}_{0}+\hat{\bar{\bm{q}}}_{1}e^{\text{i}\omega t}+\hat{\bar{\bm{q}}}_{-1}e^{-\text{i}\omega t}.$$ (3.21)

Assume that the LNS equations are subjected to a constant forcing $\bm{f}(t)=\hat{\bm{f}}_{0}$. The LNS operator exhibits frequency content at 0 and $\pm\omega$, indicating that a constant forcing will directly induce a response at these frequencies. However, higher harmonics are also indirectly stimulated, and the triggering of responses at higher harmonics follows a cascade pattern, with norms gradually decaying towards zero, i.e., $\lim_{j\to\infty}||\hat{\bm{q}}_{j}||=0$. Hence, the truncation must be set at a sufficiently high level to preserve the response norms.

In practice, by defining $\Omega_{q}$ for both input perturbation and the response, we compute accurate harmonic resolvent modes as long as the norms of higher frequency modes become less relevant within the spectrum. In other words, if we set both forcing and response frequencies sufficiently high and, upon computing harmonic resolvent modes, find that the norms of input and output modes with higher frequency content are relatively small, it indicates convergence, and including higher frequencies has negligible impact on the results. In a case, for instance, where response modes with higher frequency content remain relatively important, extending the frequency range of the output modes is necessary without changing the input frequency content. In our algorithm, the input and output frequency ranges are adjustable as needed.

3.4 Cross-frequency harmonic resolvent analysis

Harmonic resolvent analysis captures the interaction between various frequencies within $\Omega_{\bar{q}}$ and $\Omega_{q}$, and the typical objective is to identify the optimal forcing mode, potentially spanning various frequencies, that produces the optimal response across a range of frequencies. Alternatively, one can compute the unit-norm forcing at one frequency $\omega_{1}\in\Omega_{q}$ that triggers the most amplified response at another frequency $\omega_{2}\in\Omega_{q}$. This is achieved by defining

$$\bm{H}_{\omega_{2},\omega_{1}}=\bm{C}\bm{H}\bm{B}$$ (3.22)

and computing the modes and gains following the same procedure as before. Here, matrices $\bm{B}$ and $\bm{C}$ are constructed to extract the $\omega_{1}$ and $\omega_{2}$ frequencies from the forcing and response modes, respectively. This represents a special case of the modified harmonic resolvent operator in (2.6).

3.5 Subharmonic resolvent analysis

The frequency content of the forcing inputs is typically composed of harmonics of the base flow frequency. However, Padovan & Rowley (2022) have shown that the harmonic resolvent system may exhibit sensitivity to subharmonic inputs. For instance, in the case where the fundamental base flow frequency is denoted by $\omega_{f}$, the linearized system may be sensitive to forcing inputs with $\gamma<\omega_{f}$, such as $\omega_{f}/2$. More generally, if the frequency content of a periodic LNS operator is limited to the frequency range $\Omega_{\bar{q}}$, it may be desirable to compute input-output modes with frequencies outside of this range, i.e., , $\gamma\notin\Omega_{\bar{q}}$. In such cases, the subharmonic resolvent modes become relevant.

To identify unique subharmonic inputs, we can take advantage of the linear nature of the harmonic resolvent operator, which allows interaction with the fundamental frequency $\omega_{f}$ and its harmonics. It can be demonstrated that the interval $(-\omega_{f}/2,\omega_{f}/2]$ encompasses all sets of subharmonic frequencies where the input-output modes are unique. As discussed in $\lx@sectionsign$3, a forcing with a frequency $\gamma\in(-\omega_{f}/2,\omega_{f}/2]$ can only trigger responses at frequencies $\gamma+j\omega_{f}$, where $j$ is an integer. Therefore, the input and output modes exist within the frequency range $\Omega_{\gamma}=\gamma\oplus\Omega_{\bar{q}}$, where $\oplus$ denotes element-wise addition. Harmonics of $\gamma$ trigger a unique set of outputs as long as $j\gamma\leq|\omega_{f}|$. Due to the linear relationship, we can show that all harmonics are isolated sets that need to be studied separately. Note that $-\omega_{f}/2$ is excluded since $\Omega_{-\omega_{f}/2}=\Omega_{\omega_{f}/2}$, resulting in redundancy. For $\gamma\in\Omega_{\gamma}$, one can derive the subharmonic resolvent system as (Padovan & Rowley, 2022)

$$\hat{\bm{q}}=(\text{i}\gamma\bm{I}-\bm{T})^{-1}\hat{\bm{f}}.$$ (3.23)

To further elucidate the point, consider an example where the subharmonic frequency of interest is $\omega_{f}/5$ along with its harmonics. Note that the valid harmonics are the ones that fulfill the condition $j\omega_{f}/5\in(-\omega_{f}/2,\omega_{f}/2]$, which includes the set $\{-2/5,-1/5,0,1/5,2/5\}$. While $\Omega_{0}=0\oplus\{-N_{b}/2,\ldots,N_{b}/2\}$, where $N_{b}$ is the number of active frequencies within the base flow, is identical to harmonic inputs, the sets $\Omega_{-2/5}$, $\Omega_{-1/5}$, $\Omega_{1/5}$, and $\Omega_{2/5}$ represent distinct input-output modes. In general, for any two distinct subharmonic frequencies $\gamma_{1}$ and $\gamma_{2}$ within the interval $(-\omega_{f}/2,\omega_{f}/2]$, the corresponding subharmonic modes in $\Omega_{\gamma_{1}}$ are decoupled from those in $\Omega_{\gamma_{2}}$ as there exists no common frequency $\gamma_{1}-\gamma_{2}\in\Omega_{\bar{q}}$. While our algorithm described in $\lx@sectionsign$4 is described for harmonic resolvent analysis, it can be easily extended to compute subharmonic resolvent modes, as described in Appendix A.

3.6 Computing harmonic resolvent modes using RSVD-LU

The RSVD algorithm (Halko et al., 2011) is a randomized linear algebra technique developed to efficiently identify the singular vectors with the highest gains in a given matrix. In the context of harmonic resolvent analysis, the matrix of interest is the harmonic resolvent operator. This operator, similar to the resolvent operator, relies on the inverse of a sparse matrix, as shown in (3.6), making the computation of modes challenging. Modifying the original RSVD algorithm is necessary to address these challenges, enabling harmonic resolvent analysis. The modifications for resolvent analysis have been extensively documented in the literature by Ribeiro et al. (2020) and Farghadan et al. (2023), while the outline for the harmonic resolvent operator can be found in the appendix of Padovan et al. (2020).

In brief, the computation involves computing the action of $\bm{H}$ and $\bm{H}^{*}$ for sketching the range and image of the harmonic resolvent operator, respectively, and approximating the leading harmonic modes. Computing actions of $\bm{H}$ and $\bm{H}^{*}$, however, requires solving linear systems

$$\begin{gathered}\bm{T}\hat{\bm{q}}=\hat{\bm{f}},\\ \bm{T}^{*}\hat{\bm{q}}=\hat{\bm{f}},\end{gathered}$$ (3.24)

respectively, in Fourier space. These linear systems are solved via LU-decomposition of $\bm{T}$, but due to its large size and ill-conditioned nature, computing its LU decomposition can be a formidable computational obstacle, particularly for flows with three inhomogeneous directions. We propose an alternative way in the following section.

4 Computing harmonic resolvent modes using time stepping

In this section, we show how time stepping can be used to efficiently compute harmonic resolvent modes. The same concept was first used by Monokrousos et al. (2010) in the context of the resolvent analysis and further optimized by Martini et al. (2021) and Farghadan et al. (2023). Here, we introduce the extension of our algorithm to compute the harmonic resolvent modes and gains.

4.1 Computing the action of $\bm{H}$ using time stepping

Figure 1: Flowchart depicting the action of $\bm{H}$ on a forcing comprised of $N_{\omega}$ frequencies within the RSVD-LU (top route) and the RSVD-$\Delta t$ (bottom route) algorithms. Both routes produce the same result, but the bottom route is computationally advantageous for large systems.

Computing the action of $\bm{H}$ in the Fourier space poses a bottleneck within the RSVD-LU algorithm. To overcome this limitation, we propose to use time-stepping as an effective surrogate approach. Given that $\Omega_{\bar{q}}$ represents the frequency content of the base flow around which $\bm{T}$ is constructed, we define the time-domain linearized operator as

$$\bm{A}_{p}(t)=\sum_{j\omega_{f}\in\Omega_{\bar{q}}}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}.$$ (4.1)

The action of $\bm{H}$ on $\hat{\bm{f}}$ in Fourier space is expressed in (3.7), where the frequency content of the forcing lies within $\Omega_{q}$.

The steady-state solution of (3.1) is given by

$$\hat{\bm{q}}_{s}=\bm{H}\hat{\bm{f}},$$ (4.2)

when subjected to the forcing represented as

$$\bm{f}=\sum_{j\omega_{f}\in\Omega_{q}}\hat{\bm{f}}_{j}e^{\text{i}j\omega_{f}t}.$$ (4.3)

The steady-state response can be obtained through time-stepping before taking a Fourier transform. The time-domain forcing is constructed to include frequencies identical to those present in $\hat{\bm{f}}$. Since the time-domain and frequency-domain frequencies are the same, both approaches yield identical results, i.e., $\hat{\bm{q}}_{s}=\hat{\bm{q}}$. In practice, the time-stepping method is inherently discrete, and the frequencies are retained up to the maximum limit determined by the chosen time step (below the Nyquist frequency (Nyquist, 1928)). By employing a suitably small time step, we can capture frequencies up to the desired limit.

The equivalence between computing the action of $\bm{H}$ in both the RSVD-LU and RSVD-$\Delta t$ algorithms is depicted in figure 1. In the upper route, within the RSVD-LU algorithm, the LNS equations are first transformed into Fourier space before solving a coupled linear system in the frequency domain for a given forcing. In the lower route, within the RSVD-$\Delta t$ algorithm, the LNS equations are initially integrated in the time domain before being transformed into Fourier space. Both routes produce identical outputs so long as numerical artifacts are minimized.

4.2 RSVD-$\Delta t$ for harmonic resolvent analysis

By recognizing the equivalence between time stepping and solving linear systems in the frequency domain, we can effectively address the limitations of the RSVD-LU algorithm, leading to the development of RSVD-$\Delta t$. The algorithmic steps for extending RSVD-$\Delta t$ to harmonic systems are presented in Algorithm 1.

Algorithm 1 RSVD-$\Delta t$ for harmonic resolvent analysis

1:Input parameters: $\bm{A}_{p},k,q,\Omega_{q},\Omega_{\bar{q}},$TSS, $dt,T_{t}$

2:$\hat{\bm{\varTheta}}\leftarrow\;$randn$(NN_{\omega},k)$ $\triangleright$ Create random test matrices

3:$\hat{\bm{Y}}\leftarrow\;$DirectAction$(\bm{A}_{p},\hat{\bm{\varTheta}},$TSS, $dt,T_{t})$ $\triangleright$ Sample the range of $\bm{H}$

4:if $q>0$ then $\triangleright$ Optional power iteration

5:     $\hat{\bm{Y}}\leftarrow\;$PI$(\bm{A}_{p},\hat{\bm{Y}},q,$TSS, $dt,T_{t})$ $\triangleright$ Power iteration with time-stepping

6:$\hat{\bm{Q}}_{\varOmega}\leftarrow\;$qr$(\hat{\bm{Y}}_{\varOmega})$ $\triangleright$ Build the orthonormal subspace $\hat{\bm{Q}}_{\varOmega}$

7:$\hat{\bm{S}}\leftarrow\;$ AdjointAction$(\bm{A}_{p}^{*},\hat{\bm{Q}},$TSS, $dt,T_{t})$ $\triangleright$ Sample the image of $\bm{H}$

8:$(\tilde{\bm{U}}_{H},\bm{\varSigma}_{H},\bm{V}_{H})\leftarrow\;$svd$(\hat{\bm{S}}_{H})$ $\triangleright$ Obtain $\bm{\varSigma}_{H},\bm{V}_{H}$

9:$\bm{U}_{H}\leftarrow\;\hat{\bm{Q}}_{\varOmega}\tilde{\bm{U}}_{H}$ $\triangleright$ Recover $\bm{U}_{H}$

10:Output parameters: $\bm{U}_{H},\bm{\varSigma}_{H},\bm{V}_{H}$

Algorithm 1: $k,q,\Omega_{q},\Omega_{\bar{q}}$ are common parameters with RSVD-LU. $(\cdot)_{\varOmega}$ indicates all frequencies are merged into a single column, and TSS is an abbreviation for time-stepping schemes (e.g., backward Euler). DirectAction and AdjointAction are functions that solve the direct and adjoint LNS equations, respectively, with a predefined forcing.

The algorithm proposed in this study is based on the RSVD-$\Delta t$ algorithm introduced by Farghadan et al. (2023), and we will provide a concise overview of the key steps. Initially, a random matrix $\hat{\bm{\varTheta}}$ is generated (line 2) to compute the sketch of $\bm{H}$ through time stepping (line 3). DirectAction is a function that computes the action of $\bm{H}$ onto a given forcing using time-stepping and transforms the steady-state responses to Fourier space. Next, an optional power iteration is performed (lines 4 and 5) via $q$ successive applications of DirectAction and AdjointAction, enhancing the accuracy of harmonic resolvent modes. QR decomposition is performed to obtain $\hat{\bm{Q}}_{\varOmega}$ (line 6) before constructing $\bm{S}$ via time stepping (line 7). AdjointAction is the function that computes the action of $\bm{H}^{*}$ onto a given forcing using time-stepping and transforms the steady-state responses to Fourier space. An SVD (line 8) is conducted to obtain the optimal forcing modes $\bm{V}_{H}\in\mathbb{C}^{NN_{\omega}\times k}$ and gains $\bm{\varSigma}_{H}\in\mathbb{R}^{k\times k}$ of $\bm{H}$. Lastly, the optimal response modes $\bm{U}_{H}\in\mathbb{C}^{NN_{\omega}\times k}$ are recovered in line 9.

Compared to the RSVD-$\Delta t$ algorithm for resolvent analysis, two notable differences are present. In the extension of RSVD-$\Delta t$, the response modes at various frequencies are merged prior to both the QR decomposition and the SVD steps. This deviation from the original algorithm is motivated by the nature of the problem that requires interactions between frequencies, unlike in resolvent analysis where the QR decomposition and SVD are performed separately for each individual frequency of interest. The second difference is the generation of LNS operators over time, whereas in resolvent analysis, the LNS operator remains constant throughout the integration.

If $\bm{T}$ approaches machine precision singularity, it is necessary to project out the forcing along the $\bm{u}$ direction before the DirectAction and along $\bm{v}$ before the AdjointAction, as described in $\lx@sectionsign$3.2. Moreover, to mitigate numerical artifacts, it is advisable to project out the response along $\bm{v}$ and $\bm{u}$ after the DirectAction and AdjointAction, respectively.

5 Computational complexity

A brief comparison is conducted between memory requirements and CPU cost of the RSVD-LU and RSVD-$\Delta t$ algorithms. We also propose strategies to minimize the memory consumption. For more detailed information on the CPU and memory scaling of the LU decomposition of the resolvent operator, we refer readers to Farghadan et al. (2023). To offer empirical insights and confirmation of the theoretical scalings discussed below, we present a two-dimensional test case in $\lx@sectionsign$7.

5.1 Memory benefits of RSVD-$\Delta t$

In the case of resolvent analysis, the memory scaling of the LU decomposition is empirically $O(N^{1.2})$ and $O(N^{1.6})$ for two- and three-dimensional systems, respectively (Towne et al., 2022; Farghadan et al., 2023). The memory scaling of LU decomposition is expected to be worse for the harmonic resolvent operator, as the lower- and upper-triangular matrices are denser. When $\bm{T}$ is singular or nearly singular, the computational task becomes even more challenging, as discussed in Padovan et al. (2020).

On the other hand, the memory usage of RSVD-$\Delta t$ for harmonic resolvent analysis is contingent on the sizes of various matrices, specifically the sparse LNS operators $\hat{\bm{A}}_{p}\in\mathbb{C}^{N\times N\times N_{b}}$, and the dense forcing matrix $\hat{\bm{F}}\in\mathbb{C}^{NN_{\omega}\times k}$ and response matrix $\hat{\bm{Q}}\in\mathbb{C}^{NN_{\omega}\times k}$, all of which are stored in Fourier space. The scaling for resolvent analysis is $O(NN_{\omega})$, and the sole extra space is for the storage of LNS operators.

One approach to generating LNS operators for all time points involves retaining $N_{b}$ coefficients in Fourier space and subsequently constructing $\bm{A}_{\varOmega}=\{\bm{A}_{p,1},\bm{A}_{p,2},\bm{A}_{p,3},\dots,\bm{A}_{p,N_{t}}\}$ once for all time steps within a period. However, this approach can be memory-intensive for large systems, especially when $N_{t}\sim O(10^{3}-10^{5})$. An alternative strategy to reduce memory consumption is to create LNS operators on the fly, as shown in figure 2. As a result, no time-domain LNS operator is permanently stored in memory. Each LNS operator is created on-the-fly using the discrete Fourier transform (DFT) of $N_{b}$ Fourier coefficients,

$$\bm{A}_{p}(t)=\sum_{j=-N_{b}/2}^{N_{b}/2}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}.$$ (5.1)

Typically, $\bm{A}_{p}$ is very sparse, with $\texttt{nnz}(\bm{A}_{p})\sim O(N)$, where $\texttt{nnz}(\cdot)$ represents the count of non-zero elements, resulting in an overall memory consumption of $O(NN_{b})$. The methods employed for generating LNS operators, both in the streaming and memory-intensive approaches, mirror those used for creating forcing terms in resolvent analysis, as elaborated in Farghadan et al. (2023). Additionally, the CPU cost of this procedure scales as $O(\texttt{nnz}(\bm{A}_{p})N_{b})$ or $O(NN_{b})$.

Figure 2: Schematic of the action of $\bm{H}$ with streaming DFT/iDFT to transform between the Fourier and time domains.

Streaming DFT is indeed effective in reducing memory consumption; nevertheless, additional memory savings can be specifically achieved for real-valued LNS operators. In many cases, $\bm{A}_{p}$ is real-valued, i.e., confined to $\mathbb{R}^{N\times N}$, yielding $\hat{\bm{A}}_{p,j}$ = $\bar{\hat{\bm{A}}}_{p,-j}$, where $\bar{(\cdot)}$ denotes the complex conjugate. Hence, $\hat{\bm{A}}_{p,j}$ with $j\geq 0$ suffices to construct $\bm{A}_{p}$ as

$$\bm{A}_{p}(t)=\hat{\bm{A}}_{p,0}+2\mathcal{R}e\left(\sum_{j=1}^{N_{b}/2}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}\right),$$ (5.2)

resulting in halving the number of Fourier coefficients of the LNS operator to $\lfloor N_{b}/2\rfloor+1$. Here, $\mathcal{R}e$ represents the real part.

Another opportunity to reduce memory usage for real-valued LNS matrices arises from the symmetry between positive and negative frequencies. Rewriting (3.3) for $-k\omega_{f}$ yields

$$\text{i}(-k\omega_{f})\hat{\bm{q}}_{-k}-\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{-k+j}\hat{\bm{q}}_{-j}=\hat{\bm{f}}_{-k},$$ (5.3)

and the complex conjugate version of (3.3) becomes

$$-\text{i}k\omega_{f}\bar{\hat{\bm{q}}}_{k}-\sum_{j=-\infty}^{\infty}\bar{\hat{\bm{A}}}_{k-j}\bar{\hat{\bm{q}}}_{j}=\bar{\hat{\bm{f}}}_{k}.$$ (5.4)

Both equations (5.3) and (5.4) yield the same output for a given $\bar{\hat{\bm{f}}}_{k}=\hat{\bm{f}}_{-k}$ since $\bar{\hat{\bm{A}}}_{k-j}=\hat{\bm{A}}_{-k+j}$, inducing the harmonic resolvent response and forcing modes containing $\bm{U}_{H,j}=\bar{\bm{U}}_{H,-j}$ and $\bm{V}_{H,j}=\bar{\bm{V}}_{H,-j}$ for $j\neq 0$, respectively. Hence, we only keep $\lfloor N_{\omega}/2\rfloor+1$ Fourier coefficients of the forcing and response modes, reducing memory consumption by half when storing these dense matrices.

5.2 CPU benefits of RSVD-$\Delta t$

The computational cost of computing harmonic resolvent modes using RSVD-LU is predominantly determined by the LU decomposition of $\bm{T}$. In the case of resolvent analysis, the CPU scaling of the LU decomposition is $O(N^{1.5})$ and $O(N^{2})$ for two- and three-dimensional systems, respectively (Towne et al., 2022; Farghadan et al., 2023). When, in general, block matrix $\bm{T}$ contains $\hat{\bm{A}}_{p,i}$ terms, and yields a denser matrix than block diagonal with a complex sparsity pattern, the scaling is expected to be $O((N_{\omega}N)^{1.5})$ or $O((N_{\omega}N)^{2})$ or worse.

The cost of using RSVD-$\Delta t$ is directly tied to the cost of computing the actions of $\bm{H}$ and $\bm{H}^{*}$. In resolvent analysis, the total CPU cost is proportional to the state dimension $N$. This cost can be broken down into three main components: $(i)$ time-stepping, e.g., Adams-Bashforth methods, which scales as $O(N)$ when $\bm{A}$ is sparse, $(ii)$ creating forcing from Fourier space to the time domain, which scales as $O(NN_{t}N_{\omega})$, $(iii)$ taking the response from the time domain to Fourier space, which scales as $O(NN_{\omega}^{2})$. Among these components, the time-stepping part is the most significant contributor to the overall cost. The CPU cost of harmonic resolvent analysis can be assessed similarly to resolvent analysis, with one exception. For periodic flows, the LNS operator varies over time, whereas it remains constant for statistically stationary flows. We showed in $\lx@sectionsign$5.1 that the cost of creating LNS operators over time scales as $O(NN_{b})$ for sparse matrices.

In summary, the CPU and memory costs of RSVD-$\Delta t$ exhibit linear scaling with respect to the state dimension $N$, regardless of whether the base flow is steady or periodic. Moreover, the dominant terms in the CPU cost are independent of the number of retained frequencies $N_{\omega}$. Both of these properties make the RSVD-$\Delta t$ algorithm considerably more scalable than RSVD-LU for harmonic resolvent analysis.

6 Minimizing the CPU cost of RSVD-$\Delta t$

Similar to time stepping in the context of resolvent analysis, periodic systems may also experience a slow decay of transient response, which is undesirable, since we require the steady-state response in isolation to compute the action of the harmonic resolvent operator using time stepping. This phenomenon might be linked to instances when $\bm{T}$ is near singularity. Regardless, in cases where the LNS equations are absolutely stable but exhibit slowly decaying modes, we propose a strategy that can effectively reduce the CPU cost of time stepping.

6.1 Stability analysis: Floquet theorem and transient response

The steady-state response of the LNS equations, subject to forcing as described in equation (3.1), is of interest for our analysis. In the case of a linear time-invariant system where $\bm{A}$ is independent of time, the decay rate is controlled by the least-damped eigenvalue of $\bm{A}$. This decay rate can be estimated by integrating a homogeneous system over a sufficiently long period. However, for a linear time-periodic system, the decay rate is determined by the least-damped Floquet exponent.

Assuming $\bm{\Phi}(t)$ is the fundamental solution of (3.1), the monodromy matrix is constructed as

$$\bm{M}=\bm{\Phi}(0)^{-1}\bm{\Phi}(T),$$ (6.1)

where $T$ is the period of the flow. When $\bm{\Phi}(0)=\bm{I},\bm{\Phi}(t)$ is referred to as the principal fundamental matrix, and the monodromy matrix simplifies to $\bm{\Phi}(T)$. In other words, we evaluate the principal fundamental matrix at the end of the first period. The eigenvalues of $\bm{M}$, known as Floquet multipliers $\mu_{j}$, allow us to determine the Floquet exponents $\lambda_{j}=\log(\mu_{j})/T$. The least-damped Floquet exponent has the smallest real part and it determines the decay rate of the transient response. If any mode is located in the unstable half-plane, i.e., $\mathcal{R}e(\lambda_{j})>0$, the system is globally unstable.

It is worth noting that the transient length remains independent of the steady-state period. This characteristic is particularly important as it ensures that the length of integration does not rely on $T$, thereby avoiding potential deterioration in the performance of RSVD-$\Delta t$ for flows with low periodicity. In practice, both finding the principal fundamental matrix and determining the Floquet exponents can be computationally expensive, especially for large-scale systems. As with autonomous cases, it is possible to run a homogeneous simulation to analyze the long-term behavior of the transient response for periodic flows. This equivalence has been illustrated on a test case in $\lx@sectionsign$7.

If $\bm{T}$ is singular but all other modes are stable, then the Floquet exponent with smallest real part will be zero with $\widehat{d\bm{\bar{q}}/dt}$ as the corresponding Floquet mode. This can be seen from (3.12), where the phase shift direction is a non-trivial solution of the homogenous system. Leclercq & Sipp (2023) have also shown that $d\bm{\bar{q}}/dt$ is associated with a zero Floquet exponent under the assumption that the base flow satisfies the Navier-Stokes equations as in (3.11).

6.2 Efficient transient removal

Our strategy leverages the periodic nature of the steady-state part and the linear evolution of the transient part of the solution to directly compute and eliminate the undesired transient portion. For a given pair of solutions $\bm{q}_{1}$ and $\bm{q}_{2}$ at times $t_{1}$ and $t_{2}=t_{1}+T$, respectively, they can be expressed as a sum of their steady-state and transient components as

$$\bm{q}_{1} =\bm{q}_{s,1}+\bm{q}_{t,1},$$ (6.2)

$$\bm{q}_{2} =\bm{q}_{s,2}+\bm{q}_{t,2},$$

where $\bm{q}_{t,1}$ and $\bm{q}_{t,2}$ represent the transient part which decays to zero as $t\to\infty$ and $\bm{q}_{s,1}$ and $\bm{q}_{s,2}$ denote the steady-state part which evolves periodically. The time distance between two snapshots is one period, hence

$$\bm{q}_{s,2}=\bm{q}_{s,1}.$$ (6.3)

Also, the evolution of the transient part can be expressed as

$$\bm{q}_{t,2}=\bm{\Phi}(T)\bm{q}_{t,1},$$ (6.4)

where $\bm{\Phi}$ is the principal fundamental matrix of (3.1). Simplifying (6.2), (6.3), and (6.4) leads to

$$(\bm{\Phi}(T)-\bm{I})\bm{q}_{t,1}=\bm{b},$$ (6.5)

where $\bm{b}=\bm{q}_{2}-\bm{q}_{1}$ is known.

By obtaining $\bm{q}_{t,1}$ from solving (6.5), the steady-state solution can be recovered as $\bm{q}_{s,1}=\bm{q}_{1}-\bm{q}_{t,1}$. The central challenge in removing the transient is the computational cost associated with solving the linear system in (6.5). Instead, we employ Petrov-Galerkin (or Galerkin) projection to obtain an approximate solution to (6.5) in a more cost-effective manner.

A low-dimensional representation of the transient part can be expressed as

$$\bm{q}_{t,1}=\bm{\phi}\bm{\beta}_{1},$$ (6.6)

where $\bm{\phi}\in\mathbb{C}^{N\times r}$ is an orthonormal low-dimensional test basis with $r\ll N$, and $\bm{\beta}_{1}\in\mathbb{C}^{r}$ represents the coefficients describing the transient response in this basis. Substituting (6.6) into (6.5), the linear system

$$(\bm{\Phi}(T)-\bm{I})\bm{\phi}\bm{\beta}_{1}=\bm{b}$$ (6.7)

is overdetermined. We use Petrov-Galerkin projection with a low-dimensional trial basis $\bm{\psi}\in\mathbb{C}^{N\times r}$ to close (6.7), yielding

$$\tilde{\bm{M}}\bm{\beta}_{1}=\bm{\psi}^{*}\bm{b},$$ (6.8)

where

$$\tilde{\bm{M}}=\bm{\psi}^{*}(\bm{\Phi}(T)-\bm{I})\bm{\phi}$$ (6.9)

is the map between coefficients. Solving (6.8) for $\bm{\beta}_{1}$ is inexpensive due to its reduced dimension, and from (6.6) the transient estimate is obtained as

$$\bm{q}_{t,1}=\bm{\phi}(\bm{\psi}^{*}\bm{\phi}-\tilde{\bm{M}})^{-1}\bm{\psi}^{*}\bm{b}.$$ (6.10)

The reduced operator $\tilde{\bm{M}}$ is obtained by integrating the columns of $\bm{\phi}$ for one period, giving $\bm{\Phi}(T)\bm{\phi}$, and then projecting $(\bm{\Phi}(T)-\bm{I})\bm{\phi}$ against the columns of $\bm{\psi}$. This integration happens once and its CPU cost is equivalent to integrating the LNS equations for $r$ periods. The analogous process occurs for the adjoint equations.

How should the test and trial bases be selected? We have found a test basis spanning the least-damped modes to be effective. Rather than determining these modes through an exhaustive Floquet analysis, we employ a homogenous simulation as illustrated in the top row of the flowchart in figure 3. By initiating a sufficiently long simulation from a normalized random initial condition, the system asymptotically converges to the least-damped modes. The longer the integration length, the lower-dimensional the space becomes. As demonstrated in a test case in $\lx@sectionsign$7, the dimension of the test bases can, in some cases, be as small as a single column. Our test cases employ Galerkin projection with identical test and trial bases.

Figure 3: Flowchart depicting the action of $\bm{H}$ using time stepping (a) with efficient transient removal strategy and (b) with natural transient decay.

The obtained test basis effectively captures the transient response at a specific phase $\theta$. Utilizing $\bm{\phi}$ and performing Galerkin projection, we obtain $\tilde{\bm{M}}$. The time-stepping is then carried out for a sufficiently long duration to allow the initial transients to vanish before applying (6.8) to determine the transient at phase $\theta$ as shown in the middle row of the flowchart in figure 3. Once the steady-state response at the same phase is updated as $\bm{q}_{s,1}=\bm{q}_{1}-\bm{q}_{t,1}$, it is used as a new initial condition, which is synchronized with the forcing and will not excite a transient response (within the span of the bases used to construct $\tilde{\bm{M}}$). Then, integration for one period is sufficient, and $N_{\omega}$ steady-state snapshots are collected as desired. The resulting steady-state snapshots are identical to those obtained by a prolonged wait in the case of natural decay, as illustrated in the bottom row of the flowchart in figure 3. This procedure is also applicable to the adjoint LNS equations. Given that the adjoint equations differ from the LNS equations, a new basis needs to be constructed for them. Implementing this strategy can reduce the integration length by a factor of 10 or more, depending on the desired accuracy.

7 Test cases

We assess the effectiveness of RSVD-$\Delta t$ in periodic systems using two test cases. First, we use a modified Ginzburg-Landau equation to validate the RSVD-$\Delta t$ algorithm and to illustrate the transient removal strategy. Second, we consider a periodic flow around an airfoil, similar to a test case from Padovan et al. (2020), to evaluate the accuracy and cost-savings of the RSVD-$\Delta t$ algorithm compared to RSVD-LU.

7.1 Periodically varying complex Ginzburg-Landau equation

The one-dimensional complex Ginzburg-Landau equations is a widely used model for understanding and controlling the non-modal growth of instabilities in transitional and turbulent shear flows (Chomaz et al., 1988; Hunt & Crighton, 1991; Bagheri et al., 2009; Chen & Rowley, 2011; Cavalieri et al., 2019). Here, a modified system is used as an inexpensive test case to validate our algorithm. The original complex Ginzburg-Landau system follows the form of (2.1) with

$$\begin{gathered}\mathcal{A}=-\nu\frac{\partial}{\partial x}+\gamma\frac{\partial^{2}}{\partial x^{2}}+\mu(x),\\ \mu(x)=(\mu_{0}-c_{\mu}^{2})+\frac{\mu_{2}}{2}x^{2},\\ \bm{B}=\bm{C}=\bm{I}.\end{gathered}$$ (7.1)

Following Bagheri et al. (2009), we set $c_{\mu}=0.2,\mu_{2}=-0.01,\gamma=1-\text{i},$ and $\nu=2+0.2\text{i}$. This system is globally stable when $\mu_{0}<\mu_{0,cr}\approx 0.3977$ (Bagheri et al., 2009). By substituting $\mu_{0}(t)=\bar{\mu}_{0}+\mu_{p}sin(\omega_{f}t-\pi/2)$ in place of $\mu_{0}$, we transform the Ginzburg-Landau equations to a periodically varying system with fundamental frequency $\omega_{f}$, following the form of (3.1). We set $\mu_{0}=0.395$, $\mu_{p}=0.1$, and $\omega_{f}=0.1$; the mean value of $\mu(t)$ is equal to $\mu_{0}$, whose value is close to the critical value $\mu_{0,cr}$ to promote significant growth. The periodic linear operators $\bm{A}_{p,j}$ are constructed using a central finite-difference method for $x\in[-50,50]$ with $N=1000$ grid points, and we set $\bm{W}=\bm{I}$ on account of the uniform grid.

Nine operators are built within the interval $[0,T)$, with $T=2\pi/\omega_{0}\approx 62.83$, such that $\bm{A}_{p}(t)=\bm{A}_{p}(t+T)$, before conducting a DFT to obtain $\hat{\bm{A}}_{p}$. Figure 4(a) depicts the Frobenius norm of $\hat{\bm{A}}_{p}$ up to the fourth harmonic, normalized relative to the maximum norm. The spectrum reveals that only the mean and the first harmonic are active, as expected since $\mu_{0}(t)$ oscillate at a single frequency. Accordingly, $\bm{A}_{p}$ can be reconstructed at any given time using $\hat{\bm{A}}_{p,0}$, and $\hat{\bm{A}}_{p,\pm 1}$, and the harmonic resolvent operator is constructed around $\Omega_{\bar{q}}=\{-1,0,1\}\omega_{f},$ i.e., $N_{b}=3$. Notably, the first harmonic’s norm is several orders of magnitude smaller than that of the base flow. Nonetheless, it should not be disregarded, as we will elucidate its significance in subsequent discussions.

Figure 4: Ginzburg-Landau test case: (a) the normalized Frobenius norm of $\hat{\bm{A}}_{p,\omega}$ up to the fourth harmonic. (b) The norm of the transient response over time (blue) along with the expected decay rate from Floquet analysis (red).

An essential initial step involves ensuring the stability of $\bm{A}_{p}$. Typically, this is achieved by integrating the homogeneous system, i.e., integrating (3.1) with zero forcing, starting from a normalized random initial condition. However, the compact size of this system allows us to validate the time-stepping approach for computing the least-damped decaying mode against a Floquet stability analysis. The least-damped Floquet exponent, computed as $\lambda_{l}=-0.0021-0.0477\text{i}$, is expected to govern the decay rate of the transient response. Figure 4(b) displays the norm of snapshots of the homogeneous system response over time, along with a reference line representing $e^{\lambda_{l}t}$. The natural decay rate aligns with the reference line, as expected, ensuring that the transient run will ultimately converge to the least-damped Floquet mode. Additionally, the system is not singular, as the transient dynamics never reach a naturally stable state, and none of the Floquet exponents have a zero real value.

Before proceeding to the computation of harmonic resolvent modes and gains using RSVD-$\Delta t$ and RSVD-LU, the transient removal strategy is demonstrated, and computing the actions of $\bm{H}$ and $\bm{H}^{*}$ are validated for a given random forcing containing $\Omega_{q}\in\{-10,-9,...,9,10\}\omega_{f}$, i.e., $\omega\in[-1,1]$ with $\Delta\omega=0.1$. The responses have been computed up to the $20^{th}$ harmonic, i.e., $\omega\in[-2,2]$, using both time stepping and directly in Fourier space as outlined in (3.7). We used a $4^{th}$ order Runge–Kutta (RK4) integration scheme and a time step of $dt=0.001$ and $T_{t}=15000$ for this purpose.

Figure 5(a) depicts the spectrum of response norms. In order to highlight the substantial difference in output generated by the harmonic resolvent operator around the base flow ($\Omega_{\bar{q}}=\{0\}\omega_{f},N_{b}=1$), the norm of the responses for the same input forcing is also presented. Figure 5(b) displays the relative errors. The harmonic resolvent with $N_{b}=1$ produces inaccurate results across the entire spectrum and has zero response for $|\omega|>1$ since the modes at different frequencies are entirely decoupled. Therefore, neglecting $\hat{\bm{A}}_{p,\pm 1}$ is impractical, as mentioned earlier, and induces significant alterations in the harmonic resolvent modes and gains. On the other hand, the time-stepping output exhibits almost machine-precision accuracy within $|\omega|\leq 1$, and its accuracy diminishes as response norms decrease within $|\omega|>1$. The validation of $\bm{H}^{*}$ follows a similar procedure but is omitted here for brevity.

Figure 5: Ginzburg-Landau test case: (a) the norms of the Fourier space response for $N_{b}=3$ (black) and $N_{b}=1$ (blue), along with the norm of the steady-state response obtained via time-stepping (red). (b) The relative errors of the norms, with the norms from Fourier space with $N_{b}=3$ (black line in (a)) serving as the reference.

The transient norm decay rate, as shown in figure 4(b), follows $e^{-0.0021t}$. This indicates that approximately 40 periods are necessary for the transient norm to decay to 1%. However, we can significantly reduce this duration by employing our transient removal strategy. To begin, we precompute the test basis by conducting a sufficiently long homogeneous simulation, and we find that a single mode suffices. Next, we integrate the orthogonalized basis for one period to obtain $\tilde{M}$ and complete the preprocessing step.

Employing the same forcing as the validation case, we obtain the responses in time at the identical phase (initiated at $t=0$ or $\theta=0$ with a time interval of $\Delta t=T$), and compare them with the solutions after the removal of transients. Figure 6(a) presents the norm of the responses, which illustrates the effectiveness of our strategy commencing from the second period. This plot suggests that after two periods using our transient removal strategy with Galerkin projection, we obtain an accurate steady-state snapshot which can be used as the synchronized initial condition for the forced equations. Figure 6(b) demonstrates that the snapshots ultimately align perfectly with the updated snapshot. Hence, using our strategy, a total of three periods are required – two periods before transient removal and one period after the initial condition is obtained – to acquire all $N_{\omega}$ steady-state snapshots, a 10-fold speed-up compared to the natural decay. A similar observation was made regarding the adjoint equations. In brief, our strategy can significantly reduce computational costs, potentially by one or more orders of magnitude, depending on the desired level of accuracy.

Figure 6: Ginzburg-Landau test case: (a) the response norms before (black) and after (red) transient removal at the same phase, i.e., snapshots with $\Delta t=T$ interval. (b) Comparison of the normalized response at the second period in blue, the $50^{th}$ period in black, and the second period after transient removal in red.

Finally, RSVD-$\Delta t$ employing an RK4 integration scheme with $dt=0.003$ and $T_{t}=2T$, is carried out along with our efficient transient removal strategy to compute the harmonic resolvent modes for $k=5$ and $q=1$. To establish a reference, the harmonic resolvent modes are also computed using RSVD-LU with an equivalent number of test vectors and power iterations. Figure 7(a) illustrates that RSVD-$\Delta t$ computes the five leading gains identically to RSVD-LU. The gain relative error in figure 7(b) confirms that relative errors remain below $10^{-9}$ across all modes. The inner products of response and forcing modes, computed via RSVD-$\Delta t$ and RSVD-LU, demonstrate parallel directions as the relative errors remain below $10^{-8}$ across all modes. The inner-product error between two unit-norm vectors $\bm{v}_{1}$ and $\bm{v}_{2}$ is defined as

$$e_{ip}=1-\langle\bm{\bm{v}_{1}},\bm{v_{2}}\rangle.$$ (7.2)

Figure 7: Ginzburg-Landau test case: (a) the five leading gains using RSVD-$\Delta t$ (blue) and RSVD-LU (red). (b) The corresponding gain relative error (purple) and relative inner product errors (7.2) between the response modes (cyan) and the forcing modes (green) computed via RSVD-$\Delta t$ and RSVD-LU.

7.2 Flow over an airfoil

Our second test case consists of the flow over a NACA0012 airfoil at a low Reynolds number of $Re=200$ and a steep angle of attack of $\alpha=20^{\circ}$, which is dominated by periodic motions. This problem serves as a benchmark for computational cost and accuracy comparisons. A direct numerical simulation is conducted using the “CharLES” compressible flow solver. To mimic the incompressible simulations by Padovan et al. (2020), we set the Mach number to 0.05. The simulation employs a C-shaped mesh with a total of 62,000 cells. The leading edge of the airfoil is located at the origin $(x/L_{c},y/L_{c})=(0,0)$, where $L_{c}=1$ is the chord length. The computational domain spans $x/L_{c}\times y/L_{c}\in[-49,50]\times[-50,50]$. A constant time step of $\Delta tU_{\infty}/L_{c}=6.88\times 10^{-5}$ is used, corresponding to a Courant-Friedrichs-Lewy number of 0.91, where $U_{\infty}$ is the inflow streamwise velocity. Integration continues for a sufficient duration until the flow reaches a periodic state. Figure 8(a) displays the power spectral density computed from the transverse velocity at $(x/L_{c},y/L_{c})=(3.0,-0.43)$, where a pronounced vortex shedding occurs behind the trailing edge. The shedding frequency is $St_{f}=0.114$, where the Strouhal number is defined as $St=\frac{\omega L_{c}sin(\alpha)}{2\pi U_{\infty}}$. The consistency between the dominant frequency identified in the power spectral density and obtained from our linearized code has been confirmed in our previous work (Jung et al., 2023).

Figure 8: Airfoil test case: (a) the PSD spectrum based on transverse velocity at $(x/L_{c},y/L_{c})=(3.0,-0.43)$. (b) The normalized Frobenius norm of $\hat{\bm{A}}_{p,St}$ up to the eighth harmonic.

The linearized operators are constructed at 100 time points within $T=2\pi/St_{f}\approx 55$. Performing a DFT on the $\bm{A}_{i}$, we generate 100 LNS operators $\hat{\bm{A}}_{j}$ in Fourier space to be used for harmonic resolvent analysis. The Frobenius norms $||\hat{\bm{A}}_{p,St}||_{F}$, depicted in figure 8, indicate that $\bm{A}_{p}(t)$ can be accurately reconstructed using up to the 5${}^{th}$ harmonic. Given that the linearized operator $\bm{A}_{p}(t)$ is real-valued for this problem, we retained only $\lfloor N_{b}/2\rfloor+1=6$ coefficients, i.e., the zeroth and 5 positive harmonics. For time stepping, $dt=0.0045$ is chosen to ensure the stability and accuracy of the RK4 integration scheme. The input and output perturbation frequency range spans up to the $7^{th}$ harmonic, i.e., $N_{\omega}=15$. The domain of interest is $x/L_{c}\times y/L_{c}\in[-4,12]\times[-2.5,2.5]$, identical to that in Padovan et al. (2020). We seek optimal forcing and response under the Chu energy norm (Chu, 1965), which has been used in several previous studies (Towne et al., 2018; Schmidt et al., 2018; Heidt & Colonius, 2023).

A homogeneous simulation of the time-periodic system exhibits a slow decay rate, so we employ the transient removal strategy. Three snapshots constitute the test basis, obtained after a sufficiently long time interval where most initial transients are eliminated. A random forcing, encompassing the set $\Omega_{q}\in\{-7,-6,...,6,7\}St_{f}$ frequencies, is applied to the LNS equations for over 30 periods. The snapshot norms slowly approach towards the steady-state norm but remain far from converging, as illustrated in figure 9. Utilizing our strategy, the relative norm error falls below the 1% threshold after 20-30 periods in most cases with different random forcings. The efficacy of our strategy is further emphasized in figure 9(b), where the 1% relative error is achieved before 10 periods for the optimal forcing (from RSVD-$\Delta t$ output). Ultimately, we set $T_{t}=29T$ and conducted our simulations for a total duration of 30 periods to compute the actions of $\bm{H}$ and $\bm{H}^{*}$. Extrapolating the natural decay rate suggests that, in the absence of our transient-removal strategy, the norms reach the 1% threshold after around 2000 periods, demonstrating over a 60-fold acceleration in time stepping and overall algorithm performance with our strategy. This observation holds true for the adjoint equations as well.

Figure 9: Airfoil test case: (a) the response norms before (black) and after (red) transient removal, and the true steady-state norm (blue), at the same phase, i.e., snapshots with $\Delta t=T$ interval, for a random forcing. (b) The relative error between steady-state response norm and the response norms before (black) and after (red) transient removal.

The number of test vectors is set to $k=5$, and $q=2$ power iteration proves sufficient for the convergence of both gains and modes. The suboptimal forcing and response modes are shown in figure 10. The results from RSVD-$\Delta t$ closely align with the existing data obtained from RSVD-LU by Padovan et al. (2020), despite variations in the CFD solver, boundary conditions, domain setup, and energy norm. The gain plot in figure 11(a) exhibits a similar pattern to Padovan et al. (2020), with the exception that we have retained the optimal gain associated with the phase shift (the exact values differ due to differences in the problem setup). We noted that the leading mode and gain, i.e., the mode associated with the phase shift, converged without power iterations due to a substantial two-orders-of-magnitude gap between the optimal and suboptimal gains.

Figure 10: Airfoil test case: real part of the vorticity field computed from (a, c, e, g) the input mode and (b, d, f, h) the output mode associated with the first suboptimal gain.

From a computational standpoint, a fair comparison between the RSVD-LU and RSVD-$\Delta t$ algorithms is feasible when both utilize the same parameters. However, the use of $N_{b}=11$ proves to be excessively memory-intensive for the RSVD-LU algorithm, exceeding the available 3.5 TB of memory on our cluster. The CPU and memory scaling plots in figure 11(b) and (c), respectively, depict the LU decomposition of $\bm{T}$ with $N_{b}=3$, maintaining a constant $N$ while varying the number of blocks $N_{\omega}$ from 5 to 15 with $\Delta N_{\omega}=2$. On the same plots, we present the cost of computing the action of $\bm{H}$ on a vector using time-stepping. The total duration is set to 30 periods to obtain an accurate solution. We only vary $N_{\omega}$ while keeping $N_{b}=3$, the overall cost of RSVD-$\Delta t$ remains constant and unaffected by changes in $N_{\omega}$. This implies that the creation of LNS operators and the time-stepping process incur significantly higher costs compared to the transformations between Fourier space and the time domain for forcing and response. The memory of RSVD-$\Delta t$ to store LNS operators remains independent of $N_{\omega}$, as expected (see $\lx@sectionsign$5.1), while storing $\bm{Q}$ and $\bm{F}$ matrices grows linearly with number of frequencies $N_{\omega}$.

RSVD-LU is estimated to require 1543 CPU-hours for the LU decomposition of $\bm{T}$, and an additional 15 CPU-hours for solving the LU-decomposed system for each test vector, yielding a total cost of $\text{CPU}_{\text{{RSVD-LU}}}\approx 1543+15\times 5\times 2\times 2=1843$ CPU-hours for $k=5,q=1$. On the other hand, employing RSVD-$\Delta t$ with the aforementioned time-stepping parameters and $N_{b}=11$ and $N_{\omega}=15$ incurs a cost of approximately 7 CPU-hours per period for a test vector, resulting in a total cost of $\text{CPU}_{\text{{RSVD-$\Delta t$}}}\approx 30\times 7\times 5\times 2\times 2=4200$ CPU-hours. The cost of QR decomposition and the final SVD are ignored as they are orders of magnitude faster. In terms of memory consumption, RSVD-LU peaks at $\text{RAM}_{\text{RSVD-LU}}\approx 2036GB$ for the LU decomposition of $\bm{T}$. RSVD-$\Delta t$ , on the other hand, stores $\lfloor N_{b}/2\rfloor+1=6$ sparse matrices, each of size 0.4 GB, and three dense matrices of size $Nk(\lfloor N_{\omega}/2\rfloor+1)=313630\times 5\times 8$, totaling $\text{RAM}_{\text{{RSVD-$\Delta t$}}}\approx 0.4\times 6+0.185\times 3\approx 3GB$. This translates to a memory saving of approximately three orders of magnitude, even with an unbalanced number of base frequencies. Overcoming the memory consumption hurdle, which is typically the limiting factor in practice, RSVD-$\Delta t$ emerges as a viable tool for analyzing larger-scale flows compared to RSVD-LU.

Figure 11: Airfoil test case: (a) the five leading gains using RSVD-$\Delta t$ . (b) The CPU-hour, and (c) memory usage scaling of RSVD-LU (red) and RSVD-$\Delta t$ (blue) to compute the action of $\bm{H}$ onto a vector, i.e., $k=1$. The memory usage of RSVD-$\Delta t$ is decomposed into memory required to store LNS operators (solid) and forcing and response matrices (dashed).

8 Conclusions

This paper presents an extension of the RSVD-$\Delta t$ algorithm that was originally developed for resolvent analysis of steady base flows to handle harmonic resolvent analysis of periodic flows. Specifically, we demonstrate that the time stepping technique employed within RSVD-$\Delta t$ can effectively replace the actions of $\bm{H}$ and its complex conjugate $\bm{H}^{*}$ in Fourier space. In terms of CPU cost, the RSVD-$\Delta t$ algorithm exhibits a scaling of $O(N)$ for both resolvent and harmonic resolvent analyses. This offers a significant advantage, considering that the LU decomposition of $\bm{T}$ scales as $O((N_{\omega}N)^{c})$ with $c\geq 1.5$. Regarding memory usage, our algorithm only stores relevant matrices in Fourier space and exhibits $O(NN_{\omega})$ scaling; in contrast, the memory peak of LU decomposition of $\bm{T}$ empirically scales with $O((N_{\omega}N)^{1.5})$ or worse. One difference between applying RSVD-$\Delta t$ to resolvent and harmonic resolvent analyses lies in the necessity to update the LNS operators during the time-stepping process for the latter. These operators are efficiently generated on the fly from their $N_{b}$ Fourier components.

The error sources in our algorithm align with those in resolvent analysis, extensively investigated in Farghadan et al. (2023). A unique contribution of this paper is the introduction of a novel transient removal strategy for harmonic resolvent analysis. While reminiscent of our approach for resolvent operators, the challenge in harmonic resolvent analysis lies in the intertwining of all retained frequencies. Our strategy takes advantage of the differing evolution of the steady-state and transient components of the response and Petrov-Galrkin or Galerkin projections.

A validation of RSVD-$\Delta t$ against RSVD-LU is conducted using a periodic Ginzburg-Landau system. Additionally, we verify the governing role of the Floquet exponent in determining the decay rate of the least-damped mode of the system. Extending the application of our algorithm, we analyze a two-dimensional flow passing an airfoil, providing insights into the CPU and memory complexities. Despite dealing with a mid-sized flow scenario, a substantial memory gap persists between RSVD-LU and RSVD-$\Delta t$ algorithms. The computed forcing and response modes corresponding to the first suboptimal gain closely resemble those obtained by Padovan et al. (2020). The adaptation of the transient removal strategy for these test cases significantly enhanced the performance of RSVD-$\Delta t$, resulting in 10-fold speed-up for the Ginzburg-Landau problem and 60-fold speed-up for the airfoil to reach a 1% relative error. The speed-up is even greater for lower error tolerances.

Moreover, our algorithm exhibits versatility, accommodating non-identity weight, input, and output matrices. This capability extends to computing modes for the modified harmonic resolvent operator, and in particular, the analysis of cross-frequency amplification mechanisms via $\bm{H}_{\omega_{2},\omega_{1}}$. Our algorithm tackles the computation of subharmonic resolvent modes if desired. Implementation of the RSVD-$\Delta t$ algorithm within Petsc (Balay et al., 2019) and Slepc (Hernandez et al., 2005) environments leverage parallel computations for enhanced efficiency. Our time-stepping approach allows for a matrix-free application. This can be implemented using any code equipped with linear direct and adjoint capabilities, bypassing the explicit formation of $\bm{T}$ (de Pando et al., 2012; Martini et al., 2021). Finally, the time-stepping within our algorithm can also be used along with other SVD algorithms, e.g., Arnoldi and power iteration, if desired.

In turbulent flows, where the state dimension $N$ can become very large due to higher resolution requirements, CPU cost and memory requirements can swiftly impose constraints on the applicability of RSVD-LU. This limitation applies even to steady problems and is further exacerbated for harmonic resolvent analysis due to the inflated frequency-domain operators necessitated by triadic frequency coupling. Our algorithm circumvents this issue by strategically operating in the time domain, avoiding the need for limiting operations like LU decomposition, leading to linear cost scaling. Because of this key distinction, the RSVD-$\Delta t$ could enable harmonic resolvent analysis of previous intractable turbulent flows.

Appendix A RSVD-$\Delta t$ for the subharmonic resolvent operator

In $\lx@sectionsign$3.5, we provided an overview of subharmonic resolvent analysis. In this appendix, we briefly outline the application of RSVD-$\Delta t$ in computing subharmonic resolvent modes and gains. Consider the frequency of interest as $\gamma\in\Omega_{\gamma}$. This set specifies the perturbation frequency, while the base flow frequency content is confined to $\Omega_{\bar{q}}$.

To compute the actions of $\bm{H}$ and $\bm{H}^{*}$ using time-stepping, we must adhere to both the base flow frequency, enforcing a duration of $T=2\pi/\omega_{f}$, and the perturbation frequency, enforcing a duration of $T_{p}=2\pi/\gamma$, in order to obtain $N_{\omega}$ steady-state solutions to (3.1). Thus, we need to integrate for $T_{sub}$ such that $T_{sub}/T\in\mathbb{N}$ and $T_{sub}/T_{p}\in\mathbb{N}$, during which $N_{\omega}$ equidistant snapshots are saved. For instance, if we consider $\gamma=\omega_{f}/5$, requiring $T_{p}=2\pi/(\omega_{f}/5)=5T$, integrating for $T_{sub}=5T$ is sufficient to obtain the steady-state solutions. All the other steps remain the same as introduced in Algorithm 1.