Searching turbulence for periodic orbits with dynamic mode decomposition

TL;DR

This paper introduces a novel method using dynamic mode decomposition to efficiently identify unstable periodic orbits in turbulent flows, surpassing traditional recurrence-based techniques in accuracy and robustness.

Contribution

The authors develop a DMD-based approach to generate initial guesses for UPOs directly from turbulent data, avoiding the need for near recurrences and improving orbit detection.

Findings

Successfully extracts UPO periods from short data segments.

Finds more converged periodic orbits than recurrent flow analysis.

Identifies UPOs where traditional methods fail to converge.

Abstract

We present a new method for generating robust guesses for unstable periodic orbits (UPOs) by post-processing turbulent data using dynamic mode decomposition (DMD). The approach relies on the identification of near-neutral, repeated harmonics in the DMD eigenvalue spectrum from which both an estimate for the period of a nearby UPO and a guess for the velocity field can be constructed. In this way, the signature of a UPO can be identified in a short time series without the need for a near recurrence to occur, which is a considerable drawback to recurrent flow analysis, the current state-of-the-art. We first demonstrate the method by applying it to a known (simple) UPO and find that the period can be reliably extracted even for time windows of length one quarter of the full period. We then turn to a long turbulent trajectory, sliding an observation window through the time series and performing many DMD computations. Our approach yields many more converged periodic orbits (including multiple new solutions) than a standard recurrent flow analysis of the same data. Furthermore, it also yields converged UPOs at points where the recurrent flow analysis flagged a near recurrence but the Newton solver did not converge, suggesting that the new approach can be used alongside the old to generate improved initial guesses. Finally, we discuss some heuristics on what constitutes a "good" time window for the DMD to identify a UPO.

Tables

Table 1. Table 1: Periodic orbits found from the DMD results reported in figure 5 . The guessed period of the UPO, T g subscript 𝑇 𝑔 T_{g} , was determined from ω f ​ ( n = 2 ) subscript 𝜔 𝑓 𝑛 2 \omega_{f}(n=2) (equation 8 ). In the penultimate column, the UPOs found from DMD-generated guesses are connected to known UPOs converged in the slightly shorter L x = 1.75 ​ π subscript 𝐿 𝑥 1.75 𝜋 L_{x}=1.75\pi box listed on channelflow.org (see also the discussion at the start of §2). The final column shows a comparison to recurrent flow analysis, T ω R ​ F ​ A superscript subscript 𝑇 𝜔 𝑅 𝐹 𝐴 T_{\omega}^{RFA} indicates the period of any UPOs found with that method for start times close to that used in the DMD.

$t_s$	$T_g$	$T_{\omega }$	% error in $T_g$	channelflow $T_{\omega }$	$T_{\omega }^{RFA}$($t_s$, % error)
0	67.69	86.43	21.6	-	-
10	62.12	85.60	27.5	85.27	-
55	92.64	95.76	3.3	-	-
115	92.42	97.62	5.3	-	100.17(128,4.2)
175	84.35	77.25	9.1	-	-
180	81.47	76.54	6.4	-	-
200	75.04	78.23	4.1	-	-
235	61.38	71.83	14.5	-	71.83(235,12.3)
410	72.86	95.76	23.9	-	103.01(412,3.7)
415	84.19	104.28	23.9	-	103.01(412,3.7)
545	76.13	86.43	11.9	-	-
715	68.12	64.95	4.8	-	-
725	68.99	94.76	27.2	-	-
730	70.19	64.95	8.0	-	85.60(733,6.5)
735	71.31	88.61	17.3	87.89	85.60(733,6.5)
740	71.58	88.20	19.5	88.90	85.60(733,6.5)
780	88.73	90.71	3.3	90.31	-
815	76.19	78.23	2.6	-	-
825	79.71	70.12	13.7	-	-
880	95.14	104.28	8.8	-	-

Equations

$$R(t,T):=\frac{\Vert \mathbf{u}(t+T)-\mathbf{u}(t)\Vert }{\Vert \mathbf{u}(t)\Vert },$$

$${\mathbf{\Psi }}^t=\left[\begin{array}{cccc}\mathbf{\psi }(\mathbf{u}(t_1))& \mathbf{\psi }(\mathbf{u}(t_2))& \cdots & \mathbf{\psi }(\mathbf{u}(t_M))\end{array}\right],$$

$${\mathbf{\Psi }}^{t+\delta t}=\left[\begin{array}{cccc}\mathbf{\psi }(\mathbf{u}(t_1+\delta t))& \mathbf{\psi }(\mathbf{u}(t_2+\delta t))& \cdots & \mathbf{\psi }(\mathbf{u}(t_M+\delta t))\end{array}\right],$$

$$\hat{\mathbf{K}}:={\mathbf{\Psi }}^{t+\delta t}{\left({\mathbf{\Psi }}^t\right)}^{+}$$

$$\mathbf{\psi }(\mathbf{u}):=\mathbf{u}-{\mathbf{u}}_C+{\mathbf{\delta w}}_{ϵ}.$$

$$\hat{\mathbf{K}}{\mathbf{v}}_j=e^{{\lambda }_j\delta t}{\mathbf{v}}_j.$$

$${\phi }_{{\lambda }_j}(\mathbf{u})={\mathbf{w}}_j^H\mathbf{\psi }(\mathbf{u}).$$

$${\omega }_f(n):=\frac{2}{n(n+1)}\sum _{j=1}^n{\omega }_j,$$

$${\epsilon }_{\omega }(n):=\frac{1}{n|{\omega }_f{|}^2}\sum _{j=1}^n|{\omega }_j-j{\omega }_f{|}^2,$$

$${\mathbf{u}}_g(t)=z{\mathbf{e}}_x+\sum _{j=-n}^n{a}_j{\mathbf{v}}_j{e}^{{\lambda }_{j}t}.$$

$$J(\mathbf{a}):=\frac{1}{M}\sum_{m=0}^{M-1}\big{|}\boldsymbol{\psi}(\mathbf{u}(\mathbf{x},m\delta t))-\sum_{j=-n}^{n}a_{j}\mathbf{v}_{j}e^{m\lambda_{j}\delta t}\big{|}^{2}.$$

$${\phi }_{{\lambda }_j}^{DMD}(\mathbf{u})=({\mathbf{w}}_j^{DMD}{)}^H\mathbf{\psi }(\mathbf{u}),$$

$${\phi }_{{\lambda }_j}(\mathbf{u})=({\mathbf{w}}_j^{UPO}{)}^H\mathbf{\psi }(\mathbf{u}).$$

Full text

1 Introduction

Since the discovery of the first unstable period orbit (UPO) buried in a turbulent attractor (plane Couette flow) by Kawahara & Kida (2001), there has been a surge in interest and a large number of other periodic solutions found both in this configuration (e.g. Viswanath, 2007; Cvitanovic & Gibson, 2010) and in other canonical turbulent flows (Chandler & Kerswell, 2013; Willis et al., 2013; Lucas et al., 2017). The discovery of a large number of UPOs supports a perspective of turbulence in which the flow is viewed as a trajectory in a very high-dimensional dynamical system, wandering between unstable exact coherent structures (Kerswell, 2005; Kawahara et al., 2012). Individual UPOs can offer a good deal of insight into the physical processes sustaining the turbulent flow owing to their simple time dependence (Waleffe, 1997; Wang et al., 2007; Hall & Sherwin, 2010). Furthermore, periodic orbit theory (Cvitanović et al., 2016) suggests that statistics of the turbulence can be predicted from the UPOs if enough of them are found (see the attempts in Chandler & Kerswell, 2013; Lucas & Kerswell, 2015). Despite this increasing interest, the methods for finding UPOs are somewhat crude and have not changed significantly in the two decades since the first turbulent solutions were discovered.

The standard method for finding and converging UPOs begins with a search for near recurrences in data from a numerical simulation. In such a ‘recurrent flow analysis’, likely UPOs are flagged when the distance in state space (measured with an $L_{2}$ norm) between the present and past states drops below a threshold value (Kawahara & Kida, 2001; Viswanath, 2007; Cvitanovic & Gibson, 2010; Chandler & Kerswell, 2013). The resulting set of candidate orbits, each augmented with a guessed period from the time between two similar states, is then input into a Newton solver. The main downside of the approach is that it requires the turbulence to shadow a periodic orbit for at least one full cycle, and hence it can be increasingly ineffective as the Reynolds number is increased (Chandler & Kerswell, 2013). In addition, the sensitivity of the Newton solver to initial conditions can result in failures at the convergence stage when there may in fact be a UPO nearby. In this paper we introduce a new method based on dynamic mode decomposition (DMD) that goes some way to addressing these issues.

The DMD algorithm was originally invented by Schmid (2010) and can serve as an excellent data-driven method for performing global stability analyses in complex geometries, since it requires only raw data in the form of snapshot pairs. The output of the algorithm is a linear operator that (in a least squares sense over the input data) provides the best mapping between snapshots separated by a fixed time $\delta t$. An eigendecomposition of this operator yields a set of ‘dynamic modes’ which evolve exponentially in time, the time dependence being set by the associated eigenvalues. Crucially, this approach can identify frequencies of oscillation with corresponding periods that can be far longer than the time window over which the observations were recorded (Schmid, 2010; Rowley et al., 2009; Schmid et al., 2011). Since its invention, a number of variants of the algorithm have been proposed and DMD has found applications in areas of science and engineering beyond fluid mechanics (Jovanović et al., 2014; Williams et al., 2015; Kutz et al., 2016).

Despite its apparent linearity, DMD can also be a useful tool in strongly nonlinear flows. Our interest in applying DMD to turbulent data as a means to identify UPOs rests on the equivalence (under strict conditions) of the algorithm to the Koopman decomposition (Rowley et al., 2009; Tu et al., 2014; Williams et al., 2015) of a nonlinear dynamical system. The Koopman operator is an infinite-dimensional, linear operator that propagates observables of the state forward in time along trajectories of the dynamical system (Koopman, 1931; Mezić, 2005), and a Koopman decomposition of the state leads to a representation in which the nonlinear dynamics are expressed as linear superposition of Koopman modes (e.g. fixed velocity fields) which evolve exponentially in time. Koopman decompositions have been performed analytically for some simple nonlinear ODEs (Bagheri, 2013; Brunton et al., 2016; Rowley & Dawson, 2017) and PDEs (Page & Kerswell, 2018).

Applying a Koopman decomposition to a turbulent flow yields a representation of the state as a superposition of a set of harmonic averages and a broadband continuous spectrum (Mezic & Banaszuk, 2004; Mezić, 2005, 2013; Arbabi & Mezić, 2017), but individual simple invariant sets also possess their own local Koopman decompositions. For example, Mezic (2017) has shown that the Koopman eigenvalues for a nonlinear system collapsing onto a limit cycle consists of a set of repeated neutral harmonics (the limit cycle’s fundamental frequency and higher harmonics) and an infinite lattice of decaying eigenvalues which can be determined from linear combinations of the cycle’s Floquet multipliers. However, in systems with multiple exact coherent structures there are multiple local Koopman decompositions around each structure. These linear representations break down (are no longer convergent) at certain crossover points in state space (Page & Kerswell, 2018), a fact which impacts the ability of DMD to extract the Koopman modes associated with the expansions. That DMD no longer coincides with Koopman when the observation window includes the crossover point is likely the reason behind the poor performance of the algorithm when Bagheri (2013) applied it to the ‘spin up’ problem of flow past a cylinder. Page & Kerswell (2018) demonstrated the phenomenon explicitly in a Stuart-Landau equation and numerically along heteroclinic connections in the Navier-Stokes equations. In doing this, it was observed that DMD can extract some of the Koopman eigenvalues associated with the local Koopman decomposition around a periodic orbit – though only if the computation is performed within the ‘Koopman expansion zone’ where the local decomposition is convergent. Here, we go further and layout how DMD can be applied to turbulent data to identify and converge UPOs.

The remainder of this paper is organised as follows. In §2 we outline a procedure for performing and post-processing DMD calculations to generate guesses for periodic orbits and demonstrate its utility by applying it to a gently periodic edge state. In §3 the method is then used to search for periodic orbits close to a long turbulent trajectory and is compared to recurrent flow analysis. We examine some of the DMD guesses further in §4, and identify characteristics of “good” initial guesses. Finally, concluding remarks are provided in §5.

2 Methodology

In this paper we present a new method for identifying periodic orbits using DMD applied to short turbulent trajectories and contrast it to recurrent flow analysis, the current state-of-the-art. To do this we post-process numerical simulation data using both methods and attempt to converge the resulting UPO guesses via Newton-Krylov iteration. Throughout, the flow configuration is plane Couette flow in a minimal flow unit (Hamilton et al., 1995) non-dimensionalized by the upper/lower plate speed $U_{0}$ and the channel half-height $d$. The Reynolds number is fixed at $Re:=U_{0}d/\nu=400$, while the exact horizontal dimensions of the periodic computational domain are $(L_{x},L_{y})=(1.755\pi,1.2\pi)$ to match Kawahara & Kida (2001).

In our direct numerical simulations, the Navier-Stokes equations are solved using a fractional-step method with an implicit Crank-Nicholson scheme for the diffusive terms and a third-order Runge-Kutta scheme for the advection terms (Dubief et al., 2005). Second-order finite differences are used in all three spatial directions with grid stretching applied in the vertical. For the majority of results the grid resolution is $N_{x}\times N_{y}\times N_{z}=96\times 96\times 129$. Note that the slight difference in box size from the $L_{x}=1.75\pi$ box used by Cvitanovic & Gibson (2010) (and the solutions catalogued at channelflow.org) precludes a direct comparison with this extensive set of UPOs. Altering our box size to $1.75\pi$ indicates that some of UPOs we find are known solutions listed at channelflow.org. At our resolution we can match the periods of these known periodic orbits to within a relative error of $\sim 0.1$%, and have verified that this error can be reduced by further increasing our resolution. We make every effort to connect our results to previously discovered exact solutions (see table 1).

In a recurrent flow analysis (e.g. Chandler & Kerswell, 2013), guesses for periodic orbits are identified via local minima in an $L_{2}$ norm,

[TABLE]

which fall below some threshold $R_{\text{thresh}}$. If the minima occur at times $t_{i}$, then the guesses for periodic orbits are simply the corresponding states $\{\mathbf{u}(t_{i})\}$, each augmented with a guessed period $T_{g}=T_{i}$ – the future time at which the local minimum in $R(t,T)$ occurs. Note that in the form (1) we are comparing future states at $t+T$ with present states at $t$ in order to allow for comparison to results from DMD with observation windows starting at time $t$.

As discussed in the introduction, there are shortcomings to this approach for generating guesses for periodic orbits. For instance, given a near recurrence, the (future) time $T$ at which the minimum value of $R(t,T)$ occurs may differ significantly from the period of the nearby UPO. It is plausible that the sensitivity of Newton-Raphson to initial conditions could result in a lack of convergence even if there is indeed a structure nearby. In these situations, DMD may be a useful tool since it can be used to extract a set of modes associated with ‘background’ periodic motion, resulting in a more robust initial condition for the root-finder and a more informed estimate of the period. Perhaps more importantly, DMD does not require a near recurrence to occur before it can ‘sense’ periodic motion since the algorithm allows for the extraction of frequencies corresponding to oscillations far longer than the length of observation window (Schmid, 2010). Both of these scenarios are sketched in figure 1.

In DMD, vector observables of the state at discrete times $\{t_{j}\}$ are stored in a data matrix (Schmid, 2010; Tu et al., 2014),

[TABLE]

where the observable functionals that make up the elements of $\boldsymbol{\psi}$ are a design choice, often motivated by knowledge of the underlying dynamical system (e.g. Williams et al., 2015). A second data matrix is formed with observations made a time $\delta t$ (another design choice) later,

[TABLE]

where we refer to corresponding columns in $\boldsymbol{\Psi}^{t}$ and $\boldsymbol{\Psi}^{t+\delta t}$ as a ‘snapshot pair’. The DMD operator is then the best-fit (in a least-squares sense) linear operator that maps between the data matrices,

[TABLE]

where the superscript $+$ identifies the Moore-Penrose pseudo-invserve, accomplished via an SVD of the data matrix $\boldsymbol{\Psi}^{t}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{W}^{H}$. In reality the dimensionality of the problem requires that $\hat{\mathbf{K}}$ is never computed directly, but instead the $r\times r$ (where $r$ is the rank of $\boldsymbol{\Psi}^{t}$) matrix $\mathbf{U}^{H}\hat{\mathbf{K}}\mathbf{U}$ is used – the projection of the DMD matrix onto POD (Proper Orthogonal Decomposition) modes.

To search for periodic orbits we will be applying DMD to relatively short turbulent trajectories ($O(100)$ advective time units with $M\sim 100$ snapshot pairs), and hence the data matrix will likely be full rank $r=M$. This tends to lead to overfitting and unphysical DMD modes, a behaviour which can be avoided through a variety of techniques, e.g. by either performing a low-rank truncation or by introducing a sparsity constraint (Jovanović et al., 2014). Rather than selecting a threshold singular value or modifying the DMD algorithm, we instead attempt to find a DMD operator that can map between velocity fluctuations above a certain threshold by adding a small amount of white noise to our observable vector,

[TABLE]

where $\mathbf{u}_{C}:=z\mathbf{e}_{x}$ is the laminar base state and $\boldsymbol{\delta w}_{\epsilon}$ is the white noise of amplitude $\epsilon$. The minimum singular value $\sigma_{min}$ scales proportionally with $\epsilon$. Throughout, we set $\epsilon=10^{-7}$, though we have verified that our results are robust under modest changes in $\epsilon$ (plus or minus a decade).

Once the (projected) DMD operator is computed, an eigendecomposition can be performed. Under certain conditions (see Williams et al., 2015; Rowley & Dawson, 2017; Page & Kerswell, 2018) the output of DMD coincides with a Koopman analysis. In these cases, approximate Koopman modes for the observable $\boldsymbol{\psi}(\mathbf{u})$ can be obtained as right eigenvectors of $\hat{\mathbf{K}}$ (DMD modes),

[TABLE]

Note that if $\mathbf{v}^{\prime}$ is a right eigenvector of $\mathbf{U}^{H}\hat{\mathbf{K}}\mathbf{U}$, the DMD mode is obtained from $\mathbf{v}=\mathbf{U}\mathbf{v}^{\prime}$ (Schmid, 2010). Koopman eigenfunctions are related to the left eigenvectors of $\hat{\mathbf{K}}$, $\{\mathbf{w}_{j}\}$, via

[TABLE]

2.1 Spotting a periodic orbit

In a dynamical system with a single attracting periodic orbit, Mezic (2017) has shown that the spectrum of the associated Koopman operator consists of a set of repeated harmonics (multiples of the fundamental frequency of the limit cycle) along with an infinite lattice of decaying eigenvalues that can be computed from the Floquet multipliers. In the turbulent Couette flow considered in this work, all periodic orbits are unstable and each orbit will have an associated locally convergent Koopman decomposition (Page & Kerswell, 2018) of a similar form to that derived by Mezic (2017) but with unstable eigenvalues related to the set of unstable Floquet multipliers. While DMD on a trajectory shadowing the UPO for (potentially less than) one cycle is unlikely to be able to identify the stability properties of the nearby structure, the results presented in Page & Kerswell (2018) indicate that the signature of the periodic orbit – its neutral harmonics – in the DMD results can remain even for relatively short time windows.

In previous studies involving DMD on flows collapsing onto limit cycles (e.g. Bagheri, 2013; Page & Kerswell, 2018), identifying the underlying period in the DMD eigenvalue spectrum is straightfoward either because the period is known beforehand or the DMD can be run over an indefinitely long time horizon due to the stability of the exact coherent structure. However, in a turbulent flow it is unknown a-priori at any given time whether we are shadowing a periodic orbit, and if so what its period is. Therefore, given an eigenvalue spectrum from a DMD calculation on a segment of a turbulent orbit we would like to both (i) estimate the period of a nearby UPO (if one exists) and (ii) assess the quality of our “guess” to make an informed decision about whether there is indeed a nearby structure.

For step (i), the estimation of a UPO period, we proceed as follows: A maximum cut-off growth/decay rate, $\lambda_{r}^{\text{max}}$, is chosen to select a subset of the DMD modes, $\tilde{\Lambda}:=\{\lambda\;|\;|Re(\lambda)|\leq\lambda_{r}^{\text{max}}\}$. The truncated set of eigenvalues $\tilde{\Lambda}$ contains at least one mode with $Im(\lambda)=0$ and pairs of complex conjugate eigenvalues. Of the modes with $Im(\lambda)=0$, the one with the smallest real component is assumed to be the time average of a nearby UPO, while the first $n$ modes with positive imaginary components are assumed to constitute the UPO’s fundamental frequency and its $n-1$ higher harmonics. Under this assumption, a fundamental frequency can be estimated via

[TABLE]

where $\omega_{j}:=Im(\lambda_{j})$. The guess for a velocity field at a point on the UPO is then built from the time average mode and the first $n$ complex conjugate pairs (this procedure is described in more detail in §2.3 below), while the period is esimated as $T_{g}=2\pi/\omega_{f}$. For step (ii), the estimate of how good such a guess is expected to be, we determine the degree to which the finite-frequency modes approximate a series of harmonics,

[TABLE]

with $\varepsilon_{\omega}=0$ for a perfectly periodic signal. In the following analysis, DMD results for which values of $\varepsilon_{\omega}$ fall below a prescribed threshold are taken to indicate the likely existence of a nearby UPO.

2.2 Motivating example: Gently periodic edge state

To demonstrate the mechanics of the approach outlined above we first apply it to a trajectory that shadows the simplest UPO in this box (the edge state) for several complete cycles before becoming turbulent. This ‘gentle’ periodic orbit was first discovered by Kawahara & Kida (2001) and has a period $T_{\omega}\approx 86.4$ in this box (see also Kawahara, 2005). The trajectory is obtained via repeated bisection between relaminarising and turbulent trajectories (Schneider et al., 2007) and we store one full period at a resolution of $\Delta t=0.1$.

DMD spectra from two calculations on this dataset are reported in figure 2. The two calculations are identical in all respects apart from the length of the observation window, $T_{w}$, over which the DMD is performed. In the first case, the length of the observation window matches the period of the underlying UPO, $T_{w}\approx T_{\omega}$, and the DMD spectrum consists of a single neutral mode and a series of near-perfect harmonics. The fundamental frequency estimated from the first purely imaginary mode and the first harmonic matches the fundamental frequency of the UPO, while the associated value $\varepsilon_{\omega}$ is very close to zero. This result is entirely expected owing to the observation of a full period and the proximity of the trajectory to the UPO.

More interestingly, the results of the DMD calculation on the window $T_{w}=T_{\omega}/4$ also clearly show the signature of the UPO despite the short observation time. The eigenvalues are now slightly off the neutral line $\lambda_{r}=0$, but there are still clearly repeated harmonics (the third harmonic $\sim 3\omega_{f}$ is missing). The period estimated from the first two harmonics is very close to the true period of the UPO, and the corresponding value of $\varepsilon_{\omega}\lesssim 10^{-3}$.

The two DMD calculations are compared further in figure 3. The neutral DMD mode obtained over the short time window $T_{\omega}/4$ is qualitatively similar to the true neutral Koopman mode associated with the orbit (the DMD mode from the observation window $T_{w}=T_{\omega}$), though the performance is notably worse for the higher harmonics. However, initial conditions built from these modes (via the procedure described in the following section) still rapidly converge to the UPO when used as an initial guess in Newton-Krylov routine.

The performance of shortened observation windows on the edge UPO is examined for a wide range of window lengths, each averaged over a large range of possible start times $t_{s}\in[0,T_{\omega})$, in figure 4. The results indicate that short windows $T_{w}\lesssim 0.4T_{\omega}$ can perform well sometimes. On average, the DMD gets progressively worse at identifying the underlying UPO as the length of the observation window is reduced, though the accuracy of an individual computation depends subtly on the start time of observations $t_{s}$ as this controls what segment of the UPO is seen by the DMD. Importantly, the good results shown for the $T_{w}=T_{\omega}/4$ window in figure 2 would not have been obtained if an alternate start time had been selected. General rules of thumb as to what exactly the DMD needs to see in order to provide a robust estimate of the period $T_{\omega}$ of a UPO are discussed in §4 following the presentation of results from a turbulent computation.

2.3 Building a guess for a periodic orbit

While the above results were obtained by examining trajectories that shadow a known orbit, the periodic orbits close to a turbulent trajectory will not be known beforehand, and we would like to also use the output of DMD to build initial ‘guesses’ to supply to a Newton solver if the value of $\varepsilon_{\omega}$ drops below the prescribed threshold. The approach we adopt is outlined here.

For DMD results where $\varepsilon_{\omega}<\varepsilon_{\omega}^{\text{thresh}}$ for a given number of harmonics $n$, the period of the UPO is estimated from the fundamental frequency $\omega_{f}(n)$, defined in equation (8), by $T_{g}=2\pi/\omega_{f}$. The neutral mode and the first $n$ complex conjugate pairs used to estimate $\omega_{f}(n)$ are then used to construct a guess for the velocity field on the UPO,

[TABLE]

There are a variety of ways to define the unknown coefficients $\{a_{j}\}$. We set them such that the average (squared) deviation of the guess from the turbulent trajectory used in the DMD is smallest, i.e. by minimising

[TABLE]

The minimising $\hat{\mathbf{a}}$ has been presented previously in Page & Kerswell (2018). Finally, he time $t$ in (10) is selected as $\hat{t}=\arg$ $\min_{t}|\mathbf{u}_{g}(t)-\mathbf{u}(t)|$.

3 Turbulent orbits

Motivated by the success of the DMD approach at identifying properties of the edge-state UPO on short time windows, we now apply our method to turbulent data and attempt to converge UPO guesses generated as described above with a Newton solver. To allow for a direct comparison to recurrent flow analysis, we perform many DMD computations on a short time window which is passed through a continuous turbulent trajectory of length $>1000$ advective time units.

3.1 Recurrent flow analysis

The results of a recurrent flow analysis are summarised in the top panel of figure 5, where we report contours of $R(t_{s},T)$ (see equation 1). Local minima which fall below a threshold value of $R_{\text{thresh}}=0.07$ are supplied as guesses in a Newton-Krylov algorithm. Of the 13 guesses highlighted, 6 converged (green $+$, while red $\times$ indicate failures).

It is notable that the guessed period, determined from the $T$ at which the local minimum in $R(t_{s},T)$ occurs, can often differ significantly from the period of the converged solution (note the vertical distance between green $+$ and green $o$). For example, the (good) guesses flagged at $t=235$ and $t=262$ converged to periodic orbits with relative errors in the period against the guess of $12.3$% and $10.5$% respectively. As we shall see below, the combination of the state $\mathbf{u}(t)$ and period $T$ which correspond to a local minima in $R(t,T)$ can often result in failure in the Newton algorithm even if there is a nearby UPO. This is particularly evident in figure 5 for $t>800$, where recurrent flow analysis fails to identify and successful guesses, while our new approach results in three converged UPOs.

3.2 Dynamic mode decomposition

The output of the DMD-based UPO detection algorithm outlined in §2.1 is shown directly below the recurrence plot in figure 5, with guessed periods $T_{g}$ reported based on both two and three harmonics. To generate these results, DMD computations were performed with a fixed observation window of length $T_{w}=60$, which was moved through the time series in steps of $\Delta_{w}=5$. Highlighted points along the curves $T_{g}(t_{s})$ indicate that the value of $\varepsilon_{\omega}$ (equation 9) dropped below a specified threshold. These points tend to be arranged in clusters with similar guesses for the period $T_{g}$, consistent with the trajectory briefly shadowing a UPO before being flung out of a particular Koopman expansion zone (Page & Kerswell, 2018). The detection threshold was set at $\varepsilon_{\omega}=10^{-3}$ for $n=2$ and $\varepsilon_{\omega}=5\times 10^{-3}$ for $n=3$; the threshold was increased for $n=3$ since resolution of the third harmonic is more challenging. The increase in the number of harmonics used in the search for UPOs results in a drop in the number of proposed guesses, most notably for guesses with long periods. For example, the number of guesses with $T_{g}>80$ falls from 16 to to 3. However, in general there is good correspondence for the locations $\{t_{s}^{i}\}$ of the possible UPOs between the $n=2$ and $n=3$ results.

The choice of observation window length $T_{w}=60$ was motivated by two factors: (i) if the time window is too long and the turbulence visits the neighbourhoods of multiple UPOs, the spectrum from the DMD will likely be unrelated to any of them (see the discussion in the introduction and Page & Kerswell, 2018) and (ii) $T_{w}=60$ is shorter than many of the known UPOs in this box, most of which have periods $T_{\omega}>75$ (Cvitanovic & Gibson, 2010) and will test our method’s ability to identify and converge solutions without a near recurrence. One immediate consequence of the choice $T_{w}=60$ is a loss of resolution of the longest UPO found by the recurrent flow analysis at $t_{s}=412$ with period $T_{\omega}=103.01$, though notably the DMD-based algorithm does flag a pair of likely UPOs with different periods close to this point – a behaviour which is discussed in more detail below. However, in general DMD identifies guesses with similar periods to the UPOs spotted in the recurrent flow analysis. Some heuristics around characteristics possessed by ‘good’ DMD guesses are discussed in §4, though ultimately we suspect the selection of $T_{w}$ in new applications is something that will have to be determined via trial-and-error subject to the tradeoff described above.

All of the guesses in figure 5 identified by DMD in the 2-harmonic case were supplied as initial conditions in a Newton-Krylov routine. Of the 46 guesses highlighted, 20 converged to periodic orbits. This total includes UPOs found in regions where the recurrent flow analysis indicated the presence of a structure, but also in regions where no near-recurrent episodes were detected; note the guesses for $t_{s}\lesssim 50$ and for $150\lesssim t_{s}\lesssim 200$. The various UPOs converged from the DMD-based guesses identified in figure 5 are summarised in Table 1, and we will now discuss differences between these results and the recurrent flow analysis in some detail. Our discussion will explore (i) cases where a UPO is found using recurrent flow analysis but the DMD-generated guess fails to converge; (ii) cases where the DMD predictions and recurrence predictions differ and different UPOs are converged; (iii) cases where DMD and a recurrence analysis show similar predictions but only the DMD guess converges and finally (iv) cases where DMD yields a converged UPO in the absence of a near recurrence.

(i) There are two instances where both the recurrent flow analysis and our DMD method indicate a likely UPO, but only the guess from the recurrent flow analysis converges in the Newton solver ($t_{s}\approx 262$ and $t_{s}\approx 748$ in figure 5). Intriguingly, the guesses for the period generated from the DMD are more accurate than those identified in the recurrence plot. This suggests that the DMD-mode based guesses constructed using equation (10) constitute poor initial conditions for the Newton solver. In instances like these with both a near recurrence and a DMD indicator, a more robust approach could potentially combine elements of both methods, perhaps with a guess that is a concatenation of the state itself and the period coming from DMD, i.e. $[\mathbf{u}(t_{i}),T_{g}]^{T}$.

(ii) On several occasions both the DMD and recurrent flow analysis identify possible UPOs at the same approximate $t_{s}$, but the converged solutions differ. In some of these cases, it can be argued that one approach is ‘more accurate’ than the other. For example, compare the UPO with period $T_{\omega}=103.01$ found from a near recurrence at $t_{s}=412$ to the DMD predictions at $t_{s}\in\{410,415\}$. As described above, the short DMD time window is unable to accurately parameterise the long UPO with significant errors in the predicted period(s). However, both DMD-based guesses do converge – to UPOs with different periods that are commensurate with the period obtained via recurrent flow analysis. The change in period in the Newton iterations is striking and is not something we have observed in guesses from recurrent flow analysis. There are other times where the UPOs found via DMD and recurrence differ (e.g. $t_{s}\approx 128$, $t_{s}\approx 730$) but the differences in the predicted and true periods are less extreme and it is difficult to argue which guess was “better”. Here, there is a plausible connection with the project-then-search method described by Ahmed & Sharma (2017) for finding equilibria. In that method, known exact solutions are projected onto a few resolvent modes before being input as new guesses to a Newton solver, resulting in a large number of new equilibria which are qualitatively similar to the known solutions. In our approach, we initialise the Newton solver with a guess constructed from just a few (5 for the 2 harmonic cases) DMD modes.

These observations raise an interesting ambiguity: given that there appear to be groups of UPOs which are visually very similar (i.e. constitute a closed loop of the same physical processes) with very similar periods, many of which can be accessed from a Newton search with slightly different initial conditions, can we determine which solution the turbulent flow is actually shadowing? Although we don’t explore this here, more insight could potentially be gained by employing a longer observation window (e.g. taken to be slightly longer than any identified near recurrence) and by using more harmonics to generate the initial guess.

There is less ambiguity in cases where either the recurrent flow analysis guess fails (iii), or where no recurrence is observed (iv). In these situations the DMD-based approach offers significant advantages. An example of the former scenario is seen in the DMD-based successes for $t_{s}\in\{815,825\}$, where two near recurrences are found nearby but neither converges. In both of these cases, the period generated by the DMD-based approach differs non-negligibly from the time $T$ at which $R(t_{s},T)$ has a local minima. This is another scenario where the two methods can work well together to converge exact solutions where they may previously have been dismissed.

Examples of the second scenario (iv), a converged UPO from DMD without a corresponding near recurrence, are observed at $t_{s}\in\{180,200,780\}$. These results are also notable due to the quality of the guess (e.g. the relative error in the period at $t_{s}=780$ is $3.3$%). The guess at $t_{s}=780$ is particularly interesting since the DMD window corresponds to only $2/3$ of the full period, and in figures 6 and 7 we examine it further, comparing the DMD modes used to build the UPO guess to the true Koopman modes of the UPO (the coefficients of a Fourier series). There is good qualitative agreement between the modes, providing further evidence that DMD of modest-duration turbulent trajectories can really see the essence of a nearby exact coherent structure.

It is worth emphasising that this final class of guesses could not be obtained another way as there is no near recurrence in the data, and so this method appears to offer some promise at higher Reynolds numbers where the UPOs tend to be more unstable. A common feature of the various guesses from DMD without a near recurrence is that only the first two harmonics are resolved accurately (see the difference in the bottom two panels of 5). This behaviour was also observed for the simple edge state examined in §2.2, where we observed improved resolution of the higher-frequency Koopman modes as the length of the observation window was increased. However, for turbulent trajectories the observation window cannot be increased arbitrarily and would require a careful adjustment to avoid the inclusion of snapshots where the flow has left the neighbourhood around the UPO where the local Koopman decomposition holds.

4 Discussion

In this section we briefly examine some features of our initial guesses to provide some insight into what constitutes a good guess for UPO. We do this by exploring how well the short-time DMD computations can see Koopman eigenfunctions related to the structure of interest and by performing DMD on the converged turbulent UPOs themselves to identify a rough rule-of-thumb on where the best DMD time windows tend to be located in state space.

The idea that DMD can identify the signature of UPOs in a long trajectory is based on the fact that these exact coherent structures possess a local Koopman decomposition which is related to both the period of the solution itself and the eigenvalues governing the dynamics within its linear subspace (Mezic, 2017). For short time windows, DMD can typically only resolve a few eigenvalues related to the structure itself (the UPO’s Fourier harmonics) as seen in §3 (see also Page & Kerswell, 2018). The ‘quality’ of a particular DMD-based guess, and some indication as to whether the convergence of a UPO occurred due to luck, can be examined by computing approximations to the Koopman eigenfunctions (see equation 7) using left eigenvectors from the short-time-window DMD,

[TABLE]

and comparing them to those obtained from the converged UPO itself,

[TABLE]

The above procedure is performed for two periodic orbits and their original guesses in figure 8, where the velocity field used in the projections is that of the UPO itself. One of the guesses examined was selected due to its particularly poor esimate of the period of the UPO (at $t_{s}=10$, $T_{g}=62.12$ and $T_{\omega}=85.60$, a relative error of $27.5$%) and the resolution of the Koopman eigenfunctions is also poor. In particular, the projection shows sharp gradients and appears almost discontinuous in places. It does not at all resemble the true Koopman eigenfunction extracted from the UPO itself. On the other hand, the second guess was considerably more accurate and the approximation to the Koopman eigenfunctions is qualitatively much better. For other UPOs we have observed similar behaviour (not shown); that more accurate initial guesses (as measured by the predicted period) tend to produce Koopman eigenfunctions which resemble qualitatively those of the UPO itself.

Computation of the Koopman eigenfunctions along the UPO itself is only possible after the solution has been converged, and it is natural to speculate whether there are additional indicators of guess quality that can be assessed before the Newton solver is initialised. While the eigenspectrum-based measure $\varepsilon_{\omega}$ (equation 9) provides a relatively robust classifier to identify sensible initial guesses, we have observed some additional heuristics that can be used to screen a large number of initial guesses. When examinining a long time series with short DMD windows, perhaps the most obvious clue is the observation of ‘clusters’ of similar guesses (in terms of period $T_{g}$) in time. This is evident in some of the successful guesses highlighted in figure 5 and could be further emphasised by shrinking the timestep between subsequent DMD calculations (in figure 5 this is $\Delta_{w}=5$).

Another interesting feature of good guesses from short DMD windows is that they tend to share certain characteristics in terms of their location in an energy production-dissipation plane. To demonstrate this effect, we perform many DMD calculations with fixed window length $T_{w}=60$ on two of the longer converged UPOs and plot the ‘best’ and ‘worst’ time windows in figure 9, where the criteria used to rank the quality of the DMD spectra was the error in the guessed period (using the first two harmonics), $|2\pi/\omega_{f}(n=2)-T_{\omega}|/T_{\omega}$. The best guesses tend to sample the upper-right quadrant (the fast, high production, high dissipation region) at the expense of the slower, gentler region in the lower-left corner. Note that for longer orbits, good estimates (relative error $\lesssim 5$%) can also be obtained on trajectories which mostly shadow the portion of the curve $D<I$ (not shown) – the key requirement appears to be the inclusion of the most extreme values of production and dissipation.

We note that the presentation in this manuscript has focused on one particular time window length and DMD design in order to introduce the method in a clear and systematic way. There is clearly scope for the results to be improved by adjusting the DMD design based on local performance. For example, where near recurrences are observed it would make sense to increase the time window to the predicted period, while for cases where no near recurrence is observed but DMD does predict a UPO, it would be interesting to search for a ‘best’ local guess by small adjustments in window length and position.

5 Conclusion

We have presented a new method based on DMD to both identify the signature of nearby UPOs in time-series of turbulent flows, and to generate robust initial guesses that can serve as inputs to a Newton-Krylov algorithm. The approach is designed to be applied to short turbulent trajectories, with the existence of a nearby UPO deemed ‘likely’ when the DMD eigenvalue spectrum exhibits repeated harmonics of a fundamental frequency. Initial conditions for the Newton solver are then built from these near-periodic DMD modes and a guess for the period based on the estimated fundamental frequency.

The ability of the method to function without the need for a near recurrence within the observation window was demonstrated on a simple UPO before we applied it to a long turbulent trajectory. With a modest time window (shorter than most known UPOs), the new method identified many more UPOs than were found with a recurrent flow analysis of the same dataset. These solutions were found in regions where (i) a near recurrence was found and converged, (ii) a near recurrence was identified but the Newton solver failed to converge and (iii) where no near recurrence was flagged. These outcomes indicate that the method can be both a useful complementary tool and a potentially powerful alternative to recurrent flow analysis. In particular, the success of the method in regions where no near recurrence occurs suggests it may be of use at higher Reynolds numbers where the probability of shadowing a UPO for a full cycle drops.

An interesting behaviour was found at times where both the DMD based method and recurrent flow analysis flagged the probable existence of a UPO. In some cases, the DMD approach converged to a different solution than the guess from a recurrent flow analysis. Often, the error in the predicted period for these alternate solutions was comparable to the recurrent flow analysis, and given the sensitivity of the Newton method to initial conditions it is difficult to assert which one was indeed the solution being shadowed in phase space. Motivated by the variable discrepancy between the predicted period and that of the converged solution, we also presented some analysis of the features of both ‘good’ and ‘bad’ initial guesses generated by the method. Good guesses tend to require the DMD to see the most extreme dissipation events along an orbit at the expense of the slower, more gentle dynamics.

There are a number of interesting avenues open to further investigation. The most natural question is whether the method can be extended to much higher Reynolds numbers to converge UPOs, where recurrent flow analysis struggles to identify any guesses at all (e.g. see Chandler & Kerswell, 2013; Lucas & Kerswell, 2015). There are also some intriguing properties of the DMD-based approach that merit further study. For example, DMD on short time windows performs “best” (i.e. the predicted period matches that of the underlying simple invariant set) when it sees the faster, high-dissipation region of the UPO. In instances where the observation window does not have this property, can the performance of the DMD be improved, for instance by modifying how snapshot pairs are distributed in time?

Acknowledgements

We are very grateful to Yves Dubief for sharing his viscoelastic code with us.