Overcoming slowly decaying Kolmogorov n-width by transport maps: application to model order reduction of fluid dynamics and fluid–structure interaction problems

by   Monica Nonino, et al.

In this work we focus on reduced order modelling for problems for which the resulting reduced basis spaces show a slow decay of the Kolmogorov n-width, or, in practical calculations, its computational surrogate given by the magnitude of the eigenvalues returned by a proper orthogonal decomposition on the solution manifold. In particular, we employ an additional preprocessing during the offline phase of the reduced basis method, in order to obtain smaller reduced basis spaces. Such preprocessing is based on the composition of the snapshots with a transport map, that is a family of smooth and invertible mappings that map the physical domain of the problem into itself. Two test cases are considered: a fluid moving in a domain with deforming walls, and a fluid past a rotating cylinder. Comparison between the results of the novel offline stage and the standard one is presented.



page 14

page 17


A learning-based projection method for model order reduction of transport problems

The Kolmogorov n-width of the solution manifolds of transport-dominated ...

A monolithic and a partitioned Reduced Basis Method for Fluid-Structure Interaction problems

The aim of this work is to present a brief report concerning the various...

Model Reduction for Advection Dominated Hyperbolic Problems in an ALE Framework: Offline and Online Phases

Model order reduction (MOR) techniques have always struggled in compress...

Model order reduction for bifurcating phenomena in Fluid-Structure Interaction problems

This work explores the development and the analysis of an efficient redu...

A Deep Learning approach to Reduced Order Modelling of Parameter Dependent Partial Differential Equations

Within the framework of parameter dependent PDEs, we develop a construct...

Finite strain homogenization using a reduced basis and efficient sampling

The computational homogenization of hyperelastic solids in the geometric...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The reduced basis method [27, 25, 48, 10, 31, 26, 39]

is a powerful tool when requiring fast simulations of parametrized partial differential equations (PDEs): its efficiency relies on the possibility to construct an approximation of the solution for any value of the parameter in the span of a few basis functions, which are computed in the (expensive) offline phase. Despite its capability has been acclaimed in a large variety of situations, model reduction of advection dominated (or even hyperbolic) problems is still a challenging task

[41, 23, 19]. It has therefore become clear that a modification of the way the reduced basis method works is necessary, especially in order to be able to obtain a small basis set also in these more challenging situations.

Assume that is the physical domain of the problem of interest. Let , , be a compact set, and denote by the parameter. Furthermore, let be the time, for some . In the following the symbol will either stand for in parametrized stationary problems, in non-parametric unsteady problems, or the pair in parametrized unsteady problems. For any value of , we seek the solution , , to the following PDE:


where is a nonlinear operator working on functions defined over , and , or , respectively, in the three aforementioned cases. Specific examples of problems of interests in computational fluid dynamics and fluid-structure interaction will be introduced throughout the paper.

Let be a component of , and be the solution manifold, embedded in some normed linear space , defined as:

The choice of a component-wise solution manifold (rather than just one solution manifold associated to the vector

) will be motivated in section 2. One fundamental assumption of the reduced basis method is that can be approximated in an accurate way by a sequence of finite dimensional spaces: any element of can be recovered using a linear combination of solutions of (1), which are computed only once and for all. The mathematical entity that incorporates this concept is the Kolmogorov -width of , defined as


where is any linear subspace of of dimension . The Kolmogorov -width tells us the entity of the error that we commit by approximating any element of with an element of a linear space . Therefore the faster decays as we let grow, the greater possibility we have to build a good linear approximation space of of low dimension. In the majority of the problems there is no explicit analytic formula for , yet there are some situations where we can compute good bounds on the Kolmogorov -width[36, 17]. In general, since it is very difficult to provide such bounds on , we can only hope that the

-width of the solution manifold is small. A heuristic way to check that this hypothesis is satisfied by the problem of interest is to run a proper orthogonal decomposition (POD) on a set of snapshots, and check the rate of decay of the eigenvalues

returned by the POD: if the decay fast, then we can expect to be able to build a good low dimensional linear approximation space for . This assumption often fails in transport-dominated problems, which show a very slow decay of the eigenvalues of the POD, and thus the inability of the reduced basis method to reconstruct any element of by using a small number of basis functions.

A growing number of works which focus on constructing alternative (nonlinear) model order reduction techniques, that can be effectively applied to transport-dominated problems, has appeared in recent years. Abgrall et al. [3] show an norm minimization technique, to be used for the approximation of nonlinear hyperbolic equations, without anyway curing the problem of high dimensional solution manifolds. Approximated Lax pairs have been employed for model reduction of nonlinear problems arising in cardiac electrophysiology[21, 22] . Gerbeau et al [21] construct a reduced basis space at each time step: the reduced basis are the modes of a Schroedinger operator where the potential is the solution at the previous time step. The drawback of this procedure is the increase of the number of basis functions as the approximation accuracy requirement is increased. Nonlinear model reduction techniques in metric spaces are proposed by Ehrlacher et al. [19]

, relying either on a tangent principal component analysis or barycentric greedy algorithm. Adaptivity is also employed to overcome difficulties arising in model reduction of transport dominated flows, see Carlberg

[15] for the use of localized reduced bases or Peherstorfer [42]. the procedure to update the basis set. Several other different techniques rely on a composition with a suitably defined transport map (on which we will focus on in the rest of the work). Such map may be obtained as the solution of a Monge-Kantorovich optimal mass transport problem [28, 9], or provided analytically in simple cases as the one-dimensional problem considered by Cagniart et al. [14]. More complicated configurations, possibly including shocks, can still be handled relying on transport maps based on more advanced shape parametrization maps [13]. Extension to multiple transport phenomena are also possible, most notably by means of the shifted POD [45]

, as well as (possibly interacting) shocks by recent extensions of the transported snapshots interpolation

[59, 58, 60] and transported snapshots model reduction methods [38]. While in simple cases the underlying transport maps are provided by the user, registration techniques can be employed to automate their selection in more complex geometrical configurations [55]. Application of the use of transport maps for model reduction based on an embedded high fidelity method is shown by Karatzas et al. [30]. An alternative is the freezing method [40, 47, 11], in which the key tool is the identification of a Lie group acting on a frozen solution component. However, to properly develop and analyze the proposed methodology for complex problems in fluid dynamics, very involved mathematical tools and settings are needed which as of now hinders its applicability in a broader setting. This consideration, together with previous observations on other techniques, makes it clear that there is the need for a lighter, simpler and more natural framework. A methodology that satisfies these requirements and that is based on the definition of some transport maps has been introduced and applied to some toy problems by Cagniart [12] and Cagniart et al.[14]. The goal of this work is to present an application of model reduction based on transported snapshots for problems in fluid dynamics and fluid-structure interaction, focusing in particular on the comparison between results of the standard offline phase and the novel one which relies on transport. As (especially in fluid dynamics problems) further challenges are present during the online stage, such as keeping into account stabilization [57, 4] or turbulence modelling [16, 49], we limit our exposition to the offline stage; future research work will extend the results presented here to the online stage. The article is organized as follows: in Section we present the nonlinear model reduction approach by transport maps, and we explain how it helps to overcome the problem of slow decay of the Kolmogorov -width. In Section we present a first example of application: a time dependent CFD problem, without parameters. In Section we study a fluid–structure interaction problem, where the solution behaves like a travelling wave, thus featuring a slow decay of the Kolmogorov -width of the problem; also in this case, the problem is time dependent and non-parametric. In Section , we further extend the CFD test case introduced in Section by adding a physical parameter, namely the Reynolds number. Conclusions follow in section 6.

2 Nonlinear model reduction by transport maps

In this section we summarize the nonlinear approach that we will apply in the forthcoming sections to fluid dynamics and fluid-structure interaction problems. We will closely follow the presentation and the notation introduced by Cagniart et al. [14] for a parametrized viscous Burgers equation. The idea proposed therein is to “pre-process” the solution manifold by a composition with a map in a family of smooth and invertible mappings

Maps in the family are parametrized by the same appearing in (1), and are essentially problem-specific. For instance, in the work by Cagniart et al. [14], a family of translations was employed for the viscous Burgers equation. Specific choices of will be discussed alongside the problem of interest. We then introduce the “pre-processed” solution manifold

Assuming that is carefully chosen, has a smaller Kolmogorov -width than . We remark here that one may well choose different families for different components of the solution ; indeed, it is often the case, especially in multiphysics problems, that the qualitative evolution in (which usually suggests the choice of the family itself) behaves differently for different components of the solution.

The practical realization of this preprocessing procedure is incorporated in the offline phase. Given a discrete training set , we compute each solution component associated to any . The discrete approximations of the corresponding standard and preprocessed solution manifolds

provide snapshots for a compression by a POD. The compression is here applied to both and to provide a comparison between the standard offline phase and one with preprocessing, but in practical computations one would neglect the compression of , as it is understood that would result in a POD basis set of lower dimension.
As anticipated in the previous section, we will only focus on the offline part of the reduced order method; efficient evaluation of the online solution in a Galerkin projection setting is a future development of this work; nevertheless we provide for completeness some details on how the online phase of this new reduction method should be carried out. Let us assume that an accurate approximation of the solution component at time is known, as a linear combination of our basis functions :

In order to recover as an expansion:

the idea proposed by Cagniart et al.[14] is to iterate between the search for the reduced coordinates and the search for a suitable parameter . This procedure is carried out by means of a minimization problem of the norm of the residual evaluated at .

In the next two sections we will see two applications of this preprocessing procedure: the first problem we are going to study is a computational fluid dynamics (CFD) problem, and the second problem is a fluid–structure interaction problem. These two applications are quite different for what concerns the physics behind them (the latter is a coupled multiphysics problem, the former is not), but they both feature a slow decay of the eigenvalues returned by running a POD on . In the first case is characterized by the change of the direction of propagation of a vortex close to a rotating cylinder, a feature which is difficult to reproduce with a small number of basis functions. The second problem is a transport dominated problem, where the solution behaves like a wave travelling in the domain. Also in this case the travelling wave would be hard to be reconstructed with a standard model order reduction technique.

3 A CFD test case

In this Section we are going to show the results that we obtain with the preprocessing procedure on a CFD test case, namely a fluid past a rotating cylinder. The problem described in the following is inspired by the one presented in Cagniart’s thesis [12], although our formulation is slightly different: in the reference, the direction of the fluid velocity at the inlet boundary is changing in time, and the cylinder is kept fixed, whereas on the contrary in our problem the inlet direction is constant, but the cylinder is rotating.

3.1 Problem formulation

In this problem we have a flow past a rotating cylinder; this situation, with or without the rotation feature, is quite interesting and has been indeed studied in a large number of papers, both in the incompressible case [8, 50, 29, 51] and in the compressible regime [56]. It is well known that if the Reynolds number is greater than , then a vortex shedding phenomenon occurs [29]. When, in addition, the cylinder is rotating, it might happen that we see a change in the direction of propagation of this vortex; this phenomenon is quite complicated and is strictly related to a variety of different physical quantities. In fact, it is related not only to , but also to the cylinder rotation rate , which is defined as:

where is the diameter of the cylinder, is the angular velocity of a point on the surface of the cylinder and is the oncoming free stream velocity. For vortex shedding occurs, where is a critical value that is a function of . We refer to Stojković et al. [50] for different values of related to different values of the Reynolds number.

To model our problem we use incompressible Navier-Stokes equation. Figure 1 shows the physical domain for our problem. We do not take into account any physical parameter. Our problem reads as follows: for any time find and such that:

We impose homogeneous Neumann conditions on top and bottom walls and , and also on the right boundary , as all these boundaries are considered as outlets due to the fact that the vortex sheet rotates alongside with the cylinder. We impose a Dirichlet condition on , where is a fixed horizontal inflow tabulated in Table 1.

Furthermore, is the tangential velocity at the surface of the rotating cylinder. At the beginning of the simulation the cylinder is not moving, and it stays still until a fully developed vortex shedding phenomenon is reached; after that the cylinder starts to rotate counterclockwise, first with a constant acceleration , then it keeps on rotating with a constant angular velocity. Thus, denoting by be the angular velocity of the cylinder, we set


Once we have the angular velocity, we can compute thanks to the relation:

where is the radius of the cylinder, and assuming to be tangent to .

Figure 1: Physical domain of the CFD test case.

As far as the problem formulation concerns, let us mention the fact that, in preparation of a forthcoming online stage, we adopted a supremizer enrichment technique, meaning that we decided to enrich the space of the fluid velocity snapshots with some supremizers snapshots . The supremizer snapshots are needed at the reduced order level in order to have a more stable approximation of the fluid pressure. For further details on the formulation of the supremizer in the POD framework of parametrized fluid flows we refer to Ballarin et al.[5].

Data Value Data Value
s s
Table 1: Problem data for the CFD test case.
Figure 2: Decay of the eigenvalues for the fluid velocity in the simulation with a rotating cylinder.

In Table 1 we can find the problem data that we used in our simulation. After we obtain a fully developed Karman vortex, and after the cylinder starts to rotate, we have a noticeable change in the direction of propagation of the vortex. This change of direction may hinder the representation of the solution by a small number of basis functions, and this expectation is confirmed by running a POD on the fluid velocity snapshots collection, as we can see from Figure 2.

3.2 Preprocessing step

We will focus on the preprocessing of the fluid velocity , since it is more straightforward to visualize the direction of propagation of the Karman vortex and hence understand the idea beneath the deformation map.

Figure 3: Subdivision of the physical domain into (red) and (blue).

First of all, let us notice from Figure 3 that, when building the mesh, we have defined a fictitious subdomain (in red). It is very important to choose the radius of the circular subdomain in such a way that is entirely contained in the physical domain and yet is able to capture all the complex behaviour of the solution that is strictly related to the rotation of the cylinder. In fact the blue subdomain in Figure 3 deals with the small vortexes that originate in the wake of the cylinder before it starts to rotate: these vortexes propagate in the domain and their behaviour is not heavily affected by the rotation of the cylinder, and therefore we should be able to reproduce their behaviour with a coarse set of basis functions. Therefore as the solution features the most interesting phenomena in a neighborhood of the cylinder , we are not going to preprocess the entire snapshot, but we are going to focus instead just on its restriction on . This procedure is more adapted, in the sense that we are not considering the solution manifold as a global entity, but we think of it as made of two manifolds, and : this division is based on the division of the physical domain of interest into two subdomains and . Therefore:


Being far from the rotating cylinder, and therefore far from the rotating phenomenon, we expect to be a better behaved solution manifold with respect to . This means that has a small Kolmogorov -width, and does not need any preprocessing. On the contrary we will focus on , which will have a slowly decaying Kolmogorov -width. Before going any further, let us remark that the subdivision of in and is strictly dependent on the subdivision of , therefore the fictitious cylinder has to be chosen wisely, in such a way that is able to capture most of the rotating phenomenon. In particular, in our simulations, we have chosen to be a cylinder of radius , times larger than the radius of the physical cylinder.

For the preprocessing step therefore we first restrict the snapshots to the subdomain . Then, we want to build a (one-parameter) family of smooth and invertible maps

where is the parameter identifying each map. In our case it is natural to choose, at each time , to be the angle spanned by the direction of propagation of the vortex (obtained through a postprocessing of the solution ) and the horizontal axis. We choose the following preprocessing map:


The idea behind (6) is pretty straightforward: at each timestep of our simulation, we compute how much the vortex has changed its direction of propagation, and we therefore find . After that, in the preprocessing step we take all the snapshots , we restrict them to the subdomain , and then we rotate them back to a horizontal direction of propagation of the vortex with .

Figure 4: POD comparison before (black) and after (magenta) the preprocessing of .

In Figure 4 we see the improvements (in terms of POD eigenvalues decay) that we get by applying the rotation to the snapshots: as we can see, after the preprocessing we obtain an improvement e.g. of almost orders of magnitude comparing the -th eigenvalue of standard and preprocessed manifolds.

4 A multiphysics problem

We are now interested in applying the preprocessing procedure to a fluid-structure interaction problem, whose solution exhibits a transport dominated behaviour. The problem formulation features an Arbitrary Lagrangian Eulerian (ALE) formulation, see Richter [46] for more details. An extensive explanation on how to treat such a coupled problem in a reduced order modelling setting can be found in Ballarin et al.[6] (monolithic approach) and in Ballarin et al.[7] (partitioned scheme). For further details on reduced order models and applications of FSI problems we refer to Bertagna et al.[10], Lassila et al. [32, 33] and Colciago [18].

4.1 Problem formulation

We have a two dimensional rectangle of height and length , filled with a Newtonian fluid. The structure is made by the top and the bottom of the rectangle, which are considered to be deformable, and their thickness is negligible with respect to the height of the rectangle. Since the structure is thin, it is described by a one dimensional model. We further assume that the displacement of the walls in the horizontal direction is negligible, and hence the structure presents only vertical motion: the behaviour of the compliant walls is therefore described by the generalized string model [43, 44]. We want to describe the behaviour of the solution (and the domain itself) in the time interval .

Let denote the structure, and let be the fluid domain at time . The fluid equations are formulated, in this particular case, on a moving fluid domain . Let be the fluid reference domain: for convenience we take (the blue fluid domain in Figure 5).
The map that maps to at time is :

where denotes the vertical component of the displacement (or deformation) of the fluid mesh, as such displacement in the horizontal direction is negligible, since the horizontal displacement of the structure is negligible.
Thanks to the introduction of we can therefore map the fluid equations back to the reference domain , thus introducing additional terms to the classical Navier–Stokes equations; this results in an Arbitrary Lagrangian Eulerian (ALE) formulation [46, 24] of our original problem. Since we are considering the case of very small deformations, we will define as an harmonic extension of the structure displacement :


There are alternative ways of defining the fluid displacement [46] (pseudo-elasticity, bi-harmonic).

Let now be the gradient of , and let be the Jacobian. The fluid-structure interaction problem reads: find fluid velocity , fluid pressure and structure displacement such that:


Here is the fluid viscosity, is the structure viscosity, is the structure height (or thickness), and are the structure constitutive parameters.

is the Cauchy stress tensor for the fluid, and is defined as:

being the kinematic viscosity of the fluid and being the identity matrix. Finally, is the unit vector .

Figure 5: Physical domain: fluid subdomain (blue) and structure subdomain (red). The fluid-structure interface coincides with the structure in our case.

Since it is a fluid-structure interaction problem, we need some coupling conditions. Let us denote the fluid structure interface at time ; we have:


The condition on the continuity of the displacement is a geometric condition, which stems from the fact that we do not want the fluid and the solid subdomain to overlap, and which we already employed in the definition of (7). As far as the second condition concerns, let us remark a very important fact: the quantity does not coincide with the fluid velocity , as denotes the velocity of the mesh, and therefore it is not the physical velocity of the fluid particles. For this reason, in general situations in FSI like our problem, we ask for a continuity condition for the physical fluid velocity and the mesh velocity at the fluid–structure interface, that is exactly where the mesh velocity coincides with the structure velocity (see (7)). Finally, there is a third coupling condition, which does not appear together with the others, since it is included in the right hand side of the generalized string equation (8); this condition arises from the classical action-reaction principle.

Together with the coupling conditions we also have to give some boundary and some initial conditions. For the latter we assume that at time the system is at rest; the boundary conditions can be summarized in the following system:


where n is the outward normal. The third condition says that the structure is fixed at its extremities.
As far as the problem formulation concerns, let us remark that we adopted a supremizer enrichment technique also in this multiphysics test case, always to obtain, in the POD framework, a set of basis functions that allow for a stable approximation of the fluid pressure.

4.2 Transport dominated FSI problem

Data Value Data Value
g/cm dyn/cm
Table 2: Problem data for the FSI test case

Problem data used for the simulation of our test case can be found in Table 2: corresponding values are taken from the numerical results presented by Sy and Murea[54, 37].

Figure 6: Fluid pressure behaviour: the solution is pictured here at time , and . The peak of the wave is propagating into the domain, creating a transport phenomena.
Figure 7: Fluid displacement behaviour: again the solution is pictured at time , and . The peak of the wave is still very small at the beginning, it grows for some time and then it starts to propagate.
Figure 8:

Decay of the singular values for the POD on the fluid pressure (red, top left), fluid displacement (blue, top right) and fluid supremizer (black, bottom center).

The behaviour of the fluid pressure and of the extended displacement is shown in Figure 6 and Figure 7 respectively. Our problem is transport dominated: if we look at Figure 6 for example, the change in time of the position of the peak of the pressure wave will be a difficult feature to capture at the reduced order level with just a few modes. This expectation is finally confirmed at the numerical level, as we can see in Figure 8, by running a POD on , , and also . Therefore we rely on a transformation on the set of solutions, in order to compensate this transport phenomenon.

In order to make the following exposition more clear and easy to read, from now on we focus only on a particular component of the solution of our problem, namely the fluid pressure. It is anyway important to keep in mind that, based on our simulation, all the components of the solution to the FSI problem are subject to a transport phenomenon, and hence every consideration that we are going to make on could be easily applied to any of the other components of the solution.

4.3 Preprocessing step

Let us see more in detail how to apply this preprocessing procedure to the fluid pressure solution manifold . Figure 6 shows how the peak of the pressure is transported in the domain. We would like to align the peaks at all time steps in a reference configuration; in this way, we obtain a set of snapshots where the pressure wave is not moving at all. In this case therefore a low number of modes will be sufficient to give a good representation of the situation.

Starting from this observation we build a one parameter family of mappings

such that for every in , the peak of is located not at its original position, but is instead moved to the middle of the domain. In this way the new snapshots will all have the peak located at the exact same position, meaning the middle of the domain.

When building the map , another aspect to which we should pay great attention is the boundary conditions. Since we are not working in a periodic setting, we want to make sure that the preprocessed snapshots satisfy the same boundary conditions as the original snapshots. An easy way to make sure that these requirement is satisfied is to keep the points in and the points in fixed. After some computations it follows that one possible map is the following:

where is the abscissa of the position of the peak of the wave at time . We assume that the abscissa of the points on the inlet boundary is and the length of the domain is ; in addition, the position in which we are moving the peak of the wave at every time is exactly in the center of the domain. Let us remark that the map is just a stretching in the horizontal direction: this is due to the fact that we do not have any transport phenomena in the vertical direction, and hence there is no need for a transformation in the axis. So, with an abuse of notation, we can think of as .

Figure 9: Comparison on the decay of the eigenvalues for the POD on the original manifolds and on the preprocessed manifolds for the solutions components (top left), (top right) and (bottom center).
Figure 10: Original snapshots for at time , and (a), and corresponding preprocessed snapshots (b).

Our family of mappings is a one-parameter family, therefore to identify the stretching map we have to compute the parameter , which depends on time . First of all we compute all the snapshots , where , and . Once we have the snapshots, we locate the abscissa of the peak of the wave ; we then preprocess with and obtain a new snapshot ; see Figure 10 for a comparison between the original snapshots, and the preprocessed ones. The snapshot at time has not been preprocessed, and the reason for this is as follows: recall the expression of from Table 2; as we can see, represents a given pulse at the inlet boundary, which is nonzero up to time . So, up to time the pulse makes the wave grow; when the pulse is null, the wave starts to propagate into the domain. As previously mentioned, we are looking for a stretching map that moves the peak of the wave in a particular location of our choice, and at the same time keeps the points on the inlet boundary and on the outlet boundary fixed. However, in our simulation the peak of the wave is located exactly or very close to the inlet boundary for the first seconds: the transport phenomena is not immediate; this implies that, until time we do not need to apply any preprocessing step. At the reduced level this translates into the fact that, up to we adopt a standard reduction technique, using a standard reduced basis obtained with a standard POD on the first snapshots. From the snapshot to the last one, on the other hand, we first perform the preprocessing step, and then perform a POD.

In Figure 9 we can compare the decay of the eigenvalues for the POD on the pressure with and without this preprocessing procedure. As we can see, with the preprocessing technique we do actually get an improvement in the decay of the eigenvalues, and in fact with less than modes we reach a level of , which is one order of magnitude less than the one we get with less than modes in the standard case. We get the same results for the extended displacement .

5 CFD test case with a physical parameter

So far we have considered test cases where the only parameter was time, nevertheless we are interested in investigating what happens if we add another physical parameter to the problem. We then go back to the CFD test case presented in section , but now the Reynolds number will be considered as a parameter. The fluid velocity solution manifold will be defined as:

where now . Of course in this case we have to pay attention not only to the Reynolds number, but also to the rotation rate of the cylinder, since it is strictly related to the development of a vortex shedding phenomenon. For our particular CFD test case we choose the following parameter range:

and we choose . We discretize the parameter space, choosing a set of parameter samplings . For our problem we choose a Lagrange distribution sampling, i.e:

After we have obtained a set of snapshots for each parameter in the parameter sampling set, we observe that the change in the Reynolds number leads to changes in the behaviour of the fluid velocity, with the vortex shedding that may occur earlier or later, but in all the situations we see that after a while, due to the rotation of the cylinder, the direction of propagation of the vortex changes.

5.1 POD-Greedy

With the addition of a physical parameter, the exploration strategy will be carried out in a different way with respect to a standard POD on the set of snapshots; we are going to explore the parameter space with a pseudo-Greedy algorithm, and we are going to explore in time with a POD on the set of snapshots corresponding to each parameter selected by the pseudo-Greedy strategy. First of all we discretize the parameter space in order to obtain a sampling set of cardinality of our choice; in order to do so, we choose a Lagrange distribution sampling, i.e:

Once we have the parameter sampling set , we compute the truth solution for each one of these parameters. Afterwards, the POD is applied in the following way:

  1. for , we run a standard POD on the corresponding snapshots;

  2. we now have at hand a set of reduced basis ;

  3. for , , we orthogonalize each snapshot with respect to the linear space ;

  4. we then run a POD on the set of orthogonalized snapshots, and add the resulting basis functions to the already existing set of basis functions.

We remark that we are not actually using a proper Greedy algorithm, because we are choosing a priori the parameter sampling set; the reason for this “pseudo”-Greedy is that it is beyond the scope of this work to focus on an error estimator to be used for the Greedy algorithm, and also a simple Lagrange distribution will be sufficient to have an insight on what is going on before and after the preprocessing step. We also remark that different possibilities for the POD are possible: the one we are using here is based on an orthogonalization step, where we are getting rid of the superfluous information before pursuing a POD

[26]. Another possibility would be to run first a standard POD on the snapshots, for every parameter, and then, at the end, run another standard POD on the set of obtained reduced basis, again to get rid of the superfluous information [39]. We do not present here any result showing the rate of decay of the greatest eigenvalue returned by running a POD with orthogonalization on the set of snapshots corresponding to each parameter of the sampling set. The reason of this choice is the fact that, as we previously remarked, we did not actually use a Greedy procedure to select the parameters in the parameter space; in addition, a plot showing the decay of the first eigenvalues of each POD gives a good idea of how well the orthogonalization procedure is doing. The use of a proper Greedy algorithm based on a reliable error estimator and the proof that the orthogonalization procedure in the POD produces good advantages are outside the scope of this paper.

5.2 Preprocessing step

The preprocessing procedure is carried out, for every value of the parameter in the discretized parameter space, in the same way it was carried out in the case with no physical parameter: , and for every timestep we compute , which is the angle spanned by the direction of propagation of the vortex and the horizontal axis. Then the deformation map is defined as:

Also in this case the idea is to rotate back to a horizontal direction of propagation all the vortexes. We remark that the preprocessing procedure is carried out on the snapshots restricted to the domain , exactly as we did for the problem in Section .

Figure 11: Comparison of the rate of decay of the eigenvalues on the original velocity solution manifold and on the preprocessed solution manifold for (top left), (top right) and for (bottom).

Figure 11 shows the results that we obtain for three different values of the Reynolds number. As we can see, there is an improvement in the rate of decay of the eigenvalues, with results showing a difference of one order of magnitude for (right column) with just modes.

6 Conclusions and future perspectives

In our work we used a preprocessing of the snapshots during the offline stage of the reduced basis method to improve the rate of decay of the Kolmogorov -width of the solution manifold of the problem of interest. We focused on two different test cases: a fluid dynamics problem where a Karman vortex develops in the wake of a cylinder and where the vortex changes its direction of propagation due to the rotation of the cylinder, and a multiphysics problem, where the solution behaves like a travelling wave. The preprocessing procedure is different according to the problem at hand, since the definition of the deformation maps changes according to the behaviour of the solution and to the boundary conditions used. The method performs very well for one dimensional problems [12]. In this work we focused on two dimensional problems in a non-periodic setting. We can say that by adopting this preprocessing procedure of the set of snapshots, the saves from a computational point of view in terms of the dimension of the set of basis functions needed to reach a certain approximation accuracy is evident. The results that we obtained are promising: in the coupled problem, to reach a magnitude of for the eigenvalues related to the POD on the pressure, we need less modes with respect to the standard situation with no preprocessing, and we can say the same about other components (fluid displacement, fluid velocity), thus lowering the dimension of the system to be solved in the online phase of the reduced method of at least . We obtained promising results also in the fluid dynamics test case. In the non parametric problem, results for the fluid velocity show that the eigenvalues after the preprocessing decay with almost two orders of magnitude faster than the standard case: to reach a magnitude of in the standard case we need almost modes whereas in the preprocessed case we need modes. In the parametric case we analyzed the results for three different values of the physical parameter: all the three cases showed good results.
For future perspectives, we do believe that there are several possibilities to further develop and improve this preprocessing technique. In our work, both in the FSI problem and in the CFD one, we relied on our a priori knowledge of the behaviour of the solution in order to compute exactly the parameter identifying the deformation map at time : we computed exactly the position of the peak of the wave in the multiphysics problem, and we computed exactly the angle spanned by the direction of propagation of the vortex and the horizontal axis in the fluid dynamic problem. There is thus a large space to develop methods to find a good parameter , instead of using the exact one, that ensure a good performance of the preprocessing technique and that are computationally feasible. In this work we did not go in the details of the online phase of the reduced method with the preprocessing technique, as we were mainly interested in studying the performance of the new method in different situations, and comparing the rate of decay of the eigenvalues to quantify the improvements that we obtain. Future research work will include: the efficient evaluation of the online phase, and, in addition, application of this reduced order method coupled with the preprocessing technique in the framework of inverse problems[35, 20, 53, 52, 61].


We acknowledge the support by European Union Funding for Research and Innovation – Horizon 2020 Program – in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” (PI Prof. Gianluigi Rozza). We also acknowledge the INDAM-GNCS project “Advanced intrusive and non-intrusive model order reduction techniques and applications”, . The computations in this work have been performed with multiphenics [2], which is a library that aims at providing tools in FEniCS [34] for an easy prototyping of multiphysics problems on conforming meshes, and RBniCS [1], which is an implementation of several reduced order modelling techniques in FEniCS. We acknowledge developers and contributors to these three libraries.


  • [1] RBniCS - reduced order modelling in FEniCS. http://mathlab.sissa.it/rbnics, 2015.
  • [2] multiphenics - easy prototyping of multiphysics problems in FEniCS. http://mathlab.sissa.it/multiphenics, 2016.
  • [3] R. Abgrall, D. Amsallem, and R. Crisovan. Robust model reduction by -norm minimization and approximation via dictionaries: application to nonlinear hyperbolic problems. Advanced Modeling and Simulation in Engineering Sciences, 3(1):1, 2016.
  • [4] M. Azaïez, T. C. Rebollo, and S. Rubino. A cure for instabilities due to advection-dominance in pod solution to advection-diffusion-reaction equations. arXiv:1907.05614, 2019.
  • [5] F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza. Supremizer stabilization of POD–Galerkin approximation of parametrized steady incompressible Navier–Stokes equations. International Journal for Numerical Methods in Engineering, 102(5):1136–1161, 2015.
  • [6] F. Ballarin and G. Rozza. POD–Galerkin monolithic reduced order models for parametrized fluid-structure interaction problems. International Journal for Numerical Methods in Fluids, 82(12):1010–1034, 2016.
  • [7] F. Ballarin, G. Rozza, and Y. Maday. volume 17 of MS&A, chapter Reduced-order semi-implicit schemes for fluid-structure interaction problems, pages 149–167. Springer International Publishing, Cham, 2017.
  • [8] S. Behara and S. Mittal. Flow past a circular cylinder at low reynolds number: Oblique vortex shedding. Physics of Fluids, 22(5):054102, 2010.
  • [9] F. Bernard, A. Iollo, and S. Riffaud. Reduced-order model for the BGK equation based on POD and optimal transport. Journal of Computational Physics, 373:545–570, 2018.
  • [10] L. Bertagna and A. Veneziani. A model reduction approach for the variational estimation of vascular compliance by solving an inverse fluid–structure interaction problem. Inverse Problems, 30(5):055006, 2014.
  • [11] W. Beyn and V. Thümmler. Freezing solutions of equivariant evolution equations. SIAM Journal on Applied Dynamical Systems, 3(2):85–116, 2004.
  • [12] N. Cagniart. A few non linear approaches in model order reduction. 2018.
  • [13] N. Cagniart, R. Crisovan, Y. Maday, and R. Abgrall. Model order reduction for hyperbolic problems: a new framework. hal:01583224, 2017.
  • [14] N. Cagniart, Y. Maday, and B. Stamm. Model Order Reduction for Problems with Large Convection Effects, pages 131–150. Springer International Publishing, 2019.
  • [15] K. Carlberg. Adaptive -refinement for reduced-order models. International Journal for Numerical Methods in Engineering, 102(5):1192–1210, 2015.
  • [16] T. Chacón Rebollo, E. Delgado Ávila, M. Gómez Mármol, F. Ballarin, and G. Rozza. On a certified smagorinsky reduced basis turbulence model. SIAM Journal on Numerical Analysis, 55(6):3047–3067, 2017.
  • [17] A. Cohen and R. DeVore. Kolmogorov widths under holomorfic mappings. IMA Journal of Numerical Analysis, 36(1):1–12, 2015.
  • [18] C. M. Colciago and S. Deparis. Reduced order models for fluid-structure interaction problems with applications in haemodynamics. page arXiv:1801.06127, 2018.
  • [19] V. Ehrlacher, D. Lombardi, O. Mula, and F.-X. Vialard. Nonlinear model reduction on metric spaces. application to one-dimensional conservative pdes in wasserstein spaces. arXiv:1909.06626, 2019.
  • [20] F. Galarce, J.-F. Gerbeau, D. Lombardi, and O. Mula. State estimation with nonlinear reduced models. application to the reconstruction of blood flows with doppler ultrasound images. https://arxiv.org/abs/1904.13367, 2019.
  • [21] J.-F. Gerbeau and D. Lombardi. Approximated lax pairs for the reduced order integration of nonlinear evolution equations. Journal of Computational Physics, 265:246–269, 2014.
  • [22] J.-F. Gerbeau, D. Lombardi, and E. Schenone. Reduced order model in cardiac electrophysiology with approximated lax pairs. Advances in Computational Mathematics, 41(5):1103–1130, 2015.
  • [23] C. Greif and K. Urban. Decay of the kolmogorov -width for wave problems. arXiv:1903.08488, 2019.
  • [24] M. Gurtin. An Introduction to Continuum Mechanics. Mathematics in Science and Engineering. Elsevier Science, 1982.
  • [25] B. Haasdonk. Reduced Basis Methods for Parametrized PDEs – A Tutorial Introduction for Stationary and Instationary Problems, chapter 2, pages 65–136. Siam, 2017.
  • [26] B. Haasdonk and M. Ohlberger. Reduced basis method for finite volume approximations of parametrized linear evolution equations. ESAIM: M2AN, 42(2):277–302, 2008.
  • [27] J. S. Hesthaven, G. Rozza, and B. Stamm. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Springer Briefs in Mathematics. Springer International Publishing, 2015.
  • [28] A. Iollo and D. Lombardi. Advection modes by optimal mass transfer. Phys. Rev. E, 89:022923, 2014.
  • [29] S. Kang, H. Choi, and S. Lee. Laminar flow past a rotating circular cylinder. Physics of Fluids, 11(11):3312–3321, 1999.
  • [30] E. N. Karatzas, F. Ballarin, and G. Rozza. Projection-based reduced order models for a cut finite element method in parametrized domains. Computers & Mathematics with Applications, in press, 2019.
  • [31] T. Lassila, A. Manzoni, A. Quarteroni, and G. Rozza. Model Order Reduction in Fluid Dynamics: Challenges and Perspectives, volume 9 of MS&A, pages 235–273. Springer International Publishing, 2014.
  • [32] T. Lassila, A. Quarteroni, and G. Rozza. A reduced basis model with parametric coupling for fluid-structure interaction problems. SIAM J. Sci. Comput., 34(2):1187–1213, 2012.
  • [33] T. Lassila, A. Quarteroni, and G. Rozza. A reduced basis model with parametric coupling for fluid-structure interaction problems. SIAM Journal on Scientific Computing, 34(2):A1187–A1213, 2012.
  • [34] A. Logg, K. Mardal, and G. Wells. Automated Solution of Differential Equations by the Finite Element Method. Springer-Verlag, Berlin, 2012.
  • [35] Y. Maday, O. Mula, A. T. Patera, and M. Yano. The generalized empirical interpolation method: stability theory on Hilbert spaces with an application to the stokes equation. Computer Methods in Applied Mechanics and Engineering, 287:310–334, 2015.
  • [36] J. Melenk. On n-widths for elliptic problems. Journal of Mathematical Analysis and Applications, 247:272–289, 2000.
  • [37] C. Murea and S. Sy. A fast method for solving fluid–structure interaction problems numerically. International Journal for Numerical Methods in Fluids, 60:1149 – 1172, 2009.
  • [38] N. J. Nair and M. Balajewicz. Transported snapshot model order reduction approach for parametric, steady-state fluid flows containing parameter-dependent shocks. International Journal for Numerical Methods in Engineering, 117(12):1234–1262, 2019.
  • [39] N.-C. Nguyen, G. Rozza, and A. T. Patera. Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers’ equation. Calcolo, 46(3):157–185, 2009.
  • [40] M. Ohlberger and S. Rave. Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing. Comptes Rendus Mathematique, 351, 2013.
  • [41] M. Ohlberger and S. Rave. Reduced basis methods: success, limitations and future challenges. Proceedings of the Conference Algoritmy, pages 1–12, 2016.
  • [42] B. Peherstorfer. Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling. arXiv:1812.02094, 2018.
  • [43] A. Quarteroni and L. Formaggia. Mathematical modelling and numerical simulation of the cardiovascular system. Handbook of Numerical Analysis, 12:3–127, 2004.
  • [44] A. Quarteroni, M. Tuveri, and A. Veneziani. Computational vascular fluid dynamics: problems, models and methods. Computing and Visualization in Science, 2:163–197, 2000.
  • [45] J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann. The shifted proper orthogonal decomposition: a mode decomposition for multiple transport phenomena. SIAM Journal on Scientific Computing, 40(3):A1322–A1344, 2018.
  • [46] T. Richter. Fluid–structure Interactions. Model, Analysis and Finite Element, volume 118 of Lecture Notes in Computational Science and Engineering. Springer International Publishing, 2017.
  • [47] C. W. Rowley, I. G. Kevrekidis, J. E. Marsden, and K. Lust. Reduction and reconstruction for self-similar dynamical systems. Nonlinearity, 16(4):1257–1275, 2003.
  • [48] G. Rozza, D. B. P. Huynh, and A. T. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch. Comput. Meth. Eng., 15:229–275, 2008.
  • [49] G. Stabile, F. Ballarin, G. Zuccarino, and G. Rozza. A reduced order variational multiscale approach for turbulent flows. Advances in Computational Mathematics, in press, 2018.
  • [50] D. Stojković, M. Breuer, and F. Durst. Effect of high rotation rates on the laminar flow around a circular cylinder. Physics of Fluids, 14(9):3160–3178, 2002.
  • [51] D. Stojković, P. Schön, M. Breuer, and F. Durst. On the new vortex shedding mode past a rotating circular cylinder. Physics of Fluids, 15(5):1257–1260, 2003.
  • [52] M. Strazzullo, F. Ballarin, R. Mosetti, and G. Rozza. Model reduction for parametrized optimal control problems in environmental marine sciences and engineering. SIAM Journal on Scientific Computing, 40(4):B1055–B1079, 2018.
  • [53] M. Strazzullo, F. Ballarin, and G. Rozza. POD–Galerkin model order reduction for parametrized time dependent linear quadratic optimal control problems in saddle point formulation. https://arxiv.org/abs/1909.09631, 2019.
  • [54] S. Sy and C. Murea. Algorithm for solving fluid–structure interaction problem on a global moving mesh. Coupled systems mechanics, 1, 03 2012.
  • [55] T. Taddei. A registration method for model order reduction: data compression and geometry reduction. arXiv:1906.11008, 2019.
  • [56] A. R. Teymourtash and S. E. Salimipour. Compressibility effects on the flow past a rotating cylinder. Physics of Fluids, 29(1):016101, 2017.
  • [57] D. Torlo, F. Ballarin, and G. Rozza. Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs. SIAM/ASA Journal on Uncertainty Quantification, 6(4):1475–1502, 2018.
  • [58] G. Welper. h and hp-adaptive interpolation by transformed snapshots for parametric and stochastic hyperbolic pdes. arXiv:1710.11481, 2017.
  • [59] G. Welper. Interpolation of functions with parameter dependent jumps by transformed snapshots. SIAM Journal on Scientific Computing, 39(4):A1225–A1250, 2017.
  • [60] G. Welper. Transformed snapshot interpolation with high resolution transforms. arXiv:1901.01322, 2019.
  • [61] Z. Zainib, F. Ballarin, S. Fremes, P.Triverio, L. Jimeńez-Juan, and G. Rozza. Reduced order methods for parametric optimal flow control in coronary bypass grafts, towards patient–specific data assimilation. https://arxiv.org/abs/1911.01409, 2019.