{'fraction_non_alphanumeric': 0.0280091241843934, 'fraction_numerical': 0.00802786766396124, 'mean_word_length': 4.736863461148306, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CD-ROM: Complementary Deep-Reduced Order Model

Model order reduction through the POD-Galerkin method can lead to dramatic gains in terms of computational efficiency in solving physical problems. However, the applicability of the method to non linear high-dimensional dynamical systems such as the Navier-Stokes equations has been shown to be limited, producing inaccurate and sometimes unstable models. This paper proposes a closure modeling approach for classical POD-Galerkin reduced order models (ROM). We use multi layer perceptrons (MLP) to learn a continuous in time closure model through the recently proposed Neural ODE method. Inspired by Taken’s theorem as well as the Mori-Zwanzig formalism, we augment ROMs with a delay differential equation architecture to model non-Markovian effects in reduced models. The proposed model, called CD-ROM (Complementary Deep Reduced Order Model) is able to retain information from past states of the system and use it to correct the imperfect reduced dynamics. The model can be integrated in time as a system of ordinary differential equations using any classical time marching scheme. We demonstrate the ability of our CD-ROM approach to improve the accuracy of POD-Galerkin models on two CFD examples, even in configurations unseen during training.


Introduction
The simulation of complex physical processes requires solving high dimensional partial differential equations(PDE), often at great numerical cost. This constitutes a major limitation on the applicability of accurate simulation methods to engineering problems, such as model predictive control or iterative design optimisation. Despite this complexity, physical problems can often be suitably approximated by a reduced number of dominant structures, suggesting the existence of low-dimensional manifolds supporting the underlying dynamics and fueling the development of reduction methods.
The construction of representative reduced order models often relies on the identification of a suitable modal basis to project a description of the system onto. These modes are most often extracted from data, through different methods. For example, the well-known Proper Orthogonal Decomposition (POD)-Galerkin method ( [2]) has been applied to the reduction of linear systems with great results ( [38]). It was also applied to more complex problems such as the Navier Stokes equations with varying degrees of success ( [13,7]). A three dimensional Galerkin model was shown in [14] to capture the transient and oscillation regimes of the cylinder flow, provided the POD basis included a mode tailored to support the transition. Models based on Galerkin POD formulations have the advantage of reducing computational cost while preserving a high level of interpretability in terms of coherent structures and involved dynamics. Nevertheless, despite these promising results, the predictive capabilities of Galerkin models were shown to be limited and often lead to inaccurate reduced models. Indeed, while the POD is optimal for the reconstruction error in a suitable norm, the Galerkin projection might not be ( [38]).
Because of these limitations, some research efforts have been focused on deriving purely data-driven models via sparse regression methods ( [31,42]), cluster modeling ( [26]) or Koopman theory ( [23]) to cite just a few. The representation of the dynamics with the Koopman operator has garnered a lot of interest in recent years, both through the DMD algorithm ( [17]) and its extensions (e.g., EDMD and KDMD [30,29]) which leverage nonlinear state transformations to simplify the modeling problem.
These purely data-driven works are part of a wider trend towards hybridisation of simulation methods with machine learning techniques. Specifically, neural networks which have been applied to numerous physical modeling problems, from black box resolutions through graphs networks allowing the resolution of problems on ustructured meshes ( [62]), to more interpretable approaches like the physics informed neural networks proposed in [50]. During the last few years, several works proposed un-supervised learning techniques to avoid the mesh construction (mesh-free methods). These algorithms train deep neural networks to approximate PDE solutions by satisfying the differential operator, initial conditions and boundary conditions for a specific PDE. [39] used fully connected layers to approximate the solution on complex geometries. [44,45] discussed the possibility to solve high dimensional problems through neural networks. The authors of [45] proved that the neural network converges to the PDE as the number of hidden units increases. A more detailed review on machine learning/CFD hybridization can be found in [51]. At the same time, the deep learning community has developed powerful tools by adapting physical modeling concepts to deep learning, for instance the Neural ODE approach was proposed for the modeling of continuous transformations ( [40]), stable neural architectures were also developed by embedding neural networks with invariant structures ( [34,48,46]). Neural networks have also been applied in the context of model order reduction, notably to learn useful representations of a system's state ( [35,37,53,49]).
As part of this interest in hybridization, we propose a novel architecture for the development of closures for Galerkin reduced models. Despite their limited performance with nonlinear systems, the physical information embedded in their structure can be retained and used as a basis for the construction of a hybrid reduction approach. Applications of this idea have already been proposed in the literature ( [54]) and we refer the reader to [56] for a thorough review on the topic. Building on this existing body of work, we derive information about the structure of ROM closures. Specifically, theoretical results such as Taken's theorem [4] and the Mori-Zwanzig formalism [12] are used to motivate a non-Markovian closure term in reduced order models. Gathering information from the past helps to recover unobserved dynamics. This is a standard procedure among system identification techniques such as auto-regressive models ( [22,61]) or HAVOK ( [33]). Similarly, we propose a non-Markovian closure architecture, using neural networks to correct Galerkin models in a time-continuous manner. The training is carried out using a modified version of Neural ODEs [40], leveraging the adaptive checkpoint adjoint method [55] which is better adapted to complex dynamical systems.
The outline of the paper is as follows. The POD-Galerkin reduction method is introduced in Section 2 and its limitations are discussed. Section 3 details our main contribution on the derivation of a continuous in time correction architecture for reduced order models. Motivation for our work through comparisons with existing approaches is provided in Section 4, while the selected test cases and results are respectively discussed in Section 5 and 6. Section 7 concludes the paper.

Model reduction approach
In this Section we aim at describing the POD-Galerkin model reduction strategy, here illustrated in the context of fluid dynamics. In particular we focus on the incompressible Navier-Stokes equations owing to their relevance in engineering applications and since they represent a challenging test-bed for model reduction techniques. The dynamics of the velocity field u(x, t) and pressure field p(x, t) are governed by: with Re the Reynolds number. Computationally demanding operations are required to numerically solve Eq. (1). The typical approach of POD Galerkin Reduced Order Models (ROMs) is then employed to identify a low-dimensional approximation of the Navier-Stokes equations based on time invariant spatial features. The method, and its derivation for the Navier-Stokes equations, is presented in Section 2.1, and its limitations are discussed in Section 2.2. For a deeper introduction to POD Galerkin reduced order models we refer the interested reader to [20,27].

POD Galerkin models
After gathering data from experiments or the numerical simulation of Eq. (1), one can use the Proper Orthogonal Decomposition, also known as Principal Component Analysis (PCA), to identify dominant structures in the data. The POD method relies on a matrix, referred to as the snapshot matrix, U = {u(t i )|i = 1, . . . , n t } whose columns u(t) ∈ R dnx correspond to the d-dimensional state of the system (for example, the three components of the fluid velocity) indexed upon the n x points of a spatial grid at time t. In the sequel, it is supposed that n x ≥ n t .
Let us consider the thin singular value decomposition (SVD) of this real-valued matrix: U = VΣW T , where the matrix V ∈ R dnx×nt and W ∈ R nt×nt respectively hold the left and right singular vectors of U. The diagonal matrix Σ ∈ R nt×nt holds the singular values of U, usually arranged such that σ 1 ≥ σ 2 ≥ . . . ≥ σ nt ≥ 0.
Columns of V are the time-invariant spatial modes defining a POD basis. These POD modes are orthonormal v i , v j = δ ij , with δ ij the Kronecker delta and ·, · is here Euclidean: These modes are useful for dimensionality reduction because, for any r < n t , the subspace spanned by the basis V r = {v 1 , v 2 , . . . , v r } optimally approximates the data U over the set of dn x × r matrices in the sense that it minimizes the reconstruction error E r defined as: This error can be related to the sum of the discarded singular values: E 2 r = nt k=r+1 σ 2 k . This implies that the information captured by the first r modes in the POD basis can be quantified by looking at the following ratio: ( Introducing a projector P = V r V T r of R dnx onto the invariant subspace spanned by the basis V r , and Q = (I − P) its orthogonal complement, the full state vector u(t) can be decomposed as: The key idea is to select a suitable r in Eq. (3) based on an energy criterion such that Qu(t) can be neglected for the problem under consideration. The full state u(t) can then be approximated as: where a(t) := V T r u(t) is the vector containing temporal coefficients associated with the spatial modes in V r . With the approximation defined in Eq. (5) holds, model reduction arises whenever r n x . Under this modal decomposition ansatz, the original spatio-temporal problem of solving the PDE is reduced to the computation of the temporal coefficients a(t) as the spatial information is expressed in the invariant POD basis V r . However, to achieve a gain in computational cost, one needs a way to compute the reduced coordinates vector without accessing the full order solution. This is where the Galerkin projection method can be used to obtain an equation for the reduced coordinates. The method is based on the reduction of the original model equation: where f corresponds to the full order dynamics discretised in space. By projecting this equation onto the timeindependent POD basis, a dynamical model for the reduced coordinates is obtained: With this Galerkin projection, the dimension of the problem is reduced but the model still requires knowledge of the full state u(t) to be resolved. To obtain a model that solely depends on the coordinates a(t), the dynamics f of the model have to be evaluated from the approximate form of the state (see Eq. (5)). This leads to the following system of r ODEs: If f is the right hand side of the incompressible Navier-Stokes in Eq. (1) and the POD basis is retrieved from the velocity field snapshot matrix, the approximation in (8) leads to the following ROM: where tensors L and Q are reduced versions of the dissipative and advective parts of the original Navier Stokes equations. Because operators L = ∆u(x, t) and Q = ∇(u(x, t)u(x, t)) are spatial operators (here u(x, t) ∈ R d refers to the velocity vector), they can be pre-assembled from the POD modes. In contrast with this favorable situation, other PDE problems might necessitate that the dynamics be evaluated from the reconstructed solution u(t) ≈ V r a(t), which can imply a significant cost and require additional reduction work, [18]. It is worth noting that the pressure term is not present in Eq. (9). This is due to the fact that, in the full order model, the pressure term accounts for the incompressibility condition (∇ · u = 0) and can be interpreted as a Lagrange multiplier used to ensure this constraint during the simulation.
Because the reconstructed solution is a linear combination of the POD modes which are all divergence-free, the incompressibility condition is always satisfied and the pressure term can then be ignored without significant accuracy loss ( [38,14]).

Limitations of POD-Galerkin models
Because POD-Galerkin models rely on the projection of the original model equations, they retain part of the structure from the full order physical model and this makes them robust to new initial conditions comparatively to purely data-based approaches. It should be noted however that approximation errors are embedded by construction in their structure. Equation (7) shows that an exact expression for the dynamics of the reduced coordinates is obtained by projection of the model equation on the POD basis. If the orthogonal complement to V r is not neglected, equation (7) reads: with Qu(t) F Pu(t) F if r is large enough. Considering a first order approximation of Eq. (10), the temporal dynamics of the POD coefficients can be evaluated from: with J Pu(t) being the Jacobian of the discretised nonlinear system f , evaluated at the state Pu(t). The residual term R accounts for the dynamics of the complement of the POD basis in the space spanned by the POD basis. Classical POD-Galerkin reduced models ignore this error term and directly compute the dynamics of the model from the approximate solution as presented in Section 2.1.
This approximation means that small errors on the dynamics will compound over time and lead to significant discrepancies between the true and simulated trajectories. This is especially true in the case of nonlinear dynamical systems where orthogonal projection on the POD basis can suppress an important part of the dynamics. In the case of the Navier Stokes equations, POD-Galerkin models have been shown to fail to reproduce the dynamics even in simple cases like the flow over a cylinder [14].

Non-Markovianity and Taken's theorem
The residual depends on information from a subspace orthogonal to the span of the POD basis. This means that an accurate correction model cannot be directly computed from the reduced state a(t). However, we leverage the fact that the information lost by projection of the full order state can be retrieved by considering past states of the system. This hypothesis is formalized by the Takens' theorem ( [4]), which states that, under mild conditions, the dynamics of a state vector can be reproduced by constructing a time-embedding from time-lagged observables: , with k large enough. 1 This is illustrated in Figure 1, which shows the Lorenz attractor observed via its embedded X-component. By constructing a 3-dimensional embedding of the obtained time-series, an attractor is obtained, which preserves the topology of the true attractor (e.g., symmetries, correlation dimension, etc.).
This suggests that the correction term R in POD-Galerkin reduced order models is non-Markovian and should consider past states of the system a(t − τ ), a(t − 2τ ), . . . . However, such a discrete time-embedding is not well suited in the context of continuous time models such as equation (9). Indeed, using discrete time steps in combination with adaptive step time marching schemes, such as the Runge-Kutta method [3], would require specific considerations about interpolation between the simulation steps and the required embedding steps. To address this time-continuity issue while retaining a non-Markovian correction structure, we propose to use delay differential equations (DDE) with a continuous embedding of the past information: with λ ∈ R + sufficiently large for the integral above to be bounded.
These equations retain information from past states of the system in a time-continuous manner and are used for a number of modeling applications such as epidemiology or population dynamics ( [8,11]). In this formulation, the dynamics g depend both on the partially observed states of the system a(t) and a memory variables y(t), corresponding to the integral of past observables, damped in time by an exponential decay. In fact, the memory term y can be defined in many ways, depending on the problem at hand. However, the exponential decay formulation was chosen because i) it provides the model with the ability to consider recent states of the system while older observations are discarded and ii) it can be solved by directly augmenting the original system with a second ODE for the memory: The exponential decay e (τ −t)λ acts as a filter of width 1/λ on the observables evolution. Other applications of the same augmented ODE exist in literature but with different purposes such as modelling of subgrid-scales in Large-Eddy Simulations (LES) [15] or to find unstable-steady solutions of the Navier-Stokes equations [16]. Finally, the augmented dynamical model is presented on the right hand side of the image.
It is worth noting that the memory variables y have the same dimension as the observations a in (13). This limitation might introduce a significant information bottleneck in the model. Indeed, Takens theorem states that the dimension required to obtain a satisfactory embedding can go as high as twice the intrinsic dimension of the true attractor.
Although there is no similar result for the continuous case, the limited dimension of the memory y may prevent deriving an accurate correction model. Moreover, we have made the convenient choice of the exponential kernel because it can be described by a linear ODE, although more complex kernels could be considered. As a result, we define an encoding map E : R r → R n E , used to lift the observations a to a higher dimensional space to increase the dimension of the memory and relax the constraint imposed by our choice of linear memory dynamics: 2 Using this modified DDE architecture to correct reduced order models, the correction operator R acting on the memory y(t) becomes an application from memory space to phase space: R n E → R r . Finally, the following augmented reduced order model is obtained: The concepts behind this final formulation are summarised on Figure 2 to help illustrate the idea.

Data driven learning of the residual
Although the motivation for the existence of the residual R and encoding map E was outlined above, little information about their form can be derived from the previous expressions. In this context, we leverage the universal function approximator property of artificial neural networks [9] to model the missing terms.

Neural Networks
Our work relies on feed-forward neural networks with a very simple structure. These networks can be used to learn any smooth nonlinear continuous application Ψ : R r → R k by optimising the weights of a sequence of L layers: Figure 3: Superposition of sinusoidal signals of frequencies 1, 10 and 100 Hz, filtered by an exponential decay with a decay rate λ = 50. The higher frequency (100 Hz) is filtered out while the other two are reproduced in the memory signal. This implies that recent events in memory due to the 100 Hz frequency are filtered out so that the memory is mostly driven by older events associated with lower frequencies. Note that the memory signal has been scaled by a factor λ for clarity.
where σ l is a nonlinear activation function, and the dimension of η l corresponds to the number of neurons in the layer l. It has been shown that, provided the dimension of the layer is high enough, the trainable parameters W l and b l can be optimised to approximate any function, [9]. The encoder E(a) and the residual R(y) are both approximated with neural networks with parameters θ E = {W l , b l }, l = 1, .., L E and θ R respectively.

Memory time scales
Physical systems often involve a variety of phenomena each evolving at different time scales. Capturing these phenomena can be critical to accurately model these systems, which is why the memory should be able to retain information at different rates. The time scales accounted for by the memory are driven by the parameter λ of the exponential decay in Eq. (12), which acts as a low-pass filter on the encoded trajectory (see Figure 3 for an illustration).
To retain information at different rates, λ can be adjusted for each observable E i (a). Equation (15) is then modified accordingly with the single λ parameter replaced with a diagonal matrix Λ ∈ R n E + whose entries can be learned to select time scales relevant to the system at hand. This finally results in the Complementary Deep ROM architecture:

Training strategy
Optimising the parameters of the model (18) requires consideration of past states of the system as well as their impact on the dynamics in the future. In some sense, the problem is similar to the optimisation of classical recurrent neural networks for the simulation of dynamical systems [63]. The major difference is that the aim of the present work is to derive a continuous time dynamical model, while recurrent networks have traditionally been used to model transitions between discrete time instants. As a result, we can not use standard backpropagation through time to optimise the model. Instead, the NeuralODE approach [40] is used. Initially proposed as the continuous equivalent of residual networks ( [32]), Neural ODEs are based on the well established adjoint backpropagation method ( [1]). It is briefly presented here while a more detailed derivation is given in Appendix A. The adjoint backpropagation method can be used to solve constrained optimization problems of the form: The method allows for the computation of the gradient d dθ J, at roughly twice the cost of integrating forward in time from t = 0 to t = T , through the following expression: where the adjoint state µ can be computed by integrating the adjoint ODE: Computing the gradient of J with respect to the parameters thus requires integration of the adjoint ODE backwards in time, from t = T to t = 0.
If the dynamics f are expressed by a Neural Network, the required vector-Jacobian products µ T ∂f ∂x t and µ T ∂f ∂θ t can seamlessly be evaluated through automatic differentiation by any classical deep learning framework.
As explained in Section 2.1, reduced models of the incompressible Navier-Stokes equations can be expressed directly in terms of the reduced coordinates a and evaluated through simple tensorial expressions (see equation (9)). This means that the ROM dynamics can be directly computed and back-propagated through once the reduced operators L and Q are assembled. The Deep ROM architecture (18) can then be optimised within the Neural ODE framework by concatenating the reduced state and the memory into a single state vector Since the last modes in the POD basis V are typically associated to low-energy dissipative scales [56], if r is not large enough, the ROM model retrieved by V r might be unstable and diverge after a few time integration steps. This potentially unstable primal model is however embedded in the optimization framework so that the memory is trained accordingly and it results in a stable, accurate, model.

Training Data and Memory Initialisation
To optimise the parameters of equation (25) through the Neural ODE approach, one only needs knowledge of the true trajectory of the reduced coordinates a (t). As presented in Section 2.1, this trajectory data can be obtained by projecting the solutions computed with the full order solver on the POD basis: If the snapshots are sampled with a time interval ∆ t , simulating the corrected ROM (25) for n t time steps leads to the following L 2 mean squared error: for which a gradient can be computed by integrating the adjoint equation (23). It should be noted that the memory term needs to be initialised by evaluating the memory integral at time t = 0, which requires knowledge of the past of the true trajectory: To be able to compute this integral, the infinite horizon of the memory integral can be relaxed by defining a finite τ min from the longest time scale λ min of the matrix Λ and a threshold 1, chosen to be small enough such that the relative error made on the initial memory is sufficiently small: Because this initial memory term directly depends on the parameters of the encoder and the matrix Λ, it needs to be re-computed at each training epoch. However, the impact of variations in the model's parameters on y(0) is not accounted for during the optimization and could lead to small errors in the approximation of the gradient of the loss (27). This might require limiting the learning rate to minimize variations on the initial condition between epochs.

Interpretation of the model
Before presenting the selected test cases and results, further justification and insights into the model are discussed in this Section. First, the Mori-Zwanzig formalism is introduced to frame our model in the context of dynamical systems theory. Then, we study how it can be compared to purely data driven approaches, such as reservoir computing.

The Mori-Zwanzig formalism
The Mori-Zwanzig formalism provides a closed form for the dynamics of partially observed systems by distinguishing three separate terms: where a(t) are observables of a system defined as a projection of the full order state onto the observable space and Ω is the projected part of the original dynamics. These two quantities can respectively be identified as the reduced coordinates and dynamics of classical POD-Galerkin models. The remaining terms account for the impact of the non-observed coordinates of the system on the resolved dynamics. K represents the dynamical exchanges between resolved and unresolved dynamics during the simulation, while F accounts for the missing information from the initial condition. This term vanishes if the initial condition and dynamics are fully resolved at t = 0, in which case the Mori-Zwanzig equation becomes: The residual model R proposed above can then be identified with the memory integral defined by the application K, providing a strong connection with our non-Markovian correction hypothesis. Framing the CD-ROM architecture in the context of the Mori-Zwanzig formalism further justifies our modeling choices. Yet, it does not provide additional information on the form of correction model since little information is known about the convolution operator which can be infinite dimensional in certain cases. Choices about the structure of K need to be made. In this work, assumptions are made about the vanishing impact of past states of the system on the residual model, accounted for by the matrix Λ in our approach.

Deep Learning interpretation of the model
In this paragraph, the link between the Complementary Deep ROM architecture and classical deep learning models is discussed. Numerous methods have been proposed to model sequential data, each relying on a specific mechanism to extract and retain meaningful information from the past states of the system. The most popular architectures, GRU and LSTM ( [24,10,25]), both use a combination of gating mechanisms to learn long term dependencies in a sequence.
Other approaches like reservoir computing rely on an underlying dynamical system forced by the sequence data to predict the required output. In fact, strong similarities can be identified between our approach and echo state networks (ESN [19,52]), a widespread reservoir computing architecture. ESNs are based on the simulation of random dynamics described by matrices W in and W R . The matrix W in is used to encode some data x t in higher dimension, while the matrix W R holds the weights of the reservoir used to advance the state (memory) y t in time: For the sake of comparison, the nonlinear activation function σ is dropped, and we consider that equation (33) results from the Euler integration with time-step ∆ t of a continuous system describing the dynamics of the memory dy dt : = y t + ∆ t dy dt (35) which leads to the following expression: Because the spectral radius of the matrix W R is constrained to be less than unity, all the eigenvalues of the operator have negative real parts, leading the ESN to have memory dissipation properties similar to those of our model. To underline this similarity, we can express the state of the ESN at time t by diagonalising the operator W R = PΛ R P −1 : Thus, a parallel between our model and the ESN is outlined. Major differences remain in the fact that the dynamics (W R ) and encoding matrix (W in ) are not optimised during the ESN training. It should also be noted that nonlinearity is introduced in the ESN dynamics through the σ activation function, while our model is based on linear memory dynamics and on a nonlinear encoding operator. These comparisons between direct deep learning methods and the continuous correction approach help build intuition about the role of each term in the model. The encoder can be compared with the input gate of an LSTM, or the W in matrix of the ESN, the Λ matrix provides a tunable forget mechanism, while the residual term R plays the role of the output operator in our "continuous recurrent network".

Case presentation and reduced models
In this Section, we introduce two flow problems to illustrate the benefits of the Complementary Deep ROM architecture. The first case is a standard configuration of the flow over a cylinder, often used as a benchmark for reduction methods. The second case is the fluidic pinball flow, introduced in [36] for the development of new control strategies.

Flow over a cylinder
The two-dimensional incompressible flow over a cylinder has been extensively studied in the context of reduced order modeling [7,14] and model identification [31,42] which makes this test case a good initial benchmark for the proposed correction method. The flow is simulated at a Reynolds number of Re = 100 based on the cylinder diameter and the velocity of the incoming flow. In this regime, the flow is laminar and exhibits vortex shedding in the wake of the cylinder.  The flow is governed by the incompressible Navier-Stokes equations, here solved using the FEniCs finite elements solver [21,28]. The retained mesh is shown on Figure 4. It is rectangle-shaped, spanning from x = −5 to x = 15 in the streamwise direction, and from y = −5 to y = 5 in the transversal direction. The cylinder has a diameter D = 1, centered around the origin. The inflow is modeled as a uniform axial flow (u(−5, y) = [U ∞ , 0]) and a freeslip condition is used for the lower and upper boundaries of the rectangle while a no-slip condition is enforced at the cylinder surface. Finally, a stress-free condition is used for the outlet.
A reduced order modeling strategy is employed to obtain a baseline model for the correction approach. The main results are summarised here but we refer to [14] for details on the reduction strategy. The vortex shedding regime of the cylinder flow is simulated with the FEniCs solver to obtain snapshot data and compute the POD modes. The first two modes, accounting for more than 95% of the Frobenius norm of the snapshot matrix are selected. The steady solution of the system is computed with a Newton method to construct a so-called shift mode (a ∆ ). This mode is computed as a vector orthogonal to the plane described by the first two modes, pointing to the base flow solution u b , and serves as a support for the simulation of the transition of the system from its steady state to the vortex shedding limit cycle. A three-dimensional POD basis is thus finally obtained. They are shown in Fig. 5 in terms of the vorticity field.
Through Galerkin projection of the discretised Navier-Stokes equations, a system of 3 coupled ODEs is obtained, describing the dynamics of the reduced coordinates vector. To simulate the transition from the base flow to the vortex shedding regime, we chose an initial condition u 0 close to the base flow u b as suggested by Loiseau [42]: where v 1 is the first POD mode, and > 0 is a small coefficient used to perturb the unstable base flow. The results of the simulation are displayed on Figure 6. Even though the equation model is able to simulate the transient dynamics, its trajectory strongly diverges from the projected snapshot data. The transition starts much later than in the full order simulation, due to a growth rate of the ROM's transition lower than what it should be. Another significant issue with the model is its stabilisation around the limit cycle (a ∆ = 0) where an overshoot can be observed before the vortex shedding regime is established, which is not observed in the snapshot data. These discrepancies between the two trajectories can be attributed to the ignored residual term, making this model a good baseline for the correction approach illustrated in Section 6.

Fluidic pinball
The second case used to demonstrate the approach is the so-called fluidic pinball. Initially proposed as a challenging test bed for the development of control laws [36], the fluidic pinball case offers a good trade-off between complexity of its dynamics and interpretability [47,59]. The simulation domain (Figure 7) is composed of three equidistant cylinders, each generating vortices in the wake which interact to create rich dynamics. The mesh used for the simulation was provided by the authors of [41]. Displayed on Figure 7, the domain is a rectangle spanning from x = −6 to x = 20 in the streamwise direction, and y = −6 to y = 6 in the transverse direction. Three identical cylinders with diameter D = 1 are arranged in an equilateral triangle, with centers' coordinates (0, 0.75), (0, −0.75) and (−1.25, 0) respectively. The boundary conditions are identical to those of the cylinder case in Sec. 5.1. The inflow is modeled as a uniform axial flow, the upper and lower boundaries of the computational domain are modeled as free-slip, while a no-slip condition is used for the walls of the three cylinders and the outlet is modeled as stress-free.
The flow is simulated at a Reynolds number of Re = 130 for 1800 time units (t.u.) in the chaotic regime. At this Reynolds number, the flow dynamics have been shown to be chaotic [59], which makes it a challenging problem for any correction model as the smallest error on the dynamics will make the simulated trajectory diverge exponentially fast in time from the truth. This kind of chaotic problems are starting to get traction as interesting benchmarks for forecasting tasks and modeling problems [60]. It should also be noted that the pinball flow is much more complex than the cylinder, as evidenced from Figure 8 where it is seen that many POD modes are required to account for a significant part of the energy. Almost a thousand modes would be required to capture 99% of the Frobenius norm of the snapshot data, where only 8 are required to achieve the same accuracy in the cylinder case.
We chose to build a POD-Galerkin model of this flow using only the first 10 POD modes. Although this choice is somehow arbitrary, it was made to challenge the correction method, as the mean reconstruction error of about 60% has an important impact on the approximated dynamics. Indeed, the obtained reduced model quickly separates from the original trajectory, as expected from chaotic dynamics. More problematic is the fact that the reduced model is very unstable and diverges after 70 t.u. of simulation, as shown on Figure 9. Application of the correction method is presented in Section 6, a discussion on the impact of the number of modes is provided in Appendix B, while additional information on the computational cost of the approach is given in Appendix C.  6 Results and discussion

Cylinder case
The reduced model derived in Section 5 above was shown to be efficient for the simulation of the vortex shedding regime, but not suited to the simulation of transient dynamics. To apply the proposed correction procedure to this model, the modeling terms introduced in Sec. 3 are added to the reduced model. The dimension of the memory is chosen to be ten times the dimension of the reduced state. The residual R and encoder E are defined as multi layer perceptrons, using the Rectified Linear Unit activation function. Finally, the diagonal of the memory matrix Λ is initialised at random from a normal distribution. The model is trained with the Adam optimizer, and the trajectory loss introduced earlier (Eq. (27)) is used in combination with additional regularization terms. Details on the regularization of the loss are provided in Appendix D.
Since the targeted trajectory for this case is the transition trajectory, which starts close to the base flow of the system, the initial memory can be computed with minimal error through the following integral: where a b = V T r u b are the reduced coordinates of the base flow and τ min is the longest time horizon defined by the Λ matrix as in Equation (29). Finally, the parameters of the models are progressively optimised to reproduce the true transition trajectory as shown on Figure 10.
The results of the corrected model's performance on the training and test trajectories are presented on Figure 11. It can be seen that the corrected model follows the training trajectory perfectly, triggering the transition at the right time, and correcting the oscillations of the original ROM during the stabilisation on the limit cycle. Moreover, the graph shows  Finally, a test trajectory was simulated in FEniCs by taking a random initial condition in the phase space spanned by the 3 selected modes. The first 8 t.u. of DNS simulation are used to initialise the memory following Equation (30). Figure 12 presents the performance of each model on this trajectory. Because the starting point is not close to the base flow, the uncorrected ROM instantly exhibits transient dynamics, however, the growth rate of the transition is still inaccurate and the shift mode's trajectory presents the same non-physical oscillations around the limit cycle observed in the training trajectory. The corrected model does much better than the original ROM, simulating a more accurate transition, and stabilising on the limit cycle almost perfectly.

Fluidic pinball results
While the cylinder case discussed above offers a simple test bed for the presentation of the approach and its potential, it has already been shown that very parsimonious models could be used to model its dynamics [42], which suggests this configuration might not require a high dimensional neural networks to learn a correction term. The case of the fluidic pinball is more challenging and can better underline the ability of our method to handle complex physics. As presented in the previous section, the POD-Galerkin approach is not well suited to the reduction of the pinball case. The number of modes required to reconstruct the snapshot data with a satisfying accuracy is very large and using a small number of modes leads to an unstable model.
To apply the correction approach to the pinball case, a correction model is built. The encoder and residual models are multi layer perceptrons and the Sigmoid Linear Unit activation function is used as it leads to smoother integration. The weights of the residual model are initialised to be close to 0 so that the ROM is initially almost uncorrected. Diagonal Figure 12: Simulation results obtained on a test trajectory. As in Figure 11, the first graph represents the phase space trajectory, while the time evolution of the shift mode's coefficient is presented on the second graph. As presented in the previous section, the uncorrected 10-mode ROM is unstable and diverges after 70 t.u. of simulation. As a result, trying to optimise the correction model for long trajectories directly would lead to a very unstable training process. To address this issue and stabilise the training, the corrected ROM is trained on sub-trajectories of only 10 t.u. This length is chosen as it is short enough for the model to remain stable and long enough for the impact of the encoder on the memory to be accounted for. Once a good correction model has been trained for 10 t.u. long trajectories, the length is progressively augmented to attain the target horizon of 120 t.u. Besides the stabilisation of the training, using sub-trajectories also allows for parallel training. Multiple sub-trajectories can be sampled from the snapshot data and integrated in parallel on a GPU to dramatically speed-up training. The pinball correction model was trained to follow true trajectories up to 120 t.u. Results of the simulated flow fields are presented on Figure 13. It can be seen that the projection on the 10-dimensional POD basis effectively filters part of the spatial structures, leading to the divergence of the uncorrected reduced model. In contrast, the corrected ROM is able to reproduce the projected flow field accurately.
Trajectory results are also presented on Figure 14. The model starts quickly diverging from the true trajectory after the training horizon (120 t.u.) as is expected from the chaotic nature of the problem. Deriving a model to perfectly follow the DNS trajectory would here make little modeling sense. More interesting is the fact that, despite leaving the training trajectory, the corrected model does not diverge even when integrated for over 1000 t.u. This suggests that the dynamics correction learned by the model has some physical consistency, dissipating the necessary energy which would otherwise have made the simulation diverge in the uncorrected case.
The intuition that the model was able to learn a Physics-compatible correction is confirmed when looking at the statistics of the attractor spanned by the CD-ROM's trajectory. Using the nolitsa library [43], the correlation dimension ( [5]) was estimated, as well as the maximum Lyapunov exponent ( [6]) of the corrected and true trajectories. The results Figure 14: Results of the correction approach applied to the 10-mode ROM of the fluidic pinball. The initial condition is a random point selected in the training trajectory. Plots (a),(b),(c) describe the coefficients of the first three POD modes simulated with different models. Plot (d) presents the relative Euclidean distance between the true data and the trajectory simulated with the CD-ROM. Figure 15: Statistics of the 10-dimensional attractors described by the projected DNS data (black) and the CD-ROM (blue). Left : Average number of points g in a sphere of radius r, the data follows the law g = r C d where C d corresponds to the correlation dimension of the attractor. Right : Time evolution of the average distance d between trajectories starting from arbitrarily close initial conditions on the attractor. The data follows the law d = e tle where l e corresponds to the maximum lyapunov exponent of the system at hand. are shown on Figure 15 where the model is seen to reproduce well the characteristics of the true attractor. One can also look at the probability distributions of the mode's amplitudes, presented on Figure 16. Once again, the simulated trajectory reproduces the results of the true simulation.
Note that long trajectories (several hundreds of t.u.) were simulated with the corrected model to obtain these statistics. Due to the chaotic nature of the problem and the length of the integration period, the model has visited parts of the attractor different from those seen in the training trajectory. Despite this, the CD-ROM remains stable and describes an attractor with statistics that are very similar to those of the true attractor, further supporting the approach for the reduced modeling of complex dynamics.

Conclusion
This paper was concerned with an augmented reduced order modeling strategy based on the hybridisation of the classical Galerkin projection method and simple neural networks. By studying the limitations of the Galerkin projection, we established links with the field of partially observed systems, which led us to use past observables of the studied system as a critical ingredient for the correction of Galerkin models. Building on this result, the CD-ROM architecture was proposed to extract and exploit useful information from the system trajectory, by embedding the model with a delay differential equation structure. Moreover, the training strategy based on adjoint optimisation ensures a-posteriori performance of the model on the training trajectory.
The CD-ROM approach was demonstrated on two simple CFD test cases, namely: the flow over a cylinder and the fluidic pinball. Numerical experiments have shown that the corrected models were able to capture the true dynamics with a high degree of accuracy, reproducing the true transition in the case of the cylinder flow, and following the training trajectory for multiple Lyapunov times in the fluidic pinball case. Moreover, these experiments outlined the reliability of the corrected model as it performed better than its uncorrected counterpart even outside the training conditions. The results obtained with the fluidic pinball are particularly encouraging, we showed that the correction model was able to stabilise the original Galerkin model in a consistent physical manner, as the attractor simulated with the CD-ROM approach is similar to the original attractor.
On-going work investigates the application of the CD-ROM approach on a more concrete industrial use case, to demonstrate the capabilities of the method on a real-world application. As a more general perspective, one could investigate the links between the Mori-Zwanzig theory and the CD-ROM approach, as the similarities discussed in the paper suggest the existence of a formal theoretical justification for the model. Another perspective is to consider additional filter shapes in the time filtering mechanism of the CD-ROM approach, in order to improve the memory aggregation process.

A.1 Adjoint Backpropagation
The adjoint backpropagation algorithm on which Neural ODEs are based can be used to solve problems of the form 3 : Generally, this kind of problems is solved with a gradient descent method, which requires evaluation of the following gradient: Because the derivative ofẋ(t) = g(x; θ) is parameterised by θ, there is an implicit relation x = x(t; θ). The estimation of the sensitivity ∂x ∂θ | T requires the consideration of the impact of θ on the whole time-integration, from t = 0 to t = T . The adjoint backpropagation is used to avoid computing this term, which is done by evaluating the sensitivity of the following Lagrangian: The vector of Lagrangian multiplier µ is a function of time. Distributing the product in the integral and integrating the first term by parts leads to an expression where the sensitivity ∂x ∂θ T can be isolated: A so-called adjoint equation can be derived to avoid having to solve for the sensitivity ∂x ∂θ . Enforcing a vanishing variation of the Lagrangian wrt x at optimality yields: Solving the adjoint equation (45) for the values of µ(t) leads most of the terms in Equation (44) to vanish, so that: As mentioned earlier, both this integral and the adjoint equation can be easily evaluated ifẋ(t) is approximated by a neural network as the required vector-Jacobian products µ ∂g ∂x and µ(t) ∂g ∂θ can be easily computed in a deep learning framework.

A.2 Adaptive checkpoint adjoint
The original NeuralODE paper ( [40]) proposed integrating forward in time to obtain the initial condition of the adjoint equation evaluated from x(T ), while discarding intermediate values x(t). This choice was made with the goal of reducing the memory footprint of the method. However, discarding the intermediate time-steps x(t) means that they have to be recomputed during backpropagation to evaluate the adjoint ODE (45), as well as the gradient integral (46). This can be done as a single backward in time integration by concatenating the different quantities (x, µ, dL dθ ) in a single state vector: Not only does this increase the computational cost of the method, but it can also lead to erroneous gradients, as numerical errors during integration can lead to differences between the forward and backward trajectories for x(t). This was observed in [55] and is illustrated on Figure 17, which shows that, despite using the same parameters, the forward and backward time trajectories diverge due to numerical errors. To address this issue, one can use the Adaptive Checkpoint Adjoint method. This method retains the intermediate integration steps (t, x(t)) and simply evaluates the integrals for the adjoint µ and the gradient dJ dθ at the time steps selected during the forward integration, using the forward trajectory x(t).
Although this increases the memory footprint of the method, it limits the computational cost of the backward pass, as the time steps are already selected, and the trajectory x(t) does not have to be integrated another time. Most importantly, this limits the numerical errors introduced by the integration schemes, which can have a significant impact on the training, especially in the case of chaotic time-derivatives. In Section 5, we present a 10-mode reduced order model of the fluidic pinball problem. The number of modes is chosen somewhat arbitrarily to challenge the method. The more modes are used to model the flow, the better the model will be at reproducing the true dynamics, reducing the complexity of the required residual term. In this Appendix, we provide more insights into the impact of the number of modes on the modeling problem to clarify the choice of using 10 modes to model the pinball flow.  Figure 18 presents the performance of different reduced models of the pinball flow. It is clearly seen that increasing the number of modes is beneficial for the performance of the reduced models. The magnitude of the closure term, as well as the speed at which the reduced model diverges from the true trajectory, are reduced when the number of modes increases. This Figure also shows that a higher number of modes leads to more stable reduced models. However, further experiments showed that even well resolved models such as the one using 173 modes were not stable and would diverge in certain conditions.

B.2 Training convergence
To add to the argument, a comparison of the training convergence between two models using different number of modes is discussed. Two models are built, using respectively 10 and 55 POD modes, and the correction architecture is trained using the same parameters, presented in Table 1.
In both cases, the same loss is used, combining the optimisation of the distance between simulated and true trajectory with the residual regularization introduced in Appendix D: where c = 0.1 is a constant weighting the importance of the stochastic residual regularization term w.r.t. the trajectory loss. Notice that, to ease the comparison between the two models, the loss is scaled by the number of modes used in the ROM. The two models are trained in the same fashion, sub-trajectories of a hundred time steps are sampled in the training base and the loss (50) is optimised until a chosen threshold (5 10 −4 ) is reached, at which point the length of the sub-trajectories is increased by fifty time steps. This process is repeated until the model is able to reproduce sub-trajectories of a thousand time steps. Figure 19 presents the evolution of the loss for the two models, as expected, the 55-mode model is quicker to train as it first reaches the threshold in 1400 epochs, while it takes the 10-mode model more than twice the number of epochs to achieve the same performance. Similarly, the 55-mode model achieves the required precision on trajectories of a 1000 time steps in only 4700 epochs, which is again more than twice as fast as the 10-mode model.
These results confirm the interest of only using the first 10 modes to challenge our correction approach. The relatively high magnitude of the residual term to be learned, the instability embedded in the model and the low degree of reso- lution of the problem (only 42% of the snapshot information) all constitute significant complexities which could arise in real world applications. While we showed that using a higher number of modes would simplify the modeling problems, this constraining choice helps demonstrate the applicability of the CD-ROM method to challenging modeling problems.

C Computational cost
The overall goal of the approach is to obtain an accurate and inexpensive model for the simulation of complex dynamical systems. As was seen above, the POD-Galerkin approach is a very efficient way of limiting computational costs by reducing the dimension of the system of interest. However, we showed that this method often produces inaccurate models when applied to nonlinear systems limiting its applicability to real-life problems.
Using our proposed CD-ROM architecture, we showed that a correction model can be trained to complement standard POD-Galerkin models and improve their accuracy.
We compare the cost of integrating our model with a 173-mode uncorrected model, which performs much better than the 10-mode uncorrected ROM but still diverges in certain conditions. The results are presented on Figure 20 and show that, although the corrected model does diverge from the true trajectory after 1 Lyapunov time, it performs better than its uncorrected counterparts. Indeed, the 173-mode model diverges earlier than even the 10-mode model in this particular condition, and is more expensive to simulate than the CD-ROM. The fact that the CD-ROM is more stable than the 173-mode model also adds to the interest of the proposed method.

D.1 Residual regularization
We observed that, only training the model to follow the trajectory data can lead to poor local optima with large corrections applied to the original model. Since the magnitude of the correction is meant to be small when r is sufficiently large, this tends to indicate over-fitting. This is an issue as such a model does not capture the true dynamics, and will quickly diverge when evaluated on conditions different from the training trajectory. To address this, a regularization term can be added in the model, to limit the magnitude of the corrections applied to the ROM: where n t is the number of time steps of length ∆ t in the optimised trajectory, a is the training trajectory data and ρ is an hyper-parameter chosen to balance the importance of the regularization w.r.t. the rest of the loss. Although this fairly simple approach already helped guiding the training, we obtained better results by adding a custom regularization loss. This loss is based on the approximate value of the residual on the true trajectory, which can be obtained by computing the time-derivative of the true reduced coordinatesȧ through finite differences, and computing the difference with the derivative defined by the uncorrected ROM: With this correction data, the memory can be evaluated from the true trajectory integrating in time the following ODE: One can then define a regularization term L corr for the loss: The regularized loss definition was observed to lead to models with better generalization properties. Note that L corr can be computed for a small random subset of the a training batch at each epoch, to accelerate training while keeping a "stochastic" regularization for the residual.

D.2 Encoded space regularization
Regularizing the encoded space can both help smooth the training process and increase the robustness of the model to unseen conditions. Taking inspiration from existing work ( [35]), we propose to add the identity function to the encoder model. This means that useful information is embedded in memory in the form of a time convolution of the past resolved states of the system. The encoded state E(a t ) is then constituted of both the reduced state a t and a learnable nonlinear transformation MLP(·, θ E ) of it: With this structure, the encoder will only be learning additional nonlinear transformations of the state, simplifying the training and introducing a level of structure in the encoded space as we ensure its first dimensions are coherent with the phase space.

Funding
The research work of E. Menier was funded through the HSA project, at IRT SystemX, Paris-Saclay, France, and supported by the "Programme d'Investissements d'Avenir".