A Reduced-Order Shifted Boundary Method for Parametrized incompressible Navier-Stokes equations

TL;DR

This paper introduces a reduced-order model for steady incompressible Navier-Stokes equations that leverages the Shifted Boundary Method and proper orthogonal decomposition, enabling efficient simulations with complex geometries without remeshing.

Contribution

It presents a novel embedded reduced basis approach using the Shifted Boundary Method for parametrized geometries, avoiding remeshing and reference domain formulations.

Findings

Validated convergence and efficiency on 2D examples

Effective handling of geometrical parametrization

Reduced computational cost for multiple parameters

Abstract

We investigate a projection-based reduced-order model of the steady incompressible Navier-Stokes equations for moderate Reynolds numbers. In particular, we construct an "embedded" reduced basis space, by applying proper orthogonal decomposition to the Shifted Boundary Method, a high-fidelity embedded method recently developed. We focus on the geometrical parametrization through level-set geometries, using a fixed Cartesian background geometry and the associated mesh. This approach avoids both remeshing and the development of a reference domain formulation, as typically done in fitted mesh finite element formulations. Two-dimensional computational examples for one and three-parameter dimensions are presented to validate the convergence and the efficacy of the proposed approach.

Figures

Tables

Table 1. Table 1 . Relative error between the full-order solution and the reduced basis solution for velocity and pressure in the case of the geometrical parametrization with a one-dimensional parameter space. Results are reported for different dimensions of the reduced basis spaces.

Snapshots:	600
Modes	relative error
$N$	$𝒖$	$p$
20	0.0752959	19.590245
40	0.0193090	1.0666177
60	0.0128775	0.8244898
80	0.0064458	0.5823617
100	0.0039197	0.3994657
120	0.0026156	0.0731744
140	0.0024058	0.0545712

Table 2. Table 2 . (I) The relative error between full-order solution and reduced basis solution for velocity and pressure for 900 900 900 snapshots in the case with a geometrical parametrization using three-dimensional parameter space. (II) Execution time at the reduced order level. The computation time includes the assembling of the full-order matrices, their projection and the resolution of reduced problem. Times are for the resolution of one random value of the input parameter. The time execution at full-order level is equal to ≈ 50.70 absent 50.70 \approx 50.70 sec.

Snapshots:	900
Modes	(I) relative error		(II) execution time
$N$	$𝒖$	$p$	(sec)
20	0.0794779	2.4790682	8.5617990
40	0.0388933	0.5546928	8.9083982
60	0.0294700	0.4594517	8.8653870
80	0.0200465	0.3642104	8.8223756
100	0.0085595	0.0558062	8.7412506
120	0.0076601	0.0428119	8.7539559
140	0.0068780	0.0182534	8.8129689

Table 3. Table 3 . Supremizers “inf-sup” condition enrichment: Relative error between the full-order solution and the reduced basis solution for velocity and pressure in the case of the geometrical parametrization with a three-dimensional parameter space. Results are reported for different dimensions of the reduced basis spaces.

Snapshots:	900
Modes	relative error
$N$	$𝒖$	$p$
20	5.2407257	46.628631
40	0.0699833	0.3109475
60	0.0537041	0.2613186
80	0.0374247	0.2116895
100	0.0107190	0.0365011
120	0.0093516	0.0237161
140	0.0062227	0.0132555

Equations

$$\rho \nabla \cdot (\mathbf{u}(\mu )\otimes \mathbf{u}(\mu ))+\nabla p(\mu )-\nabla \cdot (2\nu \mathbf{ϵ}(\mathbf{u}(\mu )))-\rho \mathbf{g}(\mu )$$

$$\nabla \cdot \mathbf{u}(\mu )$$

$$\mathbf{u}(\mu )$$

$$-(\rho \mathbf{u}(\mu )\otimes \mathbf{u}(\mu ){\chi }_{{\Gamma }_N^{-}(\mu )}+p(\mu )\mathbf{I}-2\nu \mathbf{ϵ}(\mathbf{u}(\mu )))\mathbf{n}$$

$${\mathbf{V}}_h(\mathcal{D}(\mu ))=\left\{\mathbf{v}\in (C^0(\mathcal{D}(\mu )){)}^2:\mathbf{v}{|}_K\in (P^1(K){)}^d,\forall K\in {\mathcal{D}}_{\mathcal{T}}(\mu )\right\},$$

$$Q_{h}({\mathcal{D}}(\mu))=\left\{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varv}\in C^{0}({{\mathcal{D}}(\mu)}):{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varv}|_{K}\in{P}^{1}(K),\forall K\in{\mathcal{{{D}}_{T}}}(\mu)\right\},$$

$$0$$

$${\mathbf{M}}^h:$$

$$\tilde{\mathbf{x}}\mapsto \mathbf{x},$$

$${\mathbf{d}}_{{\mathbf{M}}^h}(\tilde{\mathbf{x}})=\mathbf{x}-\tilde{\mathbf{x}}=[{\mathbf{M}}^h-\mathbf{I}](\tilde{\mathbf{x}}).$$

$${\tilde{\Gamma }}_D=\{\tilde{E}\subseteq \tilde{\Gamma }:{\mathbf{M}}^h(\tilde{E})\subseteq {\Gamma }_D\}$$

$${\mathbf{V}}_h(\tilde{\mathcal{D}}(\mu ))=\left\{\mathbf{v}\in (C^0(\tilde{\mathcal{D}}(\mu )){)}^2:\mathbf{v}{|}_K\in (P^1(K){)}^d,\forall K\in {\tilde{\mathcal{D}}}_{\mathcal{T}}(\mu )\right\},$$

$$Q_h(\tilde{\mathcal{D}}(\mu ))=\left\{\text{varv}\in C^0(\tilde{\mathcal{D}}(\mu )):\text{varv}{|}_K\in P^1(K),\forall K\in {\tilde{\mathcal{D}}}_{\mathcal{T}}(\mu )\right\},$$

$$0$$

$$0={\mathbb{I}}_{\mathbb{NS}}[\nu ](\mathbf{u},p;\mathbf{w},q)+{\mathbb{II}}_{\mathbb{STAB}}[\nu ](\mathbf{u},p;\mathbf{w},q),$$

$$\begin{gathered} G(U(\mu ))U(\mu ):=G\left(\left[\begin{array}{c}\mathbf{u}(\mu )\\ p(\mu )\end{array}\right]\right)=\left[\begin{array}{c}\mathbf{A}\mathbf{u}(\mu )+\mathcal{C}(\mathbf{u}(\mu ))+{\mathbf{B}}^{T}p(\mu )\\ (\mathbf{B}+\hat{\mathbf{B}})\mathbf{u}(\mu )+\mathbf{C}p(\mu )\end{array}\right] \\ =\left[\begin{array}{cc}\mathbf{A}+\mathcal{C}(\mathbf{u}(\mu ))& {\mathbf{B}}^T\\ \mathbf{B}+\hat{\mathbf{B}}& \mathbf{C}\end{array}\right]\left[\begin{array}{c}\mathbf{u}(\mu )\\ p(\mu )\end{array}\right], \end{gathered}$$

$$F(\mu ):=\left[\begin{array}{c}{\mathbf{F}}_g(\mu )\\ {\mathbf{F}}_q(\mu )\end{array}\right],$$

$$R(U(\mu ))=G(U(\mu ))U(\mu )-F(U(\mu )).$$

$$\nabla G(U^{n-1}(\mu ))=\left[\begin{array}{cc}\frac{\partial (\mathbf{A}\mathbf{u}(\mu )+\mathcal{C}(\mathbf{u}(\mu )))}{\partial \mathbf{u}}& {\mathbf{B}}^T\\ \mathbf{B}& 0\end{array}\right],$$

$$\nabla G(U^{n-1}(\mu ))\delta U(\mu )=-R(U^{n-1}(\mu )).$$

$${\mathcal{S}}_{\mathbf{u}}=[\mathbf{u}({\mu }^1),\dots ,\mathbf{u}({\mu }^{N_s})]\in {\mathbb{R}}^{N_u^h\times N_s},{\mathcal{S}}_{\mathbf{p}}=[p({\mu }^1),\dots ,p({\mu }^{N_s})]\in {\mathbb{R}}^{N_p^h\times N_s}.$$

$$E_{N_{POD}}=\sum _{i=1}^{N_s}||{\mathbf{u}}_i-\sum _{k=1}^{N_{POD}}{a}_i^k{\mathbf{\phi }}_k|{|}_{L^2(\mathcal{D})}^2,$$

$${\mathcal{X}}^u{\mathbf{Q}}^u={\mathbf{Q}}^{\mathbf{u}}{\mathbf{\lambda }}^{\mathbf{u}},\text{mbox}for{\mathcal{X}}_{ij}^u=({\mathbf{u}}_i,{\mathbf{u}}_j{)}_{\mathcal{D}}\text{mbox},i,j=1,\dots ,N_s,$$

$${\mathbf{\phi }}_{\mathbf{i}}=\frac{1}{N_s{{\lambda }_{ii}^u}^{1/2}}\sum _{j=1}^{N_s}{\mathbf{u}}_j{Q}_{ij}^u.$$

$${\chi }_i=\frac{1}{N_s{{\lambda }_{ii}^p}^{1/2}}\sum _{j=1}^{N_s}{p}_j{Q}_{ij}^p.$$

$$\begin{array}{r}{\mathbf{L}}_u=[{\mathbf{\phi }}_1,\dots ,{\mathbf{\phi }}_{N_u^r}]\in {\mathbb{R}}^{N_u^h\times N_u^r},{\mathbf{L}}_p=[{\chi }_1,\dots ,{\chi }_{N_p^r}]\in {\mathbb{R}}^{N_p^h\times N_p^r}.\end{array}$$

$${\mathbf{u}}^{\mathbf{r}}\approx \sum _{i=1}^{N_u^r}{a}_i(\mu ){\mathbf{\phi }}_i(\mathbf{x})={\mathbf{L}}_u\mathbf{a}(\mu ),\text{mbox}{p}^r\approx \sum _{i=1}^{N_p^r}{b}_i(\mu ){\chi }_i(\mathbf{x})={\mathbf{L}}_p\mathbf{b}(\mu ).$$

$$\mathbf{L}\nabla G(U^{n-1}(\mu )){\mathbf{L}}^T\delta V(\mu )=-{\mathbf{L}}^{T}R(U^{n-1}(\mu )),$$

$$\nabla G^r(V^{n-1}(\mu ))\delta V(\mu )=-R^r(V^{n-1}(\mu )).$$

$$\underset{0\ne q\in Q_h}{\inf }\underset{0\ne \mathbf{v}\in {\mathbf{V}}_h}{\sup }\frac{⟨\nabla \cdot \mathbf{v},q⟩}{||\nabla \mathbf{v}||||q||}\ge \beta >0.$$

$${\tilde{\mathbf{L}}}_u=[{\mathbf{\phi }}_1,\dots ,{\mathbf{\phi }}_{N_u^r},{\mathbf{\eta }}_1,\dots ,{\mathbf{\eta }}_{N_{\text{sup}}^r}]\in {\mathbb{R}}^{N_u^h\times (N_u^r+N_{\text{sup}}^r)}.$$

$$\Delta {\mathbf{s}}_i=-\nabla {\chi }_i\text{mbox}in\mathcal{D}({\mu }^i),{\mathbf{s}}_i=0\text{mbox}on\Gamma ({\mu }^i)$$

Full text

1. Introduction

We consider a nonlinear system arising from computational fluid dynamics problems. In particular, a stationary Navier–Stokes system is examined and solved by the Shifted Boundary Method (SBM), an embedded boundary finite element method (EBM) that was recently proposed. The geometrical parametrization is described by means of level sets defined over a fixed (undeformed) background mesh. Moreover, a reduced order method (ROM), based on a proper orthogonal decomposition approach (POD-Galerkin) is tested, with the purpose of creating an appropriate embedded reduced order basis and decreasing the computational cost in the numerical solution procedure.

Starting with the pioneering work of Peskin [1], the scientific community has shown great interest on embedded or immersed methods. Some recent developments are represented by Ghost-Cell Finite Difference Methods [2], Cut-Cell Finite Volume Methods [3, 4], Immersed Interface Methods [5, 6], the Ghost Fluid Method [7, 8], the Volume Penalty Method [9, 10, 11], Cut-FEMs [12] and the Shifted Boundary Method (SBM) [13, 14, 15]. The interested reader can find additional details in [16, 17, 18, 19, 20, 21] and references therein. Later on, in the context of Navier–Stokes equations, Cut-FEMs [22, 23] and unfitted discontinuous Galerkin method [24] were employed. Recently, the authors of [13] introduced a new unfitted finite method called the Shifted Boundary Method (SBM) for the Poisson and Stokes problems. Their work circumvented the small cut-cell problem by defining a surrogate domain consisting of all uncut elements of the mesh that lie inside the computational domain. The definition of the surrogate domain, however, poses the challenge of devising effective strategies to the imposition of boundary conditions in order to maintain optimal convergence rates. To that end, weakly enforcing a Taylor extrapolation of the numerical solution at the surrogate boundary had the effect of correcting the boundary condition to account for the discrepancy between the surrogate and true boundary. Numerical results for the Poisson and Stokes problems showed optimal convergence rates and adequate conditioning of the matrix. Afterwards, the SBM was extended to the advection/diffusion equation and the laminar and turbulent incompressible Navier-Stokes equations  [14] in addition to hyperbolic systems such as wave propagation problems in acoustics and shallow water flows  [15]. The authors of [25] proposed a number of strategies to increase the accuracy of the method in elliptic problems and the authors of  [26] presented a complete numerical analysis of the SB approach for the Stokes problem. Important advances on the original method was presented in [27], which contained many simplifications and improvement to the original SBM formulation, and [28], which addresses numerical convergence and stability in the case of general domains with corners and edges. In [29], the SBM was extended to problems with internal interfaces.

Although immersed and embedded methods show improved features than fitted mesh methods in the case of geometrical design changes, there are still many cases where the approximate solutions of partial differential equations become computationally unaffordable, e.g. in real time problems, uncertainty or parametrization of the geometry etc. In these cases reduced order modeling techniques appear beneficial, [30, 31, 32, 33].

The main goal of this paper is to show how the SBM can solve geometrically parametrized nonlinear partial differential systems within the ROM framework. For this purpose, recently developed POD techniques [34, 35, 36, 37] will be applied. In these seminal works the methodology was developed on simple linear problems such as the Poisson and the Stokes setting. The main novelty of the current article is the extension of the methodology to more complex and nonlinear settings such as the Navier-Stokes one. Such an extension due to the non-linearity of the problem is not trivial. The ROM formulation had in fact to be reshaped and adapted to the underlying iterative approach used by the FOM solver. As it will be more clear in section 2 the non-linearity is resolved using a Newton’s method. Therefore, at each iteration, it is required the computation of the Jacobian of the residual vector and an incremental solution snapshot is obtained by the resolution of the resulting linear system of equations. The same iterative procedure has been implemented also at the ROM level during the online phase. The problem is therefore written in incremental form and the POD basis functions are computed starting from the incremental snapshots and not from the converged ones. The ROM is obtained by the projection of the Jacobian matrix and of the residual vector onto the POD incremental basis functions in order to obtain an incremental reduced solution. The key feature of our approach is the avoidance of a remeshing stage and/or morphing (i.e., a mapping of all the deformed geometries to a reference domain, see e.g. [30, 38, 39, 40, 41, 42, 33] for the use of this strategy in traditional body-fitted mesh finite element methods). This contribution is organized as follows:

1. (1)

In Section 2 we define the continuous strong formulation of the mathematical problem and the Nitsche weak formulation. A discrete Shifted Boundary weak formulation is also presented, for the full-order discretization, together with an incremental iterative scheme needed to solve the high fidelity problem during the offline stage. 2. (2)

We introduce the reduced order model formulation, the Proper Orthogonal Decomposition and its main features in Section 3. 3. (3)

In section 4, the proposed ROM-SBM technique is tested on a geometrically parametrized problem of the flow around an embedded rectangular domain, and convergence results, errors and execution times are also reported. 4. (4)

Finally, in Section 5, conclusions and perspectives for future improvements and developments are introduced.

2. The mathematical model and the full-order approximation

2.1. Strong formulation of the steady Navier-Stokes problem

The stationary Navier-Stokes equations for viscous incompressible flow describe the flow of a Newtonian, incompressible viscous fluid in a domain when convective forces are not negligible with respect to viscous forces. Consider an open domain $\mathcal{D}$ in ${\mathbb{R}}^{d}$, with $d=2,3$ the number of space dimensions with Lipschitz boundary $\Gamma$, decomposed into two sub-boundaries $\Gamma_{D}$, $\Gamma_{N}$. Let $\mathcal{P}\subset\mathbb{R}^{k}$ be a $k-$dimensional parameter space and $\mu\in\mathcal{P}$ a parameter vector.

The strong form of the stationary Navier-Stokes equations with Dirichlet and Neumann boundary conditions on $\Gamma_{D}$ and $\Gamma_{N}$, respectively is given by:

[TABLE]

where the variable $\mu$ introduces a geometrical parameterization. We denote by $\nu$ the viscosity, $\rho$ the density, ${\bm{\epsilon}({\bm{u}}{(\mu)})}=1/2\left(\nabla{{\bm{u}}(\mu)}+\nabla{{\bm{u}}({\mu})}^{T}\right)$ the velocity strain tensor (i.e., the symmetric gradient of the velocity ${\bm{u}}(\mu)$), $p(\mu)$ the pressure, ${\bm{g}}(\mu)$ a body force, ${\bm{g}_{D}}(\mu)$ the values of the velocity on the Dirichlet boundary and ${\bm{g}_{N}}(\mu)$ the normal stress on the Neumann boundary.

We will use for inflow and outflow boundaries the notation $\Gamma^{-}_{D}(\mu)=\{{\bm{x}}\in\Gamma_{D}(\mu)|{{\bm{g}}}_{D}(\mu)\cdot{\bm{n}}<0\}$, $\Gamma^{-}_{N}(\mu)=\{{\bm{x}}\in\Gamma_{N}(\mu)|{\bm{u}}(\mu)\cdot{\bm{n}}<0\}$, and \Gamma^{+}_{D}(\mu)={{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{D}(\mu)\setminus\Gamma^{-}_{D}(\mu)}}} and \Gamma^{+}_{N}(\mu)={{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{D}(\mu)\setminus\Gamma^{-}_{N}(\mu)}}}, where ${\mathbb{R}}\text{e}[\nu]=\rho UL/\nu$ is the Reynolds number, $U$ and $L$ are the characteristic speed and length of the problem, $h$ is the characteristic mesh size, and $\gamma_{N1}$ and $\gamma_{N2}$ are penalty parameters, with $\chi_{\Gamma^{-}_{D}}(\mu)$ and $\chi_{\Gamma^{-}_{N}}(\mu)$ the characteristic functions of the boundaries $\Gamma^{-}_{D}(\mu)$ and $\Gamma^{-}_{N}(\mu)$, respectively.

2.2. The conformal Nitsche’s weak formulation

For the sake of simplicity, in this subsection, we will omit the parameter dependency with respect to $\mu$. Let ${\bm{V}_{h}}({\mathcal{D}}(\mu))$ and $Q_{h}({\mathcal{D}}(\mu))$ be the spaces of continuous, piecewise-linear, vector- and scalar-valued functions. Namely:

[TABLE]

[TABLE]

where $\mathcal{D_{T}}(\mu)$ is the conformal mesh to the true geometry grid. We can now introduce a Nitsche’s [43, 44, 45] variational formulation of the Navier-Stokes equations:

*find ${\bm{u}}\in{\bm{V}_{h}}({\mathcal{D}}(\mu))$ and $p\in Q_{h}({\mathcal{D}}(\mu))$ such that, $\forall{\bm{w}}\in{\bm{V}_{h}}({\mathcal{D}}(\mu))$ and $\forall q\in Q_{h}({\mathcal{D}}(\mu))$, *

[TABLE]

Here we use the standard notation $(\cdot,\cdot)_{{\mathcal{D}}}$, $\langle\cdot,\cdot\rangle_{\Gamma_{D}}$ , $\langle\cdot,\cdot\rangle_{\Gamma_{N}}$ for the $L^{2}({\mathcal{D}})$, $L^{2}(\Gamma_{D})$ and $L^{2}(\Gamma_{N})$ inner products over ${\mathcal{D}}$, $\Gamma_{D}$ and $\Gamma_{N}$, respectively, and we denote the diameter of element $K\in{\mathcal{T}}$ as $h_{K}$. The overall size of the mesh is denoted by $h=\max_{K\in{\mathcal{T}}}h_{K}$.

2.3. Shifted boundary variational formulation

In this subsection, we introduce the Shifted Boundary Method for the Navier-Stokes equations [13, 14]. We define a surrogate computational domain $\tilde{\mathcal{D}}$ with surrogate boundary $\tilde{\Gamma}$ that approximate the true computational domain $\mathcal{D}$ and its boundary $\Gamma$ as seen in Figure 1 and Figure 2(a). Furthermore, $\bm{\tilde{n}}$ indicates the unit outward-pointing normal to the surrogate boundary $\tilde{\Gamma}$, which differs from the outward-pointing normal $\bm{n}$ to ${\Gamma}$ (see Figure 2(b)). $\tilde{\Gamma}$ is composed of the edges/faces of the mesh that are the closest to the true boundary $\Gamma$ in the sense of the closest-point projection, as shown in Figure 2(b).

In particular, we consider a mapping

[TABLE]

as defined in [28, Section 2.1], which associates to any point $\tilde{\bm{x}}\in\tilde{\Gamma}$ on the surrogate boundary a point $\bm{x}=\bm{M}^{h}(\tilde{\bm{x}})$ on the physical boundary $\Gamma$. Through $\bm{M}^{h}$, a distance vector function $\bm{d}_{\bm{M}^{h}}$ can be defined as

[TABLE]

For the sake of simplicity, we set $\bm{d}=\bm{d}_{\bm{M}^{h}}$ where $\bm{d}=\|\,\bm{d}\,\|\bm{\nu}$ and $\bm{\nu}$ is a unit vector.

For smooth surfaces with a single type of boundary condition, $\bm{M}^{h}$ is the closest point projection and therefore $\bm{\nu}=\bm{n}$ (see Figure 2(b)). On the other hand, when the boundary $\Gamma$ is partitioned into a Dirichlet boundary $\Gamma_{D}$ and a Neumann boundary $\Gamma_{N}$ with $\Gamma=\overline{\Gamma_{D}\cup\Gamma_{N}}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$ (see Figure 2(c)), we need to identify whether a surrogate edge $\tilde{E}\subset\tilde{\Gamma}$ is associated with $\Gamma_{D}$ or $\Gamma_{N}$. Through $\bm{M}^{h}$, we can partition $\tilde{\Gamma}$ as $\overline{\tilde{\Gamma}_{D}\cup\tilde{\Gamma}_{N}}$ with $\tilde{\Gamma}_{D}\cap\tilde{\Gamma}_{N}=\emptyset$ such that

[TABLE]

and $\tilde{\Gamma}_{N}=\tilde{\Gamma}\setminus\tilde{\Gamma}_{D}$. In Figure 2(c), for example, $\bm{M}^{h}$ associates the surrogate edge $\tilde{E}$ with $\Gamma_{D}$ whereby all the boundary information will be transported from $\Gamma_{D}$ to $\tilde{E}$. In turn, the distance $\bm{d}$ is given as the closest point projection from $\tilde{E}$ to $\Gamma_{D}$ which may not align with the true normal to $\Gamma_{D}$. For an extensive discussion on how to construct the distance $\bm{d}$ on domains with corners and edges in two and three dimensions, see [28].

Since the true surface is smooth between edges and corners, it is allowed to assume that ${\bm{M}^{h}}$ is continuous, and Lipschitz. Using ${\bm{M}^{h}}$, the unit normal vector $\bm{n}$ to the boundary $\Gamma$ can easily extend to the boundary $\tilde{\Gamma}$ as $\bm{\bar{}}{\bm{n}}(\tilde{\bm{x}})\equiv{\bm{n}}({\bm{M}^{h}}(\tilde{\bm{x}}))$. In the following, whenever we write $\bm{n}(\tilde{\bm{x}})$ we actually mean $\bm{\bar{}}{\bm{n}}(\tilde{\bm{x}})$ at a point ${\tilde{\bm{x}}}\in\tilde{\Gamma}$. The above constructions are the basis for the extension of the boundary conditions on $\Gamma$ to the boundary $\tilde{\Gamma}$ of the surrogate domain.

We can now introduce the SBM variational formulation [13, 14]. Consider the surrogate domain $\tilde{\mathcal{D}}$ and a discretization of the continuous boundary value problem with a mesh ${\mathcal{{\tilde{D}}_{T}}}$ consisting of simplexes $K$ belonging to a tessellation $\mathcal{T}$. Moreover, we introduce the discrete spaces ${\bm{V}}_{h}(\tilde{\mathcal{D}}(\mu))$ and $Q_{h}(\tilde{\mathcal{D}}(\mu))$, for velocity and pressure and we assume that a stable and convergent base formulation for the flow exists for these spaces in the case of conformal grids. Using the modified continuous piecewise-linear spaces

[TABLE]

[TABLE]

we can now introduce the Shifted boundary variational formulation:

• find ${\bm{u}}\in{\bm{V}}_{h}({\tilde{\mathcal{D}}(\mu)})$ and $p\in Q_{h}({\tilde{\mathcal{D}}(\mu)})$ such that, $\forall{\bm{w}}\in{\bm{V}}_{h}({\tilde{\mathcal{D}}(\mu)})$ and $\forall q\in Q_{h}({\tilde{\mathcal{D}}(\mu)})$,*

[TABLE]

where we indicated again by $(\cdot,\cdot)_{\tilde{\mathcal{D}}}$, $\langle\cdot,\cdot\rangle_{\tilde{\Gamma}_{D}}$ , $\langle\cdot,\cdot\rangle_{\tilde{\Gamma}_{N}}$ the $L^{2}(\tilde{\mathcal{D}})$, $L^{2}(\tilde{\Gamma}_{D})$ and $L^{2}(\tilde{\Gamma}_{N})$ inner products over $\tilde{\mathcal{D}}$, $\tilde{\Gamma}_{D}$ and $\tilde{\Gamma}_{N}$, respectively. For more details we refer to [14].

2.4. Variational multiscale stabilized finite element formulation

Because the proposed variational statement is not numerically stable, we introduce SUPG and PSPG stabilizing operators according to [46, 47], to which the reader can refer for more details. Since this is not the main focus of this paper, we refer the reader to [14, 48] for their implementation. The abstract variational form would then read:

*find ${\bm{u}}\in{\bm{V}}_{h}(\tilde{\mathcal{D}})$ and $p\in Q_{h}(\tilde{\mathcal{D}})$ such that, $\forall{\bm{w}}\in{\bm{V}}_{h}(\tilde{\mathcal{D}})$ and $\forall q\in Q_{h}(\tilde{\mathcal{D}})$, *

[TABLE]

where ${\mathbb{{II}_{STAB}}}$ $[\nu]({\bm{u}},p;{\bm{w}},q)$ denotes the SUPG/PSPG stabilization operators.

In what follows, it is useful to define the Navier-Stokes operator

[TABLE]

where $\bm{A}$ corresponds to the discrete diffusion operator, ${\bm{\mathcal{C}}}$ to the nonlinear convection operator. The rectangular matrix $\bm{B}^{T}$ represents the discrete gradient operator while $\bm{B}$ represents its adjoint, the divergence operator. $\bm{\hat{B}}$, $\bm{C}$ are instead associated with additional SUPG/PSPG stabilization operators. The right hand side

[TABLE]

is consisting of forcing and boundary data related to stabilization and Nitsche weak enforcement boundary terms. Using these definitions we can express the following residual of the algebraic system of equations:

[TABLE]

Furthermore, the Jacobian of $G(U(\mu))U(\mu)$ reads

[TABLE]

and yields the following iterative algebraic system of equations for the increment $\delta U(\mu)=U^{n}(\mu)-U^{n-1}(\mu)$:

[TABLE]

In the latter system of equations we underline some features that will play important role later in the ROM strategy described in § 3: the discretized differential operators $\bm{A}$, ${\bm{\mathcal{C}}}$, $\bm{B}$ and $\bm{\hat{B}}$ are parameter dependent and in the typical saddle point structure of the problem, the incompressibility equation is partially relaxed adding a stabilization term $\bm{C}$. This stabilization term, at full-order level, permits to the fulfillment of the “inf-sup” condition and the use of otherwise unstable pair of finite elements, such as $\mathbb{P}_{1}-\mathbb{P}_{1}$ it helps to preserve the stability of the reduced order model.

The presented formulation is used to solve the full-order problem during an offline stage and to produce the snapshots necessary for the construction of the ROM of § 3.

3. Reduced order model with a POD-Galerkin method

In this section, a POD-Galerkin approach will be analyzed as in [30, 49]. We will simplify the high fidelity model system to a reduced order one, which preserves its essential properties with the purpose of reducing computational cost in a way adapted to embedded-immersed boundary finite element methods. This approach allows flexibility with geometrical changes and to effectively and efficiently overcome several related issues that appear when using traditional FEM (see for instance [50, 34, 35, 36]).

Our interest is focused on the nonlinear system of Navier-Stokes equations with parametrized geometries and on the advantages of the SBM. We construct the ROM in a classical way starting from high dimensional SBM approximations (offline stage which involves the solution of a possibly large number of high fidelity problems). Reduced basis methods have been obtained starting from full-order approximations on fitted mesh methods for non-linear problems [51, 52, 53] – while for linear elliptic and linear parabolic equations we refer to [42, 54].

The first step in a reduced order method consists of obtaining a set of high fidelity solutions of the parametrized problem under input parameters variation. The aim of ROMs is to approximate any member of this solution set with a reduced number of basis functions. During the costly offline stage, one works out the solution set and examines its components in order to construct a reduced basis that approximates any member of the solution set to a prescribed accuracy. During a second stage, namely the online stage, after the Galerkin projection of the full-order differential operators describing the governing equations onto the reduced basis spaces, it is possible to solve a reduced problem for any new value of the input parameters with a reduced computational effort.

In general, POD-Galerkin ROMs for the incompressible Navier-Stokes equations are unstable [55, 56, 41], due to pressure instabilities, while for dynamic instabilities on transient problems the interested reader could see for instance [57, 58, 59, 60, 61]. Nevertheless, in the present approach, the SUPG and PSPG stabilization which is applied on the high fidelity solver is strongly propagating through the reduced basis construction procedure to the reduced level, and there is no need for further reduced basis stabilization and supremizers enrichment as in  [40, 38, 62, 63], see Appendix A.

The SBM unfitted/surrogate mesh Nitsche finite element method is used to apply parametrization and reduced order techniques considering Dirichlet combined with Neumann boundary conditions.

We highlight that a parametrized ROM method without the use of the transformation to reference domains will be used taking advantage of the fixed, geometrical parameter independent, background mesh [34, 36, 35].

3.1. The Proper Orthogonal Decomposition (POD)

For the projection of the high fidelity system to the reduced order one, there are several techniques. For more details about the different strategies the interested reader may see [42, 32, 64, 31, 65, 66]. In the present work, the POD is applied to the parameter-dependent full snapshots matrices. The full-order model is solved for each $\mu^{k}\in\mathcal{K}=\{\mu^{1},\dots,\mu^{N_{k}}\}\subset\mathcal{P}$ where $\mathcal{K}$ is a finite dimensional training set of parameters chosen inside the parameter space $\mathcal{P}$. The number of snapshots is denoted by $N_{s}$ and the number of degrees of freedom for the discrete full-order solution by $N_{u}^{h}$, $N_{p}^{h}$ for the velocity and pressure, respectively. The snapshots matrices $\bm{\mathcal{S}_{u}}$ and $\bm{\mathcal{S}_{p}}$, for velocity and pressure, are then given by $N_{s}$ full-order snapshots:

[TABLE]

Given a general scalar or vectorial function $\bm{u}:{\mathcal{D}}\to\mathbb{R}^{d}$, with a certain number of realizations $\bm{u}_{1},\dots,\bm{u}_{N_{s}}$, the POD problem consists in finding, for each value of the dimension of POD space $N_{POD}=1,\dots,N_{s}$, the scalar coefficients $a_{1}^{1},\dots,a_{1}^{N_{s}},\dots,a_{N_{s}}^{1},\dots,a_{N_{s}}^{N_{s}}$ and functions $\bm{\varphi}_{1},\dots,\bm{\varphi}_{N_{s}}$ that minimize the quantity:

[TABLE]

In this case the velocity field $\bm{u}$ is used as an example. It can be shown [67] that the minimization problem of equation (3.1) is equivalent to solving the following eigenvalue problem:

[TABLE]

where $\bm{\mathcal{X}^{u}}$ is the correlation matrix obtained from the parameter dependent snapshots $\bm{\mathcal{S}_{u}}$, $\bm{Q^{u}}$ is a square matrix of eigenvectors and $\bm{\lambda^{u}}$ is a diagonal matrix of eigenvalues.

The basis functions can then be obtained with:

[TABLE]

The same procedure is applied for the pressure field considering the snapshots matrix consisting of the snapshots $p_{1},p_{2},\dots,p_{N_{s}}$. The correlation matrix of the pressure field snapshots $\bm{C^{p}}$ is assembled and we solve a similar eigenvalue problem ${\bm{C}^{p}}{\bm{Q}^{p}}={\bm{Q}^{p}}{\bm{\lambda}^{p}}$. The POD modes ${\chi_{i}}$ for the pressure field can be computed as

[TABLE]

The POD spaces are constructed for both velocity and pressure using the above methodology resulting in the spaces:

[TABLE]

where $N_{u}^{r}$, $N_{p}^{r}<N_{s}$ are chosen according to the eigenvalue decay of $\bm{\lambda}^{u}$ and $\bm{\lambda}^{p}$, [42, 33].

Remark 3.1*.*

Following the ideas of [34], in the out-of-interest area (i.e., the region outside the true geometry), we prefer to use the solution values that are computed using the shifted boundary method and in particular the smooth mapping ${\mathcal{M}}$ from the true to the surrogate domain. This approach allows a smooth extension of the boundary solution to the neighboring ghost elements with values which are decreasing smoothly to zero, see for instance the zoomed image in Figure 4. This approach guarantees a regular solution in the background domain and permits therefore the construction of an effective reduced basis.

3.2. The projection stage and the generation of the ROM

Once the POD functional spaces are constructed, the reduced velocity and pressure fields can be approximated with:

[TABLE]

The reduced solution vectors $\bm{a}\in\mathbb{R}^{N^{r}_{u}\times 1}$ and $\bm{b}\in\mathbb{R}^{N^{r}_{p}\times 1}$ depend only on the parameter values and the basis functions $\bm{\varphi}_{i}$ and ${\chi}_{i}$ depend only on the physical space. Denoting by ${\bm{L}}=\begin{bmatrix}\bm{L}_{u}&\bm{0}\\ \bm{0}&\bm{L}_{p}\end{bmatrix}$ and by ${\bm{L}}^{T}=\begin{bmatrix}{\bm{L}}^{T}_{u}&{\bm{0}}\\ {\bm{0}}&{\bm{L}}^{T}_{p}\end{bmatrix}$, the unknown vectors of coefficients $V=\begin{bmatrix}\bm{a}\\ \bm{b}\end{bmatrix}$ then can be obtained through a Galerkin projection of the full-order system of equations onto the POD reduced basis spaces. The subsequent solution of a reduced iterative algebraic system of equations for the increment $\delta V(\mu)$ then becomes,

[TABLE]

which leads to the following algebraic reduced system:

[TABLE]

We remark here that at the reduced order level, we need to assemble the FOM problem in order to compute the reduced differential operator, but this expensive operation could be avoided, for example, using hyper reduction techniques [68, 69, 70, 71]. Moreover, during the online stage, also the stabilization term ${\bm{C}}$ (as well as the nonlinear term ${\bm{\mathcal{C}}}$) is projected onto the reduced basis space, which allows the inf-sup stabilization condition to propagate in an efficient way into the reduced model, as we will see in the following numerical examples.

Remark 3.2*.*

We remark here that the reduced basis spaces have been generated using the iterative solution snapshots and not only the final solutions. This is justified by the fact that also at the reduced order level an iterative procedure is used to solve the non-linear problem and to obtain the reduced basis solutions.

4. Numerical experiments

We consider two different test cases, based on the setup shown in Figure 3. The first one consists of a geometrical parametrization using a one-dimensional parameter space where both parameters $\mu_{0}$, $\mu_{1}$ are fixed. The second one consists of a geometrical parametrization with a three-dimensional parameter space where all the parameters $\mu_{0}$ $\mu_{1}$ and $\mu_{2}$ are left free. The problem domain is the rectangle $[-2,2]\times[-1,1]$, in which an embedded rectangular is immersed. The viscosity $\nu$ is set to $1$ and the Reynolds number is set to $100$. A constant velocity in the $x$ direction, $u_{\text{in}}=1$ is applied at the left side of the domain, and an open boundary condition with $\nabla\bm{u}\cdot\bm{n}=0$ on the right. In addition, a slip (no penetration) boundary condition is applied on the top and bottom edges while on the boundary of the embedded rectangle a no slip boundary condition is applied. The results for the test problems have been obtained with a mesh size of $h=0.0350$ for the background mesh, using $15022$ triangles for the discretization and $\mathbb{P}1/\mathbb{P}1$ polynomials.

4.1. Geometrical parametrization with one-dimensional parameter space

In this first experiment, the embedded domain consists of a rectangle of size $\mu_{1}\times\frac{\mu_{2}}{2}=2\times\frac{\mu_{2}}{2}$ where its aspect ratio inside the domain is parametrized with a geometrical parameter describing the size of the rectangle embedded domain with respect to its y-length with $\textit{aspect ratio}=\frac{4}{\mu_{2}}$ as in Figure 3. The horizontal x-length of the size of the box is not parametrized and the box’s center is located on the left bottom corner of the domain $[-2,-1]$. The ROM has been trained with $600$ samples and tested onto $10$ samples, for a parameter range $\mu_{2}\in[1,2]$ chosen randomly inside the parameter space. To test the accuracy of the ROM we compared its results on $10$ additional samples that were not used to create the ROM and were selected randomly within the same range. Under these considerations we record in Figure 6 the first six modes for the velocity magnitude and pressure, while in Table 1 the respective relative errors $||{\bf u}_{FOM}-{\bf{u}}_{ROM}||_{L^{2}({{\mathcal{D}}})}/||{\bf{u}}_{FOM}||_{L^{2}({{\mathcal{D}}})}$ and $||p_{FOM}-p_{ROM}||_{L^{2}({{\mathcal{D}}})}/||p_{FOM}||_{L^{2}({{\mathcal{D}}})}$ are summarized and they are depicted in Figure 9. In Figure 7 we depict the FOM and ROM solutions together with the relative error for both velocity magnitude and pressure, as well as the FOM and ROM streamlines for the velocity vector, for four different values of the input parameter.

4.2. A geometrical parametrization study with a three-dimensional parameter space

The second case considers a geometrical parametrization with a three-dimensional parameter space in the range $(\mu_{0},\mu_{1},\mu_{2})\in[-0.6,-0.5]\times[0.3,0.4]\times[1.3,1.4]$. We perform this test to examine the performances of the methodology on a more complex scenario where the box of the previous numerical test is parametrized onto the $x$ position of its center ($\mu_{0}$) and its aspect ratio ($\frac{\mu_{2}}{\mu_{1}}$) for both the $x$ and $y$ direction, namely width and height. The results are reported in Figure 6 (the first six modes for the velocity magnitude and pressure), in Table 2 I ($||{\bf u}_{FOM}-{\bf{u}}_{ROM}||_{L^{2}({{\mathcal{D}}})}/||{\bf{u}}_{FOM}||_{L^{2}({{\mathcal{D}}})}$ and $||p_{FOM}-p_{ROM}||_{L^{2}({{\mathcal{D}}})}/||p_{FOM}||_{L^{2}({{\mathcal{D}}})}$ relative errors report) and they are illustrated in Figure 10 I. In Figure 8 we visualize the FOM and ROM solutions together with the relative errors for both velocity magnitude and pressure, as well as the FOM and ROM streamlines for the velocity, for four different values of the input parameters. Finally in Table 2 II and  Figure 10 II the time savings are reported and visualized.

4.3. Some Comments

Additional testing using more snapshots and with and without supremizers yielded very similar results, in which supremizers delivered slightly worse velocity results and slightly improved pressure results. We found that the supremizers did not substantially improved the errors. In our opinion, this is a consequence of the incremental iterative formulation in the offline solver, which preserves the effects of the SUPG and PSPG stabilization in the reduced model, which allows an inf-sup stable reduced basis. Some basics related to these stability issues, and numerical results can be found in Appendix A. Both the full-order and the ROM simulation were run in serial on an Intel® CoreTM i7-4770HQ 3.70GHz CPU.

5. Conclusions and future developments

We have introduced a geometrically parametrized ROM model emulator of the two-dimensional Navier–Stokes equations, of much reduced computational cost. The ROM model evaluation through numerical tests shows good convergence properties and low errors. Comparing with high-fidelity solutions, numerical ROM errors improve with the increase of the size of training data and of the number of basis components. As future perspectives, we indicate applications to the time dependent Navier–Stokes systems, fluid structure interaction problems, and shallow water flows. Additionally, from the model reduction point of view, we will pursue further developments in hyper-reduction techniques [69, 70, 71, 68] and transportation methodologies [34].

Acknowledgments

This research has been supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research under Early Career Research Program Grant SC0012169 and the Army Research Office (ARO) under Grant W911NF-18-1-0308, Hellenic Foundation for Research and Innovation (HFRI) and the General Secretariat for Research and Technology (GSRT), under grant agreement No 1115, the National Infrastructures for Research and Technology S.A. (GRNET S.A.) under project ID pa190902, the European Research Council Executive Agency by means of the H2020 ERC Consolidator Grant project AROMA-CFD “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” - GA 681447, (PI: Prof. G. Rozza), INdAM-GNCS-2019, and by project FSE - European Social Fund - HEaD ”Higher Education and Development” SISSA operazione 1, Regione Autonoma Friuli - Venezia Giulia.

Appendix A

In this appendix we explore the experiment as described in Section 4.2. We will justify why in all previous experiments we did not use any additional RB stabilization, verifying numerically that the SUPG and PSPG stabilization which is applied on the high fidelity solver is strongly propagating through the reduced basis construction procedure to the reduced level. In contrast, we refer the interested reader to the classical works of [40, 38, 62, 63] where the reduced solution stabilization is deemed necessary. Subsequently, we introduce the basics related to the stability issues that could appear in the offline and online stage. In the full order method stage it is well known that the spaces have to satisfy the, also parametric in our case, Ladyzhenskaya-Brezzi-Babuska “inf-sup” condition see e.g. [72]. In particular, it is required that there should exist a constant $\beta>0$, independent to the discretization parameter $h$, such that

[TABLE]

In the present work this condition is fulfilled for the high fidelity solution through the SUPG and PSPG stabilization. Even if, the offline stage and the snapshots are realized and computed by a stabilized numerical method, it is not guaranteed that this stability is preserved onto the reduced basis spaces [38, 56, 40]. Next we briefly introduce the “inf-sup” condition enforcement in the reduced level using supremizers, we illustrate the relative errors results and we are comparing them with the case without supremizers enrichment. Within this approach, the velocity supremizer basis functions $\bm{L}_{\text{sup}}=[\bm{\eta}_{1},\dots,\bm{\eta}_{N_{\text{sup}}^{r}}]\in\mathbb{R}^{N_{u}^{h}\times N_{\text{sup}}^{r}}$ are constructed and added to the reduced velocity space (see Section 3.1) which is finally transformed into

[TABLE]

To obtain the latter enrichment for each pressure basis function $\chi_{i}$ the auxiliary “supremizer” problem:

[TABLE]

is solved with an SBM Poisson solver starting from the parameter value $\mu^{i}$. For each pressure basis function the corresponding supremizer element can be found and the solution ${\bm{s}_{i}}$ permits the realization of the “inf-sup” condition.

We emphasize that the above supremizer basis functions do not depend on the particular pressure basis functions and on the geometrical parameters, they are computed during the offline phase, and their calculation is based on the pressure snapshots.

Obviously, if someone compares the Table 3 with Table 2 and examines their visualization in Figure 11, supremizers drove the reduced solution to slightly worse velocity results and slightly improved pressure results. So, the supremizers did not substantially improved the errors and this is the reason that we avoided their application. In our opinion, this phenomenon is a consequence of the incremental iterative formulation in the offline solver, which preserves the effects of the SUPG and PSPG stabilization in the reduced model and allows an inf-sup stable reduced basis.