Parametric model-order-reduction development for unsteady convection

A time-averaged error indicator with POD-hGreedy is developed to drive parametric model order reduction (pMOR) for 2D unsteady natural convection in a high-aspect ratio slot parameterized with the Prandtl number, Rayleigh number, and slot angle with respect to the gravity. The error indicator is extended to accommodate the energy equation and Leray regularization. Despite being two-dimensional and laminar, the target flow regime presents several challenges: 1) there is a bifurcation in the angle parameter space; 2) the solution can be multivalued, even at steady state; and 3) the solution exhibits spatio-temporal chaos at several points in the parameter space. The authors explore several reduced-order models (ROMs) and demonstrate that Leray-regularized Galerkin ROMs provide a robust solution approach for this class of flows. They further demonstrate that error-indicated pMOR can efficiently predict several QOIs, such as mean flow, mean Nusselt number and mean turbulent kinetic energy, even in the presence of a bifurcation. Finally, they show that spatio-temporal chaos can lead to lack of reproducibility in both the full-order model and the reduced-order model and that the variance in the full-order model provides a lower bound on the pMOR error in these cases.

1 Introduction

Fluid-thermal analysis via direct numerical (DNS), large-eddy (LES), and even unsteady Reynolds-averaged Navier-Stokes (uRANS) have became tractable in geometries of ever increasing complexity due to advances in high-performance computing and modern algorithms over the past several decades. Despite these advances, when it comes to routine analysis and design of thermal hydraulic systems, which requires running hundreds of cases, the cost remains prohibitive. To overcome this issue, a rapid turn-around tool for engineering query is required; parametric model order reduction (pMOR) is one of the promising approaches.

The main idea of pMOR is to reduce the computational burden by employing reduced-order models (ROMs) built on data from full-order models (FOMs) such as DNS, LES, or RANS-based simulations. A fundamental requirement in this case is to determine how well these approaches can reproduce the flow dynamics with same input parameters as the originating FOM, which is known as the reproduction problem. For turbulent flows, FOM can require $\mathcal{N}=10^7-10^{11}$ degrees-of-freedom, while ROMs could potentially represent the flow dynamics which govern the behavior of quantities of interest (QOIs) with only N ≈ 10¹–10² basis functions. Most often, the basis is obtained from a proper orthogonal decomposition (POD) of the FOM flow snapshots (i.e., the velocity and temperature fields). The basis is then used in conjunction with the governing partial differential equation to form a low-dimensional dynamical system that, ideally, captures the principal features of the underlying flows [1–8].

As noted in [9], a successful pMOR for unsteady flows must be able to address 1) the reproduction problem and 2) the parametric problem. In the reproduction problem, the ROM and FOM are evaluated at the same $p_j^{*}\in {\mathcal{P}}_{\text{anchor}}$ and the ability of the ROM to recover the QOIs at $p_j^{*}$ is examined. In the parametric problem, ROMs are built from a small number of FOMs that are generated over a set of anchor points in the parameter space, ${\mathcal{P}}_{\text{anchor}}=\{p_1^{*},\dots ,p_m^{*}\}$, and the ability of ROMs to predict the QOIs at $p\notin {\mathcal{P}}_{\text{anchor}}$ is examined. An error indicator to assess the ROM’s fidelity at any given $p\in {\mathcal{P}}_{\text{train}}$ is usually required to efficiently construct ${\mathcal{P}}_{\text{anchor}}$. Error-indicated pMOR has the potential to be an important engineering analysis tool [10, 11].

In this work we explore the pMOR process for a surprisingly challenging 2D natural convection problem that serves as a surrogate for more difficult 3D buoyancy-driven flows encountered in a variety of mechanical and nuclear engineering applications. The geometry is the tilted slot configuration of Figure 1 and the governing equations are the unsteady Boussinesq equations. The problem is characterized by four parameters, the Rayleigh number, Ra, the Prandtl number, Pr, the slot angle with respect to gravity, θ_g, and the aspect-ratio, Γ. However, for small aspect ratios, Γ ≤ 8, we find the flow is rather simple. Hence in this work, we focus on the more challenging case of Γ = 40 with (Pr, Ra, θ_g) as the parameter space.

FIGURE 1

FIGURE 1. (A) Problem configuration. (B) FOM mean (or steady) temperature for 19 uniformly-spaced θ_g ∈ [0°, 180°] at Ra = 10⁴ and Pr = 7.2. (C) FOM mean temperature for 6 uniformly-spaced Ra ∈ [2 × 10⁵, 7 × 10⁵] at θ_g = 90° and Pr = 7.2.

Figure 1A illustrates the problem configuration. The aspect ratio is defined as Γ = L/H, where H is the width of the slot and L is the height of the slot. We take H = 1 and L = 40 throughout the study. With θ_g = 0° the flow corresponds to standard Rayleigh-Bénard convection, 90° corresponds to vertical slot convection, and θ_g = 180° leads to a pure conduction solution with the cold side on the bottom of the horizontal slot. The Rayleigh and Prandtl numbers are

$\mathrm{R}\mathrm{a}=\frac{\rho \beta \mathrm{\Delta }T{H}^{3}g}{\nu \alpha },\text{Pr}=\frac{\nu }{\alpha },(1)$

where ρ is the fluid density, β is the thermal expansion coefficient, ΔT = 2 is the wall-to-wall temperature difference, H is the gap width, g is the gravitational acceleration, ν is the kinematic viscosity, and α is the thermal diffusivity. Figure 1B illustrates representative mean (or steady) temperature fields for 19 uniformly-spaced θ_g ∈ [0°, 180°] at Ra = 10⁴ and Pr = 7.2. Figure 1C shows the mean temperature field for 6 uniformly-spaced Ra ∈ [2 × 10⁵, 7 × 10⁵] at θ_g = 90° and Pr = 7.2. The configuration has many interesting applications. For example, it represents energy-efficient double-glazed windows, in which the sealed air gap between the two panes acts as an added layer of insulation. Finding the optimum angle θ_g that enhances the heat transfer is an important question. Convection in a tilted fluid layer is also of meteorological and oceanographic interest. More information on the impact of θ_g and Γ on heat transport and flow organization for this configuration can be found in [12].

There has been significant work on pMOR development for the steady Boussinesq equations, including rigorous error estimation [10, 11, 13–15]. For pMOR, the steady problem is easier than the unsteady problem for several reasons: 1) rigorous error estimates are usually achievable, 2) there is often a well defined attractor, and 3) no temporal instability needs to be considered. Once the problem becomes unsteady many open research issues remain. To our knowledge, there are few pMOR works addressing the unsteady parameterized Boussinesq equations. In [16], the authors develop rigorous a posteriori error bounds applied to a 2D Rayleigh-Bénard problem parameterized with Gr and θ_g. However, due to exponential instability in time, the rigor is not for very high Gr and large final times. In [17], the authors overcame the high Gr issue by considering a space-time formulation which enabled effective long-time certification of a reduced basis approximation of noncoercive PDEs. However, the approach is limited due to large offline computational effort since only one snapshot is generated from one FOM solve due to the formulation. For example, to cover the parameter space, 125 FOMs are solved during the offline in their case.

Fick et al. [9] developed a POD-hGreedy pMOR to study challenging incompressible flow using a time-averaged error indicator. The authors showed that the error indicator is highly-correlated with the error in mean flow prediction and can be efficiently computed through an offline/online strategy. We view the methodology as having high potential for routine analysis and design of turbulent flows that are characteristic of thermal hydraulic systems. Hence, we explore that approach here by extending the time-averaged error indicator to accommodate the energy equation. In previous work [18] on ROM stabilization and turbulent thermal transport problems, we investigated the performance of pMOR with constrained stabilization [9], and Leray regularization [19]. Here, we extend the error indicator to support Leray regularization. For each approach, we assess the performance through the mean flow and QOIs including, mean Nusselt number (Nu), standard deviation in Nu, mean temperature fluctuation and mean turbulent kinetic energy (TKE)¹.

Even though our 2D model problem generates only laminar flows, pMOR is quite challenging in this application for several reasons: 1) there is bifurcation in θ_g parameter space; 2) the solution can be multivalued, even at steady state; and 3) the solution exhibits spatio-temporal chaos at several points in the parameter space. As our initial efforts happened to be focused in one of the spatio-temporal chaos regimes, we decided to map out a larger space to identify where pMOR could succeed, where it would have difficulty, and where it a priori could not succeed. Table 1 reflects a broad range of flow regimes identified inside the high-aspect ratio slot from hundreds of FOM simulations conducted at multiple Ra, θ_g with Pr = 0.07, 0.71, 7.2. We categorize the flow into six types: 1) motionless, 2) steady, 3) periodic, 4) quasi-periodic, 5) chaotic and 6) spatio-temporal chaotic. We identify the flow regimes by examining the (mean) solution field and the energy and Nu histories. Such analysis can readily distinguish the motionless, steady and periodic flow cases. Even though the energy and Nu analysis seem to be a reliable way to distinguish the quasi-periodic and chaotic flow, it is only a heuristic—a more rigorous analysis is through computing the power spectrum of Nu or energy [20]. Tools such as Lyapunov exponent and fractal dimension are probably the most widely used diagnostic for chaotic systems [21, 22]. The first five types of low have consistent mean flow in differing time windows, each averaged over 500 convective time units (CTUs). We define a flow to be spatio-temporal chaotic [23] if its mean solution is not consistent in at least three different time windows. This type of flow has strong irregularities in both space and time and has been observed in Rayleigh Bénard convection and in other complex dynamical systems [24]. To characterize spatio-temporal chaos, one could also consider Lyapunov exponents at each grid point. A detail analysis of spatio-temporal chaos is beyond the scope of this paper. A comprehensive review on this topic can be found in [25].

TABLE 1

TABLE 1. Distribution of six flow types with Ra and θ_g at Pr = 0.07, 0.71, 7.2.

In the present work, we start with solution reproduction problem, which represents the first step towards the development of a ROM for the parametric problem. We study the reproduction capability of the FOM, ROM and ROM with stabilization for the six types of flow reported in Table 1. We report only the cases of chaotic and spatio-temporal chaotic flow. Notice that it has been reported that ROM can often capture the first five types of flow accurately, in some cases, with the need of stabilization methods. However, for the spatio-temporal chaotic flow, to our knowledge, it has not been studied. We believe this is the first work investigate ROM’s reproducibility of spatio-temporal chaotic flow.

The pMOR development is broken into several parametric problems. We start with a problem at Ra = 10³ and Pr = 7.2 in which the solution is steady for all θ_g but nonetheless exhibits a bifurcation at θ_g = 20°². We find either h- or p-Greedy with the error indicator based on the dual norm of the residual is able to drive pMOR successfully³. We next consider two sets of parametric problems: 1) problem parameterized with θ_g at higher Ra. 2) Problem parameterized with Ra. In the first set, similar to the steady case, a bifurcation is observed and the solution space is a blend of steady and unsteady solutions. In the second set, no bifurcation is observed and the solutions are all unsteady.

By proceeding in this manner, we are able to isolate several difficulties and eventually come up with an important observation for the pMOR: Accurate prediction $(<10\%)$ with pMOR is achievable if the solution in the parametric space is either only chaotic or the spatio-temporal chaos is not significant, regardless a bifurcation exists or not. Once the spatio-temporal chaos becomes significant, the performance of the pMOR deteriorates and the maximum errors of the mean flow and QOIs are dominated by the flow chaos.

The paper is organized as follows. In Section 2, we introduce the model problem and governing equations. In Section 2.1, we introduce the Galerkin formulation for the FOM. The ROM, as well as Leray regularization is introduced in 2.2. In Section 3, we consider the solution reproduction problem and assess the numerical performance. The parametric problem is discussed starting from Section 4. We first introduce POD-hGreedy algorithms in Section 4.1 with some remarks on applying POD-pGreedy to this model problem at the end of the section. We then introduce the time-averaged error indicator with thermal extension in Section 4.2. A straight-forward integration with Leray regularization is also shown in the same section. In Sections 4.3 and 4.4, we present the pMOR results with θ_g variation and Ra variation. In Section 5, we discuss the spatio-temporal chaos and multiple states issues found in this model problem. Finally, we conclude the paper in Section 6.

2 A parametrized natural convection problem

We start with the Boussinesq equations for buoyancy-driven flow [26],

$\frac{\partial \mathbf{\text{u}}}{\partial t}+\left(\mathbf{\text{u}}\cdot \nabla \right)\mathbf{\text{u}}+\nabla p=\nu {\nabla }^2\mathbf{\text{u}}+T\mathbf{\text{g}}\left({\theta }_g\right),\nabla \cdot \mathbf{\text{u}}=0,(2)$

$\frac{\partial T}{\partial t}+\left(\mathbf{\text{u}}\cdot \nabla \right)T=\alpha {\nabla }^{2}T,(3)$

where u is the velocity, p is the pressure, T is the temperature and g(θ_g) is the unit vector represents the direction of the buoyancy force and it is defined by $\mathbf{\text{g}}({\theta }_g)=\cos ({\theta }_g)\widehat{ı}+\sin ({\theta }_g)\widehat{\mathrm{J}}$, with θ_g the angle of the slot with respect to the gravity. The velocity boundary conditions are no-slip (u = 0). The temperature boundary conditions are no-flux (insulated) on the top and bottom, heated (T = 1) on the left wall, and cooled (T = −1) on the right wall. The initial conditions are u = 0 and T = 0.

In our non-dimensional setting, we set ν = (Pr/Ra)^1/2 and α = (Pr Ra)^−1/2. The Rayleigh number, Ra = ρβgH³ΔT/(να), represents the ratio of buoyancy force to thermal and momentum diffusive force The Prandtl number Pr = ν/α, reflects the relative importance of momentum diffusivity compared to thermal diffusivity. With this nondimensionalization the characteristic velocity is $U_c=\sqrt{\beta gH\mathrm{\Delta }T}$, the characteristic length is the slot width H, and the reference time is t_r = H/U_c. The temperature is made dimensionless by subtracting the temperature on the right wall and scaling with ΔT = 2. We note that U_c is sometimes referred to as the “free-fall” velocity, indicating that one might expect ‖u‖ ≈ 1, with only a weak dependence on Ra. While that expectation is realized for θ_g = 0, we in fact see much larger velocities (‖u‖ ≈ 40) for θ_g = 90° because the domain height L = 40H in that case.

For unsteady problems, the QOIs are the mean flow, mean Nu, standard deviation in Nu, mean TKE and mean temperature fluctuation. The symbol ⟨⋅⟩ is used to indicate a time-averaged quantity. The mean velocity and temperature field are defined as:

$⟨\mathbf{\text{u}}⟩=\frac{1}{J-J_0}\sum _{j=J_0+1}^J\mathbf{\text{u}}\left(t^j\right),⟨T⟩=\frac{1}{J-J_0}\sum _{j=J_0+1}^{J}T\left(t^j\right),(4)$

with t^j = jΔt and Δt being the time step. The selection of J₀ is based on when the solution reaches it statistically steady state. The mean quantities are then averaged over 500 CTUs, with the time scale defined above. The instantaneous Nusselt number is defined as

$\text{Nu}\left(t\right)=\frac{q_w^{\prime \prime }}{k\left(\mathrm{\Delta }T\right)/H},(5)$

with $q_w^{\prime \prime }=-{\int }_{\partial {\mathrm{\Omega }}_h}k\nabla T\cdot \widehat{\mathbf{\text{n}}}dS$ being the integrated heat flux on the heated wall, ∂Ω_h. The mean Nu and the its standard deviation are then defined as

$⟨\text{Nu}⟩:=\frac{1}{J-J_0}\sum _{j=J_0+1}^J\text{Nu}\left(t^j\right),\text{Std}\left(\text{Nu}\right):=\sqrt{\frac{1}{J-J_0}\sum _{j=J_0+1}^J\left(\text{Nu}\left(t^j\right)-⟨\text{Nu}⟩\right)}.(6)$

The mean TKE and mean temperature fluctuation are defined as

$\langle \text{TKE}\rangle =\frac{1}{2\left(J-J_0\right)}\sum _{j=J_0+1}^J\Vert \mathbf{\text{u}}\left(t^j\right)-⟨\mathbf{\text{u}}⟩{\Vert }_{L^2}^2,\langle T_{\text{fluc}}\rangle =\frac{1}{J-J_0}\sum _{j=J_0+1}^J\Vert T\left(t^j\right)-⟨T⟩{\Vert }_{L^2}^2.(7)$

For steady problems, the QOIs are simply the steady solutions to Eqs. 2 and 3 and the corresponding Nu using Eq. 5.

2.1 Galerkin formulation for the full-order model

The FOM is constructed through the spectral element method (SEM) and the P_q–P_q−2 velocity-pressure coupling [27], where the velocity is represented as a tensor-product Lagrange polynomial of degree q in the reference element $\widehat{\mathrm{\Omega }}≔{[-\mathrm{1,1}]}^2$ while the pressure is of degree q − 2. The solution in Ω = ⋃_eΩ^e consists of local representations of u, p, and T that are mapped from $\widehat{\mathrm{\Omega }}$ to Ω^e for each element, e = 1, … , E. In the current FOMs for the slot problem we use E = 516 elements (an array of 6 × 86 in the H × L directions), of order q = 9, for a total of $\mathcal{N}\approx 42000$ grid points. The FOM simulations are performed using the open-source code Nek5000 [28].

For any u(x, t), we have a corresponding vector of basis coefficients $\underline{\mathbf{\text{u}}}={[{\mathbf{\text{u}}}_1\dots {\mathbf{\text{u}}}_{\mathcal{N}}]}^T$ such that

$\mathbf{\text{u}}\left(\mathbf{\text{x}},t\right)=\sum _{j=1}^{\mathcal{N}}{\mathbf{\text{u}}}_j\left(t\right){\varphi }_j\left(\mathbf{\text{x}}\right)\in {\mathbf{\text{X}}}_0^{\mathcal{N}}\subset {\mathcal{\{H\}}}_0^1,(8)$

with ϕ_j(x) the underlying spectral element basis functions spanning the FOM approximation space, ${\mathbf{\text{X}}}_0^{\mathcal{N}}$. Because the SEM is nodal-based, each u_j(t) represents the two velocity components at grid point x_j in the spectral element mesh at time t. Similarly, the temperature is given by

$T\left(\mathbf{\text{x}},t\right)=\sum _{j=1}^{\bar{\mathcal{N}}}{T}_j\left(t\right){\varphi }_j\left(\mathbf{\text{x}}\right)\in X_0^{\mathcal{N}}\subset {\mathcal{H}}_0^1.(9)$

Here, ${\mathcal{H}}^1$ is the set of square-integrable functions on Ω whose gradient is also square-integrable and $X^{\mathcal{N}}\subset {\mathcal{H}}^1$ is the finite dimensional SEM approximation space spanned by {ϕ_j(x)}. ${\mathcal{H}}_0^1$ is the set of functions in ${\mathcal{H}}^1$ that vanish wherever Dirichlet conditions associated with Eq. 3 are applied on the domain boundary ∂Ω and ${\mathcal{H}}_b^1$ is the set of functions in ${\mathcal{H}}^1$ that satisfy the prescribed Dirichlet conditions for temperature. Bold-face indicates that the space is spanned by vector-valued functions having d components (d = 2 or 3) and, in the case of ${\mathbf{\text{X}}}_0^{\mathcal{N}}\subset {\mathcal{\{H\}}}_0^1$, that the functions vanish where Dirichlet conditions are applied for Eq. 2. The pressure p is in $Y^{\mathcal{N}}\subset {\mathcal{L}}^2(\mathrm{\Omega })$, which is the space of piecewise continuous functions on Ω such that ∫_Ωp² dx < ∞. For convenience, we denote ${\overline{\mathbf{\text{Z}}}}^{\mathcal{N}}=({\mathbf{\text{X}}}^{\mathcal{N}},Y^{\mathcal{N}},X^{\mathcal{N}})$ as the collection of the relevant finite-dimensional spaces and will add a subscript 0 or b where required to explicitly indicate the imposed boundary conditions.

Both the FOM and ROM are cast within the same Galerkin framework. To begin, we introduce several inner products for elements in the FOM space, ${\overline{\mathbf{\text{Z}}}}^{\mathcal{N}}$. For any pair of scalar fields $(p,q)\in {\mathcal{L}}^2(\mathrm{\Omega })$ and d-dimensional vector fields, v(x) = [v₁(x) … v_d(x)], u(x) = [u₁(x) … u_d(x)] whose components are also in ${\mathcal{L}}^2$, let

$\left(q,p\right):={\int }_{\mathrm{\Omega }}qpd\mathbf{\text{x}},\left(\mathbf{\text{v}},\mathbf{\text{u}}\right):={\int }_{\mathrm{\Omega }}\left(v_1{u}_1+\cdots +v_d{u}_d\right)d\mathbf{\text{x}}.(10)$

Further, for S, $T\in X^{\mathcal{N}}$ and v, $\mathbf{\text{u}}\in {\mathbf{\text{X}}}^{\mathcal{N}}$, let

$a\left(S,T\right):=\left(\nabla S,\nabla T\right),a\left(\mathbf{\text{v}},\mathbf{\text{u}}\right):=\left(\nabla \mathbf{\text{v}},\nabla \mathbf{\text{u}}\right),(11)$

$c\left(S,\mathbf{\text{u}},T\right):=\left(S,\mathbf{\text{u}}\cdot \nabla T\right),c\left(\mathbf{\text{v}},\mathbf{\text{u}},\mathbf{\text{w}}\right):=\left(\mathbf{\text{v}},\mathbf{\text{u}}\cdot \nabla \mathbf{\text{w}}\right).(12)$

For the FOM, we consider the (semi-discrete) weak form of Eqs. 2 and 3 [29], Find $(\overline{\mathbf{\text{u}}},p,\bar{T})\in {\overline{\mathbf{\text{Z}}}}_b^{\mathcal{N}}$ such that, for all $(\mathbf{\text{v}},q,S)\in {\overline{\mathbf{\text{Z}}}}_0^{\mathcal{N}}$,

$\left(\mathbf{\text{v}},\frac{\partial \overline{\mathbf{\text{u}}}}{\partial t}\right)+\nu a\left(\mathbf{\text{v}},\overline{\mathbf{\text{u}}}\right)-\left(\nabla \cdot \mathbf{\text{v}},p\right)=-c\left(\mathbf{\text{v}},\overline{\mathbf{\text{u}}},\overline{\mathbf{\text{u}}}\right)+\left(\mathbf{\text{v}},\mathbf{\text{g}}\left({\theta }_g\right)\bar{T}\right),(13)$

$-\left(q,\nabla \cdot \overline{\mathbf{\text{u}}}\right)=0,(14)$

$\left(S,\frac{\partial \bar{T}}{\partial t}\right)+\alpha a\left(S,\bar{T}\right)=-c\left(S,\overline{\mathbf{\text{u}}},\bar{T}\right).(15)$

Here, we have introduced $\overline{\mathbf{\text{u}}}=\mathbf{\text{u}}+{\mathbf{\text{u}}}_0(\mathbf{\text{x}})$ and $\bar{T}=T+T_0(\mathbf{\text{x}})$ as functions that have been augmented by (potentially trivial) lifting functions, u₀ and T₀, which are functions of space only. If these functions satisfy the (time-independent) boundary conditions, then one can account for inhomogeneous boundary conditions by moving them to the right-hand side. In the case of the ROM, the lifting functions can also provide an initial approximation to the solution. In the sequel, our principal unknowns will be u and T.

Following [30], we consider a semi-implicit scheme BDFk/EXTk to discretize Eqs. 13–15 in time; kth-order backward differencing (BDFk) is used for the time-derivative term, kth-order extrapolation (EXTk) is used for the advection and buoyancy terms and implicit treatment on the dissipation terms. As discussed in [30], k = 3 is used to ensure the imaginary eigenvalues associated with skew-symmetric advection operator are within the stability region of the BDFk/EXTk time-stepper. Denoting the solution at time tⁿ = Δt ⋅ n as $({\overline{\mathbf{\text{u}}}}^n,p^n,{\bar{T}}^n)$, the full discretization of the FOM reads Find $({\mathbf{\text{u}}}^n,p^n,T^n)\in {\overline{\mathbf{\text{Z}}}}_0^{\mathcal{N}}$ such that, for all $(\mathbf{\text{v}},q,S)\in {\overline{\mathbf{\text{Z}}}}_0^{\mathcal{N}}$,

$\frac{{\beta }_0}{\mathrm{\Delta }t}\left(\mathbf{\text{v}},{\mathbf{\text{u}}}^n\right)+\nu a\left(\mathbf{\text{v}},{\mathbf{\text{u}}}^n\right)-\left(\nabla \cdot \mathbf{\text{v}},p^n\right)=\left(\mathbf{\text{v}},{\mathbf{\text{f}}}^n\right),(16)$

$-\left(q,\nabla \cdot {\mathbf{\text{u}}}^n\right)=\left(q,\nabla \cdot {\mathbf{\text{u}}}_0\right),(17)$

$\frac{{\beta }_0}{\mathrm{\Delta }t}\left(S,T^n\right)+\alpha a\left(S,T^n\right)=\left(S,Q^n\right).(18)$

Equations 16–18 represent a linear unsteady Stokes plus unsteady heat equation to be solved at each time-step tⁿ. The inhomogeneous terms comprise the BDF, advection, buoyancy and lifting terms

$\left(\mathbf{\text{v}},{\mathbf{\text{f}}}^n\right):=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}\left(\mathbf{\text{v}},{\mathbf{\text{u}}}^{n-s}\right)+{\alpha }_s\left(c\left(\mathbf{\text{v}},{\overline{\mathbf{\text{u}}}}^{n-s},{\overline{\mathbf{\text{u}}}}^{n-s}\right)-\left(\mathbf{\text{v}},\mathbf{\text{g}}\left({\theta }_g\right){\bar{T}}^{n-s}\right)\right)\right]-\nu a\left(\mathbf{\text{v}},{\mathbf{\text{u}}}_0\right),(19)$

$\left(S,Q^n\right):=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}\left(S,T^{n-s}\right)+{\alpha }_{s}c\left(S,{\overline{\mathbf{\text{u}}}}^{n-s},{\bar{T}}^{n-s}\right)\right]-\alpha a\left(S,T_0\right).(20)$

Here, the β_ss and α_ss are the respective sth-order BDF and extrapolation coefficients for the BDFs/EXTs time-stepper [30]. Note that the right-hand side of Eq. 17 will be zero if u₀ is divergence free or at least satisfies the weak divergence-free condition Eq. 14.

Under the assumption that ∇ ⋅u₀ = 0, the compact matrix form [27, 31, 32] for Eqs. 16–20 is

$\left[\begin{array}{cc}\hfill \mathbf{\text{H}}\hfill & \hfill -{\mathbf{\text{D}}}^T\hfill \\ \hfill -\mathbf{\text{D}}\hfill & \hfill 0\hfill \end{array}\right]\left(\begin{array}{c}\hfill {\underline{\mathbf{\text{u}}}}^n\hfill \\ \hfill {\underline{p}}^n\hfill \end{array}\right)=\left(\begin{array}{c}\hfill \widehat{\underline{\mathbf{\text{f}}}}\left({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n;{\theta }_g\right)\hfill \\ \hfill 0\hfill \end{array}\right),(21)$

$H_{\alpha }{\underline{T}}^n=\widehat{\underline{Q}}\left({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n\right).(22)$

Here, ${\underline{\mathbf{\text{u}}}}^n$, ${\underline{p}}^n$, and ${\underline{T}}^n$ are the vectors of spectral element basis coefficients. The corresponding block matrices are

$\mathbf{\text{H}}=\left[\begin{array}{cc}\hfill H_{\nu }\hfill & \hfill \hfill \\ \hfill \hfill & \hfill H_{\nu }\hfill \end{array}\right],\mathbf{\text{D}}=\left[D_1{D}_2\right],(23)$

with $H_{\nu }=\frac{{\beta }_0}{\mathrm{\Delta }t}M+\nu A$ and $H_{\alpha }=\frac{{\beta }_0}{\mathrm{\Delta }t}M+\alpha A$, with matrices M and A defined below. The velocity data vectors are $\widehat{\underline{\mathbf{\text{f}}}}({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n;{\theta }_g)={[{\widehat{\underline{f}}}_1{\widehat{\underline{f}}}_2]}^T$, with

${\widehat{\underline{f}}}_m:=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}M{\underline{u}}_m^{n-s}+{\alpha }_s\left(C\left({\underline{\overline{\mathbf{\text{u}}}}}^{n-s}\right){\underline{\bar{u}}}_m^{n-s}-g_{m}M{\underline{\bar{T}}}^{n-s}\right)\right]-\nu A{\underline{u}}_{m,0},m=\mathrm{1,2}.(24)$

where g₁ = cos(θ_g) and g₂ = sin(θ_g) represent the parametric forcing. The thermal load in Eq. 22 is

$\widehat{\underline{Q}}\left({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n\right):=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}M{\underline{\bar{T}}}^{n-s}+{\alpha }_{s}C\left({\underline{\overline{\mathbf{\text{u}}}}}^{n-s}\right){\underline{\bar{T}}}^{n-s}\right]-\alpha A{\underline{T}}_0.(25)$

Entries of the respective stiffness, mass, convection, and gradient matrices are

$A_{ij}={\int }_{\mathrm{\Omega }}\nabla {\varphi }_i\nabla {\varphi }_{j}d\mathbf{\text{x}},(26)$

$M_{ij}={\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_{j}d\mathbf{\text{x}},(27)$

$C_{ij}\left(\mathbf{\text{w}}\right)={\int }_{\mathrm{\Omega }}{\varphi }_i\cdot \left(\mathbf{\text{w}}\cdot \nabla \right){\varphi }_{j}d\mathbf{\text{x}},(28)$

$D_{m,ij}={\int }_{\mathrm{\Omega }}{\psi }_i\frac{\partial {\varphi }_j}{\partial x_m}d\mathbf{\text{x}},m=\mathrm{1,2}.(29)$

Note that {ϕ_i(x)} forms the spectral element velocity/temperature basis while {ψ_i(x)} constitutes the pressure basis.

2.2 Galerkin formulation for the reduced-order model

Within the Galerkin framework of the preceding section it is relatively straightforward to develop a ROM. One defines a set of functions ${\mathit{\zeta }}_j(\mathbf{\text{x}})\in {\mathbf{\text{X}}}^N\subset {\mathbf{\text{X}}}^{\mathcal{N}}$, ${\theta }_j(\mathbf{\text{x}})\in X^N\subset X^{\mathcal{N}}$ such that the coarse-space (ROM) solution, is expressed as

${\overline{\mathbf{\text{u}}}}_{\text{c}}\left(\mathbf{\text{x}}\right)=\sum _{j=0}^N{\mathit{\zeta }}_j\left(\mathbf{\text{x}}\right)u_{\text{c},j},{\bar{T}}_{\text{c}}\left(\mathbf{\text{x}}\right)=\sum _{j=0}^N{\theta }_j\left(\mathbf{\text{x}}\right)T_{\text{c},j}.(30)$

For the ROM, we insert the expansions Eq. 30 into Eqs. 13–15 and require equality for all (v, S) in ${\mathbf{\text{Z}}}_0^N$. In order to set the boundary conditions, we have augmented the trial (approximation) spaces X^N and X^N with the lifting function ζ₀: = u₀ and θ₀: = T₀. The corresponding test spaces, ${\mathbf{\text{X}}}_0^N:={\{{\mathit{\zeta }}_j\}}_{j=1}^N$ and $X_0^N:={\{{\theta }_j\}}_{j=1}^N$, satisfying homogeneous boundary conditions, as is standard for Galerkin formulation. The coarse space ${\mathbf{\text{Z}}}_0^N≔({\mathbf{\text{X}}}_0^N,X_0^N)$ is typically based on a linear combination of full spectral element solutions of Eqs. 21 and 22, such as snapshots at certain time-points, tⁿ, or solutions at various parametric values. Under these conditions and with a carefully chosen u₀, X^N is a set of velocity-space functions that are (weakly) divergence-free and the pressure terms drop out of Eqs. 13 and 14. We also note that ζ_j and θ_j are modal, not nodal, basis functions. In this work, we consider proper orthogonal decomposition (POD) to construct the reduced-basis. The N-dimensional POD-space is the space that minimizes the averaged distance between the snapshot set and the N-dimensional subspace of the snapshots set in the H¹ semi-norm. Further details of the POD basis selection are provided in [18] and references therein.

The matrix form for the ROM is readily derived by constructing a pair of rectangular basis matrices, B and B, having entries

${\mathbf{\text{B}}}_{ij}={\mathit{\zeta }}_j\left({\mathbf{\text{x}}}_i\right),B_{ij}={\theta }_j\left({\mathbf{\text{x}}}_i\right),(31)$

where the x_is are the spectral element nodal points. The coarse-system matrices are H_c,ν = B^THB and H_c,α = B^TH_αB and the governing system for the ROM becomes

$H_{\text{c},\nu }{\widehat{\underline{u}}}_{\text{c}}^n={\mathbf{\text{B}}}^T\widehat{\underline{\mathbf{\text{f}}}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n;{\theta }_g\right),H_{\text{c},\alpha }{\widehat{\underline{T}}}_{\text{c}}^n=B^T\widehat{\underline{Q}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n\right).(32)$

We refer Eq. 32 as Galerkin ROM (G-ROM) throughout the paper. The ROM coefficient vectors, ${\widehat{\underline{u}}}_{\text{c}}^n={[u_{\text{c},1}^n\dots u_{\text{c},N}^n]}^T$, ${\widehat{\underline{T}}}_{\text{c}}^n={[T_{\text{c},1}^n\dots T_{\text{c},N}^n]}^T$, are determined by solving the 2 N × N linear systems. Note that the coefficients for the lifting functions are prescribed: u_c,0 = T_c,0 = 1. The initial coefficients for the ROM are obtained by projecting the initial condition onto the coarse space ${\mathbf{\text{Z}}}_0^N$ with the H¹ semi-norm,

${\widehat{\underline{u}}}_{\text{c}}^0={\mathbf{\text{B}}}^T\mathbf{\text{A}}{\underline{\mathbf{\text{u}}}}^0,{\widehat{\underline{T}}}_{\text{c}}^0=B^{T}A{\underline{T}}^0,(33)$

which follows from the fact that the columns of B and B are, respectively, A- and A-orthonormal, where A = block-diagonal(A). To recover the spectral element representation, we simply prolong the N-length vectors ${\widehat{\underline{u}}}_{\text{c}}^n$ and ${\widehat{\underline{T}}}_{\text{c}}^n$ with the set of basis functions and add it with the lifting functions u₀ and T₀

${\underline{\overline{\mathbf{\text{u}}}}}_{\text{c}}^n=\mathbf{\text{B}}{\widehat{\underline{u}}}_{\text{c}}^n+{\underline{\mathbf{\text{u}}}}_0,{\underline{\bar{T}}}_{\text{c}}^n=B{\widehat{\underline{T}}}_{\text{c}}^n+{\underline{T}}_0.(34)$

The functional representations, ${\overline{\mathbf{\text{u}}}}_{\text{c}}^n(\mathbf{\text{x}})$ and ${\bar{T}}_{\text{c}}^n(\mathbf{\text{x}})$, are then obtained from Eqs. 8 and 9.

Next we consider the Galerkin ROM with Leray regularization using a spatial filter, following ideas presented in [19]. The approach simply requires regularizing the advecting field in the Navier-Stokes equations and energy equation through a low-pass filter function F (i.e., ${\underline{\overline{\mathbf{\text{u}}}}}_{\text{filtered}}=F(\underline{\overline{\mathbf{\text{u}}}})$). As noted in [33, 34], a small amount of regularization is sufficient to make gains in proving existence and uniqueness of the solution to the Navier-Stokes equations. Thus, Leray regularization is of interest both from a numerical (and physical) stabilization viewpoint and from a theoretical perspective.

The formulation of G-ROM with Leray regularization is shown in Eq. 35 and the only difference comparing to G-ROM Eq. 32 are the velocity data and thermal load.

$H_{\text{c},\nu }{\widehat{\underline{u}}}_{\text{c}}^n={\mathbf{\text{B}}}^T{\widehat{\underline{\mathbf{\text{f}}}}}_{\text{filtered}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n;{\theta }_g\right),H_{\text{c},\alpha }{{\widehat{\underline{T}}}_{\text{c}}}^n=B^T{\widehat{\underline{Q}}}_{\text{filtered}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n\right),(35)$

where ${\widehat{\underline{\mathbf{\text{f}}}}}_{\text{filtered}}({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n;{\theta }_g)={[{\widehat{\underline{f}}}_1{\widehat{\underline{f}}}_2]}_{\text{filtered}}^T$. We use the subscript filtered to denote the advecting field in the velocity data and thermal load are being filtered,

${\widehat{\underline{f}}}_{\text{filtered},m}:=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}M{\underline{u}}_m^{n-s}+{\alpha }_s\left(C\left({\underline{\overline{\mathbf{\text{u}}}}}_{\text{filtered}}^{n-s}\right){\underline{\bar{u}}}_m^{n-s}-g_{m}M{\underline{\bar{T}}}^{n-s}\right)\right]-\nu A{\underline{u}}_{m,0},m=\mathrm{1,2}.(36)$

and

${\widehat{\underline{Q}}}_{\text{filtered}}\left({\overline{\mathbf{\text{u}}}}^n,{\bar{T}}^n\right):=-\sum _{s=1}^k\left[\frac{{\beta }_s}{\mathrm{\Delta }t}M{\underline{\bar{T}}}^{n-s}+{\alpha }_{s}C\left({\underline{\overline{\mathbf{\text{u}}}}}_{\text{filtered}}^{n-s}\right){\underline{\bar{T}}}^{n-s}\right]-\alpha A{\underline{T}}_0.(37)$

In this work, we will be focusing on a PDE- (or differential-) based filter, which is characterized by a filter width, δ [35]. Following [19, 36], such filters are developed in a POD context as follows: Find ${\overline{\mathbf{\text{u}}}}_{\text{c},\text{filter}}\in {\mathbf{\text{X}}}^N$ such that

$\left(\left(I-{\delta }^2{\nabla }^2\right){\overline{\mathbf{\text{u}}}}_{\text{c},\text{filter}}^n,{\mathit{\zeta }}_j\right)=\left({\overline{\mathbf{\text{u}}}}_{\text{c}},{\mathit{\zeta }}_j\right),\forall j=1,\dots ,N.(38)$

Besides the differential filter, one could also consider a more economic spatial filter, namely, a POD-projection (Proj) filter as discussed in [19]. In this case, one simply truncates the higher POD mode contributions when constructing ${\overline{\mathbf{\text{u}}}}_{\text{c},\text{filtered}}^n$, just as one would do in a Fourier reconstruction [18].

Besides G-ROM and LDF-ROM, we also consider the constrained-evolution stabilization introduced in [9]. The idea behind this approach is to use information from the snapshot set to establish a priori limits on the ROM coefficients ${\widehat{\underline{u}}}_{\text{c}}$ by replacing Eq. 32 with a constrained minimization problem. At each time-step, the coefficients satisfy

${\widehat{\underline{u}}}_{\text{c}}^n=\arg \underset{{\widehat{\underline{u}}}_{\text{c}}\in {\mathcal{R}}^N}{\min }\frac{1}{2}\Vert H_{\text{c},\nu }{\widehat{\underline{u}}}_{\text{c}}-{\mathbf{\text{B}}}^T\widehat{\underline{\mathbf{\text{f}}}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n;{\theta }_g\right){\Vert }_{H^{-1}}^2,\text{s.t.}{m}_j\le u_j^n\le M_j.(39)$

where the constraints m_j and M_j on the basis coefficients $u_{\text{c},j}^n$, j = 1, … , N are derived from the observation snapshot set. A constrained minimization problem for the thermal ROM coefficients, ${\widehat{\underline{T}}}_{\text{c}}$, is derived similarly. We denote Eq. 39 as C-ROM. Further implementation details can be found in [9, 18].

3 The solution reproduction problem

In this section, we consider the solution reproduction problems for Ra = 2 × 10⁵, 7 × 10⁵ at Pr = 7.2 and θ_g = 90°, where the solutions are chaotic. We assess the performance of G-ROM, C-ROM and LDF-ROM introduced in Section 2.1 through the accuracy of the mean field and the QOIs. The mean field is computed by averaging the POD coefficients and reconstructing with the rectangular basis matrices B and B. The QOIs are the mean Nu and Std(Nu), which are estimated through Eq. 6, and the mean TKE and mean temperature fluctuation, estimated through Eq. 7. The quantities are averaged over 500 CTUs.

Although of limited practical interest, the solution reproduction problem is an important step towards the development of a MOR procedure for the parametric problem. The reproduction results for spatio-temporal chaotic flow are presented and discussed in Section 5.1.

3.1 Numerical results

Results for Ra = 2 × 10⁵ are shown in Figures 2A–D. The performances of G-ROM, C-ROM and LDF-ROM are indicated by blue, orange, green solid line. In LDF-ROM, the radius of the differential filter is δ = 0.015625. Figure 2A shows the behavior of the relative ${\mathcal{H}}^1$ error in the mean flow prediction versus N. We observe less than 1% error in C-ROM with small N. The error in G-ROM and LDF-ROM decreases as N increases and eventually reaches 1% with N = 100. Similar trends for the mean Nu, Std(Nu) and mean temperature fluctuation are observed in Figures 2B–D. The FOM data is denoted as black solid line. For mean Nu prediction, we observe around 0.1% error in C-ROM for almost all N. The error in G-ROM and LDF-ROM decreases as N increases and LDF-ROM has error around 0.01% with N = 100. For Std(Nu) and mean temperature fluctuation prediction, we observe convergence in both QOIs with G-ROM and LDF-ROM. For C-ROM, the predictions are in good agreement with the FOM data for all values of N.Behaviors of the same quantities for Ra = 7 × 10⁵ are shown in Figures 2E–H. In LDF-ROM, the radius of the differential filter is δ = 0.03125. We still observe convergence in the mean flow and QOIs but it is much slower because the flow is more chaotic. The error in the mean flow is around 10% in G-ROM and the prediction in QOIs are over-estimated and less accurate than for the other methods. With LDF-ROM, the predictions are slightly better. C-ROM is still the most effective and has 5% error in mean flow prediction. Although the prediction in mean Nu is only as good as the LDF-ROM, the Std(Nu) and mean temperature fluctuation are bounded from above for all values of N and converge to the correct value as N increases.

FIGURE 2

FIGURE 2. Performance comparison between G-ROM, C-ROM and LDF-ROM for the solution reproduction problem at θ_g = 90° and Pr = 7.2. (A–D) Behavior of the relative ${\mathcal{H}}^1$ error in mean flow prediction, predicted mean Nu, Std(Nu), and mean temperature fluctuation as a function of number of modes, N at Ra = 2 × 10⁵ (δ = 0.015625 in LDF-ROM). (E–H) Performance for the same quantities at Ra = 7 × 10⁵ (δ = 0.03125 in LDF-ROM).

Note that the differential filter radius δ selected for the two Ra yields the best accuracy in mean flow among the five differential filter radius δ = 0.25, 0.125, 0.0625, 0.03125, 0.015625. Besides, we find δ = 0.015625 yield the best results at smaller Ra and as Ra increases, results with δ = 0.03125 becomes better and are comparable with δ = 0.015625 at Ra = 7 × 10⁵. The tendency is reasonable since the flow is more chaotic as Ra increases, therefore one should expect a larger δ to stabilize the flow.

From the results, we observe the mean flow and QOIs converge with N for all ROMs. With higher Ra, convergence in those quantities is much slower and a larger N is required to reach to the same accuracy as in the lower Ra case. Note that, because of the O(N³) online costs, requiring a large value of N for convergence might require off-line resources for timely simulation, which would greatly diminish the intrinsic advantage of the ROM/pMOR framework. This potentiality highlights the importance of stabilization methods. Indeed, we find that C-ROM is able to predict the mean flow and QOIs with a better accuracy with smaller N. On the other hand, although LDF-ROM is not as effective as C-ROM, and only slightly better than G-ROM, it will play a key role in the pMOR presented in Section 4, especially in the parametric problem with θ_g variation.

4 The parametric problem

We turn now to study the performance of pMOR for the slot problem with G-ROM, C-ROM, and LDF-ROM. Two sets of parameterization are considered. With Pr = 7.2 fixed, we seek to estimate the solution and QOIs of Eqs. 13 and 15 for: 1) ${\theta }_g\in \mathcal{P}=[0°,180°]$ at multiple Ra, and 2) $\mathrm{R}\mathrm{a}\in \mathcal{P}=[2\times {10}^5,7\times {10}^5]$ with θ_g = 90°. Throughout, we take the lifting functions to be the zero velocity field and the heat conduction solution.

For efficient selection of the pMOR anchor points we consider the POD-hGreedy algorithm proposed in [9] which combines POD in time with Greedy in parameter. The term Greedy refers to the optimization strategy of basing anchor point selection on the training point that exhibits the largest value in the error estimate. Error-indicated selection of the anchor points reduces the number of FOM solves and is thus critical for the feasibility of pMOR. Here, the error indicator corresponds to the dual norm of the residual associated with the time-averaged momentum and energy equations.

The section is organized as follows. In Section 4.1 we present the POD-hGreedy algorithm. In Section 4.2, we extend the time-averaged error indicator introduced in [9] to accommodate the energy equation and Leray regularization. Finally, in Sections 4.3 and 4.4, we present the numerical results.

4.1 Proper orthogonal decomposition-hGreedy algorithm

In this section, we present the POD-hGreedy algorithm for the construction of the reduced spaces ${\{{\mathbf{\text{X}}}_{0,\ell }^N,X_{0,\ell }^N\}}_{\ell =1}^L$, and the partition ${\{{\mathcal{I}}_{\ell }\}}_{\ell =1}^L$ of $\mathcal{P}$, based on the results of L full-order simulations associated with the parameters $p_1^{*},\dots ,p_L^{*}$. The algorithm is similar to the one in [9] but with extensions for thermal fields.

To begin, we introduce the discretized parameter space ${\mathcal{P}}_{\text{train}}={\{p_i\}}_{i=1}^{n_{\text{train}}}$, $p_1\le \dots \le p_{n_{\text{train}}}$, the integers L which fix the maximum number of offline solves, the integer n_cand < L, which is the number of ROM evaluations performed online for a given value of the parameters, and an error indicator Δ. The error indicator takes as input sequences ${\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\}}_{n=0}^J$ and ${\{{\bar{T}}_{\text{c}}^n\}}_{n=0}^J$ and the value of the parameter, p*, and returns an estimate of the error in the prediction of the mean flow. We formally present the indicator in Section 4.2.

Algorithm 1 presents the computational procedure for both offline and online stage. The offline procedure starts with an anchor point that could either be selected randomly from the training space ${\mathcal{P}}_{\text{train}}$ or be user specified. At each iteration ℓ, a FOM simulation at the anchor point $p_{\ell }^{*}$ is conducted and returns a set of snapshots. The snapshot set is then processed through POD and returns the first N orthonormalized POD modes. The value of N is determined by reproduction problem at the anchor point. The ROM and the error indicator, Δ_ℓ, are then built with the reduced spaces ${\mathbf{\text{X}}}_{0,\ell }^N$ and $X_{0,\ell }^N$. The coefficients and the error estimates are then computed for each $p\in {\mathcal{P}}_{\text{train}}$ and the next anchor point is identified as the parameter that has the maximum value in the current (including previous) error estimate. The procedure starts again with the new selected anchor point. If the error indicator is sufficiently small over all points in ${\mathcal{P}}_{\text{train}}$ or the procedure reaches the maximum number of FOM solves L, the offline stage terminates.

Given the ROM/anchor point data $({\mathbf{\text{X}}}_{0,\ell }^N,X_{0,\ell }^N,p_{\ell }^{*})$ for ℓ = 1, … , L, and error indicator, Δ, the hGreedy online stage starts with finding the n_cand candidate anchor points nearest to the test parameter p. The ROMs associated to the candidate anchor points are then used to compute the coefficients and error estimate at p. The coefficients are then returned based on the ROM that has smallest error estimate. The POD-hGreedy approach is analogous to h-refinement in the finite element method in that the POD bases are not shared between anchor points. Convergence is therefore expected to be linear in the distance from the nearest anchor point.

Algorithm 1. POD-hGreedy algorithm for the construction of ${\{{\mathbf{\text{X}}}_{0,\ell }^N,X_{0,\ell }^N,{\mathcal{I}}_{\ell }\}}_l$.

Another strategy is the POD-pGreedy algorithm, following the definitions of [37], as first proposed in [38] and analyzed in [39]. The algorithm combines data from different parameters to generate a single reduced basis set that covers the entire parameter space $\mathcal{P}$. The procedure is similar to Algorithm 1 but with few differences:

1) The reduced bases are shared between anchor points. POD is still used to construct the new basis but the collected snapshot set is projected onto the orthogonal complement of the existing basis.

2) In the online/training stage, only one ROM is used instead of a set of ROMs and there is no need to check for the nearest anchor points.

3) The anchor point is selected based on the single error estimate Δ in current iteration, rather than the individual estimates for each ROM.

Although it has a better convergence rate than POD-hGreedy, POD-pGreedy can easily fail for unsteady problems. Combining modes at different anchor points, especially ones whose solution exhibits different physics, can easily lead to instability and deteriorate the performance, as noted in [9]. Moreover, stabilizations that work for POD-hGreedy can fail in the POD-pGreedy approach. For example, in C-ROM, it is not clear how to construct the constraints for the combined basis. A naive approach is to apply POD to all the snapshots at anchor points. However, this approach is inefficient and can be limited by the computer storage requirements during the offline phase. Leray regularization with the projection filter (i.e., trivially truncated basis set for the advector) is also limited since the combined basis is no longer ordered in a Fourier-like, energy-decaying, sequence. To address this, one could apply POD to all the snapshots that have been collected but this approach is again limited by the storage and therefore not practical. An alternative is to consider DF filter, denoted as LDF-ROM here. Once the radius δ is specified, it will filter right amount of energy in each basis.

4.2 A time-averaged error indicator

In this section, we extend the time-averaged error indicator proposed in [9] to accommodate the energy equation and Leray regularization. The error indicator is based on the dual norm of the discrete time-averaged residual. Given the ROM solution sequence ${\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\}}_{n=J_0+1}^J$ and ${\{{\bar{T}}_{\text{c}}^n\}}_{n=J_0+1}^J$ and the parameters of interest $\underline{p}=(\nu ,{\theta }_g,\alpha )$, the discrete time-averaged residual for velocity and temperature are defined as: ⁴

$\langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\mathbf{\text{v}};\nu ,{\theta }_g\right):=\frac{1}{J-J_0}\sum _{n=J_0+1}^J{r}_u\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g\right),\forall \mathbf{\text{v}}\in {\mathbf{V}}_{\text{div}},(40)$

$\langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,S;\alpha \right):=\frac{1}{J-J_0}\sum _{n=J_0+1}^J{r}_T\left({\bar{T}}_{\text{c}}^n,S;\alpha \right),\forall S\in X_0^{\mathcal{N}},(41)$

where $r_u({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g)\in {\mathbf{V}}_{\text{div}}^{\prime }$ (dual space of ${\mathbf{V}}_{\text{div}}$) and $r_T({\bar{T}}_{\text{c}}^n,S;\alpha )\in {X_0^{\mathcal{N}}}^{\prime }$ (dual space of $X_0^{\mathcal{N}}$) are the residual associated with Eq. 32 at time tⁿ and defined as

$\begin{array}{ll}\hfill r_u\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g\right)& =\sum _{s=1}^3{\alpha }_s\left[\left(\mathbf{\text{v}},\mathbf{\text{g}}\left({\theta }_g\right){\bar{T}}_{\text{c}}^{n-s}\right)-c\left(\mathbf{\text{v}},{\overline{\mathbf{\text{u}}}}_{\text{c}}^{n-s},{\overline{\mathbf{\text{u}}}}_{\text{c}}^{n-s}\right)\right]\hfill \\ \hfill & -\sum _{s=0}^3\frac{{\beta }_s}{\mathrm{\Delta }t}\left(\mathbf{\text{v}},{\mathbf{\text{u}}}_{\text{c}}^{n-s}\right)-\nu a\left(\mathbf{\text{v}},{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right),\hfill \end{array}(42)$

$r_T\left({\bar{T}}_{\text{c}}^n,S;\alpha \right)=-\sum _{s=1}^3{\alpha }_{s}c\left(S,{\overline{\mathbf{\text{u}}}}_{\text{c}}^{n-s},{\bar{T}}_{\text{c}}^{n-s}\right)-\sum _{s=0}^3\frac{{\beta }_s}{\mathrm{\Delta }t}\left(S,T_{\text{c}}^{n-s}\right)-\alpha \left(S,{\bar{T}}_{\text{c}}^n\right).(43)$

Note for simplicity, we assume only BDF3/EXT3 is used for time discretization in Eqs. 42 and 43. Besides, the residual is defined over ${\mathcal{\{V\}}}_{\text{div}}:=\{\mathbf{\text{v}}|\mathbf{\text{v}}\in {\mathbf{\text{X}}}_0^{\mathcal{N}},\nabla \cdot \mathbf{\text{v}}=0\}$ and $X_0^{\mathcal{N}}$, rather than ${\mathbf{\text{X}}}_0^N\subset {\mathcal{H}}_0^1$ and $X_0^N$, because we measure our reduced-basis error relative to the FOM.

We define the time-averaged error indicator, $\mathrm{\Delta }:{\otimes }_{n=J_0}^J{\mathbf{\text{X}}}^N\times X^N\times \mathcal{P}\to {\mathbb{R}}_{+}$, as follows:

$\mathrm{\Delta }\left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,{\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J;\underline{p}\right):=\sqrt{\Vert \langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\nu ,{\theta }_g\right){\Vert }_{{\mathbf{V}}_{\text{div}}^{\prime }}^2+\Vert \langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\alpha \right){\Vert }_{{X_0^{\mathcal{N}}}^{\prime }}^2}.(44)$

The residuals Eqs. 42 and 43 can be further expressed in the matrix-vector form since the spaces ${\mathcal{\{V\}}}_{\text{div}}$ and $X_0^{\mathcal{N}}$ are finite dimensional,

$r_u\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g\right)={\underline{\mathbf{v}}}^T{\underline{\mathbf{r}}}_u^n={\underline{\mathbf{v}}}^T\widehat{\underline{\mathbf{\text{f}}}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n;{\theta }_g\right)-{\underline{\mathbf{v}}}^T\mathbf{\text{H}}\mathbf{\text{B}}{\widehat{\underline{u}}}_{\text{c}}^{n}n,\forall \underline{\mathbf{v}}\in {\mathbb{R}}^{2\mathcal{N}},(45)$

$r_T\left({\bar{T}}_{\text{c}}^n,S;\alpha \right)={\underline{S}}^T{\underline{r}}_T^n={\underline{S}}^T\widehat{\underline{Q}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n\right)-{\underline{S}}^T{H}_{\alpha }B{{\widehat{\underline{T}}}_{\text{c}}}^n,\forall \underline{S}\in {\mathbb{R}}^{\mathcal{N}}.(46)$

The matrix-vector version of the discrete time-averaged residual Eqs. 40 and 41 is then expressed as

$\langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\mathbf{\text{v}};\nu ,{\theta }_g\right)={\underline{\mathbf{v}}}^T{\underline{\mathbf{R}}}_u={\underline{\mathbf{v}}}^T\left(\frac{1}{J-J_0}\sum _{n=J_0+1}^J\widehat{\underline{\mathbf{\text{f}}}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n;{\theta }_g\right)-\mathbf{\text{H}}\mathbf{\text{B}}{{\widehat{\underline{u}}}_{\text{c}}}^n\right),(47)$

$\langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,S;\alpha \right)={\underline{S}}^T{\underline{R}}_T={\underline{S}}^T\left(\frac{1}{J-J_0}\sum _{n=J_0+1}^J\widehat{\underline{Q}}\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,{\bar{T}}_{\text{c}}^n\right)-H_{\alpha }B{{\widehat{\underline{T}}}_{\text{c}}}^n\right),(48)$

$\forall \underline{\mathbf{v}}\in {\mathbb{R}}^{2\mathcal{N}}$ and $\forall \underline{S}\in {\mathbb{R}}^{\mathcal{N}}$. The norm of the residual is closely related to the error and it is tempting to use $\Vert {\underline{R}}_u{\Vert }_2$ and $\Vert {\underline{R}}_T{\Vert }_2$ to estimate the error. However, this is not correct since $⟨R_u⟩({\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\}}_{n=J_0}^J,\cdot ;\nu ,{\theta }_g):{\mathbf{V}}_{\text{div}}\to \mathbb{R}$ and $⟨R_T⟩({\{{\bar{T}}_{\text{c}}^n\}}_{n=J_0}^J,\cdot ;\alpha ):X_0^{\mathcal{N}}\to \mathbb{R}$ are bounded linear functionals whose size is appropriately measured through the dual norm:

$\Vert \langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\nu ,{\theta }_g\right){\Vert }_{{\mathcal{\{V\}}}_{\text{div}}^{\prime }}={\sup }_{\mathbf{\text{v}}\in {\mathbf{V}}_{\text{div}}}\frac{\langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\mathbf{\text{v}};\nu ,{\theta }_g\right)}{\Vert \mathbf{\text{v}}{\Vert }_{{\mathbf{V}}_{\text{div}}}},(49)$

$\Vert \langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\alpha \right){\Vert }_{X_0^{\mathcal{N}}}={\sup }_{S\in {X_0^{\mathcal{N}}}^{\prime }}\frac{\langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,S;\alpha \right)}{\Vert S{\Vert }_{X_0^{\mathcal{N}}}}.(50)$

Thanks to the Riesz representation theorem, there exist a unique $⟨{\widehat{\mathbf{R}}}_u⟩\in {\mathbf{V}}_{\text{div}}$ and $⟨{\widehat{R}}_T⟩\in X_0^{\mathcal{N}}$ such that

${\left(\langle {\widehat{\mathbf{R}}}_u\rangle ,\mathbf{\text{v}}\right)}_{{\mathbf{V}}_{\text{div}}}=\langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\mathbf{\text{v}};\nu ,{\theta }_g\right),\forall \mathbf{\text{v}}\in {\mathbf{V}}_{\text{div}},(51)$

${\left(\langle {\widehat{R}}_T\rangle ,S\right)}_{X_0^{\mathcal{N}}}=\langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,S;\alpha \right),\forall S\in X_0^{\mathcal{N}}.(52)$

It thus follows that

$\Vert \langle R_u\rangle \left({\left\{{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\nu ,{\theta }_g\right){\Vert }_{{\mathbf{V}}_{\text{div}}^{\prime }}=\Vert \langle {\widehat{\mathbf{R}}}_u\rangle {\Vert }_{{\mathbf{V}}_{\text{div}}},(53)$

$\Vert \langle R_T\rangle \left({\left\{{\bar{T}}_{\text{c}}^n\right\}}_{n=J_0}^J,\cdot ;\alpha \right){\Vert }_{{X_0^{\mathcal{N}}}^{\prime }}=\Vert \langle {\widehat{R}}_T\rangle {\Vert }_{X_0^{\mathcal{N}}}.(54)$

Equations 51 and 52 allows one to compute the Riesz representers $⟨{\widehat{\mathbf{R}}}_u⟩$ and $⟨{\widehat{R}}_T⟩$ and Eqs. 53 and 54 allows one to evaluate the dual norm of the residual through Riesz representation without computing the supremum.

In practice, determination of the Riesz representers, $⟨{\widehat{\mathbf{R}}}_u⟩$ and $⟨{\widehat{R}}_T⟩$, is relatively straightforward because the coarse (i.e., ROM) and truth (FOM) representations live in finite-dimensional spaces, meaning that there is a direct linear-algebra problem to be solved for the Riesz representers. Expanding Eqs. 51 and 52, we have the corresponding linear algebra statement,

$\left[\begin{array}{cc}\hfill \mathbf{\text{A}}\hfill & \hfill -{\mathbf{\text{D}}}^T\hfill \\ \hfill -\mathbf{\text{D}}\hfill & \hfill 0\hfill \end{array}\right]\left(\begin{array}{c}\hfill \langle {\underline{\widehat{\mathbf{R}}}}_u\rangle \hfill \\ \hfill \underline{p}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\underline{\mathbf{R}}}_u\hfill \\ \hfill 0\hfill \end{array}\right),(55)$

$A\langle {\underline{\widehat{R}}}_T\rangle ={\underline{R}}_T.(56)$

Here, A corresponds to H introduced in Eq. 23 with β₀ = 0 and ν = 1. We remark that the essential difference between the velocity and temperature representers is that the former satisfies the divergence-free constraint by virtue of the 2 × 2 block system in Eq. 55. Evaluation of the error indicator Δ entails solving Eqs. 55 and 56, computing the corresponding ${\mathcal{H}}^1$ norms of the outputs, and ultimately using these results in Eq. 44.

While use of the direct approach requires access to the FOM machinery in order to generate an error indicator, we note that such access is readily available during the pMOR training/construction phase. The advantage of this approach is that the number of Stokes/Poisson solves scales as the number of points in the training space, which is typically less than N². The other is through the offline/online computational decomposition which takes the advantage of the affine decomposition and expands the residual. By expanding the residuals ⟨R_u⟩ and ⟨R_T⟩, 2(N + 1)² + 6(N + 1) linear functionals are derived, where 2(N + 1)² is due to the convection term in the Navier-Stokes and energy equations. Applying the Riesz representation theorem to each linear functional, we end up solving 2(N + 1)² + 6(N + 1) Riesz representers, where (N + 1)² + 4(N + 1) of them are solved through Stokes problems and (N + 1)² + 2(N + 1) of them are solved through Poisson problems. Note that the Riesz representers must be stored in order to accomplish the decomposition and each is a vector of size $\mathcal{N}$ since it lives in the FOM space. For example, if N = 60, one would have to store at least 7,200 vectors of size $\mathcal{N}$ which can be prohibitive, even for large multicore workstations. Even if there is no storage limitation, the offline cost is quite high when N is large as it scales quadratically. Once it is done, the online cost is $\mathcal{O}(N^{2}J+N^4)$, where $\mathcal{O}(N^{2}J)$ is to solve Eq. 32 and $\mathcal{O}(N^4)$ is required to compute the error estimate. Further details of the decomposition are provided in [9].

4.2.1 Time-averaged error indicator with Leray regularization

The integration of the time-averaged error indicator Δ with Leray regularization is rather straightforward. Recall the difference between G-ROM Eq. 32 and LDF-ROM Eq. 35 is simply the advecting field being filtered. Hence, the residuals $r_u({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g)$ and $r_T({\bar{T}}_{\text{c}}^n,S;\alpha )$ for all n = J₀ + 1, … , J are simply modified with the filtered advecting field,

$\begin{array}{ll}\hfill r_u\left({\overline{\mathbf{\text{u}}}}_{\text{c}}^n,\mathbf{\text{v}};\nu ,{\theta }_g\right)& =\sum _{s=1}^3{\alpha }_s\left[\left(\mathbf{\text{v}},\mathbf{\text{g}}\left({\theta }_g\right){\bar{T}}_{\text{c}}^{n-s}\right)-c\left(\mathbf{\text{v}},{\overline{\mathbf{\text{u}}}}_{\text{c},\text{filtered}}^{n-s},{\overline{\mathbf{\text{u}}}}_{\text{c}}^{n-s}\right)\right]\hfill \\ \hfill & -\sum _{s=0}^3\frac{{\beta }_s}{\mathrm{\Delta }t}\left(\mathbf{\text{v}},{\mathbf{\text{u}}}_{\text{c}}^{n-s}\right)-\nu a\left(\mathbf{\text{v}},{\overline{\mathbf{\text{u}}}}_{\text{c}}^n\right),\hfill \end{array}(57)$

$r_T\left({\bar{T}}_{\text{c}}^n,S;\alpha \right)=-\sum _{s=1}^3{\alpha }_{s}c\left(S,{\overline{\mathbf{\text{u}}}}_{\text{c},\text{filtered}}^{n-s},{\bar{T}}_{\text{c}}^{n-s}\right)-\sum _{s=0}^3\frac{{\beta }_s}{\mathrm{\Delta }t}\left(S,T_{\text{c}}^{n-s}\right)-\alpha a\left(S,{\bar{T}}_{\text{c}}^n\right),(58)$

for all $\mathbf{\text{v}}\in {\mathbf{V}}_{\text{div}}$ and $S\in X_0^{\mathcal{N}}$. The corresponding time-averaged error indicator Δ is then defined based on the modified residuals Eqs. 57 and 58.

4.3 Parametric model order reduction results: θ_g variation

In this section, we consider the parametric problem parameterized with θ_g at Pr = 7.2. The problem has three characteristics: bifurcation, spatio-temporal chaos over a certain range of θ_g, and a solution manifold that is a blend of steady and unsteady solutions. To identify the major pMOR challenges for this case, three values of Ra are considered:

1) Ra = 1 × 10⁴ where the FOM is steady except at θ_g = 30°.

2) Ra = 8 × 10⁴ where the FOM is unsteady for θ_g ∈ [0°, 70°] and steady for θ_g ∈ [80°, 180°].

3) Ra = 3 × 10⁵ where the FOM is unsteady for θ_g ∈ [0°, 120°] and steady for θ_g ∈ [130°, 180°].

In each case, the ROM is constructed through Algorithm 1. In order to assess performance, we generate FOM data for θ_g = 0°, 10°, …, 180° (n_train = 19 datapoints). The FOM solution is obtained by solving Eqs. 2 and 3. For parameters where the problem is steady, the solution and the Nu are collected after the solution difference between ten time steps is less than 10^–6. For unsteady problems, the mean flow, mean Nu, Std(Nu), mean TKE and mean temperature fluctuation are averaged over 500 CTUs after the solution has reached a statistically steady state.

Although not shown here, we remark that at Ra = 1 × 10³, the problem is steady with a bifurcation at θ_g = 20°. In this case, either h- or p-Greedy with residual dual-norm base error indicator accurately estimates the solution and QOIs over the parameter space θ_g ∈ [0°, 180°].

4.3.1 Ra = 1 × 10⁴

To examine the feasibility of the pMOR (Algorithm 1) in the unsteady case, we begin with Ra = 1 × 10⁴, in which only one unsteady solution is introduced at θ_g = 30°.

Figure 3E shows the steady (or mean) velocity magnitude for 19 uniformly-space training points θ_g = 0°, 10°, …, 180°. The corresponding temperature distributions are in Figure 1B. At θ_g = 180°, we observe no flow and the temperature is simply the conduction solution. As θ_g decreases, we observe slot convection and then about θ_g = 40° there is a bifurcation to the wavy flow and rolls in the velocity. Moreover, we observe spatial-temporal chaos at θ_g = 30°. Figures 3A–D show the results of the application of Algorithm 1 for the construction of the G-ROM, C-ROM and LDF-ROM for the pMOR. The algorithm starts with ${\theta }_{g,1}^{*}=0°$ and is performed with L = 8 iterations.

FIGURE 3

FIGURE 3. POD-hGreedy performance comparison between G-ROM, C-ROM and LDF ROM for the parametric problem parameterized with θ_g at Ra = 1 × 10⁴ and Pr = 7.2. (A) Behavior of the error indicator Δ defined in Eq. 44 for eight iterations with G-ROM. (B–D) Behavior of the relative ${\mathcal{H}}^1$ error in predicted solution, mean Nu and mean TKE with θ_g based on eight anchor points. (N = 1 for all ${\theta }_g\in {\mathcal{P}}_{\text{anchor}}$ except N = 70 at ${\theta }_{g,3}^{\ast }=$ 30°.). (E) Mean (or steady) velocity magnitude for 19 uniformly-spaced θ_g ∈ [0°, 180°]. (The corresponding temperature solutions are in Figure 1B).

Figure 3A demonstrates the selection process of anchor points (denoted by red circles) for the G-ROM case. We briefly walk through the process: At the first iteration, the error estimate Δ₁(θ_g) for ${\theta }_g\in {\mathcal{P}}_{\text{train}}$ is computed (blue dashed line). With the largest error estimate, θ_g = 90° is then identified as the second anchor point. The third anchor point is then selected from ${\theta }_g\in {\mathcal{P}}_{\text{train}}$ which maximizes the error estimate Δ_1,2(θ_g): = min{Δ₁(θ_g), Δ₂(θ_g)} over ${\mathcal{P}}_{\text{train}}$. (We reiterate that minimizing over the individual error estimates is a property of the h-Greedy process—there is not a single unifying error estimate as is the case for p-Greedy.) The process continues until the error estimate reaches the desired tolerance or the number of offline solves reaches its maximum. The black solid line denotes the minimum of all error estimates computed up to current iteration. In this case, it represents min{Δ₁(θ_g), …, Δ₈(θ_g)}. Note that the error estimate at θ_g = 180° in each model Δ_ℓ(θ_g = 180°) is small due to the choice of lift function.

From Figure 3A, we observe the error estimate is small at anchor points where the problem is steady. On the other hand, although the error estimate Δ₃(θ_g) (greed line) is small at ${\theta }_{g,3}^{\ast }=30°$ compared to other points in ${\mathcal{P}}_{\text{train}}$, it is not as small as the estimate at other anchor points. Because of the unsteadiness, it can’t reach the same magnitude as in the steady cases.

Following this procedure for the other cases, we present models results for the G-ROM, C-ROM and LDF-ROM cases, denoted respectively blue, orange, green solid lines in Figures 3B–D. The behavior of the relative ${\mathcal{H}}^1$ error in the predicted solution is shown in Figure 3B. The corresponding Galerkin projection is denoted as the black dashed line. We found that G-ROM and LDF-ROM have a similar performance: the error in the solution is nearly identical to the Galerkin projection for cases where the problem is steady, including those that are not in the ${\mathcal{P}}_{\text{anchor}}$. The maximum error is at θ_g = 30° where the problem is unsteady. Both methods have around 15% error in the mean flow. Note that at θ_g = 30°, the number of modes N is carefully selected since the ROM diverges after certain N due to the spatio-temporal chaos (15% error with N = 70, 17% error with N = 80 and 23% error with N = 90). On the other hand, the mean solution prediction made by C-ROM has 73% in maximum error and in order for C-ROM to reach same accuracy as in G-ROM and LDF-ROM, two more iterations are required. Already, with this modest Ra = 1 × 10⁴, we find C-ROM is less efficient than G-ROM and LDF-ROM. The pMOR behavior for mean Nu and mean TKE are shown along with the FOM results in Figures 3C,D. Again, G-ROM and LDF-ROM are able to make accurate predictions while C-ROM has maximum 22% error in mean Nu and in particular is unable to capture the peak in mean TKE at θ_g = 30°.

Before closing this section, we highlight some observations with respect to the solution manifolds.

1) The solutions at θ_g ∈ [0°, 30°] are Rayleigh-Bénard with differing numbers of rolls, analogous to orthogonal sine and cosine functions at different wave numbers. Therefore, there is little hope in reproducing the solution except at selected anchor points. QOI’s such as mean Nu, however, are less sensitive to precise mean flow fields and are therefore more tractable.

2) At θ_g = 170° and θ_g = 180° the thermal metrics are not too different despite the $\mathcal{O}(1)$ difference in velocity solutions.

The first issue is resolved by the error indicator picking θ_g ∈ [0°, 30°] as anchor points. The second issue can be a source of error as Ra increases. With θ_g = 160° as an anchor point and the solution at θ_g = 180° as the lift function, the error at θ_g = 170° is 9% for Ra = 1 × 10⁴, 16% for Ra = 8 × 10⁴, and 19% for Ra = 3 × 10⁵.

4.3.2 Ra = 8 × 10⁴

Figure 4 shows pMOR results analogous to Figure 3 for the case Ra = 8 × 10⁴. Here, we consider only the G-ROM and LDF-ROM. The algorithm starts with ${\theta }_{g,1}^{\ast }=0°$ and terminates at L = 6 iterations. Figures 4E,F show the FOM mean (or steady) temperature and velocity solution at the training points. (In an actual pMOR, these FOMs would of course not be computed a priori.) With this Ra, the bifurcation occurs at θ_g = 40°. Moreover, we observe spatio-temporal chaos for θ_g ∈ [0°, 30°] with the lower values being more chaotic.

FIGURE 4

FIGURE 4. POD-hGreedy performance comparison between G-ROM and LDF ROM for the parametric problem parameterized with θ_g at Ra = 8 × 10⁴ and Pr = 7.2. (A) Behavior of the error indicator Δ defined in Eq. 44 for six iterations with LDF-ROM. (B–D) Behavior of the relative ${\mathcal{H}}^1$ error in predicted solution, mean Nu, and mean TKE with θ_g based on six anchor points. (N = 100 at θ_g = 0°, 20°, N = 80 at θ_g = 10°, 30°, 40°, N = 1 at θ_g = 110°, 140°, 160°). (E,F) Mean (or steady) temperature and velocity solution for 19 uniformly spaced θ_g ∈ [0°, 180°].

The anchor-point selection process with LDF-ROM is demonstrated in Figure 4A. Starting with ${\theta }_{g,1}^{\ast }=$ 0° the peak error in first iteration is at 110°, which is chosen to be ${\theta }_{g,2}^{\ast }$, and so on. Again, we find the error indicator is small at anchor points where the problem is steady and that it is larger where it is unsteady (θ_g ∈ [0°, 60°]). Nonetheless, the error indicator is still able to identify where solution changes rapidly and select most of the anchor points in the region [0°, 40°].

The behaviors of the relative ${\mathcal{H}}^1$ error in the predicted solution with θ_g using G-ROM and LDF-ROM are shown in Figure 4B. For θ_g ∈ [80°, 180°], where the solution is steady, we find the estimation is almost identical to the Galerkin projection in both models. On the other hand, for θ_g ∈ [0°, 70°], where the solution is unsteady, we find the error at anchor points θ_g = 0°, θ_g = 10° is large due to the spatio-temporal chaos. The maximum error is around 19% at θ_g = 10° with LDF-ROM while 69% at θ_g = 40° with G-ROM. Although the maximum error in G-ROM can be reduced by further iterations of the algorithm, the error will eventually be dominated by the high reproduction error arising from spatio-temporal chaos at θ_g = 0° and 10°.

The behavior for mean Nu and mean TKE are shown in Figures 4C,D for G-ROM and LDF-ROM. Despite large errors in the mean flow prediction at θ_g = 0°, 10°, the LDF-ROM is able to predict mean Nu with a maximum error around 5% whereas G-ROM has maximum error around 28%. In addition, LDF-ROM is able to more accurately predict mean TKE than G-ROM.

4.3.3 Ra = 3 × 10⁵

For (Pr,Ra) = (7.2, 3 × 10⁵) the flow is quite chaotic (similar to what is found for Pr = 0.71 at lower Ra). Figures 5E,F show the mean (or steady) temperature and velocity solution at θ_g = 0°, 10°, …, 180° (19 datapoints). This time, the bifurcation occurs at θ_g = 60°. We use this elevated Rayleigh-number case to explore the behavior of the h-Greedy pMOR convergence by considering application of the algorithm to two different training sets, ${\mathcal{P}}_1=$ [60°, 70°, … , 180°] and ${\mathcal{P}}_2=$ [0°, 10°, … , 180°]. The set ${\mathcal{P}}_1$ excludes the spatio-temporal chaotic regime while ${\mathcal{P}}_2$ spans the full range of flow phenomena.

FIGURE 5

FIGURE 5. POD-hGreedy performance comparison between G-ROM and LDF ROM for the parametric problem parameterized with θ_g ∈ [60°, 180°] at Ra = 3 × 10⁵ and Pr = 7.2. (A) Behavior of the error indicator Δ defined in Eq. 44 for five iterations with LDF-ROM. (B–D) Behavior of the relative ${\mathcal{H}}^1$ error in predicted solution, mean Nu and mean TKE with θ_g based on five anchor points. (N = 80 at θ_g = 60°, 80°, 110°, N = 50 at θ_g = 120° and N = 1 at θ_g = 130°, 140°, 160°). (E,F) Mean (or steady) temperature and velocity solution for 19 uniformly spaced θ_g ∈ [0°, 180°].

The anchor point selection process for ${\mathcal{P}}_1$ with LDF-ROM is demonstrated in Figure 5A, starting with ${\theta }_{g,1}^{\ast }=60$° and proceeding for L = 5 iterations. Again, we observe that the error estimate at the anchor points, ${\theta }_{g,1}^{\ast }=60$°, ${\theta }_{g,3}^{\ast }=80$° and ${\theta }_{g,5}^{\ast }=120$° are larger than other anchor points because of unsteadiness. The behavior of the relative ${\mathcal{H}}^1$ error in predicted solution is shown in Figure 5B. For θ_g ∈ [130°, 180°], where the solution is steady, the ROM estimates at the anchor points are almost identical to the Galerkin projection. For θ_g ∈ [60°, 120°], where the solution is unsteady, the errors at the anchor points (θ_g = 60°, 80°, 110°, 120°) are less than 10%. However, because of the irregular solution manifold, there is a 20% maximum error at θ_g = 170°, despite the ROM being based on the nearby θ_g = 160° anchor point.

The behavior of the mean Nu and mean TKE are shown in Figures 5C,D. The maximum error in the predicted mean Nu is around 5% with LDF-ROM and 8% with G-ROM. The mean TKE estimation is also reasonable but is underestimated at θ_g = 70°.

Next we examine the same problem configuration but with the full parameter space ${\mathcal{P}}_2$. The problem now includes spatio-temporal chaos for θ_g ∈ [0°, 40°] with the lower values being more chaotic. Figure 6 show the results of the application of Algorithm 1 for the construction of the LDF-ROM for the parametric problem with ${\mathcal{P}}_2$. The algorithm starts with ${\theta }_{g,1}^{*}=0°$ and is performed with L = 6 iterations. The anchor point selection process is demonstrated in Figure 6A. We observe the same issue as in the previous cases, where unsteadiness leads to larger error estimates than with the steady regimes.

FIGURE 6

FIGURE 6. POD-hGreedy performance comparison between G-ROM and LDF ROM for the parametric problem parameterized with θ_g ∈ [60°, 180°] at Ra = 3 × 10⁵ and Pr = 7.2. (A) Behavior of the error indicator Δ defined in Eq. 44 for five iterations with LDF-ROM. (B–D) Behavior of the relative ${\mathcal{H}}^1$ error in predicted solution, mean Nu and mean TKE with θ_g based on five anchor points. (N = 80 at θ_g = 10°, 20°, 40°, 60°, 90°, 100°, N = 70 at θ_g = 0°, N = 50 at θ_g = 120° and N = 1 at θ_g = 160°.)

The relative ${\mathcal{H}}^1$ error in predicted solution is shown in Figure 6B. Again we find the estimation is almost identical to the Galerkin projection for θ_g ∈ [130°, 180°] where the solution is steady. For θ_g ∈ [0°, 120°] where the solution is unsteady, the errors at anchor points θ_g = 20°, 40°, 90°, 120° are less than 10% but 35% at θ_g = 0°, which corresponds to “simple” Rayleigh-Bénard convection. Note that 35% is the error after carefully chosen N and spatial radius δ in Leray filtering. The corresponding mean Nu and mean TKE behavior are shown in Figures 6C,D. The maximum error in the predicted mean Nu is around 12%. For mean TKE, the estimation for θ_g ∈ [60°, 180°] is acceptable, while it is overestimated for θ_g ∈ [0°, 50°].

We are also aware that in some applications, the Std(Nu) could be considered as QOI. However, comparing to the mean Nu and mean TKE, we find Std(Nu) is in general a more challenging QOI. Figure 7 shows the predicted Std(Nu) in the three Ra cases. We observe accurate prediction in Ra = 1 × 10⁴ case. However, unlike the mean TKE, the Std(Nu) soon becomes intractable with Ra = 8 × 10⁴ even with Leray regularization and is even worse in Ra = 3 × 10⁵.

FIGURE 7

FIGURE 7. Behavior of the predicted Std(Nu) with θ_g. (A): G-ROM, C-ROM and LDF-ROM estimation at Ra = 1 × 10⁴ with eight anchor points. (B): G-ROM and LDF-ROM estimation at Ra = 8 × 10⁴ with six anchor points. (C): G-ROM and LDF-ROM estimation at Ra = 3 × 10⁵ with six anchor points.

4.4 Parametric model order reduction results: Ra variation

In this section, we consider the slot problem at θ_g = 90° and Pr = 7.2 with the parametric space defined by $\mathrm{R}\mathrm{a}\in \mathcal{P}=[2\times {10}^5,7\times {10}^5]$. Unlike the problem with θ_g variation, all solutions are unsteady and there is no parametric bifurcation. In order to assess performance, we generate FOM data, including mean flow, mean Nu, Std(Nu), mean TKE and mean temperature fluctuation, for Ra = 2 × 10⁵, 2.5 × 10⁵, … , 7 × 10⁵ (n_train = 11 datapoints). The quantities are averaged over 500 CTUs once the solution reaches the statistically steady state.

Figure 8 shows the results of the application of Algorithm 1 for the pMOR using G-ROM, C-ROM and LDF-ROM. The solid line denotes the performance of the reduced model which minimizes the error indicator, and thus is selected by the Greedy procedure (cf. Algorithm 1, n_cand = 2). Anchor points are denoted as red circle while FOM data is denoted as black solid line. The algorithm starts with ${\mathrm{R}\mathrm{a}}_1^{\ast }=2\times {10}^5$ and is performed with L = 5 iterations. The number of POD basis N with anchor points are listed in the figure caption.

FIGURE 8

FIGURE 8. POD-hGreedy performance comparison between G-ROM, C-ROM and LDF ROM for the parametric problem parameterized with Ra at θ_g = 90° and Pr = 7.2. (A–D) Behavior of the relative ${\mathcal{H}}^1$ error in predicted mean flow, predicted mean Nu, Std(Nu) and mean temperature fluctuation based on five anchor points. (N = 80, ${\mathrm{R}\mathrm{a}}_1^{\ast }=2\times {10}^5$, ${\mathrm{R}\mathrm{a}}_2^{\ast }=7\times {10}^5$, ${\mathrm{R}\mathrm{a}}_3^{\ast }=2.5\times {10}^5$, ${\mathrm{R}\mathrm{a}}_4^{\ast }=3.5\times {10}^5$, ${\mathrm{R}\mathrm{a}}_5^{\ast }=5\times {10}^5$ for C-ROM and LDF-ROM and ${\mathrm{R}\mathrm{a}}_5^{\ast }=4.5\times {10}^5$ for G-ROM).

Figure 8A shows the behavior of the relative ${\mathcal{H}}^1$ error in mean flow prediction with Ra. First, we observe the errors at the anchor points are less than 10% with C-ROM and LDF-ROM while G-ROM has 10% error at ${\mathrm{R}\mathrm{a}}_2^{\ast }=7\times {10}^5$. The maximum error is roughly 11% in G-ROM, 10% in LDF-ROM and 8% in C-ROM. Comparing with the Galerkin projection error (denoted by the black dashed line), the pMOR accuracy is seen to be quite satisfactory throughout $\mathcal{P}$.

Figure 8B shows the behavior of the predicted mean Nu with Ra. At the anchor points, we observe good agreement between ROMs and FOM and that stabilization does improve its accuracy. The maximum relative error is roughly 8% in G-ROM, 6.5% in C-ROM and 5% in LDF-ROM. Figures 8C,D show the behavior of the predicted Std(Nu) and mean temperature fluctuation. In both QOIs, we find C-ROM outperforms the other two models. At Ra = 7 × 10⁵, LDF-ROM is only slightly better than G-ROM.

This parametric space is in general more tractable than those involving variation in θ_g. This outcome might be anticipated by observing the mean temperature fields shown in Figure 1C, which suggests that the solution manifold with respect to Ra is quite smooth. This is also reflected in the QOIs, for example, the mean Nu, Std(Nu) and mean temperature fluctuation behave almost linearly as Ra increases.

5 Discussion

In this section, we investigate some of the flow behaviors exhibited by the FOM to better understand how they influence the relative performance of the ROMs. We note that we cannot, in general, expect a ROM to be able to predict FOM behavior if the flow itself is not predictable. Thus, variability in the FOM provides an anticipated lower bound on ROM performance for the reproduction problem.

5.1 Spatio-temporal chaos

As pointed out in the introduction, we classify a flow to be spatio-temporal chaotic by examining its consistency in mean flow over various time windows. (We use this simple metric here for convenience—we have also examined the flow fields and the time traces of multiple QOIs.) Here, we explore how lack of consistency influences four QOIs, mean Nu, Std(Nu), mean temperature fluctuation, and mean TKE, at three successive time windows, W1, W2 and W3. These quantities are used to indicate the variability in the FOM. As with the preceding cases, we consider averaging times of 500 CTUs for each of the three windows. The starting time for W1 differs with given parameters as some cases take a longer time to reach a statistically steady state. For example, for Pr = 0.71 at Ra = 1.8 × 10⁴ and θ_g = 90°, the flow is chaotic until 6,000 CTUs and then becomes periodic.

Figures 9A–D show the behavior of the four QOIs with Ra at three Pr for θ_g = 90°. Pr = 0.07 is denoted as green line, Pr = 0.71 is denoted as orange line, while Pr = 7.2 is denoted as blue line. Window W1 is denoted by a solid line, W2 by a dashed line, and W3 by a dotted line.

FIGURE 9

FIGURE 9. Parametric variability in the FOM: green–Pr = 0.07, orange–Pr = 0.71, and blue–Pr = 7.2. Each plot reveals absence/presence of chaotic effects by presenting statistics taken over three time windows, W1, W2, and W3. (A–D): Ra-dependence of mean Nu, Std(Nu), mean temperature fluctuation, and mean TKE computed over time windows W1, W2 and W3. (E–H): θ_g-dependence at fixed Ra.

From Figures 9A–D we can see the following:

1) For Pr = 0.07 the QOIs are fairly consistent except for Ra ∈ [3 × 10⁴, 5 × 10⁴].

2) For Pr = 0.71, we find large variability in Std(Nu), mean temperature fluctuation and mean TKE for Ra > 2 × 10⁴.

3) Pr = 7.2 exhibits the least variability.

For Ra where we find that the QOI variability is high, we have also examined the mean flow at multiple time windows and found that those are also inconsistent. In all cases, the mean Nu is quite repeatable.

Figures 9E–H show the behavior of the same QOIs as a function of θ_g. For Pr = 7.2, we consider three different values of Ra. We find for most of the θ_g, the variability is small except for small θ_g where we also report spatio-temporal chaotic flow. For Pr = 0.07, we find the QOIs has small variability with Ra = 1 × 10⁴. However, for Pr = 0.71, we find large variability, especially in the mean TKE. Not only it has spatio-temporal chaotic flow (for example θ_g = 100°), but also the solution manifold is not smooth. By varying θ_g = 80° to θ_g = 100°, the solution changes from steady to periodic then spatio-temporal chaotic.

From Figures 9E–H we observe the following:

1) For Pr = 0.07, Std(Nu) exhibits up to 50% variability (e.g., at θ_g = 70°) while ⟨T_fluc⟩ and ⟨TKE⟩ have orders-of-magnitude relative variability at θ_g = 0°.

2) For Pr = 0.71, ⟨T_fluc⟩ and ⟨TKE⟩ exhibit significant variability for θ_g ∈ [60°, 110°].

3) For Pr = 7.2, the most notable variation is at θ_g = 30° for Std(Nu), ⟨T_fluc⟩, and ⟨TKE⟩ at Ra = 10⁴. Remarkably, the higher Rayleigh number cases do not exhibit as much variance.

As in Figures 9A–D, the mean Nu is seen to be a repeatable QOI. It is worth noting the real challenge and sensitivity of this class of problems is illustrated in Figure 9H. Here, we observe for the (Pr, Ra) = (0.71, 10⁴) case that the flow alternates from steady to unsteady at multiple points along the one-dimensional θ_g parameter space as indicated by several distinct zeros in the TKE.In Figure 10, we further explore the influence of spatio-temporal chaos by examining the mean-flow distributions and ROM performance for (Pr, Ra) = (7.2, 3.5 × 10⁴) at θ_g = 0° and θ_g = 20°. Figures 10A–C show the behavior of the ${\mathcal{H}}^1$ error and mean Nu predictions as a function of N, along with the mean-velocity magnitude distributions over seven time windows for θ_g = 0°. Figures 10D–F show the same quantities but with θ_g = 20° and only three time windows. In Figure 10C, we observe that for θ_g = 0° the number of rolls in the mean velocity field changes across different time windows. Similar changes are observed, to a lesser extent, at θ_g = 20°. Hence, both solutions are categorized as spatio-temporal chaotic flow, but the θ_g = 0° case is more significant. Comparing the mean flow error and mean Nu, we observe that the ROM convergence for the reproduction problem is slower (or nonexistent) at θ_g = 0°, while the convergence behavior is more favorable at θ_g = 20°.

FIGURE 10

FIGURE 10. ROM performance comparison between problem having significant (θ_g = 0°) and less significant (θ_g = 20°) spatio-temporal chaotic flow at Pr = 7.2 and Ra = 3 × 10⁵. (A–C) Behavior of the relative ${\mathcal{H}}^1$ error and mean Nu as a functions of N and magnitude of the mean velocity at multiple time windows W1, … , W7 at θ_g = 0°. (D–F) performance for the same quantities at θ_g = 20°.

We have also computed the relative error between FOM mean flows across seven time windows for the two values of θ_g. The maximum relative ${\mathcal{H}}^1$ error is 34% for θ_g = 0° and 10% for θ_g = 20°. These FOM discrepancies can be considered as a bound on the predictive capabilities of the ROM. Indeed, the values of 34 and 10% are consistent with the lower bounds realized in Figures 10A,D.

In Figure 11, we examine the influence of Prandtl number by comparing results for Pr = 0.71 and 7.2 at (Ra, θ_g)= (3 × 10⁵, 90°). Figures 11A–D show the convergence behavior for the ${\mathcal{H}}^1$ error and mean Nu as well as mean temperature and x-velocity fields at three time windows for Pr = 0.71, while Figures 11E–H show the same quantities for Pr = 7.2. From the mean fields, we observe that the number of rolls and its position changes with time window at Pr = 0.71, while minimal variance is observed at Pr = 7.2. Hence the solution at Pr = 0.71 is considered to be spatio-temporal chaotic while only chaotic at Pr = 7.2. Comparing the behavior of the relative ${\mathcal{H}}^1$ error in the mean flow and mean Nu, we observe convergence issues in the ROM at Pr = 0.71, while the same metrics converge at Pr = 7.2. We further compute the relative variance between two FOM mean flows across seven time windows for the two considered θ_g. The maximum relative ${\mathcal{H}}^1$ error is 14% for Pr = 0.71 and only 1% for Pr = 7.2. This again explains why the approximation errors are different between the two cases. Considering these variance levels as a lower bound for the ROM, we can view 20% error in C-ROM as acceptable. On the other hand, the approximation error for Pr = 7.2 is able to reach below 10%.

FIGURE 11

FIGURE 11. ROM performance comparison between Pr = 0.71 and Pr = 7.2 at Ra = 3 × 10⁵ and θ_g = 90°. (A–D) Behavior of the relative ${\mathcal{H}}^1$ error and mean Nu as a functions of N, temperature and mean x-velocity at three time windows, W1, W2 and W3 at Pr = 0.71. (E–H): Performance for the same quantities at Pr = 7.2.

In summary, the results of this section show that convergence issues and variations in the QOIs in the ROM can have high a correlation with the flow being spatio-temporal chaotic. From Table 1 and Figure 9, we see that Pr = 0.71 has a more complicated solution manifold found with the other two Prandtl numbers and also exhibits spatio-temporal chaos at a relatively small Rayleigh number, Ra = 10⁴.

5.2 Multiple steady-state solutions

In this section, we report the existence of multiple steady-state solutions observed for the case (θ_g, Ra)=(0°, 10⁴). The variations are characterized by different numbers of recirculation rolls, which are induced by using different initial conditions. Figure 12 shows steady temperature solutions generated by starting with steady solutions, χ_n from other values of θ_g = n°, save for the χ₉₀ case, which corresponds to a single snapshot of the unsteady flow/temperature field at θ_g = 90°. Multiple steady states are also observed for this Prandtl number at θ_g = 10° and 20° and have been reported by other authors as well [12, 14, 40].

FIGURE 12

FIGURE 12. Different steady temperature solutions at θ_g = 0°, Ra = 1 × 10⁴ and Pr = 0.71. The solutions were computed from different initial conditions, χ_n, corresponding to FOM solutions at θ_g = n°. For n = 2, 10, 20, and 180 the solution is steady, whereas χ₉₀ is simply a snapshot from the θ_g = 90° case.

For solution reproduction, the multiplicity of the solutions is not an issue as long as the ROM uses the same initial condition as the FOM. However, for parametric problem, these multiple states could easily lead to an incorrect (or at least, unexpected or unverifiable) conclusion. For example, if the ROM anchor at θ_g = 10° is used to approximate the solution at θ_g = 0°. With the initial condition at θ_g = 10°, the approximate solution will be the third temperature solution shown in Figure 12. However, if one collects the FOM data at θ_g with zero initial condition, one will consider the first temperature solution as the truth solution. As we could consider those roll solutions as sine and cosine functions, the first and third solution are nearly orthogonal; their difference is $\mathcal{O}(1)$ and one would thus conclude that the pMOR had failed when it in fact had generated a valid solution.

5.3 Discussion summary

We have noted in Table 1 the broad range of flow regimes encountered for the tilted slot problem and in this section have illuminated a correlation between the flow states and predictive power of the MOR/pMOR framework. The cases with spatio-temporal chaos are generally the most challenging for model-order reduction and the pMOR errors are found to be (approximately) bounded from below by the variance observed in successive FOM simulations performed at the same parametric point. The development of the pMOR thus needs to be performed with care.

Two parameterizations were considered: 1) θ_g-variation, where a bifurcation exists and solution space is a blend of unsteady and steady solutions, and 2) Ra-variation, where no bifurcation exists and one finds only unsteady solutions. In the θ_g-variation problem, accurate prediction in mean flow and other QOIs by the pMOR was demonstrated in the Ra = 1 × 10⁴ and Ra = 8 × 10⁴ cases. In high Ra cases, acceptable prediction of Nu is achieved with LDF-ROM but a small mean-flow error is not realizable because of spatio-temporal chaos.

The results also indicate that the LDF-ROM is a better candidate for parametric problems with bifurcation than C-ROM. This observation is new, yet consistent with the results of [9], where the authors show that C-ROM is effective for parametric problems that do not have a bifurcation. For the parametric problem parameterized with Ra, without spatio-temporal chaos, we find that pMOR with any of the three methods, G-ROM, C-ROM, or LDF-ROM, is able to predict the mean flow quite well. In this case, C-ROM is the most accurate in mean flow prediction and other targeted QOIs. This result is not surprising given that the solution manifold does not have a bifurcation. Lastly, we remark that Std(Nu) is generally the most challenging QOI of those explored here.

From the results, we are able to make an important observation. For parametric problems where pMOR is successful (e.g., errors $<$ 10%), the solution is either only chaotic (e.g., Ra variation with θ_g = 90°) or the solution does not have significant spatio-temporal chaos (e.g., θ_g variation with Ra = 1 × 10⁴, 8 × 10⁴). Once the spatio-temporal chaos becomes significant, the predictive power of pMOR deteriorates and the maximum errors are dominated by variance in the truth solution.

Although not shown here, we have also applied POD-pGreedy to this problem. In the parametric problem parameterized with θ_g, it works only in the steady case Ra = 1 × 10³. Once the unsteady solution emerges, for example at Ra = 1 × 10⁴, combining modes associated with different values of θ_g leads to an unstable ROM even with the Leray regularization. Although no rigorous proof is given, we hypothesize that the issue is due to the bifurcation in solution behavior. This point was also suggested in [9], which empirically showed that combining modes associated with qualitatively different behaviors might lead to poor prediction. By contrast, when the current problem is parameterized with Ra we find that POD-pGreedy is more efficient than the h-refinement approach.

6 Conclusion

In this paper an error-indicated pMOR is applied to a 2D unsteady natural convection in a tilted high-aspect ratio slot. We first considered the solution reproduction problem (non-predictive case) to demonstrate the convergence of the ROMs and the effectiveness of the stabilization methods. We next addressed the parametric problem (predictive case) to validate the error indicator and, more broadly, the stabilized POD-hGreedy procedures. Principal contributions include, 1) extension of the error indicator proposed in [9] to buoyancy-driven flows; 2) demonstration that Leray-regularized Galerkin ROMs provide a robust solution approach for this class of flows; 3) identification of spatio-temporal chaos as a source of irreproducibility in both the FOM and the ROM and that the variance in the FOM provides a lower bound on the pMOR error in these cases; 4) observation that accurate prediction $(<10\%)$ with pMOR is achievable if the solution in the parametric space is either only chaotic or the spatio-temporal chaos is not significant, regardless of whether a bifurcation exists or not. Once the spatio-temporal chaos becomes significant, the performance of the pMOR deteriorates and the maximum errors of the mean flow and QOIs are dominated by the flow chaos.

We also highlight a number of challenges that are particularly relevant for buoyancy-driven flows and which should be taken into consideration in the design of pMOR strategies for 3D buoyancy driven turbulent flow. First, one needs to be aware of potential convergence issues for the mean flow and other QOI predictions when the FOM exhibits large-scale spatio-temporal chaos. Second, it is difficult to combine modes associated with different flow regimes, especially for the pGreedy case. Third, even relatively simple (e.g., steady) flows can exhibit multiple states at a given parameter. And fourth, there are large offline costs both in terms of computational time and required storage for error indicator and O(N⁴) costs for online-only error indicators⁵.

We outline potential next steps in pMOR development for this class of problems.

1) Extension to higher dimensional parameter space. In this work, we considered only one-dimensional parameter space since the pMOR behavior needed to be carefully diagnosed; however higher dimensional parameter spaces are more interesting for engineering applications.

2) hp-Greedy with a bifurcation detection technique. Although we find LDF-ROM is more efficient than C-ROM for parametric problems that have a bifurcation, the h-refinement strategy considered in this paper might require an infeasible number of offline simulations as the dimension of the parameter space increases. To tackle complex parametrizations, more advanced sampling strategies that combine h- and p-refinement [37], potentially with bifurcation detection, could be beneficial [42].

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

P-HT derived the formulation and error indicator, performed the numerical simulations, collected the data, analyzed and interpreted the results. PF supervised the project, provided critical feedback and helped on interpreted the results. Both authors prepared the manuscript.

Funding

This research is supported by the DOE Office of Nuclear Energy under the Nuclear Energy University Program (Proj. No. DE_NE0008780). Simulations were performed at the DOE Office of Science User Facility ALCF (Argonne Leadership Computing Facility).

Acknowledgments

We thank Anthony T. Patera (MIT) for his valuable comments, many contributions, insights and guidelines.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Technically, since these flows are not turbulent the TKE should be referred to as velocity variance. Because the is more widely used and the mathematical formulation is the same in either case, we prefer to use the more widely recognized appellation, TKE.

²The ROM coefficients in the steady problems are typically found through a Newton minimization over the POD approximation space [10, 11].

³The pMOR greedy strategy uses the maximal indicated error among the parametric training set to select the next anchor point. p-Greedy combines basis functions from FOMs at different anchor points to form an enriched approximation space; h-greedy builds an independent ROM for each anchor point Ngoc Cuong et al. [10].

⁴In [16], the velocity and temperature residuals are coupled because the velocity-temperature solutions are obtained through a coupled Newton’s method. Here, we do not couple the residuals because we solve Eqs. 21 and 22 separately.

⁵The O(N⁴) costs, which arise from the rank-3 advection tensor, might be mitigated by an O(N²) approximation to the advection operator, such as suggested in [41].

References

1. Grepl MA, Patera AT. A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations. ESAIM: Math Model Numer Anal (2005) 39:157–81. doi:10.1051/m2an:2005006

CrossRef Full Text | Google Scholar

2. Rozza G, Huynh DBP, Patera AT. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch Comput Methods Eng (2008) 15:229–75. doi:10.1007/s11831-008-9019-9

CrossRef Full Text | Google Scholar

3. Merzari E, Ninokata H, Mahmood A, Rohde M. Proper orthogonal decomposition of the flow in geometries containing a narrow gap. Theor Comput Fluid Dyn (2009) 23:333–51. doi:10.1007/s00162-009-0152-3

CrossRef Full Text | Google Scholar

4. Merzari E, Pointer WD, Fischer P. A POD-based solver for the advection-diffusion equation. Fluids Eng Division Summer Meet (2011) 44403:1139–47.

CrossRef Full Text | Google Scholar

5. Quarteroni A, Rozza G, Manzoni A. Certified reduced basis approximation for parametrized partial differential equations and applications. J Math Ind (2011) 1:3. doi:10.1186/2190-5983-1-3

CrossRef Full Text | Google Scholar

6. Wang Z, Akhtar I, Borggaard J, Iliescu T. Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison. Computer Methods Appl Mech Eng (2012) 237:10–26. doi:10.1016/j.cma.2012.04.015

CrossRef Full Text | Google Scholar

7. Kaneko K, Fischer P. Augmented reduced order models for turbulence. Front Phys (Forthcoming 2022).

Google Scholar

8. Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differential equations: An introduction. Switzerland: Springer (2015).

Google Scholar

9. Fick L, Maday Y, Patera AT, Taddei T. A stabilized POD model for turbulent flows over a range of Reynolds numbers: Optimal parameter sampling and constrained projection. J Comput Phys (2018) 371:214–43. doi:10.1016/j.jcp.2018.05.027

CrossRef Full Text | Google Scholar

10. Ngoc Cuong N, Veroy K, Patera AT. Certified real-time solution of parametrized partial differential equations. In: Handbook of materials modeling. Dordrecht, Netherlands: Springer (2005). p. 1529–64.

CrossRef Full Text | Google Scholar

11. Veroy K, Patera AT. Certified real-time solution of the parametrized steady incompressible Navier–Stokes equations: Rigorous reduced-basis a posteriori error bounds. Int J Numer Methods Fluids (2005) 47:773–88. doi:10.1002/fld.867

CrossRef Full Text | Google Scholar

12. Wang Q, Wan ZH, Yan R, Sun DJ. Multiple states and heat transfer in two-dimensional tilted convection with large aspect ratios. Phys Rev Fluids (2018) 3:113503. doi:10.1103/physrevfluids.3.113503

CrossRef Full Text | Google Scholar

13. Deparis S. Reduced basis error bound computation of parameter-dependent Navier–Stokes equations by the natural norm approach. SIAM J Numer Anal (2008) 46:2039–67. doi:10.1137/060674181

CrossRef Full Text | Google Scholar

14. Deparis S, Rozza G. Reduced basis method for multi-parameter-dependent steady Navier–Stokes equations: Applications to natural convection in a cavity. J Comput Phys (2009) 228:4359–78. doi:10.1016/j.jcp.2009.03.008

CrossRef Full Text | Google Scholar

15. Ballarin F, Rebollo TC, Ávila ED, Mármol MG, Rozza G. Certified Reduced Basis VMS-Smagorinsky model for natural convection flow in a cavity with variable height. Comput Mathematics Appl (2020) 80:973–89. doi:10.1016/j.camwa.2020.05.013

CrossRef Full Text | Google Scholar

16. Knezevic DJ, Nguyen NC, Patera AT. Reduced basis approximation and a posteriori error estimation for the parametrized unsteady Boussinesq equations. Math Models Methods Appl Sci (2011) 21:1415–42. doi:10.1142/s0218202511005441

CrossRef Full Text | Google Scholar

17. Yano M. A space-time Petrov–Galerkin certified reduced basis method: Application to the Boussinesq equations. SIAM J Sci Comput (2014) 36:A232–66. doi:10.1137/120903300

CrossRef Full Text | Google Scholar

18. Kaneko K, Tsai PH, Fischer P. Towards model order reduction for fluid-thermal analysis. Nucl Eng Des (2020) 370:110866. doi:10.1016/j.nucengdes.2020.110866

CrossRef Full Text | Google Scholar

19. Wells D, Wang Z, Xie X, Iliescu T. An evolve-then-filter regularized reduced order model for convection-dominated flows. Int J Numer Methods Fluids (2017) 84:598–615. doi:10.1002/fld.4363

CrossRef Full Text | Google Scholar

20. Karimi A, Paul MR. Quantifying spatiotemporal chaos in Rayleigh-Bénard convection. Phys Rev E (2012) 85:046201. doi:10.1103/physreve.85.046201

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wolf A, Swift JB, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D: Nonlinear Phenomena (1985) 16:285–317. doi:10.1016/0167-2789(85)90011-9

CrossRef Full Text | Google Scholar

22. Goldhirsch I, Sulem PL, Orszag SA. Stability and Lyapunov stability of dynamical systems: A differential approach and a numerical method. Physica D: Nonlinear Phenomena (1987) 27:311–37. doi:10.1016/0167-2789(87)90034-0

CrossRef Full Text | Google Scholar

23. Cross M, Hohenberg P. Spatiotemporal chaos. Science (1994) 263:1569–70. doi:10.1126/science.263.5153.1569

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Egolf DA, Melnikov IV, Pesch W, Ecke RE. Mechanisms of extensive spatiotemporal chaos in Rayleigh–Bénard convection. Nature (2000) 404:733–6. doi:10.1038/35008013

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cross MC, Hohenberg PC. Pattern formation outside of equilibrium. Rev Mod Phys (1993) 65:851–1112. doi:10.1103/revmodphys.65.851

CrossRef Full Text | Google Scholar

26. Busse F. Non-linear properties of thermal convection. Rep Prog Phys (1978) 41:1929–67. doi:10.1088/0034-4885/41/12/003

CrossRef Full Text | Google Scholar

27. Maday Y, Patera AT, Ronquist EM. A well-posed optimal spectral element approximation for the Stokes problem. Tech Rep (1987) 1987.

Google Scholar

28. Fischer PF, Lottes JW, Kerkemeier SG. Nek5000 web page (2008).

Google Scholar

29. Quarteroni A, Valli A. Numerical approximation of partial differential equations. In: Springer series in computational mathematics. Berlin: Springer (1994).

CrossRef Full Text | Google Scholar

30. Fischer P, Schmitt M, Tomboulides A. Recent developments in spectral element simulations of moving-domain problems. In: Recent progress and modern challenges in applied mathematics, modeling and computational science. Berlin: Springer (2017). p. 213–44.

CrossRef Full Text | Google Scholar

31. Rønquist EM, Patera AT. A Legendre spectral element method for the Stefan problem. Int J Numer Methods Eng (1987) 24:2273–99. doi:10.1002/nme.1620241204

CrossRef Full Text | Google Scholar

32. Maday Y, Patera AT, Rønquist EM. The PN × PN−2 method for the approximation of the Stokes problem. Paris: Laboratoire d’Analyse Numérique (1992).

Google Scholar

33. Guermond JL, Oden JT, Prudhomme S. Mathematical perspectives on large eddy simulation models for turbulent flows. J Math Fluid Mech (2004) 6:194–248. doi:10.1007/s00021-003-0091-5

CrossRef Full Text | Google Scholar

34. Guermond JL, Pasquetti R, Popov B. Entropy viscosity method for nonlinear conservation laws. J Comput Phys (2011) 230:4248–67. doi:10.1016/j.jcp.2010.11.043

CrossRef Full Text | Google Scholar

35. Mullen J. Development of a parallel spectral element based large eddy simulation model for the flow of incompressible fluids in complex geometries. Ph.D. thesis. Providence, Rhode Island: Brown University (1999).

Google Scholar

36. Sabetghadam F, Jafarpour A. α regularization of the POD-Galerkin dynamical systems of the Kuramoto–Sivashinsky equation. Appl Mathematics Comput (2012) 218:6012–26. doi:10.1016/j.amc.2011.11.083

CrossRef Full Text | Google Scholar

37. Eftang JL, Knezevic DJ, Patera AT. An hp certified reduced basis method for parametrized parabolic partial differential equations. Math Computer Model Dynamical Syst (2011) 17:395–422. doi:10.1080/13873954.2011.547670

CrossRef Full Text | Google Scholar

38. Haasdonk B, Ohlberger M. Reduced basis method for finite volume approximations of parametrized linear evolution equations. ESAIM: Math Model Numer Anal (2008) 42:277–302. doi:10.1051/m2an:2008001

CrossRef Full Text | Google Scholar

39. Haasdonk B. Convergence rates of the POD–Greedy method. ESAIM: Math Model Numer Anal (2013) 47:859–73. doi:10.1051/m2an/2012045

CrossRef Full Text | Google Scholar

40. Gelfgat AY, Bar-Yoseph P, Yarin A. Stability of multiple steady states of convection in laterally heated cavities. J Fluid Mech (1999) 388:315–34. doi:10.1017/s0022112099004796

CrossRef Full Text | Google Scholar

41. Barrault M, Maday Y, Nguyen NC, Patera AT. An empirical interpolation method: Application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique (2004) 339:667–72. doi:10.1016/j.crma.2004.08.006

CrossRef Full Text | Google Scholar

42. Pichi F. Reduced order models for parametric bifurcation problems in nonlinear PDEs. Ph.D. thesis. Trieste, Italy: Scuola Internazionale Superiore di Studi Avanzati (2020).