Residual-based data-driven variational multiscale reduced order models for parameter-dependent problems

Computational and Applied Mathematics (2025) 44:308 https://doi.org/10.1007/s40314-025-03273-0 Residual-based data-driven variational multiscale reduced order models for parameter-dependent problems Birgul Koc1·Samuele Rubino2·Tomás Chacón Rebollo2·Traian Iliescu3 Received: 5 January 2025 / Revised: 22 April 2025 / Accepted: 16 May 2025 © The Author(s) 2025 Abstract In this paper, we propose a novel residual-based data-driven closure strategy for reduced order models (ROMs) of under-resolved, convection-dominated problems. The new ROM closure model is constructed in a variational multiscale (VMS) framework by using the available full order model data and a model form ansatz that depends on the ROM residual. We emphasize that this closure modeling strategy is fundamentally different from the current data-driven ROM closures, which generally depend on the ROM coefficients. We investigate the new residual-based data-driven VMS ROM closure strategy in the numerical simulation of three test problems: (i) a one-dimensional parameter-dependent advection-diffusion problem; (ii) a two-dimensional time-dependent advection-diffusion-reaction problem with a small diffusion coefficient (ε=1e−4); and (iii) a two-dimensional flow past a cylinder at Reynolds number Re =1000. Our numerical investigation shows that the new residual-based data-driven VMS-ROM is more accurate than the standard coefficient-based data-driven VMS-ROM. Keywords Reduced order models ·Variational multiscale ·Data-driven modeling ·Residual Mathematics Subject Classification 65M60 ·76-04 1 Introduction Reduced order models (ROMs) have been instrumental in significantly reducing the computational cost of full order models (FOMs) (e.g., finite element or finite volume methods) BBirgul Koc [email protected] Samuele Rubino [email protected] Tomás Chacón Rebollo [email protected] Traian Iliescu [email protected] 1Departamento EDAN, Universidad de Sevilla, Sevilla, Spain 2Departamento EDAN & IMUS, Universidad de Sevilla, Sevilla, Spain 3Department of Mathematics, Virginia Tech, Blacksburg, USA 0123456789().: V,-vol 123 308 Page 2 of 28 B. Koc et al. in applications that require repeated model runs, such as, design and control, uncertainty quantification, inverse problems, and data assimilation. However, in convection-dominated problems (e.g., turbulent flows), standard ROMs generally yield inaccurate results. Indeed, to ensure a low computational cost, relatively low-dimensional ROMs are generally used in practice. Convection-dominated problems, however, usually require a large number of ROM basis functions in order to accurately represent the underlying dynamics. Thus, standard, lowdimensional ROMs usually yield spurious numerical oscillations that significantly degrade the solution accuracy. To alleviate this inaccurate behavior, the standard ROMs are generally equipped with (i) numerical stabilizations of different types (e.g., projection-based (Azaïez et al. 2021; Chacón Rebollo et al. 2022; Novo and Rubino 2021), subspace rotation (Balajewicz et al. 2016), variational multiscale (Bergmann et al. 2009; Iliescu and Wang 2013,2014; Reyes and Codina 2020), streamline-upwind Petrov-Galerkin (Parish et al. 2020), filter-based (Girfoglio et al. 2021), and physically constraints (Sanderse 2020)); or (ii) ROM closures, which are terms that are added to the standard ROM to model the effect of the unresolved ROM scales. There are several types of ROM closure strategies, which are reviewed in Ahmed et al. (2021). The variational multiscale (VMS) ROM closures are a popular class of ROM closure strategies. The VMS-ROM closures leverage the VMS framework (Hughes et al. 1998), which has been extensively used at the FOM level (see Ahmed et al. (2017) for a review). Specifically, the physical space is first decomposed in the ROM space (i.e., the space used to construct the ROM) and the space of unresolved scales. Then, the VMS-ROM equations are obtained by projecting the underlying equations onto the ROM space. We emphasize that the VMS-ROM includes both standard ROM operators (that depend on the ROM space), as well as a ROM closure term, i.e., a ROM operator that depends on the space of unresolved scales. To obtain a practical, self-contained ROM, this closure term in the VMS-ROM needs to be modeled by using exclusively the ROM space. This, in a nutshell, is the celebrated ROM closure problem. Several successful VMS-ROM closure strategies have been proposed over the years. A review of these approaches is performed in Section IV.A.5 of Ahmed et al. (2021). Next, we outline several VMS-ROM closure models. In Bergmann et al. (2009), a residual-based VMS model was proposed as a ROM stabilization strategy. In Reyes and Codina (2020), a twoscale VMS-ROM equipped with time-dependent orthogonal sub-grid scales was developed. In Iliescu and Wang (2013), the authors proposed a VMS-ROM closure that includes an artificial viscosity added only to the small resolved scales of the gradient. The numerical tests in Iliescu and Wang (2013) showed the increased numerical stability and accuracy of the VMS-ROM over the standard G-ROM and illustrated the theoretical convergence rates. In particular, a problem displaying shock-like phenomena was considered (a 2D traveling wave) at a moderate Péclet number (ν=10−4). In Iliescu and Wang (2014), the VMS-ROM was extended and studied for the incompressible Navier-Stokes equations. Recent VMS-ROM developments can be found in, e.g., Eroglu et al. (2017); Parish and Duraisamy (2017); Roop (2013); Stabile et al. (2019); Tello et al. (2019). A paradigm shift in the development of VMS-ROM closures occurred in Mou et al. (2021) (see also Xie et al. (2018) for relevant work), where the classical, physical modeling used to develop VMS-ROM closures was replaced with data-driven modeling. Instead of using traditional arguments (e.g., eddy viscosity), the novel data-driven VMS-ROM (d2-VMS-ROM) proposed in Mou et al. (2021) was constructed by leveraging available data. Specifically, the d2-VMS-ROM was built by first postulating a model form for the closure term (i.e., a linear or quadratic model), and then solving a least squares problem to determine the model parameters that yielded the closest fit between the model form and the available data. The 123 Residual-based data-driven variational multiscale… Page 3 of 28 308 original least squares formulation in the d2-VMS-ROM was replaced with a machine learning strategy in Ahmed et al. (2023) (see Xie et al. (2020) for relevant work). The d2-VMS-ROM has been extended in different directions (e.g., adding physical constraints (Mohebujjaman et al. 2019), providing mathematical support (Koc et al. 2022), and developing a stochastic framework for efficient data assimilation (Mou et al. 2023)), and has been successfully used in challenging numerical simulations (e.g., from the quasigeostrophic equations (Mou et al. 2020) to the turbulent channel flow (Mou 2021)). Despite its significant achievements, the d2-VMS-ROM has been exclusively used in its coefficient-based form. Specifically, the data-driven closure term in the d2-VMS-ROM depends exclusively on the ROM coefficients. In this paper, we propose a fundamentally different d2-VMS-ROM strategy, in which the ROM closure term is a function of the ROM residual. The main advantage of the new residual-based ROM closure term is that it is consistent: As the ROM dimension increases, the residual decreases, and thus the ROM closure term decreases as well. This behavior is consistent with the physical role of the ROM closure term: As the ROM dimension increases, more physical scales are resolved by the d2-VMS-ROM, and thus the role of the ROM closure term decreases. We emphasize that, in contrast to the new residual-based d2-VMS-ROM, the current coefficient-based d2-VMSROM are not consistent (i.e., as the ROM dimension increases, the coefficient-based ROM closures do not necessarily decrease). In this paper, we perform a numerical investigation of the new residual-based d2-VMS-ROM and show that it is more accurate and efficient than the classical coefficient-based d2-VMS-ROM. The outline of the paper is as follows: In Sect. 2, we briefly outline the standard Galerkin ROM for the incompresible Navier-Stokes equations. In Sect. 3, we describe the smallscale to large-scale decomposition that underpins the VMS-ROM framework. In Sect. 4,we introduce two strategies for the construction of the novel residual-based d2-VMS-ROM, and outline the standard coefficient-based d2-VMS-ROM. In Sect. 5, we perform a numerical investigation of the new residual-based d2-VMSROMs in the simulation of three test problems: (i) a one-dimensional parameter-dependent advection-diffusion equation; (ii) a two-dimensional time-dependent advection-diffusionreaction with small viscosity; and (iii) a two-dimensional flow past a circular cylinder at Reynolds number Re =1000. To assess the performance of the new residual-based d2-VMSROMs, we compare them with the standard Galerkin ROM, the classical coefficient-based d2-VMS-ROM, and an ideal VMS-ROM, in which the closure term is computed from the FOM data. As benchmark for our numerical investigation, we use the FOM data. In Sect. 6, we conclude the paper with a short summary and future research directions. 2 Galerkin reduced order model (G-ROM) This section provides a brief overview of the standard Galerkin ROM (G-ROM) strategy, which is one of the most common types of ROMs for fluid flows (Hesthaven et al. 2016; Noack et al. 2011; Quarteroni et al. 2015). As a mathematical model, we use the incompressible Navier-Stokes equations (NSE) (1)–(2): ∂u ∂t−Re−1u+u·∇u+∇p=f,(1) ∇·u=0,(2) where uis the velocity, pthe pressure, fthe force, and Re the Reynolds number. For clarity of presentation, we use homogeneous Dirichlet boundary conditions and zero force, i.e., f=0. 123 308 Page 4 of 28 B. Koc et al. In Algorithm 1, we outline the construction of G-ROM, which is carried out using the velocity field. In the ROM framework, we employ a divergence-free velocity basis (utilizing Scott-Vogelius elements in the finite element setting). Consequently, the pressure term is omitted in the G-ROM. For ROMs that include the pressure approximation, see, e.g., Ballarin et al. (2015); Bergmann et al. (2009); Chacón Rebollo et al. (2022); Hesthaven et al. (2016); Noack et al. (2005); Novo and Rubino (2021); Quarteroni et al. (2015); Reyes and Codina (2020); Stabile and Rozza (2018). Algorithm 1 Galerkin ROM (G-ROM) 1: Use available FOM data to construct dominant modes by using the proper orthogonal decomposition (POD), {ϕ1,...,ϕL},Ld(where dis the dimension of the input dataset), which correspond to the largest relative kinetic energy content and represent the dominant spatial structures of the given test problem; 2: Construct a ROM velocity approximation: uL= L  j=1 (aL)jϕj,(3) as a linear combination of ROM basis functions ϕjwith ROM coefficients (aL)j; 3: Replace uin the given test problem with the ROM solution uLgivenin(3); 4: Use the Galerkin projection, which projects the system obtained in step 3 onto the ROM velocity space XLspanned by {ϕ1,...,ϕL}. By using Algorithm 1for the NSE equations (1)–(2), we obtain the following G-ROM: daL dt =ALL aL+aLBLLL aL,(4) where ALL is an L×Lmatrix with entries Aij := −Re−1(∇ϕi,∇ϕj)and BLLL is an L×L×Ltensor with entries Bijk := −(ϕi,ϕj·∇ϕk).TheG-ROM(4)isanL-dimensional system of ODEs that can be used for time intervals and/or parameters different from those used in the training regime (i.e., in the construction of the G-ROM). 3 Variational multiscale reduced order model (VMS-ROM) In this section, we construct the VMS-ROM framework, which will be used in the next sections to build the d2-VMS-ROMs. First, we note that when all the available ROM modes are used to create a ROM solution, the ROM approximation becomes ud= d  j=1 (ad)jϕj.(5) In this case, udis the most accurate ROM approximation of the FOM solution with the given data in the POD sense (i.e., from the energetic point of view (Volkwein 2013)). For laminar flows, using a few (Ld) ROM basis functions is enough to capture the main dynamics of the given problem, i.e., we are in the resolved regime.Inthisregime,a 123 Residual-based data-driven variational multiscale… Page 5 of 28 308 low-dimensional ROM solution uL, with small Ld, yields an accurate approximation of the FOM solution. However, for turbulent flows, the low-dimensional ROM solution (3) with small Ld is not an accurate approximation of the FOM solution, i.e., we are in the under-resolved regime. To increase the accuracy of the L-dimensional ROM solution (3), we generally have two options: (i) increase the G-ROM dimension, L, or (ii) add numerical stabilization or a low-dimensional closure term to the G-ROM. In this paper, we aim to increase the numerical accuracy without significantly increasing the computational cost. Thus, we choose the second option. Next, we explain what ROM closure modeling is (see Ahmed et al. (2021) for a review) and how it is performed in a VMS setting. The orthogonality of the ROM basis functions (which is intrinsic to the POD framework) allows us to decompose the ROM space as follows: Xd=XL⊕XS,where Xd:= span{ϕ1, ..., ϕd},XL:= span{ϕ1, ..., ϕL},andXS:= span{ϕL+1, ..., ϕd}.Byusing this decomposition, we define the large-scale and sub-scale solutions of the most accurate (in the POD sense) ROM solution, ud, as follows: uL:= L  j=1 (aL)jϕj,(6a) uS:= d  j=L+1 (aS)jϕj.(6b) Next, we note that the most accurate ROM approximation in (5), ud, solves the following d-dimensional weak form of NSE equations (1)–(2) Dt(vd,ud)+a(vd,ud)+b(vd,ud,ud)=0∀vd∈Xd,(7) where the bilinear forms are Dt(vd,ud):= (vd,∂ tud)and a(vd,ud):= Re−1(∇vd,∇ud), and the nonlinear form is b(vd,ud,ud):= (vd,ud·∇ud). Furthermore, (·,·)represents the inner product in L2(). By using the VMS method and choosing vL=ϕL∈XLand vS=ϕS∈XS(vd=vL+vS), we can decompose (7) into two problems as follows: Dt(ϕL,uL+uS)+a(ϕL,uL+uS)+b(ϕL,uL+uS,uL+uS)=0(8a) Dt(ϕS,uL+uS)+a(ϕS,uL+uS)+b(ϕS,uL+uS,uL+uS)=0.(8b) The matrix-vector forms of (8a)–(8b) are as follows: daL dt =ALL aL+ALS aS+a LBLLLaL+a LBLLSaS+a SBLSLaL+a SBLSSaS, (9a) daS dt =ASS aS+ASL aL+a LBSLLaL+a LBSLSaS+a SBSSLaL+a SBSSSaS. (9b) The VMS-ROM idea can be explained by using the matrix formulation in (9a)–(9b). First, we note that this matrix form was obtained by using the variational formulation in (8a)–(8b). Second, we note that the matrix formulation in (9a)–(9b) is a multiscale formulation since aL and aScorrespond to the large and small scales in the system, respectively. Thus, the matrix formulation in (9a)–(9b) is truly a variational multiscale ROM formulation. The VMS-ROM rationale is that, since Ld,aLcan be computed efficiently. In contrast, since LS, we should try to avoid the expensive computation of aS. However, the challenge is that the equations for aLand aSin (9a)–(9b) are coupled. 123 308 Page 6 of 28 B. Koc et al. Fig. 1 Diagram of the algebraic form of the coupled system (9a)–(9b) The VMS-ROM strategy centers on two simple ideas, which are illustrated in the schematic in Fig. 1: The first idea is that VMS-ROMs reduce the large system of equations (9a)–(9b)to a low-dimensional system of equations for aL. The second idea is that, to obtain an accurate approximation for aL, the effect of aSneeds to be modeled (i.e., the closure problem needs to be addressed). In the next section, we present two new data-driven strategies for modeling the effect of aS(Sect. 4.1–4.2). 4 Data-driven variational multiscale ROM (d2-VMS-ROM) In this section, we explain how we build the low-dimensional closure term in d2-VMSROM (10), which aims to increase the numerical accuracy by modeling the effect of the sub-scales: daL dt =ALL aL+aLBLLL aL+(Closure-Term).(10) In the following sections, we present two fundamentally different types of Closure-Term. In Sect. 4.1–4.2, we propose novel residual-based d2-VMS-ROMs in which the sub-scale information in the closure term is modeled by using the large-scale ROM residual, ResS(aL). Specfically, in Sect. 4.1 we propose R1-ROM which is associated with a single Closure-Term, and in Sect. 4.2 we propose R2-ROM which leverages two Closure-Terms to better control sub-scale effects. In Sect. 4.3, we describe the classical coefficient-based d2-VMS-ROM, which has recently been proposed (Mou et al. 2021). In this model, the sub-scale information in the Closure-Term is modeled by using the large-scale ROM coefficient vector, aL.Our goal in the numerical investigation in Sect. 5is to show that the new residual-based d2VMS-ROMs (R1-ROM and R2-ROM) are more accurate and efficient than the standard coefficient-based d2-VMS-ROM. To build the d2-VMS-ROM (10), we first note that, to ensure that it is an efficient, Ldimensional model, its closure term should be modeled using only the large-scale ROM coefficient, aL. The next step in the d2-VMS-ROM construction is to postulate a model form for the Closure-Term, denoted Ansatz(aL). To find the Ansatz-Operators that yield the most accurate results, we use data-driven modeling. Specifically, in the offline phase, we solve a least squares problem that minimizes the difference between Ansatz(aL)and the true sub-scale term, denoted Sub-Scale-Term(aL,aS), evaluated with the FOM sub-scale and 123 Residual-based data-driven variational multiscale… Page 7 of 28 308 large-scale coefficient vectors: min Ansatz-Operators M  k=1  Sub-Scale-Term(ak L,ak S)−Ansatz(ak L)   2 L2,(11) where Mrepresents the number of snapshots. Then, in the online phase, the d2-VMS-ROM with the Ansatz computed from (11) is used for time intervals and/or parameters different from those used in the training stage. We note that alternative data-driven strategies for ROM closure have been used in, e.g., Maulik et al. (2020); Prakash and Zhang (2024), and reviewed in Sanderse et al. (2024). 4.1 Residual-based d2-VMS-ROM with one ansatz (R1-ROM) In this section, we introduce the first residual-based d2-VMS-ROM, R1-ROM. In Algorithm 2, we outline the construction of the R1-Closure-Term. Algorithm 2 Residual-Based Closure Term with One Ansatz (R1-Closure-Term) 1: Introduce the Sub-Scale-Term, which needs to be modeled, and the form of Ansatz, which appears in the Closure-Term of R1-ROM (14): Sub-Scale-Term =aS(12) ≈Ansatz =  AResS(aL)+ResS(aL) BResS(aL), where the large-scale ROM residual, derived from (9b), has the form ResS(aL):= ASL aL+a LBSLL aL. 2: Solve the least squares problem (11) to obtain the S-dimensional (the dimension of the sub-scale) Ansatz-Operators  Aand  B, using the Sub-Scale-Term and Ansatz defined in (12). 3: Substitute all aSterms in the large-scale equation (9a) with Ansatz from (12) to derive the R1-ClosureTerm (13): R1-Closure-Term ≈ALS Ansatz +a LBLLS Ansatz (13) +AnsatzBLSLaL+AnsatzBLSS Ansatz. Replacing the Closure-Term in the d2-VMS-ROM (10) with the R1-Closure-Term (13) from Algorithm 2, we derive the following R1-ROM formulation: daL dt =ALL aL+a LBLLL aL+(R1-Closure-Term).(14) We highlight that the large-scale ROM residual, ResS(aL),in(12) does not include a time-derivative term due to the orthogonality of the POD modes in the L2inner product. The idea behind the R1-ROM strategy is to avoid solving the expensive, high-dimensional equation (9b). Instead, we only leverage the information in (9b) (i.e., the fact that aSdepends on the residual ResS(aL)) to model the sub-scales in (9a). We emphasize that the closure term in the residual-based R1-ROM (14) is consistent since it depends on the large-scale ROM residual, ResS(aL). In contrast, the closure term of the standard coefficient-based C-ROM (which is presented in Sect. 4.3) is not consistent since it depends on the large-scale trajectory, aL. 123 308 Page 8 of 28 B. Koc et al. 4.2 Residual-based d2-VMS-ROM with two ansatzes (R2-ROM) In this section, we introduce the second residual-based d2-VMS-ROM, R2-ROM. The main differences between the R1-ROM and R2-ROM are the following: (i) the R2-ROM uses two ansatzes whereas the R1-ROM uses only one ansatz, which is a common ansatz, and (ii) because the R2-ROM uses more ansatzes, it has more information related to the sub-scales. In Algorithm 3, we outline the construction of R2-ROM. Algorithm 3 Residual-Based Closure Term with Two Ansatzes (R2-Closure-Term) 1: Introduce a first sub-scale term Sub-Scale-Term1, which needs to be modeled, and the form of Ansatz1. The Sub-Scale-Term1 and its corresponding Ansatz1 are modeled similarly to R1-ROM (see (12)) to approximate the sub-scale ROM coefficient, i.e., aSin (9a): Sub-Scale-Term1 =aS(15) ≈Ansatz1 =  A1ResS(aL)+ResS(aL) B1ResS(aL). 2: Solve the least squares problem (11) to compute the first set of Ansatz operators,  A1and  B1, which are of dimension S(the sub-scale dimension), utilizing the Sub-Scale-Term1 and Ansatz1 as defined in (15); 3: Replace all instances of aSin the large-scale equation (9a) with Ansatz1 from (15); 4: Introduce a second sub-scale term (Sub-Scale-Term2) and the corresponding Ansatz2 to enhance numerical accuracy: Sub-Scale-Term2 := ResL(Ansatz1)=ALS (Ansatz1)+(Ansatz1)BLSS (Ansatz1)(16) ≈Ansatz2 =  A2ResL(aS)+(ResL(aS)) B2ResL(aS); 5: Solve the least squares problem (11) to obtain the second set of Ansatz operators,  A2and  B2, of dimension L(the large-scale dimension), based on Sub-Scale-Term2 and Ansatz2, as defined in (16); 6: Replace the term ResL(aS)in the large-scale equation (9a) with Ansatz2 from (16), and substitute all remaining instances of aSin the large-scale equation with Ansatz1 from (15), thereby deriving the R2closure term: R2-Closure-Term ≈  A2ResL(Ansatz1)+(ResL(Ansatz1)) B2ResL(Ansatz1)(17) +a LBLLS Ansatz1 +Ansatz1BLSL aL. Replacing the Closure-Term in the d2-VMS-ROM (10) with the R2-Closure-Term (17) from Algorithm 3, we derive the following R2-ROM formulation: daL dt =ALL aL+a LBLLL aL+(R2-Closure-Term).(18) We emphasize that both the large-scale ROM residual, ResS(aL)in (15), and the sub-scale ROM residual, ResL(aS),in(16) do not contain a time-derivative term. This is due to the orthogonality of the POD modes in the L2inner product. Furthermore, the inclusion of a second ansatz allows us to refine the approximation of ResL(Ansatz1)using sub-scale FOM data. 123 Residual-based data-driven variational multiscale… Page 9 of 28 308 The expectation is that the residual-based R2-ROM (18) could yield more accurate results than the residual-based R1-ROM (14) since the R2-Closure-Term (17) gradually models the sub-scale and has more sub-scale information than the R1-Closure-Term (13). However, in numerical simulations, since the R2-Closure-Term (17) is more complex, it could yield inaccurate results. 4.3 Coefficient-based d2-VMS-ROM (C-ROM) In this section, we outline the coefficient-based C-ROM strategy (Mou et al. 2021). In Algorithm 4, we outline the construction of C-Closure-Term. Algorithm 4 Coefficient-Based Closure Term with One Ansatz (C-Closure-Term) 1: Introduce the sub-scale term Sub-Scale-Term, which needs to be modeled, and the form of Ansatz, which appears in the Closure-Term of C-ROM (21): Sub-Scale-Term =ALS aS+a LBLLS aS+a SBLSL aL+a SBLSS aS(19) ≈Ansatz =∗  AaL+a L  Ba L, where Sub-Scale-Term is derived from (9a). 2: Solve the least squares problem (11) to compute the Ansatz-Operators  Aand  B, which are of dimension L(the large-scale dimension), utilizing the Sub-Scale-Term and Ansatz defined in (19). 3: Substitute Sub-Scale-Term in the large-scale equation (9a) with Ansatz from (19) to derive the C-ClosureTerm (20): C-Closure-Term ≈  AaL+a L  Ba L.(20) Replacing the Closure-Term in the d2-VMS-ROM (10) with the C-Closure-Term (20) from Algorithm 4, we derive the following C-ROM formulation: daL dt =ALL aL+a LBLLL aL+(C-Closure-Term).(21) 4.4 Ideal variational multiscale ROM (I-ROM) In this section, we outline the ideal ROM (I-ROM), which is used as a benchmark model to discuss the effect of Closure-Term in (10). We emphasize that I-ROM is a purely theoretical model, used to assess the accuracy of ROM closure models. Specifically, the closure term in I-ROM is computed directly from both the large-scale and sub-scale FOM data: I-Closure-Term =ALS aS+a LBLLS aS+a SBLSL aL+a SBLSS aS.(22) Thus, the I-ROM closure term cannot be used in practical settings, where FOM data is not available. Replacing the Closure-Term in the d2-VMS-ROM (10) with the I-ClosureTerm (22), we obtain the I-ROM: aL dt =ALL aL+a LBLLL aL+(I-Closure-Term).(23) 123 308 Page 16 of 28 B. Koc et al. Fig. 7 ADR equation; L2-POD decay of eigenvalues Table 4 ADR equation; average relative L2-projection error (30) for G-ROM, I-ROM, C-ROM, R1-ROM, and R2-ROM for various Lvalues (number of POD modes) LG-ROM I-ROM C-ROM R1-ROM R2-ROM 2 2.14e+00 3.87e-01 1.75e+00 1.60e+00 1.81e+00 4 2.19e+00 4.01e-01 2.00e+00 1.73e+00 1.10e+00 6 2.13e+00 4.59e-01 2.02e+00 1.68e+00 5.25e-01 8 2.06e+00 6.02e-01 1.99e+00 1.66e+00 5.38e-01 10 1.99e+00 9.25e-01 1.95e+00 1.76e+00 5.17e-01 In Fig. 8, we plot the relative ROM error in time for L=4,6, and 10. We observe that the relative R2-ROM error recovers the I-ROM as Lincreases. Furthermore, for L=10, the relative R2-ROM error is lower than the relative I-ROM error. Finally, in Fig. 9, we plot the FOM and all ROM solutions at the final time, T=1, for L=8. We observe that the C-ROM does not diminish the oscillatory bump, whereas the R1ROM tries to decrease the magnitude of the bump. The R2-ROM decreases the oscillations the most, but it also produces another oscillatory bump. 5.3 Two-dimensional flow past a cylinder In this section, we consider a 2D channel flow past a circular cylinder at Reynolds number Re =1000. As criteria in our numerical investigation, we use the average L2ROM errors, kinetic energy, vortex shedding frequency, and Pareto plots. Computational Setting As a mathematical model, we use the NSE (1)–(2). The computational domain is a 2.2×0.41 rectangular channel with a radius =0.05 cylinder, centered at (0.2,0.2);seeFig.10. We prescribe no-slip boundary conditions on the walls and cylinder, and the following inflow and outflow profiles (John 2004; Mohebujjaman et al. 2019,2017; Rebholz and Xiao 2017): u1(0,y,t)=u1(2.2,y,t)=6 0.412y(0.41 −y), (31) u2(0,y,t)=u2(2.2,y,t)=0,(32) where u=u1,u2. There is no forcing and the flow starts from rest. Snapshot Generation For the spatial discretization, we use the pointwise divergence-free, LBB stable (P2,Pdisc 1) Scott-Vogelius FE pair on a barycenter refined regular triangular mesh of the barycenter (John 123 Residual-based data-driven variational multiscale… Page 17 of 28 308 Fig. 8 ADR equation; relative ROM error over the extrapolation testing time interval for various Lvalues et al. 2016). The mesh provides 103K(102962) velocity and 76K(76725) pressure degrees of freedom. We utilize the commonly used linearized BDF2 temporal discretization and a time step size t=0.002 for both the FOM and the ROM time discretizations. In the first time step, we use a backward Euler scheme so that we have two initial time-step solutions, as required for the BDF2 scheme. ROM Construction and Testing The FOM simulation achieves the statistically steady state after t=13 in the numerical investigation. To build the ROM basis functions and operators, we decided to use 3 time units of FOM data. Thus, we collect FOM snapshots from t=13 to t=16 and label them as training FOM data. In Fig. 11, we plot the decay of eigenvalues of the ROM basis. To train the Ansatz Operators in (11), which are used to construct the Closure-Term in (10), we use the training FOM data. Due to the periodic behavior of the flow and in order to decrease the computational cost of constructing the Ansatz Operators, we use a half-period of training FOM data, i.e., 68 FOM snapshots, from t=13 to t=13.134. Then we test all ROMs in the extrapolation testing time interval, t=16 to t=23. In Sect. 5.3.1–5.3.3, we compare the quality of all ROMs based on three criteria: L2-norm error, kinetic energy error, and vortex shedding frequency, respectively. In Table 5, we list the average L2FOM consistency error (26), computed using a reduced basis of d=22 POD modes. This choice of d, which is substantially lower than the full rank of the system, is used consistently in Sect. 5.3.1 and 5.3.2. From the results in Table 5,we observe that both R1-ROM and R2-ROM (for both ansatzes) have lower FOM consistency errors than C-ROM. Furthermore, the FOM consistency errors of R1-ROM and R2-ROM 123 308 Page 18 of 28 B. Koc et al. Fig. 9 ADR equation; ROM solutions at the final time, T=1, for L=8 Fig. 10 Geometry of the flow past a circular cylinder numerical experiment Fig. 11 Flow past a cylinder; L2-POD decay of eigenvalues 123 Residual-based data-driven variational multiscale… Page 19 of 28 308 Table 5 Flow past a cylinder; d=22; average L2FOM consistency error (26)ofC-ROM,R1-ROM,and R2-ROM for various Lvalues (number of POD modes). LC-FOM-Consistency R1-FOM-Consistency R2-FOM-Consistency Ansatz1 Ansatz2 2 2.28e-01 2.59e-02 2.05e-02 7.37e-02 5 1.49e+00 1.72e-02 2.62e-02 3.61e-03 8 4.01e-01 1.22e-04 6.60e-05 3.03e-05 11 1.35e+00 1.27e-08 1.54e-08 3.63e-09 14 2.46e-01 5.39e-03 4.24e-03 2.25e-05 17 7.29e-01 6.19e-03 6.19e-03 1.03e-04 20 1.43e-01 2.34e-03 2.52e-03 8.35e-04 22 0 0 0 0 Table 6 Flow past a cylinder; d=14; average L2FOM consistency error (26)ofC-ROM,R1-ROM,and R2-ROM for various Lvalues (number of POD modes). LC-FOM-Consistency R1-FOM-Consistency R2-FOM-Consistency Ansatz1 Ansatz2 2 2.27e-01 2.29e-02 2.29e-02 7.19e-03 3 2.59e+00 6.39e-04 3.58e-03 1.41e-03 5 1.49e+00 1.63e-02 1.63e-02 1.82e-03 7 7.63e-02 4.28e-02 4.28e-02 6.37e-03 9 2.10e+00 1.50e-02 1.50e-02 4.84e-04 11 1.33e+00 1.77e-02 1.14e-02 6.61e-03 13 1.15e+00 7.90e-03 7.90e-03 2.28e-03 14 0 0 0 0 reach the lowest value at L=11, and then increase again. Since the least squares problem (11) is sensitive based on using the size of the data, to investigate this behavior in Table 5,inTable6, we also list the average L2FOM consistency error (26), calculated using a reduced basis with d=14 POD modes. In Table 6, we observe a similar behavior as that observed in Table 5 but less significant increasing, from L=5. For the flow past a cylinder flow, due to the complex nonlinear interactions in the NavierStokes equations, the consistency error has a more complex behavior than in the linear ADR case. Specifically, the consistency error decreases as the dimension Lof the reduced space increases, up to some value L=Lopt (depending on the total number of modes taken into account to build the ROM, d), and then it approximately stagnates: Lopt =11 when d=22, and Lopt =3whend=14. 5.3.1 L2-norm error In this section, we compare the numerical accuracy of the d2-VMS-ROMs, I-ROM, and G-ROM defined in Sect. 4by using the average relative L2ROM projection error (30). 123 308 Page 20 of 28 B. Koc et al. Table 7 Flow past a cylinder; average relative L2-projection error (30) for G-ROM, I-ROM, C-ROM, R1-ROM, and R2-ROM for various Lvalues (number of POD modes). LG-ROM I-ROM C-ROM R1-ROM R2-ROM 2 1.94e+00 1.33e-02 6.96e-01 6.07e-01 3.05e-01 3 1.54e+00 4.80e-03 9.20e-01 1.03e-01 1.13e-01 4 1.19e+00 4.61e-03 3.26e-01 1.95e-01 9.55e-02 5 1.19e+00 2.15e-03 9.51e-01 5.09e-01 4.03e-01 6 5.76e-01 3.09e-02 2.41e-01 7.10e-02 3.09e-02 7 4.91e-01 2.00e-02 4.56e-01 8.60e-02 3.51e-01 8 2.58e-01 4.23e-03 1.53e-01 4.12e-02 6.40e-02 Fig. 12 Flow past a cylinder; Pareto plot of average relative L2error of d2-VMS-ROMs, i.e., R1-ROM, R2-ROM, and C-ROM In Table 7, we list the average relative L2ROM errors (27) for different ROM dimensions, L. We observe that both R1-ROM and R2-ROM yield more accurate results than C-ROM and G-ROM. Overall, the R2-ROM is more accurate than the R1-ROM, except for L =3,7,8. In Fig. 12, we present a Pareto plot of the d2-VMS-ROMs, i.e., R1-ROM and R2-ROM, and C-ROM, averaging the L2error and offline ansatz cost over the low ROM dimensions, i.e., L =2,3,4,5, and the high-ROM dimensions, i.e., L =6,7,8. For the low-ROM dimension, we observe that the R2-ROM yields the most accurate model, although it is the most expensive model. For the high-ROM dimension L =6,7,8, the R1-ROM is the most accurate and least expensive model. 5.3.2 Kinetic energy error In this section, we compare the numerical accuracy of the d2-VMS-ROMs, I-ROM, and G-ROM defined in Sect. 4by using the kinetic energy (KE) criterion: Ekin := 1 2u2 L2=1 2 |u|2d. (33) The aim of our numerical investigation is to observe how the kinetic energy of ROMs evolves in the extrapolation time interval, e.g. whether it dissipates, blows up, or stagnates. EKE =M k=1|EFOM kin (tk)−EROM kin (tk)| M k=1|EFOM kin (tk)|.(34) In Table 8, we list the average relative kinetic energy error (34) of the ROMs. We observe that R2-ROM generally yields a lower kinetic energy error than R1-ROM, and a much lower kinetic energy error than C-ROM and G-ROM for low ROM dimensions. For high ROM 123 Residual-based data-driven variational multiscale… Page 21 of 28 308 Table 8 Flow past a cylinder; average relative kinetic energy error (34) for G-ROM, I-ROM, C-ROM, R1-ROM, and R2-ROM for various Lvalues (number of POD modes) LG-ROM I-ROM C-ROM R1-ROM R2-ROM 2 4.99e-01 5.91e-03 1.40e-02 1.16e-01 4.34e-03 3 1.62e-01 1.35e-03 1.25e-01 4.57e-03 2.48e-02 4 1.08e-01 2.14e-03 2.22e-02 5.79e-02 4.52e-03 5 8.35e-02 9.35e-04 2.40e-01 1.32e-01 9.41e-02 6 2.14e-01 1.94e-03 1.17e-01 2.74e-02 2.73e-03 7 2.06e-01 2.51e-03 2.08e-01 1.44e-02 1.14e-01 8 7.53e-02 1.73e-03 3.47e-02 8.57e-03 1.31e-02 Fig. 13 Flow past a cylinder; Pareto plot of average relative KE error of d2-VMS-ROMs, i.e., R1-ROM, R2-ROM, and C-ROM dimensions, R1-ROM yields a lower kinetic energy error than R2-ROM, and a much lower kinetic energy error than C-ROM and G-ROM. In Fig. 13, we present the Pareto plot of the d2-VMS-ROMs, i.e., R1-ROM, R2-ROM, and C-ROM, averaging the average relative kinetic energy and offline ansatz cost over the low ROM dimensions, i.e., L =2,3,4,5, and the high ROM dimensions, i.e, L =6,7,8. For the low Lvalues, we observe that R2-ROM yields the most accurate results, although its computational cost is the highest. On the other hand, for the high Lvalues, the R1-ROM is the most accurate and its computational cost is the lowest. The conclusion of the kinetic energy Pareto plot in Fig. 13 is consistent with the conclusion of the average L2Pareto plot in Fig. 12. In Fig. 14, we plot the kinetic energy (33) of the FOM projection, G-ROM, I-ROM, C-ROM, R1-ROM, and R2-ROM for the ROM dimension values L =2,4,6,8 over the extrapolation time interval. These plots support the results of Table 8. 5.3.3 Vortex shedding frequency matching between FOM and d2-VMS-ROMs In this section, we compute the average vortex shedding frequency fsof d2-VMS-ROMs and FOM based on the vortex shedding period Ts, which is defined as follows: Ts=1 Ns Ns  k=1 (ts(k+1)−ts(k)), (35) where ts(k)denotes the time instances corresponding to successive peaks in the kinetic energy within the extrapolation testing time interval [18,23]. These peaks are used to estimate the dominant vortex shedding cycle of the respective models. 123 308 Page 22 of 28 B. Koc et al. Fig. 14 Flow past a cylinder; kinetic energy of FOM projection, G-ROM, I-ROM, C-ROM, R1-ROM, and R2-ROM for L =2,4,6,8 In Table 9, we list the vortex shedding frequency, i.e., fs=1/Ts, of the FOM and d2VMS-ROMs, i.e., C-ROM, R1-ROM, and R2-ROM, for various Lvalues, i.e., L=2,4,6,8. Overall, the results in Table 9show that R2-ROM yields the most accurate predictions of the Strouhal number, especially for low Lvalues. In Table 10, we list the relative vortex shedding frequency errors of d2-VMS-ROMs: Es=|fFOM s−fROM s| fFOM s .(36) Based on the relative errors in vortex shedding frequency listed in Table 10,R2-ROM demonstrates superior accuracy in capturing the vortex shedding period of the FOM among the d2-VMS-ROM variants. 123 Residual-based data-driven variational multiscale… Page 23 of 28 308 Table 9 Flow past a cylinder; vortex shedding frequency, i.e., fsfor FOM, C-ROM, R1-ROM, and R2-ROM for various L values (number of POD modes) LFOM C-ROM R1-ROM R2-ROM 2 3.77e+00 3.74e+00 3.83e+00 3.79e+00 4 3.77e+00 3.79e+00 3.73e+00 3.78e+00 6 3.77e+00 3.82e+00 3.78e+00 3.77e+00 8 3.77e+00 3.78e+00 3.77e+00 3.77e+00 Table 10 Flow past a cylinder; relative vortex shedding frequency error (36)forC-ROM, R1-ROM, and R2-ROM for various Lvalues (number of POD modes) LC-ROM R1-ROM R2-ROM 2 9.14e-03 1.40e-02 3.79e-03 4 4.21e-03 1.04e-02 4.19e-04 6 1.88e-02 8.39e-04 4.19e-04 8 1.26e-03 0 0 Fig. 15 Flow past a cylinder; vortex shedding periods of d2-VMS-ROMs, i.e., R1-ROM, R2-ROM, and CROM for L=2,6 In Fig. 15, we present the sequence of kinetic energy peaks for the FOM and the d2-VMSROMs for L=2andL=6, thus visualizing the vortex shedding periods exhibited by each model. For L=6, in Fig. 15, we observe that both R1-ROM and R2-ROM more accurately recover the vortex shedding periods of the FOM than C-ROM. Furthermore, the amplitudes of the kinetic energy peaks in C-ROM are noticeably higher than those in R1-ROM and R2-ROM. These observations are consistent with the lower ROM errors (see Table 7)and kinetic energy errors (see Table 8) exhibited by R1-ROM and R2-ROM. For L=2, R2-ROM demonstrates the best match to FOM in terms of both vortex shedding frequency and amplitude of kinetic energy peaks. Although the vortex shedding frequencies of C-ROM and R1-ROM are similar, R1-ROM significantly overestimates the peak amplitudes compared to FOM. 123 308 Page 24 of 28 B. Koc et al. Fig. 16 Flow past a cylinder; the velocity of the I-ROM and most accurate d2-VMS-ROM, i.e., R2-ROM, for L=2,4,6atT=23 123 Residual-based data-driven variational multiscale… Page 25 of 28 308 Finally, in Fig. 16, we plot the velocity field of I-ROM and the most accurate d2-VMSROMs, i.e., R2-ROM, for L=2,4,6 at the final time, T=23. We observe that R2-ROM recovers the main features of the flow with good accuracy. 6 Conclusions and outlook In this paper, we proposed a novel residual-based data-driven ROM closure model for underresolved, convection-dominated problems. The new data-driven ROM closure model was constructed by leveraging the VMS framework and the available data. To build the new datadriven VMS-ROM (d2-VMS-ROM), we first postulated a closure model form ansatz that depends on the ROM residual, and then we solved a least squares problem to find the ansatz parameters that yield the closest fit between the ansatz and the FOM data. We also considered two types of residual-based ansatzes, which yielded two types of d2-VMS-ROMs, denoted R1-ROM and R2-ROM. The main novelty of the proposed residual-based d2-VMS-ROMs is that they depend on the ROM residual instead of the ROM coefficients (which is the standard approach in current data-driven ROM closures). To assess the novel residual-based d2-VMS-ROMs, we compared them with the standard coefficient-based d2-VMS-ROM, denoted C-ROM (Mou et al. 2021; Xie et al. 2018). Finally, for comparison purposes, we also investigated a standard G-ROM in which neither stabilization nor closure was used. We investigated the new residual-based d2-VMS ROMs, i.e., R1-ROM and R2-ROM, as well as the standard C-ROM and G-ROM, in the numerical simulation of three test problems: (i) a one-dimensional parameter-dependent advection-diffusion problem (Section 5.1); (ii) a twodimensional time-dependent advection-diffusion-reaction problem with a small diffusion coefficient ε=1e−4 (Section 5.2); and (iii) a two-dimensional flow past a cylinder at Reynolds number Re =1000 (Section 5.3). Our numerical investigation yielded the following conclusions: The novel residual-based d2-VMS-ROMs, R1-ROM and R2-ROM, were significantly more accurate than the standard coefficient-based d2-VMS-ROM, C-ROM. Furthermore, R2-ROM generally yielded slightly more accurate results than R1-ROM, but there were cases in which R1-ROM was more accurate. Finally, all the d2-VMS-ROMs (i.e., R1-ROM, R2-ROM, and C-ROM) were significantly more accurate than the standard G-ROM. Since R1-ROM has a simpler formulation than R2-ROM and is easier to construct, and since the R1-ROM and R2-ROM accuracies are similar, R1-ROM appears preferable to R2-ROM in practice. There are several research directions that can be pursued next. Probably the most important is the extension of the new residual-based d2-VMS-ROM framework to more complex convection-dominated problems, such as under-resolved turbulent flows. Another important research direction is providing mathematical support for the new residual-based d2-VMSROM. For example, we plan to prove the residual-based d2-VMS-ROM’s verifiability, i.e., to show that when the closure model error decreases, the ROM error decreases at the same rate. The first step in this direction has been taken in Koc et al. (2022),whereweprovedthe verifiability of the standard coefficient-based d2-VMS-ROM, C-ROM. Acknowledgements The first author is partially supported by Project PID2021-123153OB-C21 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE and Juan de la Cierva 2022 with project number 2023/1061. The second and third authors are funded by Project PID2021-123153OB-C21 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. The fourth author is funded by ARIA MSCA-RISE EU Grant 872442 and National Science Foundation grant DMS-2012253. 123