# Discontinuous Galerkin schemes for the Stokes problem

Mayssa Mroueh, Erell Jamelot, Pascal Omnes

► To cite this version:

Mayssa Mroueh, Erell Jamelot, Pascal Omnes. Discontinuous Galerkin schemes for the Stokes problem. EPJ Web of Conferences, 2024, 302, pp.03006. 10.1051/epjconf/202430203006 . cea-04902442

HAL Id: cea-04902442

https://cea.hal.science/cea-04902442v1

Submitted on 20 Jan 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discontinuous Galerkin schemes for the Stokes problem

Mayssa Mroueh<fnref n="1" />,<fnref n="3" />,*, Erell Jamelot<fnref n="1" />,**, and Pascal Omnes<fnref n="2" />,<fnref n="3" />,***

**Abstract.** Our aim is to develop a robust and flexible code to simulate flows in nuclear core reactors. Discontinuous Galerkin schemes are thus proposed here to solve the Stokes problem, which stands as a fundamental element in fluid mechanics. A priori error estimates are provided when the solution is weakly regular. Details on how to solve the large-scale linear system are given. We moreover explain how to discretize a singular source term. Finally, we give numerical results of two test-cases for which the solution is weakly regular.

## 1 Motivation

In order to develop a robust and flexible code to simulate flows in nuclear core reactors, we develop numerical tools to solve the Stokes problem, which characterizes incompressible Newtonian flows. In numerical analysis, the study of the Stokes problem plays an essential role in the development of an appropriate numerical scheme to simulate the Navier-Stokes equations. Moreover, it allows to model numerous phenomena observed in nature, ranging from blood flows to oceanic movements. The numerical resolution of this problem thus holds great significance, demanding efficient and accurate methods.

In this article, we study the resolution of the Stokes problem with some Discontinuous Galerkin schemes, based on the Symmetric Interior Penalty method, which is well suited to a discretization using quadratic polynomials.

In Section 2, we present the Stokes problem and its mathematical framework. In Section 3, we outline the variational formulation used to discretize this problem, giving a priori error estimates when the solution is weakly regular. In Section 4, we detail how we solve the large-scale linear system. In Section 5, we explain how we discretize a singular source term. In Section 6, we present the numerical results of our simulations, comparing the performance of the different approaches.

<footnote n="1">¹Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">²Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, 91191, Gif-sur-Yvette, France.</footnote>

<footnote n="3">³Université Sorbonne Paris Nord, LAGA, CNRS UMR 7539, Institut Galilée, 93430 Villetaneuse, France</footnote>

<footnote>*e-mail: mayssa.mroueh@cea.fr</footnote>

<footnote>**e-mail: erell.jamelot@cea.fr</footnote>

<footnote>***e-mail: pascal.omnes@cea.fr</footnote>
## 2 The Stokes Problem

Let $N_D \in \{2, 3\}$ be the dimension and $\mathcal{I}_D = \{x, y\}$ for $N_D = 2$, $\mathcal{I}_D = \{x, y, z\}$ for $N_D = 3$. Let $(O, x, y)$ for $N_D = 2$ and $(O, x, y, z)$ for $N_D = 3$ be the Cartesian coordinates system, of orthonormal basis $(e_d)_{d \in \mathcal{I}_D}$. Let $\Omega \subset \mathbb{R}^{N_D}$ be a connected and bounded domain, with polytopal<fnref n="1" /> Lipschitz boundary $\partial\Omega$. The vector field $\mathbf{u}$ represents the velocity of the fluid, and the scalar field $p$ represents the pressure divided by the fluid density, which is assumed to be constant. The constant parameter $\nu > 0$ represents the kinematic viscosity. The vector field $\mathbf{f}$ represents the external force field on the fluid, divided by the fluid density. Consider the Stokes problem:

$$ \text{Find } (\mathbf{u}, p) \text{ s.t. } -\nu\Delta\mathbf{u} + \nabla p = \mathbf{f} \quad (i), \quad \nabla \cdot \mathbf{u} = 0 \quad (ii) \text{ in } \Omega, \quad (1) $$

with Dirichlet boundary conditions for the velocity: $\mathbf{u} = 0$ on $\partial\Omega$; and a normalization condition for the pressure $p$: $\int_{\Omega} p = 0$. Equation (1)-(i) (resp. (1)-(ii)) corresponds to the momentum equation (resp. the conservation of mass equation).

Let $L_0^2(\Omega) = \{q \in L^2(\Omega) \mid \int_{\Omega} q = 0\}, H^1(\Omega) = \{v \in L^2(\Omega) \mid \text{grad } v \in (L^2(\Omega))^{N_D}\}, H_0^1(\Omega) = \{v \in H^1(\Omega) \mid v|_{\partial\Omega} = 0\}$. Recall the definition of fractional order Sobolev spaces cf. [6]. For all $0 < s < 1$, we set:

$$ H^s(\Omega) := \left\{ v \in L^2(\Omega) \mid \int_{\Omega} \int_{\Omega} \frac{|v(\mathbf{x}) - v(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^{d+2s}} d\mathbf{x} d\mathbf{y} < \infty \right\}. \quad (2) $$

Note that $H_0^1(\Omega) = H^1(\Omega)$, $\forall t \in [0, \frac{1}{2}]$. Let $X \in \{H_0^1(\Omega), L_0^2(\Omega), H^s(\Omega)\}$, $\mathbf{X} = X^{N_D}$ and $\mathbb{L}^2(\Omega) := (L^2(\Omega))^{N_D \times N_D}$. For $\mathbf{v} \in \mathbf{H}^1(\Omega)$, the tensor $\text{Grad } \mathbf{v} \in \mathbb{L}^2(\Omega)$ is s.t. for all $d, d' \in \mathcal{I}_D$, $(\text{Grad } \mathbf{v})_{d,d'} = \partial_{d'}v_d$. The variational formulation of Problem (1) reads: Find $(\mathbf{u}, p) \in \mathbf{H}_0^1(\Omega) \times L_0^2(\Omega)$ s.t. for all $(\mathbf{v}, q) \in \mathbf{H}_0^1(\Omega) \times L_0^2(\Omega)$

$$ v(\text{Grad } \mathbf{u}, \text{Grad } \mathbf{v})_{\mathbb{L}^2(\Omega)} - (\text{div } \mathbf{v}, p)_{\mathbb{L}^2(\Omega)} + (\text{div } \mathbf{u}, q)_{\mathbb{L}^2(\Omega)} = f(\mathbf{v}), \quad (3) $$

where $f: \mathbf{H}_0^1(\Omega) \to \mathbb{R}$ is s.t. $f(\mathbf{v}) = \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}_0^1(\Omega)}$.

Using [3, Corollary I-2.4], one proves that Problem (3) is well-posed. When $\Omega$ is polytopal convex and $\mathbf{f} \in \mathbb{L}^2(\Omega)$, the solution of (1) is s.t. $(\mathbf{u}, p) \in (\mathbf{H}_0^1(\Omega) \cap \mathbf{H}^2(\Omega)) \times (L_0^2(\Omega) \cap \mathbf{H}^1(\Omega))$. However, when $\Omega$ is not convex and $\mathbf{f} \in \mathbb{L}^2(\Omega)$ or when $\Omega$ is convex and $\mathbf{f} \in (\mathbf{H}^1(\Omega))'$, with $t \in [0, \frac{1}{2}]$, the fields $\mathbf{u}$ and $p$ are less regular and convergence of numerical methods are less studied, and more difficult to obtain.

## 3 Discretization

Consider $(\mathcal{T}_h)_h$ a polytopal triangulation sequence of $\Omega$, where $h$ denotes the mesh size. For a triangulation $\mathcal{T}_h$, we use the following index sets:

* Let $\mathcal{I}_K$ be the set of polytope indices, s.t. $\mathcal{T}_h := \bigcup_{\ell \in \mathcal{I}_K} K_\ell$ and $N_K = \text{card}(\mathcal{I}_K)$.

* Let $\mathcal{I}_F$ be the set of facet<fnref n="2" /> indices, s.t. $\mathcal{F}_h := \bigcup_{f \in \mathcal{I}_F} F_f$ and $N_F = \text{card}(\mathcal{I}_F)$.

* Define $\mathcal{I}_F = \mathcal{I}_F^i \cup \mathcal{I}_F^b$, where $\forall f \in \mathcal{I}_F^i, F_f \in \Omega$ and $\forall f \in \mathcal{I}_F^b, F_f \in \partial\Omega$.

<footnote n="1">¹i.e. polygonal ($N_D = 2$) or polyhedral ($N_D = 3$)</footnote>

<footnote n="2">²i.e. edge ($N_D = 2$) or face ($N_D = 3$)</footnote>
Let $N_{\partial}$ be the maximum number of facets of a polytopal element of $\mathcal{T}_h$ ($N_{\partial} = d + 1$ for a simplex), $\sigma$ be the shape-regularity parameter of $\mathcal{T}_h$. Let $\mathcal{I}_F = \mathcal{I}_F^i \cup \mathcal{I}_F^b$, where $\forall f \in \mathcal{I}_F^i$, $F_f \subset \Omega$; $\forall f \in \mathcal{I}_F^b$, $F_f \in \partial\Omega$ and $N_{\partial}^{i,b} = \text{card}(\mathcal{I}_F^{i,b})$. For all $f \in \mathcal{I}_F$, we let $\mathbf{n}_f$ be its unit normal (outward oriented if $F_f \in \partial\Omega$).

Let $\mathcal{P}_h H^1 = \{v \in L^2(\Omega); \forall \ell \in \mathcal{I}_K, v_{|K_\ell} \in H^1(K_\ell)\}$ and $\mathcal{P}_h H^1 = [\mathcal{P}_h H^1]^{N_D}$ and:

$$ (v, w)_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{grad} v, \mathbf{grad} w)_{L^2(K_\ell)} \quad \|v\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \|\mathbf{grad} v\|_{L^2(K_\ell)}^2 $$

$$ (\mathbf{v}, \mathbf{w})_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{Grad} \mathbf{v}, \mathbf{Grad} \mathbf{w})_{L^2(K_\ell)} \quad \|\mathbf{v}\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \|\mathbf{Grad} \mathbf{v}\|_{L^2(K_\ell)}^2 $$

Let $f \in \mathcal{I}_F^i$ s.t. $F_f = \partial K_L \cap \partial K_R$ and $\mathbf{n}_f$ is oriented outward $K_L$. Let $v \in H^s(\Omega)$, with $s > \frac{1}{2}$, the jump (resp. average) of $v$ across $F_f$ is defined as: $[v] := v_{|K_L} - v_{|K_R}$ (resp. $\{v\} := \frac{1}{2}(v_{|K_L} + v_{|K_R})$). For $f \in \mathcal{I}_F^b$, we set: $[v] := v_{|F_f}$ and $\{v\} := v_{|F_f}$.

For all $k \in \mathbb{N}$ we denote $P_{\text{disc}}^k(\mathcal{T}_h) = \{v \in L^2(\Omega); \forall \ell \in \mathcal{I}_K, v_{|K_\ell} \in P^k(K_\ell)\}$, $\mathbf{P}_{\text{disc}}^k(\mathcal{T}_h) = [P_{\text{disc}}^k(\mathcal{T}_h)]^{N_D}$ and $N_k := \dim(P^k(K_\ell)) = \frac{(k+N)!}{k!N!}$. Let $\mathcal{I}_k = \{1, \dots, N_k\}$. Call $(\phi_{\ell,i})_{i \in \mathcal{I}_k}$ a polynomial basis of $P^k(K_\ell)^3$ We discretize Problem (1) using discontinuous Galerkin schemes [2, Chapter 6], [7]. Let's denote by $k_u \in \mathbb{N}^*$ (resp. $k_p \in \mathbb{N}$) the polynomial approximation order of the velocity (resp. pressure). We set $\mathbf{X}_h = \mathbf{P}_{\text{disc}}^{k_u}(\mathcal{T}_h)$ and $L_h := P_{\text{disc}}^{k_p}(\mathcal{T}_h) \cap L_0^2(\Omega)$. The bilinear form $a_h: \mathbf{X}_h \times \mathbf{X}_h \to \mathbb{R}$ discretizes the first term of Eq. (3):

$$ \forall (\mathbf{v}_h, \mathbf{w}_h) \in \mathbf{X}_h \times \mathbf{X}_h, \quad a_h(\mathbf{v}_h, \mathbf{w}_h) = \sum_{d \in \mathcal{I}_D} a_h^{\text{sip}}((\mathbf{v}_h)_d, (\mathbf{w}_h)_d), \quad (4) $$

where for all $k \in \mathbb{N}^*$, $a_h^{\text{sip}}: P_{\text{disc}}^k(\mathcal{T}_h) \times P_{\text{disc}}^k(\mathcal{T}_h) \to \mathbb{R}$ is s.t.

$$ a_h^{\text{sip}}(v_h, w_h) = (v_h, w_h)_h + \sum_{f \in \mathcal{I}_F} (\{\mathbf{grad} v_h\} \cdot \mathbf{n}_f, [w_h])_{L^2(F_f)} \\ + \sum_{f \in \mathcal{I}_F} (\{\mathbf{grad} w_h\} \cdot \mathbf{n}_f, [v_h])_{L^2(F_f)}, \quad (5) \\ + \eta s_h(v_h, w_h) $$

$$ \text{with } s_h(v_h, w_h) := \sum_{f \in \mathcal{I}_F^i} (h_F)^{-1} ([v_h], [w_h])_{L^2(F_f)}, \quad \forall (v_h, w_h) \in P_{\text{disc}}^k(\mathcal{T}_h) \times P_{\text{disc}}^k(\mathcal{T}_h). $$

One can prove that when $\eta > \frac{(k+1)(k+N_D)}{N_D} \sigma^{N_D} N_{\partial}$, the bilinear form $a_h^{\text{sip}}(\cdot, \cdot)$ is coercive

with respect to the following norm: $\|v_h\|_{\text{sip}} := (\|v_h\|_h^2 + \eta^{-1} s_h(v_h, v_h))^{\frac{1}{2}}$.

The first and second (resp. third, fourth) term(s) in Equation (5) ensure consistency (resp. ensures symmetry, stability) of the bilinear form $a_h^{\text{sip}}(\cdot, \cdot)$.

The bilinear form $b_h: \mathbf{X}_h \times L_h \to \mathbb{R}$ discretizes the second term of Eq. (3):

$$ b_h(\mathbf{v}_h, q_h) = - \sum_{\ell \in \mathcal{I}_K} (\text{div} \, \mathbf{v}_h, q_h)_{L^2(K_\ell)} - \sum_{f \in \mathcal{I}_F} ([\mathbf{v}_h] \cdot \mathbf{n}_f, \{q_h\})_{L^2(F_f)}. \quad (6) $$

The linear form $f_h: \mathbf{X}_h \to \mathbb{R}$ discretizes the right-hand side term of Eq. (3):

$$ f_h(\mathbf{v}_h) = \langle \mathbf{f}, \mathbf{v}_h \rangle_{H^1(\Omega)} = \sum_{\ell \in \mathcal{I}_K} \langle \mathbf{f}, \mathbf{v}_h \rangle_{H^1(K_\ell)}. \quad (7) $$

<footnote n="3">³For $N_D = 2$, $k = 1$, letting $\mathbf{x}_\ell$ be the barycenter of $K_\ell$ and $h_\ell$ be its diameter, we consider $(\phi_{\ell,i})_{i \in \mathcal{I}_k} = \{1, \frac{x_\ell - x_\ell}{h_\ell}, \frac{y_\ell - y_\ell}{h_\ell}\}$, with the convention notation $\mathbf{x}_\ell := (x_\ell, y_\ell)^T$ and $(x, y) = (x_1, x_2)$.</footnote>
The discretization of Problem (3) reads:

Find $(\mathbf{u}_h, p_h) \in \mathbf{X}_h \times L_h$ s.t. for all $(\mathbf{v}_h, q_h) \in \mathbf{X}_h \times L_h$

$$v a_h(\mathbf{u}_h, \mathbf{v}_h) + b_h(\mathbf{v}_h, p_h) - b_h(\mathbf{u}_h, q_h) + \beta s_h(p_h, q_h) = f_h(\mathbf{v}_h), \quad (8)$$

with $s_h(v_h, w_h) := \sum_{f \in \mathcal{F}_F^i} h_F([v_h], [w_h])_{L^2(F_h)}$, $\forall (v_h, w_h) \in P_{disc}^k(\mathcal{T}_h) \times P_{disc}^k(\mathcal{T}_h)$.

Let us set $\mathcal{X} := \mathbf{H}_0^1(\Omega) \times \mathbf{L}_0^2(\Omega)$ which a Hilbert space endowed with the norm:

$$\text{For } (\mathbf{u}, p) \in \mathcal{X}, \quad \|(\mathbf{u}, p)\|_{\mathcal{X},v} = \left( \| \text{Grad } \mathbf{u} \|_{L^2(\Omega)}^2 + v^{-2} \|p\|_{L^2(\Omega)}^2 \right)^{\frac{1}{2}}. \quad (9)$$

**Proposition 1** Problem (8) is well-posed for $\beta = \begin{cases} 0 & \text{when } k_{\mathbf{u}} > k_p \\ 1 & \text{when } k_{\mathbf{u}} \in \{k_p, k_p - 1\} \end{cases}$. When $\beta = 0$ and $k_{\mathbf{u}} > k_p$, the bilinear form $b_h(\cdot, \cdot)$ satisfies a discrete inf-sup condition with a constant which does not depend on $h$.

The following error estimates can be proven:

**Theorem 3.1 (Convergence)** Let $(\mathbf{u}, p)$ be the unique solution to Problem (1) and $(\mathbf{u}_h, p_h) \in \mathbf{X}_h \times L_h$ be the unique solution to Problem (8).

Let $\Omega$ be a domain with a Lipschitz boundary. There exist two constants, $C_{\mathbf{u}}(\mathbf{f})$, $C_{\mathbf{p}}(\mathbf{f}) > 0$ independent of $h$ such that:

1. If $\Omega$ is convex and there exists $s \in ]1/2, 1]$ such that $\mathbf{f} \in \mathbf{H}^{s-1}(\Omega)$. Then, by interpolation, $(\mathbf{u}, p) \in \mathbf{H}_0^1(\Omega) \cap \mathbf{H}^{1+s}(\Omega) \times \mathbf{H}^s(\Omega) \cap \mathbf{L}_0^2(\Omega)$, and we have:

$$\|\mathbf{u} - \mathbf{u}_h\|_{L^2(\Omega)} \le v^{-1} C_{\mathbf{u}}(\mathbf{f}) h^{1+s} \quad \|p - p_h\|_{L^2(\Omega)} \le C_p(\mathbf{f}) h^s.$$

2. If $\Omega$ is non-convex, and $\mathbf{f} \in \mathbf{L}^2(\Omega)$, then there exists $s \in ]1/2, 1]$, such that $(\mathbf{u}, p) \in \mathbf{H}_0^1(\Omega) \cap \mathbf{H}^{1+s}(\Omega) \times \mathbf{H}^s(\Omega) \cap \mathbf{L}_0^2(\Omega)$ and we have:

$$\|\mathbf{u} - \mathbf{u}_h\|_{L^2(\Omega)} \le v^{-1} C_{\mathbf{u}}(\mathbf{f}) h^{2s} \quad \|p - p_h\|_{L^2(\Omega)} \le C_p(\mathbf{f}) h^s.$$

Denote by $N_{\mathbf{u}} = N_{k_{\mathbf{u}}} \times N_T$ (resp. $N_p = N_{k_p} \times N_T$) the number of degrees of freedom of $P_{disc}^{k_{\mathbf{u}}}(\mathcal{T}_h)$ (resp. $P_{disc}^{k_p}(\mathcal{T}_h)$). Let's call $A \in \mathbb{R}^{N_{\mathbf{u}}} \times \mathbb{R}^{N_{\mathbf{u}}}$ the matrix s.t.:

$$\forall (i, \ell), (j, \ell') \in \mathcal{I}_{k_{\mathbf{u}}} \times \mathcal{I}_{K}, \quad A_{(i, \ell), (j, \ell')} := a_h^{\text{isp}}(\phi_{i, \ell}, \phi_{j, \ell'}). \quad (10)$$

Let $A \in \mathbb{R}^{N_D \times N_{\mathbf{u}}} \times \mathbb{R}^{N_D \times N_{\mathbf{u}}}$ be the block diagonal matrix s.t.: $A = (\delta_{d,d'}A)_{d,d' \in \mathcal{I}_D}$. For $d \in \mathcal{I}_D$, we call $B_d \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_{\mathbf{u}}}$ the matrix s.t. for all $(i, \ell) \in \mathcal{I}_{k_p} \times \mathcal{I}_K$, $(j, \ell') \in \mathcal{I}_{k_{\mathbf{u}}} \times \mathcal{I}_K$, $(B_d)_{(i,\ell),(j,\ell')} := b_h(\phi_{\ell',j}\mathbf{e}_d, \phi_{i,\ell})$ and $B = (B_d)_{d \in \mathcal{I}_D} \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_D \times N_{\mathbf{u}}}$. Let $M_p \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$ and $S_p \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$ the matrices s.t. for all $(i, \ell), (j, \ell') \in \mathcal{I}_{k_p} \times \mathcal{I}_K$, $(M_p)_{(i,\ell),(j,\ell')} = \delta_{\ell,\ell'}(\phi_{i,\ell}, \phi_{\ell,j})_{L^2(\mathcal{K}_\ell)}$ and $(S_p)_{(i,\ell),(j,\ell')} = s_h(\phi_{i,\ell}, \phi_{\ell,j})$. Notice that $A, B_d$ and $S_p$ are sparse: for $\ell, \ell' \in \mathcal{I}_K$ s.t. $\partial K_\ell \cap \partial K_{\ell'} = \emptyset$, $A_{(i,\ell),(j,\ell')} = 0$, $(B_d)_{(i,\ell),(j,\ell')} = 0$ and $(S_p)_{(i,\ell),(j,\ell')} = 0$. The matrix $A$ is symmetric, positive definite. When $k_{\mathbf{u}} > k_p$, matrix $B$ is of rank $N_p - 1$. Notice that matrix $S_p$ is also of rank $N_p - 1$. When $k_{\mathbf{u}} \le k_p$, matrix $B$ is of rank $< N_p - 1$. Let $F_{\mathbf{u}} = (F_{\mathbf{u},d})_{d \in \mathcal{I}_D} \in \mathbb{R}^{N_D \times N_{\mathbf{u}}}$ the vector s.t. for all $(i, \ell) \in \mathcal{I}_{k_{\mathbf{u}}} \times \mathcal{I}_K$, $(F_{\mathbf{u},d})_{i,\ell} = f_h(\phi_{i,\ell}\mathbf{e}_d)$. It stands: $\mathbf{u}_h = \sum_{\ell \in \mathcal{I}_K} \sum_{i \in \mathcal{I}_{k_{\mathbf{u}}}} \mathbf{u}_{i,\ell} \phi_{i,\ell}$ and $p_h = \sum_{\ell \in \mathcal{I}_K} \sum_{i \in \mathcal{I}_{k_p}} p_{i,\ell} \phi_{i,\ell}$. Let's call $U_h = (U_{h,d})_{d \in \mathcal{I}_D} \in \mathbb{R}^{N_D \times N_{\mathbf{u}}}$ (resp. $P_h \in \mathbb{R}^{N_p}$) the vector of the discrete velocity (resp. pressure) degrees of freedom: for all $(i, \ell) \in \mathcal{I}_{k_{\mathbf{u}}} \times \mathcal{I}_K$
(resp. $(i, \ell) \in \mathcal{I}_{k_p} \times \mathcal{I}_K$), $(U_{h,d})_{i,\ell} = \mathbf{u}_{i,\ell} \cdot \mathbf{e}_d$ (resp. $(P_h)_{i,\ell} = p_{i,\ell}$). The linear system related to Problem (8) reads: Find $(U_h, P_h) \in \mathbb{R}^{N_D \times N_u} \times \mathbb{R}^{N_p}$ s.t.:

$$ \left\{ \begin{array}{l} \nu A U_h + B^T P_h = F_u \\ -B U_h + \beta S_p P_h = F_p \end{array} \right. \quad \text{and} \quad \sum_{\ell \in \mathcal{I}_K} \sum_{i \in \mathcal{I}_{k_p}} p_{i,\ell} \int_{K_\ell} \phi_{i,\ell} = 0, \qquad (11) $$

where $F_p \in \mathbb{R}^{N_p}$ vanishes when considering homogeneous Dirichlet boundary conditions. To compute $p_h \in L_h$, we can either impose the last equation in (11) at each iteration of our solver or cancel a row from matrix $B$ and the same row and corresponding column from matrices $M_p$ and $S_p$. We make the abuse of keeping the same notation. Let's call $K = BA^{-1}B^T + \beta v S_p \in \mathbb{R}^{N_p-1} \times \mathbb{R}^{N_p-1}$, which is a symmetric positive definite matrix. To solve the coupled velocity-pressure problem (11), one relies usually on the three + one steps below (the fourth step being straightforward):

$$ \begin{array}{ll} \text{Prediction:} & \text{Solve in } U_h^* \text{ s.t. } \nu A U_h^* = F_u. \\ \text{Pressure solver:} & \text{Solve in } P_h \text{ s.t. } KP_h = \nu(F_p + B U_h^*). \\ \text{Correction:} & \text{Solve in } \delta U_h \text{ s.t. } \nu A \delta U_h = -B^T P_h. \\ \text{Update:} & U_h = \delta U_h + U_h^*. \end{array} \qquad (12) $$

One can prove that matrix $K$ is equivalent to $M = M_p + \beta S_p$, so that $M$ can be used as a preconditioner to solve $KP_h = \nu(F_p + B U_h^*)$. Using some iterative solver for this system, each matrix-vector product with $K$ requires the solution of linear systems such as $AX = b$ at each iteration, where $A$ is defined by Eq. (10). Let us give details on our resolution algorithm.

## 4 Domain Decomposition and Schur Complement

Let $N \in \mathbb{N}^*$, $I_N = \{1, \dots, N\}$ and $\tilde{I}_N = \{1, \dots, N-1\}$. We split $\mathcal{T}_h$ into $N$ disjoint subsets $(\mathcal{T}_{h,n})_{n \in \tilde{I}_N}$. Let's denote by $(\mathcal{F}_{h,n})_{n \in \tilde{I}_N}$ the associated sets of facets. The splitting is s.t. for $n, m \in \tilde{I}_N$, for $n \neq m$, $\mathcal{F}_{h,n} \cap \mathcal{F}_{h,m} = \emptyset$. For $n \in \tilde{I}_N$, we consider matrices $A_{n,n} = (A_{\ell,\ell'})_{\ell \in \mathcal{T}_{h,n}, \ell' \in \mathcal{T}_{h,n}}$, $A_{n,N} = (A_{\ell,\ell'})_{\ell \in \mathcal{T}_{h,n}, \ell' \in \mathcal{T}_{h,N}}$ and for $n, m \in \tilde{I}_N$, $A_{n,m} = 0$. The matrix $A$ can be rewritten by blocks as follows: $A = ((A_{n,m})_{n,m \in \tilde{I}_N \times \tilde{I}_N})$. For all $n \in \tilde{I}_N$, $A_{n,n}$ is a symmetric positive definite matrix, and for all $n \in \tilde{I}_N$, $A_{n,N}$ is the coupling matrix between $\mathcal{F}_{h,n}$ and $\mathcal{F}_{h,N}$ (we have: $A_{N,n} = (A_{n,N})^T$). Similarly, we set $X = (X_n)_{n \in \tilde{I}_N}$ and $b = (b_n)_{n \in \tilde{I}_N}$. The linear system $AX = b$ is solved as follows:

$$ \text{For all } n \in \tilde{I}_N, \text{ compute} \quad \left\{ \begin{array}{l} \tilde{b}_n = A_{n,N}^T A_{n,n}^{-1} b_n \\ \tilde{A}_{n,n} = A_{n,N}^T A_{n,n}^{-1} A_{n,N} \end{array} \right. . \\ \text{Compute} \quad \left\{ \begin{array}{l} \tilde{b}_N = b_N - \sum_{n \in \tilde{I}_N} \tilde{b}_n \\ \tilde{A}_{N,N} = A_{NN} - \sum_{n \in \tilde{I}_N} \tilde{A}_{n,n} \end{array} \right. \qquad (13) $$

Solve for $X_N$ s.t. $\tilde{A}_{N,N} X_N = \tilde{b}_N$. For all $n \in \tilde{I}_N$, solve for $X_n$ s.t. $A_{n,n} X_n = b_n - A_{n,N} X_N$.

For $n \in \tilde{I}_N$, we store the Cholesky decomposition of matrix $A_{n,n}$ to compute vector $\tilde{b}_i$ and matrix $\tilde{A}_{n,n}$. These computations can be done in parallel. After, we build $\tilde{A}_{N,N}$ and store its LU decomposition. The matrices $(\tilde{A}_{n,n})_{n \in \tilde{I}_N}$ can be calculated once and for all in parallel. For
the resolution, we can also calculate the vectors $(\tilde{b}_n)_{n \in \tilde{I}_N}$ and $(X_n)_{n \in \tilde{I}_N}$ in parallel. This work is in progress.

The code is being developed in MATLAB, based on a prototype created by Erell Jamelot, which is available on GitHub [5]. This code will be published soon on GitHub.

In order to validate our theoretical analysis and our numerical developments, we study low-regular test-cases, when $\Omega$ is non-convex and $\mathbf{f} = 0$; and when $\Omega$ is convex and $\mathbf{f} \notin \mathbf{L}^2(\Omega)$. In the next section, we give details on this last implementation.

# 5 Numerical integration of a low-regular source term

Consider here the 2D case where $\Omega = (0, 1)^2$ and the meshes $(\mathcal{T}_h)_h$ are made of triangles. Let $S_0 = (x_0, y_0)$, where $x_0 = y_0 = 0.5$ and $\mathbf{x}_0 = \overrightarrow{OS_0}$.

Call $(r, \theta)$ the polar coordinates centered in $S_0$: $r = |\mathbf{x} - \mathbf{x}_0|$, $\theta = \arctan\left(\frac{y-y_0}{x-x_0}\right) + \lambda\pi$, with $\lambda = 1$ if $x \le x_0$, and elsewhere.

Let's set $\mathcal{I}_{K,x_0} := \{\ell \in \mathcal{I}_K : \mathbf{x}_0 \in K_\ell\}$ and define $\Omega_{x_0} = \{K_\ell, \ell \in \mathcal{I}_{K,x_0}\}$.

The polar basis $(\mathbf{e}_r, \mathbf{e}_\theta)$ is s.t.: $\mathbf{e}_r = \cos \theta \mathbf{e}_x + \sin \theta \mathbf{e}_y$, $\mathbf{e}_\theta = -\sin \theta \mathbf{e}_x + \cos \theta \mathbf{e}_y$.

Let $\alpha, \beta \in \mathbb{R}^*$. We will compute the following prescribed solution to Problem (1):

$$ (\mathbf{u}, p) = (r^\alpha \mathbf{e}_\theta, r^\beta - p_0), \quad \text{where } p_0 = \frac{1}{|\Omega|} \int_{\Omega} r^\beta d\mathbf{x} \quad (14) $$

When $\alpha = 1, \beta = 2$, it is close to the Rankine vortex model [11, Annexe B]. The source term is s.t. $\mathbf{f} \in (\mathbf{H}^1(\Omega))'$, with $t = 1-s$, $s = \min\{\alpha, \beta + 1\}$ and:

$$ \mathbf{f} = v(1 - \alpha^2) r^{\alpha-2} \mathbf{e}_{\theta} + \beta r^{\beta-1} \mathbf{e}_{r}. \quad (15) $$

Notice that $\mathbf{f}_{|\Omega_{x_0}} \in \mathbf{L}^1(\Omega)$ and $\mathbf{f}_{|\Omega \setminus \Omega_{x_0}} \in \mathbf{L}^2(\Omega)$. Let $\mathbf{v}_h \in \mathbf{X}_h$. It stands:

$$ \sum_{\ell \in \mathcal{I}_K} \langle \mathbf{f}, \mathbf{v}_h \rangle_{\mathbf{H}^1(K_\ell)} = \sum_{\ell \in \mathcal{I}_K \setminus \mathcal{I}_{K,x_0}} \int_{K_\ell} \mathbf{f} \cdot \mathbf{v}_h d\mathbf{x} + \sum_{\ell \in \mathcal{I}_{K,x_0}} \langle \mathbf{f}, \mathbf{v}_h \rangle_{\mathbf{H}^1(K_\ell)}. \quad (16) $$

Hence, computing the right-hand side of Equation (8) is not straightforward.

The seven-point Gauss quadrature [8] is used to approximate the first term in the right-hand side of (16).

To approximate the second term in the right-hand side of (16), we change the variables to polar coordinates, integrating first with respect to $r$ and then using a trapezoidal rule with respect to $\theta$.

Let $K_\ell = S_0 S_1 S_2$, $\ell \in \mathcal{I}_{K,x_0}$ be a triangle s.t. $S_i = (x_i, y_i)$ (resp. $S_i = (R(\Theta_i), \Theta_i)$ in Cartesian (resp. polar) coordinates for $i = 1, 2$.

Let $(X_i, Y_i) := (x_i - 0.5, y_i - 0.5)$, $(X, Y) := (x - 0.5, y - 0.5)$. We aim to evaluate $\int_{K_\ell} r^\gamma f(\theta) X^k Y^k d\mathbf{x}$, for $\gamma \in \{\alpha - 2, \beta - 1\}$ and $k, k' \in \mathbb{N}$ and where $f(\theta) = \sin \theta$ or $\cos \theta$. We need to compute the intersection of lines $\mathcal{D}_{12} := (S_1 S_2)$ and $\mathcal{D}_\theta = (S_0, \mathbf{e}_r)$:

* In the case where $X_1 \neq X_2$, the equations of $\mathcal{D}_{12}$, $\mathcal{D}_\theta$ are s.t.:

$$ \mathcal{D}_{12} : Y = \frac{Y_1 - Y_2}{X_1 - X_2} X + \frac{X_1 Y_2 - X_2 Y_1}{X_1 - X_2}; \quad \mathcal{D}_\theta : Y = \tan \theta X. $$

$\mathcal{D}_{12}$ and $\mathcal{D}_\theta$ intersect at $(X_\theta, Y_\theta)$ s.t.:

$$ X_\theta = \frac{X_1 Y_2 - X_2 Y_1}{\mathbf{e}_\theta \cdot \overrightarrow{S_1 S_2}} \cos \theta, \quad Y_\theta = \frac{X_1 Y_2 - X_2 Y_1}{\mathbf{e}_\theta \cdot \overrightarrow{S_1 S_2}} \sin \theta, $$

$$ \text{Let } R_\theta := (X_\theta^2 + Y_\theta^2)^{\frac{1}{2}} = \frac{X_1 Y_2 - X_2 Y_1}{\mathbf{e}_\theta \cdot \overrightarrow{S_1 S_2}}. $$
• In the case where $X_1 = X_2$ the equations of $\mathcal{D}_{12}$, $\mathcal{D}_{\theta}$ are s.t.:

$$ \mathcal{D}_{12} : X = X_1 \quad \mathcal{D}_{\theta} : Y = \tan \theta X. $$

$\mathcal{D}_{12}$ and $\mathcal{D}_{\theta}$ intersect at $(X_{\theta}, Y_{\theta})$ s.t.: $X_{\theta} = X_1$, $Y_{\theta} = \tan \theta X_1$.

$$ \text{Let } R_{\theta} := (X_{\theta}^{2} + Y_{\theta}^{2})^{\frac{1}{2}} = \frac{X_{1}}{|\cos \theta|}. $$

We have to evaluate:

$$ \int_{K_{\ell}} r^{\gamma} f(\theta) X^{k} Y^{k'} d\mathbf{x} = \int_{\theta=\Theta_{1}}^{\Theta_{2}} \int_{r=0}^{R(\theta)} r^{\gamma} f(\theta) X(r, \theta)^{k} Y(r, \theta)^{k'} r dr f(\theta) d\theta. $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Triangle s.t. $X_1 \neq X_2$.</figcaption></figure>

Let $\gamma' = \gamma + k + k'$ and $\tilde{f}(\theta) = (\cos \theta)^k (\sin \theta)^{k'} f(\theta)$.

$$ \text{We get: } \int_{K_{\ell}} r^{\gamma} f(\theta) X^{k} Y^{k'} d\mathbf{x} = \frac{1}{\gamma' + 2} \int_{\Theta_{1}}^{\Theta_{2}} R(\theta)^{\gamma+2} \tilde{f}(\theta) d\theta. $$

Let $N_{\theta} \in \mathbb{N}^{*}$, $\theta_{0} = \Theta_{1}$, $\theta_{N_{\theta}+1} = \Theta_{2}$, and $\Delta\theta = \frac{\Theta_{2}-\Theta_{1}}{N_{\theta}+1}$. Using the following trapezoidal rule, it stands:

$$ \int_{\theta=\Theta_1}^{\Theta_2} R(\theta)^{\gamma+2} \tilde{f}(\theta) d\theta \approx \sum_{n=0}^{N_{\theta}} \frac{\Delta\theta}{2} \left( R(\theta_{n+1})^{\gamma+2} \tilde{f}(\theta_{n+1}) + R(\theta_n)^{\gamma+2} \tilde{f}(\theta_n) \right). $$

## 6 Numerical Results

The results presented here are based on a in-house code that will be shared soon. Denote by $\mathbf{P}_{dg}^{ku} - \mathbf{P}_{dg}^{kp}$ the DG scheme s.t. $(\mathbf{u}_h, p_h) \in \mathbf{P}_{dg}^{ku} \times \mathbf{P}_{dg}^{kp}$.

Call $\Pi_h$ (resp. $\pi_h$) representing the $L^2$ orthogonal projection onto $\mathbf{X}_h$ (resp. $L_h$). The errors are given by:

$$ \varepsilon_{0}^{\nu}(\mathbf{u}_{h}) := \frac{\|\Pi_{h}\mathbf{u} - \mathbf{u}_{h}\|_{L^{2}(\Omega)}}{\|\mathbf{u}, p\|_{X,v}} \quad \text{and} \quad \varepsilon_{0}^{\nu}(p_{h}) := \frac{v^{-1}\|\pi_{h}p - p_{h}\|_{L^{2}(\Omega)}}{\|\mathbf{u}, p\|_{X,v}} $$

The convergence rates of $\varepsilon_0^\nu(\mathbf{u}_h)$ and $\varepsilon_0^\nu(p_h)$ are denoted by $\tau_\mathbf{u}$ and $\tau_\mathbf{p}$.
## 6.1 Validation of the theoretical convergence rates

The results are given here for $\nu = 1$.

For the convex domain, consider Problem (1) with solution (14) and $(\alpha, \beta) = (0.7, -0.3)$. Then $(\mathbf{u}, p) \in \mathbf{H}^{1+0.7}(\Omega) \times \mathbf{H}^{0.7}(\Omega)$, hence $(\tau_{\mathbf{u}}, \tau_{p}) = (1 + 0.7, 0.7)$.

For the non-convex domain, consider Problem (1) with $\Omega = (0, 1)^2 \setminus [1/2, 1]^2$ and $(\mathbf{u}, p) \in \mathbf{H}^{1+\alpha}(\Omega) \times \mathbf{H}^{\alpha}(\Omega)$ given in [1]. It is studied in [9, 10], with $\alpha \approx 0.54$, $\alpha_{\pm} = \alpha \pm 1$ and $\omega = 3\pi/2$:

$$ 
\left\{ 
\begin{aligned}
\mathbf{u}(r, \theta) &= r^{\alpha} \left( \begin{array}{l} \cos(\theta)\psi'(\theta) + \alpha_{+} \sin(\theta)\psi(\theta) \\ \sin(\theta)\psi'(\theta) - \alpha_{+} \cos(\theta)\psi(\theta) \end{array} \right), \text{ where:} \\
p(r, \theta) &= -r^{\alpha_{-}} (\alpha_{-})^{-1} (\alpha_{+})^{2} \psi'(\theta) + \psi^{(3)}(\theta)
\end{aligned} 
\right. 
\quad (17) 
$$

$\psi(\theta) = (\alpha_{+})^{-1} \sin(\alpha_{+} \theta) \cos(\alpha \omega) - \cos(\alpha_{+} \theta) + (\alpha_{-})^{-1} \sin(\alpha_{-} \theta) \cos(\alpha \omega) + \cos(\alpha_{-} \theta)$. We expect $\tau_{\mathbf{u}} = 2 \alpha \approx 1.09$ and $\tau_{p} = \alpha$.

Let us consider the $\mathbf{P}_{dg}^1 - \mathbf{P}_{dg}^1$ scheme. Figures 2b and 2a represent $\varepsilon_0^\nu(\mathbf{u}_h)$ and $\varepsilon_0^\nu(p_h)$ against the mesh size. Observe that the errors are lower when the domain is convex.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Plots of $\varepsilon_0^\nu(\mathbf{u}_h)$ and $\varepsilon_0^\nu(p_h)$ against the mesh size.</figcaption></figure>

Convergence rates between two consecutive meshes are recorded in Table 1. In both test-cases, we observe a good agreement with the theoretical ones (cf. Th. 3.1).

<table><caption>Table 1: Convergence rates computed between two consecutive meshes.</caption><thead><tr><th>test-case</th><th>τ</th><th>theory</th><th>mesh 1-2</th><th>mesh 2-3</th><th>mesh 3-4</th><th>mesh 4-5</th></tr></thead><tbody><tr><td rowspan="2">(14) (Ω convex)</td><td>τ<sub>u</sub></td><td>1.7</td><td>1.63</td><td>1.66</td><td>1.67</td><td>1.68</td></tr><tr><td>τ<sub>p</sub></td><td>0.7</td><td>0.78</td><td>0.74</td><td>0.73</td><td>0.71</td></tr><tr><td rowspan="2">(17) (Ω non-convex)</td><td>τ<sub>u</sub></td><td>1.09</td><td>1.02</td><td>0.94</td><td>1.08</td><td>1.08</td></tr><tr><td>τ<sub>p</sub></td><td>0.54</td><td>0.85</td><td>0.59</td><td>0.59</td><td>0.53</td></tr></tbody></table>
## 6.2 Comparing different DG schemes

Consider test-case (14) with $v = 10^{-6}$ which is close to the kinematic viscosity of the pressurized water in nuclear reactor core and $(\alpha, \beta) = (6, 3)$. In order to compare DG schemes with different pairs of polynomial orders ($k_u, k_p$), we plot the error $\varepsilon_0^v(\mathbf{u}_h)$ (resp. $\varepsilon_0^v(p_h)$) against the mesh size on Fig. 3 or against the CPU time on Fig. 4.

* $\mathbf{P}_{dg}^1 - P_{dg}^1$ vs $\mathbf{P}_{dg}^2 - P_{dg}^1$: for a given CPU time, the error $\varepsilon_0^v(\mathbf{u}_h)$ and the error $\varepsilon_0^v(p_h)$ are similar for both schemes. However, using the $\mathbf{P}_{dg}^1 - P_{dg}^1$ scheme requires less memory footprint and seems therefore more efficient.

* $\mathbf{P}_{dg}^1 - P_{dg}^2$ vs $\mathbf{P}_{dg}^2 - P_{dg}^2$: for a given CPU time, the error $\varepsilon_0^v(p_h)$ are similar for both schemes, and the error $\varepsilon_0^v(\mathbf{u}_h)$ is better using the $\mathbf{P}_{dg}^1 - P_{dg}^2$ scheme, which therefore seems more efficient.

<figure><img src="image_5.png" /><figcaption>Figure 3: Plots of $\varepsilon_0^v(\mathbf{u}_h)$ and $\varepsilon_0^v(p_h)$ against the mesh size.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Plots of $\varepsilon_0^v(\mathbf{u}_h)$ and $\varepsilon_0^v(p_h)$ for $v = 10^{-6}$ against CPU time (s).</figcaption></figure>
Observe that for a given CPU time, the $P_{dg}^1 - P_{dg}^2$ scheme yields better results. Additionally, it exhibits the fastest convergence rate.

# 7 Conclusion

This study explored DG schemes for solving the Stokes problem. Through analysis and numerical experiments, we demonstrated their efficiency and accuracy even with low-regularity solutions. The $P_{dg}^1 - P_{dg}^2$ scheme seems more efficient, which, to our knowledge, is a new result. We provided robust a priori error estimates and we proposed to solve the linear system with the help of a domain decomposition method (DDM). Further, we plan to challenge the DDM in parallel and to implement Oseen problem, using non polynomial basis functions.

# References

[1] M. Dauge, Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. Part I: Linearized equations, SIAM J. Math. Anal. **20** (1989).

[2] D. A. Di Pietro, A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Springer, Mathématiques & Applications (2012).

[3] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations, Springer-Verlag (1986).

[4] E. Jamelot, M. Mroueh and P. Omnes, A Priori error estimates for Discontinuous Galerkin Finite Element discretization of the Stokes problem with low regularity (in progress).

[5] E. Jamelot, cea-trust-platform Stokes_NCFEM https://github.com/cea-trust-platform/Stokes_NCFEM (2022).

[6] E. Di Nezza, G. Palatucci and E. Valdinoci, Hitchhiker's guide to the fractional Sobolev spaces, Bull. Sci. Math **136**, 521–573 (2012).

[7] B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation, SIAM (2008).

[8] G. Vinsard, Modélisation/Calculs scientifiques - GME5 Leçon 3 Intégration numérique (2017).

[9] R. Verfürth, A Posteriori Error Estimators for the Stokes Equations, Numerische Mathematik Springer-Verlag (1989).

[10] M. Cermak, F. Hecht, Z. Tang and M. Vohralík, Adaptive inexact iterative algorithms based on polynomial-degree-robust a posteriori estimates for the Stokes problem, Springer-Verlag GmbH Deutschland (2017).

[11] C. Moldoveanu, Simulation des grandes échelles de tourbillons longitudinaux soumis à une turbulence extérieure intense, PhD, Institut de Mécanique des Fluide de Toulouse, France (2007).