tensor-product-stokes.tex

\documentclass{article}

\usepackage{amsmath}
%\usepackage{amsfonts}
\usepackage{amsthm}
%\usepackage{amssymb}
%\usepackage{mathrsfs}
%\usepackage{fullpage}
%\usepackage{mathptmx}
%\usepackage[varg]{txfonts}
\usepackage{color}
\usepackage[charter]{mathdesign}
\usepackage[pdftex]{graphicx}
%\usepackage{float}
%\usepackage{hyperref}
%\usepackage[modulo, displaymath, mathlines]{lineno}
%\usepackage{setspace}
%\usepackage[titletoc,toc,title]{appendix}
\usepackage{natbib}

%\linenumbers
%\doublespacing

\theoremstyle{definition}
\newtheorem*{defn}{Definition}
\newtheorem*{exm}{Example}

\theoremstyle{plain}
\newtheorem{thm}{Theorem}
\newtheorem{lem}{Lemma}
\newtheorem{prop}{Proposition}
\newtheorem{cor}{Corollary}
\newtheorem{asm}{Assumption}

\newcommand{\argmin}{\text{argmin}}
\newcommand{\ud}{\hspace{2pt}\mathrm{d}}
\newcommand{\bs}{\boldsymbol}
\newcommand{\PP}{\mathsf{P}}
\let\divsymb=\div % rename builtin command \div to \divsymb
\renewcommand{\div}[1]{\operatorname{div} #1} % for divergence
\newcommand{\Id}[1]{\operatorname{Id} #1}

\title{Tensor product elements for the Stokes equations}
\author{Daniel R. Shapero}
\date{}

\begin{document}

\maketitle

% --------------------
\section{Introduction}

The goals of this paper are (1) to prove the inf-sup stability of certain finite element discretizations of the Stokes equations on \emph{extruded domains} in $\mathbb{R}^3$ and (2) to demonstrate the use of these elements on several real problems.
An extruded domain $\Omega$ is one that is obtainable as a mapping of the Cartesian product of a \emph{footprint} domain $\Phi$ in $\mathbb{R}^2$ with the unit interval $[0, 1]$.
We might not have that $\Omega = \Phi \times [0, 1]$ exactly, but rather that $\Omega$ is obtainable by applying some vertical stretching transformation to the final coordinate.
For example, we can get such domains by using \emph{terrain-following} coordinates.

Extruded domains are very common in the geophysical sciences.
For example, flow in the ocean, atmosphere, and of large ice sheets can all be described as extruded domains.
Moreover all of these systems are \emph{thin-film} flows, where the horizontal length scales are much larger than the vertical length scales by a factor of 20 or more.
These kinds of flows are often very smooth in the vertical dimension, and it may be more computationally efficient to discretize a problem using high-degree vertical basis functions in only a few layers than it is to use many vertical layers with low-degree basis functions.

In the following, we will express the Stokes equations as a linearly-constrained quadratic minimization problem.
Let $V$ and $Q$ be Hilbert spaces, $A : V \to V^*$ and $B : V \to Q^*$ bounded linear operators with $A$ symmetric and positive-definite.
The general form of constrained optimization problem that we'll consider here is to find a critical point $u$, $p$ of the \emph{Lagrangian}
\begin{equation}
    L(u, p) = \frac{1}{2}\langle Au, u\rangle - \langle f, u\rangle + \langle Bu - g, p\rangle.
\end{equation}
For the Stokes equations, these operators are
\begin{equation}
    \langle Au, v\rangle = \int_\Omega \mu\dot\varepsilon(u):\dot\varepsilon(v)\; dx
\end{equation}
where $\mu$ is the (strictly positive) viscosity coefficient and $\dot\varepsilon(u) \equiv \frac{1}{2}(\nabla u + \nabla u^*)$ is the strain rate tensor, and
\begin{equation}
    \langle Bv, q\rangle = \int_\Omega (\nabla\cdot v)\,q\;dx.
\end{equation}
The constraint that $Bu - g = 0$ enforces the physical condition that the fluid is incompressible or has prescribed sources of mass.

We will consider shape-regular triangulations of the domain $\Omega$ with a maximum cell diameter $h$.
Let $V_h$, $Q_h$ denote a pair of finite-dimensional subspaces of $V$ and $Q$ respectively.
The crucial point we need to prove is that the approximation spaces we choose for the velocity and pressure satisfy the \emph{Ladyzhenskaya-Babu\v{s}ka-Brezzi} (LBB) or \emph{inf-sup} condition:
\begin{equation}
    \inf_{q\in Q_h}\sup_{v \in V_h}\frac{\langle Bv, q\rangle}{\|v\|_V\|q\|_Q} \ge \beta > 0
\end{equation}
for some $\beta$ that is independent of the mesh size $h$ \citep{boffi2013mixed}.
Another way of stating this condition is that
\begin{equation}
    \|B^*q\|_{V^*} \ge \beta\|q\|_Q
\end{equation}
for all $q$ in $Q_h$, where the norm on the left-hand side is an operator norm in the dual space $V^*$.
This latter condition is another way of saying that $B^*$ has a bounded inverse map and that $\|B^{-*}\| \le \beta^{-1}$.

\textbf{To our knowledge, the LBB-stability of tensor product elements with high-degree polynomial vertical basis functions has not been established in the computational fluid dynamics literature.}
\citet{canuto1984combined} established the LBB-stability of a basis using the tensor product of the MINI element in the horizontal direction and a Fourier basis with an equal number of modes for the velocity and pressure in the vertical direction.
Their proof of LBB-stability crucially relied on special properties of the Fourier basis and does not generalize to other element families.
\citet{nakahashi1989finite} studied prismatic elements for the Navier-Stokes equations but did not address LBB stability.
\citet{isaac2015solution} explored higher-order elements and multigrid solvers for the Stokes equations in very anisotropic domains.
This work used equal polynomial degree in the vertical and horizontal for the velocity space.

We propose to use velocity and pressure spaces for the 3D problem that are fashioned out of stable elements for the 2D problem.
Certain choices of 2D basis will give a stable and conforming basis for the 3D problem, like the MINI, Taylor-Hood, and Crouzeix-Raviart elements.
This list is not exhaustive and other elements may be amenable to our approach.
\textcolor{red}{Do we want to try RT or other $H(\text{div})$ elements?}
Our approach differs from that of \citet{isaac2015solution} in our use of higher degree polynomials in the vertical direction than in the horizontal.


\section{Elements}

To prove the inf-sup stability of tensor product elements for the Stokes equations on extruded domains, we will proceed through a more general route.
Rather than start from the assumption that the domain is the product of a footprint domain in 2D and the unit interval, we will instead suppose that the domain $\Omega$ is equal to a product of domains $\Omega_x$ and $\Omega_z$ of arbitrary dimensions $d_x$ and $d_z$.
We will refer to these as the ``horizontal'' and ``vertical'' domains respectively.
Our approach can be summarized as building up stable elements on the product domain through a certain combination of stable elements on the factor domains.

The key observation enabling the constructions that follow is that \textbf{the product structure of the spatial domain implies a certain product structure for the space of vector fields defined on that domain, and for the divergence operator.}
To see how the spaces and operators can be decomposed, first define the pairs of spaces $V_x = H^1(\Omega_x, \mathbb{R}^{d_x})$, $Q_x = L^2(\Omega_x)$ of vector and scalar fields on the horizontal domain and likewise $V_z$, $Q_z$ for the vertical domain.
Additionally, we let $S_x = H^1(\Omega_x)$ and $S_z = H^1(\Omega_z)$ be spaces of scalar fields on each product domain.
\begin{lem} The Sobolev space $V = H^1(\Omega_x\times\Omega_z, \mathbb{R}^{d_x + d_z})$ of velocity fields on the product domain $\Omega_x \times \Omega_z$ can be decomposed as the direct sum of tensor products as follows:
\begin{equation}
    V = \left(V_x \otimes S_z\right) \oplus \left(S_x \otimes V_z\right).
    \label{eq:velocity-space-decomposition}
\end{equation}
\end{lem}

Knowing that $V$ can be written as a direct sum of tensor products, we can then ask how the divergence operator acts on each summand and factor.
Define the operator $R_x : V \to V_x \otimes S_z$ that selects the horizontal components of the velocity field and likewise $R_z : V \to S_x\otimes V_z$ the operator that selects the vertical components.
We then define their partial inverses, the injection operators $R_x^*$, $R_z^*$ from $V_x\otimes S_z$ or $S_x\otimes V_z$ into $V$ respectively.
Let $\Id_H$ denote the identity map on a Hilbert space $H$.
The key result that helps us construct stable elements for the Stokes problem is:
\begin{lem} If $B : V \to Q^*$ is the divergence operator acting on the space of Sobolev vector fields on $\Omega_x\times\Omega_z$, then
\begin{equation}
    B = (B_x \otimes \Id_{S_z})R_x + (\Id_{S_x}\otimes B_z)R_z.
    \label{eq:b-decomposition}
\end{equation}
\end{lem}

For example, the 3D divergence operator can be written as the sum of (1) a 2D divergence operator acting only on the horizontal components of the velocity field that ignores their dependence on the vertical coordinate and (2) a 1D divergence operator acting only on the vertical components that ignores their dependence on the horizontal coordinate.

The above discussion has operated at the level of the idealized problem; we now proceed to the discretization.
We can prove that the inf-sup condition holds for a pair $V^h$, $Q^h$ of discrete spaces if there exists a \emph{compatible interpolation operator} \citep{boffi2013mixed}.
A $B$-compatible interpolation operator $\Pi^h$ is a linear operator from $V$ to $V^h$ such that, for all $v$ in $V$ and $q$ in $Q^h$,
\begin{equation}
    \langle B(v - \Pi^hv), q\rangle = 0,
\end{equation}
and $\|\Pi^hv\| \le C\|v\|$ where the constant $C$ is independent of the mesh size $h$.

We now introduce two pairs of discrete spaces $V_x^h$, $Q_x^h$ and $V_z^h$, $Q_z^h$.
We will also need two discrete spaces $S_x^h$, $S_z^h$ that discretize scalar fields in $S_x$ and $S_z$ respectively.
We can then decompose the space of discrete velocity fields as follows:
\begin{equation}
    V^h = \left(V_x^h \otimes S_z^h\right)\oplus\left(S_x^h\otimes V_z^h\right)
\end{equation}
by analogy with equation \eqref{eq:velocity-space-decomposition}.
The final pieces are the following
\begin{asm} There is a $B_x$-compatible interpolation operator $\Pi_x^h : V_x^h \to V_x$ and a $B_z$-compatible $\Pi_z^h : V_z^h \to V_z$.
    \label{asm:b-compatible-operators}
\end{asm}
\begin{asm} There is an operator $\Psi_x^h : S_x \to S_x^h$ such that
    \begin{equation}
        \langle \phi - \Psi_x^h\phi, q\rangle = 0
    \end{equation}
    for all $\phi$ in $S_x$ and $q \in Q_x^h$ and likewise $\Psi_z^h : S_z \to S_z^h$, $Q_z^h$.
    \label{asm:ortho-projection}
\end{asm}

This latter assumption holds if we can guarantee that the dimension of $S_\circ^h$ is greater than or equal to $Q_\circ^h$.
In concrete examples, we will show cases where $Q_\circ^h \subset S_\circ^h$ and the assumption is guarantee trivially by taking $\Psi_\circ^h$ to be the orthogonal projection operator.
Now all the pieces are in place to complete our construction.
\begin{thm}\label{thm:main-theorem} The operator $\Pi^h : V \to V^h$ defined by
\begin{equation}
    \Pi^h = R_x^*(\Pi_x^h\otimes\Psi_z^h)R_x + R_z^*(\Psi_x^h\otimes\Pi_z^h)R_z
\end{equation}
    is a $B$-compatible interpolation operator if $B$ has the structure of equation \eqref{eq:b-decomposition} and Assumptions \ref{asm:b-compatible-operators} and \ref{asm:ortho-projection} hold.
    \proof First, observe that $R_x^*R_x + R_z^*R_z = \Id_V$.
    We then have that
    \begin{equation}
        \Id_V - \Pi^h = R_x^*(\Id_{V_x}\otimes\Id_{S_z} - \Pi_x^h\otimes\Psi_z^h)R_x + R_z^*(\Id_{S_x}\otimes\Id_{V_z} - \Psi_x^h\otimes\Pi_z^h)R_z
    \end{equation}
    using the definition of $\Pi^h$.
    Next, we'll need the fact that $R_zR_x^* = R_xR_z^* = 0$.
    We can then compute
    \begin{align}
        B(\Id_V - \Pi^h) & = (B_x\otimes \Id_{S_z})(\Id_{V_x}\otimes\Id_{S_z} - \Pi_x\otimes\Psi_z)R_x \nonumber\\
        & \qquad + (\Id_{S_x}\otimes B_z)(\Id_{S_x}\otimes\Id_{V_z} - \Psi_x\otimes\Pi_z)R_z \\
        & = \left(B_x\otimes\Id_{S_z} - B_x\Pi_x\otimes\Psi_z\right)R_x \nonumber \\
        & \qquad + (\Id_{S_x}\otimes B_z - \Psi_x\otimes B_z\Pi_z)R_z
        \label{eq:b-id-pi}
    \end{align}
    where we have used that $(A_1\otimes A_2)\cdot(B_1\otimes B_2) = (A_1\cdot B_1)\otimes(A_2\cdot B_2)$.
    Now we let $u$ be an element of $V$ and $q$ an element of $Q^h$.
    We can assume that
    \begin{equation}
        v = \left(\begin{matrix}v_x(x)\phi_z(z) \\ \phi_x(x)v_z(z)\end{matrix}\right), \quad q = q_x(x)q_z(z)
    \end{equation}
    where $v_x \in V_x$, $\phi_z \in S_z$, $\phi_x \in S_x$, $v_z \in V_z$, and $q_x \in Q_x^h$, $q_z \in Q_z^h$.
    We can extend to the general case by linearity and taking limits.
    Substituting this into equation \eqref{eq:b-id-pi},
    \begin{align}
        \langle B(\Id_V - \Pi^h)v, q\rangle & = \langle(B_x\otimes\Id_{S_z} - B_x\Pi_x\otimes \Psi_z)(v_x\phi_z), q_x q_z\rangle \nonumber \\
        & \qquad + \langle(\Id_{S_z}\otimes B_z - \Psi_x\otimes B_z\Pi_z)(\phi_xv_z), q_xq_z\rangle.
        \label{eq:b-id-pi-vq}
    \end{align}
    Focusing just on the first term, we can distribute the tensor products of operators over those of functions:
    \begin{align}
        & \langle(B_x\otimes\Id_{S_z} - B_x\Pi_x\otimes \Psi_z)(v_x\phi_z), q_x q_z\rangle \nonumber\\
        &\qquad = \langle (B_xv_x)\phi_z - (B_x\Pi_xv_x)(\Psi_z\phi_z), q_x q_z\rangle \nonumber\\
        &\qquad = \langle B_xv_x, q_x\rangle\langle\phi_z, q_z\rangle - \langle B_x\Pi_xv_x, q_x\rangle\langle \Psi_z\phi_z, q_z\rangle \\
        &\qquad = \langle B_xv_x, q_x\rangle\langle \phi_z - \Psi_z\phi_z, q_z\rangle \\
        &\qquad = 0
    \end{align}
    where we apply the definition of inner products in the tensor product of spaces, the fact that $\Pi_x$ is $B_x$-compatible, and finally the Assumption \ref{asm:ortho-projection} for $\Psi_z$.
    The same reasoning shows that the second term of equation \eqref{eq:b-id-pi-vq} is zero as well.\qed
\end{thm}

Now we can look at the case of interest to us, where $\Omega_x$ is a domain in $\mathbb{R}^2$ and $\Omega_z$ is the unit interval $[0, 1]$.
We will denote by $CG_k$ the space of continuous Galerkin basis functions on triangles of degree $k$, $DG_k$ the discontinuous Galerkin basis, and $B_k$ the space of bubble functions.
In the following, we'll consider a few different element families:
\begin{enumerate}
    \item MINI: $V_x^h = (CG_1 \oplus B_3)^2$, $Q_x^h = CG_1$
    \item Taylor-Hood: $V_x^h = {CG_2}^2$, $Q_x^h = CG_1$
    \item Crouzeix-Raviart: $V_x^h = (CG_2\oplus B_3)^2$, $Q_x^h = DG_1$
\end{enumerate}
All of these element pairs are stable for the 2D Stokes equations and, with appropriate choices of (1) the vertical spaces $V_z$, $Q_z$ and (2) the scalar extension spaces $S_z^h$ and $S_x^h$, we can apply Theorem \ref{thm:main-theorem}.

When the vertical space is the interval $[0, 1]$, the divergence operator reduces to just differentiation $\partial_z$.
In that case, there is a $\partial_z$-compatible interpolation operator $\Pi_z^h$ for the vertical spaces $V_z^h = CG_{m + 1}$, $Q_z^h = DG_m$ for any polynomial degree $m$.
(Note that the 1D Stokes equations as such are trivial; the divergence-free constraint in 1D forces all velocity solutions to be constant.
We can still ask whether there are 1D discrete function spaces $V_z^h$, $Q_z^h$ for which there is a compatible interpolation operator for $\partial_z$.
The existence of such discrete spaces and interpolation operators is all that matters.)

The final ingredient is to pick the discrete space $S_z^h$ that we use to extend the horizontal velocities in the vertical and the discrete space $S_x^h$ that extends the vertical velocities in the horizontal.
In the examples above -- MINI, Taylor-Hood, and Crouzeix-Raviart -- the velocity space can always be expressed as $V_x^h = K^2$ for some scalar space $K$.
For simple velocity spaces like these, we can take
\begin{align*}
    S_x^h & = K \\
    S_z^h & = V_z^h
\end{align*}
because elements of the space $V_z^h$ are scalar fields.
Having made this choice, we then find that
\begin{align*}
    V^h & = (V_x^h\otimes S_z^h) \oplus (S_x^h \otimes V_z^h) \\
    & = (K^2\otimes V_z^h) \oplus (K\otimes V_z^h) \\
    & \cong (K \otimes V_z^h)^3.
\end{align*}
This choice trivially guarantees the existence of the projection operators of Assumption \ref{asm:ortho-projection}.
We can thus assert the following
\begin{cor} \label{cor:stable-pairs} The following elements are stable for the 3D Stokes problem on an extruded domain:
    \begin{enumerate}
        \item $V = ((CG_1 \oplus B_3) \otimes CG_{k + 1})^3$, $Q = CG_1\otimes DG_k$ 
        \item $V = (CG_2 \otimes CG_{k + 1})^3$, $Q = CG_1\otimes DG_k$
        \item $V = ((CG_2\oplus B_3) \otimes CG_{k + 1})^3$, $Q = DG_1\otimes DG_k$
    \end{enumerate}
\end{cor}

Two final properties are worth noting.
First, we have used $CG$ and $DG$ to denote the continuous and discontinuous Galerkin spaces.
This notation presumes that we are using the usual Lagrange finite element bases, where the degrees of freedom are the values of the basis function at interpolation points.
One can equivalently use the modal Gauss-Lobatto-Legendre (GLL) and Gauss-Legendre (GL) bases in the vertical direction, where the basis functions are mutually orthogonal.

Second, observe that the pressure space for the Crouzeix-Raviart element in 2D contains the piecewise-constant functions.
As a consequence, the velocity solution conserves mass element-wise.
Devising inf-sup stable pairs for the 3D Stokes problem on arbitrary triangular meshes is much more difficult than in 2D.
The pressure space for the last element pair in Corollory \ref{cor:stable-pairs}, defined by extruding the 2D Crouzeix-Raviart element into 3D, also has the desirable property of containing the piecewise constant functions.

There are element families that do not have a simple expression as a product of scalar spaces.
For example, \citet{cockburn2007note} discretize the velocities in the related Oseen problem using Raviart-Thomas elements, which makes the velocities exactly divergence-free.
The Raviart-Thomas space cannot be expressed as a product of scalar spaces.



% ----------------------
\section{Demonstrations}

The computational demonstrations in this section all use the Firedrake software package \citep{FiredrakeUserManual}.
Firedrake includes built-in support for tensor product meshes and elements, as well as sophisticated loop optimizations such as sum-factorization that accelerate computations on extruded meshes \citep{mcrae2016automated}.

\subsection{Empirical confirmation of inf-sup stability}

\begin{figure}
    \begin{center}
        \includegraphics[width=0.6\linewidth]{demo/lbb-stability/results-2d.pdf}
    \end{center}
    \caption{Empirical inf-sup constants for the three 2D element families on the unit square as the mesh is refined.
    These elements are known to stable for the 2D Stokes problem.}
    \label{fig:empirical-lbb-2d}
\end{figure}

Given a mesh of the spatial domain and a pair of discrete spaces $V^h$, $Q^h$, we can check empirically whether they satisfy the inf-sup condition by solving a generalized eigenvalue problem \citep{rognes2012automated}.
Let $M$ and $N$ be the Riesz representers for the canonical mapping from $V \to V^*$ and $Q \to Q^*$ respectively.
Consider the following eigenproblem: find velocity-pressure pairs $u$, $p$ and scalars $\lambda$ such that
\begin{align}
    Mu + B^*p & = 0 \\
    Bu & = -\lambda Np
\end{align}
The eigenvalues of this problem are all real and positive, and the square root of the smallest eigenvalue is equal to the inf-sup constant for $V^h$, $Q^h$ \citep{malkus1981eigenproblems}.
We can form this eigenproblem in Firedrake and then call out to the SLEPc package to solve it \citep{hernandez2005slepc}.
This empirical check is not a substitute for a formal proof of inf-sup stability, but rather serves as a ``smoke test'' that our numerical implementation and our proofs are not obviously wrong.

We performed the empirical inf-sup stability check using the unit square as our footprint domain, with a mesh spacing starting at 1/4 and going down to 1/32.
Figure \ref{fig:empirical-lbb-2d} shows the results for 2D elements which we know to be inf-sup stable.
\textcolor{red}{Do the 3D elements.}


\subsection{Lid-driven cavity flow}

\subsection{Flow around a cylinder}

The footprint domain for this problem is the rectangle $[-L_x, +L_x] \times [-L_y, +L_y]$ with a circle of radius $R$ centered at the point $(0, \ell_y)$ removed.
The footprint is extruded to a depth of 1.
We set the pressure to be equal to 1 along the inflow boundary $y = -L_y$ and equal to 0 along the outflow obundary $y = +L_y$.
On the top, bottom, side walls, and cylindrical exclusion, we enforce no-slip boundary conditions on the velocity.
The results are shown in Figure \ref{fig:cylinder-flow}.

\begin{figure}
    \begin{center}
        \includegraphics[width=0.6\linewidth]{demo/flow-around-cylinder/cylinder_flow.pdf}
    \end{center}
    \caption{Mesh for flow around cylinder problem (left) and quiver plot of the depth-averaged velocity field (right).}
    \label{fig:cylinder-flow}
\end{figure}

As a smoke test, we evaluate the total flux of material into and out of the domain.
The sum of the influx and outflux should be zero for an incompressible fluid.
The relative error between the total flux and the influx for different element choices is shown in \textcolor{red}{add a table}.
The extruded Crouzeix-Raviart element has a much smaller relative flux error than the MINI and Taylor-Hood elements.
The Crouzeix-Raviart pressure space is discontinuous and therefore preserves mass element-wise rather than only in a global sense.
Nonetheless, there is no reason a priori why it should also have much lower global flux error than the alternatives \textcolor{red}{or is there?}

\subsection{Thin-film flow with terrain-following coordinates}

The preceding examples considered cases where the spatial domain is equal to $\Omega \times [0, 1]$ for some footprint domain $\Omega$.
Here we will consider a different case: the domain is the region between two smooth surfaces $b$ and $s$, respectively the bed and surface, that both map $\Omega$ to the reals, with the thickness $h = s - b$ always positive.
Rather than solve on the physical domain, we can introduce the \emph{terrain-following coordinate}
\begin{equation}
    \zeta = \frac{z - b}{h}
\end{equation}
which now takes values in $[0, 1]$.
The computational domain is then $\Omega \times [0, 1]$ but we need to introduce a metric tensor and alter the usual differential operators.



% ------------------
\section{Discussion}

\pagebreak

\bibliographystyle{plainnat}
\bibliography{tensor-product-stokes.bib}

\end{document}