An optimisation-based domain-decomposition reduced order model for the incompressible Navi, Dispense di Matematica Computazionale

Scarica An optimisation-based domain-decomposition reduced order model for the incompressible Navi e più Dispense in PDF di Matematica Computazionale solo su Docsity! An optimisation–based domain–decomposition reduced order model for the incompressible Navier-Stokes equations Ivan Prusak∗a, Monica Nonino†b, Davide Torlo‡a, Francesco Ballarin§c, and Gianluigi Rozza ¶a aMathematics Area, mathLab, SISSA, 34136 Trieste, Italy bFakultät für Mathematik, Universität Wien, 1090 Wien, Austria cDepartment of Mathematics and Physics, Università Cattolica del Sacro Cuore, 25133 Brescia, Italy, Abstract The aim of this work is to present a model reduction technique in the frame- work of optimal control problems for partial differential equations. We combine two approaches used for reducing the computational cost of the mathematical numeri- cal models: domain–decomposition (DD) methods and reduced–order modelling (ROM). In particular, we consider an optimisation–based domain–decomposition algorithm for the parameter–dependent stationary incompressible Navier–Stokes equations. Firstly, the problem is described on the subdomains coupled at the inter- face and solved through an optimal control problem, which leads to the complete separation of the subdomain problems in the DD method. On top of that, a reduced model for the obtained optimal–control problem is built; the procedure is based on the Proper Orthogonal Decomposition technique and a further Galerkin projection. The presented methodology is tested on two fluid dynamics benchmarks: the sta- tionary backward–facing step and lid-driven cavity flow. The numerical tests show a significant reduction of the computational costs in terms of both the problem di- mensions and the number of optimisation iterations in the domain–decomposition algorithm. Keywords: domain decomposition, optimal control, reduced order modelling, computational fluid dynamics, Proper Orthogonal Decomposition ∗iprusak@sissa.it †monica.nonino@univie.ac.at ‡dtorlo@sissa.it §francesco.ballarin@unicatt.it ¶grozza@sissa.it 1 ar X iv :2 21 1. 14 52 8v 2 [ m at h. N A ] 3 A ug 2 02 3 1 Introduction In the last decades, there has been a growing interest in approximation techniques for partial differential equations (PDEs) that exploit high–performance computing within different fields of applications: industrial applications, naval engineering, aeronautics engineering, medical engineering, etc. Very often these problems have prohibitively high computational costs, and there is always the need of much more effective algorithms in order to alleviate the complexities of numerical models. Two of the most investigated and most important topics for rendering low computational costs are the reduced–order modelling for parameter–dependent PDEs [25] and domain–decomposition algorithms [42]. In the former case, equations of interest usually depend on a given set of param- eters; these parameters can describe either the physical properties of the sought quantities or the geometrical configuration of the physical domain over which the problem is posed. Model–order reduction is a technique based on the effective decoupling of the computationally expensive of- fline and usually computationally cheap online phase which provides a solution for any parameter value: for details we refer to [25]. Model order reduction has been successfully employed in different fields such as fluid dynamics [3, 10, 12, 15, 33, 43, 46, 47, 48, 49, 50, 52] and structural mechanics [6, 7, 23, 44, 51, 54]. Among the aforementioned applications a significant type of problems often emerges, namely saddle-point problems [8, 17], for which special care has to be taken in order to construct stable pairs of the reduced spaces; in particular, in fluid dynam- ics problems this is achieved by introducing so-called velocity supremisers, see, for instance, [5, 17, 33, 53]. Another very efficient way for reducing the computational complexity of numerical modes is Domain Decomposition (DD) method. Any domain decomposition method is based on the assumption that a given physical domain of interest is partitioned into subdomains; the original problem is then recast upon each subdomain yielding a family of subproblems of reduced size that are coupled to one another through the values and fluxes of the unknown solution at the subdomain interfaces [41, 42]. Very often the interface coupling is relaxed at the expense of providing an iterative process among subdomains, allowing a split of each of the subdomain solvers and making it computationally feasible. Domain decomposition methods can be extremely advantageous in the case of very complex geometries as well as in the case of multi-physics problems. The latter is even more attractive if we consider that there are often available state-of-the-art codes for a subcomponent model of a multi-physics problem which can be effectively exploited by decoupling algorithms; see, for instance, [16, 18, 27, 32]. In this paper, we bring our attention to domain–decomposition methods using an optimisation approach to ensure the coupling of the interface conditions between subdomains as it is presented, for example, in [22, 19]. In particular, we exploit both aforementioned techniques: optimisation– based domain decomposition algorithm in combination with projection–based reduced–order models. This paper is the first step towards the development of an efficient reduced–order model for an optimisation–based domain–decomposition algorithm for Fluid–Structure Interaction (FSI) problems [6]. It is even more attractive in the view of the articles [29, 30] where the authors are suggesting that this approach leads to a stable segregated model for FSI problems in the case of added-mass effect [11]; we also mention here some already successful ROM results in developing stable semi-implicit partitioned approaches, e.g., [4, 34, 35]. Very recently, authors of the paper [13] have introduced a novel partitioned approach for ROMs, where they couple either two different reduced–order models on each subdomain or a reduced–order model on one subdomain and a full–order (Finite Element) model on the other for the case of nonstationary diffusion–advection problems. In this context, the construction followed in this paper could be also applicable to the coupling presented in [13], as long as there is a way of casting functions defined on the subdomain interface onto the approximation spaces 2 Γ0Ω1 Ω2Γ𝐷,1 Γ𝑁,2 Γ𝑁,1 Γ𝐷,2 Figure 2: Domain Decomposition of the fluid domain functional J (𝑢1, 𝑢2) =: 1 2 ∫ Γ0 |𝑢1 − 𝑢2 |2 𝑑Γ. (3) Instead of (3) we can also consider the penalised or regularised functional J𝛾 (𝑢1, 𝑢2; 𝑔) =: 1 2 ∫ Γ0 |𝑢1 − 𝑢2 |2 𝑑Γ + 𝛾 2 ∫ Γ0 |𝑔 |2 𝑑Γ, (4) where 𝛾 is a constant that can be chosen to change the relative importance of the terms in (4). Thus we face an optimisation problem under PDE constraints: minimise the functional (3) (or (4)) over a suitable function 𝑔, subject to (2). 2.3 Variational Formulation of the PDE constraints For 𝑖 = 1, 2 define the following spaces and the norms with which each of them is endowed: • 𝑉𝑖 := { 𝑢 ∈ 𝐻1 (Ω𝑖 ;R2) } , | | · | |𝑉𝑖 = | | · | |𝐻1 (Ω𝑖 ) , • 𝑉𝑖,0 := { 𝑢 ∈ 𝐻1 (Ω𝑖 ;R2) : 𝑢 |Γ𝑖,𝐷 = 0 } , | | · | |𝑉𝑖,0 = | | · | |𝐻1 0 (Ω𝑖 ) , • 𝑄𝑖 := { 𝑝 ∈ 𝐿2 (Ω𝑖 ;R) } , | | · | |𝑄𝑖 = | | · | |𝐿2 (Ω𝑖 ) . Then, we define the following bilinear and trilinear forms: for i=1,2 • 𝑎𝑖 : 𝑉𝑖 ×𝑉𝑖,0 → R, 𝑎𝑖 (𝑢𝑖 , 𝑣𝑖) = 𝜈(∇𝑢𝑖 ,∇𝑣𝑖)Ω𝑖 , • 𝑏𝑖 : 𝑉𝑖 ×𝑄𝑖 → R, 𝑏𝑖 (𝑣𝑖 , 𝑞𝑖) = −(div𝑣𝑖 , 𝑞𝑖)Ω𝑖 , • 𝑐𝑖 : 𝑉𝑖 ×𝑉𝑖 ×𝑉𝑖,0 → R, 𝑐𝑖 (𝑢𝑖 , 𝑤𝑖 , 𝑣𝑖) = ((𝑢𝑖 · ∇)𝑤𝑖 , 𝑣𝑖)Ω𝑖 , where (·, ·)𝜔 indicates the 𝐿2 (𝜔) inner product. Consequently, the variational counterpart of (2) reads as follows: for 𝑖 = 1, 2, find 𝑢𝑖 ∈ 𝑉𝑖 and 𝑝𝑖 ∈ 𝑄𝑖 s.t. 𝑎𝑖 (𝑢𝑖 , 𝑣𝑖) + 𝑐𝑖 (𝑢𝑖 , 𝑢𝑖 , 𝑣𝑖) + 𝑏𝑖 (𝑣𝑖 , 𝑝𝑖) = ( 𝑓𝑖 , 𝑣𝑖)Ω𝑖 + ( (−1)𝑖+1𝑔, 𝑣𝑖 ) Γ0 ∀𝑣𝑖 ∈ 𝑉𝑖,0, (5a) 𝑏𝑖 (𝑢𝑖 , 𝑞𝑖) = 0 ∀𝑞𝑖 ∈ 𝑄𝑖 , (5b) 𝑢𝑖 = 𝑢𝑖,𝐷 on Γ𝑖,𝐷 . (5c) 5 Remark. In general, the fluxes through an interface Γ0 for the weak formulation of Navier–Stokes equation lives in the space 𝐻− 1 2 (Γ0) so that, in theory, the definition (4) of functional J𝛾 is not justified as it includes the 𝐿2 (Γ0)–norm of the function 𝑔. Although, as it will be evident in Section 3, the family of optimisation algorithms which are used to tackle the optimal–control problem in hand, in fact, define the respective approximation of 𝑔 that belongs to the space 𝐻 1 2 (Γ0). 2.4 Optimality system One of the ways to address the constrained optimisation problem is to reformulate the initial problem in terms of a Lagrangian functional by introducing the so–called adjoint variables. In this way, the optimal solution to the original problem is sought among the stationary points of the Lagrangian, see, for instance, [21, 26]. We define the Lagrangian functional as follows: L(𝑢1, 𝑝1, 𝑢2, 𝑝2, 𝜉1, 𝜉2, 𝜆1, 𝜆2; 𝑔) := J𝛾 (𝑢1, 𝑢2; 𝑔) − 2∑︁ 𝑖=1 [𝑎𝑖 (𝑢𝑖 , 𝜉𝑖) (6) +𝑐𝑖 (𝑢𝑖 , 𝑢𝑖 , 𝜉𝑖) +𝑏𝑖 (𝜉𝑖 , 𝑝𝑖) + 𝑏𝑖 (𝑢𝑖 , 𝜆𝑖)] + 2∑︁ 𝑖=1 ( 𝑓𝑖 , 𝜉𝑖)Ω𝑖 + 2∑︁ 𝑖=1 ((−1)𝑖+1𝑔, 𝜉𝑖)Γ0 . Notice that technically we should have also included Lagrange multipliers corresponding to the non–homogeneous Dirichlet boundary conditions (5c) in the definition of the functional L, but since the functional J𝛾 does not explicitly depend on 𝑢1,𝐷 and 𝑢2,𝐷 the corresponding Dirichlet boundary conditions for the adjoint equation that we are going to derive below will be homogeneous on these parts of the boundaries. We now apply the necessary conditions for finding stationary points of L. Setting to zero the first variations w.r.t. 𝜉𝑖 and 𝜆𝑖 , 𝑖 = 1, 2 yields the state equations (5a)-(5b). Setting to zero the first variations w.r.t. 𝑢1, 𝑝1, 𝑢2 and 𝑝2 yields the adjoint equations: 𝑎𝑖 (𝜂𝑖 , 𝜉𝑖) + 𝑐𝑖 (𝜂𝑖 , 𝑢𝑖 , 𝜉𝑖) + 𝑐𝑖 (𝑢𝑖 , 𝜂𝑖 , 𝜉𝑖) + 𝑏𝑖 (𝜂𝑖 , 𝜆𝑖) = ((−1)𝑖+1𝜂𝑖 , 𝑢1 − 𝑢2)Γ0 , ∀𝜂𝑖 ∈ 𝑉𝑖,0, (7a) 𝑏𝑖 (𝜉𝑖 , 𝜇𝑖) = 0, ∀𝜇𝑖 ∈ 𝑄𝑖 . (7b) Finally, setting to zero the first variations w.r.t. 𝑔 yields the optimality condition: 𝛾(ℎ, 𝑔)Γ0 + (ℎ, 𝜉1 − 𝜉2)Γ0 = 0, ∀ℎ ∈ 𝐿2 (Γ0). (8) 2.5 Sensitivity derivatives In order to obtain the expression for the gradient of the optimisation problem in hand, we will resort to the sensitivity approach, see, for instance, [21, 26]. The approach consists of finding equations for direction derivatives of the state variable with respect to control, called sensitivities. The first derivative 𝑑J𝛾 𝑑𝑔 of J𝛾 is defined through its action on variation ?̃? as follows:〈 𝑑J𝛾 𝑑𝑔 , ?̃? 〉 = (𝑢1 − 𝑢2, ?̃?1 − ?̃?2)Γ0 + 𝛾(𝑔, ?̃?)Γ0 , (9) where ?̃?1 ∈ 𝑉1,0, ?̃?2 ∈ 𝑉2,0 are the solutions to: 6 𝑎𝑖 (?̃?𝑖 , 𝑣𝑖) + 𝑐𝑖 (?̃?𝑖 , 𝑢𝑖 , 𝑣𝑖) + 𝑐𝑖 (𝑢𝑖 , ?̃?𝑖 , 𝑣𝑖) +𝑏𝑖 (𝑣𝑖 , 𝑝𝑖) = ((−1)𝑖+1?̃?, 𝑣𝑖)Γ0 ∀𝑣𝑖 ∈ 𝑉𝑖,0, (10a) 𝑏𝑖 (?̃?𝑖 , 𝑞𝑖) = 0 ∀𝑞𝑖 ∈ 𝑄𝑖 . (10b) We can make use of the adjoint equations (7) in order to find the representation of the gradient of the functional J𝛾 . Let 𝜉1 and 𝜉2 be the solutions to (7), ?̃?1 and ?̃?2 be the solutions to (10). By setting 𝜂𝑖 = ?̃?𝑖 in (7a), 𝜇𝑖 = 𝑝𝑖 in (7b), 𝑣𝑖 = 𝜉𝑖 in (10a) and 𝑞𝑖 = 𝜆𝑖 in (10b) we obtain: (𝑢1 − 𝑢2, ?̃?1 − ?̃?2)Γ0 = (?̃?, 𝜉1 − 𝜉2)Γ0 , so that it yields the explicit formula for the gradient of J𝛾 : 𝑑J𝛾 𝑑𝑔 (𝑢1, 𝑢2; 𝑔) = 𝛾𝑔 + (𝜉1 − 𝜉2) |Γ0 , (11) where 𝜉1 and 𝜉2 are determined from 𝑔 through (7). Notice that the gradient expression (11) is consistent with the optimality condition (8) derived in the previous section. 3 Gradient–based algorithm for PDE–constraint opti- misation problem In view of being able to provide a closed–form formula for the gradient for the objective functional J𝛾 , the natural way to proceed is to resort to a gradient–based iterative optimisation algorithm. In order to keep the exposition simple, we consider the following simple gradient method with a constant step size 𝛼 > 0: given a starting guess 𝑔 (0) , let 𝑔 (𝑛+1) = 𝑔 (𝑛) − 𝛼 𝑑J𝛾 𝑑𝑔 ( 𝑢 (𝑛) 1 , 𝑢 (𝑛) 2 ; 𝑔 (𝑛) ) . (12) Combining this with (11) we obtain: 𝑔 (𝑛+1) = 𝑔 (𝑛) − 𝛼 ( 𝛾𝑔 (𝑛) + (𝜉 (𝑛)1 − 𝜉 (𝑛) 2 ) |Γ0 ) , (13) or 𝑔 (𝑛+1) = (1 − 𝛼𝛾) 𝑔 (𝑛) − 𝛼(𝜉 (𝑛)1 − 𝜉 (𝑛) 2 ) |Γ0 , (14) where 𝜉 (𝑛) 1 and 𝜉 (𝑛) 2 are determined from (7) with 𝑔 replaced by 𝑔 (𝑛) . In summary, we have the following algorithm: Algorithm 1. 1. Choose 𝑔 (0) , 𝛼 > 0. 2. For n=0,1,2,... until convergence (a) Determine 𝑢 (𝑛) 1 ∈ 𝑉1, 𝑢 (𝑛)2 ∈ 𝑉2 by solving (5a)–(5b) with 𝑔 = 𝑔 (𝑛) . (b) Determine 𝜉 (𝑛) 1 ∈ 𝑉1,0, 𝜉 (𝑛)2 ∈ 𝑉2,0 by solving (7) with 𝑢1 = 𝑢 (𝑛) 1 , 𝑢2 = 𝑢 (𝑛) 2 . (c) Update 𝑔 (𝑛+1) by setting 𝑔 (𝑛+1) := (1 − 𝛼𝛾) 𝑔 (𝑛) − 𝛼 ( 𝜉 (𝑛) 1 − 𝜉 (𝑛) 2 ) |Γ0 . In practice, the typical methods used to solve problems like the one considered in this paper are Broyden–Fletcher–Goldfarb–Shanno (BFGS) and Newton Conjugate Gradient (CG) algorithms which tend to show much faster convergence and higher efficiency with respect to the steepest-decent algorithm. 7 where 𝑝𝑖,ℎ, 𝑖 = 1, 2 are the finite-element pressure solutions of the Navier-Stokes problem and the left-hand side is the scalar product which defines a norm on the space𝑉𝑖,0,ℎ. For more details, we refer to [5, 17]. Another way to apply the supremiser is to apply it directly to the reduced basis of the velocity spaces, but this might lead to parameter–dependent reduced spaces [5]. Other simplifications may work in a similar fashion, we might compare them in future works. 5.2 Reduced Basis Generation Once we obtain the homogenised snapshots 𝑢𝑖,0,ℎ and the pressure supremisers 𝑠𝑖,ℎ for 𝑖 = 1, 2, we are ready to construct a set of reduced basis functions. A very common choice when dealing with Navier-Stokes equations is to use the Proper Orthogonal Decomposition (POD) technique, which is based on the Singular Value Decomposition of the snapshot matrices; see, for instance, [25]. In order to implement this technique we will need two main ingredients: the matrices of the inner products and the snapshot matrices. First, we define the basis functions for the FE element spaces used in the weak formulation (17), (18) and (20) as follows: U𝑖,0,ℎ = { 𝜙 𝑢𝑖 1 , ..., 𝜙 𝑢𝑖 N𝑢𝑖 ℎ } − the FE basis of the space 𝑉𝑖,0,ℎ, 𝑖 = 1, 2, P𝑖,ℎ = { 𝜙 𝑝𝑖 1 , ..., 𝜙 𝑝𝑖 N𝑝𝑖 ℎ } − the FE basis of the space 𝑄𝑖,ℎ, 𝑖 = 1, 2, Ξ𝑖,0,ℎ := U𝑖,0,ℎ, N 𝜉𝑖 ℎ := N𝑢𝑖 ℎ , 𝑖 = 1, 2, G𝑖,ℎ = { 𝜙 𝑔 1 , ..., 𝜙 𝑔 N𝑔 ℎ } − the FE basis of the space 𝑋ℎ, where N∗ ℎ , ∗ ∈ {𝑢1, 𝑝1, 𝑢2, 𝑝2, 𝑔} denotes the dimension of the corresponding FE space. We proceed by building the snapshot matrices. In doing so we sample a parameter space and draw a discrete set of 𝑀 parameter values; there are various sampling techniques, among which we point out the uniform sampling. Then, the snapshots are taken as a high–fidelity, i.e. Finite Element, solutions at each parameter value in the sampling set. We proceed by building the snapshot matrices S𝑢𝑖 ∈ RN 𝑠 ℎ ×4𝑀 , S𝑠𝑖 ∈ RN 𝑠 ℎ ×4𝑀 , S𝑝𝑖 ∈ RN 𝑠 ℎ ×4𝑀 , S𝜉𝑖 ∈ R N𝑎 ℎ ×2𝑀 for 𝑖 = 1, 2 and S𝑔 ∈ RN 𝑔 ℎ ×𝑀 defined as follows: S𝑢1 = [𝑢1 1,0,ℎ, ..., 𝑢 𝑀 1,0,ℎ, 0, ..., 0, 0, ..., 0, 0, ..., 0], S𝑠1 = [𝑠1 1,ℎ, ..., 𝑠 𝑀 1,ℎ, 0, ..., 0, 0, ..., 0, 0, ..., 0], S𝑝1 = [0, ..., 0, 𝑝1 1,ℎ, ..., 𝑝 𝑀 1,ℎ, 0, ..., 0, 0, ..., 0], S𝑢2 = [0, ..., 0, 0, ..., 0, 𝑢1 2,0,ℎ, ..., 𝑢 𝑀 2,0,ℎ, 0, ..., 0], S𝑠2 = [0, ..., 0, 0, ..., 0, 𝑠1 2,ℎ, ..., 𝑠 𝑀 2,ℎ, 0, ..., 0], S𝑝2 = [0, ..., 0, 0, ..., 0, 0, ..., 0, 𝑝1 2,ℎ, ..., 𝑝 𝑀 2,ℎ], S𝜉1 = [𝜉1 1,ℎ, ..., 𝜉 𝑀 1,ℎ, 0, ..., 0], S𝜉2 = [0, ..., 0, 𝜉1 2,ℎ, ..., 𝜉 𝑀 2,ℎ], S𝑔 = [𝑔1 ℎ , ..., 𝑔𝑀 ℎ ], where N 𝑠 ℎ = N𝑢1 ℎ + N 𝑝1 ℎ + N𝑢2 ℎ + N 𝑝2 ℎ , N𝑎 ℎ = N 𝜉1 ℎ + N 𝜉2 ℎ and 𝑀 is the number of snapshots. Notice that since all the snapshots of the variables 𝜉1,ℎ and 𝜉2,ℎ are divergence-free on the domain of definition, the reduced spaces constructed for those variables will already contain this information, so that it allows us not to store the snapshots of the variables 𝜆1,ℎ and 𝜆2,ℎ, which 10 are playing the role of the Lagrange multipliers relative to the divergence free-conditions, as they do not contain any important information. The next step is to define the inner-product matrices 𝑋𝑢𝑖 , 𝑋𝑝𝑖 , 𝑋𝜉𝑖 for 𝑖 = 1, 2 and 𝑋𝑔. These matrices have the block diagonal structure as follows: 𝑋𝑢1 = diag ( 𝑥𝑢1 , 0𝑝1 , 0𝑢2 , 0𝑝2 ) , 𝑋𝑝1 = diag ( 0𝑢1 , 𝑥𝑝1 , 0𝑢2 , 0𝑝2 ) , 𝑋𝑢2 = diag ( 0𝑢1 , 0𝑝1 , 𝑥𝑢2 , 0𝑝2 ) , 𝑋𝑝2 = diag ( 0𝑢1 , 0𝑝1 , 0𝑢2 , 𝑥𝑝2 ) , 𝑋𝜉1 = diag ( 𝑥𝑢1 , 0𝜉2 ) , 𝑋𝜉2 = diag ( 0𝜉1 , 𝑥𝑢2 ) , 𝑋𝑔 = 𝑥𝑔 . Above, we used the following notations: 0∗ ∈ RN ∗ ℎ ×N∗ ℎ is a zero square matrix of dimension 𝑁∗ ℎ × N∗ ℎ , where ∗ ∈ {𝑢1, 𝑝1, 𝑢2, 𝑝2, 𝜉1, 𝜉2, 𝑔} and (𝑥𝑢𝑖 ) 𝑗𝑘 = ( ∇𝜙𝑢𝑖 𝑘 ,∇𝜙𝑢𝑖 𝑗 ) Ω𝑖 , for 𝑗 , 𝑘 = 1, ...,N𝑢𝑖 ℎ , 𝑖 = 1, 2, (𝑥𝑝𝑖 ) 𝑗𝑘 = ( 𝜙 𝑝𝑖 𝑘 , 𝜙 𝑝𝑖 𝑗 ) Ω𝑖 , for 𝑗 , 𝑘 = 1, ...,N 𝑝𝑖 ℎ , 𝑖 = 1, 2, (𝑥𝑔) 𝑗𝑘 = ( 𝜙 𝑔 𝑘 , 𝜙 𝑔 𝑗 ) Γ0 , for 𝑗 , 𝑘 = 1, ...,N𝑔 ℎ . We are now ready to introduce the correlation matrices C𝑢𝑖 , C𝑠𝑖 , C𝑝𝑖 , C𝜉𝑖 for 𝑖 = 1, 2 and C𝑔, all of dimension 𝑀 × 𝑀 , as: C∗ := S𝑇 ∗ 𝑋∗𝑆∗ for every ∗ ∈ {𝑢1, 𝑝1, 𝑢2, 𝑝2, 𝜉1, 𝜉2, 𝑔} and C𝑠𝑖 := S𝑇 𝑠𝑖 𝑋𝑢𝑖 𝑆𝑠𝑖 , 𝑖 = 1, 2. Once we have built the correlation matrices, we are able to carry out a POD compression on the sets of snapshots. This can be achieved by solving the following eigenvalue problems: C∗Q∗ = Q∗Λ∗ (22) where ∗ ∈ {𝑢1, 𝑠1, 𝑝1, 𝑢2, 𝑠2, 𝑝2, 𝜉1, 𝜉2, 𝑔}, Q∗ is the eigenvectors matrix and Λ∗ is the diagonal eigenvalues matrix with eigenvalues ordered by decreasing order of their magnitude. The 𝑘-th reduced basis function for the component ∗ is then obtained by applying the matrix S∗ to 𝑣∗ 𝑘 – the 𝑘-th column vector of the matrix Q∗: Φ∗ 𝑘 := 1√︃ 𝜆∗ 𝑘 S∗𝑣∗𝑘 , where 𝜆∗ 𝑘 is the 𝑘-th eigenvalue from (22). Therefore, we are able to form the set of reduced basis as A𝑠 := ⋃ ∗∈{𝑢1 ,𝑠1 , 𝑝1 ,𝑢2 ,𝑠2 , 𝑝2 } { Ψ∗ 1 , ...,Ψ ∗ 𝑁∗ } , A𝑎 := ⋃ ∗∈{ 𝜉1 , 𝜉2 } { Ψ∗ 1 , ...,Ψ ∗ 𝑁∗ } , 11 A𝑔 := { Φ 𝑔 1 , ...,Φ 𝑔 𝑁𝑔 } , where the integer numbers 𝑁∗ indicate the number of the basis functions used for each component and Ψ 𝑢1 𝑘 = ©­­­« Φ 𝑢1 𝑘 0 0 0 ª®®®¬ , Ψ 𝑠1 𝑘 = ©­­­« Φ 𝑠1 𝑘 0 0 0 ª®®®¬ , Ψ 𝑝1 𝑘 = ©­­­« 0 Φ 𝑝1 𝑘 0 0 ª®®®¬ , Ψ 𝑢2 𝑘 = ©­­­« 0 0 Φ 𝑢2 𝑘 0 ª®®®¬ , Ψ 𝑠2 𝑘 = ©­­­« 0 0 Φ 𝑠2 𝑘 0 ª®®®¬ , Ψ 𝑝2 𝑘 = ©­­­« 0 0 0 Φ 𝑝2 𝑘 ª®®®¬ , Ψ 𝜉1 𝑘 = ( Φ 𝜉1 𝑘 0 ) , Ψ 𝜉2 𝑘 = ( 0 Φ 𝜉2 𝑘 ) . We note that the first and the third blocks include both the 𝑢1, 𝑠1 and the 𝑢2, 𝑠2 basis functions - it is here that we use the pressure supremiser enrichment of the velocities spaces discussed at the beginning of this section. We provide the following renumbering of the functions for further simplicity: Φ 𝑢𝑖 𝑁𝑢𝑖 +𝑘 := Φ 𝑠𝑖 𝑘 , Ψ 𝑢𝑖 𝑁𝑢𝑖 +𝑘 := Ψ 𝑠𝑖 𝑘 , for 𝑘 = 1, ..., 𝑁𝑠𝑖 , 𝑖 = 1, 2, and we redefine 𝑁𝑢𝑖 := 𝑁𝑢𝑖 + 𝑁𝑠𝑖 , 𝑖 = 1, 2. Finally, we introduce three separate reduced basis spaces - for the state, the adjoint and the control variables, respectively: 𝑉 𝑠 𝑁 = spanA𝑠 , dim𝑉 𝑠 𝑁 = 𝑁𝑢1 + 𝑁𝑝1 + 𝑁𝑢2 + 𝑁𝑝2 , 𝑉𝑎 𝑁 = spanA𝑎 , dim𝑉 𝑠 𝑁 = 𝑁𝜉1 + 𝑁𝜉2 , 𝑉 𝑔 𝑁 = spanA𝑔, dim𝑉 𝑠 𝑁 = 𝑁𝑔 . 5.3 Online Phase Once we have introduced the reduced basis spaces we can define the reduced function expansions 𝑈𝑁 = (𝑢1,0,𝑁 , 𝑝1,𝑁 , 𝑢2,0,𝑁 , 𝑝2,𝑁 ) ∈ 𝑉 𝑠 𝑁 ,Ξ𝑁 = (𝜉1,𝑁 , 𝜉2,𝑁 ) ∈ 𝑉𝑎 𝑁 , 𝑔𝑁 ∈ 𝑉 𝑔 𝑁 as 𝑢𝑖,0,𝑁 := 𝑁𝑢𝑖∑︁ 𝑘=1 𝑢𝑖,0,𝑘Φ 𝑢𝑖 𝑘 , 𝑖 = 1, 2, 𝜉𝑖,𝑁 := 𝑁𝜉𝑖∑︁ 𝑘=1 𝜉 𝑖,𝑘 Φ 𝜉𝑖 𝑘 , 𝑖 = 1, 2, 𝑝𝑖,𝑁 := 𝑁𝑝𝑖∑︁ 𝑘=1 𝑝 𝑖,𝑘 Φ 𝑝𝑖 𝑘 , 𝑖 = 1, 2, 𝑔𝑁 := 𝑁𝑔∑︁ 𝑘=1 𝑔 𝑘 Φ 𝑔 𝑘 . In the previous equations, the underscore indicates the coefficients of the basis expansion of the reduced solution. Then the online reduced problem reads as follows: minimise over 𝑔𝑁 ∈ 𝑉 𝑔 𝑁 the functional J𝛾,𝑁 (𝑢1,𝑁 , 𝑢2,𝑁 ; 𝑔𝑁 ) := 1 2 ∫ Γ0 𝑢1,𝑁 − 𝑢2,𝑁 2 𝑑Γ + 𝛾 2 ∫ Γ0 |𝑔𝑁 |2 𝑑Γ (23) where 𝑢1,𝑁 = 𝑢1,0,𝑁 + 𝑙1,𝑁 , 𝑢2,𝑁 = 𝑢2,0,𝑁 + 𝑙2,𝑁 for (𝑢1,0,𝑁 , 𝑝1,𝑁 , 𝑢2,0,𝑁 , 𝑝2,𝑁 ) ∈ 𝑉 𝑠 𝑁 satisfy the following reduced equations ∀𝑣𝑁 = (𝑣1,𝑁 , 𝑞1,𝑁 , 𝑣2,𝑁 , 𝑞2,𝑁 ) ∈ 𝑉 𝑠 𝑁 : 12 Physical parameters 2 : 𝜈, ?̄? Range 𝜈 [0.5, 2] Range ?̄? [0.5, 6.5] Resulting 𝑅𝑒 number [0.75, 40] FE velocity order 2 FE pressure order 1 Total number of FE dofs 27,890 Number of FE dofs at the interface 130 Optimisation algorithm L-BFGS-B 𝐼𝑡𝑚𝑎𝑥 40 𝑇𝑜𝑙𝑜𝑝𝑡 10−5 𝑀 900 𝑁𝑚𝑎𝑥 50 Table 1: Computational details of the offline stage. preserve all the necessary physical information in the reduced model. In particular, we can see that a higher number of modes is needed to correctly represent the adjoint variables 𝜉1 and 𝜉2. Figures 6–9 represent the first four POD modes for each of the variables 𝑢1, 𝑢2, 𝑠1, 𝑠2, 𝑝1, 𝑝2, 𝜉1 and 𝜉2. We stress that the POD modes were obtained separately for each component and the resulting figures are obtained by gluing the subdomain function just for the sake of visualisation. (a) The first mode (b) The second mode (c) The third mode (d) The fourth mode Figure 6: The first POD modes for the velocities 𝑢1 and 𝑢2 (subdomain functions are glued together for visualisation purposes). 15 (a) POD singular values as a function of number 𝑛 of POD modes (log scaling in 𝑦-direction) (b) Energy retained by the first 𝑁𝑚𝑎𝑥 POD modes (log scaling in 𝑥-direction) Figure 5: Results of the offline stage: POD singular eigenvalue decay (a) and retained energy (b) of the first 𝑁𝑚𝑎𝑥 POD modes (a) The first mode (b) The second mode (c) The third mode (d) The fourth mode Figure 7: The first POD modes for the pressure supremisers 𝑠1 and 𝑠2 (subdomain functions are glued together for visualisation purposes). Figure 6 shows the first modes for the fluid velocities 𝑢1 and 𝑢2: in particular, notice that the modes corresponding to 𝑢1 (on the left section of the domain) are zero at the inlet boundary due to the use of lifting function. In Figure 7, we can see the first four modes for 𝑠1 and 𝑠2: here, the corresponding functions are mostly localised inside the domains Ω1 and Ω2 thanks to the homogeneous conditions at the boundaries and the non-zero forcing term coming from the pressure. 16 (a) The first mode (b) The second mode (c) The third mode (d) The fourth mode Figure 8: The first POD modes for the pressures 𝑝1 and 𝑝2 (subdomain functions are glued together for visualisation purposes). (a) The first mode (b) The second mode (c) The third mode (d) The fourth mode Figure 9: The first POD modes for the adjoint velocities 𝜉1 and 𝜉2 (subdomain functions are glued together for visualisation purposes). Figure 8 represents the first modes for the pressures 𝑝1 and 𝑝2: we point out the signs of the oscillation behaviour, which suggests that the supremiser enrichment might be needed to assure stability of the reduced–order solution. Finally, Figure 9 shows the first four modes for the adjoint variables 𝜉1 and 𝜉2: note that they are concentrated only around the interface Γ0 because the only nonzero contribution in the adjoint equations is coming from the source terms, which are defined solely on the interface Γ0. 17 Abs. error 𝑝ℎ := | |𝑝𝑖,ℎ − 𝑝ℎ | |𝐿2 (Ω𝑖 ) on domain Ω𝑖 , Rel. error 𝑝ℎ := | |𝑝𝑖,ℎ − 𝑝ℎ | |𝐿2 (Ω𝑖 ) | |𝑝ℎ | |𝐿2 (Ω𝑖 ) on domain Ω𝑖 , for 𝑖 = 1, 2. (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 (d) Iteration 40 Figure 13: High–fidelity solution for the pressures 𝑝1 and 𝑝2. Values of the parameters ?̄? = 4, 𝜈 = 0.75 and 𝑅𝑒 ≈ 19 Iteration Functional Value Gradient norm 0 7.902 2.213 5 1.956 1.210 10 0.403 2.132 40 0.007 0.069 Table 4: Functional values and the gradient norm for the FOM optimisation solution at parameter values ?̄? = 4, 𝜈 = 0.75 and 𝑅𝑒 ≈ 19 Iteration Abs. error 𝑢ℎ Rel. error 𝑢ℎ Abs. error 𝑝ℎ Rel. error 𝑝ℎ Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 0 0.2520 11.9830 0.0181 1.0000 31.6121 21.1630 0.5859 1.0000 5 0.6639 5.0075 0.0478 0.4179 20.7060 10.2359 0.3838 0.4837 10 0.2704 1.3722 0.0195 0.1145 6.7317 2.8262 0.1248 0.1335 40 0.0865 0.2566 0.0062 0.0214 1.4498 0.6443 0.0269 0.0304 Table 5: Absolute and relative errors of the FOM optimisation solution with respect to the monolithic solution at the parameter values ?̄? = 4, 𝜈 = 0.75 and 𝑅𝑒 ≈ 19 20 Iteration Functional Value Gradient norm 0 4.8 · 10−1 0.391 5 5.4 · 10−3 0.047 10 3.6 · 10−4 0.015 Table 6: Functional values and the gradient norm for the ROM optimisation solution at parameter values ?̄? = 1, 𝜈 = 1 and 𝑅𝑒 = 3 Iteration Abs. error 𝑢𝑁 Rel. error 𝑢𝑁 Abs. error 𝑝𝑁 Rel. error 𝑝𝑁 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 0 0.0284 2.9935 0.0083 1.0000 10.9522 7.0679 0.5198 1.0000 5 0.0746 0.1956 0.0217 0.0653 0.8548 0.5672 0.0406 0.0803 10 0.0135 0.0357 0.0039 0.0119 0.1714 0.1186 0.0081 0.0168 Table 7: Absolute and relative errors of the ROM optimisation solution with respect to the monolithic solution at the parameter values ?̄? = 1, 𝜈 = 1 and 𝑅𝑒 = 3 Figures 14 – 17 represent the reduced–order solutions for two different values of the param- eters (?̄?, 𝜈) = (1, 1) and 𝑅𝑒 = 3 and (?̄?, 𝜈) = (4, 0.75) and 𝑅𝑒 ≈ 19. In each of the cases, we choose the following number of the reduced basis functions: 𝑁𝑢1 = 𝑁𝑠1 = 𝑁𝑝1 = 𝑁𝑢2 = 𝑁𝑠2 = 𝑁𝑝2 = 𝑁𝑔 = 10 and 𝑁𝜉1 = 𝑁𝜉2 = 30. As was previously anticipated, we use a higher number for the adjoint variables 𝜉1 and 𝜉2 since they show much slower decay of the singular values (see Figure 5a). The solutions were obtained by carrying out 10 optimisation iterations of L–BFGS–B algorithm. Figures 14 and 16 show the intermediate solutions at iteration 0, 5 and 10 for the fluid velocities 𝑢1 and 𝑢2, whereas Figures 15 and 17 show the corresponding pressures 𝑝1 and 𝑝2. The final solution, at the 10th iteration, shows continuity between subdomain solutions at the interface Γ0. (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 Figure 14: Reduced order solution for the velocities 𝑢1 and 𝑢2. Values of the parameters ?̄? = 1, 𝜈 = 1 and 𝑅𝑒 = 3. Number of POD modes: 10 - for each state variable, each supremiser and the control, 30 – for both adjoint velocities 21 (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 Figure 15: Reduced order solution for the pressures 𝑝1 and 𝑝2. Values of the parameters ?̄? = 1, 𝜈 = 1 and 𝑅𝑒 = 3. Number of POD modes: 10 - for each state variable, each supremiser and the control, 39 – for both adjoint velocities Iteration Functional Value Gradient norm 0 7.869 2.120 5 0.107 0.401 10 0.060 0.555 Table 8: Functional values and the gradient norm for the ROM optimisation solution at parameter values ?̄? = 4, 𝜈 = 0.75 and 𝑅𝑒 ≈ 19 (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 Figure 16: Reduced order solution for the velocities 𝑢1 and 𝑢2. Values of the parameters ?̄? = 4, 𝜈 = 0.75 and 𝑅𝑒 ≈ 19. Number of POD modes: 10 - for each state variable, each supremiser and the control, 39 – for both adjoint velocities 22 Γ𝑙𝑖𝑑 Γ𝑤𝑎𝑙𝑙 Γ𝑤𝑎𝑙𝑙 Γ𝑤𝑎𝑙𝑙Ω (a) Physical domain Ω1 Ω2 Ω2 (b) Domain splitting Figure 18: Lid-driven cavity flow geometry our test case. The best performance has been achieved by using the limited-memory Broy- den–Fletcher–Goldfarb–Shanno (L-BFGS-B) optimisation algorithm, and two stopping criteria are applied: either the maximal number of iteration 𝐼𝑡𝑚𝑎𝑥 is reached or the gradient norm of the target functional is less than the given tolerance 𝑇𝑜𝑙𝑜𝑝𝑡 . Physical parameters 2 : 𝜈, ?̄? Range 𝜈 [0.05, 2] Range ?̄? [0.5, 10] Resulting 𝑅𝑒 number [0.25, 200] FE velocity order 2 FE pressure order 1 Total number of FE dofs 14,867 Number of FE dofs at the interface 138 Optimisation algorithm L-BFGS-B 𝐼𝑡𝑚𝑎𝑥 100 𝑇𝑜𝑙𝑜𝑝𝑡 10−6 𝑀 300 𝑁𝑚𝑎𝑥 100 Table 11: Computational details of the offline stage. Snapshots are sampled from a training set of 𝑀 parameters uniformly distributed in the 2-dimensional parameter space, and the first 𝑁𝑚𝑎𝑥 POD modes have been retained. Figure 19a 25 0 20 40 60 80 100 n 10 15 10 12 10 9 10 6 10 3 100 n/ m ax Singular values u1 sp1 p1 u2 sp2 p2 1 2 g (a) POD singular values as a function of number 𝑛 of POD modes (log scaling in 𝑦-direction) 100 101 n 0.6 0.7 0.8 0.9 1.0 E n Retained energy u1 sp1 p1 u2 sp2 p2 1 2 g (b) Energy retained by the first 𝑁𝑚𝑎𝑥 POD modes (log scaling in 𝑥-direction) Figure 19: Results of the offline stage: POD singular eigenvalue decay (a) and retained energy (b) of the first 𝑁𝑚𝑎𝑥 POD modes shows POD singular values for all the state, the adjoint and the control variables. As it can be seen, the POD singular values corresponding to the adjoint velocities 𝜉1 and 𝜉2 feature a slower decay compared to the one for the other variables. In Figure 19b, we can see the behaviour of the energy 𝐸𝑛 retained by the first 𝑁 modes for different components of the solution. Note that, as it was in the previous numerical case, a higher number of modes is needed to correctly represent the adjoint variables 𝜉1 and 𝜉2. Figures 20–23 represent first three POD modes for the variables 𝑢1, 𝑢2, 𝑠1, 𝑠2, 𝑝1, 𝑝2 and 𝜉1, 𝜉2. We stress that the POD modes were obtained separately for each component and the resulting figures are obtained by gluing the subdomain functions just for the sake of visualisation. Figure 20 shows the first modes for the fluid velocities 𝑢1 and 𝑢2. In particular, we notice that the modes corresponding to 𝑢2 (on the upper section of the domain) are zero at the lid boundary due to the use of lifting function. Figure 23 shows the first three modes for the adjoint variables 𝜉1 and 𝜉2: note that they are concentrated only around the interface Γ0 because the only nonzero contribution in the adjoint equations is coming from the source terms, which are defined solely on the interface Γ0. Figure 20: The first POD modes for the velocities 𝑢1 and 𝑢2 (subdomain functions are glued together for visualisation purposes). 26 Figure 21: The first POD modes for the supremiser variables 𝑠1 and 𝑠2 (subdomain functions are glued together for visualisation purposes). Figure 22: The first POD modes for the pressures 𝑝1 and 𝑝2 (subdomain functions are glued together for visualisation purposes). Figure 23: The first POD modes for the adjoint velocities 𝜉1 and 𝜉2 (subdomain functions are glued together for visualisation purposes). 27 (a) Iteration 0 (b) Iteration 3 (c) Iteration 10 Figure 26: Reduced-order solution for the velocities 𝑢1 and 𝑢2. Values of the parameters ?̄? = 5, 𝜈 = 0.05 and with 𝑅𝑒 = 100. Number of POD modes: 10 - for each state variable, each supremiser and the control, 15 – for both adjoint velocities Iteration Functional Value Gradient norm 0 4.8 · 10−1 3.153 3 2.4 · 10−2 1.634 10 7.2 · 10−5 0.023 Table 16: Functional values and the gradient norm for the ROM optimisation solution at parameter values ?̄? = 5, 𝜈 = 0.05 and with 𝑅𝑒 = 100 Iteration Abs. error 𝑢𝑁 Rel. error 𝑢𝑁 Abs. error 𝑝𝑁 Rel. error 𝑝𝑁 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 0 0.3411 0.1796 1.0000 0.1523 0.2431 0.2519 1.0000 0.1864 3 0.0512 0.0552 0.1501 0.0468 0.0531 0.0646 0.3634 0.0478 10 0.0050 0.0056 0.0147 0.0047 0.0139 0.0139 0.0956 0.0103 Table 17: Absolute and relative errors of the ROM optimisation solution with respect to the monolithic solution at the parameter values ?̄? = 5, 𝜈 = 0.05 and with 𝑅𝑒 = 100 We present additional details in Tables 16 - 19. In particular, in Tables 16 and 18 we list the values for the functional J𝛾 and the 𝐿2 (Γ0)-norm of the gradient 𝑑J𝛾 𝑑𝑔 at the different iteration of the optimisation procedure, while Table 17 and Table 19 contain the 𝐿2-relative errors with respect to the monolithic (the entire–domain) solutions 𝑢ℎ, 𝑝ℎ. Analyzing the results, we are able to see that the reduced basis method gives us a solution as accurate as the high–fidelity model. The reduced–order approximation of the optimisation problem at hand allowed us to reduce the dimension of the high-fidelity optimisation functional by more than 10-20 times and enabled us to use half optimisation algorithm iterations (each optimisation step requires at least one solve of the state and the adjoint equations). 30 (a) Iteration 0 (b) Iteration 3 (c) Iteration 10 Figure 27: Reduced-order solution for the velocities 𝑢1 and 𝑢2. Values of the parameters ?̄? = 1, 𝜈 = 0.1 and with 𝑅𝑒 = 10. Number of POD modes: 10 - for each state variable, each supremiser and the control, 15 – for both adjoint velocities Iteration Functional Value Gradient norm 0 2.6 · 10−2 2.6 · 10−1 3 1.5 · 10−5 1.0 · 10−2 10 7.1 · 10−7 1.2 · 10−3 Table 18: Functional values and the gradient norm for the ROM optimisation solution at the parameter values ?̄? = 1, 𝜈 = 0.1 and with 𝑅𝑒 = 10 Iteration Abs. error 𝑢𝑁 Rel. error 𝑢𝑁 Abs. error 𝑝𝑁 Rel. error 𝑝𝑁 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 Ω1 Ω2 0 0.0668 0.0591 1.0000 0.2424 0.0349 0.0403 1.0000 0.0936 3 0.0010 0.0019 0.0155 0.0076 0.0024 0.0020 0.0752 0.0047 10 0.0004 0.0004 0.0066 0.0017 0.0020 0.0019 0.0621 0.0045 Table 19: Absolute and relative errors of the ROM optimisation solution with respect to the monolithic solution at the parameter values ?̄? = 1, 𝜈 = 0.1 and with 𝑅𝑒 = 10 In order to provide more visually representable results (the scale of the solution on the subdomains Ω1 and Ω2 has a few orders of the difference in the magnitude), we provide the graphs of the velocities 𝑢1 and 𝑢2 separately in Figures 28 and 29. Additionally, in Table 20 we provide a comparison between full–order and reduced–order models in terms of the relative errors between ROM solutions with respect to the corresponding FOM solutions. The considerations drawn in the previous section are valid also for this test case. Remark. In both numerical cases presented above, it might seem that due to the fact that the non–homogeneous Dirichlet boundary condition is present only on the boundary of one of the subdomains only a few corrections are needed on this subdomain. On the other hand, this is true only for the velocity field, as it can be seen in the tables listing the errors (for instance in Table 3). Indeed, the errors for the pressure on those subdomains are higher than on the other one. Regarding the cavity flow, our original idea was to split the domain vertically, but in that case, the convergence even at full–order level was much slower, hence, we opted for the horizontal split. 31 (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 (d) Iteration 25 Figure 28: Reduced–order solution for the velocity 𝑢1. Values of the parameters ?̄? = 5, 𝜈 = 0.05 and with 𝑅𝑒 = 100 (a) Iteration 0 (b) Iteration 5 (c) Iteration 10 (d) Iteration 25 Figure 29: Reduced–order solution for the velocity 𝑢2. Values of the parameters ?̄? = 5, 𝜈 = 0.05 and with 𝑅𝑒 = 100 32 References [1] multiphenics - easy prototyping of multiphysics problems in FEniCS, http://mathlab.sissa.it/multiphenics, 2016. [2] RBniCS - reduced order modelling in FEniCS, http://mathlab.sissa.it/rbnics, 2015. [3] S. Ali, F. Ballarin, and G. Rozza. Stabilized reduced basis methods for parametrized steady Stokes and Navier–Stokes equations. Computers & Mathematics with Applications, 80(11):2399–2416, 2020. High-Order Finite Element and Isogeometric Methods 2019. [4] M. Astorino, F. Chouly, and M. A. Fernández. Robin Based Semi-Implicit Coupling in Fluid-Structure Interaction: Stability Analysis and Numerics. SIAM Journal on Scientific Computing, 31(6):4041–4065, 2010. [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. Reduced-order semi-implicit schemes for fluid- structure interaction problems. In Model Reduction of Parametrized Systems, pages 149– 167. Springer, 2017. [8] M. Benzi, G. H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numerica, 14:1–137, 2005. [9] G. Caloz and J. Rappaz. Numerical analysis for nonlinear and bifurcation problems. Handbook of numerical analysis, 5:487–637, 1997. [10] G. Carere, M. Strazzullo, F. Ballarin, G. Rozza, and R. Stevenson. A weighted POD- reduction approach for parametrized PDE-constrained Optimal Control Problems with ran- dom inputs and applications to environmental sciences. Computers & Mathematics with Applications, 102:261–276, 2021. [11] P. Causin, J. Gerbeau, and F. Nobile. Added-mass effect in the design of partitioned algorithms for fluid–structure problems. Computer Methods in Applied Mechanics and Engineering, 194(42):4506–4527, 2005. [12] R. Crisovan, D. Torlo, R. Abgrall, and S. Tokareva. Model order reduction for parametrized nonlinear hyperbolic problems as an application to uncertainty quantification. Journal of Computational and Applied Mathematics, 348:466–489, 2019. [13] A. de Castro, P. Kuberry, I. Tezaur, and P. Bochev. A novel partitioned approach for reduced order model – Finite Element model (ROM-FEM) and ROM-ROM coupling, 2022. [14] M. B. del Pino, T. C. Rebollo, and M. G. Mármol. A boundary-oriented reduced Schwarz do- main decomposition technique for parametric advection-diffusion problems. arXiv preprint arXiv:2305.19199, 2023. [15] S. Deparis and G. Rozza. Reduced basis method for multi-parameter-dependent steady Navier-Stokes equations: Applications to natural convection in a cavity. Journal of Com- putational Physics, 228(12):4359–4378, 2009. [16] V. Ervin, E. Jenkins, and H. Lee. Approximation of the Stokes–Darcy system by optimiza- tion. Journal of Scientific Computing, 59, 06 2014. 35 [17] A.-L. Gerner and K. Veroy. Certified reduced basis methods for parametrized saddle point problems. SIAM Journal on Scientific Computing, 34(5):A2812–A2836, 2012. [18] P. Gosselet, V. Chiaruttini, C. Rey, and F. Feyel. A monolithic strategy based on an hy- brid domain decomposition method for multiphysic problems. application to poroelasticity. Revue Européenne des Éléments Finis, 13, 04 2012. [19] M. Gunzburger and H. K. Lee. A domain decomposition method for optimization problems for partial differential equations. Computers & Mathematics with Applications, 40(2):177– 192, 2000. [20] M. Gunzburger, J. Peterson, and H. Kwon. An optimization based domain decomposition method for partial differential equations. Computers & Mathematics with Applications, 37(10):77–93, 1999. [21] M. D. Gunzburger. Perspectives in Flow Control and Optimization. Society for Industrial and Applied Mathematics, 2002. [22] M. D. Gunzburger and H. K. Lee. An optimization-based domain decomposition method for the Navier–Stokes equations. SIAM Journal on Numerical Analysis, 37(5):1455–1480, 2000. [23] B. Haasdonk. Reduced Basis Methods for Parametrized PDEs—A Tutorial Introduction for Stationary and Instationary Problems, chapter 2, pages 65–136. Society for Industrial and Applied Mathematics, 2017. [24] M. Hess, A. Alla, A. Quaini, G. Rozza, and M. Gunzburger. A localized reduced-order modeling approach for PDEs with bifurcating solutions. Computer Methods in Applied Mechanics and Engineering, 351:379–403, 2019. [25] J. S. Hesthaven, G. Rozza, and B. Stamm. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Springer Briefs in Mathematics. Springer, Switzerland, 1 edition, 2015. [26] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich. Optimization with PDE Constraints, volume 23 of Mathematical modelling. Springer, 2009. [27] T. T. P. Hoang and H. Lee. A global-in-time domain decomposition method for the coupled nonlinear Stokes and Darcy flows. Journal of Scientific Computing, 87, 04 2021. [28] A. Iollo, G. Sambataro, and T. Taddei. A one-shot overlapping Schwarz method for component-based model reduction: application to nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, 404:115786, 2023. [29] P. Kuberry and H. K. Lee. A decoupling algorithm for fluid-structure interaction prob- lems based on optimization. Computer Methods in Applied Mechanics and Engineering, 267:594–605, 2013. [30] P. Kuberry and H. K. Lee. Analysis of a fluid-structure interaction problem recast in an optimal control setting. SIAM Journal on Numerical Analysis, 53(3):1464–1487, 2015. [31] Y. A. Kuznetsov, I. A. Kuznetsov, and Y. Kuznetsov. Elements of applied bifurcation theory, volume 112. Springer, 1998. [32] J. E. Lagnese, G. Leugering, and G. Leugering. Domain in Decomposition Methods in Optimal Control of Partial Differential Equations. Number 148 in International Series of Numerical Mathematics. Springer Science & Business Media, 2004. [33] T. Lassila, A. Manzoni, A. Quarteroni, and G. Rozza. Model order reduction in fluid dynamics: challenges and perspectives. In A. Quarteroni and G. Rozza, editors, Reduced Order Methods for Modeling and Computational Reduction, volume 9, pages 235–274. Springer MS&A Series, 2014. 36 [34] M. Nonino, F. Ballarin, and G. Rozza. A monolithic and a partitioned reduced basis method for fluid–structure interaction problems. Fluids, 6(6), 2021. [35] M. Nonino, F. Ballarin, G. Rozza, and Y. Maday. Projection based semi–implicit partitioned reduced basis method for non parametrized and parametrized Fluid–Structure Interaction problems. arXiv preprint arXiv:2201.03236, 2022. [36] T. Panitz and D. Wasan. Flow attachment to solid surfaces: the Coanda effect. AIChE Journal, 18(1):51–57, 1972. [37] F. Pichi and G. Rozza. Reduced basis approaches for parametrized bifurcation problems held by non-linear Von Kármán equations. Journal of Scientific Computing, 81:112–135, 2019. [38] F. Pichi, M. Strazzullo, F. Ballarin, and G. Rozza. Driving bifurcating parametrized nonlin- ear pdes by optimal control strategies: application to Navier–Stokes equations with model order reduction. ESAIM: Mathematical Modelling and Numerical Analysis, 56(4):1361– 1400, 2022. [39] M. Pintore, F. Pichi, M. Hess, G. Rozza, and C. Canuto. Efficient computation of bifurcation diagrams with a deflated approach to reduced basis spectral element method. Advances in Computational Mathematics, 47:1–39, 2021. [40] A. Quaini, R. Glowinski, and S. Čanić. Symmetry breaking and preliminary results about a Hopf bifurcation for incompressible viscous flow in an expansion channel. International Journal of Computational Fluid Dynamics, 30(1):7–19, 2016. [41] A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equations. Springer Series in Computational Mathematics. 23. Springer Berlin Heidelberg, Heidelberg, DE, 1994. Written for: Numerical analysts, applied mathematicians. [42] A. Quarteroni and A. Valli. Domain decomposition Methods for Partial Differential Equa- tions. Oxford University Press, Oxford, UK, 1999. [43] G. Rozza. Reduced basis methods for stokes equations in domains with non-affine parameter dependence. Computing and Visualization in Science, 12:23–35, 2009. [44] G. Rozza, D. Huynh, and A. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Archives of Computational Methods in Engineering, 15:1–47, 09 2007. [45] R. Seydel. Practical bifurcation and stability analysis, volume 5. Springer Science & Business Media, 2009. [46] G. Stabile, F. Ballarin, G. Zuccarino, and G. Rozza. A reduced order variational multiscale approach for turbulent flows. Advances in Computational Mathematics, 45(5):2349–2368, 2019. [47] G. Stabile and G. Rozza. Finite volume POD-Galerkin stabilised reduced order methods for the parametrised incompressible Navier–Stokes equations. Computers & Fluids, 173:273– 284, 2018. [48] 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. Journal of Scientific Computing, 83(3):1–35, 2020. [49] M. Strazzullo, M. Girfoglio, F. Ballarin, T. Iliescu, and G. Rozza. Consistency of the full and reduced order models for evolve-filter-relax regularization of convection-dominated, marginally-resolved flows. International Journal for Numerical Methods in Engineering, 123(14):3148–3178, 2022. 37