#  Data-driven Koopman operator predictions of turbulent dynamics in models of
shear flows

C. Ricardo Constante-Amores1,2 Andrew J. Fox1 Carlos E. Pérez De Jesús1 and
Michael D. Graham1<EMAIL_ADDRESS>1Department of Chemical and Biological
Engineering, University of Wisconsin-Madison, Madison WI 53706, USA
2Department of Mechanical Science and Engineering, University of Illinois,
Urbana Champaign, IL 61801, USA

###### Abstract

The Koopman operator enables the analysis of nonlinear dynamical systems
through a linear perspective by describing time evolution in the infinite-
dimensional space of observables. Here this formalism is applied to shear
flows, specifically the minimal flow unit plane Couette flow and a nine-
dimensional model of turbulent shear flow between infinite parallel free-slip
walls under a sinusoidal body force (Moehlis et al., New J. Phys. 6, 56
(2004)). We identify a finite set of observables and approximate the Koopman
operator using neural networks, following the method developed by Constante-
Amores et al., Chaos, 34(4), 043119 (2024). Then the time evolution is
determined with a method here denoted as “Projected Koopman Dynamics”. Using a
high-dimensional approximate Koopman operator directly for the plane Couette
system is computationally infeasible due to the high state space dimension of
direct numerical simulations (DNS). However, the long-term dynamics of the
dissipative Navier-Stokes equations are expected to live on a lower-
dimensional invariant manifold. Motivated by this fact, we perform dimension
reduction on the DNS data with a variant of the IRMAE-WD (implicit rank-
minimizing autoencoder-weight decay) autoencoder (Zeng et al., ML&ST 5(2),
025053 (2024)) before applying the Koopman approach. We compare the
performance of data-driven state space (neural ordinary differential equation)
and Projected Koopman Dynamics approaches in manifold coordinates for plane
Couette flow. Projected Koopman Dynamics effectively captures the chaotic
behavior of Couette flow and the nine-dimensional model, excelling in both
short-term tracking and long-term statistical analysis, outperforming other
frameworks applied to these shear flows.

###### keywords:

shear flows, Koopman, data-driven models, chaos

##  1 Introduction

The study of turbulence, particularly in wall-bounded flows, holds paramount
importance in understanding fluid dynamics as it has profound implications for
various engineering applications (Pope, 2000; Avila et al., 2023).
Approximately $25\%$ of the energy consumed by industry is dedicated to
transporting fluids through pipes and channels, and about one-quarter of this
energy is dissipated due to turbulence occurring near walls  (Jiménez & Moin,
1991; Waleffe, 2001; Jiménez, 2013). To understand the underlying physical
mechanisms of turbulence, researchers have focused on canonical flows such as
plane Couette flow, which serves as an ideal test case due to its simplicity
and well-defined boundary conditions, as evidenced by numerous studies, such
as Nagata (1990); Jiménez & Moin (1991); Waleffe (1998); Kawahara & Kida
(2001); Waleffe (2001); Gibson (2002); Itano & Generalis (2009).

The focus of the present work is the extent to which modern data-driven
methods can capture the dynamics of turbulent shear flows near transition.
Because of its geometric and dynamical simplicity, the ‘minimal flow unit’
(MFU) of plane Couette flow (PCF) has been the object of much study Jiménez &
Moin (1991); Hamilton et al. (1995); Gibson (2002); Gibson et al. (2008,
2009); Halcrow et al. (2009). The MFU is the minimum domain size in which
turbulence will self-sustain; thus contains the key components of turbulent
dynamics, i.e. the ‘self-sustaining process’ (SSP) (Hamilton et al., 1995) in
which low-speed streaks near the wall become wavy, triggering their breakdown
to generate rolls, which in turn rolls lift fluid from the walls, regenerating
the streaks. Given that the state space dimension (number of mesh points times
number of velocity components) for fully resolved direct numerical simulation
(DNS) even of an MFU may be $\mathcal{O}(10^{5})$ (Jiménez & Moin, 1991;
Hamilton et al., 1995), significant efforts have been made to develop low-
dimensional models capable of capturing the SSP (Waleffe, 1997; Moehlis et
al., 2004; Gibson, 2002; Kaszás et al., 2022). Waleffe (1997) proposed a
8-mode model for turbulent shear flow between infinite parallel free-slip
walls subjected to a sinusoidal body force. The modes were selected based on
intuition from the self-sustaining process (SSP), and a Galerkin projection
was then performed onto these modes. Moehlis et al. (2004) proposed the
Moehlist, Faisst, and Eckhardt (MFE) model, which contains an additional mode
that enables modification of the mean profile by the turbulence; as the
fluctuations change the mean profile for turbulent shear flows. The MFE model
is attractive because it captures key features of shear flow transition near
onset, including interaction between roll and streak modes, intermittency and
ultimate relaminarization while maintaining computational tractability. The
present work addresses data-driven modeling of the plane Couette MFU both in
the form of full DNS results and the MFE model.

In data-driven modelling of the dynamics, two different approaches can be
taken: the state space approach and the function space or observable space
approach. For the latter, the notion that any finite-dimensional non-linear
dynamical system can be transformed into an equivalent linear system by
lifting it to an infinite-dimensional space of observables was introduced
almost a century ago by Koopman (1931). However, its significance emerged over
a decade ago with its connection to dynamic mode decomposition (DMD), as
outlined by Schmid (2010), and connected to the Koopman operator by Rowley et
al. (2009). The non-linear dynamical system is transformed into a linear one
using a set of observables, with the infinite-dimensional Koopman operator
applying to all observables. Then, the entire dynamics is represented by a
linear expansion of the modes (eigenfunctions) and their corresponding
eigenvalues (Schmid, 2022; Mezić, 2013). The challenge lies in lifting the
state space to the observable space. In practice, the infinite-dimensional
Koopman operator is approximated by a matrix-valued operator; this
approximation aims to maintain the essential properties of the operator while
operating only on a discrete finite set of observables. Several strategies
have been proposed to ensure that the input function space, and subsequently
the selection of observables, is sufficiently ‘rich’. For a given method,
there is one dictionary of observables, with many elements. Williams et al.
(2014) proposed the extended dynamic mode decomposition (EDMD), which uses
Legendre polynomials of the state space coordinates as a dictionary to create
the observable space. Their approach primarily focused on low-dimensional
systems, highlighting the difficulty of choosing appropriate dictionaries
(Korda & Mezić, 2018). Even for a simple dictionary such as orthogonal
polynomials, the number of dictionary elements grows rapidly with the number
of dimensions.

Recently, attempts have been made to use neural networks to construct a good
dictionary of observables, resulting in a good approximation of the matrix-
valued Koopman operator – see Li et al. (2017); Yeung et al. (2019); Lusch et
al. (2018); Otto & Rowley (2019); Alford-Lago et al. (2022); Constante-Amores
et al. (2024). Li et al. (2017) presented extended dynamic mode decomposition
with dictionary learning (EDMD-DL) which uses neural networks to adapt the
dictionary of observables to the data.  Recently, Constante-Amores et al.
(2024) proposed a variant of the EDMD-DL that leverages automatic
differentiation to simultaneously determine the dictionary of observables and
the corresponding approximation of the Koopman operator. In their study, the
authors contrast two different ways to handle temporal evolution, referred to
as KDLoo (Koopman Dictionary Learning, Observable-Only evolution) and KDLso
(State-Observable alternation). The former approach maintains the linearity of
the Koopman formalism but suffers from poor performance in time evolution. The
latter approach behaves substantially better in time evolution but is
nonlinear. We use the latter approach here, denoting it as ‘Projected Koopman
Dynamics’. Finally, Pan & Duraisamy (2024) emphasized the importance of
considering symmetries when learning the Koopman operator for low-dimensional
systems.

The accurate simulation of MFU PCF requires a large state space to resolve all
the relevant spatial and temporal scales. For example, Linot & Graham (2023);
Hamilton et al. (1995) required approximately $\mathcal{O}(10^{5})$ degrees of
freedom to capture the complex nonlinear turbulent dynamics. Performing
Koopman modeling in this state space poses computational challenges, as it
involves lifting this high-dimensional state space into even a higher-
dimensional space of observables.  However, the Navier-Stokes equations (NSE)
are dissipative PDEs; as a result, the long-term dynamics of the system are
confined within a finite-dimensional invariant manifold with fewer degrees of
freedom than the full state-space dimension(Temam, 1989; Zelik, 2022). Thus if
a coordinate representation (i.e. a dimension-reduction process) can be found
for data on an invariant manifold, then one can  lift the coordinates of this
invariant manifold into its equivalent observable-space representation, as was
shown by Constante-Amores & Graham (2024) for the Koopman modelling of
coherent state dynamics on invariant manifolds. Since in the present work, we
are considering time evolution with Koopman operators on invariant manifolds,
it is worth noting that Nakao & Mezic (2020) and Mezic (2020) showed that
these correspond to joint zero level sets of Koopman eigenfunctions.

The dimension of an invariant manifold for long-time dynamics is not known a
priori, so must be estimated from data. Various methods for estimating
manifold dimensions have been introduced in the literature. For the Kuramoto-
Sivashinsky equation (KSE), the manifold dimension has been estimated using
covariant Lyapunov vectors and monitoring the drop of the Lyapunov spectrum
(Yang & Radons, 2012; Yang et al., 2009). These methods require access to the
governing equations and high-precision solutions, which are generally not
available for the NSE. When the governing equations of motion are unknown and
only data is available, the most common method is Principal Component Analysis
(PCA, also known as Proper Orthogonal Decomposition, POD, in the fluid
dynamics community) which is based on the idea of the projection of the state
onto the set of orthogonal modes that capture the maximum energy of the data
(Jolliffe, 1986; Abdi & Williams, 2010; Holmes et al., 2012). PCA, being a
linear reduction method, projects data onto a flat manifold, which is not
generally appropriate for complex nonlinear problems, and subsequently, it
struggles to estimate their true dimension. Thus, nonlinear methods can be
advantageous. One of the prevalent methods for nonlinear dimension reduction
is the use of an undercomple autoencoder, a pair of neural networks in which
one network maps from a high-dimensional space to a low-dimensional space, and
the other maps back (Kramer, 1991; Hinton & Salakhutdinov, 2006; Milano &
Koumoutsakos, 2002).  For very high-dimensional systems, it can be
advantageous to perform an initial linear dimension reduction step with PCA,
then use an autoencoder for further dimension reduction (Linot & Graham, 2020;
Young et al., 2023) It can also be advantageous to apply PCA and an
autoencoder in parallel to learn a low-dimensional nonlinear representation as
a correction to a linear representation Linot & Graham (2020). Recently, Jing
et al. (2020) introduced a variation of a standard autoencoder to learn low-
rank representations for image-based classification. This method was extended
to data from dynamical systems by Zeng et al. (2024), demonstrating that this
architecture can yield robust and precise estimates of manifold dimension, as
well as an orthogonal manifold coordinate system. The present work uses a
variation of this approach.

Once, we have access to the manifold coordinates, i.e. a low-dimensional
representation of the data, we can learn an evolution equation in these
coordinates. Some examples where a state space modelling in manifold
coordinates for chaotic systems include the KSE equation (Otto & Rowley, 2019;
Linot & Graham, 2020; Liu et al., 2024), Kolmogorov flow Pérez-De-Jesús &
Graham (2023) and Couette flow Kaszás et al. (2022); Linot & Graham (2022).
Linot & Graham (2020) introduced an invariant-manifold-based data-driven
approach known as ‘DManD’ (data-driven manifold dynamics). DManD is based on
the idea of performing a nonlinear reduction of the full space with
autoencoders to find a mapping to the manifold coordinates, and then use a
neural ODEs (NODE) to learn an evolution equation in manifold coordinates – in
NODE, a neural network learned from data parameterises the vector field
(right-hand side of the ODE) determining the dynamics (Linot & Graham, 2022;
Chen et al., 2019). This framework applied to MFU PCF only requires 18 degrees
of freedom to build reduced-order models that faithfully capture the turbulent
dynamics over both short and long timescales, including the streak breakdown
and regeneration process and Reynolds stresses (Linot & Graham, 2023). We
remark that NODE is an advantageous approach as the dynamics are Markovian
(consistent with dynamics that arise from ODEs) and continuous in time.

We conclude our introduction by discussing several works in data-driven
modeling of the MFE system. Racca & Magri (2022) modeled MFE using Echo State
Networks, which expand the nonlinear input state space to create a high-
dimensional system (up to two orders of magnitude higher). Their predictions
account for history in their time evolution, making the model non-Markovian,
with temporal evolution based on discrete time maps. Eivazi et al. (2021)
employed a Koopman-based framework, proposed by Khodkar et al. (2019),
incorporating a nonlinear forcing term. They construct the forcing term as a
library of candidate nonlinear functions, allowing the model to identify the
most significant nonlinearities from the data. Here, the inclusion of
predetermined forcing candidates is needed beforehand. Their method can
capture accurate short-term predictions and reproduce long-term statistics,
but the temporal evolution includes delay embedding, thus rendering the model
non-Markovian. Similarly, approaches based on transformer architectures have
been successful for prediction of flow problems, but are also non-Markovian
(Solera-Rico et al., 2024). Finally, Fox et al. (2023) utilise a state space
approach denoted as Charts and Atlases for Nonlinear Data-Driven Dynamics on
Manifolds, or CANDyMan (Floryan & Graham, 2022). In their approach, the time
series of MFE data is divided into overlapping clusters or charts, and an
evolution equation is learned for each chart to construct a final global
model. This clustering of the data significantly improved the performance of
the models.

|
---|---
(a) | (b)
Figure 1: Schematic representation of the shear flows considered in this
study. (a) Evolution of three amplitudes, $a_{1},a_{6},a_{9}$ and
corresponding kinetic energy from one time series of the MFE data set. (b)
Three dimensional snapshot of turbulent MFU Plane Couette Flow at $Re=400$,
coloured by the streamwise velocity.

In this study, we aim to build data-driven Koopman-operator-based predictions
of turbulent dynamics in models of shear flows using the approach presented in
Constante-Amores et al. (2024). To achieve this, we begin with the nine-
dimensional MFE model, and then we extend our framework to the high-
dimensional MFU PCF (see Figure 1, for a schematic representation of the MFE
model and MFU PCF) . We remark that the MFE model, being inherently low-
dimensional, does not require dimension reduction. In contrast, the plane
Couette flow requires dimension reduction to manifold coordinates prior to
lifting the system to the observable space. To find these manifold
coordinates, we introduce a variant of IRMAE-WD. Then, we employ a Koopman
approach to learn an effective observable-space representation and dynamical
system using the framework that we termed as ‘Koopman-Data-Driven Manifold
Dynamics’ (KDManD). We also make comparisons with DManD (state space approach
using NODEs). The paper is organized as follows, Section 2 provides a detailed
explanation of the framework for finding good observables for an approximate
Koopman representation with neural networks as well as the time-evolution of
the observables, along with a description of the approach to find the manifold
coordinates for PCF. Section 3 presents the description of the MFE model and
PCF system, the results concerning the reconstruction of short-time and long-
time statistics. Finally, concluding remarks are summarized in Section 4.

##  2 Framework

Figure 2 shows a diagram of the framework applied to a full DNS plane Couette
flow, where the evolution of states and observables are done in manifold
coordinates, which encapsulate the long-time dynamics as the full DNS.

---
Figure 2: A cartoon of the data-driven framework applied to the full DNS of
turbulent plane Couette flow. The top panel represents the methodology to find
a coordinate transformation from the full state space
($\mathbf{x}\in\mathbb{R}^{d_{N}}$) to a low-dimensional representation
($\boldsymbol{h}\in\mathbb{R}^{d_{h}}$) via $\boldsymbol{\chi}$ and
$\tilde{\boldsymbol{\chi}}$ In manifold coordinates, we can use a state space
modeling approach to update the state $\boldsymbol{h}\in\mathcal{M}$ using the
state space time-evolution operator $\boldsymbol{G}_{\delta t}$, i.e.,
$\boldsymbol{h}(t+\delta t)=\boldsymbol{G}_{\delta t}(\boldsymbol{h}(t))$.
Alternatively, the manifold representation can be lifted to the high-
dimensional space of observables where the dynamics behave linearly. The
bottom path represents the observable space in which the updates the
observables, $\mathcal{G}(\boldsymbol{h})\in\mathcal{F}$, is done using the
Koopman operator, $\mathcal{K}$, i.e., $\mathcal{G}(\boldsymbol{h}(t+\delta
t))=\mathcal{K}\mathcal{G}(\boldsymbol{h}(t))$. As a result, either the state
space or function space modeling leads to the same result.

###  2.1 Finding manifold coordinates

In this work, we consider systems that are characterized by deterministic,
Markovian dynamics, so if $\boldsymbol{x}\in\mathbb{R}^{d_{N}}$ represents the
full space state, then the dynamics can be represented by an ODE as

$\dfrac{d\boldsymbol{x}}{dt}=\boldsymbol{f}(\boldsymbol{x}).$ (1)

Here, $\boldsymbol{x}$ represents the full state space from highly resolved
direct numerical simulation (DNS). The long-time dynamics of dissipative
systems collapse onto an invariant manifold of embedding dimension $d_{h}$
(i.e., $d_{h}\ll d_{N})$.  When $\boldsymbol{x}$ is mapped into the manifold
coordinates, represented by $\boldsymbol{h}\in\mathbb{R}^{d_{h}}$, an
equivalent evolution equation can be expressed as

$\frac{d\boldsymbol{h}}{dt}=\boldsymbol{g}(\boldsymbol{h}).$ (2)

With the full state, we use an undercomplete autoencoder architecture to find
a mapping for the low-dimensional representation as
$\boldsymbol{h}=\boldsymbol{\chi}(\boldsymbol{x})$, along with the inverse
$\check{\boldsymbol{x}}=\check{\boldsymbol{\chi}}(\boldsymbol{h})$ so that the
full state space may be recovered (i.e.,
$\boldsymbol{x}\approx\tilde{\boldsymbol{x}}$).

To find $\boldsymbol{\chi}$ and $\check{\boldsymbol{\chi}}$, we use a
variation of IRMAE-WD (Zeng et al., 2024). In the original IRMAE-WD, the
autoencoder is formed by standard nonlinear encoder, $\mathcal{E}$, and
decoder, $\mathcal{D}$, networks with $n$ additional linear layers with weight
matrices $\mathcal{W}_{1},\mathcal{W}_{2}\cdots\mathcal{W}_{n}$ (of size
$d_{z}\times d_{z}$). The encoder finds a compact representation
$\boldsymbol{z}\in\mathbb{R}^{d_{z}}$, and the decoder performs the inverse
operation. In IRMAE-WD, a predetermined number of linear layers is specified
in advance. The variant that we introduce in this research allows the forward
pass to propagate through different ‘branches’ with varying numbers of linear
layers (e.g., an upper and lower range of linear layers are specified). We
refer to this architecture as IRMAE-WD-B. This variant improves the robustness
of dimension estimates for more complex systems as presented by Pérez-De-Jesús
(2024).

A schematic representation of the IRMAE-WD-B framework is shown in Figure 3. A
set of linear networks
$\mathcal{W}_{i}\left(\cdot;\theta_{\mathcal{W}}\right)$ is included between
the encoder and the various decoders $\mathcal{D}_{i}$ (i.e.,
$\boldsymbol{z}=\mathcal{W}_{i}\mathcal{W}_{i-1}\ldots\mathcal{W}_{1}(\mathcal{E}(\boldsymbol{x};\theta_{\mathcal{E}});\theta_{\mathcal{W}})$).
These decoders map the low-dimensional representation back to the full space
(e.g.,
$\tilde{\boldsymbol{x}}_{i}=\mathcal{D}_{i}\left(\mathcal{W}_{i}\mathcal{W}_{i-1}\ldots\mathcal{W}_{1}(\mathcal{E}(\boldsymbol{x};\theta_{\mathcal{E}});\theta_{\mathcal{W}});\theta_{\mathcal{D}_{i}}\right)$).
Therefore, the loss function for this architecture is represented by

$\mathcal{L}_{i}\left(\tilde{\boldsymbol{x}};\theta_{\mathcal{E}},\theta_{\mathcal{W}},\theta_{\mathcal{D}_{i}}\right)=\left\langle\left\|\boldsymbol{x}-\mathcal{D}_{i}\left(\mathcal{W}_{i}\mathcal{W}_{i-1}\ldots\mathcal{W}_{1}\left(\mathcal{E}\left(\boldsymbol{x};\theta_{\mathcal{E}}\right);\theta_{\mathcal{W}}\right);\theta_{\mathcal{D}_{i}}\right)\right\|_{2}^{2}\right\rangle+\frac{\lambda}{2}\|\theta\|_{2}^{2}.$
(3)

where $\langle\cdot\rangle$ is the average over a training batch,
$\theta_{\mathcal{E}}$ the weights of the encoder, $\theta_{\mathcal{D}_{i}}$
the weights of the decoder $i$, and $\theta_{\mathcal{W}}$ the weights of the
linear network which contains the linear layers used in the branch. We also
include a $L_{2}$ (“weight decay”) regularization to the weights with
prefactor $\lambda$. In this work $\lambda$ is fixed to $10^{-6}$, following
the work from (Zeng et al., 2024).

We note that this architecture will produce a number of loss values
corresponding to the number of linear layers. To train the network, we select
a branch and backpropagate through it to update the weights. The branch
selection process involves randomly selecting a path for each epoch of
training, resulting in an improvement in the autoencoder performance.

To find the low-dimensional representation, we optimize the loss function
described in equation 3 using an Adam optimiser in pythorch (Paszke, 2019). We
allow the model to choose among 10 different paths. Each path is characterized
by an increasing number of linear layers of size $d_{z}\times d_{z}$ (e.g.,
path 1 has one linear layer, while path 10 has ten linear layers). We ran for
1000 epochs, and to further improve the model at the end of the training, we
used the path with the lowest mean squared error (MSE) during the last $10\%$
of the total number of epochs.

Post-training, a singular value decomposition (SVD) is applied to the
covariance matrix of the latent data matrix z yielding matrices of singular
vectors $\mathsfbi{U}$, and singular values $\mathsfbi{S}$. Then, we can
project z onto $\mathsfbi{U}^{T}$ to obtain
$\mathsfbi{U}^{T}\textbf{z}=\boldsymbol{h}^{+}\in\mathbb{R}^{d_{z}}$ in which
each coordinate of $\textbf{h}^{+}$ is orthogonal and ordered by contribution
(here, $\boldsymbol{h}^{+}$ refers to the projection of latent variables onto
manifold coordinates). This framework reveals the manifold dimension $d_{h}$
as the number of significant singular values, indicating that a coordinate
representation exists in which the data spans $d_{h}$ directions. Thus, the
encoded data avoids spanning directions associated with nearly zero singular
values (i.e., $\mathsfbi{U}\mathsfbi{U}^{T}{\bf
z}\approx\mathsfbi{\hat{U}}\mathsfbi{\hat{U}}^{T}{\bf z}$, where
$\mathsfbi{\hat{U}}$ are the singular vectors truncated corresponding to
singular values that are not nearly zero). Leveraging this insight, we extract
a minimal, orthogonal coordinate system by projecting ${\bf z}$ onto
$\mathsfbi{\hat{U}}$, resulting in a minimal representation
$\mathsfbi{\hat{U}}^{T}{\bf z}={\boldsymbol{h}}\in\mathbb{R}^{d_{h}}$.

Figure 3: Implicit Rank Minimizing autoencoder with weight decay and branching (IRMAE-WD-B) framework. (a) Network architecture with regularization mechanisms. Note the different paths between the encoder, $\mathcal{E}$, and decoders, $\mathcal{D}_{i}$. (b) Singular value decomposition of the covariance of the learned latent data representation $\textbf{z}_{i}$. Table 1:  Here we list the architecture and parameters utilized in the studies of this paper. Between the $\mathcal{E}$ and $\mathcal{D}$, there are $n$ sequential linear layers of shape $d_{z}\times d_{z}$. The superscripts in $\boldsymbol{\Psi}$ and $\boldsymbol{g}$ denote the number of sequential layers with the specified number of node units. Case | Function | Shape | Activation | Learning Rate
---|---|---|---|---
MFE | $\boldsymbol{\Psi}$ | $9/(200)^{4}/\Psi+9\quad$ | $(\text{elu})^{4}/\text{elu}$ | $[10^{-4},10^{-5},10^{-6}]$
Plane Couette Flow | $\boldsymbol{\chi}$ | 502/1004/256/$d_{z}$ | ReLU/ReLU/lin | $[10^{-3},10^{-4}]$
| $\check{\boldsymbol{\chi}}$ |  $d_{z}$/256/1004/502 | ReLU/ReLU/lin | $[10^{-3},10^{-4}]$
| $\boldsymbol{\Psi}$ | $d_{h}/(200)^{6}/\Psi+d_{h}\quad$ | $(\text{elu})^{6}/\text{elu}$ | $[10^{-4},10^{-5},10^{-6}]$
| $\boldsymbol{g}$ | $d_{h}/(500)^{2}/d_{h}\quad$ | $(\text{ReLU})^{2}/\text{lin}$ | $[10^{-3},10^{-4}]$

###  2.2 function space modelling: approximation of the Koopman operator

We apply a function space approach to learn an evolution equation either in
the low-dimensional representation for plane Couette Flow or full space for
the MFE model. We note that the nomenclature used in this section is based on
using the low-dimensional representation. The function space approach is based
on the infinite-dimensional linear Koopman operator $\mathcal{K}$, which
describes the evolution of an arbitrary observable
$\mathcal{G}(\boldsymbol{h})$ from time $t$ to time $t+\delta t$:
$\mathcal{G}(\boldsymbol{h}(t+\delta
t))=\mathcal{K}\mathcal{G}(\boldsymbol{h}(t))$ (Koopman, 1931; Lasota &
Mackey, 1994; Mezic, 2023). The tradeoff for gaining linearity is that
$\mathcal{K}_{\delta t}$ is also infinite-dimensional, requiring for
implementation of some finite-dimensional truncation of the space of
observables.

Here we use a variant of the EDMD-DL approach  to determine both the
dictionary of observables and the corresponding approximation of the Koopman
operator,  described in detail in Constante-Amores et al. (2024). Given a
vector of observables
$\boldsymbol{\Psi}(\boldsymbol{h}(t)){{\color[rgb]{0,0,0}\in\mathbb{R}^{d_{K}}}}$
(each element of the observable $\boldsymbol{\Psi}$ is a neural network except
for the first $d_{h}$ elements, which correspond to the manifold coordinates
themselves), now there is a matrix-valued approximate Koopman operator
$\mathsfbi{K}$ such that the evolution of observables is approximated by

$\boldsymbol{\Psi}(\boldsymbol{h}(t+\delta
t))=\mathsfbi{K}\boldsymbol{\Psi}(\boldsymbol{h}(t)).$ (4)

Given a time series of observables assembled into a matrix whose columns are
the vector of observables at different times,

$\boldsymbol{\psi}(t)=\left[\boldsymbol{\Psi}(\boldsymbol{h}(t_{1})),~{}\boldsymbol{\Psi}(\boldsymbol{h}(t_{2}))\ldots\right]$
(5)

and its corresponding matrix at $t+\delta t$,

$\boldsymbol{\psi}(t+\delta
t)=\left[\boldsymbol{\Psi}(\boldsymbol{h}(t_{1}+\delta
t)),~{}\boldsymbol{\Psi}(\boldsymbol{h}(t_{2}+\delta t))\ldots\right],$ (6)

the approximate matrix-valued Koopman operator is defined as the least-squares
solution

$\mathsfbi{K}=\boldsymbol{\psi}({t+\delta t})\boldsymbol{\psi}(t)^{\dagger},$
(7)

where $\dagger$ superscript denotes the Moore-Penrose pseudoinverse. We aim to
minimise the difference between the observable prediction
$\tilde{\boldsymbol{\psi}}(\boldsymbol{h},t+\delta t;\theta)$ (e.g.,
$\tilde{\boldsymbol{\psi}}(\boldsymbol{h},t+\delta
t;\theta)=\mathsfbi{K}\boldsymbol{\psi}(\boldsymbol{h};\theta)$), and the
corresponding observable matrix at $t+\delta t$ (equation 6). Then, the loss
function is defined by

$\mathcal{L}_{K}(\boldsymbol{h},\theta_{K})=\left\|\boldsymbol{\psi}(\boldsymbol{h},t+\delta
t;\theta_{K})-\tilde{\boldsymbol{\psi}}(\boldsymbol{h},t+\delta
t;\theta_{K})\right\|_{F}$ =
---
$=\left\|\boldsymbol{\psi}_{\text{NN}}(\boldsymbol{h},t+\delta
t;\theta_{K})-\boldsymbol{\psi}_{K}(\boldsymbol{h},t+\delta
t;\theta_{K})\boldsymbol{\psi}_{K}(\boldsymbol{h},t;\theta_{K})^{+}\boldsymbol{\psi}_{\text{NN}}(\boldsymbol{h},t;\theta_{K})\right\|_{F}$
(8)

where $\left\|\cdot\right\|_{F}$ refers to the Frobenius norm. Here,
$\theta_{K}$ stands for the weights of the neural networks. We take advantage
of automatic differentiation, and we compute the gradient of
$\partial\mathcal{L}_{K}/\partial\theta_{K}$ directly, enabling to find
$\mathsfbi{K}$ and the set of observables $\boldsymbol{\Psi}(\boldsymbol{h})$
simultaneously using the Adam optimizer (Kingma & Ba, 2014). For more details,
we refer to Constante-Amores et al. (2024), in which this Koopman methodology
is extensively explained.

For time evolution predictions, a ‘pure’ Koopman approach would simply apply
Eq. 4. (In Constante-Amores et al. (2024), this method was denoted KDLoo).
While respecting the linearity of the exact infinite-dimensional formulation,
this approach struggles to produce accurate results over long time horizons
even for simple systems, not to mention complex chaotic ones (Constante-Amores
et al., 2024). Essentially, small errors accumulate and a drift of the
dynamics are observed. Nevertheless, it has been found to work well for
systems featuring steady and oscillatory behaviours – see Constante-Amores et
al. (2024); Constante-Amores & Graham (2024).

A slight modification of Eq. 4 has been found in past work (Li et al., 2017;
Junker et al., 2022; Van Goor et al., 2023; Constante-Amores et al., 2024) to
generate substantially better time-evolution predictions for simple systems,
and we show below that it also works extremely well for complex chaotic
systems. In this approach, at each time step the state $\boldsymbol{h}(t)$ is
used to construct $\boldsymbol{\Psi}(\boldsymbol{h}(t))$, which is then
evolved according to Eq. 4 to yield $\boldsymbol{\Psi}(\boldsymbol{h}(t+\delta
t))$. Finally $\boldsymbol{h}(t+\delta t)$ is obtained by projection – recall
that the first $d_{h}$ elements of $\boldsymbol{\Psi}$ are simply the state
(in manifold coordinates) $\boldsymbol{h}$. In summary, one time step of this
approach is given as

$\boldsymbol{h}(t+\delta
t)=\mathsfbi{P}\mathsfbi{K}\boldsymbol{\Psi}(\boldsymbol{h}(t))$ (9)

where $\mathsfbi{P}$ represents a projection operator that maps down from the
full vector of observables to just the state. We refer to this approach as
‘projected Koopman dynamics’. (In Constante-Amores et al. (2024), this method
was denoted KDLso). Because $\boldsymbol{\Psi}(\boldsymbol{h})$ is nonlinear,
this is a nonlinear evolution equation for $\boldsymbol{h}$. (It is worth
noting that other Koopman-based formalisms for time evolution such as that of
Khodkar et al. (2019) are also nonlinear.)

Although this approach originated in the space of observables, in the end it
is simply a nonlinear discrete-time evolution equation in the state space of
$\boldsymbol{h}$ with a particular structure. It is a composition of a
nonlinear mapping to a high dimensional observable space given by
$\boldsymbol{\Psi}$, a linear evolution of $\boldsymbol{\Psi}$ with
$\mathsfbi{K}$ according to the Koopman formalism, and projection back from
the observable space to the manifold coordinates with $\mathsfbi{P}$. In the
deep learning context, it is a nonlinear neural network followed by two linear
layers. In combination with dimension reduction onto an invariant manifold, we
refer to this framework as ‘Koopman Data-driven Manifold Dynamics’ (KDManD).

To find an effective architecture for lifting the state space to the
observable space, we use a trial-and-error approach, varying the network
architecture, the number of dictionary elements, and the activation functions.
Table 1 summarises the network parameters that yielded the best results. In
this work, we select 200 dictionary elements to lift the state space to the
observable space. To optimize the loss function described in equation 8, we
use an Adam optimiser in pythorch (Paszke, 2019). We select a batch size of
5000, and we train the models for 5000 epochs using a learning rate scheduler
that systematically decreases its value by three orders of magnitude, from
$10^{-4}$ to $10^{-6}$ over an equivalent span of three epochs. At this point,
we see no improvement in the reconstruction error.

###  2.3 state space modelling: Neural ODEs

We use a stabilized neural ODE framework for state space modeling in the
latent space. Rather than Eq. (2), we use a slight modification:

$\frac{d\boldsymbol{h}}{dt}=\boldsymbol{g}(\boldsymbol{h})-\mathsfbi{A}\boldsymbol{h}.$
(10)

Here ${A}_{ij}=-\beta\delta_{ij}\sigma_{i}(\boldsymbol{h})$,
$\sigma_{i}(\boldsymbol{h})$ is standard deviation of the $i$th component of
$\boldsymbol{h}$, $\beta$ is a fixed parameter and $\delta_{ij}$ is the
Kronecker delta. This modification, with a small value of $\beta$ it has been
found to stabilize the dynamics against spurious growth of fluctuations
without compromising accuracy of predictions Linot & Graham (2022); Linot et
al. (2023a); Floryan (2024). Representing the vector field $\boldsymbol{g}$ on
the manifold as a neural network with weights $\theta_{f}$, we can time-
integrate Eq. 10 between $t$ and $t+\delta t$ to yield a prediction
$\tilde{\boldsymbol{h}}(t+\delta t)$, such as

$\tilde{{\boldsymbol{h}}}(t+\delta t)=\boldsymbol{h}(t)+\int_{t}^{t+\delta
t}\left(\boldsymbol{g}(\boldsymbol{h}(t^{\prime});\theta_{f})+\mathsfbi{A}\boldsymbol{h}(t^{\prime})\right)dt^{\prime}.$
(11)

Given data for $\boldsymbol{h}(t)$ and $\boldsymbol{h}(t+\delta t)$ for a long
time series we can train $\boldsymbol{g}$ to minimize the $L_{2}$ difference
between the prediction $\tilde{\boldsymbol{h}}(t+\delta t)$ and the known data
$\boldsymbol{h}(t+\delta t)$. Then the loss function is defined by

$\mathcal{L}_{g}=\frac{1}{dN}\sum_{i=1}^{N}||\boldsymbol{h}(t_{i}+\delta
t)-\tilde{\boldsymbol{h}}(t_{i}+\delta t)||_{2}^{2}.$ (12)

We use automatic differentiation to determine the derivatives of
$\boldsymbol{g}$ with respect to $\theta_{f}$. To optimize the loss function
described in equation 12, we use an Adam optimiser in pythorch (Paszke, 2019).

##  3 Results

###  3.1 MFE model

####  3.1.1 Description of MFE data

The MFE model, developed by Moehlis et al. (2004), is a low-dimensional model
that describes the dynamics of an incompressible fluid between infinite, free-
slip walls driven by a spatially sinusoidal body force. The domain is periodic
of size $L_{x}\times L_{y}\times L_{z}$, with infinite, parallel walls at
$y=\pm L_{y}/2$ and periodic boundaries in $x$ and $z$. Here, $x$, $y$, and
$z$ are the streamwise, wall-normal, and spanwise coordinates, respectively.
The model displays features consistent with turbulence in the transition
regime, with long periods of turbulent behavior with infrequent quasi-
laminarization events (also called hibernating (Xi & Graham, 2010) or
quiescent (Hamilton et al., 1995) intervals) and ultimately full
laminarization at long times (Hamilton et al., 1995; Xi & Graham, 2010; Hof et
al., 2006). The model is an approximation to the Navier-Stokes equations,
generated through a severe truncation of the Fourier Galerkin method, that
represents the flow as combinations of spatial Fourier modes
$\boldsymbol{v}_{i}(\boldsymbol{x})$. The velocity field $\boldsymbol{v}$ is
given by a superposition of nine models

$\boldsymbol{v}(\mathbf{x},t)=\sum\limits^{9}_{i=1}a_{i}(t)\boldsymbol{v}_{i}(\mathbf{x})$
(13)

with modes describing the laminar flow profile and turbulent deviations caused
by streaks, vortices, and the interactions between them. The dynamics of the
system are captured through the evolution of the nine modal amplitudes
$a_{i}(t)$, governed by nine ordinary differential equations given by the
Galerkin projection, which can be found explicitly in Moehlis et al. (2004).
For this case, we perform no dimension reduction, so
$\boldsymbol{x}=\boldsymbol{h}=[a_{1},a_{2},\ldots,a_{9}]^{T}$.

Domain dimensions of $L_{x}=4\pi$, $L_{y}=2,L_{z}=2\pi$, with a channel
Reynolds number of $Re=400$, have been shown to produce turbulent flow
dynamics with sufficient self-sustainment to be captured through data-driven
modeling (Srinivasan et al., 2019). We generated 100 unique time series as
training data by integrating the MFE equation using a fourth-order Runge-Kutta
method with a time step of 0.5 (i.e., $\delta t=0.5$). Each time series
captures the transient turbulent state, featuring intervals of turbulence
interrupted by quasi-laminarization events, ultimately leading to full
laminarization at long time. The flow is often characterized using the total
kinetic energy (KE), defined as $KE=\frac{1}{2}\sum_{i=1}^{9}a_{i}^{2}$. In
this context, the turbulent state corresponds to low energy, while the laminar
state is high energy. After a long turbulent transient, every time series
eventually evolves to the laminar fixed point $a_{i}=\delta_{i1}$. The time
series were generated with initial conditions for eight of the amplitudes set
as follows:
$(a_{1},a_{2},a_{3},a_{5},a_{6},a_{7},a_{8},a_{9})=(1,0.07066,-0.07076,0,0,0,0,0)$.
The initial value of $a_{4}$ was randomly chosen within the range
$a_{4}\in[-0.1,0.1]$. These initial conditions have been shown to produce
chaotic dynamical data with quasi-laminarization events (Eivazi et al., 2021).
The results are presented with time normalized using the Lyapunov time,
$\tau_{L}$, for the system (e.g., $\tau_{L}\approx 61$ according to Racca &
Magri (2022)). The first 1000 time steps of each time series were discarded to
eliminate dependence on initial conditions. Each time series was evolved until
the laminar state was reached, requiring a variable number of time steps
depending on the initial conditions. The resulting training data set comprised
approximately 2 million snapshots of MFE amplitudes, with a mean lifetime of
$164\tau_{L}$. An example of the evolution of three of the amplitudes from a
trajectory, as well as the total kinetic energy of the trajectory,
$E_{k}=\frac{1}{2}\sum_{1}^{9}a_{i}^{2}$, is shown in Figure 1a.

####  3.1.2 Koopman MFE model

|
---|---
(a) | (b)
Figure 4: (a) Eigenvalues of the Koopman operator for MFE model (b) Temporal
evolution of three amplitudes $a_{1}$, $a_{6}$ and $a_{9}$ for the true MFE
solution (solid lines) and Koopman predictions (dashed lines) for the same
initial condition.

First, we apply Projected Koopman Dynamics to the MFE model described in
section 3.1.1 with $\delta t=0.5$. Figure 4a shows the eigenvalues of the
Koopman operator with 200 dictionary elements. We observe eigenvalues within
the unit circle (i.e., $|\lambda_{k}|\approx 0.99$). If we use the linear
evolution described by equation 4, given that the $|\lambda_{k}|<1$, the
dynamics will decay. But, by using the projection, described by Eq. 9, we
achieve a good approximation of the true dynamics.  Figure 4b shows the
temporal evolution of three amplitudes $a_{1}$, $a_{6}$ and $a_{9}$ for the
true MFE solution and Koopman predictions for the same initial condition up to
$10\tau_{L}$. The predictions show an exceptional agreement with the true
solution up to 3$\tau_{L}$, then the dynamics diverge due to their chaotic
nature.

The success of our Koopman-based model is first assessed in its ability to
reconstruct the time-averaged velocity statistics of the turbulent regime of
the MFE model, demonstrating a base level of competence vital for any data-
driven modeling technique. Using 100 randomized initial conditions, we used
our Koopman model to forecast the evolution of the modal amplitudes over a
period of 100 $\tau_{L}$. These same initial conditions were similarly evolved
using the MFE equations over the same time period. The true and predicted
amplitudes were then projected onto the corresponding spatial modes, and the
time-averaged velocity statistics of the corresponding model were calculated.

Figure 5a shows the mean streamwise velocity and Reynolds shear stress from
the resulting true and reconstructed flow fields, as well as the difference
between these two measurements. The Koopman-based model captures the structure
and magnitude of the velocity statistics extremely well, with the deviation
from the true values remaining comparatively small throughout the spatial
domain. The Koopman model outperforms the model from Fox et al. (2023),
reducing the error in the forecasts of both the velocity profile and the
Reynolds shear stress.

The Koopman-based model was next assessed on its ability to precisely forecast
the short-time evolution of the modal amplitudes in the turbulent regime.
Using 1000 random and unique initial conditions, we forecasted the evolution
of the MFE modes using our Koopman model for a time period of 10 $\tau_{L}$,
with the same initial conditions evolved separately by the MFE equations over
the same period. The accuracy of our forecasts was then evaluated by
calculating a time-dependent, relative ensemble-averaged error over the
testing time period
$E(t)=\left<||a(t)-\tilde{a}(t)||_{2}\right>/\left<||a(t)||_{2}\right>$ (see
figure 5b). The Koopman-based model performs very well, accurately forecasting
the evolution of the modal amplitudes for nearly one $\tau_{L}$ with minimal
error. The ensemble-averaged error of the forecasts grows very slowly over the
testing period, requiring several Lyapunov times to completely deviate from
the true amplitude evolution, further demonstrating the short-time accuracy of
the Koopman-based model. The Koopman model again outperforms the results of
Fox et al. (2023), with the new data-driven model reducing the error growth
significantly in the testing window.

|
---|---
(a) | (b)
Figure 5: (a) The top left and right panels display the mean streamwise
velocity and Reynolds shear stresses for the full field of both the test and
training data, respectively. The bottom panels show the differences between
the true values and the predicted values. Results from Fox et al. (2023) are
also included. (b) Ensemble averaged short-time error tracking of the
reconstruction of the Koopman model, with comparison to prior study by Fox et
al. (2023).

We next determined the ability of the Koopman-based model to predict the long-
time structure and occurrence of quasi-laminarization events in turbulent
regimes. In the MFE model, quasi-lamizarization is characterized by an
increase in the first modal amplitude, associated with the mean shear flow
profile, and a decrease in the amplitude of the remaining eight modes, which
correspond with the turbulent fluctuations. After a conditionally-dependent
time, these values collapse back towards their original range as the flow
evolves back into the turbulent regime. To observe the frequency and structure
of these events, we generated a test data set containing 100 trajectories
evolved with the MFE model for 100 Lyapunov times, separately evolving the
initial conditions of these trajectories with the Koopman-based model; as we
are here interested solely in the turbulent behavior, all full re-
laminarizations were removed from both data sets. Shown in Figure 6 is the
joint probability density function (PDF) of $a_{1}$ and $a_{3}$ for the true
and reconstructed trajectories, where the quasi-laminarization events can be
seen in the narrowing tail on the right side of the plots. The reconstruction
by the Koopman mode closely matches that of the true data, adequately
capturing the turbulent and quasi-laminar states. The accuracy is further
quantified by the relative mean squared error,
$rMSE=MSE(PDF_{true}-PDF_{pred})/\mathcal{P}-1$, calculated by taking the mean
squared error (MSE) between the two joint PDFs and normalizing it by the MSE
between two joint PDFs generated from true trajectories of the same length
with different initial conditions, $\mathcal{P}$, and subtracting one; here,
the relative MSE is 0.1, shown inset in the middle panel of Figure 6, a very
close match and an improvement over a previous study using the state-space
based CANDyMan approach (Fox et al., 2023).

Figure 6:  Joint probability density function (jPDF) of $a_{1}$ and $a_{3}$
for the true data, the Koopman model, and the work by Fox et al. (2023),
corresponding to left-to-right panels, respectively. For the models, we have
added the corresponding value of the relative mean squared error (rMSE). Note
the logarithmic scale.

Finally, we assessed the ability of the data-driven model to forecast the
occurrence of a quasi-laminarization event, here defined as an increase in the
kinetic energy of the system such that $E_{k}>0.1$. Using 1000 initial
conditions, we used the MFE equations and the Koopman model to evolve the
modal amplitudes separately for a period of 5.5 $\tau_{L}$. We then segmented
the trajectories into 0.5 $\tau_{L}$ time windows, where the window at time
$\tilde{t}=t/\tau_{L}$ refers to the time between $\tilde{t}$ and
$\tilde{t}+0.5\tau_{L}$. Each window was then investigated for a quasi-
laminarization event in both the true and forecasted data sets for each
trajectory, and the two observations were compared. If an extreme event
occurred in both the true and reconstructed trajectories, this was labeled as
a true positive ($TP$). If the true trajectory exhibited an extreme event, but
the Koopman model failed to predict one, this was labeled as a false negative
($FN$). If the model forecasted an extreme event when the true trajectory
showed none, it was identified as a false positive ($FP$) (Racca & Magri,
2022). The total number of each identification type in each window was
tabulated and used to calculate the precision $p$, recall $r$, and F-score
$F$, where $p=TP/(TP+FP)$, $r=TP/(TP+FN)$, and $F=2/(p^{-1}+r^{-1})$.
Therefore, $p=1$ indicates that all predicted events did occur, $r=1$
indicates that all events were predicted, and $F=1$ indicates that both are
true. Figure 7 displays the value of each metric at each time window over the
testing period, as well as a comparison to the results of Racca & Magri
(2022), which used an echo state network to generate a data-driven model of
the system, as well as to our previous study Fox et al. (2023). The Koopman
model performs very well, bettering the forecasting ability of the model by
Fox et al. (2023), and equaling the forecasting ability of Racca & Magri
(2022) at all time windows across the testing period. Remarkably, the Koopman
model can achieve this despite being Markovian, while the model of Racca &
Magri (2022) is notably non-Markovian (as they used an echo state network).

Finally, when evaluating long-term full laminarization, the Koopman models do
not transition to the laminar state, instead maintaining the dynamics in the
turbulent state. We note that the initial conditions of the training dataset
described in section 3.1.1 are far from the laminar state, and subsequently,
only a small portion of the data reflects the dynamics in the laminar state.
To address this issue, a new training dataset with initial conditions closer
to the laminar state should be included during training of the Koopman model.
We note that Pan & Duraisamy (2024) discussed the challenges of applying
Koopman methods to systems with multiple attractors, which will also apply to
systems with long chaotic transients (‘chaotic saddles’).

Figure 7: Precision, recall, and F-score of extreme event forecasting of the
MFE model by the Koopman model, with comparison to prior studies by Fox et al.
(2023) and Racca & Magri (2022).

###  3.2 Plane Couette Flow

####  3.2.1 Description of PCF data

We perform DNS of an incompressible viscous fluid inside of two plates moving
in opposite directions at the same speed. A schematic representation of the
numerical setup is shown in Figure 1b. The dimensionless form of the
continuity and momentum equations are respectively expressed as

$\nabla\cdot\boldsymbol{v}=0,$ (14) $\frac{\partial\boldsymbol{v}}{\partial
t}+\boldsymbol{v}\cdot\nabla\boldsymbol{v}=-\nabla
p+Re^{-1}\nabla^{2}\boldsymbol{v},$

here $t$, $\boldsymbol{v}=[v_{x},v_{y},v_{z}]$ and $p$ stand for the time, the
velocity and the pressure, respectively. Here, $Re$ is the Reynolds number,
defined as $Re=Ud/\mu$, where $U,d$ and $\mu$ stand for the plate velocity,
the distance between plates and the kinematic viscosity, respectively. We
consider a minimal flow unit (MFU) of the PCF system, and the computational
cell follows the domain presented by Hamilton et al. (1995) with $Re=400$ and
a domain size $[L_{x},L_{y},L_{z}]=[1.75\pi,2,1.2\pi]$; where,
$x\in[0,L_{x}]$, $y\in[-1,1]$ and $z\in[0,L_{z}]$ for the streamwise, wall-
normal and spanwise directions, respectively. The data was generated using the
pseudo-spectral Channelflow code on a grid of dimensions
$[N_{x},N_{y},N_{z}]=[32,35,32]$ in the $x$, $y$, and $z$ directions,
respectively (Gibson, 2012; Linot & Graham, 2023). We impose periodic boundary
conditions in $x$ and $z$ and no-slip and no-penetration boundary conditions
in $y$.

Equation 3.2.1, under the boundary conditions described above, is invariant
(and its solutions equivariant) under the discrete symmetries of point
reflections about $[x,y,z]=[0,0,0]$:

$\mathcal{P}[(v_{x},v_{y},v_{z},p)(x,y,z,t)]=(-v_{x},-v_{y},-v_{z},p)(-x,-y,-z,t),$
(15)

reflection about the $z=0$ plane:

$\mathcal{R}[(v_{x},v_{y},v_{z},p)(x,y,z,t)]=(v_{x},v_{y},-v_{z},p)(x,y,-z,t),$
(16)

and rotation by $\pi$ about the z-axis:

$\mathcal{RP}[(v_{x},v_{y},v_{z},p)(x,y,z,t)]=(-v_{x},-v_{y},v_{z},p)(-x,-y,z,t).$
(17)

In addition to the discrete symmetries, there are also continuous translation
symmetries in $x$ and $z$.

$\mathcal{T}_{\sigma_{x},\sigma_{z}}[(v_{x},v_{y},v_{z},p)(x,y,z,t)]=(v_{x},v_{y},v_{z},p)(x+\sigma_{x},y,z+\sigma_{z},t).$
(18)

We incorporate all these symmetries into the POD representation of the data
(as done by Smith et al. (2005); Linot & Graham (2023)). Then, we apply the
method of slices to achieve phase alignment of the flow in the spanwise ($z$)
direction (Budanur et al., 2017). This phase alignment in $z$ allows us to fix
the position of the low-speed streak in the center of the domain, enhancing
the predictive capability of the reduced order models.

Different initial conditions are chosen, and after removing the transient data
and the last time to avoid laminar states, we accumulate a total of 91,562
time units of data which live in the turbulent attractor, stored at intervals
of one time unit (i.e., $\delta t=1$). The data was then divided into $80\%$
for training and $20\%$ for testing purposes.

|
---|---
(a) | (b)

---
(c)

Figure 8: Identification of the manifold coordinates. (a) Normalised singular
values of the covariance of the latent space data for five different trials.
(b) Reconstruction of $\left\langle\left\|a_{n}\right\|^{2}\right\rangle$
(mean-squared POD coefficient amplitudes) for the test data from 502 POD modes
and the IRMAE-WD-B with $d_{h}=36$. (c) Reconstruction of the Reynolds
stresses for the DNS and IRMAE-WD-B with $d_{h}=36$.

####  3.2.2 Finding Manifold Coordinates

The initial step in building the state space or function space models involves
reducing the dimensionality of the full state ($\mathcal{O}(10^{5})$) to its
manifold coordinates ($\mathcal{O}(10^{1})$). This step is crucial because
lifting the state space to a high-dimensional observable space, would be
computationally intractable using the full space. The first step in building
this low-dimensional model is to perform a linear reduction with POD. We
select 256 modes, resulting in a vector of POD coefficients $a(t)$, which
contains $99.8\%$ of the energy. Projecting onto these modes gives rise to a
502-dimensional system. For more details about the POD reduction, we refer to
(Linot & Graham, 2023), in which the POD framework was extensively explained
for MFU PCF. After converting the data into POD coefficients and eliminating
the low-energy modes, the next step involves training IRMAE-WD-B for nonlinear
dimension reduction.

Figure 8a shows the normalized singular values of the covariance of the latent
space data for five different trials. It shows that the singular values drop
several orders of magnitude in the range of $32<d_{h}<36$, with some drops
reaching magnitudes as low as $\mathcal{O}(10^{-14})$. The relative error in
the POD coefficients across all the trials, measured as
$\left\langle\left\|a-\tilde{a}\right\|\right\rangle/\left\langle\left\|a\right\|\right\rangle\approx
0.06$. To demonstrate that the latent space spanned for IRMAE-WD-B is capable
of reconstructing accurately the POD coefficients, Figure 8b shows the
$\left\langle\left\|a_{n}\right\|^{2}\right\rangle$ (mean-squared POD
coefficient amplitudes) for the test data from 502 POD modes and the IRMAE-
WD-B with $d_{h}=36$. The choice of $d_{h}=36$ reflects the highest predicted
latent dimension among the five trials, representing a conservative choice.
Finally, Figure 8c shows the Reynolds stresses of the test data, and those
from the truncated latent dimension at $d_{h}=36$, and we observe an
exceptional agreement between them.

Across the five independent trials, we observed that IRMAE -WD-B does not
follow a consistent path; instead, it demonstrates an alternation of varying
numbers of linear layers across all trials. This variability suggests that a
degree of randomness may contribute to improving the dimension estimate. We
note that Linot & Graham (2023) found that a latent space with a
dimensionality of $d_{h}=18$ is capable of faithfully capturing the short-time
tracking and long-time statistics of the system, including the leading
Lyapunov exponents, Reynolds stresses, and Power-Input dissipation. To
determine this dimension, they train autoencoders until the MSE plateaus, and
then train NODE at different dimensions. They did not observe quantitative
improvements beyond $d_{h}=18$.

####  3.2.3 State and function space modelling in manifold coordinates

|
---|---
(a) | (b)

---
(c)

Figure 9: (a) Eigenvalues of the Koopman operator. (b) Kinetic energy ($E_{k}$) of the flow for the DNS, KDManD and DManD with $d_{h}=36$ up to $t/\tau_{L}=3$. (c) Two-dimensional representation of the dynamics in a $z$-$x$ plane ($y=0$) with $u_{x}$ for the DNS and the models. Time is displayed at the top of each panel and corresponds to the marked points in panel (b). |
---|---
(a) | (b)
Figure 10: (a) Normalised short-time-tracking error for 10 random initial
conditions using KDManD and DManD up to 150 time units. (b) Temporal
autocorrelation of the fluctuations of the kinetic energy for KDManD, DManD
and DNS.

Now, we apply the state space and function space approaches described in
section 2 (with $\delta t=1$) to learn a linear mapping of the dynamics in the
manifold coordinates, and a nonlinear function, respectively. We select 200
dictionary elements to find the observable space. Similarly to what was
observed for the Koopman MFE model, by using the projection for temporal
predictions, described by Eq. 9, we achieve a good approximation of the true
dynamics .

We start the discussion of the predictive capabilities of KDManD and DManD by
selecting an initial condition that displays the streak breakdown process
characteristic of wall turbulence. To evaluate how close the models can
capture the dynamics of the system, Figure 9b represents the kinetic energy of
the flow:

$E_{k}=\frac{1}{2L_{x}L_{z}}\int_{0}^{L_{x}}\int_{-1}^{1}\int_{0}^{L_{z}}\boldsymbol{v}^{2}~{}dx~{}dy~{}dz.$
(19)

We evolve the same initial condition with the DNS solver and in the models. We
observe that KDManD can faithfully capture the dynamics, with a good match up
to $3\tau_{l}$ Lyapunov times (here the Lyapunov time corresponds to around
$\tau_{l}\approx 48$ time units, reported by (Inubushi et al., 2015)). Figure
9c displays field snapshots in the $x-z$ plane ($y=0$) up to $3\tau_{L}$ for
the DNS (top panels) and KDManD (bottom panels). We observe that streaks
(velocity functions below the mean) at the center of the domain (see
$t/\tau_{L}=0$) are affected by spanwise wavy disturbances, leading to their
breakdown and the formation of rolls (see $t/\tau=0.52-1.04$). These rolls, in
turn, induce the re-emergence of streaks (see $t/\tau=2.08$), which break down
(see panels $t/\tau=2.39-3.12$) again to restart the sustaining process of
turbulence.

While Figure 9b,c depicts the trajectory for a single selected initial
condition, Figure 10a illustrates t ensemble relative tracking error in the
POD coefficients for 500 randomly chosen initial trajectories with KDManD and
DManD. Here, the relative error of the models remains below one until one
$\tau_{L}$, indicating that the models can reasonably capture the dynamics for
approximately one Lyapunov time with approximately $80\%$ of relative error at
one $\tau_{L}$. We highlight the small error at $t/\tau_{L}=0$, attributable
to the effective coordinate transformation that maps initial conditions from
the full space to the manifold coordinates. Figure 10b illustrates the
temporal correlation of the fluctuating part of the kinetic energy, defined as
$k(\tau)=E(\tau)-\langle E\rangle$. To generate this figure, we utilized 1000
different initial conditions evolved up to $3\tau_{L}$. We observe that KDManD
matches well for one $\tau_{L}$, with small discrepancies persisting up to
$3\tau_{L}$.

Next, we focus on the performance of KDManD and DManD to capture long-time
dynamics of the system. Thus, we will evaluate their performance in the
reconstruction of the Reynolds stresses, and the prediction of the energy
balance at long times. For the energy balance, we calculate the power input
($I$) injected due to its moving plates and the dissipation due to viscosity
($D$), defined as

$I=\frac{1}{2L_{x}L_{z}}\int_{0}^{L_{x}}\int_{0}^{L_{z}}\left|\frac{\partial
v_{x}}{\partial y}\right|_{y=-1}+\left|\frac{\partial v_{x}}{\partial
y}\right|_{y=1}~{}dx~{}dz,$ (20)
$D=\frac{1}{2L_{x}L_{z}}\int_{0}^{L_{x}}\int_{-1}^{1}\int_{0}^{L_{z}}|\boldsymbol{\nabla}\times\boldsymbol{v}|^{2}~{}dx~{}dy~{}dz.$
(21)

On an energy balance, the energy input rate, the dissipation rate, and the
kinetic energy $E=(I-D)/Re$, must average to zero over long times.

Figure 11 summarises the long-time statistics of the DNS and the models for an
initial condition (IC) that has been evolved for approximately $40\tau_{L}$.
Figure 11a displays four components of the Reynolds stresses, KDManD
demonstrates remarkable accuracy in predicting the DNS results, showing a very
small deviation from the true solutions, some small discrepancies are observed
in DManD model for the $\left\langle v^{\prime 2}_{z}\right\rangle$ component.
To ensure that the system effectively captures the energy balance, Figure 11b
displays the $D-I$ projection of the states predicted by KDManD and DManD,
with black shading indicating a long trajectory from the DNS solver. It is
observed that KDManD and DManD models effectively capture the core region of
these projections, with some excursions occurring at low dissipation rates
that are not present in the DNS results. Finally, Figure 11c shows the PDF
from the DNS (generated from a single trajectory evolved for $100\tau_{L}$),
the PDF of KDManD and DManD models. Overall, we can conclude that KDManD can
effectively predict accurately the long-time statistics of this complex system
in manifold coordinates.

Finally, we note that at very long times (beyond $40\tau_{L}$), most KDManD
models encounter stability problems. This is a current limitation of the
framework presented in section 2.2. We note that this issue might be resolved
with a stabilization term, as employed in the DManD framework, but addressing
this is beyond the scope of the present research. We did not observe this
stability issue in the MFE model. We believe this is due to the low
dimensionality of the system’s state space.

|
---|---
(a) | (b)

---
(c)

Figure 11: Long-time statistics of the models and the DNS. The same IC has
been evolved for $40\tau_{l}$ for the models: (a) Components of the Reynolds
stress. (b) Energy dissipation vs power input. (c) PDF plots of the
dissipation vs power input.

##  4 Conclusions

Simulations of shear flows play a crucial role in the study of turbulence as
they provide a simplified setting to explore fundamental mechanisms of
turbulence, such as the self-sustaining process of wall turbulence (Hamilton
et al., 1995). Most of the works applied for the modelling of shear flows are
based on the state space approach. In the present work, we investigate shear
flows using both state space and function space approaches. For the latter, we
lift the state space to a higher dimensional observable space where  the
temporal evolution of the observables can be captured through the Koopman
operator. To find this set of observables (we use neural networks as
‘universal approximators’) and the corresponding Koopman operator, we employ a
variant of EDMD-DL which leverages automatic differentiation to perform
gradient descent through the Moore-Penrose pseudoinverse (Constante-Amores et
al., 2024).

First, we applied the Koopman methodology to the MFE model (Moehlis et al.,
2004). Since the model is inherently low-dimensional, and only described by
nine coefficients, we applied the the Koopman framework directly in the state
space. Our data-driven Koopman model accurately predicts both the short-time
and long-time statistics of the system. We compared our results with (Fox et
al., 2023) (which used the CANDyMan framework presented in Floryan & Graham
(2022)), and the echo state results presented by Racca & Magri (2022). Our
Koopman framework yielded better short-time predictions (accurate for almost 3
Lyapunov times) and long-time predictions, including the reconstruction of
Reynolds stresses, the mean velocity profile, and the joint PDF of
coefficients $a_{1}$ and $a_{3}$.

Next, we moved to the high dimensional Minimal Flow Unit of plane Couette flow
at $Re=400$. To ensure all the relevant spatial and temporal scales are fully
resolved, the initial dimensionality of the full state space is
($\mathcal{O}(10^{5})$. Lifting this space to a higher observable space for
Koopman predictions renders this approach computationally intractable. But,
leveraging the dissipative nature of the Navier-Stokes equations, where
viscosity attenuates small-scale perturbations, results in the long time
dynamics living in an invariant manifold of lower dimension. To identify a
coordinate transformation to a manifold representation of
($\mathcal{O}(10^{1})$), we introduced the IRMAE-WD-B framework, which
provides a consistent estimate of the manifold dimension ranging between
$d_{h}=32$ and $d_{h}=36$; we selected the latter dimension for the function
space and state space modeling. Note that this dimension is larger than the
one predicted by Linot & Graham (2023), who found $d_{h}=18$ by varying the
latent dimension and evaluating the mean squared error of the reconstruction
using another type of advanced autoencoder (Linot & Graham, 2020). Then, we
lifted the manifold state space to a high-dimensional observable space where
the dynamics can be described by the linear Koopman operator. Our Koopman
framework successfully captured the streak breakdown process, tracked short-
term dynamics up to nearly one Lyapunov time, and accurately reconstructed the
long-term statistics evaluated in the Reynolds stresses and the power input-
dissipation energy balance. Additionally, we included results from state space
modeling using neural ODEs to learn the evolution equation in manifold
coordinates. We observed comparable results between the function space and
state space frameworks.

A promising direction for future research is leveraging the linear
representation afforded by the Koopman operator. This linearization simplifies
control, enabling the use of established techniques such as the Linear
Quadratic Regulator (LQR) for optimal control in nonlinear systems. We are
actively developing model-based control strategies within manifold dynamics
and extending them to the full state space (similar to the work of Linot et
al. (2023b)). Integrating LQR with the Koopman framework could significantly
enhance control precision for nonlinear systems by transforming their dynamics
into a higher-dimensional linear framework. This approach balances performance
and energy usage effectively and could advance applications from fluid
dynamics to robotic motion planning.

## Declaration of Interests

The authors report no conflict of interest.

## Acknowledgments

This work was supported by ONR N00014-18-1-2865 (Vannevar Bush Faculty
Fellowship). We gratefully acknowledge Alec J. Linot for helpful discussions.

## Data availability statement

Code and sample data that support the findings of this study for finding
manifold coordinates with IRMAE-WD-B are openly available at
https://github.com/mdgrahamwisc/IRMAE_WD_B

## References

  * Abdi & Williams (2010) Abdi, H. & Williams, L. J 2010 Principal component analysis. Wiley interdisciplinary reviews: computational statistics 2 (4), 433–459.
  * Alford-Lago et al. (2022) Alford-Lago, D. J., Curtis, C. W., Ihler, A. T. & Issan, O. 2022 Deep learning enhanced dynamic mode decomposition. Chaos 32 (3), 033116.
  * Avila et al. (2023) Avila, M., Barkley, D. & Hof, B. 2023 Transition to turbulence in pipe flow. Annu. Rev. Fluid Mech. 55 (1), 575–602.
  * Budanur et al. (2017) Budanur, N. B., Short, K. Y., Farazmand, M., Willis, A. P. & Cvitanović, P. 2017 Relative periodic orbits form the backbone of turbulent pipe flow. J. Fluid Mech. 833, 274–301.
  * Chen et al. (2019) Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. 2019 Neural ordinary differential equations. arXiv:1806.07366 .
  * Constante-Amores & Graham (2024) Constante-Amores, C. R. & Graham, M. D. 2024 Data-driven state-space and Koopman operator models of coherent state dynamics on invariant manifolds. J. Fluid Mech. 984, R9.
  * Constante-Amores et al. (2024) Constante-Amores, C. R., Linot, A. J. & Graham, M. D. 2024 Enhancing predictive capabilities in data-driven dynamical modeling with automatic differentiation: Koopman and neural ODE approaches. Chaos 34 (4), 043119.
  * Eivazi et al. (2021) Eivazi, H., Guastoni, L., Schlatter, P., Azizpour, H. & Vinuesa, R. 2021 Recurrent neural networks and Koopman-based frameworks for temporal predictions in a low-order model of turbulence. Int J Heat Fluid FL 90, 108816.
  * Floryan (2024) Floryan, D. 2024 On instabilities in neural network-based physics simulators. arXiv: 2406.13101 .
  * Floryan & Graham (2022) Floryan, D. & Graham, M. 2022 Data-driven discovery of intrinsic dynamics. Nat. Mach. Intell 4 (12), 1113–1120.
  * Fox et al. (2023) Fox, A. J., Constante-Amores, C. R. & Graham, M. D. 2023 Predicting extreme events in a data-driven model of turbulent shear flow using an atlas of charts. Phys. Rev. Fluids 8, 094401.
  * Gibson (2002) Gibson, J. F. 2002 Dynamical systems models of wall-bounded, shear-flow turbulence. PhD thesis, Cornell University.
  * Gibson (2012) Gibson, J. F 2012 Channelflow: A spectral Navier-Stokes simulator in c++. New Hampshire 900.
  * Gibson et al. (2008) Gibson, J F, Halcrow, J & Cvitanović, P 2008 Visualizing the geometry of state space in plane Couette flow. J. Fluid Mech. 611, 24.
  * Gibson et al. (2009) Gibson, J. F., Halcrow, J. & Cvitanović, P. 2009 Equilibrium and travelling-wave solutions of plane Couette flow. J. Fluid Mech. 638, 243–266.
  * Halcrow et al. (2009) Halcrow, J., Gibson, J. F., Cvitanović, P. & Viswanath, D. 2009 Heteroclinic connections in plane couette flow. J. Fluid Mech. 621, 365–376.
  * Hamilton et al. (1995) Hamilton, J M, Kim, J & Waleffe, F 1995 Regeneration mechanisms of near-wall turbulence structures. J. Fluid Mech. 287, 317–348.
  * Hinton & Salakhutdinov (2006) Hinton, G. E. & Salakhutdinov, R. R. 2006 Reducing the dimensionality of data with neural networks. Science 313 (5786), 504–507.
  * Hof et al. (2006) Hof, B, Westerweel, J, Schneider, T M & Eckhardt, B 2006 Finite lifetime of turbulence in shear flows. Nature 443, 59–62.
  * Holmes et al. (2012) Holmes, P., Lumley, J. L., Berkooz, G. & Rowley, C. W. 2012 Turbulence, Coherent Structures, Dynamical Systems and Symmetry, 2nd edn. Cambridge Monographs on Mechanics . Cambridge University Press.
  * Inubushi et al. (2015) Inubushi, M., Takehiro, S. & Yamada, M. 2015 Regeneration cycle and the covariant lyapunov vectors in a minimal wall turbulence. Phys. Rev. E 92, 023022.
  * Itano & Generalis (2009) Itano, T. & Generalis, S. C. 2009 Hairpin vortex solution in planar couette flow: A tapestry of knotted vortices. Phys. Rev. Lett. 102, 114501.
  * Jiménez & Moin (1991) Jiménez, J. & Moin, P. 1991 The minimal flow unit in near-wall turbulence. J. Fluid Mech. 225, 213–240.
  * Jiménez (2013) Jiménez, K. 2013 Near-wall turbulence. Phys. Fluids 25 (10), 101302.
  * Jing et al. (2020) Jing, L., Zbontar, J. & others 2020 Implicit rank-minimizing autoencoder. Advances in Neural Information Processing Systems 33, 14736–14746.
  * Jolliffe (1986) Jolliffe, I. T. 1986 Principal Component Analysis and Factor Analysis, pp. 115–128. New York, NY: Springer New York.
  * Junker et al. (2022) Junker, A., Timmermann, J. & Trächtler, A. 2022 Data-driven models for control engineering applications using the koopman operator. In 2022 3rd International Conference on Artificial Intelligence, Robotics and Control (AIRC), pp. 1–9.
  * Kaszás et al. (2022) Kaszás, B., Cenedese, M. & Haller, G. 2022 Dynamics-based machine learning of transitions in couette flow. Phys. Rev. Fluids 7, L082402.
  * Kawahara & Kida (2001) Kawahara, G. & Kida, S. 2001 Periodic motion embedded in plane couette turbulence: regeneration cycle and burst. J. Fluid Mech. 449, 291–300.
  * Khodkar et al. (2019) Khodkar, M. A., Hassanzadeh, P. & Antoulas, A. 2019 A koopman-based framework for forecasting the spatiotemporal evolution of chaotic dynamics with nonlinearities modeled as exogenous forcings. arXiv:1909.00076 .
  * Kingma & Ba (2014) Kingma, D. & Ba, J. 2014 Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, ICLR 2015, 1–15 .
  * Koopman (1931) Koopman, B. O 1931 Hamiltonian systems and transformation in Hilbert space. PNAS 17 (5), 315–318.
  * Korda & Mezić (2018) Korda, M. & Mezić, I. 2018 On convergence of extended dynamic mode decomposition to the koopman operator. Journal of Nonlinear Science 28, 687–710.
  * Kramer (1991) Kramer, Mark A. 1991 Nonlinear principal component analysis using autoassociative neural networks. AIChE Journal 37 (2), 233–243.
  * Lasota & Mackey (1994) Lasota, A. & Mackey, M. C 1994 Chaos, Fractals and Noise: stochastic aspects of dynamics, 2nd edn. Springer . New York: Springer.
  * Li et al. (2017) Li, Q., Dietrich, F., Bollt, E. M. & Kevrekidis, I. G. 2017 Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator. Chaos 27 (10), 103111.
  * Linot et al. (2023a) Linot, A. J., Burby, J. W., Tang, Q., Balaprakash, P., Graham, M. D. & Maulik, R. 2023a Stabilized neural ordinary differential equations for long-time forecasting of dynamical systems. J. Comput. Phys. 474, 111838.
  * Linot & Graham (2020) Linot, A. J. & Graham, M. D. 2020 Deep learning to discover and predict dynamics on an inertial manifold. Phys. Rev. E 101, 062209.
  * Linot & Graham (2022) Linot, A. J. & Graham, M D. 2022 Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations. Chaos 32 (7), 073110.
  * Linot & Graham (2023) Linot, A. J. & Graham, M. D. 2023 Dynamics of a data-driven low-dimensional model of turbulent minimal Couette flow. J. Fluid Mech. 973, A42.
  * Linot et al. (2023b) Linot, A. J., Zeng, K. & Graham, M. D. 2023b Turbulence control in plane couette flow using low-dimensional neural ODE-based models and deep reinforcement learning. Int J Heat Fluid FL 101, 109139.
  * Liu et al. (2024) Liu, A., Axås, J. & Haller, G. 2024 Data-driven modeling and forecasting of chaotic dynamics on inertial manifolds constructed as spectral submanifolds. Chaos 34 (3), 033140.
  * Lusch et al. (2018) Lusch, B., Kutz, J. & Brunton, S. 2018 Deep learning for universal linear embeddings of nonlinear dynamics. Nat. Commun 9.
  * Mezic (2020) Mezic, I. 2020 Koopman operator, geometry, and learning, arXiv: 2010.05377.
  * Mezic (2023) Mezic, I 2023 Operator is the model, arXiv: 2310.18516.
  * Mezić (2013) Mezić, I. 2013 Analysis of fluid flows via spectral properties of the koopman operator. Annu. Rev. Fluid Mech. 45 (Volume 45, 2013), 357–378.
  * Milano & Koumoutsakos (2002) Milano, M & Koumoutsakos, P 2002 Neural network modeling for near wall turbulent flow. J. Comput. Phys. 182 (1), 1–26.
  * Moehlis et al. (2004) Moehlis, J., Faisst, H. & Eckhardt, B. 2004 A low-dimensional model for turbulent shear flows. New J. Phys 6 (1), 56.
  * Nagata (1990) Nagata, M. 1990 Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity. J. Fluid Mech. 217, 519–527.
  * Nakao & Mezic (2020) Nakao, H. & Mezic, I. 2020 Spectral analysis of the Koopman operator for partial differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (11), 113131.
  * Otto & Rowley (2019) Otto, S. E. & Rowley, C. W. 2019 Linearly recurrent autoencoder networks for learning dynamics. SIAM J. Appl. Dyn. Syst 18 (1), 558–593.
  * Pan & Duraisamy (2024) Pan, S. & Duraisamy, K. 2024 On the lifting and reconstruction of nonlinear systems with multiple invariant sets. Nonlinear Dynamics 112 (12), 10157–10165.
  * Paszke (2019) Paszke, et al. 2019 Pytorch: An imperative style, high-performance deep learning library. In NIPS (ed. H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox & R. Garnett), , vol. 32. Curran Associates, Inc.
  * Pérez-De-Jesús & Graham (2023) Pérez-De-Jesús, C. E & Graham, M. D. 2023 Data-driven low-dimensional dynamic model of Kolmogorov flow. Phys. Rev. Fluids 8, 044402.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.
  * Pérez-De-Jesús (2024) Pérez-De-Jesús 2024 Learning data-driven reduced-order models of complex flows. PhD thesis, University of Wisconsin- Madison.
  * Racca & Magri (2022) Racca, A & Magri, L 2022 Data-driven prediction and control of extreme events in a chaotic flow. Phys. Rev. Fluids 7, 104402.
  * Rowley et al. (2009) Rowley, C. W., Mezic, I., Bagheri, S., Schlatter, P & Henningson, D. S. 2009 Spectral analysis of nonlinear flows. J. Fluid Mech. 641, 115–127.
  * Schmid (2010) Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28.
  * Schmid (2022) Schmid, P. J 2022 Dynamic mode decomposition and its variants. Annual Review of Fluid Mechanics 54 (1), 225–254.
  * Smith et al. (2005) Smith, T. R, Moehlis, J. & Holmes, P. 2005 Low-dimensional modelling of turbulence using the proper orthogonal decomposition: a tutorial. Nonlinear Dynamics 41, 275–307.
  * Solera-Rico et al. (2024) Solera-Rico, A., Vila, C. S., Gómez-López, M., Wang, Y., Almashjary, A., Dawson, S. T. M. & Vinuesa, R. 2024 Variational autoencoders and transformers for reduced-order modelling of fluid flows. Nature Communications 15 (1), 1361.
  * Srinivasan et al. (2019) Srinivasan, P A, Guastoni, L, Azizpour, H, Schlatter, P & Vinuesa, R 2019 Predictions of turbulent shear flows using deep neural networks. Phys. Rev. Fluids 4, 054603.
  * Temam (1989) Temam, R. 1989 Do inertial manifolds apply to turbulence? Phys. D: Nonlinear Phenom. 37 (1), 146–152.
  * Van Goor et al. (2023) Van Goor, P., Mahony, R., Schaller, M. & Worthmann, K. 2023 Reprojection methods for koopman-based modelling and prediction. In 2023 62nd IEEE Conference on Decision and Control (CDC), pp. 315–321. IEEE.
  * Waleffe (1997) Waleffe, F. 1997 On a self-sustaining process in shear flows. Phys. Fluids 9 (4), 883–900.
  * Waleffe (1998) Waleffe, F. 1998 Three-dimensional coherent states in plane shear flows. Phys. Rev. Lett. 81, 4140–4143.
  * Waleffe (2001) Waleffe, F. 2001 Exact coherent structures in channel flow. J. Fluid Mech. 435, 93–102.
  * Williams et al. (2014) Williams, M., Kevrekidis, I. & Rowley, C. 2014 A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition. J. Nonlinear Sci. 25.
  * Xi & Graham (2010) Xi, L & Graham, M D 2010 Active and hibernating turbulence in minimal channel flow of Newtonian and polymeric fluids. Phys. Rev. Lett. 104, 218301.
  * Yang & Radons (2012) Yang, H. & Radons, G. 2012 Geometry of inertial manifolds probed via a lyapunov projection method. Phys. Rev. Lett. 108, 154101.
  * Yang et al. (2009) Yang, H., Takeuchi, K. A., Ginelli, F., Chaté, H. & Radons, G. 2009 Hyperbolicity and the effective dimension of spatially extended dissipative systems. Phys. Rev. Lett. 102, 074102.
  * Yeung et al. (2019) Yeung, E., Kundu, S. & Hodas, N. 2019 Learning deep neural network representations for Koopman operators of nonlinear dynamical systems pp. 4832–4839.
  * Young et al. (2023) Young, C. D., Corona, P. T., Datta, A., Helgeson, M. E. & Graham, M. D. 2023 Scattering-Informed Microstructure Prediction during Lagrangian Evolution (SIMPLE)—a data-driven framework for modeling complex fluids in flow. Rheologica Acta 62 (10), 587–604.
  * Zelik (2022) Zelik, S. 2022 Attractors. then and now, arXiv: 2208.12101.
  * Zeng et al. (2024) Zeng, K., Pérez-De-Jesús, C. E, Fox, A. J & Graham, M. D 2024 Autoencoders for discovering manifold dimension and coordinates in data from complex dynamical systems. MLST 5 (2), 025053.