Projection-based reduced order models for a cut finite element method in parametrized domains

TL;DR

This paper introduces a reduced order modeling approach based on CutFEM that efficiently handles large domain deformations without remeshing, enabling fast evaluations for parametrized elliptic and Stokes problems.

Contribution

It combines embedded finite element methods with reduced order models to improve computational efficiency in parametrized domain problems, avoiding remeshing and reference domain complexities.

Findings

Effective in linear elliptic problems

Applicable to Stokes flow simulations

Reduces computational costs significantly

Abstract

This work presents a reduced order modelling technique built on a high fidelity embedded mesh finite element method. Such methods, and in particular the CutFEM method, are attractive in the generation of projection-based reduced order models thanks to their capabilities to seamlessly handle large deformations of parametrized domains. The combination of embedded methods and reduced order models allows us to obtain fast evaluation of parametrized problems, avoiding remeshing as well as the reference domain formulation, often used in the reduced order modelling for boundary fitted finite element formulations. The resulting novel methodology is presented on linear elliptic and Stokes problems, together with several test cases to assess its capability. The role of a proper extension and transport of embedded solutions to a common background is analyzed in detail.

Equations

$$a_{\gamma }(u(\mu ),v;\mu )={\ell }_{\gamma }(v;\mu )\text{ in }V(\mu ),$$

$$\{\begin{array}{ll}-\Delta u(\mu )=g(\mu ),& \text{in D(μ)},\\ u(\mu )=g_D(\mu ),& \text{on ΓD(μ)},\\ {\mathbf{n}}_{\mathbf{\Gamma }}\cdot \nabla u(\mu )=g_N(\mu ),& \text{on ΓN(μ)},\end{array}$$

$$\left(\nabla u(\mu ),\nabla \upsilon (\mu )\right)=\left(g(\mu ),\upsilon (\mu )\right)+{\left(g_N(\mu ),\upsilon (\mu )\right)}_{{\Gamma }_N(\mu )},\forall v(\mu )\in V_0(\mu ).$$

$$V_h(\mu )=\left\{\upsilon \in C^0({\overline{\mathcal{D}}}_{\mathcal{T}}(\mu )):\upsilon {|}_K\in P^1(K),\forall K\in {\mathcal{T}}_h(\mu )\right\},$$

$${\left(\nabla u_h(\mu ),\nabla {\upsilon }_h(\mu )\right)}_{\mathcal{D}(\mu )}-{\left({\mathbf{n}}_{\Gamma }\cdot \nabla u_h(\mu ),{\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}-{\left(u_h(\mu ),{\mathbf{n}}_{\Gamma }\cdot \nabla {\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}$$

$$+{\left({\gamma }_D{h}^{-1}{u}_h(\mu ),{\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}+{\left({\gamma }_{N}h{\mathbf{n}}_{\Gamma }\cdot \nabla u_h(\mu ),{\mathbf{n}}_{\Gamma }\cdot \nabla {\upsilon }_h(\mu )\right)}_{{\Gamma }_N(\mu )}+j(u_h(\mu ),{\upsilon }_h(\mu );\mu )$$

$$\displaystyle\qquad=\left(g{(\mu)},\upsilon_{h}(\mu)\right)_{{\mathcal{D}}{(\mu)}}+\left(g_{D}(\mu),\gamma_{D}h^{-1}\upsilon_{h}(\mu)-{\bf n}_{\Gamma}\cdot\nabla\upsilon_{h}(\mu)\right)_{\Gamma_{D}{(\mu)}}+\big{(}g_{N}(\mu),\upsilon_{h}(\mu)$$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\gamma_{N}h{\bf n}_{\Gamma}\cdot\nabla\upsilon_{h}(\mu)\big{)}_{\Gamma_{N}{(\mu)}},$$

$$\text{ where}j(u_h(\mu ),{\upsilon }_h(\mu );\mu )=\sum _{F\in {\mathcal{F}}_G(\mu )}{\left({\gamma }_{1}h[[{\mathbf{n}}_F\cdot \nabla u_h(\mu )]],[[{\mathbf{n}}_F\cdot \nabla {\upsilon }_h(\mu )]]\right)}_F,$$

$$b_{\gamma }(u_h(\mu ),{\upsilon }_h(\mu );\mu ):=-{\left(u_h(\mu ),{\mathbf{n}}_{\Gamma }\cdot \nabla {\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}+{\left({\gamma }_D{h}^{-1}{u}_h(\mu ),{\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}$$

$$+{\left({\gamma }_{N}h{\mathbf{n}}_{\Gamma }\cdot \nabla u_h(\mu ),{\mathbf{n}}_{\Gamma }\cdot \nabla {\upsilon }_h(\mu )\right)}_{{\Gamma }_N(\mu )},$$

$$a_{\gamma }(u_h(\mu ),{\upsilon }_h(\mu );\mu )={\ell }_{\gamma }({\upsilon }_h(\mu );\mu ),\forall {\upsilon }_h(\mu )\in V_h(\mu ),$$

$$a_{\gamma }(u_h(\mu ),{\upsilon }_h(\mu );\mu ):=a(u_h(\mu ),{\upsilon }_h(\mu );\mu )+b_{\gamma }(u_h(\mu ),{\upsilon }_h(\mu );\mu )+j(u_h(\mu ),{\upsilon }_h(\mu );\mu ),$$

$${\ell }_{\gamma }({\upsilon }_h(\mu );\mu )$$

$$a(u_h(\mu ),{\upsilon }_h(\mu );\mu )={\left(\nabla u_h(\mu ),\nabla {\upsilon }_h(\mu )\right)}_{\mathcal{D}(\mu )}-{\left({\mathbf{n}}_{\Gamma }\cdot \nabla u_h(\mu ),{\upsilon }_h(\mu )\right)}_{{\Gamma }_D(\mu )}.$$

$$\{\begin{array}{ll}-\nabla \cdot (2\nu \epsilon (u(\mu ))-p(\mu )I)=g(\mu ),& \text{ in }\mathcal{D}(\mu ),\\ \nabla \cdot u(\mu )=0,& \text{ in }\mathcal{D}(\mu ),\\ u(\mu )=g_D(\mu ),& \text{ on }{\Gamma }_D(\mu ),\\ (2\nu \epsilon (u(\mu ))-p(\mu )I)\cdot {\mathbf{n}}_{\Gamma }=0,& \text{ on }{\Gamma }_N(\mu ),\end{array}$$

$$V_h(\mu )$$

$$Q_h(\mu )$$

$$\nu (\nabla u,\nabla \psi )-(p,\nabla \cdot \psi )+(\nabla \cdot u,\xi )-(\nu n_{\Gamma }\cdot u-pn,\psi {)}_{{\Gamma }_D}-(u,\nu n_{\Gamma }\cdot \psi +\xi n,\psi {)}_{{\Gamma }_D}$$

$$+\sum _{K\in G_h}{\int }_{{\Gamma }_K}\nu {\gamma }_D{h}_K^{-1}u\psi =(g,\psi )-(g_D,\nu n_{\Gamma }\cdot \nabla \psi +\xi n{)}_{{\Gamma }_D}+\sum _{K\in G_h}{\int }_{{\Gamma }_K}\nu {\gamma }_D{h}_K^{-1}{g}_D\psi ,$$

$$b_{\gamma }(\mathbf{u},\mathbf{\varphi })$$

$$a(\mathbf{u},\mathbf{\varphi })$$

$${\ell }_{\gamma }(\mathbf{\varphi })$$

$$a_{\gamma }(\mathbf{u}(\mu ),\mathbf{\varphi }(\mu );\mu )={\ell }_{\gamma }(\mathbf{\varphi }(\mu );\mu ),\forall \mathbf{\varphi }(\mu )=(\psi (\mu ),\xi (\mu ))\in {\mathbf{V}}_h(\mu ).$$

$$a_{\gamma }(\mathbf{u},\mathbf{\varphi }):=a(\mathbf{u},\mathbf{\varphi })+b_{\gamma }(\mathbf{u},\mathbf{\varphi })+j_u(u,\psi ;\mu )-j_p(p,\xi ;\mu ),$$

$$j_u(u,\psi ;\mu )$$

$$j_p(p,\xi ;\mu )$$

$$\mathcal{S}=[\hat{u}({\mu }^1),\dots ,\hat{u}({\mu }^M)]\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}.$$

$${\hat{V}}_h=\left\{\upsilon \in C^0(\overline{\mathcal{B}}):\upsilon {|}_K\in P^1(K),\forall K\in {\mathcal{B}}_h\right\},$$

$$\mathcal{C}Q=Q\Lambda ,\text{being }{\mathcal{C}}_{ij}=⟨\hat{u}({\mu }^i),\hat{u}({\mu }^j){⟩}_{L^2({\mathcal{T}}_h)},\text{mbox}fori,j=1,\dots ,M,$$

$${\hat{V}}_N=\text{span}\{{\phi }_1,\dots ,{\phi }_N\}$$

$$u_N(\mu )=\sum _{i=1}^N{\alpha }_i(\mu ){\phi }_i,$$

$$\sum _{i=1}^N{\alpha }_i(\mu )a_{\gamma }({\phi }_i,{\phi }_j;\mu )={\ell }_{\gamma }({\phi }_j;\mu ),\forall j=1,\dots ,N.$$

$${\mathcal{S}}_u=[\hat{u}({\mu }^1),\dots ,\hat{u}({\mu }^{N_s})]\in {\mathbb{R}}^{{\mathcal{N}}_h^u\times M},{\mathcal{S}}_p=[\hat{p}({\mu }^1),\dots ,\hat{p}({\mu }^{N_s})]\in {\mathbb{R}}^{{\mathcal{N}}_h^p\times M},$$

$${\hat{V}}_h$$

Full text

*Corresponding author

1. Introduction and Motivation

A wide variety of numerical methods and computational libraries for the solution of problems governed by partial differential equations is nowadays available. However, there are still many cases in which either the solution of the governing equations or the solution of associated inverse problems becomes impractical or unfeasible using standard discretization techniques, such as the Finite Element Method (FEM). For instance, complicated topology of the problem or a complex geometry may pose a challenge in the discretization of complex phenomena, and ultimately affect the quality of the resulting simulation. Moreover, repeated queries to the underlying solver in the context of an iterative solution of inverse problems may result in unbearably large computational times. Such situations occur, for example, when a large number of different configurations are in need of being tested, such as in uncertainty quantification, optimal control and shape optimization. The overall objective of this manuscript is to investigate how the recently introduced unfitted mesh finite element methods may be used combined with reduced order modeling techniques for parametrized partial differential problems. Indeed, unfitted methods are very useful in cases characterized by complex geometrical configurations; however, as their Finite Element counterparts, they usually require large computational efforts. A combination with reduced order methods, able to widely decrease the overall computational time, would result in a very compelling methodology to be possibly applied in several different fields.

Classical embedded/immersed methods provide simple, efficient, and robust numerical schemes for solving PDE in general domains [1, 2, 3]. Since these early works, several improvements have been made, for instance for what concerns the rate of spatial accuracy near embedded boundaries. Recent improvements go under the names of Ghost-Cell finite difference methods, Cut-Cell finite volume approach, Immersed Interface, Ghost Fluid, Volume Penalty methods, for which we refer to the review paper [1] and references within. In particular, for what concerns incompressible flows in arbitrary smooth domains, the Immersed Boundary Smooth Extension method has shown high-order convergence for the incompressible Navier-Stokes equations [4].

More in detail, extended mesh finite element methods using cut elements are examined in [5, 6] for stationary Stokes flow systems, as well as for Navier-Stokes. An analysis for high Reynolds numbers, independent of the local Reynolds, has been carried out in [7, 8, 9]. XFEM approaches in 2D and 3D Navier-Stokes are reported in [10]. Higher Reynolds number aerodynamics problems in unbounded domains, thin vortex ring and Lattice Green function Immersed Boundary methods are studied in [11, 12]. Furthermore, embedded and immersed methods have been used in solving fluid structure interaction problems, see e.g. [13, 14, 15, 16, 17, 18, 19, 3, 20, 21, 22, 23]. More recently a new embedded boundary method called Shifted Boundary method was introduced in [24], with good approximation properties and experimental results in several fields, including Navier-Stokes systems and fluid structure interaction applications [24, 25, 26]. The main idea of this method is that boundary conditions are imposed on the boundary of a surrogate domain whose geometry is properly chosen avoiding cut cells.

First attempts to apply Reduced Order Methods (ROMs) in the context of viscous flows and Stokes / Navier-Stokes systems can be found in [27, 28]. ROMs based on FEM full order approximations have been used to treat several problems based on linear elliptic equations [29], linear parabolic equations [30] and even non-linear problems [31, 32]. Stability of the resulting reduced order systems is often an issue, see for more details [33, 34, 35, 36, 37, 38], and for transient systems see for instance [39, 40, 41, 42, 43].

The main novelty of this work is the combination of the cut element finite element method (CutFEM), which we will consider as high fidelity method, to projection-based model reduction techniques. In particular, this allows to overcome a reference domain formulation customarily employed in ROMs for parametrized geometries built on conforming discretizations, see e.g. [44, 45, 46, 47, 48, 49, 36, 29, 37] and references therein. Even though such reference domain formulation avoids remeshing when updating the parametric domain, the choice of the transformation map is usually problem-dependent [50, 44, 47], suitable parameter space dimensionality reduction techniques [51, 52, 53, 54] must be employed to identify the most relevant shape parameters and preserve good quality meshes, and often limited only to small parametric deformations. In contrast, the combination of ROMs with a CutFEM approach results in a novel methodology capable of handling large parametric deformations as well.

To the best of the authors’ knowledge this combination is not particularly investigated. We mention here the approach proposed in [55] for classical embedded methods and model reduction, in which two regions are separated by an evolving in time interface, where a snapshot compression problem as a weighted low-rank approximation is formulated. We also mention our on-going work on combination of the Shifted Boundary method and model reduction techniques for small parametric variations [56, 57, 58, 59]. The present work touches several key points of growing interest in the model reduction community, especially since it shares challenges to be tackled for the efficient reduction of advection dominated problems. In particular, a snapshots preprocessing is advocated by [60, 61, 62, 63, 64, 65] in order to improve the efficiency of the resulting reduced model for advection dominated problems by providing a better (smaller) representation of the reduced basis space. Such snapshots preprocessing procedure, based in our case on extension and transportation, will be pivotal in this work in order to deal with large deformations.

The work is organized as it follows: in Section 2 the abstract formulation of the problem is introduced, as well as the embedded method used for the high fidelity problem of Darcy flow pressure model and of the steady Stokes equations. The reduced order methodology is discussed in details in Section 3, whilst in Section 4 the proposed ROM technique is tested on several numerical benchmarks. Finally in Section 5 conclusions and perspectives are drawn, highlighting the directives for future improvements and developments.

2. High fidelity CutFEM approximation

We introduce in this section the abstract parametrized model problem, which has the form: find $u(\mu)\in V(\mu)$ such that

[TABLE]

where $a_{\gamma}(\cdot,\cdot;\mu)$ is the weak form of an operator defined on a domain $\mathcal{D}(\mu)\subset\mathbb{R}^{d}$, for $d=2,3$, while $\ell_{\gamma}(\cdot,\mu)$ is the right hand side of the system of equations related to the forcing term. The forms $a_{\gamma}(\cdot,\cdot;\mu)$, $\ell_{\gamma}(\cdot,\mu)$, the domain $\mathcal{D}(\mu)$, as well as the Sobolev space $V(\mu)$, depend on the parameter $\mu\in\mathcal{K}$, being $\mathcal{K}\subset\mathbb{R}^{K}$ the set of possible outcomes, which we assume to be a compact set in $\mathbb{R}^{K}$, $K\in\mathbb{N}$. Since CutFEM discretizations usually require suitable penalty and stabilization procedures, for the sake of exposition we will denote by $a_{\gamma}(\cdot,\cdot;\mu)$ and $\ell_{\gamma}(\cdot,\mu)$ the forms with penalty and stabilization, where $\gamma$ represents the penalty and stabilization coefficients, see e.g. [5] and references therein.

We start with a sketch description of the continuous strong form and the weak formulation used for the problems under consideration. The CutFEM formulation (2.1) will be used for the high fidelity simulation employed during the training of the reduced order model, which will be introduced in Section 3.

2.1. The Darcy flow pressure model

For any $\mu\in\mathcal{K}$, let $\mathcal{D}(\mu)\subset{\mathbb{R}}^{d}$ be a bounded domain depending on $\mu$, with boundary $\Gamma(\mu)$. We consider the following model problem in $\mathcal{D}(\mu)$: for any $\mu\in\mathcal{K}$, find $u(\mu):{\overline{\mathcal{D}}(\mu)}\to\mathbb{R}$ such that

[TABLE]

and $\Gamma_{D}(\mu)$, $\Gamma_{N}(\mu)$ are non-overlapping parts of $\Gamma(\mu)=\Gamma_{D}(\mu)\cup\Gamma_{N}(\mu)$ where Dirichlet and Neumann boundary conditions are applied, and $g(\mu)$, $g_{N}(\mu)$, $g_{D}(\mu)$ are given functions in $\mathcal{D}(\mu)$ and on the boundaries $\Gamma_{N}(\mu)$, $\Gamma_{D}(\mu)$ respectively, [66].

We start with the classical weak formulation in the spatial domain: for any $\mu\in\mathcal{K}$, find $u(\mu)\in V_{g_{D}}(\mu)=\left\{w\in H^{1}\left(\mathcal{D}(\mu)\right)\text{ with }w|_{\Gamma_{D}(\mu)}=g_{D}(\mu)\right\}$ such that

[TABLE]

We denote by $\mathcal{B}$ the background domain, and by $\mathcal{B}_{h}$ its corresponding mesh. The continuous boundary value problem is next formulated on a domain ${\mathcal{D}}_{\mathcal{T}}(\mu)$ that contains ${\mathcal{D}}(\mu)\subset{\mathcal{D}}_{\mathcal{T}}(\mu)$, while its mesh $\mathcal{T}_{h}(\mu)$ is not fitted to the domain boundary and ${\mathcal{D}}_{\mathcal{T}}(\mu)\subset\mathcal{B}$ and $\mathcal{T}_{h}(\mu)\subset\mathcal{B}_{h}$ for all $\mu\in\mathcal{K}$. Let also ${G}_{h}(\mu):=\{K\in\mathcal{T}_{h}(\mu):K\cap\Gamma(\mu)\neq\emptyset\}$ be the set of elements that are intersected by the interface. We remark that $\mathcal{T}_{h}(\mu)$, ${G}_{h}(\mu)$ and ${\mathcal{D}}_{\mathcal{T}}(\mu)$ depend on $\mu$ through $\mathcal{D}(\mu)$ (or its boundary), while the background domain $\mathcal{B}$ and its mesh $\mathcal{B}_{h}$ do not depend on $\mu$.

Furthermore, the set of element faces $\mathcal{F}_{G}(\mu)$ associated with ${G}_{h}(\mu)$, is defined as follows: for each face $F\in\mathcal{F}_{G}(\mu)$, there exist two simplices $K\neq K^{\prime}$ such that $F=K\cap K^{\prime}$ and at least one of the two is a member of ${G}_{h}(\mu)$. Note that boundary faces of $\mathcal{T}_{h}(\mu)$ are excluded from $\mathcal{F}_{G}(\mu)$. On a face $F\in\mathcal{F}_{G}(\mu)$, $F=K\cap K^{\prime}$, the jump of the gradient of $v\in C^{0}(\overline{\mathcal{D}}_{\mathcal{T}})$ is defined by $\mathopen{[\![}{{\bf n}_{F}}\cdot\nabla v\mathclose{]\!]}={\bf n}_{F}\cdot\nabla v|_{K}-{\bf n}_{F}\cdot\nabla v|_{K^{\prime}}$, where ${\bf n}_{F}$ denotes the outward pointing unit normal vector to $F$. Finally, let $h_{k}$ be the diameter of $K$, $h=\max_{K\in\mathcal{T}_{h}(\mu)}h_{k}$.

The CutFEM discretization is as follows. We seek a discrete solution $u_{h}(\mu)$ in the finite element space

[TABLE]

such that

[TABLE]

[TABLE]

for all test functions $\upsilon_{h}(\mu)\in V_{h}(\mu)$ and for \gamma_{D}\color[rgb]{0,0,0}{>0}, \gamma_{N}\color[rgb]{0,0,0}{\geq 0}, and \gamma_{1}\color[rgb]{0,0,0}{>0} penalty parameters, see for instance [67]. The stabilization term $j(u_{h}(\mu),\upsilon_{h}(\mu);\mu)$, which depends on $\gamma_{1}$, extends the coercivity from the physical domain ${\mathcal{D}}(\mu)$ to the extended mesh domain ${\mathcal{D}}_{\mathcal{T}}(\mu)$, while the penalty terms involving coefficients $\gamma_{D}$ and $\gamma_{N}$ account for a Nitsche weak imposition of boundary conditions. Assuming that for the Nitsche terms we use the notation

[TABLE]

the variational form (2.1), using the classical linear and bilinear forms, can be equivalently expressed as: find $u_{h}(\mu)\in V_{h}(\mu)$ such that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

We refer to [67] for more details concerning the CutFEM discretization of elliptic problems.

2.2. Steady Stokes problem

The strong form of the stationary Stokes flow system of equations with Dirichlet and Neumann boundary conditions, geometrically parametrized by $\mu$, is given by: for any $\mu\in\mathcal{K}$, find velocity $u(\mu):{\overline{\mathcal{D}}(\mu)}\to\mathbb{R}^{d}$ and pressure $p(\mu):{\overline{\mathcal{D}}(\mu)}\to\mathbb{R}$ such that

[TABLE]

where ${\varepsilon({u})}=1/2(\nabla{u}+\nabla{u}^{T})$ is the velocity strain tensor (i.e., the symmetric gradient of the velocity), $\nu$ is the viscosity, $\nabla\cdot u$ denotes the divergence of $u$, $g$ a body force, and $g_{D}$ the value of the velocity on the Dirichlet boundary. The first equation represents the conservation of the linear momentum of the fluid, while the second equation is the incompressibility condition and enforces the mass conversation, see e.g. [24] and references therein.

As in the Darcy case, the high fidelity CutFEM formulation is based on three ingredients: (i) discrete FE spaces, (ii) Nitsche weak imposition of Dirichlet boundary conditions, and (iii) suitable stabilization terms. For what concerns the first topic, using similar notation as in the elliptic case for geometrical quantities, we introduce the discrete spaces ${V}_{h}(\mu)$ and $Q_{h}(\mu)$, for the velocity and the pressure, respectively, as piecewise linear spaces as in [67]:

[TABLE]

Alternative choices are possible, such as, for example, discontinuous Galerkin spaces, see e.g. [68, 69, 70]. For shortness in the remainder of this Section we will drop the suffix $h$ when referring to discrete functions $\upsilon_{h}$, as well as sometimes omit the parameter $\mu$ in our notation.

The Nitsche discrete weak formulation with focus on the Dirichlet boundary conditions of the embedded geometry, takes the form

[TABLE]

being $\gamma_{D}$ a penalty coefficient. For the sake of notation, let us introduce

[TABLE]

where ${\bm{u}}=(u,p)$ and ${\bm{\phi}}=(\psi,\xi)$ denote the pair of velocity and pressure functions.

Finally, introducing suitable stabilization terms, an extended mesh weak form can be expressed using a symmetric form for the velocity and pressure, as follows: find ${\bm{u}(\mu)}=(u(\mu),p(\mu))\in\bm{V}_{h}(\mu):={V}_{h}(\mu)\times Q_{h}(\mu)$ such that

[TABLE]

where

[TABLE]

The terms $j_{u}$ and $j_{p}$ account for the stabilization due to the equal order FE spaces, and are expressed as follows:

[TABLE]

for all test functions $\bm{\phi}=(\psi,\xi)\in{V}_{h}\times Q_{h}$, and for $\gamma_{1,u}$, $\gamma_{1,p}$ positive penalty parameters. We refer to [67] for further details on CutFEM discretization for Stokes problems.

3. POD–Galerkin reduced order method

In this Section we obtain a projection-based reduced order model (ROM) built on the high fidelity discretization introduced in Section 2. Following the ROM procedure of [71], a two stage procedure will be used, the offline and the online. During the offline stage, one examines the solution manifold to construct a reduced basis that well approximates the manifold through a small number of basis functions. This may involve the solution of a large number of high fidelity problems, which may be expensive to query. In contrast, during the online stage, a Galerkin projection of the problem (2.1) onto the space spanned by the reduced basis is required. During this stage each solution for a new value of $\mu$ entails substantially reduced cost.

3.1. The Darcy flow pressure model

Let us start from the elliptic case introduced in Section 2.1. In order to generate the reduced basis space we will employ a compression by means of a Proper Orthogonal Decomposition (POD) [71], even though several other options (e.g. based on greedy procedures, such as the certified Reduced Basis method (RB) [72, 71, 73, 74, 29] or the Proper Generalized Decomposition (PGD) [72, 75, 76]) are available.

The offline phase of the POD–Galerkin ROM consists in an exploration of the solution manifold, obtained by querying the CutFEM high fidelity solver for $M\gg 0$ values of the parameter, and a successive compression to a basis of size $N<M$. The procedure starts by collecting the $M$ high fidelity solutions in the so-called snapshots matrix $\mathcal{S}$, defined as

[TABLE]

Here $\mu^{1},\ldots,\mu^{M}$ are randomly selected values in $\mathcal{K}$ and ${u}(\mu^{i})$ corresponding solutions of (2.6), $i=1,\ldots,M$. We note that, since the solution ${u}(\mu^{i})$ of (2.6) is sought in the $\mu$-dependent FE space (2.4), at least a suitable extension $\widehat{u}(\mu^{i})$ should be carried out to provide snapshots defined on the (common) background mesh $\mathcal{B}_{h}$. Such extensions define a snapshot on the ($\mu$-independent) background FE space

[TABLE]

denoting by $\mathcal{N}_{h}$ its dimension. We will return to this topic in Sections 3.3 and 3.4.

Afterwards, following e.g. [71, 77], a compression by POD is carried out. In practice, this derives the following eigenvalue problem:

[TABLE]

where ${\mathcal{C}}$ is the correlation matrix obtained starting from the snapshots ${\mathcal{S}}$, $\Lambda$ is a diagonal matrix collecting eigenvalues on the diagonal, and ${Q}$ is an orthogonal matrix of the corresponding eigenvectors. The resulting reduced space is then spanned (possibly after a $L^{2}$ normalization) by the columns of the matrix obtained as product between $\mathcal{S}$ and the first $N$ columns of $Q$. The basis functions are denoted by $\varphi_{i}$, $i=1,\ldots,N$, as well as

[TABLE]

denotes the obtained $N$-dimensional space, which will replace the high fidelity space $V_{h}(\mu)$ in all online computations. We remark that $\widehat{V}_{N}$ is a parameter independent space, composed of basis functions defined on the whole background mesh; the former assumption will be relaxed in Section 3.4.

During the online stage, a reduced solution $u_{N}(\mu)$ of the form

[TABLE]

is sought111As the basis functions are global, ${u_{N}}(\mu)$ is defined on the whole background mesh, understanding that its value outside of ${\mathcal{D}}_{\mathcal{T}}(\mu)$ is not interesting and can be discarded during the analysis of the numerical results.. The unknown coefficients $\bm{\alpha}=[\alpha_{1},\ldots,\alpha_{N}]\in\mathbb{R}^{N}$ are obtained through a Galerkin projection of the governing equations onto the reduced basis space, as follows: for any $\mu\in\mathcal{K}$, find $\bm{\alpha}=\bm{\alpha}(\mu)\in\mathbb{R}^{N}$ such that

[TABLE]

Provided that $N\ll\mathcal{N}_{h}$, this linear system is usually very inexpensive to solve. The overall efficiency of the online stage usually comes also from an inexpensive assembly of the left-hand and right-hand side of this linear system, that can be obtained owing to affinity assumptions. Such assumptions, although natural in many applicative cases, can be approximately regained (if lacking) through suitable hyper-reduction procedures [78]. We do mention that numerical examples in Section 4 do not fulfill affinity assumptions. However, being this a preliminary work on the combination of unfitted methods and ROMs, we do not aim at the utmost efficiency at this point. Thus, we will allow inefficient assembly of the left-hand and right-hand side by projecting $\mu$-dependent forms during the online stage. Further perspectives and improvements in this direction will be discussed in the conclusion in Section 5.

3.2. Steady Stokes problem

Let us now discuss the reduction of the Stokes problems presented in Section 2.2. During the offline stage, as in the previous discussion, we collect snapshots from the CutFEM high fidelity discretization. Since the problem is characterized by two unknowns, following e.g. [33, 79] we define two snapshots matrices, namely

[TABLE]

and ultimately perform two separate compressions by means of POD. We remark that two different extension operators, both denoted by $\widehat{\cdot}$ above, may be used for velocity and pressure snapshots, respectively. Here $\mathcal{N}_{h}^{u}$ and $\mathcal{N}_{h}^{p}$ denote the dimension of the background FE spaces

[TABLE]

respectively. Let us denote by $\varphi_{i}^{u}$, $i=1,\ldots,N$ the first $N$ velocity basis functions, as well as by $\varphi_{j}^{p}$, $j=1,\ldots,N$ the first $N$ pressure basis functions.

However, as well known in the reduced basis approximation of saddle point problems, the resulting basis may lead to inaccurate results, especially for what concerns the pressure approximation. Therefore, we resort to a supremizer enrichment procedure [37, 35, 33], which we summarize in the following. From a practical standpoint, this requires the solution of $M$ additional Laplace problems during the offline phase which, in strong form, read as: find $s(\mu^{i})\in V_{h}(\mu^{i})$ such that

[TABLE]

for all $i=1,\ldots,M$. We note that, with the only exception of being a vector problem, (3.4) is a special instance of (2.2). Therefore, we can use the CutFEM discretization introduced in Section 2.1 for its approximation, component by component. The obtained snapshots are then collected in

[TABLE]

where the same extension (and the same extended FE space) used for velocity is also used for supremizers; correspondingly, $N$ basis functions are generated, denoted by $\varphi_{i}^{s}$, $i=1,\ldots,N$ in the following. The resulting reduced basis spaces are therefore obtained as

[TABLE]

Using the combined velocity-pressure notation, we can express the combined reduced basis space as

[TABLE]

as in the following we will not be interested in differentiating the notation between velocity, supremizer and pressure, we equivalently write this space as

[TABLE]

owing to a suitable renumbering of the basis functions. Thus, during the online stage we seek a $3N$-dimensional solution $\bm{u}_{N}(\mu)$ of the form

[TABLE]

such that

[TABLE]

Finally, we remark that alternative approaches to supremizer enrichment are available, see e.g. [80, 34, 81, 79]. In particular, one that relies on the underlying $P^{1}/P^{1}$ (Cut) FEM stabilization [80, 81] seemed attractive in view of providing a less intrusive reduced order model, allowing to neglect the supremizer enrichment stage. However, numerical experiments carried out in the preparation of this work have shown a sensible deterioration of the accuracy when relying only on the high-fidelity stabilization without supremizer enrichment due to the weaker stabilization employed. Although, stronger CutFEM stabilizations are available in the literature, [7, 13].

3.3. Snapshots extension to the background mesh

Let us now discuss a few practical options for what concerns snapshot extension. The first two possibilities that we investigate are to carry out trivial extensions. The first option, which will be called zero extension in the following, is to extend the snapshot to zero in $\mathcal{B}\setminus\mathcal{D}(\mu)$. A slightly different option is to extend snapshots to zero in $\mathcal{B}\setminus\mathcal{D}_{\mathcal{T}}(\mu)$, as the CutFEM solution is naturally defined up to the boundary of $\mathcal{D}_{\mathcal{T}}(\mu)\supset\mathcal{D}(\mu)$. This second options will be named natural smooth extension, as this name is usually employed in the CutFEM literature [82].

The third method we propose is an harmonic extension, see e.g. [83], which, for the scalar case of Section 3.1, reads: for any $\mu\in\mathcal{K}$, assuming the corresponding solution $u(\mu)$ of (2.2) to be known, find $u_{c}(\mu):\mathcal{B}\to\mathbb{R}$ such that

[TABLE]

and then define

[TABLE]

Being another elliptic problem, the harmonic extension problem (3.5) is discretized as described in Section 2.1, up to a suitable replacement of domain, boundaries and problem coefficients. A similar approach may be used for the Stokes problem in Section 2.2, acting component by component.

3.4. Snapshots transportation on the background mesh

Similarly to problems characterized by a strong hyperbolic nature, the methodology described up to now may suffer from a slowly decreasing Kolmogorov $n$-width. To exemplify, let us consider the natural smooth extension, and let $\mu^{1}\neq\mu^{2}$ be two parameters such that $\overline{\mathcal{D}}_{\mathcal{T}}(\mu^{1})\cap\overline{\mathcal{D}}_{\mathcal{T}}(\mu^{2})=\emptyset$. Then, the snapshot corresponding to $\mu^{1}$ bears no useful information for the reduced basis representation of the phenomena associated to $\mu^{2}$, as the support of the solutions do not intersect and the extension is trivial. Even though a non-trivial extension (such as the harmonic one) may alleviate this drawback to some extent, we propose in this Section to combine extension with an additional snapshots preprocessing based on a transportation, as advocated in [60, 61, 62, 63, 64, 65].

Let $\bm{\tau}(\mu):\mathcal{B}\to\mathcal{B}$ be a bijective map chosen such that $\mathcal{D}(\mu)$ is mapped into $\mathcal{D}(\overline{\mu})$, and let $\bm{\tau}^{-1}(\mu)$ denote its inverse with respect to the spatial coordinates. Here $\overline{\mu}$ is a fixed parameter, e.g. the barycenter of $\mathcal{K}$. We will discuss how to devise $\bm{\tau}(\mu)$ from $\mathcal{D}(\mu)$ later on in this Section.

Let $u(\mu)$ be the snapshot obtained having solved (2.2) (the same reasoning goes for Stokes), and extend it to $\widehat{u}(\mu)$ according to the one of the methods introduced in the previous Section. By composition with $\bm{\tau}(\mu)$ the following transported (and extended) snapshot is defined:

[TABLE]

Such transported snapshots $\widehat{u}_{\bm{\tau}}(\mu^{i})$ can then be employed in the definition of the snapshots matrix (3.1) in place of the extended ones $\widehat{u}(\mu^{i})$, $i=1,\ldots,M$. Since all transported snapshots are centered around $\mathcal{D}(\overline{\mu})$, it is expected that a more effective reduction (i.e., faster decay for POD singular values) is obtained.

However, as a result of this transformation, the reduced basis space needs now to depend on $\mu$. Indeed, for any online query associated to a new parameter $\mu$, basis functions need to be transported back from a neighborhood of $\mathcal{D}(\overline{\mu})$ to the corresponding neighborhood of $\mathcal{D}(\mu)$. In particular, the $\mu$-dependent reduced basis space is now defined as

[TABLE]

and inverse transported basis functions $\varphi_{j}\circ\bm{\tau}^{-1}(\mu)$, $j=1,\ldots,N$ are to be employed in the solution of the reduced linear system (3.3), as well as in the representation (3.2) of the reduced order solution. We do remark that such inverse transport introduces an additional non-affinity in the problem formulation; as discussed in Section 3.1, such issue will be a topic of future developments.

Let us now go back to the choice of $\bm{\tau}(\mu)$. As customary in CutFEM methods, $\mathcal{D}(\mu)$ is obtained through a level set method222Although, we remark that the use of a level-set domain description is not an inherent limitation of CutFEMs, see for example [84] where a parametric description of the domain boundary is used.. Being the level sets functions employed in the numerical experiments simple algebraic expressions, it is a matter of simple algebraic manipulations to devise the associated maps $\bm{\tau}(\mu)$. In future, in case of lacking explicit representation of $\bm{\tau}(\mu)$, we seek to resort to optimal transportation procedures [85, 60, 86], as proposed in [62, 87].

Before concluding this Section, we finally stress the difference with a reference domain formulation, as both methods revolve around the idea of introducing a fixed parameter $\overline{\mu}$. In the reference domain formulation, the domain $\mathcal{D}(\overline{\mu})$ associated to fixed parameter $\overline{\mu}$ is the only one employed in the computations, and the differential problems associated to any $\mu$ is recast on $\mathcal{D}(\overline{\mu})$ through suitable pull backs, which are often limited to small deformations. In contrast, the proposed transport procedure requires only a postprocessing of the CutFEM solution, each obtained on its own domain $\mathcal{D}(\mu)$, without any change to the high fidelity solver. The introduction of $\overline{\mu}$ is thus mandatory in the reference domain approach and it is strongly intertwined with the snapshots computation, whilst in our case it is associated only to a desirable postprocessing carried out in order to improve the decay of the Kolmogorov $n$-width.

4. Numerical experiments

4.1. The Darcy flow pressure model

The proposed reduced order technique is numerically tested in this Section on a case characterized by large deformations. The domain $\mathcal{D}(\mu)\subset\mathbb{R}^{2}$ is a parametrized ellipse, defined through the level set function

[TABLE]

where the reference radius $R=0.05$, the length of the axes of the ellipse is parametrized by $(\mu_{1},\mu_{2})\in[0.3,1.8]^{2}$, while the position of the center of the ellipse is parametrized by $(\mu_{3},\mu_{4})\in[-0.85,0.85]^{2}$. A corresponding background domain $\mathcal{B}=[-1.2,1.2]^{2}$ is chosen so that the ellipse is strictly contained in $\mathcal{B}$ for any $\mu=(\mu_{1},\mu_{2},\mu_{3},\mu_{4})$ in the parametric range $\mathcal{K}=[0.3,1.8]^{2}\times[-0.85,0.85]^{2}$. The value $\overline{\mu}=(1,1,0,0)$, corresponding to a circle of radius $R$ centered in the origin, is chosen for what concerns the transport method introduced in Section 3.4, and the transportation mapping333We do note that this expression does not map $\mathcal{B}$ in itself. However, as the goal of the mapping is to transport snapshots in a $R$-neighborhood of the origin, we are satisfied with the proposed mapping as it provides properly mapped values in such neighborhood. Mapped snapshots will be extended to zero for points which are not in $\bm{\tau}(\mu)(\mathcal{B})$. This is trivially through in the zero and natural smooth extensions, while it is in agreement with boundary conditions on $\partial\mathcal{B}$ for the harmonic extension. As an alternative, one could flip the role of extension and transportation, first transporting snapshots and then extending the transported snapshots; we omit this case, as it only affects the harmonic extension case with negligible impact on the overall results and discussion. $\bm{\tau}(x,y;\mu)=(\mu_{1}x+\mu_{3},\mu_{2}y+\mu_{4})$ together with the inverse mapping $\bm{\tau}^{-1}(x,y;\mu)=\left(\frac{x-\mu_{3}}{\mu_{1}},\frac{y-\mu_{4}}{\mu_{2}}\right)$ is considered444In particular, we transport the elliptical geometries (corresponding to the snapshots for various parameters $\mu$ during the offline stage) to circular ones with fixed radius. Indeed, the transformation $x\mapsto\mu_{1}x+\mu_{3}$, $y\mapsto\mu_{2}y+\mu_{4}$ guarantees that the ellipse $\mu^{2}_{2}(x-\mu_{3})^{2}+\mu^{2}_{1}(y-\mu_{4})^{2}-\mu^{2}_{1}\mu^{2}_{2}R^{2}=0$ is mapped into the circle $x^{2}+y^{2}-R^{2}=0$. Afterwards, we collect the transported snapshots and we construct the basis. Then, for every new (in the online stage) geometry we apply the inverse trasportation $x\mapsto({x-\mu_{3}})/{\mu_{1}}$, $y\mapsto({y-\mu_{4}})/{\mu_{2}}$, i.e., we transport basis functions defined on the circular geometry to the original elliptical geometry. , see Figure 2.

This example is characterized by large parametric deformations (e.g., translation of any coordinate of the center of the ellipse can reach values up to $1700\%$ its radius) which a boundary fitted method with reference domain formulation (obtained e.g. by creating a boundary fitted mesh to $\mathcal{D}(\overline{\mu})$ on $\mathcal{B}$ and applying a pull back to the differential problem) would hardly handle. The strong formulation of the problem is as in (2.2) for $g(x,y;\mu)=20$, $g_{D}(x,y;\mu)=0.5+xy$, $\Gamma_{D}(\mu)\equiv\partial\mathcal{D}(\mu)$. For the background domain discretization we used mesh step size $h=0.05$, corresponding to 2806 degrees of freedom. The parameter values $\gamma_{D}$, $\gamma_{N}$, $\gamma_{1}$ under consideration are $10$, [math], $0.1$ respectively.

The results of the offline stage, run over a training set of 400 snapshots, are summarized in Figure 3, where POD eigenvalues (normalized to the maximum eigenvalue) are plotted against the number of modes. The three extension methods introduced in Section 3.3 are considered, as well as their combination to the snapshots transportation procedure introduced in Section 3.4. Results show that small reduced basis spaces cannot be obtained relying only on extension procedures; indeed, Figure 3 shows a very slow eigenvalue decay, with negligible differences among the different extension options (zero extension faring slightly worse than the others). Moreover, the training set is not large enough to provide a representative reduced basis, as the smallest eigenvalue is only 4 orders of magnitude less than the largest one. A faster decay is obtained instead thanks to the transportation. In particular, comparing to the cases with extension only, an improvement of almost three orders of magnitude is obtained for $N=100$, more than four for $N=200$; furthermore, the transport POD drops below numerical precision at $N=300$. Such improvement is instrumental in keeping a low space dimension in the ROM, say $N\approx 100$. Finally, we plot the first six basis functions for the case of natural smooth extension without transportation (Figure 4) and with transportation (Figure 5). From a qualitative point of view, the improvement of the transportation procedure is justified by looking at the extent of the support of the resulting basis functions. Indeed, as expected, basis functions generated through the transportation procedure are compactly supported in a circular neighborhood of the origin, while POD modes obtained without transport are characterized by larger (and varying) support.

Figure 6 shows the results of an error analysis between the reduced order and high fidelity approximations over a testing set of $30$ parameter values, in the $L^{2}$ norm. In particular, the average of the relative error over the testing set is plotted against the reduced basis size. The reduced solution obtained from the zero extension is inaccurate even for $N=140$, being affected by relative errors of the order of $10^{-1}$. A non-zero extension is beneficial, resulting in relative errors that are of the order of $10^{-2}$ for the maximum value of $N$. Furthermore, the combination with inverse transportation allows to further improve results, up to errors of $10^{-4}$ for $N=140$ in the case of POD basis obtained from transport and natural smooth extension. Thus, the pivotal role of snapshots transportation can be inferred from these results, being capable of improving the results of almost three orders of magnitude compared to the simplest zero extension. Nonetheless, all methods reach a plateau after which no further improvements are shown. We claim that this is due to integration errors occurring on $\partial\mathcal{D}(\mu)$ and Nitsche weak imposition of Dirichlet boundary conditions, as the maximum values of the error are consistently attained on parts of the the boundary (see Figure 7 for a representative case).

Further investigation on this issue shows slightly improved relative errors for some specific Nitsche parameters values. In particular, we tested our algorithm for several values of $\gamma_{D}$ in the range $[0.1,40]$. For each coefficient $\gamma_{D}$ we trained the ROM using $400$ snapshots and we queried the reduced solver for a fixed number of 120 basis components and 30 random geometrical parameter test cases. We employed the best extension, namely the natural smooth extension, without transportation (in order to separate the beneficial effects of transportation from those of changing $\gamma_{D}$). As seen in Figure 8, the average relative error shows a mild dependence on the chosen value of $\gamma_{D}$, at least in the range $[5,40]$. Yet, slightly improved results (compared to those of Figure 6) can be obtained for the choice of $\gamma_{D}=8$, resulting in a relative error of $0.0056$, instead of the relative error $0.0112$ for $\gamma_{D}=10$ as in Figure 6. Even though a careful tuning of $\gamma_{D}$ allows to improve the accuracy, such issue should be further investigated in future publications.

4.2. Steady Stokes problem

This experiment considers the Stokes flow around a circular cylinder, embedded in a cavity ${\mathcal{B}}=[-2,2]\times[-1,1]$. The parametrization, which allows to change the $y$ coordinate of the center of the cylinder, is inspired by the numerical example in [57, Section 4]. The embedded domain is thus parametrized through $\mu\in\mathcal{K}=[-0.5,0.5]$ according to the following level set expression:

[TABLE]

where the parameter $\mu$ describes the $y$ coordinate of the center of the embedded circular domain. The radius $R=0.2$, while the viscosity $\nu$ is set to $1$. A constant velocity in the $x$ direction, $u_{\text{in}}=1$ is applied (strongly) at the left side of the domain, and a Neumann boundary condition is applied on the right. No slip, i.e. homogeneous Dirichlet boundary conditions are applied (strongly) on the top and bottom edges, as well as (through Nitsche weak imposition) on the boundary of the embedded cylinder. The results for the test problem have been obtained with a mesh size of $h=0.0350$ for the background mesh ${{\mathcal{B}}_{h}}$, using 15022 triangles for the discretization and $P^{1}/P^{1}$ finite elements in space with stabilization terms as described in Section 2.2.

In this test case we will only employ the natural smooth extension among all the possible extension options presented in Section 3.3, as it was the one which gave best results in Section 4.1. For what concerns the transportation method of Section 3.4, we choose the reference parameter $\overline{\mu}=0$, which corresponds to the center of the parameter range $\mathcal{K}$. The expression $\tilde{y}=\bm{\tau}(y;\mu)=y+\mu(1-y^{2})$ is chosen so that the center-line $y=0$ of the domain is mapped to $\tilde{y}=\mu$, while keeping top and bottom walls fixed. Parameter range $\mathcal{K}$ has been chosen so that the inverse mapping being $\bm{\tau}^{-1}(y;\mu)=\frac{1}{2\mu}(1-(4\mu^{2}-4\mu y+1)^{\frac{1}{2}})$ is well defined.

The results of the offline stage, run over a training set of 600 snapshots, are summarized in Figure 9, where POD eigenvalues (normalized to the maximum eigenvalue) are plotted against the number of modes, for the velocity $u$, supremizer $s$ and pressure $p$. As in the scalar Darcy flow experiments, results show that smaller reduced basis spaces can be obtained by combining extension and transportation procedures; indeed, Figure 9 shows a slower eigenvalue decay for the case without transportation. The difference among the two options are less marked here than in the scalar case due to the simpler parametrization, which only affects one geometrical quantity rather than four. Nonetheless, the qualitative difference between the resulting POD modes is still very appreciable comparing Figures 10-12 (velocity and pressure, respectively, without transportation) to Figures 11-13 (velocity and pressure, respectively, with transportation). In particular, close to the inlet, POD basis functions without transportation are characterized by peaks/sinks in the $y$ direction, their number increasing with the basis index. Such behavior is required for such basis to capture the moving circle. In contrast, no such phenomenon is present when applying transportation to the snapshots.

Figure 14 shows the $L^{2}$ relative errors results, for velocity and pressure, of an error analysis between the reduced order and high fidelity approximations over a testing set of $100$ parameter values. In particular, the average of the relative error over the testing set is plotted against the reduced basis size. In the case of natural smooth extension, the reduced pressure solution obtained is still inaccurate even for $N=50$ even with the addition of supremizers, being affected by relative errors of the order of $10^{-1}$. Reduced velocity solution is more accurate, resulting in errors of the order of $10^{-2}$ for $N\geq 20$. In contrast, the combination to the transportation procedure results in relative errors that is of the order of $10^{-2}$ starting from $N=3$ for both velocity and pressure. Thus, again the role of snapshots transportation can be inferred from these results, being capable of improving the results of almost an order of magnitude compared to the natural smooth extension for a very low number of basis functions. Unfortunately, both methods again reach a plateau after which no further improvement is shown. As in the elliptic case, we claim that this is due to weak imposition of boundary conditions, as errors tend to concentrate in a neighborhood of the embedded circle (especially for the pressure, which will thus affect the downstream velocity; see Figures 15-16 for a representative case); further work on this topic is forthcoming.

5. Conclusions and future developments

In this work a POD-Galerkin ROM based on CutFEM high fidelity simulations was presented for linear PDE problems (elliptic and Stokes), characterized by a geometrical parametrization with (possibly) large variations. A CutFEM discretization naturally allows to use a level set description of the parametrized geometry. In our opinion, this results in a simpler and more versatile high fidelity method when compared to a FE formulation with pull back to a reference domain.

A separation between construction and evaluation phases is sought. While the evaluation of the ROM follows a standard Galerkin projection, a careful adaptation to embedded methods has been necessary for the construction stage, especially for what concerns the definition of snapshot on a common background mesh through a suitable combination of extension to the background mesh (in order to have all snapshots defined on a common mesh) and transportation on the background mesh itself (in order to enforce a rapid decrease of the Kolmogorov $n$-width). Indeed, thanks in particular to the transportation step, the developed ROM is able to reproduce the high fidelity solution in an accurate manner, with relative errors of the order of $10^{-4}$ for the elliptic case and $10^{-2}$ for the Stokes case.

Even though the proposed ROM follows a construction-evaluation paradigm, it is not offline-online separable in the usual sense [71], yet. Such separability, which is required for the efficient evaluation of the reduced order system, is a perspective to be thoroughly investigated in future. Suitable hyper reduction techniques, such as the empirical interpolation method [78, 88], should be employed for this goal. Certification of the error is another topic of remarkable interest, in view of the application of greedy algorithms during the generation of the reduced basis space. Furthermore, a deep investigation on the role of Nitsche terms on the accuracy of the resulting ROM is needed, as numerical results show a plateau due to the weak imposition of Dirichlet boundary conditions.

As a further future development, we mention the extension of the proposed ROM to nonlinear problems in fluid dynamics, such as the Navier-Stokes equations. Moreover, a prototypical case for the application of embedded methods are fluid-structure interaction problems, which is another future perspective of this work. This will introduce additional complexities such as additional equations to account for the structural dynamics, as well as fluid structure-interaction coupling.

Acknowledgments

We acknowledge the support by European Union Funding for Research and Innovation – Horizon 2020 Program – in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” and FSE project - European Social Fund - HEaD ”Higher Education and Development” SISSA operazione 1, Regione Autonoma Friuli Venezia-Giulia. We also acknowledge INdAM-GNCS projects “Tecniche di Riduzione di Modello per Applicazioni Mediche”, “Metodi numerici avanzati combinati con tecniche di riduzione computazionale per PDEs parametrizzate e applicazioni” and HFRI and GSRT under grant agreement No 1115. Numerical simulations were carried out using the ngsxfem extension of ngsolve [89, 90] for the high fidelity simulations, and RBniCS [91] for the reduced order ones. We acknowledge developers and contributors to both libraries. We thank in particular Prof. Christoph Lehrenfeld, Institute for Computational and Applied Mathematics, Göttingen, for useful instructions regarding the ngsxfem library.