Generalized parametric solutions in Stokes flow

TL;DR

This paper develops a PGD-based method for efficiently solving parametric Stokes flow problems, enabling rapid evaluations for design optimization and uncertainty quantification in computational mechanics.

Contribution

It introduces a PGD formulation tailored for the saddle point structure of the Stokes problem, analyzing different separated forms and demonstrating effectiveness through numerical examples.

Findings

Successful application to Stokes and Brinkman models

Efficient parametric solutions with reduced computational cost

Enhanced capability for design optimization and uncertainty quantification

Abstract

Design optimization and uncertainty quantification, among other applications of industrial interest, require fast or multiple queries of some parametric model. The Proper Generalized Decomposition (PGD) provides a separable solution, a \emph{computational vademecum} explicitly dependent on the parameters, efficiently computed with a greedy algorithm combined with an alternated directions scheme and compactly stored. This strategy has been successfully employed in many problems in computational mechanics. The application to problems with saddle point structure raises some difficulties requiring further attention. This article proposes a PGD formulation of the Stokes problem. Various possibilities of the separated forms of the PGD solutions are discussed and analyzed, selecting the more viable option. The efficacy of the proposed methodology is demonstrated in numerical examples for both Stokes and Brinkman models.

Equations

$$\{\begin{array}{rlrl}-\text{boldmath∇unboldmath}\cdot \nu \text{boldmath∇unboldmath}\text{boldmathuunboldmath}+\text{boldmath∇unboldmath}p& =\text{boldmathbunboldmath}& & \text{in Ω}\\ \text{boldmath∇unboldmath}\cdot \text{boldmathuunboldmath}& =0& & \text{in Ω}\\ \text{boldmathuunboldmath}& ={\text{boldmathuunboldmath}}_D& & \text{on ΓD}\\ -p\text{boldmathnunboldmath}+\nu \text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}\text{boldmathuunboldmath}& =\text{boldmathtunboldmath}& & \text{on ∂Ω∖ΓD}.\end{array}$$

$$\left\{\begin{aligned} \text{\boldmath$u$\unboldmath}&\in\text{\boldmath$\mathcal{V}$\unboldmath}{\color[rgb]{0,0,0}:=[\mathcal{V}\otimes\mathcal{L}_{2}(\mathcal{I}_{1})\otimes\mathcal{L}_{2}(\mathcal{I}_{2})\otimes\dotsb\otimes\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}})]^{\texttt{n}_{\texttt{sd}}}\color[rgb]{0,0,0}}\text{, and}\\ p&\in\mathcal{L}_{2}(\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath})=\mathcal{L}_{2}(\Omega)\otimes\mathcal{L}_{2}(\mathcal{I}_{1})\otimes\mathcal{L}_{2}(\mathcal{I}_{2})\otimes\dotsb\otimes\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}}),\end{aligned}\right.$$

$$A\bigl{(}\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath}\bigr{)}+B\bigl{(}\text{\boldmath$v$\unboldmath},p\bigr{)}+B\bigl{(}\text{\boldmath$u$\unboldmath},q\bigr{)}=L\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)},\;\forall(\text{\boldmath$v$\unboldmath},q)\in\text{\boldmath$\mathcal{S}$\unboldmath}\times\mathcal{L}_{2}(\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}),$$

$$\begin{split}A\bigl{(}\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath}\bigr{)}&:=\int_{\mathcal{I}_{1}}\!\int_{\mathcal{I}_{2}}\dotsi\int_{\mathcal{I}_{\texttt{n}_{\texttt{pa}}}}a\bigl{(}\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath}\bigr{)}\,d\mu_{\texttt{n}_{\texttt{pa}}}\dotsm d\mu_{2}\,d\mu_{1},\\ B\bigl{(}\text{\boldmath$u$\unboldmath},q\bigr{)}&:=\int_{\mathcal{I}_{1}}\!\int_{\mathcal{I}_{2}}\dotsi\int_{\mathcal{I}_{\texttt{n}_{\texttt{pa}}}}b\bigl{(}\text{\boldmath$u$\unboldmath},q\bigr{)}\,d\mu_{\texttt{n}_{\texttt{pa}}}\dotsm d\mu_{2}\,d\mu_{1},\\ L\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)}&:=\int_{\mathcal{I}_{1}}\!\int_{\mathcal{I}_{2}}\dotsi\int_{\mathcal{I}_{\texttt{n}_{\texttt{pa}}}}\ell\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)}\,d\mu_{\texttt{n}_{\texttt{pa}}}\dotsm d\mu_{2}\,d\mu_{1},\end{split}$$

$$\begin{split}a\bigl{(}\text{\boldmath$v$\unboldmath},\text{\boldmath$w$\unboldmath}\bigr{)}=\int_{\Omega}2\nu\,\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$v$\unboldmath}:\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$w$\unboldmath}\,d\Omega\,,&\quad b\bigl{(}\text{\boldmath$v$\unboldmath},q\bigr{)}=-\int_{\Omega}q\,\text{\boldmath$\nabla$\unboldmath}\cdot\text{\boldmath$v$\unboldmath}\,d\Omega,\quad\textrm{ and}\\ \ell\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)}=\int_{\Omega}\text{\boldmath$s$\unboldmath}\cdot\text{\boldmath$v$\unboldmath}\,d\Omega+&\int_{\partial\Omega\setminus\Gamma_{D}}\text{\boldmath$v$\unboldmath}\cdot\text{\boldmath$t$\unboldmath}\,d\Gamma.\end{split}$$

$$\left\{\begin{aligned} A\bigl{(}\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath}\bigr{)}+B\bigl{(}\text{\boldmath$v$\unboldmath},p\bigr{)}&=L\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)}&&\forall\text{\boldmath$v$\unboldmath}\in\text{\boldmath$\mathcal{S}$\unboldmath},\\ B\bigl{(}\text{\boldmath$u$\unboldmath},q\bigr{)}&=0&&\forall q\in\mathcal{L}_{2}(\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}),\end{aligned}\right.$$

$$\{\begin{array}{rl}{u}_i(\text{boldmathxunboldmath},\mu )\approx u_{i_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& =u_{i_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+F_{\text{boldmathuunboldmath},i}^n(\text{boldmathxunboldmath})L_{\text{boldmathuunboldmath},i,1}^n({\mu }_1)L_{\text{boldmathuunboldmath},i,2}^n({\mu }_2)\cdots L_{\text{boldmathuunboldmath},i,{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}}),\\ p(\text{boldmathxunboldmath},\mu )\approx p_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& =p_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+F_p^n(\text{boldmathxunboldmath})L_{p,1}^n({\mu }_1)L_{p,2}^n({\mu }_2)\cdots L_{p,{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}})\end{array}$$

$$\{\begin{array}{rl}\text{boldmathuunboldmath}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})\approx {\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\text{boldmathμunboldmath})& ={\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})+{\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}^n(\text{boldmathxunboldmath})L_{\text{boldmathuunboldmath},1}^n({\mu }_1)L_{\text{boldmathuunboldmath},2}^n({\mu }_2)\cdots L_{\text{boldmathuunboldmath},{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}}),\\ p(\text{boldmathxunboldmath},\text{boldmathμunboldmath})\approx p_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\text{boldmathμunboldmath})& =p_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})+F_p^n(\text{boldmathxunboldmath})L_{p,1}^n({\mu }_1)L_{p,2}^n({\mu }_2)\cdots L_{p,{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}}).\end{array}$$

$$b\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}^{s},q\bigr{)}=0\quad\forall q\in\mathcal{L}_{2}(\Omega)\text{ and $s=0,\dotsc,n$}\Longrightarrow b\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n},q\bigr{)}=0\quad\forall q\in\mathcal{L}_{2}(\Omega),$$

$$\{\begin{array}{rl}\text{boldmathuunboldmath}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})\approx {\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\text{boldmathμunboldmath})& ={\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})+{\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}^n(\text{boldmathxunboldmath})L_1^n({\mu }_1)L_2^n({\mu }_2)\cdots L_{{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}}),\\ p(\text{boldmathxunboldmath},\text{boldmathμunboldmath})\approx p_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\text{boldmathμunboldmath})& =p_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})+F_p^n(\text{boldmathxunboldmath})L_1^n({\mu }_1)L_2^n({\mu }_2)\cdots L_{{\text{n}}_{\text{pa}}}^n({\mu }_{{\text{n}}_{\text{pa}}}).\end{array}$$

$$\{\begin{array}{rl}\text{boldmathuunboldmath}(\text{boldmathxunboldmath},\mu )\approx {\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& ={\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+{\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}(\text{boldmathxunboldmath})L_{\text{boldmathuunboldmath}}(\mu ),\\ p(\text{boldmathxunboldmath},\mu )\approx p_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& =p_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+F_p(\text{boldmathxunboldmath})L_p(\mu ).\end{array}$$

$$A\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},\text{\boldmath$v$\unboldmath}\bigr{)}+B\bigl{(}\text{\boldmath$v$\unboldmath},F_{\!p}\,L_{p}\bigr{)}+B\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},q\bigr{)}=R\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1},\text{\boldmath$v$\unboldmath}\bigr{)}-B\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},q\bigr{)},$$

$$R\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1},\text{\boldmath$v$\unboldmath}\bigr{)}=L\bigl{(}\text{\boldmath$v$\unboldmath}\bigr{)}-A\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},\text{\boldmath$v$\unboldmath}\bigr{)}-B\bigl{(}\text{\boldmath$v$\unboldmath},p_{{}_{\texttt{PGD}}}^{n-1}\bigr{)}.$$

$$\text{boldmathvunboldmath}=\delta {\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}{L}_{\text{boldmathuunboldmath}}+{\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}\delta L_{\text{boldmathuunboldmath}}\text{ and }q=\delta F_p{L}_p+F_p\delta L_p.$$

$$\displaystyle[left={\empheqlbrace\,}]A\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}}\bigr{)}{+}B\bigl{(}\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},F_{\!p}\,L_{p}\bigr{)}$$

$$\displaystyle B\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},\delta\!F_{\!p}\,L_{p}\bigr{)}$$

$$\displaystyle A\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,\delta\!L_{\text{\boldmath$u$\unboldmath}}\bigr{)}{+}B\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,\delta\!L_{\text{\boldmath$u$\unboldmath}},F_{\!p}\,L_{p}\bigr{)}$$

$$\displaystyle B\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},F_{\!p}\,\delta\!L_{p}\bigr{)}$$

$$\begin{split}A\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}L_{\text{\boldmath$u$\unboldmath}},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}L_{\text{\boldmath$u$\unboldmath}}\bigr{)}&=\int_{\mathcal{I}}\!L_{\text{\boldmath$u$\unboldmath}}L_{\text{\boldmath$u$\unboldmath}}a\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\bigr{)}\,d\mu=\left[\int_{\mathcal{I}}L_{\text{\boldmath$u$\unboldmath}}^{2}\,d\mu\right]a\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\bigr{)},\\ B\bigl{(}\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}},F_{\!p}\,L_{p})&=\left[\int_{\mathcal{I}}L_{\text{\boldmath$u$\unboldmath}}\,L_{p}\,d\mu\right]\,b\bigl{(}\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,F_{\!p})\end{split}$$

$${\alpha }_{\mu }={\int }_{\mathcal{I}}{L}_{\text{boldmathuunboldmath}}^{2}d\mu ,{\beta }_{\mu }={\int }_{\mathcal{I}}{L}_{\text{boldmathuunboldmath}}{L}_{p}d\mu ,$$

$$\displaystyle[left={\empheqlbrace\,}]\alpha_{\mu}a\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\bigr{)}+\beta_{\mu}b\bigl{(}\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,F_{\!p})$$

$$\displaystyle=:\mathcal{R}_{\text{\boldmath$u$\unboldmath}}\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1},\delta\!\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\,L_{\text{\boldmath$u$\unboldmath}}\bigr{)}$$

$$\displaystyle\beta_{\mu}b\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},\delta\!F_{\!p}\bigr{)}$$

$$\displaystyle=:\mathcal{R}_{p}\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},\delta\!F_{\!p}\,L_{p}\bigr{)}$$

$$\displaystyle\begin{split}\mathcal{R}_{\text{\boldmath$u$\unboldmath}}\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1},\text{\boldmath$w$\unboldmath}\,\omega\bigr{)}=\biggr{[}\int_{\mathcal{I}}\omega\,d\mu\biggl{]}\,\ell\bigl{(}\text{\boldmath$w$\unboldmath}\bigr{)}-\sum_{s=0}^{n-1}&\biggr{[}\int_{\mathcal{I}}\omega\,L_{\text{\boldmath$u$\unboldmath}}^{s}\,d\mu\biggl{]}\,a\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}^{s},\text{\boldmath$w$\unboldmath}\bigr{)}\\ &-\sum_{s=0}^{n-1}\biggr{[}\int_{\mathcal{I}}\omega\,L_{p}^{s}\,d\mu\biggl{]}\,b\bigl{(}\text{\boldmath$w$\unboldmath},F_{\!p}^{s}\bigr{)},\end{split}$$

$$\displaystyle\mathcal{R}_{p}\bigl{(}\text{\boldmath$u$\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},q\,\rho\bigr{)}=-\sum_{s=0}^{n-1}\biggr{[}\int_{\mathcal{I}}L_{\text{\boldmath$u$\unboldmath}}^{s}\,\rho\,d\mu\biggl{]}b\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}^{s},q\bigr{)}.$$

$$\left(\begin{array}{cc}\mathbf{K}& \mathbf{G}\\ {\mathbf{G}}^T& 0\end{array}\right).$$

$$[left=\text{empheqlbrace}]{\alpha }_{\text{boldmathuunboldmath}}{\int }_{\mathcal{I}}\delta L_{\text{boldmathuunboldmath}}{L}_{\text{boldmathuunboldmath}}^{n}d\mu +{\beta }_{\text{boldmathuunboldmath}}{\int }_{\mathcal{I}}\delta L_{\text{boldmathuunboldmath}}{L}_p^{n}d\mu$$

$${\beta }_{\text{boldmathuunboldmath}}{\int }_{\mathcal{I}}\delta L_p{L}_{\text{boldmathuunboldmath}}^{n}d\mu$$

$$\alpha_{\text{\boldmath$u$\unboldmath}}=a\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}}\bigr{)}\;,\qquad\beta_{\text{\boldmath$u$\unboldmath}}=b\bigl{(}\text{\boldmath$F$\unboldmath}_{\!\!\text{\boldmath$u$\unboldmath}},F_{\!p}\bigr{)}\;,$$

$$\left(\begin{array}{cc}{\mathbf{M}}_{\text{boldmathuunboldmath},\text{boldmathuunboldmath}}& {\mathbf{M}}_{\text{boldmathuunboldmath},p}\\ {\mathbf{M}}_{\text{boldmathuunboldmath},p}^T& 0\end{array}\right).$$

$$\left(\begin{array}{cc}\mathbf{M}& \mathbf{M}\\ \mathbf{M}& 0\end{array}\right)\left(\begin{array}{c}{\text{boldmathLunboldmath}}_u\\ {\text{boldmathLunboldmath}}_p\end{array}\right)=\left(\begin{array}{c}{\text{boldmathRunboldmath}}_u\\ \text{boldmath0unboldmath}\end{array}\right).$$

$$[left=\text{empheqlbrace}]{\alpha }_{\text{boldmathuunboldmath}}{L}_{\text{boldmathuunboldmath}}^n(\mu )+{\beta }_{\text{boldmathuunboldmath}}{L}_p^n(\mu )$$

$${\beta }_{\text{boldmathuunboldmath}}{L}_{\text{boldmathuunboldmath}}^n(\mu )$$

$$\{\begin{array}{rl}\text{boldmathuunboldmath}(\text{boldmathxunboldmath},\mu )\approx {\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& ={\text{boldmathuunboldmath}}_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+{\text{boldmathFunboldmath}}_{\text{boldmathuunboldmath}}(\text{boldmathxunboldmath})L(\mu ),\\ p(\text{boldmathxunboldmath},\mu )\approx p_{{}_{\text{PGD}}}^n(\text{boldmathxunboldmath},\mu )& =p_{{}_{\text{PGD}}}^{n-1}(\text{boldmathxunboldmath},\mu )+F_p(\text{boldmathxunboldmath})L(\mu ).\end{array}$$

Full text

Keywords: Reduced order model, Parametric Solution, Stokes flow, Proper Generalized Decomposition

1 Introduction

Standard discretization techniques in computational mechanics have reached an amazing level of maturity and efficiency. Nonetheless, the systematic exploration of parametric solutions arising from optimization (where the best choice for the parameters is unknown) or uncertainty quantification (where the parameters have stochastic features) is often computationally unaffordable. The Proper Generalized Decomposition (PGD), among other Reduced Order Models, provides a generalized solution with an explicit parametric dependence. This compact expression containing the analytical dependence on the free parameters is also known as computational vademecum and allows an expedited exploration of the parametric space, with the computational cost of a simple interpolation, i.e. post-processing.

The PGD has been successfully employed in different problems in the broad field of mathematical and computational modelling. Essentially, PGD consists in finding a separable approximation, that is a sum of terms, each of them being a product of modal functions depending on one of the parameters. This approximation is usually computed with a greedy algorithm (obtaining the terms sequentially) and, for each term, an alternated directions iterative scheme is, in general, employed to find the different parametric modes.

PGD was already used in the framework of Stokes and Navier–Stokes governed problems to obtain separated solutions in terms of the different spatial dimensions for Cartesian domains [1, 2, 3] and also for space-time separation [4, 5]. PGD for space-space separation and space-time separation is significantly increasing the computational efficiency in solving problem with complex flow patterns in simple cartesian domains.

Here, the focus is on solving parametric problems in complex domains with arbitrary geometries (not assumed to be Cartesian). Therefore the space coordinates are treated together and separated from the different (independent) parametric dimensions.

This paper aims at analyzing the application of PGD to problems with saddle point structure, taking the Stokes problem as one of the simplest. In particular, special attention is paid to the selection of the form of parametric separation in hybrid formulations. In other words, in a velocity-pressure formulation, the question is: must the parametric modes be independent for velocity and pressure, or just one for both?

Thus, the remainder of the paper is structured as follows. The parametric version of the Stokes problem is stated in Section 2. Section 3 describes the general formulation of PGD in this case and analyzes the possible alternatives for the parametric separation, concluding that the same parametric mode must be used for both velocity and pressure. Section 4 presents examples demonstrating the viability of the devised approach.

2 The parameterized Stokes flow

The strong form of the Stokes problem can be written as

[TABLE]

User-prescribed data are the computational domain $\Omega\subset\mathbb{R}^{\texttt{n}_{\texttt{sd}}}$ ($\texttt{n}_{\texttt{sd}}$ being the number of spatial dimensions) whose boundary $\partial\Omega$ is partitioned into Dirichlet, $\Gamma_{D}$, and Neumann frontiers, the body forces $s$, the Dirichlet, \text{\boldmathu\unboldmath}_{D}, and Neumann, $t$, boundary conditions, and the kinematic viscosity $\nu$.

Any of these user-prescribed data could be a function of a set of parameters \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}} (with $\texttt{n}_{\texttt{pa}}$ number of parameters). Those affecting the right-hand-side of the resulting equations (viz. $s$, \text{\boldmathu\unboldmath}_{D} and $t$) are easy to handle. On the contrary those affecting the differential operator (viz. viscosity or domain) cannot be treated trivially. The set \text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}}, which characterizes the admissible range for parameters $\mu$, can be defined as the cartesian product of the range for each parameter, namely, \text{\boldmath\mathcal{I}\unboldmath}:=\mathcal{I}_{1}\times\mathcal{I}_{2}\times\dotsb\times\mathcal{I}_{\texttt{n}_{\texttt{pa}}} with $\mu_{i}\in\mathcal{I}_{i}$ for $i=1,\dotsc,\texttt{n}_{\texttt{pa}}$.

This can be interpreted as taking $\mu$ as additional independent variables (or parametric coordinates) instead of problem parameters. Hence, the unknown velocity-pressure pair (\text{\boldmathu\unboldmath},p) can be seen as functions in a larger dimensional space and can be expressed as \text{\boldmathu\unboldmath}(\text{\boldmathx\unboldmath},\text{\boldmath\mu\unboldmath}) and p(\text{\boldmathx\unboldmath},\text{\boldmath\mu\unboldmath}) with (\text{\boldmathx\unboldmath},\text{\boldmath\mu\unboldmath})\in\Omega\times\text{\boldmath\mathcal{I}\unboldmath}.

Consequently, formally $u$ and $p$ lie in tensor product spaces, namely

[TABLE]

where [\mathcal{V}]^{\texttt{n}_{\texttt{sd}}}:=\{\text{\boldmathu\unboldmath}\in[\mathcal{H}^{1}(\Omega)]^{\texttt{n}_{\texttt{sd}}}:\text{\boldmathu\unboldmath}=\text{\boldmathu\unboldmath}_{D}\text{ on }\Gamma_{D}\}. Note that in the definition of $\mathcal{V}$ all the spaces in the tensorial product are raised to the power of $\texttt{n}_{\texttt{sd}}$, also the parametric ones. A standard weighted residuals approach, with integrals in \Omega\times\text{\boldmath\mathcal{I}\unboldmath} and the usual integration by parts only in $\Omega$ produces a (spatially) weak form in this multi-dimensional setup. Namely, find (\text{\boldmathu\unboldmath},p)\in\text{\boldmath\mathcal{V}\unboldmath}\times\mathcal{L}_{2}(\Omega\times\text{\boldmath\mathcal{I}\unboldmath}) such that

[TABLE]

where the test function space for velocities is \text{\boldmath\mathcal{S}\unboldmath}:=[\mathcal{H}^{1}_{\Gamma_{D}}]^{\texttt{n}_{\texttt{sd}}}\otimes[\mathcal{L}_{2}(\mathcal{I}_{1})]^{\texttt{n}_{\texttt{sd}}}\otimes[\mathcal{L}_{2}(\mathcal{I}_{2})]^{\texttt{n}_{\texttt{sd}}}\otimes\dotsb\otimes[\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}})]^{\texttt{n}_{\texttt{sd}}} and [\mathcal{H}^{1}_{\Gamma_{D}}]^{\texttt{n}_{\texttt{sd}}}:=\{\text{\boldmathu\unboldmath}\in[\mathcal{H}^{1}(\Omega)]^{\texttt{n}_{\texttt{sd}}}:\text{\boldmathu\unboldmath}=\text{\boldmath0\unboldmath}\text{ on }\Gamma_{D}\}.

The following definitions of the bilinear and linear forms are necessary:

[TABLE]

where

[TABLE]

Obviously, the number of dimensions of the solution domain increases with the number of parameters. To circumvent the curse of dimensionality, the PGD approach [6, 7, 8, 9] is employed here. This approach assumes a separable structure in the approximation to (\text{\boldmathu\unboldmath},p). Note that the tensor product spaces $\mathcal{V}$ and \mathcal{L}_{2}(\Omega\times\text{\boldmath\mathcal{I}\unboldmath}) inherit the multidimensional complexity of the problem and, in principle, do not assume separability of the functions.

Remark 1** (Saddle point structure).**

Note that equation (3) is often written as

[TABLE]

to evidence the saddle point problem at hand.

3 The Proper Generalized Approximation

3.1 Three alternative forms of the separated approximation

As usual in a PGD strategy, a separated representation (\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n}) is imposed to approximate the solution of (3) in each subdomain. The couple (\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n}) stands for the PGD approximation with $n$ terms (or modes) of the velocity-pressure couple and it is defined as a sum of separated terms. Each term (mode) is the product of functions depending only on one of the arguments (\text{\boldmathx\unboldmath},\mu_{1},\mu_{2},\dotsc,\mu_{\texttt{n}_{\texttt{pa}}}). Note that, in some of the PGD implementations the separated modal functions are normalized and therefore a scalar coefficient affects each mode and characterizes its amplitude. The first mode, (\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{0},p_{{}_{\texttt{PGD}}}^{0}), is arbitrarily chosen (for instance accounting for Dirichlet boundary conditions). Then, a greedy algorithm is implemented to compute successively the last one, that is to compute term $n$ assuming that term $n-1$ is available [8, 9].

Three alternatives can be considered for this separation, see also [10], depending on how the modal functions for parameters $\mu$ are considered.

Independent component-wise separation (case #0):

A distinct parametric modal function is considered for each term of the PGD expansion providing (\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n}), each component of the velocity and the pressure ($\texttt{n}_{\texttt{sd}}+1$ components) and each parameter ($\texttt{n}_{\texttt{pa}}$). Hence the total number of parametric modal functions is $n(\texttt{n}_{\texttt{sd}}+1)\texttt{n}_{\texttt{pa}}$. This is the most general separation because a different parameter function is considered for each velocity component and for pressure, namely, for $i=1,\dotsc,\texttt{n}_{\texttt{sd}}$

[TABLE]

where $F$ denote the spatial functions (depending on $x$) characterizing the mode (subscripts indicate wether they refer to some component of the velocity or the pressure) and $L$ denote the parametric functions. All these are real-valued scalar fields. In the following the vector fields are boldfaced. This strategy imposes different parameter–dependent functions for each spatial velocity component. Consequently, the spatial differential operators are affected differently for each spatial component. This induces enormous difficulties in order to express (at least well approximately) the different forms by the sum of products of parameter-dependent functions and parameter-independent operators.

The major drawback of such an approach is that the incompressibility constrain will not be trivially enforced. That is, LBB or incompressibility stabilization must be specifically studied and it is not trivial due to the variable weightings introduced by the parameter functions.

Moreover, the implementation of the alternated directions scheme in the PGD for this separation form is highly intrusive. This is because it requires a distinct treatment of the different directions, depending on the value of the assumed parametric modes. Thus, the use of non-intrusive strategies with commercial codes becomes much more involved.

Unique parameter function for velocity independent

from the ones for pressure (case #1): With respect to the previous formulation, the functions affecting the $\texttt{n}_{\texttt{sd}}$ different components of the velocity are taken to be the same. Thus, the total number of parametric modal functions is equal to $n\,2\,\texttt{n}_{\texttt{pa}}$. Correspondingly,

[TABLE]

This approach uncouples naturally the parameter functions from the spatial ones when the spatial divergence is computed to impose incompressibility. This will have a major effect in simplifying the choice of the spatial spaces pairs for velocity and pressure.

Remark 2** (Divergence-free modes).**

Note that, given the separated representation of the PGD approximation, the point-wise divergence-free velocity, \text{\boldmath\nabla\unboldmath}\cdot\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n}=0, is guaranteed by imposing \text{\boldmath\nabla\unboldmath}\cdot\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{s}=0 for $s=0,\dotsc,n$. Likewise

[TABLE]

which imposes the usual weak divergence-free condition on $\Omega$ for any set of parameters \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}.

Unique parameter functions for velocity and pressure (case #2):

Following with the simplification in the number of parameter functions at each mode, the next step is to employ the same function for every component of the velocity and also for pressure, with a total number of parametric modal functions equal to $n\,\texttt{n}_{\texttt{pa}}$, namely.

[TABLE]

Obviously, this case also benefits from the separation of the divergence of the velocity as in the previous approximation.

As stated above, the first form of the separation (with different parameter modes affecting every component of the velocity and the pressure) leads to a cumbersome formulation requiring a highly intrusive implementation. Consequently, in the following, only the two latter alternatives are taken into consideration. For the sake of a simpler presentation and without any loss of generality, the subsequent developments are done for the particular case of only one parameter $\mu$ ($\texttt{n}_{\texttt{pa}}=1$). Thus, the first alternative under consideration (case #1) uses two parameter functions (one for $u$ and one for $p$), and the second alternative (case #2) uses just one parameter function (the same for $u$ and $p$). In the general case of $\texttt{n}_{\texttt{pa}}\geq 1$, the number of parameter functions in each case are $2\,\texttt{n}_{\texttt{pa}}$ and $\texttt{n}_{\texttt{pa}}$, respectively.

3.2 Case #1: Two parameter functions (one for $u$ and one for $p$)

For $\texttt{n}_{\texttt{pa}}=1$, the PGD approximation is written in this case as

[TABLE]

Note that, in order to shorten the writing, superscript $n$ is omitted in the notation of the unknown functions \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{n}(\text{\boldmathx\unboldmath}), L_{\text{\boldmathu\unboldmath}}^{n}(\mu), F_{\!p}^{n}(\text{\boldmathx\unboldmath}) and $L_{p}^{n}(\mu)$.

3.2.1 Solving for each mode

The approximation defined in (7) is substituted in (3) and tested in a tangent manifold. That is, the unknowns to be determined are \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}\in[\mathcal{V}]^{\texttt{n}_{\texttt{sd}}}, L_{\text{\boldmathu\unboldmath}}\in\mathcal{L}_{2}(\mathcal{I}), $F_{\!p}\in\mathcal{L}_{2}(\Omega)$, and $L_{p}\in\mathcal{L}_{2}(\mathcal{I})$ such that

[TABLE]

for all $v$ and $q$ in the tangent manifold and being the residual $R(\cdot,\cdot,\cdot)$ defined by

[TABLE]

The tangent manifold is readily characterized by choosing $v$ and $q$ as variations of \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}\,L_{\text{\boldmathu\unboldmath}} and $F_{\!p}\,L_{p}$ respectively, that is

[TABLE]

for all \delta\!\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}\in[\mathcal{H}^{1}_{\Gamma_{D}}]^{\texttt{n}_{\texttt{sd}}}, \delta\!L_{\text{\boldmathu\unboldmath}}\in\mathcal{L}_{2}(\mathcal{I}), $\delta\!F_{\!p}\in\mathcal{L}_{2}(\Omega)$, and $\delta\!L_{p}\in\mathcal{L}_{2}(\mathcal{I})$.

Following Remark 1, this problem can also be equivalently rewritten as

[TABLE]

Note that (9) is a nonlinear system of functional equations for the four unknowns \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}, L_{\text{\boldmathu\unboldmath}}, $F_{\!p}$, and $L_{p}$. In the PGD framework, (9) is iteratively solved using an alternated directions scheme. That is, first solving (9a) and (9b) for unknowns \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$, assuming that L_{\text{\boldmathu\unboldmath}} and $L_{p}$ are known. This first stage is denoted spatial iteration because it has the same structure of a standard (non parametric) Stokes problem.

Then, equations (9c) and (9d) are solved for unknowns L_{\text{\boldmathu\unboldmath}} and $L_{p}$ assuming that \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$ are known. This step is denoted parameter iteration, and it consists in iterating for every parametric direction (just once for $\texttt{n}_{\texttt{pa}}=1$, in general $\texttt{n}_{\texttt{pa}}$ steps are needed). The process is iterated between subsystem (9a) and (9b) and subsystem (9c) and (9d) until a stationary solution is reached.

Remark 3** (Solving groups of two equations).**

In other PGD formulations, the alternated direction schemes for the nonlinear systems take the modes one by one, solving for one and assuming that the rest are known. Here, the two couples of unknowns (\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}},F_{\!p}) and (L_{\text{\boldmathu\unboldmath}},L_{p}) are solved together. This is due to the Saddle Point structure inherited by the groups of equations (9a) & (9b) and (9c) & (9d). In particular, the natural unknown for (9b) would be $F_{\!p}$ and it is not appearing explicitly in the equation. Thus, it is not possible solving (9b) to find $F_{\!p}$ assuming that \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}, L_{\text{\boldmathu\unboldmath}} and $L_{p}$ are known. The same happens with $L_{p}$ in (9d).

3.2.2 The spatial iteration.

As stated above, the spatial iteration consists in solving (9a) and (9b) for unknowns \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$, assuming that L_{\text{\boldmathu\unboldmath}} and $L_{p}$ are known.

The simplest separable form of the bilinear operators introduced in (4) is, for $\texttt{n}_{\texttt{pa}}=1$

[TABLE]

In general, the separation of the bilinear form may require a sum of different terms. For the sake of a simple notation, this one-term separation is assumed to hold. The general case does not introduce additional conceptual complexity.

Thus, introducing the computable scalar quantities

[TABLE]

the system of equations (9a) and (9b) reads

[TABLE]

where the residual character of the left-hand-sides of (12a) and (12b) is emphasized introducing the notations \mathcal{R}_{\text{\boldmathu\unboldmath}} and $\mathcal{R}_{p}$ such that for any \text{\boldmathw\unboldmath}\in[\mathcal{H}^{1}_{\Gamma_{D}}]^{\texttt{n}_{\texttt{sd}}} and $\omega\in\mathcal{L}_{2}(\Omega)$,

[TABLE]

Note that problem (12) is linear for \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$ and has the same structure of a standard (nonparametric) Stokes problem.

Once the discrete subspaces approximating $[\mathcal{H}^{1}_{\Gamma_{D}}]^{\texttt{n}_{\texttt{sd}}}$ and $\mathcal{L}_{2}(\Omega)$ are chosen, the functional equation (12) results in a linear system of algebraic equations. The matrix associated with the system in the Stokes model is symmetric, with $2\times 2$ blocks and a null submatrix on the diagonal, namely

[TABLE]

A necessary condition to guarantee unicity of the solution is that the kernel of the gradient matrix $\mathbf{G}$ reduces to the trivial space, that is $\ker\mathbf{G}=\{\mathbf{0}\}$, where $\ker\mathbf{G}:=\{\mathbf{q}:\mathbf{q}\in\mathbb{R}^{\hat{n}}\text{ and }\mathbf{G}\mathbf{q}=\mathbf{0}\}$, $\hat{n}$ being the number of pressure unknowns in the spatial domain. This implies that the standard finite element approaches for incompressibility can readily be applied in the context of the PGD parameterized Stokes problem. That is, the user-preferred choice of LBB spatial elements or incompressible stabilization can be directly used in this context.

Remark 4** (Divergence-free solution).**

Note that if the first term is weakly divergence free (for instance, this is trivial for homogeneous Dirichlet boundary conditions), namely b\bigl{(}\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{0},\delta\!F_{\!p}\bigr{)}=0 for all $\delta\!F_{\!p}\in\mathcal{L}_{2}(\Omega)$, every mode will be weakly divergence free and consequently, following Remark 2, in this case \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n} is weakly divergence-free.

3.2.3 The parameter iteration

Recall that this substep is made to determine the parameter functions for each mode and consists in solving (9c) & (9d) for L_{\text{\boldmathu\unboldmath}} and $L_{p}$ assuming that the spatial functions \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$ are known. For the particular case of $\texttt{n}_{\texttt{pa}}=1$, the problem is rewritten as

[TABLE]

where

[TABLE]

and \mathcal{R}_{\text{\boldmathu\unboldmath}} and $\mathcal{R}_{p}$ defined in (13) are the known separated expressions of the residuals for velocity and pressure at the previous PGD approximation (\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1}) tested now with \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}\,\delta\!L_{\text{\boldmathu\unboldmath}} and $F_{\!p}\,\delta\!L_{p}$, respectively.

Remark 5** (Algebraic nature of (15)).**

Note that equations (15) for L_{\text{\boldmathu\unboldmath}} and $L_{p}$ are integral equations that do not derive from any differential equation but from algebraic ones. This can be readily shown by realizing that weighting function \delta\!L_{\text{\boldmathu\unboldmath}} and $\delta\!L_{p}$ could be taken (in a point collocation fashion) as a set of Dirac deltas, ensuring that the algebraic equation is fulfilled at all the points included in the collocation (the expressions in (15) do not contain any derivative of the unknowns L_{\text{\boldmathu\unboldmath}} and $L_{p}$).

In many PGD implementations, the parametric modes L_{\text{\boldmathu\unboldmath}} and $L_{p}$ are represented as Finite Element (FE) functional approximations (using the nodal values as degrees of freedom and the shape functions as functional basis) and equations (15) are solved with a Galerkin approach (taking \delta\!L_{\text{\boldmathu\unboldmath}} and $\delta\!L_{p}$ equal to the shape functions), which in this case results in a standard Least Squares functional approximation. This is typically done in order to preserve in the computation of the parametric modes L_{\text{\boldmathu\unboldmath}} and $L_{p}$ the same coding structure as for the velocity and pressure modes \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$ when solving equation (12).

Thus, once the discrete subspaces are chosen, a symmetric matrix is obtained, namely

[TABLE]

Similarly as in the previous case, the condition that \ker\mathbf{M}_{\text{\boldmathu\unboldmath},p}=\{\mathbf{0}\} ensures uniqueness of the solution.

Note that, as stated in Remarks 2 and 4, the velocity modes are divergence-free (\text{\boldmath\nabla\unboldmath}\cdot\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{s}=0 for $s=1,\ldots,n-1$) and therefore the right-hand-side of (14b) is zero. This is provoking an inconsistency that is clearly demonstrated for the particular, but not at all unusual, case of using the same discrete subspace of $\mathcal{L}_{2}(\mathcal{I})$ for both L_{\text{\boldmathu\unboldmath}} and $L_{p}$. Under such an assumption, there is only one mass matrix \mathbf{M}=\mathbf{M}_{\text{\boldmathu\unboldmath},\text{\boldmathu\unboldmath}}=\mathbf{M}_{\text{\boldmathu\unboldmath},p}, symmetric and positive definite (i.e. its kernel is zero). Thus problem (14) has a unique solution and the system to solve for each parameter substep has the following structure:

[TABLE]

This is obviously leading to an inconsistent solution of L_{\text{\boldmathu\unboldmath}}=0 and therefore \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}=0.

This shows that the second approach for the separated representation is not viable.

Note that the same conclusion is reached by following a point collocation approach as described in Remark 5. Taking \delta\!L_{\text{\boldmathu\unboldmath}}=\delta\!_{\mu} and $\delta\!L_{p}=\delta\!_{\mu}$ in (14b), that is particularizing the algebraic equations for a given value of the parameter $\mu$, (14b) results in

[TABLE]

And, being \mathcal{R}_{p}\bigl{(}\text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},F_{\!p}\,\delta\!_{\mu}\bigr{)}=0 in (16b), the solution is always L_{\text{\boldmathu\unboldmath}}^{n}(\mu)=0.

Thus, this drawback associated with the option taken in Case #1 is independent of the choice of the functional description of the parametric modes and also of the approximation criterion to compute them.

3.3 Case #2: One parameter function (same for $u$ and $p$)

For $\texttt{n}_{\texttt{pa}}=1$, using a similar notation as in (7), and with $L(\mu)$ replacing both L_{\text{\boldmathu\unboldmath}}(\mu) and $L_{p}(\mu)$

[TABLE]

The approximation defined in (17) is substituted in (3) and tested in a tangent manifold. Thus, the problem becomes, find \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}\in[\mathcal{V}]^{\texttt{n}_{\texttt{sd}}}, $F_{\!p}\in\mathcal{L}_{2}(\Omega)$ and $L\in\mathcal{L}_{2}(\mathcal{I})$ such that

[TABLE]

for all $v$ and $q$ in the tangent manifold. Now, the space of unknowns and the tangent manifold have one dimension less with respect to the previous case, since L_{\text{\boldmathu\unboldmath}} and $L_{p}$ have been replaced by $L$. The corresponding expressions for the test functions are

[TABLE]

Thus, the equation corresponding to (9) is derived by replacing also \delta\!L_{\text{\boldmathu\unboldmath}} and $\delta\!L_{p}$ by $\delta\!L$,

[TABLE]

note that suppressing one unknown suppresses also one equation, because equations (9c) and (9d) have been summed up into (20d).

Again, equations (20a) and (20b) have to be solved together, due to the saddle point structure, as noted in Remark 3.

The spatial iteration described in Section 3.2.2 is similar in this case. It consists in solving (20a) and (20b) for unknowns \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$, assuming that $L$ is known. Recalling (10) and introducing a new definition for the computable scalar quantity

[TABLE]

the system of equations (20a) and (20b) reads

[TABLE]

Note that problem (22) has the same structure as problem (12) and therefore the spatial iterations are equivalent for cases #1 and #2.

The structure of the parameter iteration, however, changes considerably in case #2. It consists in solving (20d) for $L$ assuming that \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}} and $F_{\!p}$ are known. Recalling the definitions in (15), (20d) becomes a unified version of the equation in (14) (the sum of (14a) and (14b)), namely

[TABLE]

Note that this problem results in a simple system of equations with just a mass matrix, with the right-hand-side accounting for the effect of all the residuals. The system is easily solvable and provides a single parametric mode, affecting both velocity and pressure modes, with no particular restrictions.

Consequently, the alternative analyzed as case #2 appears to be viable and, as confirmed in the numerical examples shown in the next section, is the right approach to define a parametric separation of the saddle point problems.

3.4 Least-squares PGD projection and PGD compression.

Often, the PGD separated solution is post-processed with a compression algorithm based on a least squares projection in order to reduce the number of PGD modes. This is standard in the PGD practice, because the PGD terms may contain some intrinsic redundancy that is alleviated with this post-process. When compared with an SVD separation of the complete parametric solution, the redundancy is associated with the nonorthogonality of the different terms (or, conversely, the optimality of the SVD representation is associated with their orthogonality).

The projection strategy is described in [11] where its superior performance is demonstrated when compared with standard SVD (for a 2D separation) and HOSVD (High-Order SVD, for a larger number of parameters). It consists in finding with a PGD like algorithm (that is, greedy and with an alternated directions iterative scheme) the separated functions \tilde{\text{\boldmathu\unboldmath}}_{{}_{\texttt{PGD}}} and $\tilde{p}_{{}_{\texttt{PGD}}}$ that better approximate \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}} and $p_{{}_{\texttt{PGD}}}$ with a least squares criterion. Namely, \tilde{\text{\boldmathu\unboldmath}}_{{}_{\texttt{PGD}}} and $\tilde{p}_{{}_{\texttt{PGD}}}$ are sought such that

[TABLE]

for all $v$ and $q$ ranging in some suitable spaces. Equations (24) are solved with a PGD strategy, exactly as described in the previous sections for the Stokes problem and equation (3). Typically, this operation is performed selecting \tilde{\text{\boldmathu\unboldmath}}_{{}_{\texttt{PGD}}} and $\tilde{p}_{{}_{\texttt{PGD}}}$ in the same functional spaces as \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}} and $p_{{}_{\texttt{PGD}}}$. Bur here, we consider also projecting into richer functional spaces.

In the present context, this technology deserves a particular attention because it will allow addressing a concern that is naturally raised after the considerations introduced in the previous sections. A clear conclusion of the above analysis is that the PGD solution of the Stokes problem must adopt the formulation labelled as case #2 (same parametric mode for all the velocity components and the pressure), while case #1 (same parametric mode for all the velocity components and a different one for the pressure) is not viable. Moreover, the first idea announced in section 3.1 that we may label now as case #0 (all parametric modes different) was also discarded because of its implementation complexity in commercial codes.

The discarded forms (cases #0 and #1) are richer descriptions of the solution, in the sense that, for the same number of terms in the PGD sum, the number of degrees of freedom used to describe \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}} and $p_{{}_{\texttt{PGD}}}$ is much larger in cases #0 and #1 than in case #2. Roughly speaking, the number of d.o.f. describing the parametric modes is in case #1 the double of case #2 and the factor goes to $d+1$ (being $d$ the number of spatial dimensions) when compared to case #0. Considering the total amount of d.o.f., one may think that the richer descriptions would require less terms in the PGD sum, with the subsequent computational savings. Thus, the question is: does the proposed PGD formulation of the Stokes problem requires an excessive number of PGD terms?

The PGD least squares projection is used here to answer this question. Indeed, the three alternatives are admissible and viable to solve equation (24) with a PGD approach. Thus, once \text{\boldmathu\unboldmath}_{{}_{\texttt{PGD}}} and $p_{{}_{\texttt{PGD}}}$ are computed as described in case #2, \tilde{\text{\boldmathu\unboldmath}}_{{}_{\texttt{PGD}}} and $\tilde{p}_{{}_{\texttt{PGD}}}$ may be computed solving equation (24) with any of the three formulations. This will allow checking if having more d.o.f. per term in the sum results in having a shorter PGD sum (less terms). In other words, it will indicate if the restriction of having the same parametric mode for all the components of velocity and the pressure is an artificial constraint or if, on the contrary, this additional condition fits the form of the parametric solution. The numerical examples presented in the next section demonstrate that, in a pretty general situation, the case #2 formulation does not significantly increase the number of terms in the PGD sum with respect to the alternatives corresponding to cases #0 and #1.

In short, case #1 uses L_{\text{\boldmathu\unboldmath}}\neq L_{p} and case #2 L_{\text{\boldmathu\unboldmath}}=L_{p} and consequently one would expect that one term of case #1

[TABLE]

requires two terms of case #2 (for instance, taking \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{1}=\text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}, $F_{\!p}^{1}=0$, L^{1}=L_{\text{\boldmathu\unboldmath}}, \text{\boldmathF\unboldmath}_{\!\!\text{\boldmathu\unboldmath}}^{2}=0, $F_{\!p}^{2}=F_{\!p}$ and $L^{2}=L_{p}$).

In other words, we claim that the solutions of the examples analyzed in the next section are optimally represented by the case #2 option, in the sense that reducing the number of degrees of freedom per PGD term does not increase the number of PGD terms required.

Since the alternative solutions (with cases #0 and #1) of the original problem are discarded, this assertion is demonstrated by representing (using the PGD compression) the solution obtained with case #2 in the forms of cases #0 and #1. Instead of saving a significant number of PGD terms, the compression obtained is not relevant: the number of terms required is very similar. Thus, we conclude the solution fits naturally with the functional structure provided in case #2 and the fact of reducing the number of degrees of freedom per PGD term it is not introducing extra terms in the PGD solution.

4 Numerical examples

This section presents three examples, the first two are a parametric Stokes problems, with a set of parameters determining the geometry of the computational domain. The first example is a very simple backward facing step with a single geometric parameter. This example is used to discuss the different alternatives for the PGD representation (cases #1 and #2 described in sections 3.2 and 3.3) and to analyze the effect of the corresponding PGD compression techniques.

The second example considered in this section describes a Stokes flow around a NACA airfoil, where the geometry depends on four independent parameters. The analysis includes a discussion on the how the PGD compression techniques are affected by the choice of the hybrid formulation. First, the particular treatment of geometric parameters in the PGD formulation is briefly recalled.

The third example addresses the Brinkman problem with a free parameter stating the relative weight of the Stokes and Darcy models (Brinkman is seen as a combination of both). As a consequence of the conclusion of the first example and the previous section, among the PGD formulations discussed above, only the alternative #2 (see section 3.3) is considered in the second example.

4.1 Stokes flow in domains with parametric geometry

4.1.1 Accounting for geometric parameters

The strategy to deal with geometric parameters in the PGD solver was devised in [12] for Poisson problems and then combined with material parameters in [13] and [14] for heat and wave propagation problems. The fundamental idea is using a reference domain $\mathcal{T}$ and a parametric mapping to the physical domain \Omega({\text{\boldmath\mu\unboldmath}}). Thus, the (physical) problem is stated in the reference domain. The problem in the reference domain includes some fictitious parametric properties accounting for the mapping. The mapping between $\mathcal{T}$ and \Omega({\text{\boldmath\mu\unboldmath}}) is described with a coarse FE mesh, much coarser than the computational mesh because it is only required to resolve the parametric variations of the geometry.

In the case of a Stokes problem (1), the bilinear form in the left-hand-side of the weak equation reads

[TABLE]

where {\bf J}({\text{\boldmath\mu\unboldmath}}) is the Jacobian of the mapping and \hat{\text{\boldmathx\unboldmath}} are coordinates in the reference domain $\mathcal{T}$. The main goal of using this mapping is to transform the parametric dependence of the integration domain \Omega({\text{\boldmath\mu\unboldmath}}) into a parametric dependence of material-like properties, the parametric fictitious conductivity {\bf D}({\text{\boldmath\mu\unboldmath}}). In order to use PGD, the parametric dependence of {\bf D}({\text{\boldmath\mu\unboldmath}}) has to be expressed (often it has to be approximated) in a separable form. This requires a further step computed via SVD or higher-order SVD as described in detail in [13].

4.1.2 Backward facing step

A simplified backward facing step problem based on the Stokes equation, see Figure 1, is considered where the parameter $H$ is the height of the step. The evolution of the error of the PGD in terms of the number of modes is also shown in Figure 1. The error is taken as compared with a Finite Element solution for specific values of the parameter, measured in the supreme norm, that is the infinite norm, for the space dimension. Note that the error decreases with the number of terms and is larger for the pressure than for the velocity where it stagnates at a relative error of $10^{-6}$ for 40 PGD terms.

The same curve is repeated in Figure 2, this time accompanied by three others representing the errors of the compressed PGD solutions. The additional three curves correspond to Least Squares projections using the three different alternatives, as described in Section 3.4. Recall that the three alternatives correspond to cases #0, #1 and #2 with different parametric modes in the solution. The case #0 defines the larger functional space (number of parametric modes equal to number of space dimensions plus one) which includes the functional space of case #1 (number of parametric modes equal to two). Case #2 defines a smaller functional space with a single parametric mode. In consequence, the lower projection error must correspond to case #0, then the error of case #1 is larger (or equal) and even larger (or equal) for case #2.

The results demonstrate that the three projection behave almost equally. The inclusion of the functional spaces is indeed translated into the expected inequality of the errors (error #0 $\leq$ error #1 $\leq$ error #2 ) but the difference is very small. This reveals that selecting option #2 is not practically increasing the required number of PGD terms with respect to a richer functional description (being case #0 the richest). The structure of the parametric description enforced in case #2 seems to fit the nature of the actual solution and its parametric dependence.

4.1.3 Stokes flow around a NACA airfoil

The flow around a 4-digit NACA airfoil is considered, see Figure 3. The four geometric parameters (digits, in this contex) defining the geometry of the airfoil are: 1) length, $c$, 2) thickness, $t$, 3) max camber, $m$ and 4) max camber position, $p$.

The PGD solution provides a computational vademecum containing the flow solutions for any possible NACA-4 geometry. For the sake of illustration, the evaluation for four particular values of the parameters is shown in Figure 4. These four particular parametric values are also used to check convergence with the number of PGD modes. Figure 5 presents the evolution of the errors in velocity and pressure as the number of modes increases.

Note that also in this complex example (with four parametric dimensions) the PGD solution behaves correctly when compared to standard FE solutions for specific (and representative) values of the parameters. None of this four sampling points of the 4D parametric space is a grid point (the values of the parameters selected do not coincide with the discrete grid of each parametric dimension). Note that integration in the parametric space to compute a $L_{2}$ norm of the error is avoided because of the associated computational burden.

4.2 Flow in fractured media

Describing flow in fractured, vuggy and porous media with a unique model is important for reservoir engineering. A standard approach is using Darcy’s law in the porous domain and Stokes law in the parts of the domain containing voids and fractures. The Brinkman model merges both Darcy and Stokes, see [16, 17, 18], linearly combining the effect of the two constitutive models. Thus, the Brinkman problem reads: find velocity $u$ and pore pressure $p$ such that

[TABLE]

where $\eta$ and $\tilde{\eta}$ are the dynamic and effective viscosity respectively and $\mathbf{K}$ is the permeability tensor. Note that the dynamic viscosity $\eta$ allows tuning the weight of the Darcy model in the constitutive description of the flow. For small values of $\eta$, a free-flow Stokes-like pattern is obtained and for large values of $\eta$ the solution tends to behave as porous flow (Darcy). In reservoir modelling it is usual to assume incompressibility (25b), to neglect gravity (right-hand-side of (25a)) and to set $\eta=\tilde{\eta}$ for the bulk material, [17, 18].

A parametrized version of the Brinkman problem is built upon the background permeability field provided as a test case (SPE10) by the Society of Petroleum Engineers [15]. The original three-dimensional (3D) field is restricted to a 2D domain following Könnö and Stenberg [18]. Here, we aim at analyzing the effect of adding a fracture to the layer 68 of the SPE10 model. Note that SPE10 provides a non-uniform isotropic permeability and, therefore, the matrix $\mathbf{K}$ becomes a scalar field k_{\text{SPE10}}(\text{\boldmathx\unboldmath}). In the framework of the Brinkman model, the fracture is accounted for in a natural way by significantly increasing the permeability in the region where the fracture is located. Figure 6 shows the SPE10 permeability with and without the perturbation that accounts for the fracture.

The value of the perturbed permeability in the fractured zone is not easy to set. The modeller knows only that permeability has to be significantly larger in the zones where the free-flow pattern is expected, with respect to the reference values for porous flow. Thus, it is particularly interesting having a tuning parameter that allows enforcing a gradual transition between the two regimes. This is especially relevant in flow simulations of karst reservoirs where vugs and caves are embedded in a porous rock and are connected via fracture networks at multiple scales. For example the work of Popov and coauthors [17] based on a Brinkman model, assumes a continuous permeability ranging six orders of magnitude and investigates the effect of permeability in a filled fracture.

Here, the model of Könnö and Stenberg is extended including one additional parameter that controls the permeability of the fracture, occupying a subdomain $\Lambda$ in $\Omega$ . This parameter is denoted by $\mu$ and can be interpreted as a measure of the degree of saturation of some filling inside the fracture. The influence of the input parameter $\mu$ ranging in $[0,\mu_{\text{max}}]$ in the resulting permeability distribution is given by the following expression:

[TABLE]

where n(\text{\boldmathx\unboldmath}):=\log_{10}\left(k_{\text{SPE10}}(\text{\boldmathx\unboldmath})\right). Note that for the extreme values of $\mu$, k(\text{\boldmathx\unboldmath},0)=k_{\text{SPE10}}(\text{\boldmathx\unboldmath}) and k(\text{\boldmathx\unboldmath},\mu_{\text{max}})=10^{\mu_{\text{max}}}.

In order to implement the PGD, input data must be expressed in a separated form. The only term in (25) that it is not trivially separable is the second term of (25a), involving the inverse of the permeability k(\text{\boldmathx\unboldmath},\mu). Thus, the inverse of k(\text{\boldmathx\unboldmath},\mu) as defined in (26) has to be expressed as a separated expression, that is in terms of functions that depend only on n(\text{\boldmathx\unboldmath}) and functions that depend only on $\mu$. Note that this is only needed for \text{\boldmathx\unboldmath}\in\Lambda. The part of k(\text{\boldmathx\unboldmath},\mu)^{-1} which is not trivially separable is 10^{\frac{\mu}{\mu_{\text{max}}}n(\text{\boldmathx\unboldmath})}. A SVD of a dense sampling of this function is used to separate it. The separated approximation is a description in terms of modes \Theta^{m}(n(\text{\boldmathx\unboldmath})) and $\phi^{m}(\mu)$, $m=1,\ldots,M$, and reads

[TABLE]

Thus, for \text{\boldmathx\unboldmath}\in\Lambda,

[TABLE]

Note that the separation is performed in terms of the variables n(\text{\boldmathx\unboldmath}) and $\mu$ (instead of $x$ and $\mu$). This simplifies the function to be separated and reduces the number of terms required to reach some prescribed accuracy. In this case, using $M=16$ terms provides a relative error smaller than $10^{-12}$ (for any value of n(\text{\boldmathx\unboldmath}) and $\mu$).

The PGD solution is sought with the form defined in (17) (same parameter functions for velocity and pressure). The evolution with the number of PGD terms of the error with respect to a standard FE solution is shown in Figure 7. It can be observed that the relative error is of order $10^{-5}$ with only 40 PGD terms.

Despite being relatively simple problem (with only one scalar parameter $\mu$), the convergence is faster than in other PGD solutions. The resulting velocity and pressure fields for the extreme values of $\mu$ are shown in Figure 8, and the errors with respect to the corresponding finite element solutions are in Figure 9. Note that for these two particular values, the errors are lower than for the worst case scenario depicted in the convergence curve of Figure 7.

5 Conclusions

The analysis of the different forms for the parametric separation of Stokes problems reveals that the only viable option is having a unique parametric mode for each independent parameter, affecting all the velocity components and the pressure. This choice (denoted as case #2) has less degrees of freedom than the alternative cases #1 (one parametric mode for all velocity components and a different one for the pressure) and #0 (different parametric modes for every velocity component and the pressure). The proposed formulation corresponding to case #2 is the simplest alternative, guarantees incompressibility and is not affected by any stability concerns.

The abundance of degrees of freedom may be a desirable feature of the PGD formulation because a richer functional space could compensate the computational overhead for each mode with a reduced number of modes. However, in this case increasing the unknowns leads to unsolvable problems and therefore the alternatives #0 and #1 have to be discarded.

Moreover, the PGD compression algorithm based on a Least Squares projection can be performed for the three alternative parametric representations. The analysis of the three compressions demonstrates that the parametric structure of the solutions does not require the multiplicity of parametric modes: the number of modes required to reach some prescribed accuracy for case #2 is almost the same as for the other two alternatives.

Acknowledgements

This work has been supported by Ministerio de Economía y Competitividad, grant number DPI2014-51844-C2-2-R and by Generalitat de Catalunya, grant number 2014-SGR-1471.