Certified Reduced Basis VMS-Smagorinsky model for natural convection flow in a cavity with variable height

TL;DR

This paper introduces a certified reduced basis model for natural convection in a variable height cavity, incorporating physical and geometric parameters, and demonstrates its effectiveness through numerical tests.

Contribution

The work develops a novel reduced basis VMS-Smagorinsky model with an empirical interpolation and error estimation for parametrized natural convection problems.

Findings

Effective reduced basis model for variable height cavities

Accurate approximation of sub-grid eddy viscosity terms

Successful numerical validation across parameter configurations

Abstract

In this work we present a Reduced Basis VMS-Smagorinsky Boussinesq model, applied to natural convection problems in a variable height cavity, in which the buoyancy forces are involved. We take into account in this problem both physical and geometrical parametrizations, considering the Rayleigh number as a parameter, so as the height of the cavity. We perform an Empirical Interpolation Method to approximate the sub-grid eddy viscosity term that let us obtain an affine decomposition with respect to the parameters. We construct an a posterior error estimator, based upon the Brezzi-Rappaz-Raviart theory, used in the greedy algorithm for the selection of the basis functions. Finally we present several numerical tests for different parameter configuration.

Figures

Tables

Table 1. Table 1: Computational time for FE and RB solutions, with the speedup and the error, for problem ( 1 ), R ​ a ∈ [ 10 3 , 10 5 ] 𝑅 𝑎 superscript 10 3 superscript 10 5 Ra\in[10^{3},10^{5}] (top) and R ​ a ∈ [ 10 5 , 10 6 ] 𝑅 𝑎 superscript 10 5 superscript 10 6 Ra\in[10^{5},10^{6}] (bottom), μ g = 1 subscript 𝜇 𝑔 1 \mu_{g}=1 .

Data ${\mu }_{ph}=4060$ ${\mu }_{ph}=17808$ ${\mu }_{ph}=53778$ ${\mu }_{ph}=93692$ $T_{FE}$ 633.65s 585.83s 553.25s 677.86s $T_{online}$ 0.55s 0.5s 0.46s 0.49s speedup 1133 1151 1189 1367 $\frac{{\Vert {\text{u}}_h-{\text{u}}_N\Vert }_1}{{\Vert {\text{u}}_h\Vert }_1}$ $4.01\cdot {10}^{-9}$ $4.06\cdot {10}^{-9}$ $3.56\cdot {10}^{-9}$ $3.22\cdot {10}^{-9}$ $\frac{{\Vert {\theta }_h-{\theta }_N\Vert }_1}{{\Vert {\theta }_h\Vert }_1}$ $5.59\cdot {10}^{-9}$ $4.29\cdot {10}^{-9}$ $5.36\cdot {10}^{-9}$ $4.76\cdot {10}^{-9}$ $\frac{{\Vert p_h-p_N\Vert }_0}{{\Vert p_h\Vert }_0}$ $5.75\cdot {10}^{-10}$ $1.71\cdot {10}^{-10}$ $1.96\cdot {10}^{-10}$ $2.30\cdot {10}^{-10}$

Table 2. Table 2: Computational time for FE and RB solutions, with the speed-up and the error, for Boussinesq VMS-Smagorinsky model with μ g ∈ [ 0.5 , 2 ] subscript 𝜇 𝑔 0.5 2 \mu_{g}\in[0.5,2] .

Data ${\mu }_g=0.64$ ${\mu }_g=1.08$ ${\mu }_g=1.44$ ${\mu }_g=1.87$ $T_{FE}$ 808.91s 810.16s 866.1s 851.82s $T_{online}$ 2.68s 2.55s 2.61s 2.52s speedup 301 317 331 337 ${\Vert {\text{u}}_h-{\text{u}}_N\Vert }_1/{\Vert {\text{u}}_h\Vert }_1$ $3.4\cdot {10}^{-9}$ $4.12\cdot {10}^{-9}$ $5.41\cdot {10}^{-9}$ $5.68\cdot {10}^{-9}$ ${\Vert {\theta }_h-{\theta }_N\Vert }_1/{\Vert {\theta }_h\Vert }_1$ $3.75\cdot {10}^{-8}$ $4.66\cdot {10}^{-9}$ $4.86\cdot {10}^{-9}$ $4.91\cdot {10}^{-9}$ ${\Vert p_h-p_N\Vert }_0/{\Vert p_h\Vert }_0$ $2.51\cdot {10}^{-9}$ $3.25\cdot {10}^{-9}$ $5.51\cdot {10}^{-9}$ $4.48\cdot {10}^{-9}$

Table 3. Table 3: Computational time for FE and RB solutions, with the speedup and the error, for Boussinesq VMS-Smagorinsky model with 𝝁 ∈ [ 10 3 , 10 4 ] × [ 0.5 , 2 ] 𝝁 superscript 10 3 superscript 10 4 0.5 2 \boldsymbol{\mu}\in[10^{3},10^{4}]\times[0.5,2] .

Data $Ra=2143$ $Ra=3506$ $Ra=5922$ $Ra=9618$ ${\mu }_g=1.95$ ${\mu }_g=0.71$ ${\mu }_g=1.13$ ${\mu }_g=1.63$ $T_{FE}$ 600.96s 914.18s 684.95s 630.94s $T_{online}$ 11.08s 15.73s 14.52s 11.46s speedup 54 58 47 55 ${\Vert {\text{u}}_h-{\text{u}}_N\Vert }_1/{\Vert {\text{u}}_h\Vert }_1$ $2.07\cdot {10}^{-7}$ $3.91\cdot {10}^{-7}$ $6.52\cdot {10}^{-8}$ $4.18\cdot {10}^{-8}$ ${\Vert {\theta }_h-{\theta }_N\Vert }_1/{\Vert {\theta }_h\Vert }_1$ $1.29\cdot {10}^{-7}$ $3.75\cdot {10}^{-7}$ $3.17\cdot {10}^{-8}$ $1.19\cdot {10}^{-8}$ ${\Vert p_h-p_N\Vert }_0{\Vert p_h\Vert }_0$ $1.56\cdot {10}^{-8}$ $4.21\cdot {10}^{-8}$ $5.57\cdot {10}^{-9}$ $1.52\cdot {10}^{-9}$

Equations

$$\left\{\begin{array}[]{ll}\textbf{u}_{o}\cdot\nabla_{o}\textbf{u}_{o}-Pr\Delta_{o}\textbf{u}_{o}-\nabla_{o}\cdot(\nu_{T}(\textbf{u}_{o})\nabla_{o}\textbf{u}_{o})\\ +\nabla_{o}p_{o}-Pr\,\mu_{ph}\,\theta_{o}\,\textbf{e}_{d}=\textbf{f}&\mbox{ in }\Omega_{o}(\mu_{g})\vspace{0.2cm}\\ \nabla_{o}\cdot\textbf{u}_{o}=0&\mbox{ in }\Omega_{o}(\mu_{g})\vspace{0.2cm}\\ \textbf{u}_{o}\cdot\nabla_{o}\theta_{o}-\Delta_{o}\theta_{o}-\nabla_{o}\cdot(K_{T}(\textbf{u}_{o})\nabla_{o}\theta_{o})=Q&\mbox{ in }\Omega_{o}(\mu_{g})\vspace{0.2cm}\\ \textbf{u}_{o}=0&\mbox{ on }\Gamma\vspace{0.2cm}\\ \theta_{o}=\theta_{D}&\mbox{ on }\Gamma_{D}\vspace{0.15cm}\\ \partial_{n}\theta_{o}=0&\mbox{ on }\Gamma_{N}\vspace{0.15cm}\\ \end{array}\right.$$

$$K_T(\text{u})=\frac{1}{Pr}{\nu }_T(\text{u}).$$

$$\Vert \text{v}{\Vert }_{0,4,\Omega }\le C_u\Vert \nabla \text{v}{\Vert }_{0,2,\Omega },\forall \text{v}\in Y,$$

$$\Vert \theta {\Vert }_{0,4,\Omega }\le C_{\theta }\Vert \nabla \theta {\Vert }_{0,2,\Omega },\forall \theta \in \Theta .$$

$$\Vert U{\Vert }_X^2=\Vert \nabla \text{u}{\Vert }_{0,2,\Omega }^2+\Vert \nabla \theta {\Vert }_{0,2,\Omega }^2+\Vert p{\Vert }_{0,2,\Omega }^2,\forall U=(\text{u},\theta ,p)\in X.$$

$$\left\{\begin{array}[]{l}\mbox{Find }(\textbf{u}_{o},\theta_{o}^{u},p_{o}^{u})=(\textbf{u}_{o}(\boldsymbol{\mu}),\theta^{u}_{o}(\boldsymbol{\mu}),p^{u}_{o}(\boldsymbol{\mu}))\in X^{o}\mbox{ such that}\vspace{0.3cm}\\ \begin{array}[]{ll}\tilde{a}_{u}(\textbf{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})+\tilde{b}(\textbf{v}_{o},p^{u}_{o};\boldsymbol{\mu})+\tilde{a}_{Su}^{\prime}(\textbf{u}_{o};\textbf{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})\\ +\tilde{c}_{u}(\textbf{u}_{o},\textbf{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})+\tilde{f}(\theta_{o}^{u},\textbf{v}_{o};\boldsymbol{\mu})=\tilde{F}(\textbf{v}_{o};\boldsymbol{\mu})&\quad\forall\textbf{v}_{o}\in Y^{o},\vspace{0.15cm}\\ \tilde{b}(\textbf{u}_{o},p^{v}_{o};\boldsymbol{\mu})=0&\quad\forall p^{v}_{o}\in M^{o},\vspace{0.15cm}\\ \tilde{a}_{\theta}(\theta^{u}_{o},\theta^{v}_{o};\boldsymbol{\mu})+\tilde{c}_{\theta}(\textbf{u}_{o},\theta_{o}^{u},\theta^{v}_{o};\boldsymbol{\mu})\\ +\tilde{a}_{S\theta,n}^{\prime}(\textbf{u}_{o};\theta^{u}_{o},\theta^{v}_{o};\boldsymbol{\mu})=\tilde{Q}(\theta^{v}_{o};\boldsymbol{\mu})&\quad\forall\theta^{v}_{o}\in\Theta^{o}.\end{array}\end{array}\right.$$

$$\begin{array}[]{l}\tilde{a}_{u}(\textbf{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})=Pr\displaystyle\int_{\Omega_{o}(\mu_{g})}\nabla_{o}\textbf{u}_{o}:\nabla_{o}\textbf{v}_{o}\;d\Omega_{o},\\ \tilde{a}_{\theta}(\theta_{o}^{u},\theta_{o}^{v};\boldsymbol{\mu})=\displaystyle\int_{\Omega_{o}(\mu_{g})}\nabla_{o}\theta_{o}^{u}\cdot\nabla_{o}\theta_{o}^{v}\;d\Omega_{o},\vspace{0.2cm}\\ \tilde{f}(\theta^{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})=-Pr\,\mu_{ph}\displaystyle\int_{\Omega_{o}(\mu_{g})}\theta^{u}_{o}\,\textbf{v}_{do}\;d\Omega_{o},\\ \tilde{b}(\textbf{u}_{o},p^{v}_{o};\boldsymbol{\mu})=-\displaystyle\int_{\Omega_{o}(\mu_{g})}(\nabla_{o}\cdot\textbf{u}_{o})p_{o}^{v}\;d\Omega_{o};\end{array}$$

$$\begin{array}[]{l}\tilde{c}_{u}(\textbf{z}_{o},\textbf{u}_{o},\textbf{v}_{o};\boldsymbol{\mu})=\displaystyle\int_{\Omega_{o}(\mu_{g})}(\textbf{z}_{o}\cdot\nabla_{o}\textbf{u}_{o})\textbf{v}_{o}\;d\Omega_{o},\\ \tilde{c}_{\theta}(\textbf{u}_{o},\theta^{u}_{o},\theta^{v}_{o};\boldsymbol{\mu})=\displaystyle\int_{\Omega_{o}(\mu_{g})}(\textbf{u}_{o}\cdot\nabla_{o}\theta^{u}_{o})\theta^{v}_{o}\;d\Omega_{o};\end{array}$$

$${\tilde{a}}_{Su}^{\prime }({\text{z}}_o;{\text{u}}_o,{\text{v}}_o;\mathbf{\mu })={\int }_{{\Omega }_o({\mu }_g)}{\nu }_T({\text{z}}_o)\nabla {\text{u}}_o:\nabla {\text{v}}_{o}d{\Omega }_o,$$

$${\nu }_T(\text{u})=(C_S{h}_K{)}^2|\nabla {\text{u}}_{{|}_K}|.$$

$${\varphi }_n(x)=\frac{n}{\Vert \varphi {\Vert }_{0,1,\mathbb{R}}}\varphi (nx).$$

$${\tilde{a}}_{S\theta ,n}^{\prime }({\text{u}}_o;{\theta }_o^u,{\theta }_o^v;\mathbf{\mu })={\int }_{{\Omega }_o({\mu }_g)}{\nu }_{T,n}({\text{u}}_o){\nabla }_o{\theta }_o^u\cdot {\nabla }_o{\theta }_o^{v}d{\Omega }_o,$$

$${\nu }_{T,n}(\text{u})=(C_S{h}_K{)}^2({\varphi }_n\ast |\nabla {\text{u}}_{{|}_K}|),$$

$${\tilde{a}}_{S\theta }^{\prime }({\text{u}}_o;{\theta }_o^u,{\theta }_o^v;\mathbf{\mu })={\int }_{{\Omega }_o({\mu }_g)}{\nu }_T({\text{u}}_o)\nabla {\theta }_o^u\cdot \nabla {\theta }_o^{v}d{\Omega }_o.$$

$$\tilde{F}({\text{v}}_o;\mathbf{\mu })=⟨\text{f},{\text{v}}_o⟩,\tilde{Q}({\theta }_o^v;\mathbf{\mu })=⟨Q,{\theta }_o^v⟩,$$

$$\Psi((x,y);\mu_{g})=\left(\begin{array}[]{cc}1&0\\ 0&\mu_{g}\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right),\quad\forall(x,y)\in\Omega_{r}.$$

$$\textbf{J}((x,y);\mu_{g})=\left(\begin{array}[]{cc}1&0\\ 0&\mu_{g}\end{array}\right),\mbox{ and }|\textbf{J}((x,y);\mu_{g})|=\mu_{g}.$$

$$\left\{\begin{array}[]{l}\mbox{Find }(\textbf{u}_{h},\theta_{h}^{u},p_{h}^{u})=(\textbf{u}_{h}(\boldsymbol{\mu}),\theta^{u}_{h}(\boldsymbol{\mu}),p^{u}_{h}(\boldsymbol{\mu}))\in X_{h}\mbox{ such that}\vspace{0.3cm}\\ \begin{array}[]{ll}a_{u,x}(\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})+a_{u,y}(\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})+b_{x}(\textbf{v}_{h},p^{u}_{h};\boldsymbol{\mu})\\ +b_{y}(\textbf{v}_{h},p^{u}_{h};\boldsymbol{\mu})+a_{Su,x}^{\prime}(\textbf{u}_{h};\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})\\ +a_{Su,y}^{\prime}(\textbf{u}_{h};\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})+c_{u,x}(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})\\ +c_{u,y}(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})+f(\theta_{h}^{u},\textbf{v}_{h};\boldsymbol{\mu})=F(\textbf{v}_{h};\boldsymbol{\mu})&\quad\forall\textbf{v}_{h}\in Y_{h},\vspace{0.15cm}\\ b_{x}(\textbf{u}_{h},p^{v}_{h};\boldsymbol{\mu})+b_{y}(\textbf{u}_{h},p^{v}_{h};\boldsymbol{\mu})=0&\quad\forall p^{v}_{h}\in M_{h},\vspace{0.15cm}\\ a_{\theta,x}(\theta^{u}_{h},\theta^{v}_{h};\boldsymbol{\mu})+a_{\theta,y}(\theta^{u}_{h},\theta^{v}_{h};\boldsymbol{\mu})+c_{\theta,x}(\textbf{u}_{h},\theta_{h}^{u},\theta^{v}_{h};\boldsymbol{\mu})\\ +c_{\theta,y}(\textbf{u}_{h},\theta_{h}^{u},\theta^{v}_{h};\boldsymbol{\mu})+a_{S\theta,nx}^{\prime}(\textbf{u}_{h};\theta^{u}_{h},\theta^{v}_{h};\boldsymbol{\mu})\\ +a_{S\theta,ny}^{\prime}(\textbf{u}_{h};\theta^{u}_{h},\theta^{v}_{h};\boldsymbol{\mu})=Q(\theta^{v}_{h};\boldsymbol{\mu})&\quad\forall\theta^{v}_{h}\in\Theta_{h},\end{array}\end{array}\right.$$

$$\Vert {\Pi }_h^{\ast }{\text{u}}_h{\Vert }_{1,2,\Omega }\le C_f\Vert {\text{u}}_h{\Vert }_{1,2,\Omega }{\text{u}}_h\in Y_h,$$

$$a_{u,x}({\text{u}}_h,{\text{v}}_h;\mathbf{\mu })=Pr{\mu }_g{\int }_{{\Omega }_r}({\partial }_x{u}_1{\partial }_x{v}_1+{\partial }_x{u}_2{\partial }_x{v}_2)d{\Omega }_r,$$

$$a_{u,y}({\text{u}}_h,{\text{v}}_h;\mathbf{\mu })=\frac{Pr}{{\mu }_g}{\int }_{{\Omega }_r}({\partial }_y{u}_1{\partial }_y{v}_1+{\partial }_y{u}_2{\partial }_y{v}_2)d{\Omega }_r,$$

$$b_x({\text{u}}_h,p_h^u;\mathbf{\mu })=-{\mu }_g{\int }_{{\Omega }_r}{p}_h^u{\partial }_x{u}_{1}d{\Omega }_r,b_y({\text{u}}_h,p_h^u;\mathbf{\mu })=-{\int }_{{\Omega }_r}{p}_h^u{\partial }_y{u}_{2}d{\Omega }_r,$$

$$f({\theta }_h^u,{\text{v}}_h;\mathbf{\mu })=Pr{\mu }_{ph}{\mu }_g{\int }_{{\Omega }_r}{\theta }_h^u{v}_{2}d{\Omega }_r,$$

$$a_{\theta ,x}({\theta }_h^u,{\theta }_h^v;\mathbf{\mu })={\mu }_g{\int }_{{\Omega }_r}{\partial }_x{\theta }_h^u{\partial }_x{\theta }_h^{v}d{\Omega }_r,a_{\theta ,y}({\theta }_h^u,{\theta }_h^v;\mathbf{\mu })=\frac{1}{{\mu }_g}{\int }_{{\Omega }_r}{\partial }_y{\theta }_h^u{\partial }_y{\theta }_h^{v}d{\Omega }_r,$$

$$c_{u,x}(\textbf{w}_{h},\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})=\mu_{g}\displaystyle\int_{\Omega_{r}}\big{[}(w_{1}\partial_{x}u_{1})v_{1}+(w_{1}\partial_{x}u_{2})v_{2}\big{]}\;d\Omega_{r},$$

$$c_{u,y}(\textbf{w}_{h},\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})=\displaystyle\int_{\Omega_{r}}\big{[}(w_{2}\partial_{y}u_{1})v_{1}+(w_{2}\partial_{y}u_{2})v_{2}\big{]}\;d\Omega_{r},$$

$$c_{\theta ,x}({\text{u}}_h,{\theta }_h^u,{\theta }_h^v;\mathbf{\mu })={\mu }_g{\int }_{{\Omega }_r}(u_1{\partial }_x{\theta }_h^u){\theta }_h^{v}d{\Omega }_r,$$

$$c_{\theta ,y}({\text{u}}_h,{\theta }_h^u,{\theta }_h^v;\mathbf{\mu })={\int }_{{\Omega }_r}(u_2{\partial }_y{\theta }_h^u){\theta }_h^{v}d{\Omega }_r,$$

$$a_{Su,x}^{\prime}(\textbf{w}_{h};\textbf{u}_{h},\textbf{v}_{h};\boldsymbol{\mu})=\mu_{g}\int_{\Omega_{r}}{\nu_{T}(\Pi_{h}^{*}\textbf{w})\big{[}\partial_{x}(\Pi_{h}^{*}u_{1})\partial_{x}(\Pi_{h}^{*}v_{1})}$$

$${+\partial_{x}(\Pi_{h}^{*}u_{2})\partial_{x}(\Pi_{h}^{*}v_{2})d\Omega_{r}\big{]}},$$

Full text

**Keywords. ** Reduced basis method, Empirical interpolation method, a posteriori error estimation, Boussinesq equations, Smagorinsky turbulence model.

1 Introduction

Nowadays, several industrial processes need numerical simulations, which are usually performed with the widespread high-fidelity approximation techniques such as finite element (FE), finite volumes or spectral methods, and they usually take very long time for computation (several hours, even days). In many situations, the model that represents the behavior of an industrial process is given by a Partial Differential Equation (PDE) depending on parameters. Reduced-order modeling (ROM) is used in parametrized PDE in order to try to reduce the high computational time required by its numerical solution, when large number of simulations with different parameter values are needed [1, 2, 3, 4, 5, 6, 7, 8].

In the context of fluids dynamics, even with the reduction of the computational cost provided by turbulent models, such as the Variational Multi-Scale (VMS) models (cf. [9]), with respect to the DNS, it is still expensive to compute accurately the real flows that commonly appear in industry problems, specially in cases where parameters play important roles. When a high number of computations for a fluid flow depending on parameters is required, ROM becomes useful. Several works for Reduced Basis (RB) model have been presented for Stokes equations [10, 6, 11] and Navier-Stokes equations [12, 13, 14, 15, 16]. Most recently, there have been developed works for RB turbulent models, such as the Smagorinsky model [17], or the VMS-Smagorinsky model [18]. In those last works, the ROM is constructed from the turbulent model. A different way to build ROM model for turbulent flows is the one in which initially we would build a ROM for Navier-Stokes equations, and then model the unresolved scales (either by eddy diffusion or other techniques) to build the turbulence model. This approach has been followed in [19, 20], for instance. We address here the construction of a RB model for the Boussinesq equations, that includes turbulent diffusivity for both momentum and energy equations. The turbulent diffusivities are modeled by a VMS-Smagorinsky approach, in such a way that eddy diffusion effects only act on the small resolved scales. One of the main advantages of using the VMS-Smagorinsky is that the eddy viscosity and eddy diffusivity only affect the small resolved scales, avoiding over-diffusive effects, that can be maintained in the ROM setting thanks to the Reduced Basis framework that we have developed. On the counterpart, we deal with non-linear terms related with the VMS-Smagorinsky model, that force us to use techniques such as the Empirical Interpolation Method, resulting a feasible ROM model.

In this work we consider the application of a RB Boussinesq VMS-Smago-

rinsky model to simulate a natural convection in a variable height cavity. In applications to architecture, this cavity represents a courtyard inside a building. The study of the heat exchange between the air and the walls inside the courtyard is of high interest to minimize the energy needs of the building. The variability of the cavity height is considered through a geometrical parametrization of the domain. Since we are interested in solving efficiently the parameter-dependent problem, we need to reformulate the Boussinesq VMS-Smagorinsky model in a parameter-independent domain with a change of variables. With this change of variables, we obtain operators that depend on both physical and geometrical parameters. This setting, besides the Empirical Interpolation Method (EIM) (cf. [21, 22]) for the non-linear eddy diffusivities, lets us approximate these non-linear operators by operators that depend affinely with respect to the parameters, both of physical and geometrical type. Then it is possible to store parameter-independent matrices and tensors in the offline phase.

We tackle in this work some of the intermediate difficulties which should necessarily be solved when dealing with RB modelling of natural convection turbulent flows. Particularly, we analyze how to deal with the temperature when constructing the a posteriori error estimator. From the numerical analysis point of view, we present the development of an a posteriori error bound based upon the Brezzi-Rappaz-Raviart (BRR) theory [23], used in the snapshot selection in the greedy algorithm. This a posteriori error estimator is an extension of the ones presented for the Navier-Stokes equations [12, 13, 14] and the Smagorinsky model [17]. The main difference for the a posteriori error bound presented in this paper with the previous ones, is the necessity of considering a mollifier for the thermal eddy diffusivity term, due to the fact this term is no longer Lipschitz-continuous. Thanks to the consideration of this regularized term, we are able to develop the a posteriori error estimator.

We present four different numerical tests for the buoyancy-driven cavity problem. In the first two tests, we consider a fixed height, i.e. we consider the geometrical parameter $\mu_{g}=1$, for different ranges of the Rayleigh number, ranging from moderate Rayleigh numbers values $Ra\in[10^{4},10^{5}]$, to high Rayleigh numbers values $Ra\in[10^{5},10^{6}]$. In the third one, we fix the Rayleigh number, with a moderate value $Ra=10^{5}$, and only the geometrical parameter changes. This test intends to represent a situation in which the environmental conditions are fixed and we are only interested in simulating the flow in cavities. In the last test, both the Rayleigh and the geometric parameter are taken into account. This test is more complex since two parameters are considered, thus the number of basis functions to include in our RB spaces increases with respect to the previous one. At the same time, the speed-up ratio decreases somewhat, although it remains on values around 50.

The paper is structured as follows: in section 2 we define the high fidelity problem in the reference domain from the one defined in the original domain depending on the geometric parameter. Then, in section 3, we present the Reduced Basis problem, with the EIM approximation for the eddy viscosity and eddy diffusivity terms. In section 4 we construct the a posteriori error estimator. Finally in section 5, we present the numerical results for the tests previously described, programmed in FreeFem++ (cf. [24]). Conclusions are then summarized in section 6.

2 Continuous problem and full order discretization

The aim of this work is to present a reduced order model for natural convection problems over domains with variable geometry. For this purpose, we propose a turbulence Smagorinsky model in which the buoyancy forces are modeled by the Boussinesq approach. Let $\boldsymbol{\mu}=(\mu_{ph},\mu_{g})\in\mathbb{R}^{p_{ph}\times p_{g}}$ and let $\Omega_{o}(\mu_{g})$ be a bounded polyhedral domain in $\mathbb{R}^{d}$ $(d=2,3)$, depending on the geometrical parameters, commonly called original domain in the RB framework. Let $\Gamma(\mu_{g})=\Gamma_{D}(\mu_{g})\cup\Gamma_{N}(\mu_{g})$ the Lipschitz-continuous boundary of $\Omega_{o}(\mu_{g})$, where $\Gamma_{D}$ is the part of the boundary with Dirichlet conditions and $\Gamma_{N}$ the part of the boundary with Neumann conditions.

We next present the continuous Boussinesq-Smagorinsky model that we consider in this work. Although the Smagorinsky approach is intrinsically discrete, we present it in a continuous form in order to clarify its relationship with the standard Boussinesq model:

[TABLE]

Here, $\textbf{u}_{o}$ is the velocity field, $p_{o}$ is the pressure and $\theta_{o}$ is the temperature. In addition, $\textbf{e}_{d}$ is the last vector of the canonical basis of $\mathbb{R}^{d}$, while $\mu_{ph}$ and $Pr$ are the Rayleigh and Prandtl dimensionless numbers respectively. Both the external body forces f, and the heat source term $Q$, are given data for the problem. In (1), $\nu_{T}(\textbf{u})$ is the eddy viscosity term, and $K_{T}(\textbf{u})$ is the eddy diffusivity given by

[TABLE]

In (1), we represent by $\theta_{D}$ a given temperature over the boundary $\Gamma_{D}$. For simplicity of the analysis, we further consider that $\theta_{D}=0$. In the case of considering non-homogeneous boundary conditions, it is enough to define a lift function $\theta_{g}$ such that $\theta_{g}|_{\Gamma_{D}}=\theta_{D}$. In [17], an analysis with the lift function is already done for a RBM Smagorinsky model.

Let us consider the spaces $Y^{o}=(H_{0}^{1}(\Omega_{o}))^{d}$, $M^{o}=L_{0}^{2}(\Omega_{o})$, $\Theta^{o}=H_{0}^{1}(\Omega_{o})$. We denote by $\|\cdot\|_{k,p,\Omega}$ the norm of the Sobolev space $W^{k,p}(\Omega)$. We consider the $H^{1}$-seminorm for the velocity and temperature spaces, and the $L^{2}$-norm for the pressure space. In addition, let us define the Sobolev embedding constants $C_{u}$ and $C_{\theta}$, associated to these norms, such that

[TABLE]

and

[TABLE]

Moreover, we consider the tensor space $X^{o}=Y^{o}\times\Theta^{o}\times M^{o}$, with the following associated norm:

[TABLE]

The variational formulation of problem (1), over the parameter-dependent original domain is

[TABLE]

Here, the bilinear forms $\tilde{a}_{u}(\cdot,\cdot;\boldsymbol{\mu})$, $\tilde{a}_{\theta}(\cdot,\cdot;\boldsymbol{\mu})$, $\tilde{b}(\cdot,\cdot;\boldsymbol{\mu})$ and $\tilde{f}(\cdot,\cdot;\boldsymbol{\mu})$ are defined by

[TABLE]

the trilinear forms $\tilde{c}_{u}(\cdot,\cdot,\cdot;\boldsymbol{\mu})$ and $\tilde{c}_{\theta}(\cdot,\cdot,\cdot;\boldsymbol{\mu})$ are defined by

[TABLE]

and the non-linear Smagorinsky term for eddy viscosity, $\tilde{a}_{Su}^{\prime}(\cdot;\cdot,\cdot;\boldsymbol{\mu})$, is given by

[TABLE]

with

[TABLE]

For the thermal eddy diffusivity term $\tilde{a}^{\prime}_{S\theta,n}$, let us first introduce a mollifier $\phi\in C_{c}^{\infty}(\mathbb{R})$, with supp($\phi$)$\subset B(0,1),\phi\geq 0,\|\phi\|_{0,1,\mathbb{R}}>0$, and $\phi$ is even, i.e., $\phi(-x)=\phi(x)$. Let us consider the mollifier sequence $\{\phi_{n}(x)\}_{n\geq 1}$, with $\phi_{n}\in C_{c}^{\infty}(\mathbb{R})$, supp($\phi_{n}$)$\subset B(0,1/n),\phi_{n}\geq 0,\|\phi_{n}\|_{0,1,\mathbb{R}}=1$, defined by

[TABLE]

Thus, the VMS-Smagorinsky eddy diffusivity term is defined as

[TABLE]

with

[TABLE]

where $*$ denotes the convolution. The choice of the eddy viscosity is the one suggested in [25]. For thermal eddy diffusivity, we have chosen a modified form, considering a mollifier for the thermal eddy diffusivity in [25]. The choice of this thermal eddy diffusivity with the mollifier, assure the Lipschitz continuity of the eddy diffusivity operator.

Thanks to mollifier properties (see [26] for details), it holds that $\tilde{a}_{S\theta,n}^{\prime}$ converges uniformly to $\tilde{a}_{S\theta}^{\prime}$, with

[TABLE]

Remark 1

*The mollified eddy diffusivity $\tilde{a}_{S\theta,n}^{\prime}$ is considered for the

well-possedness analysis and the development of a posteriori error estimator. In practice, we consider the eddy diffusivity term in the Boussinesq-Smagorinsky model as $\tilde{a}_{S\theta}^{\prime}$.*

Finally, the linear forms $\tilde{F}$ and $\tilde{Q}$ are given by

[TABLE]

where here $\left<\cdot,\cdot\right>$ stands either the duality paring between $Y$ and $Y^{\prime}$, and between $\Theta$ and $\Theta^{\prime}$; being $Y^{\prime}$ and $\Theta^{\prime}$ the dual spaces of $Y$ and $\Theta$, respectively.

We present a cavity domain, whose height is varied through a geometrical parameter. This geometrical parameter, that we denote by $\mu_{g}$, varies the aspect ratio of the cavity. In Fig. 1 we show the original cavity domain considered, with the geometrical parameter considered.

To be able to store parameter independent matrices in the offline phase (see section 3) of the RB method, we need to compute all the integrals in a reference domain trough a transformation of the original domain. Thus, we set $\mu_{g}^{ref}=1$, and we define the reference domain $\Omega_{r}=\Omega_{o}(\mu_{g}^{ref})$. The parameter-dependent original domain can be recovered by a transformation map, $\Psi:\Omega_{r}\times\mathcal{D}\rightarrow\mathbb{R}^{2}$, defined as

[TABLE]

As this map is linear, its Jacobian matrix and its determinant are given by

[TABLE]

Let $\{\mathcal{T}_{h}\}_{h>0}$ a regular family of triangulations on the reference domain. We denote $Y_{h}=(H_{0}^{1}(\Omega_{r})\cap V_{h}^{l}(\Omega_{r}))^{d}$, $M_{h}=L_{0}^{2}(\Omega_{r})\cap V_{h}^{m}(\Omega_{r})$, $\Theta_{h}=H_{0}^{1}(\Omega_{r})\cap V_{h}^{q}(\Omega_{r})$, with $l,m,q\in\mathbb{N}$, respectively the discrete velocity, pressure and temperature spaces on the reference domain. Although the analysis can be performed for a regular mesh, for simplicity in the analysis we suppose that the mesh we are considering is uniform. Denoting $X_{h}=Y_{h}\times\Theta_{h}\times M_{h}$, we rewrite problem (6) with respect to the reference domain, applying the change of variables of the transformation map $T$, as

[TABLE]

where the subscripts $x$ and $y$ denote the addend of the corresponding operator, relative to the partial derivative with respect to $x$ or $y$, respectively. For both the eddy viscosity and eddy diffusivity terms, we consider a VMS small-small setting approach (cf. [9]). For that, we consider a uniformly $H^{1}$-norm interpolator operator from $Y_{h}$ to $\overline{Y}_{h}=(H_{0}^{1}(\Omega_{r})\cap V_{h}^{l-1}(\Omega_{r}))^{d}$, denoted by $\Pi_{h}$. Thus, we assume that $\Pi_{h}$ satisfies that there exists a constant $C_{f}>0$ independent of $h$ such that

[TABLE]

where we denote $\Pi_{h}^{*}=Id-\Pi_{h}$. See [27, 28] for more details.

The operators in (14) have the following form:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These integrals are derived applying the well-known change of variable formula (see e.g. [5]). With this geometrical parametrization, the eddy viscosity $\nu_{T}(\cdot)$ (analogously $\nu_{T,n}(\cdot)$) also depends on the geometrical parameter, and is defined as

[TABLE]

Here we are supposing that we consider an uniform mesh in the reference domain $\Omega_{r}$, with $N_{h}$ partitions on each side. Since the mesh size, $h_{K}$, in the VMS-Smagorinsky eddy viscosity and eddy diffusivity appears in (16) in terms of the parameter-dependent original domain, we map it to the reference domain, by applying the change of variable map $\Psi$ defined in (12).

3 Reduced Basis formulation

In this section we present the RB problem derived from the discrete problem presented in section 2. We construct the low-dimensional spaces for the RB problem with the Greedy algorithm. Both the pressure and the temperature reduced basis spaces are defined with the corresponding snapshots computed solving the FE problem (14).

The RB velocity space is constructed with the velocity snapshot of the FE velocity solution, and the inner pressure supremizer (cf. [29, 30]), $T_{p}^{\boldsymbol{\mu}}:M_{h}\rightarrow Y_{h}$, defined for this problem as

[TABLE]

Thus, the reduced basis spaces are given by

[TABLE]

[TABLE]

[TABLE]

Denoting $X_{N}=Y_{N}\times\Theta_{N}\times M_{N}$, the RB problem is

[TABLE]

The eddy viscosity $\nu_{T}(\textbf{u};\mu_{g})$ must be tensorized in problem (21), for the efficient solve in the online phase. For this purpose, we consider the use of EIM (cf. [21, 22]). The eddy viscosity and eddy diffusivity terms are approximated as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\sigma_{k}(\boldsymbol{\mu})$ and $q_{k}$ are computed by the EIM algorithm (see [22] for further details).

The parameter independent matrices and tensors to store during the offline phase in order to efficiently solve problem (21), are given in this case by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here we are representing the reduced basis velocity, temperature and pressure solutions as a linear combination of the velocity, temperature and pressure snapshots, respectively, of the reduced spaces, i.e.,

[TABLE]

4 A posteriori error estimator

In order to develop the a posteriori error estimator for the Greedy algorithm, we rewrite problem (14) in a more compact form as

[TABLE]

The a posteriori error estimator is based upon the BRR theory (cf. [23]). For this purpose we define the Gateaux derivative of $A(\cdot,\cdot;\boldsymbol{\mu})$ with respect to the first variable, in the direction $Z\in X$, denoted by $\partial_{1}A(U,V;\boldsymbol{\mu})(Z)$. For this problem, denoting $Z=(\textbf{z},\theta^{z},p^{z})$, the derivative is defined by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

and

[TABLE]

The Gateaux derivative satisfies the following continuity and inf-sup conditions:

[TABLE]

[TABLE]

The existence of $\gamma_{0}\in\mathbb{R}$ and $\beta_{0}>0$ satisfying (23) and (24), respectively, are given by the following results, whose proofs can be found in A.1 and A.2 respectively.

Proposition 1

There exists $\gamma_{0}\in\mathbb{R}$ such that $\forall\mu\in\mathcal{D}$

[TABLE]

Proposition 2

Let $C(\boldsymbol{\mu},\|\phi_{n}^{\prime}\|_{0,1,\mathbb{R}})$ a constant depending on $\boldsymbol{\mu}$ and $\|\phi_{n}^{\prime}\|_{0,1,\mathbb{R}}$. Then if

[TABLE]

and

[TABLE]

then there exists $\tilde{\beta}(\boldsymbol{\mu})>0$ such that

[TABLE]

Since operator $b(\textbf{v}_{h},p_{h};\boldsymbol{\mu})$ satisfies the discrete inf-sup condition

[TABLE]

we can prove that the inf-sup condition (27) is satisfied thanks to Prop. 2.

With the following result, we prove that the Gateaux derivative of the Boussinesq-Smagorinsky operator is locally Lipschitz-continuous. The proof can be found in A.3.

Lemma 1

Let $U_{h}^{1},U_{h}^{2}\in X_{h}$. Then, in a neighborhood of $U_{h}^{1}$ and $U_{h}^{2}$, there exists a positive constant $\rho_{n}(\mu_{g})$ such that, $\forall Z_{h},V_{h}\in X_{h}$,

[TABLE]

We define the following continuity and inf-sup constants:

[TABLE]

[TABLE]

where the supremizer operator $T_{N}$ is defined as

[TABLE]

such that

[TABLE]

The existence of these constants can be proved in the same way that the existence of the constants (23)-(24). Thus, we can define the a posteriori error estimator as

[TABLE]

where $\tau_{N}(\boldsymbol{\mu})$ is given by

[TABLE]

with $\epsilon_{N}(\boldsymbol{\mu})$ the dual norm of the residual. The a posteriori error estimator is stated by the following result, whose proof can be found in A.4.

Theorem 1

Let $\boldsymbol{\mu}\in\mathcal{D}$, and assume that $\beta_{N}(\boldsymbol{\mu})>0$. If problem (22) admits a solution $U_{h}(\boldsymbol{\mu})$ such that

[TABLE]

then this solution is unique in the ball $B_{X}\left(U_{N}(\boldsymbol{\mu}),\dfrac{\beta_{N}(\boldsymbol{\mu})}{\rho_{n}(\mu_{g})}\right)$.

Moreover, assume that $\tau_{N}(\boldsymbol{\mu})\leq 1$ for all $\boldsymbol{\mu}\in\mathcal{D}$. Then there exists a unique solution $U_{h}(\boldsymbol{\mu})$ of (22) such that the error with respect $U_{N}(\boldsymbol{\mu})$, solution of (21), is bounded by the a posteriori error estimator, i.e.,

[TABLE]

with effectivity

[TABLE]

5 Numerical results

In this section we present some numerical results for the Boussinesq VMS-Smagorinsky RB model. We perform three different configurations for the parametrical set. The first configuration corresponds to the consideration only of physical parametrical set, fixing the value of the geometrical parameter $\mu_{g}=1$. Here we consider two different scenarios depending on the Rayleigh number range. In order to get error levels small enough for taking into account the a posteriori error estimator, we split the Rayleigh number ranger considered, $\mu_{ph}\in[10^{3},10^{6}]$, into two ranges, $\mu_{ph}\in[10^{3},10^{5}]$ and $\mu_{ph}\in[10^{5},10^{6}]$.

Then, we suppose that the Rayleigh number is fixed with $\mu_{ph}=10^{5}$, and we consider the geometrical parameter ranging in $\mu_{g}\in\mathcal{D}=[0.5,2]$. Finally, we consider both the geometrical parameter and the Rayleigh number. For this test, we consider the Rayleigh number, $\mu_{ph}$, ranging in $[10^{3},10^{4}]$, and the geometrical parameter, $\mu_{g}$, ranging in $\mu_{g}\in[0.5,2]$. Thus we are considering that the parameter domain is $\mathcal{D}=[10^{3},10^{4}]\times[0.5,2]$. For all cases, the Prandtl number considered is $Pr=0.71$, that corresponds to the air Prandtl number.

In all tests, we consider no-slip boundary conditions for velocity. We also consider homogeneous Neumann boundary condition for temperature at the top and bottom of the cavity, and Dirichlet conditions for the vertical walls: $\theta=1$ for the left vertical wall and $\theta=0$ for the right vertical wall. Moreover, we consider $\textbf{f}=0$ and $Q=0$ both for the momentum equation and the energy equation, respectively.

The FE solution is computed through a semi-implicit evolution approach, considering that the steady state solution is reached when the error between two iterates is below $\varepsilon_{FE}=10^{-10}$. The FE solution has been computed considering $\mathbb{P}2-\mathbb{P}2-\mathbb{P}1$ finite elements for velocity, temperature and pressure, respectively.

With respect to the constants that appear in the *a posteriori * error estimator, we explain in the following the numerical approximation of those constants. For the constant $\beta$, we use the Radial Basis Function (RBF) algorithm in order to compute efficiently $\beta(\mu)$ for all $\mu\in\cal{D}$. We follow the technique suggested in [31]. Although it can not be proved that a lower bound of the constant $\beta$ as the Successive Constraint Method (SCM) is provided, the computational time is much lower for the RBF, mainly when more than one parameter is considered, with quite good accuracy. The constant $\rho_{n}$, and more precisely the Sobolev embedding constants, are built once in the reference domain following the algorithm proposed in [12]. Moreover, we compute the Lipschitz constant $\rho_{n}(\mu_{g})$ without taking into account the term in which the mollifier takes part of it. This term in $\rho_{n}(\mu_{g})$ that comes from the mollifier is multiplied by a factor of $(C_{S}h)^{2}$, that has no relevance in the value of the Lipschitz constant. Thus, we consider accurate the approximation of the Lipschitz constant without considering the mollifier.

5.1 Physical parametrization

In this test, we consider two different scenarios, one for the Rayleigh number range $[10^{3},10^{5}]$, an the other for the Rayleigh number range $[10^{5},10^{6}]$. In both cases we consider the geometrical parameter fixed, with $\mu_{g}=1$. For the first scenario, which corresponds to the lower Rayleigh number values, the heat transfer is principally in form of diffusion, i.e., the diffusion term in the energy equation is predominant, leading to an almost vertical linear contouring for the temperature, and a recirculating motion in the core of the region is observed. As we increase the value of the Rayleigh number in $\mathcal{D}$, the flow is stretched to the walls, especially to the vertical walls; and the heat transfer starts to be driven mainly by convection. The isotherms become horizontal in a domain inside the cavity, far from the walls, that increases as the Rayleigh number increases. When we consider the second scenario, where the Rayleigh number range is higher, the velocity in the center of the cavity is practically zero, and presents large and normal gradients near the vertical walls. The temperature isolines are horizontal in a large domain inside the cavity, except near the vertical walls. This behavior agrees with the results presented in several works, e.g. [32, 33, 34]. In Fig. 2 we show the FE velocity magnitude and temperature for $\mu_{ph}=4363$, $\mu_{ph}=53778$ and $\mu_{ph}=667746$, with a fixed value of the geometrical parameter of $\mu_{g}=1$.

We consider different meshes depending the scenario. For $\mu_{ph}\in[10^{3},10^{5}]$ we consider a uniform mesh, with 50 divisions in each square side, i.e., $h=0.02\sqrt{2}$. For $\mu_{ph}\in[10^{5},10^{6}]$ we consider a finer mesh, with 70 divisions in each square side, i.e., $h=1/70\cdot\sqrt{2}$, in order to reproduce efficiently the eddies near the vertical walls appearing in this Rayleigh number range.

Concerning the time step in the evolution semi-implicit approach, we have considered a time step $\Delta t=0.01$ for the case of $\mu_{ph}\in[10^{3},10^{5}]$, and $\Delta t=2\cdot 10^{-3}$ for the case $\mu_{ph}\in[10^{5},10^{6}]$.

In the Reduced Basis framework, we perform an EIM for both the eddy viscosity and eddy diffusivity. Although for the numerical analysis performed in this work we have considered a regularized eddy diffusivity, the numerical tests are done with the eddy diffusivity defined in (2). Since the eddy diffusivity is proportional to the eddy viscosity, we only need to perform one EIM. With the EIM we are able to decouple the parameter dependence of the non-linear eddy viscosity and eddy diffusivity terms. For this test, we need $M=42$ basis until reaching a prescribed tolerance of $\varepsilon_{EIM}=5\cdot 10^{-3}$, when we consider that $\mu_{ph}\in[10^{3},10^{5}]$, and $M=150$ basis functions when we consider the second scenario where $\mu_{ph}\in[10^{5},10^{6}]$. In this last case, the Smagorinsky eddy viscosity and eddy diffusivity terms become more relevant, and for this reason, we take a lower tolerance for this test with respect to the previous one, considering $\varepsilon_{EIM}=10^{-4}$. In Fig. 3 we show the evolution of this error for both scenarios.

For the Greedy algorithm we prescribe a tolerance of $\varepsilon_{RB}=10^{-4}$ for both scenarios. For the first scenario, when $\mu_{ph}\in[10^{3},10^{5}]$, we need $N_{\max}=22$ basis to reach this tolerance. When $N=15$, holds the condition of Theorem 1 and $\tau_{N}(\boldsymbol{\mu})<1$ for all $\boldsymbol{\mu}$ in $\mathcal{D}$. In the second scenario, when $\mu_{ph}\in[10^{5},10^{6}]$, we need $N=N_{\max}=64$ basis functions to reach the tolerance previously prescribed, becoming $\tau_{N}(\boldsymbol{\mu})$ smaller than one when we get $N=52$ basis functions. In both cases, when $\tau_{N}(\boldsymbol{\mu})>1$ and the a posteriori error bound is not defined, we use as a posteriori error bound the proper $\tau_{N}(\boldsymbol{\mu})$. In Fig. 4 we show the convergence for the greedy algorithm.

In Fig. 5 (left) we show the comparison between the true error and the a posteriori error bound when $\mu_{ph}\in[10^{3},10^{5}]$, in which we can observe that the efficiency of the a posteriori error bound is between two and three orders of magnitude. Moreover, in Fig. 5 (right) we represent the comparison between the true error and the dual norm of the residual, $\epsilon_{N}(\boldsymbol{\mu})$. For this test, we can observe how the error correlates quite better with the dual norm of the residual than the a posteriori error bound, but in any case, it does not give us an upper bound of the error.

In Fig. 6 (left) we show the comparison between the true error and the a posteriori error bound when $\mu_{ph}\in[10^{5},10^{6}]$, in which we can observe that the efficiency of the a posteriori error bound is between one and two orders of magnitude. Moreover, in Fig. 6 (right) we represent the comparison between the true error and the dual norm of the residual, $\epsilon_{N}(\boldsymbol{\mu})$. For this test case, we can see how the dual norm of the residual is much lower than the error. Thus, the a posteriori error estimator gives us a better estimation for the error than the dual norm of the residual.

Finally, in Table 1, we show a comparison between the FE and RB solutions for several Rayleigh values in both scenarios. We show the computational time for solving a FE solution and a RB solution in the online phase. As can be observed, the speed-up rate of the computational time is larger than three orders of magnitude when $\mu_{g}\in[10^{3},10^{5}]$, while when $\mu_{g}\in[10^{5},10^{6}]$ the speed-up rate is close to three hundred. The difference in the speed-up magnitude between both cases is due to the longer number of EIM and RB functions computed in each case. In addition, we show the relative errors in $H^{1}$-norm for velocity and temperature, and in $L^{2}$-norm for pressure; for which we observe that the RB solution is close enough to the FE solution, with relative errors about $10^{-9}$ in both test cases. For this test, the offline phase when $\mu_{ph}\in[10^{3},10^{5}]$ took approximately 2 days in being performed. For the case when $\mu_{ph}\in[10^{5},10^{6}]$, the offline phase took approximately 3 weeks in be performed. In this offline computational time we consider either the EIM and the Greedy algorithm with the computation of the a posteriori error estimator

5.2 Geometrical parametrization

In this test, we consider a moderate Rayleigh number value $Ra=10^{5}$, and we consider the geometrical parameter ranging in $\mu_{g}\in\mathcal{D}=[0.5,2]$. The difference in the height of the cavity affects to the buoyancy force, making it more relevant when we increase the parameter value. This behavior is observed in Fig. 7, in which we show four solutions for different values of the geometrical parameter.

Firstly in the offline phase, we construct the reduced-basis space corresponding to the EIM, in which we approximate properly the eddy viscosity and eddy diffusivity terms. In this test, we need $M=73$ basis functions in order to reach a prescribed tolerance of $\varepsilon_{EIM}=10^{-4}$. In Fig. 8 (left) we show the evolution of the infinity norm of the error between the eddy viscosity $\nu_{T}(\mu_{g})$ and its EIM approximation.

For the Greedy algorithm, in this test, we prescribe a tolerance for the a posteriori error bound of $\varepsilon_{RB}=10^{-4}$. We need $N=23$ basis functions until to guarantee the condition of Theorem 1, and get $\tau_{N}(\boldsymbol{\mu})<1$. Then, we reach the prescribed tolerance when $N=N_{\max}=32$. In Fig. 8 (right), we show the maximum value for all $\mu_{g}\in\mathcal{D}$ of the a posteriori error estimator, and $\tau_{N}(\boldsymbol{\mu})$, in each iteration of the Greedy algorithm.

In Fig. 9 (left), we show a comparison between the *a posteriori * error bound and the true error for all $\mu_{g}\in\mathcal{D}$, in the last iteration of the Greedy algorithm, i.e., when $N=32$. Here we can see how the efficiency of the a posteriori error bound is about one order of magnitude. In Fig. 9 (right), we show the comparison between the a posteriori error bound and $\epsilon_{N}(\boldsymbol{\mu})$. In this case, we can observe how the dual norm of the residual is much lower than the error, and how the a posteriori error bound fits better the true error than the dual norm of the residual.

Finally, in Table 2, we summarize the results for several parameter values. We show the comparison between the time for computing a FE solution, and the online phase computational time. We obtain a speed-up rate of several hundreds in the computational time. The RB solution accuracy is fairly good, since the relative error is approximately of order $10^{-9}$ for velocity, temperature, and pressure. For this test, the offline phase took approximately 5 days in being performed. In this offline computational time we consider either the EIM and the Greedy algorithm with the computation of the a posteriori error estimator

5.3 Physical and geometrical parametrization

In this test, we perform a RB model in which a physical parameter (the Rayleigh number), and a geometric parameter are taken into account. Due to the increasing complexity in the flux with the consideration of this two parameters, we consider low range of Rayleigh number.

Thus, we consider that $\boldsymbol{\mu}=(\mu_{ph},\mu_{g})\in\mathcal{D}=[10^{3},10^{4}]\times[0.5,2]$. If we wanted to increase the Rayleigh number, we would have to consider a smaller interval for the geometric parameter. Indeed, as shown in sect. 5.1, the flow for high Rayleigh values is quite complex, thus the consideration of geometric parameter joint with the physical parameter is only possible if both intervals are not too big. If a big parameter set is required, a possible strategy is to split it in subsets of smaller amplitude.

For the EIM, in this test, we prescribe a tolerance of $\varepsilon_{EIM}=10^{-3}$. The error between $\nu_{T}(\textbf{u}_{h};\boldsymbol{\mu})$ and its interpolant fits this tolerance when $M=138$ basis functions are included in the EIM reduced-basis space. In Fig. 10 (left) we show the evolution of that error along the Greedy algorithm in the EIM.

For the Greedy algorithm in the offline phase we prescribe a tolerance of $\varepsilon_{RB}=10^{-3}$. This tolerance is reached when $N=N_{\max}=54$ basis functions are considered. We need $N=46$ basis functions to get $\tau_{N}(\boldsymbol{\mu})<1$, satisfying the conditions of Theorem 1, and having defined the a posteriori error bound $\Delta_{N}(\boldsymbol{\mu})$. In Fig. 10 (left) we show the evolution of the maximum value of $\tau_{N}(\boldsymbol{\mu})$ and $\Delta_{N}(\boldsymbol{\mu})$ in the Greedy algorithm. On the other hand, in Fig. 11 we show the value of the a posteriori error estimator, when $N=N_{\max}=54$, for all $\boldsymbol{\mu}\in\mathcal{D}$.

Finally, in Table 3 we sumarize some results obtained for some values of $\boldsymbol{\mu}\in\mathcal{D}$. There we show that the relative error between the FE solution and the RB solution is between order $10^{-7}$ and $10^{-9}$ for velocity and temperature, and between order $10^{-8}$ and $10^{-9}$ for pressure. For this test, the speedup rate obtained in the computation of the RB solution in the online phase with respect the computation of the FE solution is around fifty, due to the large number of EIM and RB basis functions needed, and the lower computational effort in the FE computation. Again, we have obtained a good accuracy in the RB solution with respect to the FE solution, with a considerable decrease of the computational time. For this test, the offline phase took approximately 2 weeks in being performed. In this offline computational time we consider either the EIM and the Greedy algorithm with the computation of the a posteriori error estimator.

6 Conclusions

In this work, we have developed a reduced turbulence model for buoyant flows in domains with geometrical variability. Specifically, we have dealt with the RB Boussinesq VMS-Smagorinsky model, for a variable height cavity. To represent this variability in the cavity height, we have parametrized the domain. Thus, we needed to reformulate our problem in a reference domain, which does not depend on the geometric parameter. As main technical tool, we have developed an a posteriori error estimator for the greedy algorithm involved in the reduced basis space construction. This construction is based upon the Brezzi-Rappaz-Raviart theory. We had to regularize the eddy viscosity for temperature, in order to ensure that the Boussinesq-Smagorinsky operator is locally Lipschitz-continuous.

Moreover, we have presented three different tests, considering geometrical parameters, physical parameters, or both. For each test, we obtained an accurate RB solution with a speedup rate going from one thousand in the simplest case, to fifty in the most complex case from one thousand for variability of only the physical parameter with diffusion-dominant effects, to nearly fifty for both geometrical and physical parameter variability.

Acknowledgements

This work has been supported by Spanish Government Project MTM2015-64577-C2-1-R and RTI2018-093521-B-C31, and COST Action TD1307

EU-MORNET. Francesco Ballarin and Gianluigi Rozza acknowledge project H2020 ERC CoG AROMA-CFD (GA 681447).