Mode decomposition of flows with shocks

Guillaume Chauvat

, Peter J. Schmid

and Ardeshir Hanifi

1

KTH Royal Institute of Technology, Engineering Mechanics, FLOW Centre and Swedish e-Science Research Centre (SeRC)

2

Imperial College London, Department of Mathematics

To be submitted

In the present work, we introduce a modification of the Proper Orthogonal Decomposition (POD) to allow the decomposition of data containing moving discontinuities without the typical drawbacks of POD such as the appearance of Gibbs’ phenomena near the interface. We use a shock-fitting approach and interpolate the data onto a grid following the discontinuities. Solutions of a model problem based on the Burgers’ equation and unsteady supersonic flow around a cylinder are used to illustrate the new method. The shock-fitting approach produces a faster decay of singular values and sharper discontinuities in the identified modes than the standard POD. The Gibbs’ phenomenon has been largely eliminated.

1. Introduction

In many technological areas, flow of different type of fluids, separated by an interface or discontinuity, is an essential part of the system. Understanding the dynamics of flow structures in those applications is key to improving their effi- ciency. The increasing performance of computers has allowed the computation of ever larger and more complex flow configurations. Nonetheless calculating the solution to the Navier–Stokes equations remains prohibitively expensive in many cases. For this reason, reduced-order models (ROMs) have been used in areas such as flow control, decreasing the dimensionality of the problem to a set of physically important modes.

The choice of the reduced-order basis defines what physical phenomena can be captured by those methods. The Proper Orthogonal Decomposition (POD) (Lumley 1970; Sirovich 1987; Rowley & Dawson 2017) is widely used for this purpose. It is usually calculated as the Singular Value Decomposition (SVD) of a matrix of snapshots (Sirovich 1987) that can be taken from numerical simulations or experimental measurements. Among applications are flow con- trol (Rowley & Dawson 2017) or design optimisation (Iuliano & Quagliarella 2013). However, the use of POD with advection-dominated equations can be problematic: the POD, as a purely spatial method, is not properly equipped to efficiently represent travelling waves. An infinite travelling sine wave can be represented with two modes, but more modes are necessary for a more complex

123

function shape even if it is simply being advected at a constant velocity since in this case the POD is given by the Fourier transform of the wave form. Variants of the POD applied in the frequency domain (Towne et al. 2018) or combined with Multiresolution Analysis (Mendez et al. 2019) can improve on this issue by taking into account some temporal information. Specifically for advection- diffusion equations, Lagrangian formulations have also been proposed (Lu &

Tartakovsky 2020), in which the data points move along the characteristics of the equations.

Difficulties become more severe when considering travelling discontinuities, due to the linear nature of POD. A POD can be used but an accurate representation of the shock motion may require the superposition of many oscillatory modes, similarly to the representation of a discontinuous function in terms of Fourier modes. In addition this approach is unable to reconstruct a shock extending beyond its range in the set used for generating the POD basis (Li & Zhang 2016). It has been proposed to determine models using a POD basis for the smooth regions and a full-order model of the regions containing moving shocks (Lucia et al. 2001, 2003). This may be effective when the shock motion is small, and requires a-priori knowledge of its range. Another possible approach, when the data on only one side of the discontinuity is of interest, is the masking of the data in the adjacent region Balajewicz & Farhat (2014). This in effect makes the discontinuity disappear. However, in most cases, especially in aero- nautical applications when a shock wave is present around a body or inside an engine inlet for example, the flow on both sides of the discontinuity is of interest. Several discontinuities can also be present at the same time, and they do not necessarily divide the domain into separate subdomains, which makes this approach less practical. As an alternative we develop a shock-fitting POD in which the data aggregated in the snapshot matrix is not sampled at fixed points but on a dynamic grid fitted to the shocks, allowing the data to remain continuous across snapshots.

In the following sections, we first describe and discuss the method in relation to the POD. We then demonstrate the advantages and feasibility of the method in a synthetic scenario using Burgers’ equation with a single shock, before illustrating it on a more realistic two-dimensional solution of the compressible Navier–Stokes equations with several intersecting shocks. Finally, we summarize our technique and present possible extensions and generalizations.

2. Shock-fitting POD

2.1. Proper Orthogonal Decomposition and Singular Value Decomposition

Among the tools for a modal decomposition of numerical or experimental data,

the Proper Orthogonal Decomposition (POD) ranks among the most widely

used. Its origin can be traced to early attempts at decorrelating structures in

sample sequences, and various forms of this decomposition are known, under its

respective names, throughout the quantitative sciences and many application fields.

From a mathematical point of view, the Proper Orthogonal Decomposition is equivalent to a Singular Value Decomposition (SVD) of a snapshot matrix A. While often linked to an equivalent eigenvalue problem for the associated correlation matrices A

A or AA

, we will continue with the SVD viewpoint.

The generally rectangular data matrix A, consisting of m rows and n columns, is decomposed into a product of three matrices according to

A = U ΣV

. (1)

where

denotes the transpose operation. With A composed of n temporal snapshots, each consisting of m spatial data points, the decomposition produces an orthogonal matrix U ∈ R

, a diagonal matrix Σ ∈ R

, and an orthog- onal matrix V ∈ R

. The entries of the diagonal matrix Σ, referred to as the singular values, are non-negative and ordered in descending order on the diagonal of Σ. The matrix U contains the spatial modes as its columns, while the columns of V contain their associated temporal evolution. The singular values can be thought of as the amplitude of the modes, since the respective columns of U and V are normalized. This type of decomposition allows the expression of the data matrix as a superposition of rank-1 subprocesses according to

A =

n

X

i=1

σ

u

v

(2)

where u

and v

denote the i-th column of U and V, respectively. This repre- sentation also allows the approximation of the data matrix A in reduced-order form, by truncating the sum in the above expression and, hence, neglecting modes that have insignificant amplitudes. In matrix form, we have a rank-k representation of the data matrix

A

= U

Σ

V

(3)

where k < min(m, n), U ∈ R

, Σ

∈ R

and V ∈ R

.

It is worth mentioning that the spatial arrangement of the data points in the columns of matrix A is arbitrary, but maintained from column to column. In other words, the decomposition can process data on any kind of mesh, structured or unstructured, equispaced or clustered. The choice of mesh, however, carries over into the spatial modes contained in U, as they inherited the same spatial arrangement as the columns of the data matrix A.

When processing flow data with interfaces, moving boundaries or shocks

on fixed grids, we encounter the situation where data points in the vicinity of

the interface finds themselves on different sides of the interface and therefore

show rapid changes in the physical variables they track. For example, in the

case of a shock sweeping across grid points, we may record pre-shock conditions

for some time and then abruptly change to post-shock conditions. This will

lead to artifacts in the modal structures when processed by the singular value

decomposition. More precisely, oscillatory features appear in the spatial modes which, via superposition, represent the rapid changes in flow regime. This is akin to the Gibbs phenomenon encountered when representing functions with discontinuities by oscillatory basis functions.

To remedy this shortcoming for the decomposition of interface-dominated flows we revert back to a Lagrangian viewpoint and consider a set of data points divided into unstructured logical grids delineated by the various interfaces. By moving the points within each deforming domain and re-interpolating the flow field variable onto them, we circumvent the crossing of data points through interfaces and thus the artifacts mentioned above. To this end, we need to detect the current location of the interface and use it to triangulate the resulting subdomains.

2.2. Shock detection and tracking

In our special case, the interfaces consist of shocks, and a variety of shock- detection techniques are readily available from the literature on computational techniques for compressible fluid dynamics or shock dynamics. While any physical quantity that exhibits a discontinuity across the shock, is a viable criterion for shock detection, we choose the pressure and progressively determine the points along the shock, searching for new points in the vicinity of already identified points. These identified points are then interpolated by a parametric cubic spline to produce a continuous representation, which in turn can be sampled at will. We stress that we determine these points until a wall or the edge of the computational domain is reached, or until the discontinuity ceases to exist. This technique combines efficiency and flexibility to treat a variety of flow problems that is characterized by unsteady shock motion.

For more complex shock dynamics, for example, with intersecting shocks, each branch is followed and parameterized separately, starting from the intersec- tion point. As long as the overall shock topology does not change qualitatively, the interpolation points can be used to drive the domain division and meshing.

For detection of special intersection points, we employ the maximum of the Q

criterion (Ando 2000). At the moment, we cannot accommodate shock motion that results in the creation and annihilation of shock-bounded domains.

Techniques for these additional complexities will be developed in a future effort.

2.3. Points generation and tracking

Based on the identified and parameterized interface location, we generate an

unstructured mesh, issuing from the detected shocks. A starting grid is first

generated by distributing points in a regular fashion along the average location of

the shocks and domain boundaries, as well as in the interior of the computational

domain. Inspired by the mesh generation algorithm DistMesh (Persson & Strang

2004), these points are then iteratively Delaunay triangulated and moved. At

each iteration step, an attractive force, linearly proportional to the distance of

nearest-neighbor points, is applied and a global equilibrium solution is sought.

A zero resultant force between the grid points yields a steady solution and a well- structured grid. Formulating the adjacency matrix of the grid, and dimensionally splitted it into Cartesian coordinate components leads to a system of equations of the form

A

x = B

x

(4)

A

y = B

y

(5)

where (x, y) = (x

, y

)

are the coordinates of the n grid points and (x

, y

) are the coordinates of the corresponding initial points. The entries in A

and A

can be derived from the adjacency matrix A, and represent the x- and y -component of the graph Laplacian of the grid. Points on the interface and boundaries of the computational domain are restricted from moving, and their corresponding rows have to be replaced by a no-force condition (or row of an identity matrix) to reflect this constraint. The matrices B

and B

accomplish this detail in the above equations. Rather than iterating a dynamical system to steady state (as in DistMesh), we solve directly for the equilibrium solution using a Krylov-subspace GMRES technique. Matrix multiplications are performed efficiently in sparse matrix arithmetic.

The x and y coordinates are treated separately to allow the pinned bound- ary points to move freely along the domain boundaries, in our case, either horizontally for the top and bottom boundaries vertically for the inlet and outlet boundaries. A nonlinear force model may have been used for a better grid point distribution; our choice of a linear force model has been motivated by the applicability of a time-efficient linear solver. Once the points have settled in their equilibrium position, the triangulation is regenerated, and the process is repeated until convergence of the point coordinates to a prescribed tolerance is achieved.

The resulting reference grid can then be modified as the shocks move by solving (4–5) with associated boundary conditions defined by the shock locations for each snapshot. The triangulation is kept the same as in the reference grid such that all points move continuously with shock displacements. Crucially, the moving grid points generally do not cross the shocks, and hence the data measured at a specific grid point avoids discontinuities in time.

A low resolution example of a deformed grid is displayed in figure 1. The points stay on the same side of the shocks in all snapshots, allowing the data, sampled at their location, to evolve continuously and smoothly.

2.4. Advantages and limitations

The present method allows the natural representation of sharp discontinuities

moving in time from one snapshot to another. Instead of requiring a large

number of modes to approximate the discontinuities, they will instead already

be reflected in each mode. Moreover, we avoid the contamination of smooth

parts of the flow field by oscillations stemming from a nearby discontinuity. As

Figure 1: Deformed interpolation grid at two different times. The shocks are printed in blue. For clarity the grid shown here only uses 3% of the points of the full resolution case.

a result, the thus identified modes contain far fewer oscillations close to the discontinuities, and the values of the modes at points close to a discontinuity are not contaminated by structures on the other side of the discontinuity, as would be the case for using a standard POD.

Despite these advantages, the present method has restrictions: it assumes that the topology of the shocks shows some similarity in all snapshots. While the shocks can change shape, even violently so, every tracked shock must be recognizable in all snapshots. If shocks intersect each other, they must do so in a similar fashion. Regardless of these restrictions, we can still analyze a wide range of flow fields governed by moving interfaces, such as encountered in complex shock motion, in fluid-structure interactions and in flow over membranes, among many others. Importantly, the computational method does not rely on a subdivision of the full domain into sub-domains. It can readily cope with a shock that ceases to exist within the the computational domain instead of being attached to a boundary.

In the present variant of the method, the spatial information is the grid is not taken into account in the modal decomposition. The effect of this omission is left for a future effort. In this fully Lagrangian approach we embed the shock motion into the snapshot matrix not only by its flow variables on a moving grid, but also by the geometry of the moving grid itself. In this case, the coordinates of the interpolation points will be concatenates with the associated flow fields. In this manner, we track the flow field and its evaluation points together in time. This technique, as well as the treatment of topology changes in the regions limited by the interfaces (such as finite region enclosed by shocks disappearing due a collision of these shocks), will be the topic of a future study.

The reduced variant considered here already furnishes many advantages over

standard techniques, as we will see below.

3. Results

We showcase the computational technique and present results of our analysis for two different cases. The first case is a model problem where we artificially introduce a time-dependent discontinuity. The second case is a full flow case exhibiting a rather complex and unsteady behaviour which includes shock motion and vortex shedding.

3.1. Burgers’ equation

In order to validate and assess our new method, we start by analyzing solutions of the two-dimensional Burgers’ equation

∂u

∂t + u · ∇u = ν∇

u − 1.2 ≤ x ≤ 1.2, −0.5 ≤ y ≤ 0.5 (6) with ν = 5 × 10

and Dirichlet conditions on the boundaries of the domain.

The boundary conditions were chosen such that the solution of the above equation has a strong discontinuity within the computational domain. The results presented here correspond to the following boundary conditions:

u = 0.3x + (0.4 + 0.8y) cos (7x + 8y + 2t), (7) v = −2y cos (y) + 0.05 cos (6x + 0.4t). (8) The equations are solved using an explicit Runge–Kutta time integration scheme and a second-order centered finite-difference discretization in space. We use a uniform Cartesian mesh composed of 160 × 100 points.

The simple structure of the discontinuity allows us to use a simplified approach for re-meshing and interpolation of data. For each x-location, the coordinate y = y

(x) of the shock is detected using the velocity gradient. Then, the domains [−0.5, y

(x)] and [y

(x), 0.5] are mapped, respectively, to [−0.5, 0]

and [0, 0.5] leaving the shock in the mapped domain at the fixed location y = 0 . The x-coordinates of these grid points coincide with the ones used in the simulation of Burgers’ equation. For this reason, the data are interpolated only in thee y-direction using a one-dimensional cubic interpolation. In what follows, we have used 666 snapshots of the two-dimensional field for the POD computations.

The singular values obtained by standard POD and shock-fitted POD are presented in figure 2. As can be seen, the singular values decay rapidly in both cases. Still, the decrease is more pronounced in the case of the shock-fitted POD for which the singular values decay by three orders of magnitude over the first 20 values. This drop is one order of magnitude larger than for the standard POD.

We recognize that the behavior of the singular values is not a true or the

only reflection of the improved performance of the modal decomposition; we

also have to compare the computed POD modes and evaluate their ability to

capture the essential features of the data sequence. A set of selected modes is

0 200 400 600 10

10

10

k σ

standard POD shock-fitting POD

Figure 2: Decay of singular values for the Burgers equation case. Results correspond to computations based on 666 snapshots of the filed.

given in figure 3. The left column shows the result of the standard POD, while the right column displays modes obtained from the shock-fitting POD. The plots in the top row correspond to the zeroth mode which represents the mean value over the snapshot sequence. Already in this mode, traces of the shock motion are visible in the results of the standard POD, while the shock-fitting POD produces a smooth field. In the subsequent modes, the shock is captured well by shock-fitted POD modes while the results of the standard POD include nonphysical structures which cannot be related to any relevant feature of the flow. We also notice that the structures off the moving shock are overwhelmed by the smeared out shock motion in the standard POD, whereas the shock-fitted POD clearly depicts the modal structures on both sides of the interface.

Furthermore, we assess the reconstructive potential of the computed modes by considering a superposition of modes to recover the original fields. In figure 4, reconstructions of the u-component of the data fields using twenty modes are presented. The shock in the reconstructed solution using standard POD modes is rather diffuse, while the shock-fitted POD modes create a much sharper structure with the same number of modes.

3.2. Supersonic flow around a cylinder

With the encouraging results of the previous section, we continue to evaluate

the performance of our shock-fitting POD method on a physical flow case: an

unsteady supersonic flow past a circular cylinder. The characteristic parameters

of the flow are given by the free-stream Mach number of M

= 2.0 and the

Reynolds number, based on the free-stream velocity and cylinder diameter, of

Re = 5000 . The cylinder has a non-dimensional radius r = 0.5. The extent of

the computational domain is x ∈ [−2, 8] and y ∈ [−4, 4]. Symmetry boundary

(a) (b)

(c) (d)

(e) (f)

Figure 3: Comparison of the u-component of standard POD modes, showing the (a) zeroth, (c) first and (e) 16-th mode, and of the shock-fitting POD modes, showing the (b) zeroth, (d) first and (f) 15-th mode, for the Burgers’ equation.

(a) (b)

Figure 4: Comparison of the reconstructed solutions with 20 modes of (a)

standard POD and (b) shock-fitting POD.

conditions are applied at y = ±4.

With a constant inlet condition, the flow is unsteady mainly due to the shedding of vortices in the wake of the cylinder. In order to accentuate the motion of the upstream bow shock, the inlet velocity is modulated in time and space according to

U

(y, t) = U



1 + a exp



− y

2σ

sin(ωt)

 (9)

where U

= 686.229 m s

, a = 0.2, σ

= 0.5 m and ω = 1000 s

. The inlet temperature and pressure are kept constant at 293 K and 1.013 × 10

Pa , respectively.

Solutions to the compressible Navier–Stokes equations are computed with OpenFOAM using a finite-volume discretization based on 721920 cells. To simplify data processing, the flow fields are interpolated onto a 300 × 200 Cartesian grid. Data for grid points inside the cylinder are set to zero. The shocks are detected using the norm of the gradient of the pressure logarithm. A typical snapshot of the flow is shown in figure 5. Plot (a) shows the local Mach number values and in (b) k∇ log pk is plotted in pseudo-colors. The locations of the detected shocks are marked in black in (a) in orange in (b). The flow exhibits a relatively complex shock structure which is driven by and interacts with the vortices shed from cylinder.

The POD computations below are based on 600 snapshots. For simplicity, the surface of the cylinder is treated as an interface, even though its location is imposed as a geometric constraint. For reference, a standard POD is performed from an SVD of the snapshot matrix. The shock-fitting POD is then computed as described above. The number of spatial points considered the computations is identical between the two methods.

We start by comparing the singular values. Figure 6 shows the decay of the singular values of the regular snapshot matrix and the corresponding matrix from the shock-fitted approach. As can be seen, the latter approach results in a slightly faster decay of the singular values, even though the difference to the standard POD is not as pronounced as in the earlier Burgers’ case.

Nonetheless, it implies that fewer modes are necessary to approximate the data with shock-fitted POD modes.

As mentioned before, the decay behavior of the singular values does not convey the full potential of the shock-fitting approach. To draw a closer comparison between the two decomposition methods, we consider the shape of the computed modes. We choose to present the modes using the fluid density (ρ) component. The structure of three selected modes is shown in figures 7 and 8.

The standard POD modes tend to contain most of their energy in spatial feature

linked to the motion of the shocks, while discounting the flow structures away

from the shocks. Figures 7 (a) and (b) correspond to the average fields. For the

case of standard POD, the shocks are far more blurred, while the equivalent

results of the shock-fitted POD exhibit rather sharp discontinuities and reflect

(a)

(b)

Figure 5: Instantaneous flow filed around the cylinder. (a) Flow filed presented in terms of the local Mach number. (b) Pressure-gradient-based shock detection.

The quantity k∇ log pk is used to isolate the structures from the flow field (a).

The locations of the detected shocks are plotted in black in (a) and in orange

in (b).

0 200 400 600 10

10

10

10

10

k σ

standard POD shock-fitting POD

Figure 6: Decay of the singular values of a data matrix sampled from supersonic flow about a circular cylinder. Results are based on 600 snapshots of the two-dimensional flow field.

(a) (b)

Figure 7: Comparison of the ρ-component of zeroth POD modes from (a) the standard POD and (b) the shock-fitting POD. The shock-fitting POD mode has been interpolated onto the base grid built around the average location of the shocks.

the correct shock motion. In figure 8 we present selected higher-order modes to

draw the attention at the representation of dynamic flow features by the two

methods.

(a) (b)

(c) (d)

(e) (f)

Figure 8: Comparison of the ρ-component of (a) the first, (c) the fourth and

(e) the eleventh standard POD modes with the equivalent first (b), fourth (d)

and eighth (f) shock-fitting POD modes. The flow fields for the shock-fitting

POD modes have been interpolated onto the base grid built around the average

location of the shocks.

0 2 4

−1 0 1

·10

x

ρ

standard POD mode 11 shock-fitting POD mode 8

(a)

−2 −1 0

1 1.5 2 2.5

x ρ/ρ

exact standard POD shock-fitting POD

(b)

Figure 9: Slice of density at y = 1.3, in (a) selected modes and (b) reconstructed solutions.

Whereas the order of the computed modes may not coincide between the two computations, we have tried to identify modes with similar structures to allow a fair comparison of the higher-order modes.

Some of the modes from the standard POD, however, do not have any equivalents among the shock-fitted modes. Examples of those with similar structures are given in figures 8 (a), (b), (e) and (f). We observe that the main difference between the two sets of modes lies in the poor capture of the unsteady shock motion by the standard POD, while the new approach succeeds at producing sharp discontinuities in the modes. One exception to this observation is concerned with the region near x = 2.5 which is caused by a complex pattern of intersecting shocks which changes across snapshots.

Outside of the shock regions, modes from the standard and shock-fitting POD

describe rather similar flow information. In figure 8 (c) an example of a standard

POD mode without a counterpart among the shock-fitted modes is given: we

juxtapose the fourth mode from the standard POD with the fourth mode of the

shock-fitted POD (see figure 8 (d)). Minimal commonality can be detected. A

final comparison is made by evaluating the reconstructed flow fields at certain

locations within the computational domain. The density along the line y = 1.3

in the selected modes and in the reconstructed fields using twenty modes is

plotted in figure 9. In subfigure (a) modes with similar structures from the two

sets of solutions are considered (the same as those in figures 8 (e) and (f)). As

can be seen, the POD modes show significant amplitude oscillations around

the shocks that superimpose to approximate the discontinuity. In figure 9 (b)

spurious oscillations around the shocks are visible in reconstruction based on

standard POD, while the shock-fitting POD solutions coincide with the original

data.

4. Summary and conclusions

We have developed a shock-fitting approach to improve the quality of POD analysis of fields with moving discontinuities such as those encountered in compressible flows with shock-shock and shock-boundary-layer interactions.

Standard POD can be seen as a singular value decomposition of a snapshot matrix sampled from fields on a fixed grid. Its use in such cases is problematic due to the difficulty of linear superposition of smooth fields to reconstruct discontinuous features. This issue is caused by the fact that some points of the modes are used to approximate the solution on either side of the discontinuities, as the interfaces sweep across the domain. By broadening the concept to a matrix of snapshots sampled at points moving in sync with the discontinuities, and thus always remaining on the same side of them, the data sampled at each point is rendered continuous across snapshots. A method to track shocks between snapshots and displace solution points continuously has been devised for two-dimensional fields.

A first application to the Burgers’ equation and a more complex application to the compressible Navier-Stokes equations have demonstrated the advantages of this method. In addition to singular values of the snapshot matrix decaying slightly faster, the modes are void of typical oscillations around shocks and focus more directly on flow features that are not linked to the interface motion, but rather with feature in the smooth regions of the flow. Instead of the unsteady discontinuities being distributed across many modes and the full interface motion being represented as a superposition of these modes, the discontinuities are incorporated into a single mode. As a result the shocks in the reconstructed solutions are also sharper and far less oscillatory.

While the shock-fitting approach is showing first and promising results, it suggests various extensions and generalizations for treating modal decoom- positions with moving interfaces. An incorporation of grid locations into the data snapshots should further improve the reconstruction capabilities of the method and result in an even more faithful representation of the flow dynamics, especially for sizeable interface motions. Furthermore, a formulation that allows for topology changes in the layout of the interface-bounded subdomains is desirable in order to address even more complex shock and interface motions.

These extensions will be topic of future research efforts.

Acknowledgement

This project has been funded through the European Union’s Horizon 2020

research and innovation programme under the Marie Skłodowska-Curie grant

agreement No 675008 and the European Research Council under grant agree-

ment 694452-TRANSEP-ERC-2015-AdG. Computer resources provided by the

Swedish National Infrastructure for Computing (SNIC) at the Center for High

Performance Computing (PDC) at the Royal Institute of Technology (Stock-

holm), the High Performance Computing Center North (HPC2N) at Umeå

University, and the National Supercomputer Centre at Linköping University are gratefully acknowledged.

References

Ando, S. 2000 Image field categorization and edge/corner detection from gradient covariance. IEEE Transactions on Pattern Analysis and Machine Intelligence

Balajewicz, M. & Farhat, C. 2014 Reduction of nonlinear embedded boundary models for problems with evolving interfaces. Journal of Computational Physics

Iuliano, E. & Quagliarella, D. 2013 Proper orthogonal decomposition, surrogate modelling and evolutionary optimization in aerodynamic design. Computers &

Fluids 84, 327–350.

Li, J. & Zhang, W. 2016 The performance of proper orthogonal decomposition in discontinuous flows. Theoretical and Applied Mechanics Letters 6 (5), 236–243.

Lu, H. & Tartakovsky, D. M. 2020 Lagrangian dynamic mode decomposition for construction of reduced-order models of advection-dominated phenomena.

Journal of Computational Physics 407, 109229.

Lucia, D., King, P., Oxley, M. & Beran, P. 2001 Reduced order modeling for a one-dimensional nozzle flow with moving shocks. In 15th AIAA Computational

Lucia, D. J., King, P. I. & Beran, P. S. 2003 Reduced order modeling of a two-dimensional flow with moving shocks. Computers & Fluids 32 (7), 917–938.

Lumley, J. L. 1970 Stochastic tools in turbulence.

Mendez, M. A., Balabane, M. & Buchlin, J.-M. 2019 Multi-scale proper or- thogonal decomposition of complex fluid flows. Journal of Fluid Mechanics 870, 988–1036.

Persson, P.-O. & Strang, G. 2004 A simple mesh generator in MATLAB. SIAM

Rowley, C. W. & Dawson, S. T. M. 2017 Model reduction for flow analysis and control. Annual Review of Fluid Mechanics 49 (1), 387–417.

Sirovich, L. 1987 Turbulence and the dynamics of coherent structures. i. coherent structures. Quarterly of Applied Mathematics 45 (3), 561–571.

Towne, A., Schmidt, O. T. & Colonius, T. 2018 Spectral proper orthogonal

decomposition and its relationship to dynamic mode decomposition and resolvent

analysis. Journal of Fluid Mechanics 847, 821–867.