3.1. Description of PCF data

In the following sections we apply DManD to DNS of turbulent PCF in a MFU domain. Specifically, we consider the well-studied Hamilton, Kim and Waleffe (HKW) domain (Hamilton et al. Reference Hamilton, Kim and Waleffe1995). We made this selection to compare our DManD results with the analysis of the SSP in this domain, to compare our DManD results with other Galerkin-based ROMs and to compare our DManD results with known unstable periodic solutions in this domain.

For PCF we solve the Navier–Stokes equations (NSEs)

(3.1) \begin{equation} \dfrac{\partial \boldsymbol{v}}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}=-\boldsymbol{\nabla} p+Re^{-1}\nabla^2 \boldsymbol{v}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0 \end{equation}

for a fluid confined between two plates moving in opposite directions with the same speed. Equation (3.1) is the non-dimensionalised form of the equations with velocities in the streamwise $x\in [0,L_x]$, wall-normal $y\in [-1,1]$ and spanwise $z\in [0,L_z]$ directions defined as $\boldsymbol {v}=[v_x,v_y,v_z]$ and pressure $p$. We solve this equation for a domain with periodic boundary conditions in $x$ and $z$ ($\boldsymbol {v}(0,y, z)=\boldsymbol {v}(L_x,y, z), \boldsymbol {v}(x,y, 0)=\boldsymbol {v}(x,y, L_z)$) and no-slip, no-penetration boundary conditions in $y$ ($v_x(x,\pm 1, z)=\pm 1, v_y(x,\pm 1, z)=v_z(x,\pm 1, z)=0$). The complexity of the flow increases as the Reynolds number $Re$ increases and the domain size $L_x$ and $L_z$ increase. Here we use the HKW cell, which is at $Re=400$ with a domain size $[L_x,L_y,L_z]=[1.75{\rm \pi},2,1.2{\rm \pi} ]$ (Hamilton et al. Reference Hamilton, Kim and Waleffe1995). The HKW cell is one of the simplest flows that sustains turbulence for extended periods of time before relaminarising. We chose to use this flow because it is well studied (refer to § 1), it isolates the SSP (Hamilton et al. Reference Hamilton, Kim and Waleffe1995) and a library of ECS exists for this domain (Gibson et al. Reference Gibson, Halcrow, Cvitanović and Viswanath2008b). Here we are interested in modelling the turbulent dynamics, so we will remove data upon relaminarisation as detailed in the following.

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

(3.2)\begin{equation} \mathcal{P}[(v_x,v_y,v_z,p)(x,y,z,t)]=(-v_x,-v_y,-v_z,p)(-x,-y,-z,t) \end{equation}

reflection about the $z=0$ plane

(3.3)\begin{equation} \mathcal{R}[(v_x,v_y,v_z,p)(x,y,z,t)]=(v_x,v_y,-v_z,p)(x,y,-z,t) \end{equation}

and rotation by ${\rm \pi}$ about the $z$-axis

(3.4)\begin{equation} \mathcal{RP}[(v_x,v_y,v_z,p)(x,y,z,t)]=(-v_x,-v_y,v_z,p)(-x,-y,z,t). \end{equation}

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

(3.5)\begin{equation} \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). \end{equation}

We incorporate all these symmetries in the POD representation (Smith et al. Reference Smith, Moehlis and Holmes2005), as we discuss further in § 3.2. Then, we use the method of slices (Budanur et al. Reference Budanur, Borrero-Echeverry and Cvitanović2015a) to phase align in the $z$ direction. By phase aligning in $z$ we fix the location of the low-speed streak. Without the alignment in $z$, models performed poorly because the models needed to learn how to represent every spatial shift of every snapshot. In what follows, we only consider phase alignment in $z$, but we note that extending this work to phase alignment in $x$ is straightforward. To phase align the data, we use the first Fourier mode method of slices (Budanur et al. Reference Budanur, Borrero-Echeverry and Cvitanović2015a). First, we compute a phase by taking the Fourier transform of the streamwise velocity in $x$ and $z$ ($\hat {v}_x(k_x,y,k_z)=\mathcal {F}_{x,z}(v_x)$) at a specific $y$ location ($y_1$) to compute the phase

(3.6)\begin{equation} \phi=\text{atan2}(\text{imag}(\hat{v}_x(0,y_1,1)),\text{real}(\hat{v}_x(0,y_1,1))). \end{equation}

The variables $k_x$ and $k_z$ are the streamwise and spanwise wavenumbers. We select $y_1$ to be one grid point off the bottom wall. Any $y$ location should work for the phase alignment, but we found choosing the point in the viscous sublayer resulted in a more gradual change of the phase. Rapid changes in the phase make prediction difficult and can require rescaling time, as in Budanur et al. (Reference Budanur, Cvitanović, Davidchack and Siminos2015b), which we want to avoid. Then we compute the pattern dynamics by using the Fourier shift theorem to set the phase to 0 (i.e. move the low-speed streak to the centre of the channel)

(3.7)\begin{equation} \boldsymbol{v}_p=\mathcal{F}_{x,z}^{-1}(\hat{\boldsymbol{v}}\exp(-{\rm i}k_z\phi)). \end{equation}

We generate turbulent PCF trajectories using the pseudo-spectral Channelflow code developed by Gibson (Reference Gibson2012) and Gibson et al. (Reference Gibson2021). In this code, the velocity and pressure fields are projected onto Fourier modes in $x$ and $z$ and Chebyshev polynomials of the first kind in $y$. These coefficients are evolved forward in time first using the multistage SMRK2 scheme (Spalart, Moser & Rogers Reference Spalart, Moser and Rogers1991), then, after taking multiple timesteps, using the multistep Adams–Bashforth backward-differentiation 3 scheme (Peyret Reference Peyret2002). At each timestep, a pressure boundary condition is found such that incompressibility is satisfied at the wall (${\rm d}v_y/{{\rm d}\kern 0.05em y}=0$) using the influence matrix method and tau correction developed by Kleiser & Schumann (Reference Kleiser and Schumann1980).

Data were generated with $\Delta t=0.0325$ on a grid of $[N_x,N_y,N_z]=[32,35,32]$ in $x$, $y$ and $z$ for the HKW cell. Starting from random divergence-free initial conditions, we ran simulations forward for either $10\,000$ time units or until relaminarisation. Then we dropped the first 1000 time units as transient data and the last $1000$ time units to avoid laminar data, and repeated with a new initial condition until we had $91\,562$ time units of data stored at intervals of one time unit. We split these data into $80\,\%$ for training and $20\,\%$ for testing. Finally, we preprocess the data by computing the mean over snapshots and the $x$ and $z$ directions $\bar {\boldsymbol {v}}(y)$ from the training data and subtracting it from all data $\boldsymbol {v}'=\boldsymbol {v}-\bar {\boldsymbol {v}}$, and then we compute the pattern $\boldsymbol {v}'_p$ and the phase $\phi$ as described previously. The pattern $\boldsymbol{u}_p$ described in § 2 is $\boldsymbol {v}'_p$ flattened into a vector (i.e. $d=3N_xN_yN_z$). The POD and NN training use only the training data, and all comparisons use test data unless otherwise specified.

3.2. Dimension reduction and dynamic model construction

3.2.1. Linear dimension reduction with POD: from ${O}(10^5)$ to ${O}(10^3)$

The first task in DManD for this Couette flow data is finding a low-dimensional parameterisation of the manifold on which the long-time dynamics lie. We parameterise this manifold in two steps. First, we reduce the dimension down from $d= {O}(10^5)$ to $d_{{POD}}={O}(10^3)$ with the POD and, second, we use an autoencoder to reduce the dimension down to $d_h$, which ideally is equal to $d_{\mathcal {M}}$. The first step is simply a preprocessing step to reduce the size of the data, which reduces the number of parameters in the autoencoder. Due to Whitney's embedding theorem (Whitney Reference Whitney1936, Reference Whitney1944), we know that as long as $d_{\mathcal {M}}< d_{POD}/2$, then this POD representation is diffeomorphic to the full state. As we show later, the manifold dimension $d_{\mathcal {M}}$ appears to be far lower than $d_{POD}/2$, so no information of the full state should be lost with this first step.

In other words, there are four different spaces for representing the state in this process. There is the (infinite-dimensional) solution space of the NSEs, the $d$-dimensional data space generated by the DNS, which we will refer to as the ‘full state’, the $d_{POD}$-dimensional POD representation of the full state and the $d_h$-dimensional representation in the manifold coordinate system. We assume that our high-resolution DNS solution accurately represents a true solution of the NSEs, and further assume that this solution lies on a low-dimensional manifold of dimension $d_{\mathcal {M}}$. Then, given a sufficient number of dimensions, the POD, followed by the autoencoder, should provide a representation diffeomorphic to the full state. However, we cannot know that a sufficient number of dimensions has been chosen until after training the models.

The POD originates with the question of what function $\boldsymbol {\varPhi }$ maximises

(3.8)\begin{equation} \dfrac{\langle |(\boldsymbol{v}',\boldsymbol{\varPhi})_E|^2\rangle }{\|\boldsymbol{\varPhi}\|_E^2}, \end{equation}

where $\langle {\cdot }\rangle$ is the ensemble average and the inner product is defined to be

(3.9)\begin{equation} (\boldsymbol{q}_1,\boldsymbol{q}_2)_E=\frac{1}{2L_xL_z}\iiint_V \boldsymbol{q}_1 \boldsymbol{\cdot} \boldsymbol{q}_2 \,{\rm d}\boldsymbol{x}, \end{equation}

with the corresponding energy norm $\|\boldsymbol {q}\|_E^2=(\boldsymbol {q},\boldsymbol {q})_E$. Solutions $\boldsymbol {\varPhi }^{(n)}$ to this problem satisfy the eigenvalue problem

(3.10)\begin{equation} \sum_{j=1}^3 \int_0^{L_x}\int_{-1}^1 \int_0^{L_z} \langle v'_i(\boldsymbol{x}, t) v_j^{\prime*}(\boldsymbol{x}^{\prime}, t)\rangle \varPhi_j^{(n)}(\boldsymbol{x}^{\prime})\, {\rm d}\boldsymbol{x}^{\prime}=\lambda_i \varPhi_i^{(n)} (\boldsymbol{x}) \end{equation}

(Smith et al. Reference Smith, Moehlis and Holmes2005; Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012). Unfortunately, upon approximating these integrals, with the trapezoidal rule for example, this becomes a $d \times d$ matrix, making computation intractable. Furthermore, computing the average in (3.10), without any modifications, results in POD modes that fail to preserve the underlying symmetries described previously.

In order to make this problem computationally tractable, and preserve symmetries, we apply the POD method used in Smith et al. (Reference Smith, Moehlis and Holmes2005), with the slight difference that we first subtract off the mean of the state before performing the analysis. The first step in this procedure is to treat the POD modes as Fourier modes in both the $x$ and $z$ directions. Holmes et al. (Reference Holmes, Lumley, Berkooz and Rowley2012) showed that for translation-invariant directions, Fourier modes are the optimal POD modes. This step transforms the eigenvalue problem into

(3.11)\begin{equation} L_x L_z \sum_{j=1}^3 \int_{-1}^1\langle \hat{v}'_i(k_x,y^{\prime},k_z,t) \hat{v}_j^{\prime*}(k_x,y^{\prime},k_z,t)\rangle \varphi_{jk_xk_z}^{(n)}(y^{\prime}) \,{{\rm d}\kern0.05em y}^{\prime}=\lambda_{k_xk_z}^{(n)} \varphi_{ik_xk_z}^{(n)}(y), \end{equation}

which reduces the $d \times d$ eigenvalue problem down to a $3 N_y \times 3 N_y$ eigenvalue problem for every wavenumber pair $(k_x,k_z)$ of Fourier coefficients. We used $5000$ snapshots evenly sampled over the training data to compute the POD modes. Then, to account for the discrete symmetries, the data are augmented such that the mean in (3.11) is computed by adding all the discrete symmetries of each snapshot, i.e. we compute $\langle \hat {v}'_i \hat {v}_j^{\prime *}\rangle$ for $\boldsymbol {v}$, $\mathcal {P}\boldsymbol {v}$, $\mathcal {R}\boldsymbol {v}$ and $\mathcal {RP}\boldsymbol {v}$ and average the results.

This analysis gives us POD modes

(3.12)\begin{equation} \boldsymbol{\varPhi}_{k_xk_z}^{(n)}(\boldsymbol{x})=\frac{1}{\sqrt{L_x L_z}} \exp \left(2 {\rm \pi}{\rm i}\left(\frac{k_x x}{L_x}+\frac{k_z z}{L_z}\right)\right) \boldsymbol{\varphi}_{k_xk_z}^{(n)}(y), \end{equation}

and eigenvalues $\lambda _{k_xk_z}^{(n)}$. The projection onto these modes results in complex values unless $k_x=k_z=0$. We sort the modes from largest eigenvalue to smallest eigenvalue ($\lambda _i$) and project the phase-aligned data (i.e. $\boldsymbol {v}'_p$, as described in § 3.1) onto the leading $256$ modes, giving us a vector of POD coefficients $\boldsymbol{a}(t)$. A majority of these modes are complex (i.e. 2 degrees of freedom), so projecting onto these modes results in a $502$-dimensional system, i.e. $d_{POD}=502$. In figure 1(a) we plot the eigenvalues, which show a rapid drop and then a long tail that contributes little to the energy content. By dividing the eigenvalues of the leading $256$ modes by the total, we find these modes contain $99.8\,\%$ of the energy. To further illustrate that $256$ modes is sufficient to represent the state in this case, we consider the reconstruction of statistics after projecting onto the POD modes. In figure 1(b) we show the reconstruction of four components of the Reynolds stress, $\langle v^{\prime 2}_x\rangle$, $\langle v^{\prime 2}_z\rangle$, $\langle v^{\prime 2}_y\rangle$ and $\langle v^{\prime }_xv^{\prime }_y\rangle$. The projection onto POD modes matches all of these quantities extremely well.

Figure 1. (a) Eigenvalues of POD modes sorted in descending order. (b) Components of the Reynolds stress for data generated by the DNS and this data projected onto 256 POD modes. In (b) the curves are, from top to bottom, $\langle v^{\prime 2}_x\rangle$, $\langle v^{\prime 2}_z\rangle$, $\langle v^{\prime 2}_y\rangle$ and $\langle v^{\prime }_xv^{\prime }_y\rangle$.

3.2.2. Nonlinear dimension reduction with autoencoders: from ${O}(10^3)$ to ${O}(10^1)$

Now that we have converted the data to POD coefficients and filtered out the low-energy modes, we next train an autoencoder to perform nonlinear dimension reduction. However, we first ‘preprocess’ the POD coefficients by normalising the coefficients before using them in the autoencoder. It is common practice before NN training to subtract the mean and divide by the standard deviation of each component. We do not take this approach here because the standard deviation of the higher POD coefficients, which contribute less to the reconstruction, is much smaller than the lower POD coefficients. In order to retain the important information in the magnitudes, but put the values in a more amenable form for NN training, we instead normalise the POD coefficients by subtracting the mean and dividing by the maximum standard deviation.

With the above preprocessing step completed, we now turn to the reduction of dimension with the nonlinear approach enabled by the autoencoder structure. We consider three approaches to reducing the dimension of the normalised POD coefficient vector $\boldsymbol{a}$: (1) training a hybrid autoencoder; (2) training a standard autoencoder; and (3) linear projection onto a small set of POD modes. We describe the first two approaches in § 2, and note that the POD projection onto $256$ (complex) modes can be written as $\boldsymbol{a}=\boldsymbol{\mathsf{U}}_r^{\rm T}\boldsymbol{u}$. The third approach corresponds to setting $\mathcal {E}$ and $\mathcal {D}$ to zero in 2.5 and 2.6.

To optimise the parameters of the hybrid and standard autoencoders, we use an Adam optimiser (Kingma & Ba Reference Kingma and Ba2015) in Keras (Chollet Reference Chollet2015) to minimise the loss in (2.7) with the normalised POD coefficients as inputs. We train for 500 epochs with a learning rate scheduler that drops the learning rate from $10^{-3}$ to $10^{-4}$ after 400 epochs. At this point, we see no improvement in the reconstruction error. For the hybrid autoencoder approach, we set $\alpha =0.01$, and for the standard autoencoder $\alpha =0$ (in practice, this term is not included). Table 1 includes additional NN architecture details. We determined these network parameters through a manual parameter sweep by varying the shape and activation functions of the NNs to achieve the lowest error without excessive computational cost.

Table 1. Architectures of NNs. ‘Shape’ indicates the dimension of each layer, ‘Activation’ the corresponding activation functions and ‘sig’ is the sigmoid activation. ‘Learning rate’ gives the learning rate at multiple times during training. We dropped the learning rate at even intervals.

In figure 2(a) we show the relative reconstruction error for the three approaches over a range of dimensions $d_h$. We compute this error by computing the energy norm of the difference between the state in the ambient space and the predicted state. Then, we divide that quantity by the energy norm of the state in the ambient space and average over the test data. We use NNs with the same architectures for both the standard and the hybrid autoencoder approaches. Due to the variability introduced into autoencoder training by randomly initialised weights and stochasticity in the optimisation, we show the error for four separately trained autoencoders, at each $d_h$. We see that the autoencoders perform much better than POD in the range of dimensions considered here.

Figure 2. (a) Relative error in the energy norm on test data for POD, standard autoencoders and hybrid autoencoders at various dimensions $d_h$. At each dimension there are four standard and four hybrid autoencoders. (b) The centreline streamwise velocity for three examples in the DNS and reconstructed using the hybrid autoencoder with $d_h=18$. The relative error of reconstruction from top to bottom is $0.021$, $0.071$ and $0.048$.

The performance of the standard and hybrid autoencoder approaches is very similar, and the MSE begins to level off for $d_h\gtrsim 18$. Given a perfect autoencoder (i.e. one that has latent dimension $d_h\geq d_{\mathcal {M}}$ and sufficient capacity to determine an exact mapping between the ambient space and the manifold), the error in this figure ideally should vanish once $d_h$ becomes high enough; however, there is always some error in the autoencoder, either due to finite capacity or limitations in the training process, that prevents this from happening. Figure 2(b) shows three examples of the flowfields and their reconstruction for a hybrid autoencoder at $d_h=18$. This autoencoder captures the key structures of these snapshots. The results in the following use the hybrid approach. We took this approach because the low-dimensional representations of the hybrid autoencoders displayed less variability in predictions from trial to trial than those of the standard autoencoders. This led, in turn, to less variability when training neural ODEs to predict the time dynamics of the different low-dimensional representations of these autoencoders.

Before discussing how to use the low-dimensional representations from these autoencoders to train stabilised neural ODEs, we briefly mention convolutional neural networks (CNNs). These could be a reasonable choice for the autoencoder in the framework outlined previously. We chose not to use them because empirically we found lower errors using dense NNs in the POD basis for simple chaotic systems such as the KSE. Once in the POD basis, a CNN holds no advantage over a densely connected network because there are no longer spatial connections between the POD coefficients. In addition, although standard CNNs are invariant to small translations of the input they are not equivariant to inputs, which is what we desire. Separating out pattern and phase variables ensures equivariance for any function approximation of $\boldsymbol{\chi}$ and $\boldsymbol{\check {\chi}}$, making it an appealing choice over trying to constrain the structure of the NN to ensure equivariance.

3.3. Short-time tracking

In the following two sections, we evaluate the performance of the DManD models at reconstructing short-time trajectories and long-time statistics. Figure 3 shows snapshots of the streamwise velocity at the centre plane of the channel, $y=0$, for the DNS and DManD at $d_h=18$. We choose to show results for $d_h=18$ because the autoencoder error begins to level off around this dimension and, as we will show, the error in statistics levels off before this dimension. The value $d_h=18$ is not necessarily the minimal dimension required to model this system. In figure 3, both the DNS and the DManD model show key characteristics of the SSP: (1) low-speed streaks become wavy; (2) the wavy low-speed streaks break down generating rolls; (3) the rolls lift fluid from the walls, regenerating streaks.

Figure 3. Snapshots of the streamwise velocity at $y=0$ from the DNS and from the DManD model at $d_h=18$.

Not only does DManD capture the qualitative behaviour of the SSP, but figure 3 also shows good quantitative agreement as well. To further illustrate this, in figure 4 we show the modal root-mean-squared (RMS) velocity

(3.13) \begin{equation} M(k_x,k_z)=\left(\int_{-1}^1 (\hat{v}^2_x(k_x,y,k_z)+\hat{v}^2_y(k_x,y,k_z)+\hat{v}^2_z(k_x,y,k_z))\,{{\rm d}\kern0.05em y}\right)^{1/2}, \end{equation}

which Hamilton et al. (Reference Hamilton, Kim and Waleffe1995) used to identify the different parts of the SSP. Specifically, we consider the $M(0,1)$ mode, which corresponds to the low-speed streak, and the $M(1,0)$ mode, which corresponds to the $x$-dependence that appears when the streak becomes wavy and breaks up. In this example, the two curves match well over a cycle of the SSP and only start to move away after $\sim$150 time units, which is about three Lyapunov times.

Figure 4. Modal RMS velocity from the DNS ($M$) and from the DManD model at $d_h=18$ ($\tilde {M}$). The markers correspond to the times in figure 3.

While the previous result shows a single example, we also consider ensembles of initial conditions. Figure 5 shows the tracking error $\|\boldsymbol{a}(t_i+t)-\boldsymbol{\tilde{a}}(t_i+t)\|$ of 10 trajectories, starting at $t_i$, for a model at $d_h=18$. Here we normalise the tracking error by the error between solutions at random times $t_i$ and $t_j$:$\mathcal {N}=\langle \|\boldsymbol{a}(t_i)-\boldsymbol{a}(t_j)\|\rangle$. Here we note that, because the 256 complex modes we use capture $99.8\,\%$ of the energy in the flow, this normalised tracking error will be approximately equal to the normalised tracking error in the full space. In this case, the darkest line corresponds to the flow field in figures 3 and 4. When considering the other initial conditions in figure 5, there tends to be a relatively slow rise in the error over $\sim$50 time units and then a more rapid increase after this point. To better understand how this tracking varies with the dimension of the model we next consider the ensemble-averaged tracking error.

Figure 5. Normalised tracking error for 10 random initial conditions (different shades) using DManD with $d_h=18$.

In figure 6(a) we show the normalised ensemble-averaged tracking error for model dimensions between $d_h=3$ and $18$. For $d_h=3\unicode{x2013}5$ there is a rapid rise in the error until $\sim$40 time units after which the error levels off. This behaviour often happens due to trajectories quickly diverging and landing on stable fixed points or POs that do not exist in the true system. For $d_h=6\unicode{x2013}10$ there is an intermediate behaviour where lines diverge more quickly than on the higher-dimensional models, but tend to approach the same tracking error at $\sim$100 time units. Then, for the remaining models $d_h=11\unicode{x2013}18$, there is a smooth improvement in the tracking error over this time interval. As the dimension increases in this range the trends stay the same, but the error continues to decrease, which is partially due to improvement in the autoencoder performance.

Figure 6. (a) Ensemble-averaged short-time tracking and (b) temporal autocorrelation of the kinetic energy for DManD models of increasing dimension. In (b), odd numbers above $d_h=5$ are omitted for clarity.

The instantaneous kinetic energy of the flow is given by $E(t)=1/2\|\boldsymbol {v}\|_E^2$ (3.9) and we denote its fluctuating part as $k(t)=E(t)-\langle E\rangle$. In figure 6(b) we show the temporal autocorrelation of $k$. Again, for $d_h=3\unicode{x2013}5$ we see clear disagreement between the true autocorrelation and the prediction. Above $d_h>5$ all of the models match the temporal autocorrelation well, without a clear trend in the error as dimension changes. All these models match well for $\sim$40 time units, with $d_h=18$ (the darkest line) matching the data extremely well for two Lyapunov times.

Finally, before investigating the long-time predictive capabilities of the model, we show the tracking of phase dynamics for $d_h=18$. As mentioned in § 2, we decouple the phase and pattern dynamics such that the time evolution of the phase only depends upon the pattern dynamics. Here we take the $d_h=18$ model and used it to train an ODE for the phase dynamics. For training we repeat the process used for training $\boldsymbol{g}_{NN}$ to train $\boldsymbol{g}_\phi$ with the loss in (2.14). Table 1 contains details on the NN architecture.

In figure 7(a) we show an example of the model phase evolution over 200 time units. In this example, the model follows the same downward drift in phase despite not matching exactly. Then, to show the statistical agreement between the DNS and the model, we show the mean squared phase displacement $\text {MSD}=\langle (\phi (t_i)-\phi (t_i+t))^2\rangle$ for both the DNS and the model in figure 7(b), as was done for Kolmogorov flow by Pérez De Jesús & Graham (Reference Pérez De Jesús and Graham2023). The curves are in good agreement. All of the remaining long-time statistics we report are phase invariant, so the remaining results use only models for the pattern dynamics.

Figure 7. (a) Example of the phase evolution and (b) the MSD of the phase evolution for the DNS and the DManD model at $d_h=18$.

3.4. Long-time statistics

Next, we investigate the ability of the DManD model to capture the long-time dynamics of PCF. An obvious prerequisite for models to capture long-time dynamics is the long-time stability of the models. As mentioned in § 2, the long-time trajectories of standard neural ODEs often become unstable, which led us to use stabilised neural ODEs with an explicit damping term. We quantify this observation by counting, of the $16$ models trained at each dimension $d_h$, how many become unstable with and without the presence of an explicit damping term. From our training data, we know where $\boldsymbol{h}$ should lie, so if it falls far outside this range after some time we can classify the model as unstable. In particular, we classify models as unstable if the norm of the final state is two times that of the maximum in our data ($\|\boldsymbol{\tilde{h}}(T)\|>2\max _t \|\boldsymbol{h}(t)\|$), after $T=10^4$ time units. In all of the unstable cases $\|\boldsymbol{\tilde {h}}(t)\|$ follows the data over some short time range before eventually growing indefinitely.

In figure 8 we show the number of unstable models with and without damping. With damping, all of the models are stable, whereas without damping almost all models become unstable for $d_h=5\unicode{x2013}16$, and around half become unstable in the other cases. In addition, with longer runs or with different initial conditions, many of the models without damping labelled as stable here also eventually become unstable. This lack of stability happens when inaccuracies in the neural ODE model push trajectories off the attractor. Once off the attractor, the model is presented with states unlike the training data leading to further growth in this error. We show more results highlighting this behaviour in Linot & Graham (Reference Linot and Graham2022) and Linot et al. (Reference Linot, Burby, Tang, Balaprakash, Graham and Maulik2023a). Thus, although some standard neural ODE models do provide reasonable statistics, using these models presents challenges due to this lack of robustness. As such, all other results we show use stabilised neural ODEs.

Figure 8. Fraction of unstable DManD models with standard neural ODEs and with stabilised neural ODEs at various dimensions.

While figure 8 indicates that stabilised neural ODEs predict $\boldsymbol{\tilde {h}}$ in a similar range to that of the data, it does not quantify the accuracy of these predictions. In fact, with few dimensions many of these models do not remain chaotic, landing on fixed points or POs. The first metric we use to quantify the long-time performance of the DManD method is the mean-squared POD coefficient amplitudes ($\langle \|a_n\|^2\rangle$). We consider this quantity because Gibson reports it for POD-Galerkin in Gibson (Reference Gibson2002) at various levels of truncation. In figure 9 we show how well the DManD model, with $d_h=18$, captures this quantity, in comparison to the POD-Galerkin model in Gibson (Reference Gibson2002). The two data sets slightly differ because we subtract the mean before applying POD and Gibson did not. The DManD method, with only 18 degrees of freedom, matches the mean-squared amplitudes to high accuracy, far better than all of the POD-Galerkin models. It is not until POD-Galerkin keeps 1024 modes that the results become comparable, which corresponds to $\sim$2000 degrees of freedom because most coefficients are complex. In addition, our method requires only data, whereas the POD-Galerkin approach requires both data for computing the POD and knowledge of the equations of motion for projecting the equations onto these modes.

Figure 9. Comparison of $\langle \|a_n\|^2\rangle$ (mean-squared POD coefficient amplitudes) from the DNS to (a) $\langle \|a_n\|^2\rangle$ from the DManD model at $d_h=18$ and (b) $\langle \|a_n\|^2\rangle$ from POD-Galerkin with $N$ POD modes (reproduced with permission from Gibson Reference Gibson2002). In (b), the quantity $\lambda$ is equivalent to $\langle \|a_n\|^2\rangle$ from the DNS.

We now investigate how the Reynolds stress and the power input versus dissipation vary with dimension. Figure 10 shows four components of the Reynolds stress at various dimensions. For $\langle v_x^{\prime 2}\rangle$ and $\langle v'_xv'_y\rangle$, nearly all the models match the data, with relatively small deviations only appearing for $d_h\sim 3\unicode{x2013}6$. For $\langle v_y^{\prime 2}\rangle$ and $\langle v_z^{\prime 2}\rangle$, this deviation becomes more obvious, and the lines do not converge until around $d_h > 10$, with all models above this dimension exhibiting a minor overprediction in $\langle v_z^{\prime 2}\rangle$.

Figure 10. Components of the Reynolds stress with increasing dimension for DManD models at various dimensions. Odd numbers above $d_h=5$ are omitted for clarity.

To evaluate how accurate the models are at reconstructing the energy balance, we look at joint PDFs of power input and dissipation. The power input is the amount of energy required to move the walls:

(3.14) \begin{equation} I=\dfrac{1}{2L_xL_z}\int_0^{L_x}\int_0^{L_z} \left.\dfrac{\partial v_x}{\partial y}\right|_{y=-1}+\left.\dfrac{\partial v_x}{\partial y}\right|_{y=1}\,{\rm d}x\,{\rm d}z, \end{equation}

and the dissipation is the energy lost to heat due to viscosity:

(3.15)\begin{equation} D=\dfrac{1}{2L_xL_z}\int_0^{L_x}\int_{-1}^1\int_0^{L_z}| \boldsymbol{\nabla} \times \boldsymbol{v}|^2 \,{\rm d}x\,{\rm d}y\,{\rm d}z. \end{equation}

These two terms define the rate of change of energy in the system $\dot {E}=(I-D)/Re$ (Kaszás, Cenedese & Haller Reference Kaszás, Cenedese and Haller2022), which must average to zero over long times. Checking this statistic is important to show the DManD models correctly balance the energy.

Figure 11(a–c) show the PDF from the DNS, the PDF for $d_h=6$ and the PDF for $d_h=18$, generated from a single trajectory evolved for $5000$ time units, and figures 11(e,f) show the absolute difference between the true and model PDFs. With $d_h=6$ the model overestimates the number of low dissipation states, while $d_h=18$ matches the density well. In figure 11(d) we compare the joint PDFs at all dimensions with the true PDF using the earth movers distance (EMD) (Rubner, Tomasi & Guibas Reference Rubner, Tomasi and Guibas1998). The EMD determines the distance between two PDFs as a solution to the transportation problem by treating the true PDF as ‘supplies’ and the model PDF as ‘demands’ and finding the ‘flow’ that minimises the work required to move one to the other. Specifically, we find the flow $f_{ij}$ that minimises $\sum _{i=1}^m\sum _{j=1}^n f_{ij}d_{ij}$ subject to the constraints:

(3.16)$$\begin{gather} f_{i j} \geq 0, \quad 1 \leq i \leq m,\ 1 \leq j \leq n, \end{gather}$$

(3.17)$$\begin{gather}\sum_{j=1}^n f_{i j} = p_i, \quad 1 \leq i \leq m, \end{gather}$$

(3.18)$$\begin{gather}\sum_{i=1}^m f_{i j} = q_j, \quad 1 \leq j \leq n, \end{gather}$$

where $p_i$ is the probability density at the $i$th bin in the model PDF and $q_j$ is the probability density at the $j$th bin in the true PDF, for PDFs with $n$ and $m$ bins (in this case $n=m$). In addition, $d_{ij}$ is the cost to move between bins, which we take to be the $L_2$ distance between bins $i$ and $j$ (for $i=j$ $d_{ij}=0$). After solving this minimisation problem for the optimal flow $f^*_{i j}$ the EMD is defined as

(3.19) \begin{equation} \text{EMD}=\frac{\displaystyle\sum_{i=1}^m \displaystyle\sum_{j=1}^n f_{i j}^* d_{i j}}{\displaystyle\sum_{i=1}^m \displaystyle\sum_{j=1}^n f_{i j}^*}. \end{equation}

For a more detailed explanation of the EMD we refer the reader to Levina & Bickel (Reference Levina and Bickel2001).

Figure 11. (a–c) Examples of joint PDFs for the true system, the DManD model at $d_h=6$ and the DManD model at $d_h=18$. (d) EMD between the PDF from the DNS and the PDFs predicted by the DManD model at various dimensions. ‘DNS’ is the error between two PDFs generated from DNS trajectories of the same length with different initial conditions. (e,f) The error associated with the DManD model PDFs at $d_h=6$ and $d_h=18$.

We compute the distance between PDFs using the EMD because it is a cross-bin distance, meaning the distance accounts for the density in neighbouring bins. This is in contrast to bin-to-bin distances, such as the Kullback–Leibler divergence, which only uses the error at a given bin. Bin-to-bin distances can vary significantly with small shifts in one PDF (misalignment) and when changing the number of bins used to generate the PDF (Ling & Okada Reference Ling and Okada2007). We choose the EMD because it does not suffer from these issues. In figure 11(d) we see a steep drop in the EMD at $d_h=5$ and after $d_h\gtrsim 15$ the joint PDFs are in excellent agreement with the DNS. The dashed line corresponds to the EMD between two different trajectories from the DNS. An interesting observation is that the DManD models appear to accurately capture the statistics of the DNS for $d_h\gtrsim 15$, and Inubushi et al. (Reference Inubushi, Takehiro and Yamada2015) reported a Lyapunov dimension of 14.8. The close agreement between this fractal dimension and the dimension at which our models perform well further supports that, in theory, for the models with $d_h\gtrsim 15$ we have a sufficient number of dimensions to ‘exactly’ parameterise the manifold on which the dynamics lie.

Finally, we investigate the Lyapunov exponents of the DManD models. In figure 12, we report the spectrum of Lyapunov exponents $\lambda ^{LE}_n$ for the DManD model with $d_h=18$, as well as the leading values computed by Inubushi et al. (Reference Inubushi, Takehiro and Yamada2015) from DNS. Figure 12(a) shows the prediction of the Lyapunov exponents through time in the DManD model, averaged over 10 initial conditions using the methods described in Sandri (Reference Sandri1996) with the code made available by Rozdeba (Reference Rozdeba2017). At around 2000 time units we see that the Lyapunov exponents have nearly converged, with relatively small discrepancies across trials, as seen by the standard deviation. In figure 12(b), we report the Lyapunov spectra for the DManD model and the DNS. The DManD model correctly predicts that there are four positive exponents, and the computed values are in very good agreement with those from the DNS. Along with the positive exponents, Inubushi et al. (Reference Inubushi, Takehiro and Yamada2015) identified that there were three zero exponents associated with the time translation symmetry and the spatial translational symmetries in $x$ and $z$. We see that the DManD model reasonably predicts two of these zero Lyapunov exponents, but lacks the third due to the phase alignment of the model. The negative values are not in such good agreement, probably for one or both of two reasons. First, the DManD model considers only the dynamics on the low-dimensional invariant manifold of the long-time dynamics, so cannot capture Lyapunov exponents corresponding to directions off of this manifold. Second, the stabilisation term included in the neural ODE representation may lead to some additional damping. Nevertheless, the quantitative prediction of both the number and magnitude of the unstable and neutrally stable Lyapunov exponents indicates the fidelity of the DManD model with regard to the chaotic dynamics of the flow.

Figure 12. (a) The mean (black) and standard deviation (blue) of the Lyapunov exponents computed with the DManD model through time for 10 random initial conditions. The grey dashed line identifies $\lambda =0$. (b) Lyapunov exponents computed from the model and from the DNS (Inubushi et al. Reference Inubushi, Takehiro and Yamada2015).

As with the other statistics, we may again investigate the exponents as we vary the DManD model dimension. Figure 13 shows the four leading Lyapunov exponents (all positive in the DNS) for the models as we vary the dimension. Again, we compute these Lyapunov exponents from 10 initial conditions evolved forward 2000 time units. The positive Lyapunov exponents may be used to compute the Kolmogorov–Sinai entropy $\text {KS}=\varSigma _{\lambda ^{{LE}}_n>0} \lambda ^{{LE}}_n$. This is a measure of the information created by a dynamical system (Eckmann & Ruelle Reference Eckmann and Ruelle1985). As we increase the model dimension, the estimate of these Lyapunov exponents improves in two different ways. First, from $d_h=4\unicode{x2013}8$ the number of positive Lyapunov exponents increases from one to three, recovering the correct number of four positive Lyapunov exponents at $d_h=9$. Second, from $d_h=9\unicode{x2013}15$ we see a gradual decrease in the two highest Lyapunov exponents and a gradual increase in the other two positive Lyapunov exponents. After $d_h\gtrsim 15$ the accuracy of the Lyapunov exponents only shows mild improvement, further supporting the conclusion that near $d_h=15$ we have a sufficient number of dimensions. At $d_h=18$ (the most accurate model) the Kolmogorov–Sinai entropy of the DManD model is $\text {KS}=0.054$ which is close to the DNS value of $\text {KS}=0.048$ (Inubushi et al. Reference Inubushi, Takehiro and Yamada2015).

Figure 13. The four leading Lyapunov exponents for the DManD models at various dimensions and for the DNS.

3.5. Finding ECS in the model

Now that we know that the DManD model quantitatively captures many of the key characteristics of MFU PCF, we now want to explore using the model to discover ECS. In particular, we first investigate whether known POs of the DNS exist in the DManD model, and then we use the DManD model to search for new POs. Finding good initial conditions for use in a PO search is a challenging task, for example, Page & Kerswell (Reference Page and Kerswell2020) compared using DManD with recurrent flow analysis for this task. Here we show that the DManD model offers a rapid method for finding useful initial conditions to input into the full DNS ECS solver. As our model predicts phase-aligned dynamics, the POs of the DManD model are either POs or RPOs, depending on the phase evolution, which we have not tracked. In the following, we omit all $\tilde {{\cdot }}$ for clarity, so all functions should be assumed to come from a DManD model.

Here we outline the approach we take to find POs, which follows Cvitanović et al. (Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2016). When searching for POs we seek an initial condition to a trajectory that repeats after some time period. This is equivalent to finding the zeros of

(3.20) \begin{equation} \boldsymbol{H}(\boldsymbol{h},T)=\boldsymbol{G}_T(\boldsymbol{h})-\boldsymbol{h}, \end{equation}

where $\boldsymbol{G}_T(\boldsymbol{h})$ is the flow map forward $T$ time units from $\boldsymbol{h}$, i.e. $\boldsymbol{G}_T(\boldsymbol{h}(t))=\boldsymbol{h}(t+T)$. We compute $\boldsymbol{G}_T(\boldsymbol{h})$ from (2.9). Finding zeros to (3.20) requires that we find both a point $\boldsymbol{h}^*$ on the PO and a period $T^*$ such that $\boldsymbol{H}(\boldsymbol{h}^*,T^*)=0$. One way to find $\boldsymbol{h}^*$ and $T^*$ is by using the Newton–Raphson method.

By performing a Taylor series expansion of $\boldsymbol{H}$ we find near the fixed point $\boldsymbol{h}^*, T^*$ of $\boldsymbol{H}$ that

(3.21) \begin{equation} \left.\begin{gathered} \boldsymbol{H}(\boldsymbol{h}^*, T^*)-\boldsymbol{H}(\boldsymbol{h}, T) \approx \boldsymbol{\mathsf{D}}_{h} \boldsymbol{H}(\boldsymbol{h}, T) \Delta \boldsymbol{h}+\boldsymbol{\mathsf{D}}_T \boldsymbol{H}(\boldsymbol{h}, T) \Delta T, \\ -\boldsymbol{H}(\boldsymbol{h}, T) \approx \boldsymbol{\mathsf{D}}_{h} \boldsymbol{H}(\boldsymbol{h}, T) \Delta \boldsymbol{h}+\boldsymbol{g}(\boldsymbol{G}_T(\boldsymbol{h})) \Delta T, \end{gathered}\right\} \end{equation}

where $\boldsymbol{\mathsf{D}}_{h}$ is the Jacobian of $\boldsymbol{H}$ with respect to $\boldsymbol{h}$, $\boldsymbol{\mathsf{D}}_T$ is the Jacobian of $\boldsymbol{H}$ with respect to the period $T$, $\Delta \boldsymbol{h}=\boldsymbol{h}^*-\boldsymbol{h}$ and $\Delta T=T^*-T$. To have a complete set of equations for $\Delta \boldsymbol{h}$ and $\Delta T$, we supplement (3.21) with the constraint that the updates $\Delta \boldsymbol{h}$ are orthogonal to the vector field at $\boldsymbol{h}$: i.e.

(3.22)\begin{equation} \boldsymbol{g}(\boldsymbol{h})^{\rm T} \Delta \boldsymbol{h}=0. \end{equation}

With this constraint, at Newton step $(i)$, the system of equations becomes

(3.23) \begin{equation} \left[\begin{array}{@{}cc@{}} \boldsymbol{\mathsf{D}}_{h^{(i)}} \boldsymbol{H}(\boldsymbol{h}^{(i)}, T^{(i)}) & \boldsymbol{g}(\boldsymbol{G}_{T^{(i)}}(\boldsymbol{h}^{(i)})) \\ \boldsymbol{g}(\boldsymbol{h}^{(i)})^{\rm T} & 0 \end{array}\right] \left[\begin{array}{@{}c@{}} \Delta \boldsymbol{h}^{(i)} \\ \Delta T^{(i)} \end{array}\right]=-\left[\begin{array}{@{}c@{}} \boldsymbol{H}(\boldsymbol{h}^{(i)},T^{(i)}) \\ 0 \end{array}\right], \end{equation}

which, in the standard Newton–Raphson method, is used to update the guesses $\boldsymbol{h}^{(i+1)}=\boldsymbol{h}^{(i)}+\Delta \boldsymbol{h}^{(i)}$ and $T^{(i+1)}=T^{(i)}+\Delta T^{(i)}$.

Typically, a Newton–Krylov method is used to avoid explicitly constructing the Jacobian (Viswanath Reference Viswanath2007). However, with DManD, computing the Jacobian is simple, fast and requires little memory because the state representation is dramatically smaller in the DManD model than in the DNS. We compute the Jacobian $\boldsymbol{\mathsf{D}}_{h} \boldsymbol{H}(\boldsymbol{h}, T)$ directly, with the same automatic differentiation tools used for training the neural ODE. Furthermore, if we had chosen to represent the dynamics in discrete, rather than continuous time, computation of general POs would not be possible, as the period $T$ can take on arbitrary values and a discrete-time representation would limit $T$ to multiples of the time step. When finding POs of the DManD model we used the SciPy ‘hybr’ method, which uses a modification of the Powell hybrid method (Virtanen et al. Reference Virtanen2020), and for finding POs of the DNS we used the Newton GMRES-Hookstep method built into Channelflow (Gibson et al. Reference Gibson2021). In the following trials we only consider DManD models with $d_h=18$.

For the HKW cell there exists a library of POs made available by Gibson et al. (Reference Gibson, Halcrow, Cvitanović and Viswanath2008b). To investigate whether the DManD model finds POs similar to existing solutions, we took states from the known POs, encoded them and used this as an initial condition in the DManD Newton solver to find POs in the model. In figure 14 we show projections of 12 known POs, which we identify by the period $T$, and compare them with projections of POs found using the DManD model. This makes up a majority of the POs made available by Gibson et al. (Reference Gibson, Halcrow, Cvitanović and Viswanath2008b). Of the other known solutions, three are RPOs with phase shifts in the streamwise direction that our model, with the current set-up, cannot capture. The other two have short periods of $T=19.02$ and $T=19.06$. Most of the POs found with DManD land on initial conditions near that of the DNS and follow similar trajectories to the DNS.

Figure 14. Power input versus dissipation of known POs (period reported in bottom right) from the DNS and POs found in the DManD model at $d_h=18$. The blue line is a long trajectory of the DNS for comparison.

How close many of these trajectories are to the true PO is surprising and encouraging for many reasons. First, the data used for training the DManD model does not explicitly contain any POs. Second, this approach by no means guarantees convergence on a PO in the DManD model. Third, starting with an initial condition from a PO does not necessarily mean that the solution the Newton solver lands on will be the closest PO to that initial condition, so there may exist POs in the DManD model closer to the DNS solutions than what we present here.

Now that we know the DManD model can find POs similar to those known to exist for the DNS, we now use it to search for new POs. First, we searched for POs in three of the $d_h=18$ models by randomly selecting 20 initial conditions and selecting 4 different periods $T=[20,40,60,80]$. We then took the initial conditions and periods for converged POs and decoded and upsampled them onto a $48 \times 49 \times 48$ grid. We performed this upsampling because Viswanath (Reference Viswanath2007) reported that solutions on the coarser grid can be computational artifacts. Finally, we put these new initial conditions into Channelflow and ran another Newton search for 100 iterations. This procedure resulted in us finding nine new RPOs and three existing POs, the details of which we include in table 2.

Table 2. Details on the RPOs and POs found using initial conditions from the DManD model. The first nine solutions are new and the last three had previously been found. ‘Label’ indicates the label in figure 15, $\sigma _z$ corresponds to the phase-shift in $z$, $T$ is the period of the orbit and ‘Error’ is $\|\textrm {shifted final state} - \textrm {initial state}\|/\|\textrm {initial state}\|$, which is the same error as in Viswanath (Reference Viswanath2007).

In figure 15 we show the new RPOs in the DManD model and the RPOs they converged to after putting them into the Channelflow Newton solver as initial guesses. Again, many of the RPOs end up following a similar path through this state space, with the biggest exceptions being RPO1 and RPO6, which converged to low-power input solutions. It is worth noting that this worked well, considering that the DManD initial conditions are POD coefficients from a model trained using data on a coarser grid than used to search for these solutions and considering Newton–Raphson methods can converge to solutions far from the initial condition (as in RPO1 and RPO6). We have described a new method to rapidly find new ECS, wherein an accurate low-dimensional model, such as the DManD model presented here, is used to quickly perform a large number of ECS searches in the model, and then these solutions can be fine-tuned in the full simulation to land on new solutions.

Figure 15. Power input versus dissipation of new POs found in the DManD model at $d_h=18$ and used to find solutions in the DNS (periods reported in table 2). The blue line is a long trajectory of the DNS for comparison.