A geometric multigrid method for isogeometric compatible discretizations of the generalized Stokes and Oseen problems

TL;DR

This paper introduces a geometric multigrid method for efficiently solving isogeometric discretizations of the generalized Stokes and Oseen problems, ensuring divergence-free velocity fields and robustness across grid resolutions.

Contribution

The paper develops a multigrid solver with Schwarz-style smoothers that guarantees divergence-free solutions and maintains convergence rates regardless of grid size or flow parameters.

Findings

Convergence rates are invariant to grid resolution.

Method produces divergence-free velocity fields.

Robust for non-advection-dominated flows.

Abstract

In this paper, we present a geometric multigrid methodology for the solution of matrix systems associated with isogeometric compatible discretizations of the generalized Stokes and Oseen problems. The methodology provably yields a pointwise divergence-free velocity field independent of the number of pre-smoothing steps, post-smoothing steps, grid levels, or cycles in a V-cycle implementation. The methodology relies upon Scwharz-style smoothers in conjunction with specially defined overlapping subdomains that respect the underlying topological structure of the generalized Stokes and Oseen problems. Numerical results in both two- and three-dimensions demonstrate the robustness of the methodology through the invariance of convergence rates with respect to grid resolution and flow parameters for the generalized Stokes problem as well as the generalized Oseen problem provided it is not advection-dominated.

Tables

Table 1. Table 1 : Number of V(1,2) cycles required for convergence for the unit square generalized Stokes problem using the multiplicative Schwarz smoother and p = 2 , 3 𝑝 2 3 p=2,3 .

	$p=2$			$p=3$
$n_{\ell }$	DOFs	$Da=1$	$Da=1000$	DOFs	$Da=1$	$Da=1000$
1	16	3	2	25	6	6
2	36	5	4	49	12	14
3	100	5	5	121	13	13
4	324	6	5	361	12	12
5	1156	6	5	1225	13	13
6	4356	7	6	4489	13	13
7	16900	7	6	17161	13	13
8	66564	7	6	67081	13	13
9	264196	7	7	265225	13	13
10	1052676	7	7	1054729	13	13

Table 2. Table 2 : Number of V(1,2) cycles required for convergence for the unit square generalized Stokes problem using the additive Schwarz smoother and p = 2 𝑝 2 p=2 .

$n_{\ell }$	DOFs	$Da=1$	$Da=1000$
1	16	5	4
2	36	12	9
3	100	15	14
4	324	16	15
5	1156	16	16
6	4356	16	16
7	16900	16	16
8	66564	16	16
9	264196	16	16
10	1052676	16	16

Table 3. Table 3 : Number of V(1,2) cycles required for convergence for the quarter annulus generalized Stokes problem using the multiplicative Schwarz smoother and p = 2 𝑝 2 p=2 .

$n_{\ell }$	DOFs	$Da=1$	$Da=1000$
1	16	2	2
2	36	3	3
3	100	8	7
4	324	16	15
5	1156	24	24
6	4356	26	26
7	16900	31	31
8	66564	34	34
9	264196	34	34
10	1052676	34	34

Table 4. Table 4 : Number of V(1,2) cycles required for convergence for the unit cube generalized Stokes problem using the multiplicative Schwarz smoother and p = 2 , 3 𝑝 2 3 p=2,3 .

	$p=2$			$p=3$
$n_{\ell }$	DOFs	$Da=1$	$Da=1000$	DOFs	$Da=1$	$Da=1000$
1	144	1	1	300	2	4
2	540	1	1	882	1	1
3	2700	2	2	3630	1	1
4	16524	1	1	19494	1	1
5	114444	1	1	124950	1	1
6	849420	1	1	888822	1	1

Table 5. Table 5 : Number of V(1,2) cycles required for convergence for the hollow cylinder generalized Stokes problem using the multiplicative Schwarz smoother and p = 2 𝑝 2 p=2 .

$n_{\ell }$	DOFs	$Da=1$	$Da=1000$
1	144	1	1
2	540	1	2
3	2700	2	2
4	16524	1	1
5	114444	1	1
6	849420	1	1

Table 6. Table 6 : Number of V(1,2) cycles required for convergence for the unit square generalized Oseen problem using the multiplicative Schwarz smoother and p = 2 𝑝 2 p=2 .

		$Re=1$		$Re=100$
$n_{\ell }$	DOFs	$Da=1$	$Da=1000$	$Da=1000$
1	16	3	2	2
2	36	5	4	3
3	100	5	5	4
4	324	6	5	9
5	1156	6	5	13
6	4356	7	6	15
7	16900	7	6	11
8	66564	7	6	7
9	264196	7	7	7
10	1052676	7	7	7

Table 7. Table 7 : Number of V(1,2) cycles required for convergence for the unit cube generalized Oseen problem using the multiplicative Schwarz smoother and p = 2 𝑝 2 p=2 .

		$Re=1$		$Re=100$
$n_{\ell }$	DOFs	$Da=1$	$Da=1000$	$Da=1000$
1	144	1	2	2
2	540	1	2	2
3	2700	2	3	3
4	16524	1	2	2
5	114444	1	1	3
6	849420	1	1	3

Equations

$$\displaystyle\begin{array}[]{rll}\frac{\partial{\bf u}}{\partial t}+{\bf u}\cdot\nabla{\bf u}-\nu\Delta{\bf u}+\nabla p&={\bf f}&\textup{ for }({\bf x},t)\in\Omega\times(0,T)\\ \nabla\cdot{\bf u}&=0&\textup{ for }({\bf x},t)\in\Omega\times(0,T)\\ {\bf u}&={\bf 0}&\textup{ for }({\bf x},t)\in\Gamma\times(0,T)\\ \left.{\bf u}\right|_{t=0}&={\bf u}_{0}&\textup{ for }{\bf x}\in\Omega\end{array}$$

$$\frac{\partial \mathbf{u}}{\partial t}(t_{n+1/2})$$

$$\Delta \mathbf{u}(t_{n+1/2})$$

$$\nabla p(t_{n+1/2})$$

$$\left(\mathbf{u}\cdot \nabla \mathbf{u}\right)(t_{n+1/2})$$

$$\displaystyle\begin{array}[]{rll}\sigma{\bf u}^{(n+1)}-\nu\Delta{\bf u}^{(n+1)}+\nabla p^{(n+1)}&={\bf f}_{GS}^{(n+1)}&\textup{ for }{\bf x}\in\Omega\\ \nabla\cdot{\bf u}^{(n+1)}&=0&\textup{ for }{\bf x}\in\Omega\\ {\bf u}^{(n+1)}&={\bf 0}&\textup{ for }{\bf x}\in\Gamma\end{array}$$

$${\mathbf{f}}_{GS}^{(n+1)}=\mathbf{f}(t_{n+1/2})+\sigma {\mathbf{u}}^{(n)}-3{\mathbf{u}}^{(n)}\cdot \nabla {\mathbf{u}}^{(n)}+{\mathbf{u}}^{(n-1)}\cdot \nabla {\mathbf{u}}^{(n-1)}+\nu \Delta {\mathbf{u}}^{(n)}-\nabla p^{(n)}$$

$$\left(\mathbf{u}\cdot \nabla \mathbf{u}\right)(t_{n+1/2})$$

$$\displaystyle\begin{array}[]{rll}\sigma{\bf u}^{(n+1)}+{\bf a}^{(n+1)}\cdot\nabla{\bf u}^{(n+1)}-\nu\Delta{\bf u}^{(n+1)}+\nabla p^{(n+1)}&={\bf f}_{GO}^{(n+1/2)}&\textup{ for }{\bf x}\in\Omega\\ \nabla\cdot{\bf u}^{(n+1)}&=0&\textup{ for }{\bf x}\in\Omega\\ {\bf u}^{(n+1)}&={\bf 0}&\textup{ for }{\bf x}\in\Gamma\end{array}$$

$${\mathbf{f}}_{GO}^{(n+1)}=\mathbf{f}(t_{n+1/2})+\sigma {\mathbf{u}}^{(n)}-{\mathbf{a}}^{(n+1)}\cdot \nabla {\mathbf{u}}^{(n)}+\nu \Delta {\mathbf{u}}^{(n)}-\nabla p^{(n)}$$

$${\mathbf{H}}_0^1(\Omega ):=\left\{\mathbf{v}\in {\mathbf{H}}^1(\Omega ):\mathbf{v}=0\text{ on }\Gamma \right\}$$

$$L_0^2(\Omega ):=\left\{q\in L^2(\Omega ):{\int }_{\Omega }qd\Omega =0\right\}$$

$$a(\mathbf{v},\mathbf{u})-b(\mathbf{v},p)+b(\mathbf{u},q)=\ell (\mathbf{v})$$

$$a({\bf v},{\bf u}):=\left\{\begin{array}[]{cl}\displaystyle\int_{\Omega}\sigma{\bf v}\cdot{\bf u}\ d\Omega+\int_{\Omega}\nu\nabla{\bf v}:\nabla{\bf u}\ d\Omega&\text{generalized Stokes}\\ \displaystyle\int_{\Omega}\sigma{\bf v}\cdot{\bf u}\ d\Omega+\int_{\Omega}{\bf v}\cdot\left({\bf a}\cdot\nabla{\bf u}\right)\ d\Omega+\int_{\Omega}\nu\nabla{\bf v}:\nabla{\bf u}\ d\Omega&\text{generalized Oseen}\\ \end{array}\right.$$

$$b(\mathbf{v},p):={\int }_{\Omega }\left(\nabla \cdot \mathbf{v}\right)pd\Omega$$

$$\ell (\mathbf{v}):={\int }_{\Omega }\mathbf{v}\cdot \mathbf{f}d\Omega$$

$$a({\mathbf{v}}_h,{\mathbf{u}}_h)-b({\mathbf{v}}_h,p_h)+b({\mathbf{u}}_h,q_h)=\ell ({\mathbf{v}}_h)$$

$$a_h({\mathbf{v}}_h,{\mathbf{u}}_h)-b({\mathbf{v}}_h,p_h)+b({\mathbf{u}}_h,q_h)=\ell ({\mathbf{v}}_h)$$

$$a_h({\mathbf{v}}_h,{\mathbf{u}}_h):=a({\mathbf{v}}_h,{\mathbf{u}}_h)-{\int }_{\Gamma }\nu {\mathbf{v}}_h\cdot {\nabla }_{\mathbf{n}}{\mathbf{u}}_{h}d\Gamma -{\int }_{\Gamma }\nu {\nabla }_{\mathbf{n}}{\mathbf{v}}_h\cdot {\mathbf{u}}_{h}d\Gamma +{\int }_{\Gamma }\frac{C_I\nu }{h}{\mathbf{v}}_h\cdot {\mathbf{u}}_{h}d\Gamma$$

$$\left[\begin{array}[]{cc}{\bf A}&-{\bf B}\\ {\bf B}^{T}&{\bf 0}\end{array}\right]\left(\begin{array}[]{c}{\bf u}\\ {\bf p}\end{array}\right)=\left(\begin{array}[]{c}{\bf f}\\ {\bf 0}\end{array}\right)$$

$$[\mathbf{A}{]}_{ij}$$

$$[\mathbf{B}{]}_{ij}$$

$$[\mathbf{f}{]}_i$$

$$\mathbf{K}\mathbf{U}=\mathbf{F}$$

$$\text{setcounter}MaxMatrixCols11\begin{array}{ccccccccccc}0& \to & \Phi & \stackrel{\vec{\nabla }}{\to }& \mathbf{\Psi }& \stackrel{\vec{\nabla }\times }{\to }& \mathbf{V}& \stackrel{\vec{\nabla }\cdot }{\to }& \text{Q}& \to & 0\end{array}$$

$$\Phi$$

$$\mathbf{V}$$

$$\text{setcounter}MaxMatrixCols11\begin{array}{ccccccccccc}0& \to & \Phi & \stackrel{\vec{\nabla }}{\to }& \mathbf{\Psi }& \stackrel{\vec{\nabla }\times }{\to }& \mathbf{V}& \stackrel{\vec{\nabla }\cdot }{\to }& \text{Q}& \to & 0\\ \downarrow {\Pi }_{\varphi }& & \downarrow {\Pi }_{\psi }& & \downarrow {\Pi }_v& & \downarrow {\Pi }_q& & & & \\ 0& \to & {\Phi }_h& \stackrel{\vec{\nabla }}{\to }& {\mathbf{\Psi }}_h& \stackrel{\vec{\nabla }\times }{\to }& {\mathbf{V}}_h& \stackrel{\vec{\nabla }\cdot }{\to }& {\text{Q}}_h& \to & 0\end{array}$$

$$\displaystyle\begin{array}[]{rl}\hat{N}_{i,p}(\xi)&:=\frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}\hat{N}_{i,p-1}(\xi)+\frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}\hat{N}_{i+1,p-1}(\xi)\hskip 5.0pt\textup{ for }p>0\\ \\ \hat{N}_{i,0}(\xi)&:=\left\{\begin{array}[]{rl}1&\xi_{i}\leq\xi<\xi_{i+1}\\ 0&\text{elsewhere}\end{array}\right.\end{array}$$

$${\hat{N}}_{\mathbf{i},\mathbf{p}}(\mathbf{\xi }):=\prod _{k=1}^d{\hat{N}}_{i_k,p_k}({\xi }_k)$$

$$S^{p_{1},p_{2},\ldots,p_{d}}_{\bm{\alpha}_{1},\bm{\alpha}_{2},\ldots,\bm{\alpha}_{d}}\left(\mathcal{M}_{h}\right):=\left\{f:\hat{\Omega}\rightarrow\mathbb{R}\ \Big{|}\ f(\bm{\xi})=\sum_{\bf i}a_{\bf i}\hat{N}_{{\bf i},{\bf p}}(\bm{\xi})\right\},$$

$${\hat{\Psi }}_h$$

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

$${\hat{\text{Q}}}_h$$

Full text

Keywords: Geometric multigrid; Isogeometric compatible discretizations; Isogeometric divergence-conforming discretizations; Generalized Stokes flow; Generalized Oseen flow; Overlapping Schwarz smoothers

1 Introduction

Isogeometric compatible discretizations444Depending on context, isogeometric compatible discretizations may also be referred to as isogeometric discrete differential forms, structure-preserving discretizations, divergence-conforming discretizations, or curl-conforming discretizations. have recently arisen as an attractive candidate for the spatial discretization of fluid flow problems [1, 2, 3, 4, 5, 6, 7]. These discretizations comprise a discrete Stokes complex [8, 9] and may be interpreted as smooth generalizations of Raviart-Thomas-Nédélec finite elements [10, 11]. When applied to incompressible flow problems, isogeometric compatible discretizations produce pointwise divergence-free velocity fields and hence exactly satisfy mass conservation. As a result, they preserve the balance law structure of the incompressible Navier-Stokes equations, and in particular, they properly conserve mass, linear and angular momentum, energy, vorticity, enstrophy (in the two-dimensional setting), and helicity (in the three-dimensional setting) in the inviscid limit [5]. Isogeometric compatible discretizations have recently been applied to Cahn-Hilliard flow [12], turbulent flow [13], and fluid-structure interaction [14] where improved results were attained in comparison with state-of-the-art discretization procedures.

Despite the promise of isogeometric compatible discretizations, very little research has been conducted so far in the area of efficient linear solvers. In fact, only the performance of Krylov subspace methods in conjunction with block preconditioners has been investigated in prior work [15, 16]. The objective of the current work is to introduce an optimally efficient linear solution procedure for isogeometric compatible discretizations of the generalized Stokes and Oseen problems. It should be noted that there are many different candidates in this regard. For instance, there exist efficient physics-based splitting methods such as the inexact Uzawa algorithm [17]. However, these techniques rely on suitable Schur complement approximations which can be difficult to design in the context of generalized Oseen flow. Alternatively, one can employ a multigrid method in conjunction with a Vanka smoother [18], a Uzawa smoother [19], or a Braess-Sarazin smoother [20]. While these techniques generally do not require accurate Schur complement approximations, they typically involve specially tuned relaxation parameters. Perhaps more concerning is the fact that all of the aforementioned procedures do not return a pointwise-divergence free velocity field unless the linear solver is fully converged.

To overcome the issues associated with the aforementioned linear solution procedures, we present a geometric multigrid methodology which relies upon Schwarz-style smoothers [21] in conjunction with specially defined overlapping subdomains that respect the underlying topological structure of the generalized Stokes and Oseen problems. This methodology is inspired by multigrid and auxiliary space preconditioning methodologies for divergence-conforming discontinuous Galerkin formulations of Stokes flow [22, 23, 24, 25] and multigrid methodologies for compatible finite element discretizations of Darcy and Maxwell problems [26, 27]. We prove that our methodology yields a pointwise divergence-free velocity field independent of the number of pre-smoothing steps, post-smoothing steps, grid levels, or cycles in a V-cycle implementation. We also demonstrate by numerical example that our methodology is optimally efficient and robust in that it exhibits convergence rates independent of the grid resolution and flow parameters for the generalized Stokes problem as well as the generalized Oseen problem provided it is not advection-dominated. It should be mentioned that the only user-defined constants in our methodology are the number of pre-smoothing steps and post-smoothing steps as well as the scaling factor if one elects to use an additive Schwarz smoother rather than a multiplicative Schwarz smoother. However, we have found that our method is optimally efficient irregardless of the number of pre- and post-smoothing steps selected.

An outline of the remainder of the paper is as follows. In Section 2, we inspire the need for efficient linear solvers for the generalized Stokes and Oseen problems through a discussion of temporal discretization of the Navier-Stokes equations. In Section 3, we discuss spatial discretization of the generalized Stokes and Oseen problems. In Section 4, we introduce the Stokes complex and demonstrate how to construct isogeometric compatible discretizations which commute with this complex. In Section 5, we present our structure-preserving geometric multigrid methodology, and we prove that this methodology indeed yields discrete velocity fields which are divergence-free. In Section 6, we apply the proposed multigrid method to a selection of generalized Stokes and Oseen problems. Finally, in Section 7, we provide concluding remarks.

2 Temporal Discretization of the Navier-Stokes Equations and the Generalized Stokes and Oseen Problems

To motivate the need for efficient linear solvers for the generalized Stokes and Oseen problems, we first demonstrate how such problems arise through semi-implicit temporal discretization of the incompressible Navier-Stokes equations subject to homogeneous Dirichlet boundary conditions. For $d\in\mathbb{Z}_{+}$, let $\Omega\subset\mathbb{R}^{d}$ denote an open, Lipschitz bounded domain, let $\Gamma$ denote the boundary of $\Omega$, and let $T\in\mathbb{R}_{+}$. Given $\nu\in\mathbb{R}_{+}$, ${\bf f}:\Omega\times(0,T)\rightarrow\mathbb{R}^{d}$, and ${\bf u}_{0}:\Omega\rightarrow\mathbb{R}^{d}$, the strong form of the Navier-Stokes problem then reads as follows: Find ${\bf u}:\overline{\Omega}\times[0,T]\rightarrow\mathbb{R}^{d}$ and $p:\Omega\times(0,T)\rightarrow\mathbb{R}$ such that:

[TABLE]

Above, ${\bf u}$ denotes the velocity field, $p$ denotes the pressure field, $\nu$ denotes the kinematic viscosity, ${\bf f}$ denotes the force per unit mass, and ${\bf u}_{0}$ denotes the initial velocity field. The velocity field is uniquely specified by the Navier-Stokes problem while the pressure field is unique up to a constant. To discretize in time, we first define a sequence of time instances $t_{0}<t_{1}<t_{2}<\ldots<t_{N}$ such that $t_{0}=0$ and $t_{N}=T$, and we denote the velocity and pressure solutions at the $n^{\textup{th}}$ time instance as ${\bf u}^{(n)}$ and $p^{(n)}$ respectively for $n=0,\ldots,N$. We further define $t_{n+1/2}=\frac{t_{n}+t_{n+1}}{2}$. Without loss of generality, we assume that the time instances are equi-spaced, and we define the time step size to be $\Delta t=t_{n+1}-t_{n}$. The velocity solution at $N=0$ is given by ${\bf u}^{(0)}={\bf u}_{0}$, while to find the velocity and pressure solutions at each subsequent time instance, we must discretize the Navier-Stokes problem in time. We discuss two demonstrative semi-implicit temporal discretization schemes herein, though the proceeding discussion also applies to other semi-implicit temporal discretization schemes555With a fully implicit time discretization scheme, one must turn to a nonlinear solution procedure such as Newton’s method. However, with Newton’s method, one solves a sequence of generalized Oseen problems. Hence, the multigrid methodology discussed here can also be employed to solve these problems.

Let us first consider the standard Crank-Nicolson/Adams-Bashforth scheme [28]. In this approach, a central difference approximation of the unsteady term and linear interpolation approximations of the diffusive and pressure force terms are employed:

[TABLE]

The advection term at $t_{n+1/2}$ is alternatively approximated using Taylor-series expansions involving time instances $t_{n-1}$ and $t_{n}$, resulting in666This approximation is not properly defined for $n$ = 0. Consequently, the approximation is replaced by ${\bf u}^{(n)}\cdot\nabla{\bf u}^{(n)}$ for $n$ = 0 in practice.:

[TABLE]

Collecting the above approximations, we find that the resulting generalized Stokes system holds for each $n=0,\ldots,N-1$:

[TABLE]

where $\sigma=\sigma=\nicefrac{{2}}{{\Delta t}}$ and:

[TABLE]

The above generalized Stokes system is reaction-dominated for small time step sizes and diffusion-dominated for large time step sizes. We demonstrate later that our geometric multigrid methodology is robust for both of these regimes.

The advantage of the Crank-Nicolson/Adams-Bashforth scheme is that the advection term is handled in a purely explicit manner. After spatial discretization, this leads to a symmetric matrix problem. However, the disadvantage of the scheme is that it is stable only if the time-step is chosen sufficiently small as to satisfy a CFL condition. With this in mind, we next consider an unconditionally stable semi-implicit scheme introduced by Guermond [29]. In this scheme, the unsteady, diffusive, and pressure force terms are approximated as before, but the advection term is approximated as follows777This approximation is also not properly defined for $n$ = 0. Consequently, the approximation is replaced by ${\bf u}^{(n)}\cdot\nabla\left(\left({\bf u}^{(n)}+{\bf u}^{(n+1)}\right)/2\right)$ for $n$ = 0 in practice.:

[TABLE]

Note that the advection velocity is approximated in an explicit manner while the gradient is approximated in an implicit manner. Collecting the above approximations, we find that the resulting generalized Oseen system holds for each $n=0,\ldots,N-1$:

[TABLE]

where $\sigma=\nicefrac{{2}}{{\Delta t}}$, ${\bf a}^{(n+1)}=\frac{3}{2}{\bf u}^{(n)}-\frac{1}{2}{\bf u}^{(n-1)}$, and:

[TABLE]

In opposition with the generalized Stokes system obtained earlier, the above system admits different behavior based on not only the scalars $\sigma$ and $\nu$ but also the advection velocity ${\bf a}^{(n+1)}$. We demonstrate later that our geometric multigrid methodology is robust for this system provided it is not advection-dominated. This holds if a CFL-like condition is satisfied.

3 Spatial Discretization of the Generalized Stokes and Oseen Problems

Now that we have motivated the need for efficient linear solvers for the generalized Stokes and Oseen problems, we turn to the question of spatial discretization. In this section, we present the basic ingredients associated with a mixed Galerkin discretization. Later, we will specialize to the setting of isogeometric compatible discretizations.

3.1 Weak Formulation of the Generalized Stokes and Oseen Problems

To begin, we must state a weak formulation for the generalized Stokes and Oseen Problems. We strictly consider the case of homogeneous Dirichlet boundary conditions without loss of generality. Before proceeding, we must first define suitable velocity and pressure test spaces:

[TABLE]

We also assume that $\sigma,\nu\in\mathbb{R}_{+}$, ${\bf a}\in{\bf H}^{1}_{0}(\Omega)$, and ${\bf f}\in{\bf L}^{2}(\Omega)$, and we assume that the advection velocity is divergence-free, that is, $\nabla\cdot{\bf a}\equiv 0$. With these assumptions in hand, the weak form of the generalized Stokes or Oseen problem is stated as follows: Find ${\bf u}\in{\bf H}^{1}_{0}(\Omega)$ and $p\in L^{2}_{0}(\Omega)$ such that:

[TABLE]

for all ${\bf v}\in{\bf H}^{1}_{0}(\Omega)$ and $q\in L^{2}_{0}(\Omega)$ where:

[TABLE]

[TABLE]

[TABLE]

3.2 Mixed Galerkin Approximation of the Generalized Stokes and Oseen Problems

To discretize in space using a mixed Galerkin formulation, we first must specify finite-dimensional approximation spaces for the velocity and pressure fields. We denote these spaces as ${\bf V}_{h}\subset{\bf H}^{1}_{0}(\Omega)$ and $\text{Q}_{h}\subset L^{2}_{0}(\Omega)$ respectively, but we defer the discussion of suitable approximation spaces to Section 4. With approximation spaces defined, the mixed Galerkin formulation of the generalized Stokes or Oseen problem is stated as follows: Find ${\bf u}_{h}\in{\bf V}_{h}$ and $p_{h}\in\text{Q}_{h}$ such that

[TABLE]

for all ${\bf v}_{h}\in{\bf V}_{h}$ and $q\in\text{Q}_{h}$. It should be noted that the velocity and pressure approximation spaces may not be arbitrarily selected. Instead, they should be chosen such that the Babuška-Brezzi inf-sup condition is satisfied [30, 31]. We later select isogeometric divergence-conforming discretizations for spatial discretization which indeed satisfy such a condition.

3.3 Weak Enforcement of No-Slip Boundary Conditions

The no-slip boundary condition, ${\bf u}\times{\bf n}={\bf 0}$ where ${\bf n}$ is the outward facing normal to $\Gamma$, leads to the formation of boundary layers for wall-bounded flows. High mesh resolution is required near boundary layers to accurately represent associated sharp layers, so when the no-slip condition is strongly enforced in a mixed Galerkin formulation, inaccurate flow field approximations are obtained for insufficiently-resolved boundary layer meshes. It has recently been shown that superior results can be achieved by imposing the no-penetration boundary condition strongly and the no-slip boundary condition weakly using a combination of upwinding and Nitsche’s method [32, 33]. With such an approach, we first specify finite-dimensional velocity and pressure approximation spaces as before, but we only require that the corresponding discrete velocity fields satisfy ${\bf v}\cdot{\bf n}=0$. That is, we specify ${\bf V}_{h}\subset{\bf H}^{1}_{n}(\Omega)=\left\{{\bf v}\in{\bf H}^{1}(\Omega):{\bf v}\cdot{\bf n}=0\textup{ on }\Gamma\right\}$ and $\text{Q}_{h}\subset L^{2}_{0}(\Omega)$. The corresponding formulation for the generalized Stokes or Oseen problem is then stated as: Find ${\bf u}_{h}\in{\bf V}_{h}$ and $p_{h}\in\text{Q}_{h}$ such that:

[TABLE]

for all ${\bf v}_{h}\in{\bf V}_{h}$ and $q_{h}\in\text{Q}_{h}$ where:

[TABLE]

Above, $h$ is the wall-normal element mesh size and $C_{I}$ is a positive constant that must be chosen sufficiently large to ensure coercivity of the bilinear form $a_{h}(\cdot,\cdot)$. Appropriate values for the constant $C_{I}$ can be obtained by solving element-wise eigenvalue problems or by appealing to analytical upper bounds for the trace inequality [34]. We choose to weakly enforce no-slip boundary conditions throughout the remainder of this work. This not only leads to more accurate numerical results, but it also ensures proper solution behavior in the limit of zero viscosity [4, 5].

3.4 The Matrix Problem

The formulation given by (18) yields a linear matrix system when the discrete velocity and pressure spaces are provided basis functions. Let $\left\{{\bf N}^{v}_{i}\right\}_{i=1}^{n_{v}}$ denote a set of vector basis functions for ${\bf V}_{h}$ where $n_{v}=\text{dim}\left({\bf V}_{h}\right)$, and let $\left\{N^{q}_{i}\right\}_{i=1}^{n_{q}}$ denote a set of scalar basis functions for $\text{Q}_{h}$ where $n_{q}=\text{dim}\left(\text{Q}_{h}\right)$. Then the resulting matrix system takes the form:

[TABLE]

where:

[TABLE]

Moreover, this matrix system can be written concisely as:

[TABLE]

where the matrix ${\bf K}$ has the block structure in (19) and the vectors ${\bf U}$ and ${\bf F}$ are vectors representing the group variables in (19).

4 The Stokes Complex and Isogeometric Compatible Discretizations

It remains to specify suitable velocity and pressure approximation spaces for the generalized Stokes and Brinkman problems. In this section, we present a particular selection of velocity and pressure approximation spaces which is not only inf-sup stable but also yields pointwise divergence-free discrete velocity fields. Before doing so, however, we first introduce the so-called Stokes complex which succinctly captures the fundamental theorem of calculus and expresses the differential relationships between potential, velocity, and pressure fields.

4.1 The Stokes Complex

The Stokes complex [8, 9] is a cochain complex of the form:

[TABLE]

in the three-dimensional setting where:

[TABLE]

are infinite-dimensional spaces of scalar potential fields, vector potential fields, velocity fields, and pressure fields. The Stokes complex is a smoothed version of the classical $L^{2}$ de Rham complex, and when the domain $\Omega\subset\mathbb{R}^{3}$ is simply connected with simply connected boundary, the Stokes complex is exact. This means that every pressure field may be represented as the divergence of a velocity field, every divergence-free velocity field may be represented as the curl of a vector potential field, and every curl-free vector potential field may be represented as the gradient of a scalar potential field. An analogous two-dimensional Stokes complex also exists, though for brevity, the interested reader is referred to [8] for more details.

It has been shown in previous works that the Stokes complex endows the generalized Stokes and Oseen problems with important underlying topological structure. In particular, the infinite-dimensional inf-sup condition may be derived from the complex [8]. As such, there is impetus for developing finite-dimensional approximations of the Stokes complex. Such discrete complexes are referred to as discrete Stokes complexes, and when these complexes are endowed with special commuting projection operators, they form the following commuting diagram with the Stokes complex:

[TABLE]

where $\Phi_{h}$, $\bm{\Psi}_{h}$, ${\bf V}_{h}$, and $\text{Q}_{h}$ are discrete scalar potential, vector potential, velocity, and pressure spaces and $\Pi_{\phi}:\Phi\rightarrow\Phi_{h}$, $\Pi_{\psi}:\bm{\Psi}\rightarrow\bm{\Psi}_{h}$, $\Pi_{v}:{\bf V}\rightarrow{\bf V}_{h}$, and $\Pi_{q}:\text{Q}\rightarrow\text{Q}_{h}$ are the aforementioned commuting projection operators. Remarkably, when ${\bf V}_{h}$ and $\text{Q}_{h}$ are selected as velocity and pressure approximation spaces in a mixed Galerkin formulation of the generalized Stokes or Oseen problem, the resulting approximation scheme is inf-sup stable and free of spurious oscillations and the returned discrete velocity solution will be pointwise divergence-free [3, 4]. Both of these properties are a direct consequence of the commuting diagram above, and for the sake of completeness, we prove the second property below.

Lemma 1**.**

Assume that the discrete velocity and pressure spaces ${\bf V}_{h}$ and $\text{Q}_{h}$ are associated with a discrete complex which commutes with the Stokes complex. Suppose ${\bf v}_{h}\in{\bf V}_{h}$ satisfies $b({\bf v}_{h},q_{h})=0$ for every $q_{h}\in\text{Q}_{h}$. Then $\nabla\cdot{\bf v}_{h}=0$ pointwise.

Proof.

Let $q_{h}=\nabla\cdot{\bf v}_{h}$. Then $\|\nabla\cdot{\bf v}_{h}\|^{2}_{L^{2}(\Omega)}=b({\bf v}_{h},q_{h})=0$ and the desired result follows. ∎

While we have demonstrated the benefit of using velocity and pressure spaces coming from a discrete Stokes complex, we have not yet described how to arrive at such spaces. In this paper, we turn to the use of so-called isogeometric compatible B-spline discretizations which are the focus of the next two subsections.

4.2 Univariate and Multivariate B-splines

The basic building blocks of isogeometric compatible B-spline discretizations, like any isogeometric analysis technology, are B-splines [35]. B-splines are piecewise polynomial functions, but unlike $C^{0}$-continuous finite elements, B-splines may exhibit high levels of continuity. Univariate B-splines are constructed by first specifying a polynomial degree $p$888The notation $p$ is used for both the pressure field as well as the polynomial degree. Thus, the reader should discern what term $p$ refers to in various portions of the paper by context., a number of basis functions $n$, and an open knot-vector $\Xi=\left\{\xi_{0},\xi_{1},\ldots,\xi_{n+p+1}\right\}$, a non-decreasing vector of knots $\xi_{i}$ such that the first and last knot are repeated $p+1$ times. We assume without loss of generality that the first and last knot are 0 and 1 respectively such that the domain of the knot vector is $(0,1)$. With a knot vector in hand, univariate B-spline basis functions are defined recursively through the Cox-deBoor formula:

[TABLE]

Figure 1 shows example sets of unvariate B-spline basis functions. We can alternatively define B-splines not from the knot vector itself, but instead a vector of unique knot values $\bm{\zeta}=\left\{\zeta_{1},\zeta_{2},\ldots,\zeta_{n_{k}}\right\}$ and a regularity vector $\bm{\alpha}=\left\{\alpha_{1},\alpha_{2},\ldots,\alpha_{n_{k}}\right\}$ such that the B-splines have $\alpha_{j}$ continuous derivatives across $\zeta_{j}$. By construction, $\alpha_{1}=\alpha_{n_{k}}=-1$. We will later employ the convention $\bm{\alpha}-1=\left\{-1,\alpha_{2}-1,\ldots,\alpha_{n_{k}-1}-1,-1\right\}$.

Given a set of knot-vectors and polynomial degrees, multivariate B-spline basis functions are obtained through a tensor-product of unvariate B-spline basis functions:

[TABLE]

where ${\bf i}=(i_{1},i_{2},\ldots,i_{d})$ and ${\bf p}=(p_{1},p_{2},\ldots,p_{d})$ We denote the corresponding space of multidimensional B-splines over the parametric domain $\hat{\Omega}=(0,1)^{d}$ as:

[TABLE]

where $\bm{\alpha}_{j}$ is the regularity vector associated with the $j^{\text{th}}$ direction where $j=1,\ldots,d$ and $\bm{\mathcal{M}}_{h}$ is the parametric mesh defined by the vectors of unit knot values in each parameteric direction. Note that the space is fully characterized by the polynomial degrees, regularity vectors, and parametric mesh as indicated by the notation. For ease of notation, however, we drop the dependence on the parameteric mesh and instead use $S^{p_{1},p_{2},\ldots,p_{d}}_{\bm{\alpha}_{1},\bm{\alpha}_{2},\ldots,\bm{\alpha}_{d}}=S^{p_{1},p_{2},\ldots,p_{d}}_{\bm{\alpha}_{1},\bm{\alpha}_{2},\ldots,\bm{\alpha}_{d}}\left(\mathcal{M}_{h}\right)$ in what follows.

4.3 Isogeometric Compatible B-splines

We are now in a position to define isogeometric compatible B-splines. Their definition is made possible through the observation that the derivative of univariate B-splines of degree $p$ are univariate B-splines of degree $p-1$. Since multivariate B-splines are tensor-products of univariate B-splines, the aforementioned property naturally generalizes to higher dimension, allowing us to build a discrete Stokes complex of B-spline spaces [8, 36]. We first define such a discrete Stokes complex in the parametric domain $\hat{\Omega}=(0,1)^{d}$ for both $d=2$ and $d=3$ before constructing a discrete Stokes complex in the physical domain of interest using a set of structure-preserving push-forward/pull-back operators.

In the two-dimensional setting, we define the following B-spline spaces over the unit square:

[TABLE]

where $\hat{\Psi}_{h}$ is the B-spline space of streamfunctions, $\hat{\bf V}_{h}$ is the B-spline space of flow velocities, and $\hat{\text{Q}}_{h}$ is the B-spline space of pressures. These discrete spaces are endowed with B-spline basis functions $\{\hat{N}^{\psi}_{i}\}_{i=1}^{n_{\psi}}$, $\{\hat{\bf N}^{v}_{i}\}_{i=1}^{n_{v}}$, and $\{\hat{N}^{p}_{i}\}_{i=1}^{n_{q}}$, respectively, where $n_{\psi}$ is the number of streamfunction basis functions, $n_{v}$ is the number of velocity basis functions, and $n_{q}$ is the number of pressure basis functions, all of which can be inferred from the chosen polynomial degrees and knot vectors. One can readily show that these spaces form the following discrete Stokes complex:

[TABLE]

and provided the functions in the B-spline pressure space are at least $C^{0}$-continuous, there exist a set of commuting projection operators that make the above discrete complex commute with the Stokes complex. Thus, we refer to the spaces $\hat{\Psi}_{h}$, $\hat{{\bf V}}_{h}$, and $\hat{\text{Q}}_{h}$ as compatible B-spline spaces. As mentioned previously, if we select $\hat{{\bf V}}_{h}$ and $\hat{\text{Q}}_{h}$ as velocity and pressure approximation spaces in a mixed Galerkin formulation of the generalized Stokes or Oseen problems, then the resulting scheme yields a pointwise divergence-free velocity field. The degrees of freedom associated with compatible B-splines are associated with the geometrical entries of the underlying control mesh. This is graphically illustrated in Figure 2 which shows that streamfunction degrees of freedom are associated with control points, velocity degrees of freedom are associated with (and aligned normal to) control edges, and pressure degrees of freedom are associated with control cells. Each degree of freedom corresponds to a particular basis function, and to visualize these basis functions, we have selected four degrees of freedom in Figure 2 and visualized the respective basis functions in Figure 3999Note that the pressure basis function we have highlighted does not have zero average over the parametric domain. In practice, we enforce this constraint using a Lagrange multiplier rather than to the individual pressure basis functions..

In the three-dimensional setting, we define the following B-spline spaces over the unit cube:

[TABLE]

where $\hat{\Phi}_{h}$ is the B-spline space of scalar potentials, $\hat{\bf\Psi}_{h}$ is the B-spline space of vector potentials, $\hat{\bf V}_{h}$ is the B-spline space of flow velocities, and $\hat{\text{Q}}_{h}$ is the B-spline space of pressures. These discrete spaces are endowed with the basis functions $\{\hat{N}^{\phi}_{i}\}_{i=1}^{n_{\phi}}$, $\{\hat{\bf N}^{\psi}_{i}\}_{i=1}^{n_{\psi}}$, $\{\hat{\bf N}^{v}_{i}\}_{i=1}^{n_{v}}$, and $\{\hat{N}^{p}_{i}\}_{i=1}^{n_{q}}$, respectively, where $n_{\phi}$ is the number of scalar potential basis functions, $n_{\psi}$ is the number of vector potential basis functions, $n_{v}$ is the number of velocity basis functions, and $n_{q}$ is the number of pressure basis functions, all of which can be inferred from the chosen polynomial degrees and knot vectors. Once again, one can show that the above spaces form the following discrete Stokes complex:

[TABLE]

and provided the functions in the B-spline pressure space are at least $C^{0}$-continuous, there exist a set of commuting projection operators that make the above discrete complex commute with the Stokes complex.

Heretofore, we have discussed how to construct compatible B-splines in the parametric domain. To define compatible B-splines in the physical domain $\Omega$, we need to first define a piece-wise smooth bijective mapping ${\bf F}:\hat{\Omega}\rightarrow\Omega$. This mapping can be defined using Non-Uniform Rational B-splines (NURBS), for instance, as is commonly done in the isogeometric analysis community [35]. With this mapping in hand, we define two-dimensional compatible B-spline spaces in the physical domain via the relations:

[TABLE]

and three-dimensional compatible B-spline spaces via the relations:

[TABLE]

where ${\bf J}=\partial_{\bm{\xi}}{\bf F}$ is the Jacobian of the parametric mapping. Corresponding basis functions in the physical domain are defined via push-forwards of the basis functions in the parametric domain, and we denote the discrete velocity basis functions as $\{{\bf N}^{v}_{i}\}_{i=1}^{n_{v}}$ and the basis functions for other quantities in analogous fashion. It is easily shown that the compatible B-spline spaces in the physical domain also comprise a discrete complex which commutes with the Stokes complex. The compatible B-splines in the physical domain are referred to as isogeometric compatible B-splines as they are built from B-splines, the basis building blocks of geometric modeling, and they are defined on the exact geometry of the problem of interest.

4.4 B-spline Refinement

One more concept needs to be introduced before proceeding forward, namely the concept of B-spline refinement. For a fixed set of polynomial degrees, B-spline refinement is carried out by a process referred to as knot insertion [37]. In the univariate setting, we start with a particular knot vector $\Xi$ and then insert a sequence of knots to arrive at a refined knot vector $\tilde{\Xi}$ such that $\Xi\subset\tilde{\Xi}$. The B-spline basis functions associated with the original knot vector, denoted as $\left\{\hat{N}_{i_{p}}(\xi)\right\}_{i=1}^{n}$, can be represented as linear combinations of the basis functions associated with the refined knot vector, denoted as $\left\{\tilde{N}_{i,p}(\xi)\right\}_{i=1}^{\tilde{n}}$, using a transformation matrix ${\bf T}$. This relationship is expressed mathematically as:

[TABLE]

for $i=1,\ldots,n$. Consequently, if a B-spline function takes the form:

[TABLE]

it can be alternately be represented as:

[TABLE]

where:

[TABLE]

for $j=1,\ldots,\tilde{n}$. Figure 4 depicts the action of knot insertion for univariate quadratic B-splines. B-spline refinement in the multivariate setting (including the compatible B-spline setting) is carried out in a tensor-product fashion, and the transformation matrix ${\bf T}$ takes the same form in both the parametric domain and the physical domain. There exist a variety of algorithms capable of performing knot insertion [37, 38, 39] which can be used to construct the transformation matrix ${\bf T}$, so we do not discuss this construction further in this paper.

5 A Structure-Preserving Geometric Multigrid Methodology

At last, we are ready to present our geometric multigrid methodology for isogeometric compatible discretizations of the generalized Stokes and Oseen problems. We begin this section by reviewing the basics of the geometric multigrid approach as well as the required ingredients in the setting of an isogeometric compatible discretization. Then, we introduce the Schwarz-style smoothers which our methodology leans upon. We then show that our methodology preserves the divergence-free constraint on the velocity field and that it effectively ellipticizes the underlying system of interest. We limit our discussion to the V-cycle algorithm, though our approach can also be applied within a W-cycle or Full Multigrid framework [40].

5.1 Nested B-spline Stokes Complexes, Intergrid Transfer Operators, and the V-Cycle Algorithm

Assume that we have a sequence of nested B-spline Stokes complexes that have been obtained through knot insertion. We denote the discrete velocity and pressure spaces associated with this sequence as $\left\{{\bf V}_{\ell}\right\}_{\ell=0}^{n_{\ell}}$ and $\left\{\text{Q}_{\ell}\right\}_{\ell=0}^{n_{\ell}}$ respectively where $n_{\ell}$ is the number of levels, and we note that:

[TABLE]

and:

[TABLE]

for each $\ell=0,\ldots,n_{\ell}$. Level $\ell=0$ corresponds to the coarsest mesh while level $\ell=n_{\ell}$ corresponds to the finest mesh. The action of knot insertion not only allows for B-spline refinement, but it also provides the intergrid transfer operators associated with a geometric multigrid method. Namely, we can build prolongation operators:

[TABLE]

for $\ell=0,\ldots,n_{\ell}-1$ using the construction provided in Subsection 4.4. We encode the action of these prolongation operators in the matrices ${\bf P}^{v}_{\ell}$ and ${\bf P}^{q}_{\ell}$ such that the following refinement operations hold:

[TABLE]

for $\ell=0,\ldots,n_{\ell}-1$ where $\left\{{\bf N}^{v}_{i,\ell}\right\}_{i=1}^{n_{v,\ell}}$ and $\left\{N^{q}_{i,\ell}\right\}_{i=1}^{n_{q,\ell}}$ denote the velocity and pressure B-spline basis functions associated with level $\ell$. Moreover, the degrees of freedom associated with pressure and velocity fields on the $\ell^{\textup{th}}$ level can be transferred to the the $(\ell+1)^{\textup{st}}$ level via the expressions:

[TABLE]

As is standard with a Galerkin formulation, restriction operators are constructed as the adjoint or transpose of the prolongation operators, namely $R^{v}_{\ell+1}=\left(P^{v}_{\ell}\right)^{*}$, $R^{q}_{\ell+1}=\left(P^{q}_{\ell}\right)^{*}$, ${\bf R}^{v}_{\ell+1}=\left({\bf P}^{v}_{\ell}\right)^{T}$, and ${\bf R}^{q}_{\ell+1}=\left({\bf P}^{q}_{\ell}\right)^{T}$ for $\ell=0,\ldots,n_{\ell}-1$. Finally, we define a prolongation matrix ${\bf P}_{\ell}$ for the full group variable such that:

[TABLE]

for $\ell=0,\ldots,n_{\ell}-1$. The corresponding restriction matrix for level $\ell$ is given by ${\bf R}_{\ell+1}=\left({\bf P}_{\ell}\right)^{T}$.

We need a few more ingredients before stating the multigrid V-cycle algorithm for our discretization scheme. First of all, we need to form the matrix system associated with the finest level, ${\bf K}{\bf U}={\bf F}$. We then form the system matrices associated with coarser levels via the relation ${\bf K}_{\ell}={\bf R}_{\ell+1}{\bf K}_{\ell+1}{\bf P}_{\ell}$ for $\ell=0,\ldots,n_{\ell}-1$ where ${\bf K}_{n_{\ell}}={\bf K}$. Second of all, we need to choose a smoother for each level $\ell$ which we encode in a smoothing matrix ${\bf S}_{\ell}$, and and we need to select a number of pre-smoothing steps $\nu_{1}$ and post-smoothing steps $\nu_{2}$. Third of all, we need to choose a suitable initial guess ${\bf U}$ for the solution on the finest level. Then, one V-cycle corresponds to a single call of the form $\text{MGV}(n_{\ell},{\bf U},{\bf F})$ to the recursive function defined below [40].

Note that the solution ${\bf U}$ is updated within the algorithm stated above. Hence, additional V-cycles simply correspond to additional calls of the form $\text{MGV}(n_{\ell},{\bf U},{\bf F})$.

5.2 Overlapping Schwarz Smoothers on Compatible Subdomains

At this juncture, we have not yet determined what smoother to employ. We turn to the use of overlapping Schwarz smoothers [21] with specially chosen overlapping subdomains which respect the underlying topological structure of the generalized Stokes and Oseen problems [24]. Namely, for each level $\ell$, we define a collection of subdomains $\{\Omega_{i,\ell}\}_{i}$ where each individual subdomain is defined as the support of a discrete streamfunction basis function in the two-dimensional setting:

[TABLE]

and a discrete vector potential basis function in the three-dimensional setting:

[TABLE]

It is easily seen that the subdomains form a cover of the physical domain, that is:

[TABLE]

For each subdomain, we define discrete velocity and pressure subspaces ${\bf V}_{i,\ell}\subset{\bf V}_{\ell}$ and $\text{Q}_{i,\ell}\subset\text{Q}_{\ell}$, respectively, as

[TABLE]

In the two-dimensional setting, we define a discrete streamfunction subspace for each subdomain as:

[TABLE]

and in the three-dimensional setting, we define a discrete vector potential subspace for each subdomain as:

[TABLE]

The degrees of freedom associated with all of the aforementioned subspaces are illustrated in the two-dimensional case in Figure 5 for two separate subdomains.

By construction, $\dim{\left(\Psi_{i,\ell}\right)}=1$ in the two-dimensional setting and $\dim{\left(\bm{\Psi}_{i,\ell}\right)}=1$ in the three-dimensional setting. Moreover, $\dim{\left({\bf V}_{i,\ell}\right)}=4$ and $\dim{\left(\text{Q}_{i,\ell}\right)}=3$ in both the two- and three-dimensional settings101010From Figure 5, it appears that $\dim{\left(\text{Q}_{i,\ell}\right)}=4$. However, the functions in $\text{Q}_{i,\ell}$ must satisfy a zero average constraint, so the dimension is one less than what is observed from the figure.. Thus, the subspaces associated with each subdomain form the following exact discrete Stokes complex in the two-dimensional setting:

[TABLE]

and the following exact discrete Stokes complex in the three-dimensional setting:

[TABLE]

Thus, as previously suggested, our choice of subdomains indeed respects the underlying topological structure of the generalized Stokes and Oseen problems.

With our subdomains defined, we can now describe our choice of smoothers, namely additive and multiplicative Schwarz smoothers using our prescribed subdomains. In this direction, let ${\bf E}^{v}_{i,\ell}$ and ${\bf E}^{q}_{i,\ell}$ denote the velocity and pressure subdomain restriction matrices for a given level $\ell$ and subdomain $i$ that take the full set of velocity and pressure degrees of freedom associated with level $\ell$ and map them to the set of pressure and velocity degrees of freedom associated with the subdomain $\Omega_{i,\ell}$. Additionally, let:

[TABLE]

denote the subdomain restriction matrix for a given level $\ell$ and subdomain $i$ for the full group variable. Then, the action of the additive Schwarz smoother is defined through:

[TABLE]

where $\eta\in(0,1)$ is a suitably chosen scaling factor [23], while the action of the multiplicative Schwarz smoother is defined through:

[TABLE]

The additive and multiplicate Schwarz smoothers are generalizations of the classical Jacobi and Gauss-Seidel smoothers, and indeed they can be implemented in an efficient, iterative manner. For both of these smoothers, a sequence of local matrix problems of the form:

[TABLE]

where ${\bf K}_{i,\ell}={\bf E}_{i,\ell}{\bf K}_{\ell}{\bf E}_{i,\ell}^{T}$ must be solved. It is easily seen that:

[TABLE]

where:

[TABLE]

Thus, with both the additive and multiplicative Schwarz smoothers, a discrete generalized Stokes or Oseen problem is solved for each subdomain. In the next subsection, we further clarify this interpretation in a variational setting. With the additive Schwarz smoother, the subdomain problems are solved independently, and their respective solutions are summed together and multiplied through by a scaling factor as indicated above. With the multiplicative Schwarz smoother, the subdomain problems are solved in a sequential fashion in analogy with the Gauss-Seidel smoother.

5.3 Preservation of the Divergence-free Constraint

Now that we have presented our choice of smoother, we demonstrate that our geometric multigrid methodology preserves the divergence-free constraint on the velocity field. Provided that the initial guess for the V-cycle algorithm satisfies the divergence-free constraint, it is sufficient to show that each smoothing step provides velocity updates that are divergence-free. One application of either the additive or multiplicative Schwarz smoother at level $\ell$ is akin to solving a collection of local subdomain problems of the form: Find $\delta{\bf u}_{i,\ell}\in{\bf V}_{i,\ell}$ and $\delta p_{i,\ell}\in\text{Q}_{i,\ell}$ such that:

[TABLE]

for all ${\bf v}_{h}\in{\bf V}_{i,\ell}$ and $q_{h}\in\text{Q}_{i,\ell}$ where ${\bf u}_{h}$ and $p_{h}$ are the approximate discrete velocity and pressure solutions. For the additive Schwarz smoother, the approximate discrete velocity and pressure solutions are updated following the solution of all of the local problems according to:

[TABLE]

while for the multiplicative Schwarz smoother, the approximate discrete velocity and pressure solutions are updated following the solution of each individual local problem according to:

[TABLE]

For each subdomain problem, if the approximate discrete velocity solution is divergence-free, it holds that:

[TABLE]

for all $q_{h}\in\text{Q}_{i,\ell}$. Since $\text{Q}_{i,\ell}=\nabla\cdot{\bf V}_{i,\ell}$, we can select $q_{h}=\nabla\cdot\delta{\bf u}_{i,\ell}$ to find:

[TABLE]

and thus the solution to the local problem is also divergence-free. Thus, if the initial guess for the V-cycle algorithm satisfies the divergence-free constraint, each subsequent application of the Schwarz smoother at any given level $\ell$ preserves the divergence-free constraint as well.

5.4 Efficacy of the Structure-Preserving Geometric Multigrid Methodology

We conclude here with a short discussion of the efficacy of our geometric multigrid methodology. We restrict our discussion to the three-dimensional setting without loss of generality. Recall that the spaces $\bm{\Psi}_{i,\ell}$, ${\bf V}_{i,\ell}$, and $\text{Q}_{i,\ell}$ form a discrete Stokes complex for a given level $\ell$ and subdomain $i$. Thus, we can express the velocity solution $\delta{\bf u}_{i,\ell}\in{\bf V}_{i,\ell}$ to (43) in terms of the curl of a vector potential $\delta\bm{\psi}_{i,\ell}\in\bm{\Psi}_{i,\ell}$ provided the velocity is divergence-free, and this vector potential can be obtained via the reduced subdomain problem: Find $\delta\bm{\psi}_{i,\ell}\in\bm{\Psi}_{i,\ell}$ such that:

[TABLE]

for all $\bm{\zeta}_{h}\in\bm{\Psi}_{i,\ell}$. This is precisely the subdomain problem associated with the global semi-elliptic generalized Maxwell problem with hyperresitivity [41]: Find $\bm{\psi}\in\bm{\Psi}_{h}$ such that:

[TABLE]

for all $\bm{\zeta}^{h}\in\bm{\Psi}_{h}$. It is known that a geometric multigrid methodology based on the use of Schwarz smoothers posed on structure-preserving subdomains is optimally convergent for Maxwell problems [27]. Consequently, we can expect that at least the discrete velocity solutions will converge in our approach.

6 Numerical Results

We now present a series of numerical tests illustrating the effectiveness of our proposed geometric multigrid methodology. Each of the tests correspond to problems with homogeneous Dirichlet boundary conditions applied along the entire domain boundary. In our discretization scheme, no-penetration boundary conditions are enforced strongly and no-slip boundary conditions are enforced weakly using a penalty constant of $C_{I}=4(p-1)$ where $p$ is the polynomial degree which is taken to be equal in each parameteric direction. It should be noted that $p$ refers to the polynomial degree of the discrete streamfunction space in the two-dimensional case and the discrete scalar potential space in the three-dimensional case. Hence, for $p=2$, the discrete pressure fields are piecewise bilinear/trilinear B-splines rather than piecewise biquadratic/triquadratic B-splines. Maximally smooth B-splines defined on uniform knot vectors are utilized throughout.

For all of the following tests, we define convergence as the number of V-cycles required to reduce the initial residual by a factor of $10^{6}$. We always initialize the V-cycle algorithm using a random initial guess which satisfies the divergence-free constraint on the velocity field. For each V-cycle, one pre-smoothing and two post-smoothing steps are employed using either the multiplicative or additive Schwarz smoother. For the additive Schwarz smoother, a scaling factor of $\eta=0.5$ is employed.

For all the problems presented here, a single element is used for the coarsest mesh and we investigate the convergence behavior for various levels of refinement. We report on the convergence behavior of our method for both the generalized Stokes and Oseen problems as well as a selection of different problem parameters, polynomial degrees, domain geometries, and number of spatial dimensions (2D and 3D). With respect to problem parameters, we consider the ratios between reaction and diffusion and advection and diffusion, which we express through a Damköhler number ($Da$) and a Reynolds number ($Re$). We define these numbers as:

[TABLE]

where $L$ is a characteristic length scale which is taken to be one throughout.

6.1 Two-dimensional generalized Stokes flow in a square domain

We first consider a two-dimensional generalized Stokes problem posed on the square domain $(0,1)^{2}$. In particular, we consider a forcing:

[TABLE]

corresponding to the manufactured solution [1]:

[TABLE]

[TABLE]

The velocity field associated with this exact solution is plotted in Figure 6.

We first present convergence results for our multigrid method using the multiplicative Schwarz smoother and polynomial degrees $p=2$ and $p=3$ in Table 1. It is clear that for a given polynomial order, the convergence behavior is robust with respect to both the number of levels of refinement and the problem parameters. We also observe that as polynomial order is increased, the convergence behavior deteriorates, albeit slightly. This is consistent with previously observed behavior for isogeometric analysis [42, 43, 44].

We next present convergence results for our multigrid method using the additive Schwarz smoother and polynomial degree $p=2$ in Table 2. As expected, overall convergence is slower with additive Schwarz than with multiplicative Schwarz, although in this case the method is still robust with regard to both the number of levels of refinement and the problem parameters.

6.2 Two-dimensional generalized Stokes flow in a quarter annulus

We next consider a two-dimensional generalized Stokes problem posed on a quarter annulus. The domain is described in Figure 7 where $r_{i}=0.075$ and $r_{o}=0.225$. We consider a manufactured solution achieved by mapping the solution presented in (47)-(48) to the quarter-annulus domain using a quadratic rational Bézier parametric mapping and appropriate push-forward operators. The velocity field associated with the exact solution are also plotted in Figure 7.

We present convergence results for our multigrid method using the multiplicative Schwarz smoother and polynomial degree $p=2$ in Table 3. Compared with the square domain, the number of V-cycles required for convergence is larger. However, the method is still robust with respect to both the number of levels of refinement and the problem parameters.

6.3 Three-dimensional generalized Stokes flow in a cube domain

We next consider a three-dimensional generalized Stokes problem posed on the unit cube $(0,1)^{3}$. In particular, we consider a forcing:

[TABLE]

corresponding to the manufactured solution [3]:

[TABLE]

[TABLE]

[TABLE]

Streamlines colored by velocity magnitude associated with the exact solution are plotted in Figure 8.

We present convergence results for our multigrid method using the multiplicative Schwarz smoother and polynomial degrees $p=2$ and $p=3$ in Table 4. Convergence appears to be much quicker in the three-dimensional setting. Notably, one V-cycle appears to be sufficient to reduce the residual by six orders of magnitude for a sufficient number of levels for both $p=2$ and $p=3$ and irrespective of the Damköhler number. We believe this may be due to the fact each velocity degree of freedom is updated twice as many times in each iteration of the smoother in the three-dimensional case as compared to the two-dimensional case.

6.4 Three-dimensional generalized Stokes flow in a hollow cylinder section

The final generalized Stokes problem considered in this paper is a three-dimensional problem posed on a hollow cylinder section. The domain for this problem is simply the quarter annulus from before extruded in the $z$-direction by a depth of $d=0.1$. We consider a manufactured solution achieved by mapping the solution presented in (50)-(52) to the hollow cylinder section domain using a quadratic rational Bézier parametric mapping and appropriate push-forward operators. Streamlines colored by velocity magnitude associated with the exact solution are plotted in Figure 9.

We present convergence results for our multigrid method using the multiplicative Schwarz smoother and polynomial degree $p=2$ in Table 5. Incredibly, one V-cycle again appears to be sufficient to reduce the residual by six orders of magnitude for a sufficient number of levels irrespective of the Damköhler number.

6.5 Two-dimensional generalized Oseen flow in a square domain

We now turn our attention to the generalized Oseen problem. We first consider a two-dimensional generalized Oseen problem posed on the square domain $(0,1)^{2}$. We manufacture a solution with a forcing:

[TABLE]

where ${\bf u}$ and $p$ are defined as in (47)-(48) such that the resulting solution is the same as the unit square generalized Stokes problem. Note that the advection velocity is taken to be the manufatured velocity field.

We present convergence results for our multigrid method using the multiplicative Schwarz smoother, polynomial degree $p=2$, and various Reynolds and Damköhler numbers in Table 6. When Reynolds number is low, the advection terms become negligible, and thus the method performs as it did on the the 2D generalized Stokes problem. As the Reynolds number is increased, the advection term becomes more significant. In this case, we have observed favorable convergence behavior as long as the Damköhler number is at least as large as the Reynolds number. When the system becomes advection-dominated, on the other hand, the multigrid method fails to converge. We expect that improved results may be obtained through the use of an alternative smoother which respects the directionality of the advection velocity.

6.6 Three-dimensional generalized Oseen flow in a cube domain

We conclude by considering a three-dimensional generalized Oseen problem posed on the unit cube $(0,1)^{3}$. We manufacture a solution with a forcing:

[TABLE]

where ${\bf u}$ and $p$ are defined as in (50)-(52) such that the resulting solution is the same as the unit cube generalized Stokes problem. As with the two-dimensional generalized Oseen problem, the advection velocity is take to be the manufactured velocity field.

We present convergence results for our multigrid method using the multiplicative Schwarz smoother, polynomial degree $p=2$, and various Reynolds and Damköhler numbers in Table 7. The same trends that were observed for the two-dimensional case are observed here as well. Namely, when the system is not advection dominated, we achieve excellent convergence behavior. Also, as was the case with the generalized Stokes flow, the three-dimensional case exhibits improved convergence as compared to the two-dimensional case.

7 Conclusions

In this paper, we presented a structure-preserving geometric multigrid methodology for isogeometric compatible discretizations of the generalized Stokes and Oseen problems which relies upon Schwarz-style smoothers in conjunction with specially chosen subdomains. We proved that our methodology yields a pointwise divergence-free velocity field independent of the number of pre-smoothing steps, post-smoothing steps, grid levels, or cycles in a V-cycle implementation, and we demonstrated the efficiency and robustness of our methodology by numerical example. Specifically, we found that our methodology exhibits convergence rates independent of the grid resolution and flow parameters for the generalized Stokes problem as well as the generalized Oseen problem provided it is not advection-dominated. We also discovered that, somewhat surprisingly, our methodology exhibits improved convergence rates in the three-dimensional setting as compared with the two-dimensional setting.

We envision several avenues for future work. First of all, we plan to conduct a full mathematical analysis of our methodology. We anticipate that this analysis will largely follow the same program of work as laid out in a recent geometric multigrid paper for divergence-conforming discontinuous Galerkin formulations of Stokes flow [23]. Second, we would like to extend the applicability of our methodology to advection-dominated Oseen problems. We anticipate the need for upwind-based line smoothers in such a setting [45]. Third, we plan to extend our methodology to multi-patch geometries and adaptive isogeometric compatible discretizations [46, 47]. Initial results in this area are quite encouraging. Finally, we plan to extend our methodology to multi-physics problems, including coupled flow transport, fluid-structure, and magnetohydrodynamics.

8 Acknowledgement

This material is based upon work supported by the Air Force Office of Scientific Research under Grant No. FA9550-14-1-0113.