Study of dynamics in post-transient flows using Koopman mode decomposition

TL;DR

This paper applies Koopman Mode Decomposition to analyze complex post-transient flow dynamics in a lid-driven cavity, revealing spectral properties and flow regime transitions with improved flow reconstruction capabilities.

Contribution

It adapts spectral analysis algorithms for Koopman operator spectra and links spectral properties to flow regime changes, advancing flow analysis methods.

Findings

Koopman spectrum characterizes flow regimes from steady to aperiodic.

Koopman modes are robust across different flow regimes.

KMD outperforms POD in reconstructing flows with quasi-periodic components.

Abstract

The Koopman Mode Decomposition (KMD) is a data-analysis technique which is often used to extract the spatio-temporal patterns of complex flows. In this paper, we use KMD to study the dynamics of the lid-driven flow in a two-dimensional square cavity based on theorems related to the spectral theory of the Koopman operator. We adapt two algorithms, from the classical Fourier and power spectral analysis, to compute the discrete and continuous spectrum of the Koopman operator for the post-transient flows. Properties of the Koopman operator spectrum are linked to the sequence of flow regimes occurring between $Re=10000$ and $Re=30000$, and changing the flow nature from steady to aperiodic. The Koopman eigenfunctions for different flow regimes, including flows with mixed spectra, are constructed using the assumption of ergodicity in the state space. The associated Koopman modes show remarkable robustness even as the temporal nature of the flow is changing substantially. We observe that KMD outperforms the Proper Orthogonal Decomposition in reconstruction of the flows with strong quasi-periodic components.c features are present in the flow.

Tables

Table 1. Table 1: Koopman frequencies with highest kinetic energy in the lid-driven cavity flow.

$Re$=13000		$Re$=16000		$Re$=19000
(periodic)		(quasi-periodic)		(mixed spectra)
$St$	$k$	$St$	$𝐤$	$St$	$𝐤$
0.000000	0	0.000000	(0,0)	0.000000	(0,0)
0.159826	1	0.155375	(1,0)	0.243822	(1,1)
0.319653	2	0.252287	(1,1)	0.150046	(1,0)
0.479479	3	0.096911	(0,1)	0.206310	(2,-1)
0.639303	4	0.193823	(0,2)	0.093778	(0,1)
0.799130	5	0.310751	(2,0)	0.300089	(2,0)
0.958954	6	0.213839	(2,-1)	0.187556	(0,2)
		0.349197	(1,2)	0.056266	(1,-1)
		0.407661	(2,1)	0.450133	(8,-8)
		0.058463	(1,-1)	0.393868	(7,-7)
		0.116927	(2,-2)	0.337602	(6,-6)
		0.135359	(-1,3)	0.225068	(4,-4)
		0.446109	(1,3)	0.487643	(2,2)
		0.290733	(0,3)	0.018756	(2,-3)
		0.038447	(-1,2)	0.431379	(6,-5)
		0.504573	(2,2)	0.543910	(8,-7)
		0.020017	(2,-3)	0.375111	(0,4)
		0.232270	(-1,4)	0.468886	(0,5)
		0.466125	(3,0)	0.693955	(9,-7)
		0.272303	(3,-2)	0.581423	(2,3)

Equations

$${\partial }_t\mathbf{u}(\mathbf{x},t)=\mathbf{F}(\mathbf{u}(\mathbf{x},t)).$$

$$g^{\tau }(\mathbf{u}(\mathbf{x},t)):=U^{\tau }g(\mathbf{u}(\mathbf{x},t))=g(\mathbf{u}(\mathbf{x},t+\tau )),$$

$${\varphi }^{\tau }(\mathbf{u}):=U^{\tau }{\varphi }_j(\mathbf{u})=e^{{\lambda }_j\tau }{\varphi }_j(\mathbf{u})$$

$${\lambda }_j=i{\omega }_j,{\omega }_j\in \mathbb{R}.$$

$$g(\mathbf{u})=\sum _{j=1}^{\infty }{g}_j{\varphi }_j(\mathbf{u}),$$

$$g^{\tau }(\mathbf{u})=U^{\tau }g(\mathbf{u})=\sum _{j=1}^{\infty }{g}_j{\varphi }_j(\mathbf{u})e^{i{\omega }_j\tau }.$$

$${\mathbf{g}}^{\tau }(\mathbf{u}):=U^{\tau }\mathbf{g}(\mathbf{u})=\sum _{j=1}^{\infty }{\mathbf{g}}_j{\varphi }_j(\mathbf{u})e^{i{\omega }_j\tau },$$

$$<f,g{>}_{\mathcal{H}}={\int }_{A}f{g}^{\ast }d\mu .$$

$$U^{\tau }\mathbf{g}\stackrel{\mu }{=}\sum _{\mathbf{k}\in {\mathbb{Z}}^m}{\mathbf{g}}_{\mathbf{k}}{\varphi }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{\Omega }\tau }.$$

$$U^{\tau }\mathbf{g}({\mathbf{u}}_0)=\sum _{\mathbf{k}\in {\mathbb{Z}}^m}{\mathbf{g}}_{\mathbf{k}}{\varphi }_{\mathbf{k}}({\mathbf{u}}_0)e^{i\mathbf{k}\cdot \mathbf{\Omega }\tau }.$$

$$U^{\tau }g\stackrel{\mu }{=}\sum _{\mathbf{k}\in {\mathbb{Z}}^m}{g}_{\mathbf{k}}{\varphi }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{\Omega }\tau }+{\int }_{-\infty }^{\infty }{e}^{i\alpha \tau }d{E}_{\alpha }(g).$$

$$<g,U^{\tau }g{>}_{\mathcal{H}}=\sum _{k=1}^{\infty }|g_k{|}^2{e}^{i{\omega }_k\tau }+{\int }_{-\infty }^{\infty }{e}^{i\alpha \tau }{\rho }_g(\alpha )d\alpha .$$

$${\mu }_g(I)={\int }_I{\rho }_g(\alpha )d\alpha .$$

$$<g,U^{\tau }g{>}_{\mathcal{H}}=r_g(\tau ):=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^{T}g(t)g^{\ast }(t+\tau )dt.$$

$$\{X_t{\}}_{t\in \mathbb{R}}.$$

$$\mathbb{E}(X_t)=\mathbb{E}(X_{t+\tau })=m,\text{for all }\tau \in \mathbb{R},$$

$$cov(X_t,X_{t+\tau })=\mathbb{E}((X_t-m)(X_{t+\tau }-m))=cov(\tau ).$$

$$\{U^{t}g{\}}_{t\in \mathbb{R}},$$

$$\mathbf{P}(B)=\mu\big{(}(U^{t}g)^{-1}(B)\big{)},\quad B\subset\mathbb{R}.$$

$$\mathbb{E}(U^{t}g)=<U^{t}g,1{>}_{\mathcal{H}}=<g,U^{-t}1{>}_{\mathcal{H}}=<g,1{>}_{\mathcal{H}}=\mathbb{E}(g),$$

$$cov(U^{t}g,U^{t+\tau }g)$$

$$cov(g,U^{\tau }g)$$

$${\mathbf{u}}_j={\nabla }^{\perp }{\psi }_j,j=1,2,3,\dots .$$

$${\mathbf{\zeta }}_j=\nabla \times {\mathbf{u}}_j,j=1,2,3,\dots .$$

$$u_{lid}=(1-x^2{)}^2,x\in [-1,1].$$

$$\frac{\partial }{\partial t}{\nabla }^2\psi +\frac{\partial \psi }{\partial y}\frac{\partial }{\partial x}{\nabla }^2\psi -\frac{\partial \psi }{\partial x}\frac{\partial }{\partial y}{\nabla }^2\psi =\frac{1}{Re}{\nabla }^4\psi ,$$

$$\psi\bigg{|}_{\partial\Omega}=0\quad and\quad\frac{\partial\psi}{\partial n}\bigg{|}_{\partial\Omega}=u_{w},$$

$$(x_i,y_j)=\left(\cos (\frac{i\pi }{N}),\cos (\frac{j\pi }{M})\right)$$

$$i=0,1,\dots ,N,j=0,1,\dots ,M.$$

$$\psi (x,y)=(1-x^2)(1-y^2)q(x,y).$$

Full text

††preprint: APS/123-QED

I Introduction

In 1931, Bernard Koopman offered a new formulation of Hamiltonian mechanics based on the theory of Hilbert spaces Koopman (1931). In his formulation, the central object was a linear transformation, called the Koopman operator, which described the time evolution of observations on a Hamiltonian system. The potential applications of Koopman’s work to general theory of dynamical systems went mostly unrecognized for a long time, until it was brought to attention, and further developed in the context of spectral analysis and dissipative systems Mezić and Banaszuk (2004); Mezić (2005). In particular, Ref. Mezić, 2005 presented a linear expansion to describe the temporal evolution of observables on a nonlinear dynamical system in terms of the eigenfunctions and eigenvalues of the Koopman operator. This expansion, known as the Koopman Mode Decomposition (KMD), introduced the concept of Koopman modes, which are projection of observable evolution onto eigenfunctions of the Koopman operator, and give spatial “shapes” that evolve exponentially in time - with exponents which could be complex numbers. Since then, KMD has been used as a tool of data-driven analysis with applications ranging from power network stability analysis to pattern detection in neural networks Susuki and Mezić (2014); Giannakis et al. (2015); Georgescu and Mezić (2015); Erichson et al. (2015); Brunton et al. (2016); Mann and Kutz (2016)

In high-dimensional systems like fluid flows, direct analysis of the state space is computationally prohibitive and one often relies on the data-driven methods to extract the underlying spatio-temporal features from data obtained by numerical simulation or experiment. For this purpose, the KMD was introduced to fluid mechanics by Rowley et al. (2009), with an application to the problem of jet in cross flow. Rowley et al. showed how KMD can extract various oscillation frequencies and their associated spatial structures from the data. They also discovered the connection between the KMD and the numerical algorithm called Dynamic Mode Decomposition (DMD), which had been proposed by P.J. Schmid Schmid and Sesterhenn (2008). Since then, KMD and its numerical counterpart DMD, have become popular tools for extraction of physically-relevant time scales and their associated spatial structures in complex flows (Schmid, 2010; Schmid et al., 2011; Pan et al., 2011; Seena and Sung, 2011; Muld et al., 2012; Hua et al., 2016). Furthermore, the dynamical-systems origin of KMD has motivated a number of its applications in study of dynamical phenomena such as stability, bifurcation and transition in flows. Some examples are study of bifurcation and transition in flow past a cylinder (Chen et al., 2012; Bagheri, 2013), detection of significant structures in boundary layer transition (Sayadi et al., 2014; Subbareddy et al., 2014), bifurcation analysis using parametric DMD (Sayadi et al., 2015) and identification of flow regimes in thermo-fluid systems using sparse sensing (Kramer et al., 2015).

In this work, we discuss KMD as a versatile tool to identify and analyze the state space dynamics of post-transient flows (i.e. state space trajectory evolving on the attractor). The traditional approach to determine the dynamic regime of fluid systems is to look at the Fourier or power spectrum of time series . The existence of sharp peaks in those spectra is deemed to indicate periodic motion while broadband spectrum is often interpreted as a sign of chaos (e.g. Swinney and Gollub, 1978; Farmer et al., 1980; Karniadakis and Triantafyllou, 1992; Tomboulides and Orszag, 2000; Peng et al., 2003; Basley et al., 2011). We point out the relationship between this classical viewpoint and the Koopman spectral analysis, and describe how the Koopman spectrum of data can be used to determine the geometry of the attractor, using the results in Mezić (2017). The Koopman viewpoint generalizes the classical spectral analysis through the notion of Koopman eigenfunctions and modes. In particular, the Koopman eigenfunctions provide linearly evolving coordinates in the state space, and enable the use of state-space analysis and control techniques which are not reachable by classical spectral analysis. For example, the knowledge of Koopman eigenfunctions can be used to construct nonparametric predictors Giannakis (2017); Korda and Mezić (2016), state estimators Surana and Banaszuk (2016); Surana (2016) and nonlinear controllers Korda and Mezić (2016); Kaiser et al. (2017) using linear system strategies. The utility of Koopman eigenfunctions in flow prediction and control is further discussed in Taira et al. (2017). In this paper, we present a new way to construct (and visualize) the Koopman eigenfunctions in post-transient flows which are ergodic in the state space.

We apply the Koopman analysis to the two dimensional lid-driven cavity flow with regularized lid velocity. This flow provides a good benchmark for our analysis, since it shows a wide range of dynamic behavior over various Reynolds numbers. The dynamics of each flow regime is discussed in terms of the Koopman spectral properties: the Koopman spectrum determines the type of the attractor, the Koopman eigenfunctions indicate the oscillatory directions of motion in the state space, and the Koopman modes describe the evolution of velocity field in the flow domain. In particular, we use the Koopman modes to study the traveling waves that appear as a result of bifurcation from steady solution to periodic and quasi-periodic flow.

The lid-driven cavity flow becomes fully chaotic at ultimately high Reynolds numbers. In such flows, the Koopman spectrum is continuous and does not contain any (non-trivial) eigenvalues. Using the properties of the Koopman operator and plausible assumptions on the post-transient dynamics, one can show how the measurements of observable on the chaotic flow can be interpreted as a realization of a wide-sense stationary stochastic process. This observation allows us to use the techniques from random signal processing to compute the continuous spectrum of the Koopman operator. We also study the flow regimes with mixed spectra, i.e., flows that have both discrete and continuous spectrum. In those flows, the evolution of flow observables is a mixture of quasi-periodic and chaotic motion, and the Koopman eigenfunctions help us distinguish and extract the quasi-periodic components of motion in the state space.

As of now, DMD-type algorithms are the most popular methods for computation of Koopman modes and eigenvalues. A notable extension of the original algorithms in Schmid (2010); Rowley et al. (2009), is the so-called Extended DMD Williams et al. (2015) which approximates the Koopman operator as a matrix using different trajectories and a dictionary of observables. Other early works have discussed the connection of DMD with other data analysis tools such as Discrete Fourier Transform (DFT) Chen et al. (2012) and linear system identification methods Tu et al. (2014). More recently, other variants of DMD have been proposed to study problems that exhibit a large range of time scales Kutz et al. (2016) and systems with external input Proctor et al. (2016). Given the increasingly large and complex data sets that are generated by simulations and experiments, DMD has also been extended to handle larger data sets (Hemati et al., 2014; Guéniat et al., 2015), different sampling techniques (Brunton et al., 2013; Tu et al., 2014; Guéniat et al., 2015) and noise (Dawson et al., 2016; Hemati et al., 2017).

In this paper, we use a different approach for computation of Koopman spectral properties. For post-transient flows, the spectrum of the Koopman operator (including both continuous and discrete components) lies on the imaginary axis, and the problem of estimating the Koopman spectrum reduces to the classical spectral estimation of signals. This problem is challenging for flows with mixed spectra where there is no a priori model for the continuous spectrum. Our methodology for Koopman spectral estimation consists of three steps: first, we apply a high-resolution algorithm - adapted from Laskar Laskar (1990); Laskar et al. (1992)- to detect the candidate discrete Koopman frequencies and modes. The Laskar algorithm provides a controllable balance between accuracy and computational efficiency which makes it suitable for large data sets like high-resolution flow snapshots. Moreover, it makes direct use of the harmonic averaging Mezić and Banaszuk (2004); Mezić (2005) which has proven convergence properties for computation of Koopman modes. In the second step, we use the ergodic properties of the attractor to discard the spurious frequencies that are artifacts of the continuous spectra. After extracting the periodic components of the flow, we estimate the continuous Koopman spectrum by applying the Welch method Welch (1967) to the chaotic residual. Our computational approach is advantageous over DMD-type algorithms since it can handle the continuous spectrum, and it is related to the well-studied techniques and notions in spectral analysis of signals.

A key objective of modal decomposition techniques is to obtain low-dimensional representation of the data from experiments or numerical simulations. Therefore, an important question regarding any decomposition is how efficiently it can capture the flow evolution. Several authors have already proposed variations of DMD algorithm to obtain low-dimensional description of the the flow features in an optimal manner Chen et al. (2012); Wynn et al. (2013); Jovanović et al. (2014). DMD is also used in a data assimilation approach to obtain a low-dimensional dynamic model of the cylinder wake flow Tissot et al. (2014). Here, we study the efficiency of the Koopman modes by considering the error in the low-dimensional truncations of KMD in representing the cavity flow dynamics. We also compare the performance of Koopman modes with the modes obtained by Proper Orthogonal Decomposition (POD).

The outline of this paper is as follows. In §II, we review the basics of the Koopman operator theory and describe how the Koopman spectrum is related to the geometry of the attractor. We also point out the connection between the Koopman mode decomposition of different observables such as stream function, velocity field and vorticity. In §III, the flow settings for the lid-driven cavity and its numerical solution are described. Section IV discusses the problem of estimating the Koopman spectrum and modes. In §V, we present the results of Koopman spectral analysis for the cavity flow. The Koopman modes and eigenfunctions are discussed in §V.2 and §V.3, and the comparison with POD is presented in §V.4. We summarize the results and conclude in §VI.

II Koopman Operator Theory

The Koopman operator theory is a mathematical formalism that relates the observations on a system to its underlying state-space dynamics. For a viscous incompressible flow, the state space is infinite-dimensional, i.e., it consists of all divergence-free smooth vector-fields defined on the flow domain. Studying the trajectories in the state space of such system via classical tools, like Poincaré maps, is nonviable, since they involve computation or visualization in an appropriate truncation of an infinite-dimensional space. The Koopman operator viewpoint circumvents this problem by focusing on the time evolution of observables rather than state variables. For example, the pressure magnitude at a certain point in the flow domain or the total kinetic energy of the fluid are are two observables on a flow that can be analyzed through the Koopman operator framework. Observables could be multiple-valued as well, like a vector containing the values of velocity at multiple grid points, or even a field of observables such as the vorticity field. One prominent outcome of the Koopman operator theory is the Koopman mode decomposition (Mezić, 2005), which describes the evolution of such observables as a linear combination of Koopman modes, Koopman eigenfunctions and Koopman frequencies, which are all explained below. We briefly review the basic formalism of the theory but the interested reader is referred to the review article Mezić (2013) and the references therein for a more detailed exposition.

Consider the state space of a flow including all the smooth divergence-free velocity fields defined on the flow domain. The state of the flow, realized by the velocity field $\mathbf{u}$, evolves in time according to the Navier-Stokes equation written as

[TABLE]

We let $g$ be a complex-valued function on the state space of the flow, i.e., for every state $\mathbf{u}$, the observable returns the complex value $g(\mathbf{u})$. We call $g$ an observable on the flow. The Koopman operator describes how this observable changes with time. More precisely, if we assume that the solution to above equation exists and it is unique, then the Koopman operator at time $\tau\in[0,\infty)$, denoted by $U^{\tau}$, maps the function $g$ to a new function $g^{\tau}$ such that

[TABLE]

The Koopman operator is a linear operator by definition ($U(\alpha_{1}g+\alpha_{2}h)=\alpha_{1}U^{\tau}g+\alpha_{2}U^{\tau}h$ for scalars $\alpha_{1,2}$) and therefore analyzing its spectrum and eigenfunctions gives a comprehensive understanding of its action on observables. An eigenfunction of the Koopman operator is a function on the state space of the flow, similar to $g$, which evolves linearly with time. Let us denote by $\phi_{j}$ the Koopman eigenfunction associated with the Koopman eigenvalue $\lambda_{j}$. Then

[TABLE]

In this work, we consider the post-transient flow dynamics, for which, the Koopman eigenvalues are known to lie on the imaginary axis (Mezić, 2005). Therefore, we will be interested in Koopman frequencies, $\omega_{j}$, related to Koopman eigenvalues through the following,

[TABLE]

Let us temporarily assume that all the observables lie in the linear span of the Koopman eigenfunctions. Then any observable like $g$ can be expanded in the Koopman eigenfunctions,

[TABLE]

where the scalar coefficient $g_{j}$ is given by the projection of observable $g$ onto the Koopman eigenfunction $\phi_{j}$. Since the Koopman operator is linear, we can use Eq. (3) and find a new expression for evolution of observable $g$ in terms of the Koopman eigenfunctions,

[TABLE]

If we replace the single-valued observable $g$ with a vector-valued observable such as $\mathbf{g}$, and follow the above procedure, the coefficient $g_{j}$ turns into the vector of coefficients $\mathbf{g}_{j}$, and we obtain a similar expansion in the vector form,

[TABLE]

This expansion of observables in terms of Koopman eigenfunctions is the so-called Koopman Mode Decomposition (KMD). The vector $\mathbf{g}_{j}$, called the Koopman mode associated with the Koopman frequency $\omega_{j}$, describes the components of the observable $\mathbf{g}$ obtained by projection of the observable onto the Koopman eigenfunction $\phi_{j}$. As a result, the evolution of $\mathbf{g}$ in time could be described as a linear combination of Koopman modes with oscillating coefficients. We will further explain the nature of this expansion in §II.1 and Koopman modes in §V.

The above expansion can be applied to fields of observables as well, in which case, the Koopman modes become fields of coefficients. For example, projecting the velocity field observable onto a Koopman eigenfunction returns a field of coefficients which can be thought of as a steady velocity field. Note that the velocity field undergoes nonlinear time evolution described by Navier-Stokes equations, but at the same time, the Koopman mode decomposition of the velocity field as an observable, offers a linear expansion in Koopman modes. This seeming paradox between the nonlinear evolution and linear expansion of KMD is resolved once we recall that the expansion in (7) is essentially infinite dimensional and therefore it can describe the nonlinear time evolution. In the following sections, we explain how the expansion above is related to the asymptotic dynamics of the trajectories in the state space and also remark on the choice of observable for study of the cavity flow.

II.1 Flow bifurcations and Koopman mode decomposition

It is interesting to see how the expansion in Eq. (7) changes as the flow undergoes bifurcation. We will be interested in detecting the post-transient flow dynamics which is directly related to the type of attractor on which the flow trajectory is evolving. The bifurcations affect the Koopman eigenvalues, eigenfunctions and modes in the KMD, but here, we only discuss how the change in the distribution of Koopman eigenvalues can be traced back to qualitative changes in the attractor of the flow.

First, we recall some standard notions from dynamical systems theory. We call a compact invariant subset of the state space, denoted by $A$, an attractor of the dynamical system, if for many initial conditions the systems evolves toward $A$. Moreover $A$ is a minimal set in the sense that it cannot be split into smaller attractors (see e.g. Wiggins (2003)). Simple examples of attractors in the state space of flows include stable fixed points and periodic orbits which correspond to steady and time-periodic flows, respectively. The attractor could be more complicated and exhibit chaos such as the butterfly-shaped attractor of the chaotic Lorenz system. We assume that the dynamics on the attractor preserves a physical measure (i.e. a distribution), which we denote by $\mu$. Roughly speaking, this implies that the time-averages (and therefore statistical properties) of continuous observables on the flow are well-defined. Now we let $\mathcal{H}:=L^{2}(A,\mu)$ be the Hilbert space of square-integrable observables defined on the attractor ($A$) with respect to measure $\mu$. In this work, we are interested in observables that belong to $\mathcal{H}$ (which includes continuous observables as well). It turns out that Koopman operator defined in (2) is a unitary operator in $\mathcal{H}$ (i.e. its adjoint and inverse are the same), which implies that its spectrum lies on the unit circle. In the following, we use the symbol $\overset{\mu}{=}$ to describe the functional equalities, i.e., the functions on different sides of $\overset{\mu}{=}$ are equal everywhere on $A$ except on a set with zero $\mu$-measure. We also use $<f,g>_{\mathcal{H}}$ to denote the inner products in $\mathcal{H}$, i.e.,

[TABLE]

Let us revisit the key assumption that led to derivation of (7), that is, the Koopman eigenfunctions span the space of observables, in this case, $\mathcal{H}$. For simple attractors like limit cycles and torus, this assumption holds and the expansion in (7) can be used to explain the behavior of all observables. We first describe this case in more detail and then turn to the more general form of Koopman expansion for more complex dynamics.

When the trajectory in the state space of the flow evolves on a limit cycle or a torus, the post-transient flow shows (quasi-)periodic time dependence. Let $\bm{\Omega}=[\omega_{1},\omega_{2},\ldots,\omega_{m}]^{T}$ denote the vector of basic frequencies for the for the motion of state variable on an $m-$dimensional torus (for limit cycles $m=1$). The Koopman spectral expansion for the flow is given by (Mezić (2017))

[TABLE]

We have dropped the dependence of $\mathbf{g}$ and $\phi_{k}$ on the state $\mathbf{u}$ to simplify the notation. The above equation is a functional equality which holds almost everywhere on the attractor. We can evaluate it for a single trajectory starting from the initial state $\mathbf{u}_{0}$ to obtain the vector expansion

[TABLE]

The term $U^{\tau}\mathbf{g}(\mathbf{u}_{0})$ is the signal generated by observing $\mathbf{g}$ over the trajectory starting at $\mathbf{u}_{0}$. If the attractor is a limit cycle, this signal is time-periodic and (10) is simply the Fourier series expansion in time. If the attractor is a torus, this expansion is a generalized Fourier expansion for the quasi-periodic signal that is generated by measuring $\mathbf{g}$. We observe that the Koopman frequencies in the above expansion form a lattice on the frequency axis. For limit cycling systems, this lattice consists of multiples of the basic frequency $\omega_{1}$, while for the torus attractors, it is given linear combinations of the basic frequencies in $\bm{\Omega}$ with integer coefficients. Hence, a bifurcation from a limit cycle to a torus can be easily detected by counting the number of basic frequencies in the lattice of Koopman frequencies obtained from the data.

The point-evaluated expansion in (10) is more suitable for the study of fluid flows than the function expansion in (9). This is due to the fact that each flow simulation or experiment provides us with only a single trajectory in the state space and direct evaluation of the Koopman eigenfunctions on arbitrary regions of state space is not practical. In case of post-transient flows, however, the ergodicity condition - which is discussed later - allows us to construct and visualize the Koopman eigenfunctions on the attractor using the signals coming from as few as one trajectory. We will use this fact to construct and visualize the eigenfunctions in §V.

The converse of the above statements is also true, that is, if the Koopman spectrum of observables has only a countable number of frequencies, then the flow trajectory must be evolving on a torus-shaped attractor in the state space. In fact, the representation theorem from the ergodic theory states that if the post-transient flow dynamics is ergodic and smooth, the Koopman operator having only discrete spectrum implies that the motion in the state space is topologically equivalent to rotation on a torus Neumann (1932). This classic result combined with numerical KMD algorithm gives a practical framework for detecting motion on tori in high-dimensional systems.

For post-transient flows with chaotic behavior, the Koopman eigenfunctions do not span $\mathcal{H}$ and evolution of observables cannot be described based on them. In fact, the Koopman operator spectra, in addition to eigenvalues, includes continuous spectrum which is related to the chaotic component of the flow. The spectral expansion for the Koopman operator takes a more general form (see e.g. MacCluer (2008)), however, as first stated in (Mezić, 2005), we can still represent it in a way that distinguishes the quasi-periodic and chaotic components of the evolution. For the scalar observable $g$, it can be written as

[TABLE]

The first term on the right-hand-side is the contribution of discrete spectrum and describes the quasi-periodic part of the flow (similar to (9)). The second term is the contribution of the continuous spectrum. Informally speaking, $i\alpha$ with $\alpha\in(-\infty,\infty)$, denotes a continuum of eigenvalues distributed along the imaginary axis. The term $dE_{\alpha}(\cdot)$ is the spectral measure of the Koopman operator, that is, for each interval of frequencies such as $I=[\alpha_{1},\alpha_{2}]$, $\int_{\alpha\in I}dE_{\alpha}(g)$ is the projection of the observable $g$ onto the eigen-subspace of $\mathcal{H}$ associated with $I$. The above expansion in the functional form is not suitable for flow applications, and it can be converted to a scalar equality by taking its inner product with the same observable $g$. That is

[TABLE]

where we have assumed that eigenfunctions are normalized, i.e., $\|\phi_{k}\|_{\mathcal{H}}=1$. In passing from (11) to (12), we have made a technical assumption that the spectral measure of the Koopman operator for the chaotic part is absolutely continuous. The Koopman spectral density $\rho_{g}$ denotes the contribution of the continuous spectrum, such that the contribution of the frequency interval $I$ to the evolution of $g$ is given by

[TABLE]

In order to compute the spectrum of the Koopman operator from the flow data, we need to assume that the post-transient dynamics is ergodic. This implies that the statistics of the flow is independent of the initial condition, and the trajectories starting almost everywhere provide a perfect sampling of observables (in the sense defined in (14) below). The ergodicity assumption holds for post-transient evolution of typical dynamical systems, including systems with periodic and quasi-periodic attractors and many chaotic systems like Lorenz Luzzatto et al. (2005). Under this condition, we can use the pointwise ergodic theorem Petersen (1989) to approximate the inner product in (12) from the data,

[TABLE]

where $r_{g}(\tau)$ is the autocovariance function of $g$ at time $\tau$. Therefore, we can approximate the spectral density of the Koopman operator by first extracting the chaotic component of $g$, then approximating $r_{g}$ using finite-time observations (i.e. finite $T$ in (14)), and finally applying inverse Fourier transform to $r_{g}$. We will discuss the practical aspects of this computation in §IV.

II.2 Stochastic processes and Koopman representation of deterministic chaos

In analyzing the chaotic data from experiments and simulations, it is customary to use the tools from applied probability theory even in the case that underlying dynamical systems is fully deterministic. The reasoning behind this approach is the duality between the post-transient evolution of dynamical systems which is measure-preserving and the stationary stochastic processes. A classic formalism of this duality can be found e.g. in Doob and Doob (1953). In this section, we reiterate this connection in the framework of the Koopman operator theory with an emphasis on the spectral expansion of observables.

Recall that a continuous-time stochastic process is a collection of real random variables that are indexed by time, and denoted by

[TABLE]

where $X_{t}$ is the random variable at time $t$ with a specified distribution over real line. A stochastic process is wide-sense stationary if it satisfies two conditions. First, its expected value should not change with time, i.e.,

[TABLE]

and second, its autocovariance function only depends on the lag time, i.e.,

[TABLE]

Now we consider the deterministic flow evolving on the attractor $A$ which preserves the normalized measure $\mu$ ($\mu(A)=1$). We see that the collection of observables

[TABLE]

is a stochastic process defined on the probability space $(A,\mu)$. Each observable $U^{t}g$ is a random variable and it assigns a probability distribution on the real line which is given by

[TABLE]

where $\mathbf{P}(B)$ is the probability of the interval $B$. Because of the measure-preserving property of the dynamics, this probability is independent of $t$, and the stochastic process in (18) is identically distributed (but not independent). Moreover, it is a wide-sense stationary process; in view of (8), we can write

[TABLE]

and

[TABLE]

where we have used the unitary property of the Koopman operator, i.e., $(U^{t})^{*}=U^{-t}$. Using the measure-preserving property, one can show (18) is strictly stationary as well Doob and Doob (1953), but that is not required for the spectral expansion.

According to the Wiener-Khintchine theorem (e.g. Peebles ), the covariance of any wide-sense stationary process, such as (18), has a spectral expansion in the following form,

[TABLE]

where $F$ is the power spectral distribution of the process. Note that this expansion holds for the general post-transient dynamics including both chaotic and quasi-periodic behavior. In case that there are no quasi-periodic components in the flow, and $F$ is absolutely continuous similar to (12), we can rewrite the above expansion as

[TABLE]

where $\rho$ is called the Power Spectral Density (PSD) of the stochastic process. Despite the deterministic nature of our system, we observe that we can treat the chaotic component of the data as a a realization of a stationary process, and consequently, the notion of the Koopman spectral density coincides with that of PSD for random signals. This observation enables us to use the spectral estimation techniques of stochastic signals for computation of Koopman continuous spectrum.

II.3 Choice of observables and the relationship between their Koopman modes

In this section, we consider the choice of observables for KMD and the relation between their modal decomposition. This question is important since applying KMD to an observable reveals only the Koopman eigenvalues that are present in the expansion of that observable. Furthermore, one can use the relationship between the Koopman modes of different observables to reduce the computational cost of the analysis. The propositions in Appendix A assert that if two observables are related through a linear operator, then their modes are also related via the same linear operator. For example, consider the field of stream function $\psi$ and the velocity field $\mathbf{u}$ in an incompressible 2D flow. These two observables are related thorough the linear operator $\nabla^{\perp}:=[\partial/\partial y,-\partial/\partial x]^{T}$, that is, $\mathbf{u}=\nabla^{\perp}\psi$. Let $\psi_{j}$ and $\mathbf{u}_{j}$ denote the Koopman modes of these two observable fields associated with Koopman eigenvalue $\lambda_{j}$, then

[TABLE]

A similar relationship could be established between the Koopman modes of the vorticity field, denoted by $\bm{\zeta}_{j}$, and those of the velocity field,

[TABLE]

This further implies that applying KMD to either of the these observable fields yields the same Koopman eigenvalues as long as none of the modes lie in the null space of the linear operator.

The knowledge of any of the above observable fields, i.e., stream function, velocity field or vorticity, uniquely determines the state of the system and therefore it can be used to elicit the Koopman spectrum of all other observables of interest. Thus, we conclude that applying KMD to any of these fields would give us the information which is sufficient to detect the flow bifurcations. In the dynamical analysis of the cavity flow, we choose the stream function as the primary observable for the application of KMD since its Koopman modes and eigenvalues are least expensive to compute. The Koopman modes of velocity and vorticity can be computed using (21) and (22).

III The lid-driven cavity flow

The 2D lid-driven cavity flow is a simple model of an incompressible viscous fluid confined to a rectangular box with a moving lid. This flow is usually used as a benchmark for numerical simulations, and represents a simplified model of geophysical flows driven by shear Tseng and Ferziger (2001); Gildor et al. (2010), and the flow inside a common type of mixer in polymer engineering Chella and Ottino (1985). The 2D cavity flow is also realized in experiments using soap films Gharib and Derango (1989) (for experiments on 3D flow see Koseff and Street (1984)). This flow is particularly interesting for the Koopman analysis because it shows a wide range of dynamic behavior depending on the increase of the top lid velocity Shen (1991); Poliashenko and Aidun (1995); Cazemier et al. (1998); Auteri et al. (2002); Peng et al. (2003); Balajewicz et al. (2013).

Our computational model of the flow consists of a square domain $[-1,1]^{2}$, with solid stationary boundaries, except the top lid (at $y=1$) which moves with a regularized velocity profile,

[TABLE]

This boundary condition has a low-order polynomial form which satisfies the continuity and incompressibility in the top corners (as opposed to the uniform velocity profile), and it is frequently used in numerical studies on cavity flow (see e.g. Shen, 1991; Botella, 1997; Balajewicz et al., 2013).

The incompressibility of the flow allows us to use the stream function formulation of the Navier-Stokes equation,

[TABLE]

subject to two types of boundary condition on the stream function,

[TABLE]

where the wall velocity $u_{w}$ is zero everywhere except at the top wall, where $u_{w}(y=1)=u_{lid}$. The solution of the cavity flow as described above is known to exist and be unique, and moreover, the flow trajectory asymptotically converges onto a universal attractor in the state space (Temam, 1988).

For numerical solution, we have used the Chebyshev-spectral collocation method described in Ref. Trefethen, 2000. The stream function is approximated by a polynomial of order $N$ in spatial directions. This polynomial is determined by its values at the Chebyshev points,

[TABLE]

Given the polynomial approximation and the prescribed boundary condition in (25), we use the transformed variable $q(x,y)$ defined by

[TABLE]

which satisfies the Dirichlet boundary condition identically, and turns the Neumann boundary condition into Dirichlet boundary condition, i.e.,

[TABLE]

For the temporal discretization of the ordinary differential equations on $q(x_{i},y_{i})$, we have used the second-order Crank-Nicholson scheme for the diffusion terms and second-order Adams-Bashforth discretization for the convection terms. The flow solutions studied in this work are computed using zero initial velocity. The numerical solutions of the steady flow obtained by our method agree with the results reported in Ref. Shen, 1991. There is also agreement on the time periods of the periodic flows between the two studies. To the best of our knowledge, however, there are no reported benchmark solutions for quasi-periodic or aperiodic flow.

IV Numerical computation of Koopman spectrum and modes

As discussed in section II, the Koopman spectrum of post-transient flows lies on the imaginary axis, and its estimation reduces to the classical spectral analysis of flow signals. In this work, we are specially interested in flows that possess a continuous spectrum in addition to discrete frequencies. Reliable estimation of each of these two components from data has a rich history in the context of signal processing and is still a subject of ongoing research. The DFT algorithm, by itself, gives a good approximation for the location of the discrete frequencies and there are a large number of the so-called high- or super-resolution algorithms, based on DFT or otherwise, that improve the accuracy of such estimation. For continuous spectra, however, DFT is a poor estimator. Application of DFT to the autocovariance function in (14) produces an estimate of spectral density with high fluctuations that do not diminish with the increase of data samples Stoica et al. (2005). Therefore, the algorithms developed to resolve continuous spectrum use some type of local averaging over frequency domain to reduce this variance. Conversely, this averaging process reduces the frequency resolution and makes these algorithms ill-suited for detection of discrete spectra Ghil et al. (2002). As a result a judicious combination of these methods should be used for computation of mixed spectrum (i.e. including both continuous and discrete parts).

Our strategy for computing the Koopman spectrum is to first detect and extract the discrete frequencies using a high-resolution algorithm, and then apply a continuous spectra estimator to the remainder. Note that most of the developed methods for accurate estimation of mixed spectrum are parametric, in the sense that they are based on specific models for the continuous spectrum such as colored or auto-regressive noise Ghil et al. (2002); Li and Stoica (1996); Stoica et al. (1997), which are not valid for typical chaotic dynamical systems. Our methodology here is non-parametric, and besides absolute continuity (discussed in §II.1), we don’t make any assumptions on the shape of the continuous spectrum. Instead, we connect our analysis to the theory of dynamical systems through the ergodicity assumption. Namely, given that the dynamics on the attractor is ergodic, we use the fact that the Koopman modes are unique (i.e. depend only on the observable and the flow parameters) which allows us to identify and discard the spurious discrete frequencies that are not robust with respect to the choice of initial condition or the time interval of integration.

The succession of ideas in this section are as follows: first, we describe the idea of harmonic averaging from classical ergodic theory which has proven convergence properties for computation of Koopman modes given the knowledge of Koopman frequencies. Then, we discuss the Laskar algorithm for computation of the discrete spectrum and benchmark its numerical performance against other high-resolution algorithms. In the last subsection, we discuss our procedure for approximation of Koopman continuous spectrum from the chaotic component of the data, and test its performance for two well-known chaotic dynamical systems.

IV.1 Harmonic averaging and DFT

For post-transient flows, the Koopman eigenfunctions are orthogonal Mezić (2005) and the Koopman modes can be computed via direct projection of the observables onto the Koopman eigenfunctions. Let $\phi_{j}$ be the normalized Koopman eigenfunctions ($\|\phi_{j}\|=1$) associated with the frequency $\omega_{j}$. We observe that the Koopman eigenfunction evolves as $\phi^{\tau}(\mathbf{u}_{0})=e^{i\omega\tau}$ over a single trajectory of the system. Using the pointwsie ergodic theorem, we can compute the Koopman modes using the harmonic average,

[TABLE]

The above limit is known to exist for almost every initial condition under the assumption that the dynamics on the attractor is preserving a measure (Wiener and Wintner, 1941) - which is less restrictive than ergdocity. The time series obtained by experiments and simulations consists of time-discrete samples over finite intervals. Assuming uniform sampling at time instants $\{\tau_{0}=0,\tau_{1},\ldots,\tau_{N-1}=T\}$, we can approximate the harmonic average as

[TABLE]

where $\mathbf{g}(k)$ is the value of observable at the sampling time $\tau_{k}$. For any $\omega_{j}$ that is a Koopman frequency, we have $\mathbf{g}^{N}_{j}\rightarrow\mathbf{g}_{j}$ as $N\rightarrow\infty$, and otherwise $\mathbf{g}^{N}_{j}\rightarrow 0$. For periodic and quasi-periodic attractors, the rate of convergence is proportional to $N$ Mezić and Sotiropoulos (2002), but for typical chaotic systems it scales with $\sqrt{N}$ Krengel and Brunel (1985).

Given a uniform sampling in time, we can use DFT frequencies as a rough approximation of the Koopman frequencies. Let the number of samples $N$ be even, and denote the sampling interval by $\Delta\tau:=T/(N-1)$. The DFT grid of frequencies is

[TABLE]

Accordingly, computing the harmonic average in (31) reduces to computing the DFT amplitude of the observations,

[TABLE]

DFT is already shown to be equivalent to DMD when applied to a linearly independent sequence of snapshots with zero mean (Chen et al., 2012). The advantage of using DFT to find the Koopman modes lies in its relative simplicity and the availability of Fast Fourier Transform (FFT) algorithms for its implementation. On the other hand, it suffers from two basic shortcomings. The first one is the picket fencing, i.e., the Koopman frequencies depend on the dynamics and may occupy arbitrary locations on the real interval, whereas DFT frequencies are determined by the sampling rate and observation interval. The second phenomenon, known as spectral leakage, refers to the spillage of energy from a frequency to its neighborhood, and it is due to the finite length of observation interval which leads to errors in approximation of the modes Oppenheim and Buck (1999); Stoica et al. (2005). In the following, we discuss some of the methods developed to overcome these problems.

IV.2 Estimation of Koopman frequencies

The problem of detecting discrete frequencies from noisy signals is often called line spectral estimation. The general goal of line spectral estimation methods is to compute estimates of frequencies with errors smaller than the DFT frequency resolution. Many of such methods use DFT as a preliminary step because of its computational efficiency, and often utilize a combination of windowing and interpolation to reduce the leakage and fencing problem (see e.g. Jain et al. (1979); Andria et al. (1989); Agrez (2002)). Some other techniques, including Prony analysis De Prony (1795) and Nonlinear Last-Squares (NLS) method Stoica et al. (2005), treat the line spectral estimation as a data fitting problem to find the frequencies and associated amplitudes that represent the time series with least error. These methods do not face the the DFT shortcomings, but they are more costly for computation and suffer vulnerability to noise (e.g. Prony analysis) or the choice of initial guess for the optimal values of fitting (e.g. NLS) Stoica et al. (2005). We note that using Prony analysis to compute the Koopman modes is explored in Susuki and Mezić (2015).

There are also the so-called subspace techniques which exploit the linear algebraic properties of matrices that embed the signal information. The two most popular algorithm in this class are the MUSIC Schmidt (1986) and ESPRIT Roy and Kailath (1989) which use eigen-decomposition of the data covariance matrix. These methods circumvent the obstacles of the DFT-based methods by posing the frequency estimation as an eigenvalue problem, which leads to accurate estimates at a higher computational complexity due to the embedding of time-series in large matrices. Moreover, these methods are parametric and their good performance is only guaranteed when the noise follows a pre-determined model (which is usually white noise) Stoica et al. (2005). A more recent class of super-resolution algorithms recast the line spectral estimation as a convex optimization of measures on the frequency domain Candès and Fernandez-Granda (2014); Fernandez-Granda (2016). Under the two conditions of spectral sparsity and minimum separation between the frequencies, this framework recovers the exact values of frequencies from a finite number of time samples. Unfortunately, this framework has a high computational complexity and it is only suitable for discrete spectrum identification in presence of little noise.

In this work, we adapt the algorithm suggested by Laskar Laskar (1990); Laskar et al. (1992) to compute the discrete Koopman frequencies and the associated modes. This algorithm is attractive for two reasons: first, it makes explicit use of harmonic averaging which allows us to assess its convergence based on the theory. In fact, this algorithm was invented to detect chaotic motion from data in Hamiltonian systems with a moderate number of degrees of freedom, like the solar system (see Laskar (1990)). Secondly, this algorithm is related to a popular sparse approximation technique known as Orthogonal Matching Pursuit (OMP)Tropp and Gilbert (2007). OMP efficiently approximates a sparse vector (i.e. vector with few non-zero elements) given a relatively small number of linear measurements on the sparse vector through an iterative greedy algorithm. Different variants of this algorithm are frequently used in decomposition of signals and images into sinusoids, wavelets etc. (see e.g. Mallat and Zhang (1993); Wang (2006); Fannjiang and Liao (2012); Mamandipoor et al. (2016)).

The main idea in Laskar algorithm is to discretize the frequency domain and use OMP (implemented as FFT and harmonic averaging) to find the frequencies and associated amplitudes that best explain the time-sampled values of the observables (which are the linear measurements in the sense of OMP). The algorithm also uses windowing and adaptive refinement of the initial grid to diminish the effect of spectral leakage and picket fencing. The computational steps are outlined in algorithm 1, and below we describe the structure of data matrix used as the input.

Let $\{\mathbf{g}(0),\mathbf{g}(1),\ldots,\mathbf{g}({N-1})\}$ be the set of observations on the vector-valued observable $\mathbf{g}$, made on uniformly-spaced time instants $\{\tau_{0}=0,\tau_{1},\ldots,\tau_{N-1}=T\}$. The snapshot data matrix $G$ is defined as

[TABLE]

Also let $\|\cdot\|$ denote an appropriate vector norm on $\mathbf{g}(\cdot)$, and $\sigma$ denote the expected $\|\cdot\|$-norm of measurement or computation noise in the data. We denote by $S(\omega)$ the sinusoid associated with frequency $\omega$, that is

[TABLE]

We also make use of windowing functions in the general form of a weight vector:

[TABLE]

Choice of data matrix $G$ and appropriate norm: We have applied the above algorithm to the vector of the stream function ($\psi$) values at the computational grid points given in (III). The Koopman modes of velocity and vorticity are subsequently computed using (21) and (22). The results reported in this paper are computed using the sampling rate $\omega_{s}=10~{}sec^{-1}$, where $sec$ is the unit of the the characteristic time given by

[TABLE]

and $L_{R}$ and $U_{R}$ are the half of the cavity side length and maximum velocity on the top lid, respectively. Our numerical experiments show that the computed frequencies are independent of the sampling frequency $\omega_{s}\in[10,200]$. We have chosen the vector-norm in the above algorithm such that it reflects the kinetic-energy norm of the Koopman modes, that is,

[TABLE]

with $U_{R}$ denoting the maximum velocity on the top lid. We choose $\sigma^{2}=10^{-6}\|\mathbf{u}_{0}\|^{2}$. This is a heuristic choice and reflects how strong a periodic component we want to resolve. Other factors that might be considered are the accuracy of the numerical simulation and the computational cost.

Dictionary of frequencies for real-valued data and choice of the window function: Given that the spectrum is symmetric for real-valued data, we can effectively reduce the computational cost by doing the search and refinement (step 1 and 2) for a non-zero frequency $\omega_{k}$ and then add the pair $(-\omega_{k},\omega_{k})$ to the dictionary in step 3. To evaluate the filtered harmonic average in step 2, we use the Hann window given by

[TABLE]

Using the window function is not necessary but improves the detection of frequencies that are close to each other - in the case of quasi-periodic flow - as it reduces the local spectral leakage. Alternative window functions can be chosen based on the proximity and relative strength of the frequencies (see e.g. Stoica et al. (2005)).

Least-square projection and harmonic average: The least-square problem in step 4 is equivalent to orthogonal projection of observables onto the Koopman eigenfunctions. In fact, in the limit of $N\rightarrow\infty$, the computation of Koopman modes in step 4 reduces to the harmonic average in (30). To see this, note that the solution to (40) is given by $W=D^{\dagger}G^{T}$, however as $N\rightarrow\infty$, the columns of $D$ become orthogonal and $D^{\dagger}\rightarrow(1/N)D^{*}$. It is easy to check that $(1/N)D^{*}G^{T}$ yields the harmonic average of columns of $G^{T}$, i.e., the Koopman modes.

Benchmark: We compare the performance of Laskar algorithm to two other high-resolution algorithms. The first one is the Newtonized Orthogonal Matching Pursuit (NOMP) Mamandipoor et al. (2016). The benchmark study in Mamandipoor et al. (2016) suggests that NOMP is a near-optimal algorithm in the sense that its accuracy is close to the theoretical limit. This algorithm is similar to Laskar, but one of its distinctive features is the refinement of all frequency estimates after the detection of each new frequency. The computational run-time of NOMP is approximately $O(K^{2})$ longer than Laskar, where $K$ is the number of detected frequencies. We also implement the MUSIC algorithm using the rootmusic() function in MATLAB. This algorithm estimates the frequencies based on the eigen-decomposition of the data covariance matrix and has a high computational complexity which does not allow its implementation on large number of observables. In case of large data sets where these two algorithms cannot be applied to the all data, we implement them in the compact form, i.e., we choose the sampling of the stream function at 10 random points, and apply MUSIC/NOMP to compute the Koopman frequencies and then compute the Koopman modes using an orthogonal projection onto the dictionary of sinusoids.

Spectral estimation by MUSIC and NOMP do not suffer from the fencing problem since they don’t use the discretization of the frequency domain, and therefore provide accurate estimates of Koopman frequencies for (purely) periodic and quasi-periodic flows. The accuracy of Laskar algorithm, on the other hand, depends on the choice of adaptive frequency grid, and higher resolutions requires evaluating (37) over finer grids. However, comparison (in FIG. 1) shows that Laskar algorithm is better suited for computing the Koopman frequencies and modes in case of large flow data with mixed spectrum. The reason for this is two-fold: Many super-resolution algorithms are designed based on a special model for the noise spectrum (e.g., the MUSIC algorithm relies on the white noise model), which does not hold for the continuous spectrum of dynamical systems. (In the context of detecting discrete frequencies, we treat the continuous spectrum as noise.) The second reason is related to the computational cost of the algorithm. The super-resolution algorithms that are not based on FFT are computationally expensive, and they can only be applied in the compact form, (i.e. applying to one or few observables simultaneously). In such cases, they may fail to capture many of low-energy frequencies in the presence of noise, as in the case of mixed spectra shown in FIG. 1. (We have verified those frequencies by the criteria introduced in the next section and the fact that they lie on the lattice of frequencies described in (9)). In contrast, the computational parsimony of the Laskar algorithm allows estimation of the frequencies using a larger set of (and possibly all) the observables. This, in turn, increases the effective signal-to-noise ratio for the low-energy periodic components which leads to detection of their associated frequencies.

IV.3 Detection of the spurious frequencies in flows with mixed spectra

In the case of mixed spectra, application of high-resolution methods to the data might produce peaks that are not genuine Koopman frequencies, but artifacts of the continuous spectrum. In such flows, we need a criteria to distinguish such peaks from the actual frequencies. On the other hand, the assumption of the ergodicity implies that the Koopman modes are unique (i.e. don’t depend on the initial condition), and therefore Koopman modes computed over different (and sufficiently long) intervals should be the same. To use this notion we run Laskar’a algorithm on different (overlapping or non-overlapping) chunks of the snapshot matrix. We discard the frequencies whose associated modes show too much variability depending on the time interval of computation. In the results to be presented, we have discarded the modes that show more than %5 variability in the kinetic energy norm, while the modes are computed over intervals of 1000 $sec$ and longer.

IV.4 Estimation of Koopman continuous spectrum

Recall from §II.2, that the spectral density of the Koopman operator appearing in (12) coincides with the PSD of the chaotic component of signals generated by measuring observables. If $g$ is a real-valued observable with purely continuous spectrum, then we can approximate the autocovariance function in (14) using the time series data

[TABLE]

and compute the correlogram approximation of PSD as

[TABLE]

Such a direct evaluation produces a highly fluctuating estimate of $\rho_{g}$, and the computation must be modified to get a more reliable estimate Stoica and Nehorai (1989).

We use the Welch method Welch (1967) to approximate the continuous part of the Koopman spectrum. The idea behind this algorithm is simple: it approximates the spectral density of the signal over small (and possibly overlapping) subsamples of the data using FFT, and averages the computed densities over all those subsamples. The averaging process reduces the variance of PSD estimation by a factor that is equal to the number of subsamples Welch (1967). This reduction in the variance comes at the price of low spectral resolution (which increases with the length of subsamples), and using too many subsamples may result in over-averaging and getting a flat spectrum. Therefore, the number and length of windows should be chosen carefully to maintain an accurate estimate while resolving the distribution of energy over frequencies.

The Welch method reduces to the Bartlett method Bartlett (1950) in case of non-overlapping subsamples. Bartlett was among the first to realize that different subsamples of the data can be averaged to find a better estimate of the spectral density given that the autocovariance decays rapidly enough and the subsamples are sufficiently large. This methodology can be interpreted as special case in the well-known class of Balckman-Tukey estimators Blackman and Tukey (1958) and the general class of filter bank approaches. We refer the reader to Ref. Stoica et al., 2005 for a discussion of connections between these methods.

We apply the Welch method, outlined in algorithm 2, to the chaotic component of the velocity field. This component is computed by extracting the contribution of Koopman modes from the original data matrix. To measure the contribution of the continuous spectrum to the whole flow field, we compute the kinetic energy density of the continuous spectrum given by

[TABLE]

This definition would allows us to compute the kinetic energy content of each frequency interval via integrating $p(\omega)$ over that interval, i.e.,

[TABLE]

and we will recover the average kinetic energy of chaotic fluctuations by calculating $P([0,2\pi))$.

We test the Welch algorithm using two well-known chaotic dynamical systems. The first one is a discrete-time map on a periodic 2D domain, known as Arnold’s cat map, and given by

[TABLE]

For the choice of observables

[TABLE]

the Koopman spectral density is known in analytical form GOVINDARAJAN et al. (2017),

[TABLE]

where $\theta\in[0,2\pi)$ is the discrete-time frequency.

The second system that we consider is the chaotic Lorenz system:

[TABLE]

This system is known to have only continuous spectrum (except the zero frequency) Luzzatto et al. (2005), but no analytical expression exists for the spectral density of non-trivial observables. We compare the Welch estimation of the Lorenz spectrum with the recent results in Korda et al. (2017) which is based on the approximation of Fourier moments of the spectral measure and the Christoffel-Darboux kernel. In particular, we consider the observable

[TABLE]

and compute its spectral density $\rho(g_{3};\omega)$ where $\omega\in[0,\omega_{s}/2)$ is the continuous time frequency and sampling frequency $\omega_{s}$ is $10\pi$. The comparison in FIG. 2 shows great agreement between the results of Welch method, analytic densities of the cat map and the numerical results of Korda et al. (2017) on Lorenz system.

V Cavity flow dynamics and KMD

V.1 Koopman spectrum

Figure 3 shows the distribution of kinetic energy in the discrete and continuous spectrum of the Koopman operator. The energy contained in Koopman modes (black bars in the figure) is simply the kinetic energy contained in each mode, but the representation of energy over the continuous spectrum is slightly different: the amount of energy contained at each frequency interval is the integral of the kinetic energy density of continuous spectrum (defined in (47) and shown as the dashed blue curve) over that interval.

The evolution of the Koopman spectrum in FIG. 3 indicates the sequence of the bifurcations as follows: for $Re\leq 10000$, the cavity flow induced by regularized lid velocity converges to a steady laminar solution which corresponds to a fixed point in the state space of the flow. The Koopman mode expansion for steady flow is trivial (hence not shown) and consists of zero frequency with an associated mode which is the steady flow. At a Reynolds number slightly above 10000, the steady solution becomes unstable and the numerical solution converges to a time-periodic flow which maintains stability up to $Re=15000$. The kinetic energy in this range is fully distributed in the Koopman modes. The basic frequency of periodic flow decreases with the Reynolds number, until at $Re\geq 15000$, another bifurcation occurs and the solution converges to a quasi-periodic flow. The basic frequencies of the quasi-periodic flow also decrease with the $Re$, but around $Re=18000$ another bifurcation occurs and the level of kinetic energy lying in continuous spectrum quickly rises to a few percent. This kinetic energy of continuous spectra keeps rising such that at $Re\geq 22000$ we cannot detect any robust Koopman modes which indicates there are no quasi-periodic components in the state space dynamics.

The discrete Koopman frequencies obtained for periodic, quasi-periodic, and interestingly, the mixed-spectra flow match the lattice structure of the frequencies in the KMD of quasi-periodic flow (9). That is, every frequency is accurately described by the integer combination of one or two basic frequencies (see table 1 in Appendix B). From the representation theorem mentioned in §II, we recall that this means the attractor is shaped like a limit cycle or a torus in the state space. For the flow with the mixed spectra, however, no such theorem exists but we can speculate that it consists of both a quasi-periodic factor and a chaotic factor. This type of attractor is called skew-periodic in the literature of dynamical systems theory Broer and Takens (1993).

The evolution of the Koopman spectrum in FIG. 3 offers a picture of transition to chaos that is consistent with the theory of Ruelle and Takens (1971). According to this theory, the chaotic state of the flow can be reached after one or two Hopf bifurcations from an initially stable steady flow. Physical evidence for this theory appeared in the experiments on rotating Couette flow and natural convection by Swinney and Gollub (1978). In particular, they detected the flow bifurcations using the power spectrum of velocity measurements at a single point in the flow domain. The transition to chaos was marked by the sudden growth of “background noise” in the power spectrum of the quasi-periodic flow. The above results show that the Koopman spectrum can be used as a generalized spectral tool for study of bifurcations; it offers a clear quantification of the energy in terms of true periodicity and contribution of continuous spectra for deterministic flows, and furthermore it connects the discrete spectrum to the state space geometry and flow domain.

V.2 Koopman eigenfunctions

In this section, we discuss the relationship between the Koopman spectrum and the state space dynamics. This relationship is realized through the notion of Koopman eigenfunctions, which are associated with the Koopman eigenvalues. For post-transient flows, the eigenfunctions provide an intrinsic coordinate on the state space along which the time evolution is linear oscillation. First, we construct the Koopman eigenfunctions for the quasi-periodic cavity flow using the theory presented in Mezić (2017); Mezić and Banaszuk (2004), and then we discuss its application to the flow with mixed spectra. Note that for flows with ultimately high Reynolds (e.g. $Re=30000$), the only Koopman frequency is zero and there are no oscillatory Koopman eigenfunctions and modes.

Consider the quasi-periodic flow at $Re=16000$. This flow possesses a torus-shaped attractor and the state space trajectory evolves on this torus with two basic frequencies $\omega_{1}$ and $\omega_{2}$ (whose non-dimensional values are given in Appendix B). We can parameterize the torus using two time-linear coordinates, that is $(\theta,\gamma)\in[0,2\pi)^{2}$ with the linear evolution equation

[TABLE]

The evolution of the trajectories on the actual torus in the state space is nonlinear, but the tuple $(\theta,\gamma)$ are angular coordinates along a torus with uniform flow which is dynamically equivalent to the actual torus in the state space Mezić (2017). On this time-linear coordinates the Koopman eigenfunctions are the same as the Fourier functions, i.e., the Koopman eigenfunction $\phi_{k,l}$ associated with frequency $\omega_{k,l}=k\omega_{1}+l\omega_{2}$ is

[TABLE]

Although there is no analytical formula for transformation from the actual attractor to the time-linear coordinates defined above, we can construct the eigenfunctions in the state space using information on a trajectory, given that the trajectory is ergodic and sufficiently long. We can normalize an the eigenfunctions such that $\phi^{t=0}(\mathbf{u})=1$, and hence the value of the eigenfunction along the trajectory is given by setting $\theta=\omega_{1}t$ and $\gamma=\omega_{2}t$ in (49).

Figure 4 shows the construction of the eigenfunctions in the state space of the flow. The state space is realized by delay embedding of some typical observables Takens (1981) - in this case the stream function at random points in the flow domain. The attractor of the periodic flow is a limit cycle (top row in the figure) and the Koopman eigenfunctions (shown as color field) correspond to the one dimensional linear time coordinate ($\theta\in[0,2\pi)$). For the quasi-periodic flow, the attractor is a 2-torus and the Koopman eigenfunctions show the directions on the torus where the evolution is linear and periodic, e.g., the eigenfunction $\phi_{0,1}$, shown in the rightmost panel of the second row shows the coordinate $\gamma$ along which the trajectories oscillate with frequency $\omega_{2}$.

For skew-periodic attractors (i.e. flow with mixed spectra) the eigenfunctions are even more interesting because they provide coordinates on an attractor which is not exactly a torus, but possesses directions with periodic motion. The embedded attractor of the flow at $Re=19000$, for example, is similar to a torus which is related to the fact that this flow possesses a strong discrete spectrum (in the energy sense) with two basic frequencies and a relatively weak continuous spectrum. In fact, using the Koopman eigenfunctions, we can compute the factors (i.e. geometric slice) of such an attractor, on which, the motion is purely quasi-periodic. The existence of such factorization for systems with discrete Koopman eigenvalues was shown in Ref. Mezić and Banaszuk, 2004. Here, we use this idea to reconstruct the quasi-periodic component of the attractor at $Re=19000$. Let $E$ be the observable whose embedding is used to construct the attractor. According to the (11), this observable can be split to two components:

[TABLE]

where $E_{qp}$ denotes the component of $E$ that lies in the span of Koopman eigenfunctions (including the eigenfunction at zero frequency), and $E_{c}$ is the chaotic component that belongs to subspace associated with the continuous spectrum. By doing KMD on $E$, we can extract its Koopman modes and reconstruct the evolution quasi-periodic component $E_{qp}$ over the trajectory which is given by the first two terms in (11). The embedding of $E_{qp}$ constructs the torus which corresponds to the quasi-periodic part of the motion. As such, the general motion on the skew-periodic attractor (third row of the figure) can be decomposed into rotational motion along its quasi-periodic component (bottom row) superposed with chaotic motion in an unknown direction. We stress that the above constructions are valid for any type of state space realization as long as the data on an ergodic trajectory is available.

V.3 Koopman modes

The Koopman modes of the vorticity field associated with the basic frequencies of each flow are shown in FIG. 6. Each mode can be interpreted as the component of the vorticity field along the eigenfunction coordinates in the state space (the color field in FIG. 4). For the eigenfunction at zero frequency, this component is the mean flow (Koopman mode associated with zero frequency) and does not change in time. The oscillatory modes however are components of the vorticity field that linearly oscillate along the eigenfunction directions.

The major share of kinetic energy in all the unsteady regimes is contained within the mean flow. This mode is essentially composed of the central vortex in the flow and the corner eddies in the left corners. For fully chaotic flows, the mean flow is the only Koopman mode, and its structure is similar to the mean flow of periodic and quasi-periodic flow except for the intensification of the downstream eddy in the bottom right corner.

The oscillatory Koopman modes, on the other hand, describe the flow oscillations around the edge of the central vortex in the mean flow. To be more precise, the evolution of the mode $\mathbf{u}_{k}$ can be written as

[TABLE]

A careful examination of the figure shows that $Re(\mathbf{u}_{k})$ and $Im(\mathbf{u}_{k})$ are similar for each $k$ but appear to be shifted in the direction along the shear layer of the mean flow. This observation is related to the fact that the unsteady motion in periodic and quasi-periodic regimes corresponds to wave(s) that travel along the downstream edge of the central vortex. This traveling wave structure is observed in the previous studies of cavity flow (see e.g. the sequence of flow snapshots in Poliashenko and Aidun (1995); Auteri et al. (2002)) but never characterized.

The Koopman modes provide a straightforward framework to characterize the traveling waves from the data. Let’s consider a simple example first: let $[0,1]$ be a periodic domain, over which, the general form of traveling wave is

[TABLE]

with $f$ being $2\pi$-periodic. Using the Fourier series expansion, we have

[TABLE]

Clearly, the last expression is the KMD of $f$ with the Koopman modes given as $f_{j}(x)=\exp{-2\pi jkx}$. Having this example in mind, we can compute the wave numbers (and phase velocity) of traveling waves in the cavity flow through the following steps: First, we compute the phase of each Koopman mode given by

[TABLE]

Then, we sample the values of $\theta_{k}$ along the direction of travel, denoted by $\hat{x}$, and compute the average local slope of $\theta_{k}(\hat{x})$ to get the wave number of the mode. This process is summarized in FIG. 7 (a) for the Koopman mode $\mathbf{u}_{1}$ of the periodic flow at $Re=13000$. The results of this computation for different Koopman modes (shown in FIG. 7(b)) indicates that the Koopman modes associated with higher frequencies have proportionally higher wave numbers. This is expected from the KMD expansion of traveling wave (53). Moreover, it suggests that the wave numbers rarely change with the Reynolds number, however the phase velocity (slope of the lines in FIG. 7(b)) slightly decreases due to the decrease in the Koopman frequencies. We note that the above methodology based on Koopman modes is already used, in a concealed form, to extract the wave numbers from experimental data on nonlinear waves in thermally-driven flows Mukolobwiez et al. (1998); Garnier et al. (2003), and internal waves in stratified flows Mercier et al. (2008).

An interesting observation is that the oscillatory Koopman modes show a remarkable structural robustness in the range $Re=11000-19000$, despite the flow undergoing bifurcations from periodic to quasi-periodic and then skew-periodic. An examples of this robustness are the modes in panel 2, 7 and 12 in FIG. 6 (counting from top left to right) which correspond to the frequency trace of $St_{1,0}$ in periodic, quasi-periodic and mixed-spectra regimes, and panel 8 and 13 associated with the frequency $St_{0,2}$ in quasi-periodic and mixed-spectra regimes. This observation suggests that Koopman modes may provide a suitable basis for reduced modeling of flows (e.g. Noack et al. (2003)) over wide range of Reynolds number.

V.4 Spectral Projections and Proper Orthogonal Decomposition (POD)

We study the the efficiency of Koopman modes in representing the flow, by computing the error of the spectral-projection models. An $n-$dimensional spectral projection model, is an $n$-term truncation of the KMD where the modes are sorted based on their kinetic energy. The error defined as

[TABLE]

gives the kinetic-energy norm of the difference between the spectral projection model (the sum in the above expression) and the actual flow field. Given the finite-dimensional nature of these models, they are essentially quasi-periodic approximations of the flow. The time-averaged kinetic energy of the error for the spectral projections of order 1-10 is shown in FIG. 8. In the periodic and quasi-periodic flows, the bulk of the motion is readily captured by a few Koopman modes and the low-order projections approximate the flow with great accuracy. As the flow becomes less periodic with the increase of Reynolds number, the approximation error increases as well. For fully chaotic flows, the only Koopman mode is the mean flow and therefore there are no low-dimensional spectral projections except the steady one-dimensional model which is the mean flow itself. For purpose of comparison, however, we have plotted the error of spectral projections using Fourier modes (computed via FFT) at $Re$=30000. The kinetic energy of unsteady motion in this flow is spread in the continuous spectrum and any low-dimensional approximation using oscillatory components would involve large errors.

Figure 8 also shows an instructive comparison between Koopman mode decomposition and the Proper Orthogonal Decomposition (POD). POD is a decomposition of the flow field into spatially-orthogonal modes such that the POD-truncated models have the minimum energy error among all choices of orthogonal decompositions (Berkooz et al., 1993). Due to its optimality and advantageous numerical properties, POD has been the keystone of many studies on coherent structures and low-order modeling of complex flows, inclduing the lid-driven cavity flow Cazemier et al. (1998); Balajewicz et al. (2013).

In case of periodic and quasi-periodic flow, the low-dimensional spectral projection model gives a better approximation than POD-truncated models (first two panels in FIG. 8). This observation is not contradictory to the optimality of POD-truncated models, but due to the fact that Koopman modes are complex-valued for oscillating systems which allows for better representation of dynamics evolving on limit cycles or torus as previously suggested in Ref. Mezić, 2005 (see example 4 therein). We speculate that this observations would be valid for other flows that may exhibit strong quasi-periodic behavior over some parameter range, for example, periodically-driven flows (e.g. Lopez and Marques (2000)), rotating Couette flow Takeda (1999) and Rayleigh-Bernard convection Swinney and Gollub (1978), flow in converging-diverging channel Guzmán and Amon (1994) and behind bluff objects Karniadakis and Triantafyllou (1992).

For flows with mixed spectra, the error of approximation with spectral projection and POD is generally comparable but as the Reynolds number increases and the flow becomes more chaotic, POD-truncated models perform better. These observations support the alternative decomposition for model reduction proposed in Mezić (2005), in which, the KMD is used to extract the periodic components of the flow, and then POD is used for representation of the remaining chaotic component.

VI Conclusion

In this work, we used the spectral properties of the Koopman operator for analysis of a post-transient flow. Understanding of the asymptotic dynamics in the state space of the flow is achieved by inspecting the Koopman spectrum and eigenfunctions obtained by KMD of numerical data. In case of periodic or quasi-periodic attractors, the dominant Koopman frequencies possess a lattice-type structure whereas chaotic flows are associated with the continuous part of the Koopman spectrum. In between these two regimes, we observed flows that exhibit a combination of chaotic and quasi-periodic behavior (i.e. skew-periodic attractor) which is associated with mixed Koopman spectrum. The Koopman eigenfunctions, on the other hand, determine the directions of linear evolution on quasi-periodic and skew-periodic attractors.

We employed a new approach for computation of Koopman modes and frequencies which comprises of a high-resolution algorithm to detect the eigenvalues, and an averaging estimation of the continuous spectrum. This approach, in contrast to DMD-type algorithms, is capable of computing the continuous spectrum of the Koopman operator and detecting genuine Koopman frequencies for flows with mixed spectra.

Study of unsteady lid-driven cavity flow suggests that in case of periodic and quasi-periodic flow, a handful of Koopman modes are sufficient to represent the spatio-temporal features of the flow in a low-dimensional model. In fact, the Koopman modes offer a more efficient representation of such flows than POD modes. This is due to the fact that complex-valued Koopman modes are more suitable for describing the evolution on limit cycles and tori.

Acknowledgement and supplemental material

This research was partially supported by the ONR grant N00014-14-1-0633. H.A. is grateful to Ati Sharma and Babak Mamandipoor for valuable discussions on spectral estimation and POD analysis. The authors also thank Peter Schmid for suggesting the computation of wave numbers using Koopman modes.

The time series data and MATLAB codes for KMD can be found at https://mgroup.me.ucsb.edu/resources.

Appendix A Relationship between Koopman modes of different observables

We present two propositions regarding the Koopman mode decomposition of observables which are related through a linear operator or spatial differentiation.

Proposition 1

Let $g,h:\Omega\times M\rightarrow\mathbb{R}$ be two fields of observables defined on the flow domain $\Omega$ and the state space of the flow $M$. Assume $g(.,z),~{}h(.,z)\in\mathcal{F}(\Omega)$ for every state $z\in M$ where $\mathcal{F}(\Omega)$ denotes the space of fields on $\Omega$. Moreover, the observables are related through a linear operator $\mathcal{D}:\mathcal{F}(\Omega)\rightarrow\mathcal{F}(\Omega)$, such that $\mathcal{D}g=h$ and $\mathcal{D}$ is bounded in $\mathcal{F}$. Let $g_{i}$ and $h_{i}$ denote the Koopman modes of observable $g$ and $h$, respectively, associated with the Koopman eigenvalue $\lambda_{j}$. Then those modes are related to each other via $\mathcal{D}g_{i}=h_{i}$. In words, the Koopman mode decomposition commutes with the linear bounded operator $\mathcal{D}$.

Proof:

Let us assume for now that the dynamical system is measure-preserving which implies $\lambda_{j}=i\omega_{j},~{}\omega_{j}\in\mathbb{R}$. The Koopman modes are computed via the harmonic averaging and it follows that

[TABLE]

We have used the linearity and continuity (i.e. boundedness) of $\mathcal{D}$ in the 3rd and 4th equalities, respectively. If the system is dissipative, the Koopman modes are given by Generalized Laplace Analysis formula Mezić (2013). The above argument could be used, along with induction, to prove the statement for that case. $\square$

For the purposes of this paper, we are mostly interested in linear operators that involve spatial derivatives, such as gradient or curl. The derivative operator, however, is not bounded and therefore needs a special treatment which is given by the following proposition.

Proposition 2

Let $g,h:\Omega\times M\rightarrow\mathbb{R}$ be two fields of observables defined on the flow domain $\Omega$ and the state space of the flow denoted by $M$. Assume the observables are bounded over their domain and related to each other through

[TABLE]

where $D$ denotes the partial derivative with respect to a spatial coordinate in the flow domain $\Omega$. Let $g_{k}$ and $h_{k}$ denote the Koopman modes of observables associated with the Koopman frequencies $\omega_{k},~{}k=1,2,\ldots$. Then

[TABLE]

for every $x\in\Omega$. In words, the Koopman mode decomposition commutes with the spatial differentiation in the flow domain.

Proof:

Define the finite-time harmonic average of $g$ at the frequency $\omega_{k}$ by

[TABLE]

The Koopman mode of $g$ associated with $\omega_{k}$ is then given by Mezić (2005),

[TABLE]

The Leibniz rule implies that the spatial derivative commutes with the finite-time averaging in (59), i.e.,

[TABLE]

for any finite $T$. Fix $x$. If both $h(x)$ and $g(x)$ are integrable functions, then the existence of $g_{k}(x)$ implies the existence of $h_{k}(x)$ (Wiener and Wintner, 1941), and since the scalar functions $h^{T}_{k}(x)$ and $\mathcal{D}g_{k}^{T}(x)$ are equal up to any finite $T$, their limits as $T\rightarrow\infty$ must be equal, i.e., $h_{k}(x)=\mathcal{D}g_{k}(x)$. $\square$

Appendix B Discrete spectra of Koopman operator

Table 1 shows the Koopman frequencies of cavity flow, computed by algorithm 1 and ranked in order of decreasing kinetic energy. The basic frequency vector is $St=0.159826$ for the flow at $Re$=13000, $St=[0.155375,0.096911]$ for $Re$=16000 and $St=[0.150046,0.093778]$ for $Re=19000$. The rest of the frequencies can be described (to the sixth digit of accuracy) as a linear combination of the basic frequencies. The coefficients of the combination are given by $\mathbf{k}$ in the table.