text:
\section{Introduction}

Stokes flow describes the motion of an incompressible viscous fluid at slow speeds, or small scales, and can be used to model a wide range of intricate phenomena, including mantle dynamics, the swimming of micro-organisms, the sedimentation of particulates, and the flotation of water droplets in clouds. In the steady-state case, the corresponding governing equations of motion are given by the Stokes equations, which generally take on one of two forms: either
\begin{equation} \label{eq:stokes1} \begin{aligned} -\mu \nabla^2 \mathbf u + \nabla p &= {\mathbf f} \\ \nabla \cdot \mathbf u &= 0, \end{aligned} \end{equation}
or, alternatively,
\begin{equation} \label{eq:stokes2} \begin{aligned} -\nabla \cdot \bigl( \mu (\nabla \mathbf u + \nabla \mathbf u^\mathsf{T}) \bigr) + \nabla p &= {\mathbf f} \\ \nabla \cdot \mathbf u &= 0, \end{aligned} \end{equation}
where $\mu$ specifies the viscosity of the fluid, $\mathbf u$ and $p$ describes its velocity and pressure fields, and $\mathbf f$ specifies the net external forces acting on the fluid. The form given in \eqref{eq:stokes1} is here referred to as the \textit{standard form} of the Stokes equations, and \eqref{eq:stokes2} as the \textit{viscous-stress form}, with the particular choice depending on the end application. For example, the viscous-stress form is generally applicable when $\mu$ is variable or if boundary conditions on stress are imposed.

Our motivation in this work is to develop fast multigrid solvers for computing high-order accurate solutions of the Stokes systems \eqref{eq:stokes1} or \eqref{eq:stokes2}, with extension also to multi-phase variants involving interfacial jump conditions in velocity and stress, as well as to time-dependent problems. In particular, we consider a framework based on local discontinuous Galerkin (LDG) methods \cite{ldg}, and build on prior work developing efficient multigrid algorithms for LDG discretisations of elliptic interface (Poisson-like) problems \cite{dgmg,fluxx}. We show that standard and simple-to-implement geometric multigrid algorithms can be applied to the resulting multi-phase Stokes problems---in the steady-state case, results show that the solver can match the speeds of fast geometric multigrid methods for Poisson problems; in the time-dependent case, convergence rates further accelerate as the effective Reynolds number is increased.

In particular, the presented multigrid algorithm consists of a standard V-cycle using an element-wise block Gauss-Seidel smoother---individual blocks correspond to individual mesh elements, such that the elemental degrees of freedom of both velocity and pressure are collected into the same block.
Key to the rapid convergence of this approach is a suitable choice of the pressure penalty stabilisation parameter underlying the LDG framework---if the parameter is chosen well, then a highly efficient multigrid algorithm is obtained.
We discuss how to choose this parameter for steady-state Stokes problems, and develop a simple strategy for generalising this choice to time-dependent Stokes problems.
Extensive tests of the multigrid methods are presented in this paper, including problems which impose Dirichlet or stress boundary conditions, variable-viscosity problems, test cases with curved geometry using semi-unstructured meshes, and multi-phase embedded interface problems with viscosity and density coefficients exhibiting discontinuities several orders in magnitude.

\subsection*{Previous work}

A vast amount of work in computational science and engineering has been devoted to the efficient solution of Stokes systems and saddle point problems in general; for an in-depth review of the correspondingly wide array of different approaches and their applications, see Benzi, Golub, and Liesen \cite{BenziGolubLiesen2005}. These approaches include, among others: block preconditioner methods, which operate on the viscosity, gradient, and divergence operator block structure of the Stokes equations; Schur complement methods, which manipulate, and usually approximate, the Schur complement of the saddle point system; and stationary iterative methods, such as the well-known Uzawa method, which alternates between updates of velocity and pressure, holding the other fixed, through the two governing equations in \eqref{eq:stokes1}. Here, we briefly review work on multigrid-style methods, particularly schemes based on the coupled solution of velocity and pressure, as is relevant to the present work; see also the reviews \cite{Wesseling2001,WesselingOosterlee2001,OosterleeGaspar2008}.

As mentioned, the multigrid algorithm developed here uses a block Gauss-Seidel relaxation method, with each block collecting the velocity and pressure degrees of freedom on each mesh element. This approach is similar in essence to the ``symmetric coupled Gauss-Seidel method'' of Vanka \cite{Vanka1986,Vanka1986b} and can also be considered as a kind of ``box relaxation'' scheme \cite{BrandtLivne}. Vanka-type smoothers, originally devised for staggered-grid finite difference methods, visit each grid cell, solve for the velocity and pressure unknowns simultaneously via a local Stokes-like problem, and then move onto the next cell. Typically, a damping/under-relaxation parameter is needed to ensure convergence. In the original Vanka method, the Stokes system is restricted to the local variables in each grid cell and off-diagonal entries of the viscous operator are zeroed-out to facilitate a simpler update for the unknowns \cite{Wesseling2001,OlshanskiiTyrtyshnikov}. Variations have led to schemes which include the off-diagonal terms \cite{ThompsonFerziger1989,BorzacchielloLericheBlottiereGuillet2017}, line-based sweeping methods \cite{ThompsonFerziger1989,OosterleeWesseling1993b,OosterleeWesseling1993}, and have been examined with local mode analyses \cite{Sivaloganathan1991}. Vanka-type smoothers may also be considered as iterative Schwarz solvers, whereby the subdomains of the Schwarz method corresponds to the collection of degrees of freedom in each grid cell. Schemes building on this idea have since been developed for finite volume and finite element methods with much of the attention devoted to the choice of relaxation parameters, the choice of subdomains (e.g., whether to use one cell, or patches of cells), and on theoretical proofs of convergence in a multigrid setting, see, e.g., \cite{SchoberlZulehner2003,Manservisi2006,OosterleeGaspar2008,GasparNotayOosterleeRodrigo2014,Chen2015,HeMacLachlan2017,ColeyBenzakenEvans2017}; they have also found application in variable-viscosity Stokes problems \cite{BorzacchielloLericheBlottiereGuillet2017} and in computational solid mechanics \cite{WobkerTurek2009,HongKrausXuZikatanov2016}. As an example, in very recent work, Farrel, He, and MacLachlan \cite{FarrellHeMacLachlan2019} demonstrated the application of local Fourier analysis on these smoothers and found that smaller patches result in better convergence per floating point operation. Meanwhile, solvers specific to discontinuous Galerkin methods of the Stokes equations have also been devised; here, one possible approach is to exactly enforce the divergence constraint across the multigrid hierarchy through manipulation of the DG spaces, e.g., through $H(\text{div},\Omega)$-conforming discretisations. In the associated multigrid solvers, the divergence constraint is built into the coarse and fine mesh approximations, see, e.g., \cite{Chen2015,KanschatMao2015,HongKrausXuZikatanov2016,AdlerBensonMacLachlan2016,ColeyBenzakenEvans2017}.

In many of these works, satisfactory multigrid convergence rates are reported, but they generally do not match the speeds of an efficient geometric multigrid method designed for scalar elliptic equations. In some cases, performance degrades as the mesh is refined, or as viscosity ratios increase, or in the case of time-dependent Stokes problems, as the Reynolds number changes. In contrast, for the LDG schemes devised here, we found that a simple block Gauss-Seidel relaxation method, which does not use any under- or over-relaxation parameters, can result in rapid multigrid convergence across a variety of challenging Stokes problems. Another advantage to a block Gauss-Seidel method is its simple implementation and the possibility of parallelism; for example, Gmeiner \textit{et al} \cite{GmeinerHuberJohnRudeWohlmuth2016} and Bauer \textit{et al} \cite{BauerKlementOberhuberZabka2016} demonstrate massively parallel and GPU implementations, respectively. For example, some of the three-dimensional tests in this paper used half a billion degrees of freedom and scaled to several hundred computing nodes, though we do not report on scaling performance here.

\subsection*{Outline}

In the main article, the central ideas and results are presented, while nonessential details of the LDG discretisation, grid convergence analyses, and implementation possibilities are deferred to the \iffalse Supplementary Material\fi Appendices. \iffalse The remainder of this article is organised as follows. \fi First, we outline the essential components of the LDG framework for the multi-phase Stokes equations. Second, the design of a standard multigrid V-cycle is outlined, afterwhich the role of pressure penalty stabilisation on multigrid efficiency is examined. In the remaining two sections, results are presented for a variety of problems for the steady-state and time-dependent Stokes equations, respectively. We then conclude, summarising the key observations made in this work along with a discussion of future research avenues.

\section{Local Discontinuous Galerkin Methods for Multi-Phase Stokes Problems}

In this work, we build on the LDG schemes developed by Cockburn \textit{et al} \cite{CockburnKanschatSchotzauSchwab2002} and extend them to the variable-viscosity multi-phase Stokes problem. The governing equations are written as follows: we seek to determine a velocity field $\mathbf u : \Omega \to {\mathbb R}^d$ and pressure field $p : \Omega \to {\mathbb R}$ such that
\begin{equation} \label{eq:govern1} \left. \begin{aligned} -\nabla \cdot \bigl(\mu_i (\nabla \mathbf u + \gamma \nabla \mathbf u^\mathsf{T}) \bigr) + \nabla p &= \mathbf f \\ -\nabla \cdot \mathbf \mathbf u &= f \end{aligned} \right\} \text{ in $\Omega_i$,} \end{equation}
subject to the interfacial jump conditions,
\begin{equation} \label{eq:govern2} \left. \begin{aligned} \jump{\mathbf u} &= {\mathbf g}_{ij} \\ \jump{\mu (\nabla \mathbf u + \gamma \nabla \mathbf u^\mathsf{T}) \mathbf n - p \mathbf n} &= {\mathbf h}_{ij} \end{aligned} \right\} \text{ on $\Gamma_{ij}$,} \end{equation}
and boundary conditions,
\begin{equation} \label{eq:govern3} \begin{aligned} \mathbf u &= {\mathbf g}_{\partial} && \text{on } \Gamma_D, \\ \mu (\nabla \mathbf u + \gamma \nabla \mathbf u^\mathsf{T}) \mathbf n - p \mathbf n &= {\mathbf h}_{\partial} && \text{on } \Gamma_N, \end{aligned} \end{equation}
where $\Omega$ is a domain in $\mathbb R^d$ divided into one or more subdomains $\Omega_i$ (denoted ``phases''), $\Gamma_{ij} := \partial \Omega_i \cap \partial \Omega_j$ is the interface between phase $i$ and $j$, and $\Gamma_D$ and $\Gamma_N$ denote the parts of $\partial \Omega$ on which velocity Dirichlet or stress boundary conditions are imposed, respectively. Here, either $\gamma = 0$ or $\gamma = 1$ depending on whether the Stokes equations are in standard form or viscous-stress form, respectively. The operator $\jump{\cdot}$ denotes the jump in a quantity across an interface and $\mathbf n$ is to be understood from context---on $\partial \Omega$, $\mathbf n$ denotes the outward unit normal to the domain boundary, whereas for an interface $\Gamma_{ij}$, $\mathbf n$ denotes the unit normal to $\Gamma_{ij}$, oriented consistently with the definition of the jump operator. Finally, $\mu_i$ is a phase-dependent viscosity coefficient, while $\mathbf f$, $f$, $\mathbf g$, and $\mathbf h$ provide the data to the multi-phase Stokes problem and are given functions defined on $\Omega$, its boundary, and internal interfaces.

Here, we mainly consider meshes arising from Cartesian grids along with semi-unstructured quadtree/octree-based implicitly defined meshes of more complex curved domains. In this setting, it is natural to adopt a tensor-product piecewise polynomial space. Let $\mathcal E = \bigcup_i E_i$ denote the set of mesh elements, let $p \geq 1$ be an integer,\footnote{The meaning of $p$, whether as pressure or polynomial degree, should be clear from context.} and define ${\mathcal Q}_p(E)$ as the space of tensor-product polynomials of (one-dimensional) degree $p$ on element $E$. For example, ${\mathcal Q}_2$ is the space of biquadratic or triquadratic polynomials, with dimension $9$ or $27$ in 2D or 3D, respectively. Define the corresponding space of discontinuous piecewise polynomial functions as
\[ V_h = \bigl\{ u : \Omega \to \mathbb R \enskip \bigl| \enskip u|_E \in \mathcal Q_p(E) \text{ for every } E \in \mathcal E \bigr\},\]
with analogous definitions for the space of piecewise polynomial vector-valued fields, $V_h^d$, and the space of matrix-valued fields, $V_h^{d \times d}$. As discussed in ref.~\cite{CockburnKanschatSchotzauSchwab2002}, it is possible to build LDG methods for the Stokes equations wherein the discrete pressure field has either the same polynomial degree as the discrete velocity field or is in a lower degree space. In this work, we focus on the case the two have the same degree, i.e., we seek a discrete solution such that $\mathbf u_h \in V_h^d$ and $p_h \in V_h$.

In one possible construction of the LDG framework, the governing set of equations \eqref{eq:govern1}--\eqref{eq:govern3} can be discretised in a three-step process: (i) define a discrete stress tensor $\bm{\tau}_h \in V_h^{d \times d}$ equal to the discretisation of $\nabla \mathbf u_h + \gamma \nabla \mathbf u_h^\mathsf{T}$, taking into account Dirichlet source data $\mathbf g$; (ii) define $\bm{\sigma}_h \in V_h^{d \times d}$ as the viscous stress $\mu \bm{\tau}_h - p_h {\mathbb I}$ via an $L^2$ projection of $\mu \bm{\tau}_h$ onto $V_h^{d \times d}$; (iii) compute a discrete divergence of $\bm{\sigma}_h$, taking into account Neumann-like data $\mathbf h$, and add penalty stabilisation parameters for both velocity and pressure, setting the result equal to the $L^2$ projection of the given right hand side $\mathbf f$. Details of this construction, along with the associated treatment of the divergence constraint, are provided in \iffalse the Supplementary Material \fi \cref{sec:ldg}; here, we summarise the main outcomes of essential relevance. \iffalse to the remainder of this the main article.\fi The LDG discretisation results in a symmetric linear system for $(\mathbf u_h,p_h)$ having the form
\begin{equation} \label{eq:blockform} \begin{pmatrix} A & M {\mathcal G} \\ {\mathcal G}^\mathsf{T} M & -E \end{pmatrix} \begin{pmatrix} \mathbf u_h \\ p_h \end{pmatrix} = \begin{pmatrix} {\mathbf b}_{\mathbf u} \\ b_p \end{pmatrix} \end{equation}
where $M$ is the block-diagonal mass matrix, $\mathcal G$ is a discrete gradient operator, and $({\mathbf b}_{\mathbf u}, b_p)$ collects the entire influence of the source data $\mathbf f$, $f$, $\mathbf g$, and $\mathbf h$ onto the right hand side. Here, $A$ implements the viscous part of the Stokes momentum equations, and can be written in $d \times d$ block form corresponding to its action on the $d$ components of $\mathbf u_h$, with the $(i,j)$th block given by
\[ A_{ij} = \delta_{ij} \bigl( \textstyle{\sum_{k=1}^d} G_k^\mathsf{T} M_\mu G_k \bigr) + \gamma\, G_j^\mathsf{T} M_\mu G_i + \delta_{ij} \tilde{E}, \]
where $M_\mu$ is a $\mu$-weighted mass matrix, $G$ is a second discrete gradient operator closely related to the adjoint of ${\mathcal G}$, and $\tilde E$ is the operator associated with velocity penalty stabilisation. Note that if $\gamma = 0$, then $A$ is block diagonal with identical blocks corresponding to a discretisation of the Laplacian operator $-\nabla \cdot (\mu \nabla)$. Meanwhile, noting that the adjoint of $\mathcal G$ is given by $M^{-1} {\mathcal G}^\mathsf{T} M$, one may observe that the divergence constraint of the Stokes equations is implemented in the $(p,\mathbf u)$ block of \eqref{eq:blockform} through an effective discrete divergence operator which is equal to the negative adjoint of $\mathcal G$.

There is one last operator to define in \eqref{eq:blockform}, whose presence is of key importance to multigrid efficiency: the pressure stabilisation operator $E$, which weakly enforces continuity of the pressure field. The symmetric positive semidefinite matrix $E$ is defined such that\footnote{In a convenient abuse of notation, a piecewise polynomial function (e.g., in $V_h$) may carry the same notation as its corresponding coefficient vector in the basis of $V_h$, with the precise meaning understood from context. For example, in the identity $\smash{u^\mathsf{T} M v = \scalebox{0.8}{$\int_\Omega$} u\,v}$, the left hand side employs vectors and matrices relative to the chosen basis of $V_h$, whereas the right hand side employs the functional form.}
\begin{equation} \label{eq:eop} u^\mathsf{T} E v = \int_{\Gamma_0} \tau_p \,\jump{u} \jump{v} \end{equation}
holds for every $u, v \in V_h$. Here, the integral is taken over the union of every non-interfacial interior mesh face, $\jump{\cdot}$ denotes the jump across the face, and $\tau_p$ is a \textit{pressure stabilisation penalty parameter} which scales proportional to the element size $h$ and inversely proportional to the (local) viscosity coefficient:
\begin{equation} \label{eq:taueqn} \tau_p = \tau\,h / \mu, \end{equation}
where $\tau$ is a user-defined constant prefactor. Provided $\tau$ is positive, Cockburn \textit{et al} \cite{CockburnKanschatSchotzauSchwab2002} (see also extensions \cite{CockburnKanschatSchotzau2003,CockburnKanschatSchotzau2004}) prove the well-posedness of the single-phase symmetric saddle point problem in \eqref{eq:blockform}, including satisfaction of the inf-sup conditions.
In these cited works, however, the particular choice of $\tau$ in \eqref{eq:taueqn} is not extensively discussed. One of the main results in this work is to demonstrate that $\tau$ can be chosen so as to achieve excellent multigrid solver efficiency when computing solutions to the Stokes equations.

In the remainder of this paper, we will refer to the Stokes problem both in its operator block form \eqref{eq:blockform} and through the more succinct notation
\[ {\mathcal A}_h x_h = b_h, \]
where ${\mathcal A}_h$ is the symmetric saddle point operator and $x_h$ collects $\mathbf u_h$ and $p_h$ into one set of unknowns.

\section{Multigrid Design}

Prior work on designing geometric multigrid methods for LDG discretisations of \iffalse elliptic \fi Poisson-like equations \cite{ImplicitMeshPartOne,ImplicitMeshPartTwo,dgmg,fluxx} shows that one can build an efficient solver through standard multigrid concepts: a V-cycle applied to a mesh hierarchy using straightforward interpolation and restriction operators, together with standard relaxation methods, such as block Gauss-Seidel in which each block corresponds to the collective set of unknowns on each mesh element. Here, we show the same can be done for the discretised Stokes problem \eqref{eq:blockform}. (In the following, it is assumed the reader is familiar with the general design of multigrid methods; see, e.g., the books \cite{Briggs_00_01,BrandtLivne,Wesseling2001,OlshanskiiTyrtyshnikov} for reviews and applications.)

The multigrid methods designed here may be considered as a ``purely-geometric'' approach, wherein the Stokes problem is discretised on each level of the mesh hierarchy. (A convenient strategy for constructing the coarse-mesh problems---without having to explicitly form the coarse meshes themselves---is discussed shortly.) Three preliminary ingredients are needed to specify its design:
\begin{itemize}
\item \textit{Mesh hierarchy.} In this work, quadtrees and octrees are used to define the finest mesh. The tree structure naturally defines a hierarchical procedure by which to agglomerate elements to create a nested mesh hierarchy,  coarsening by a factor of two in each dimension down each level. Regarding the multi-phase case, element agglomeration is permitted only between elements of the same phase---as such, the interface is sharply preserved throughout the entire multigrid hierarchy.
\item \textit{Interpolation operator.} Owing to the presence of a nested mesh hierarchy, the interpolation operator $I_{2h}^h$, which transfers coarse mesh corrections to a fine mesh, is naturally defined via injection. In particular, we define $(I_{2h}^h u)|_{E_f} = u|_{E_c}$, where $E_f$ is a fine mesh element and $E_c \supseteq E_f$ is its corresponding coarse mesh element.
\item \textit{Restriction operator.} The restriction operator $R_h^{2h}$, which transfers the residual of a fine-mesh problem to the coarse mesh, is defined as the $L^2$ projection onto the coarse mesh (or, equivalently, as the adjoint of the interpolation operator). It is related to the interpolation operator via $R_h^{2h} = M_{2h}^{-1} (I_{2h}^h)^\mathsf{T} M_h$, where $M_h$ and $M_{2h}$ are the mass matrices of the two meshes.
\end{itemize}
The last essential multigrid ingredient, and perhaps most important, is the relaxation method. As mentioned earlier, we have used a simple block Gauss-Seidel method, where each block corresponds to the collective set of degrees of freedom (i.e., velocity and pressure combined) on each mesh element. Specifically, consider a repartitioning of $\mathcal A x = b$ according to these blocks, such that $x_i$ denotes the set of velocity and pressure values on element $i$, and ${\mathcal A}_{ij}$ denotes the $(i,j)$th block of ${\mathcal A}$, whence $b_i = \sum_j {\mathcal A}_{ij} x_j$. Then, the block Gauss-Seidel method simply sweeps over the elements, in some particular order, replacing $x_i \leftarrow {\mathcal A}_{ii}^{-1}(b_i - \sum_{j \neq i} {\mathcal A}_{ij} x_j)$. Here, ${\mathcal A}_{ii}$ is the $i$th diagonal block of ${\mathcal A}$ and takes on the form of a miniature Stokes operator; referring to \eqref{eq:blockform}, we have
\[ {\mathcal A}_{ii} = \begin{pmatrix} A_{ii} & M_{ii} {\mathcal G}_{ii} \\ {\mathcal G}_{ii}^\mathsf{T} M_{ii} & -E_{ii} \end{pmatrix}. \]
Note that ${\mathcal A}_{ii}$ needs to be inverted in the Gauss-Seidel update of element $i$. Assuming that the global Stokes saddle-point problem ${\mathcal A} x = b$ satisfies the inf-sup conditions, it is straightforward to show that so too does ${\mathcal A}_{ii}$, and hence the local element-wise problem is well-posed; this has also been confirmed through numerous and extensive numerical tests. In our specific implementation, we precompute a symmetric indefinite factorisation of ${\mathcal A}_{ii}$ for every $i$, and use this factorisation as a direct solver for each of these mini-Stokes problems in the Gauss-Seidel sweep.
Regarding the element ordering, we have opted for a multi-coloured Gauss-Seidel method. The primary reason for this choice is that a multi-coloured sweep affords a simpler parallel implementation of the method, both in terms of multi-threading and in a distributed environment (e.g., through standard domain decomposition methods using MPI).

\begin{algorithm}[t]
	\caption{\sffamily\small Multigrid V-cycle $V({\mathcal E}_h, x_h, b_h)$ with $\nu_1$ pre- and $\nu_2$ post-smoothing steps on mesh ${\mathcal E}_h$ of the hierarchy}
	\begin{algorithmic}[1]
		\If{${\mathcal E}_h$ is the bottom level}
			\State Solve ${\mathcal A}_h x_h = b_h$ with bottom solver
		\Else
			\State Apply block Gauss-Seidel relaxation $\nu_1$ times
			\State $r_{2h} := (I_{2h}^h)^\mathsf{T} (b_h - {\mathcal A}_h x_h)$
			\State $x_{2h} := V({\mathcal E}_{2h}, 0, r_{2h})$
			\State $x_h \leftarrow x_h + I_{2h}^h x_{2h}$
			\State Apply block Gauss-Seidel relaxation $\nu_2$ times
		\EndIf
		\State{\textbf{return} $x_h$}
	\end{algorithmic}%
	\label{algo:vcycle}%
\end{algorithm}%

Using the defined interpolation and restriction operators and the block Gauss-Seidel relaxation method, the construction of a multigrid V-cycle is relatively standard and is outlined in \cref{algo:vcycle}.\footnote{Note that $(I_{2h}^h)^\mathsf{T}$ appears on line 5, rather than the restriction operator $R_h^{2h}$; this follows from a convenient simplification common to many finite element methods: briefly, viewed as an operator which maps $\smash{V_h^d \otimes V_h}$ to $\smash{V_h^d \otimes V_h}$, the discrete Stokes operator is given by $\smash{M^{-1} {\mathcal A}}$. Therefore, the residual of the fine mesh problem, as a piecewise polynomial function, is $\smash{M_h^{-1} (b_h - {\mathcal A}_h x_h)}$. This residual is then multiplied by $\smash{R_h^{2h}}$ to define the source data for the coarse mesh problem $\smash{M_{2h}^{-1} {\mathcal A}_{2h} = R_h^{2h} M_h^{-1} (b_h - A_h x_h)}$. Rearranging, one obtains line 5.} In this algorithm, $\mathcal A_h$ is assumed to be pre-computed on every level of the mesh; a particularly convenient method for doing so---without having to explicitly mesh each level, build quadrature rules for coarse mesh elements, or build LDG operators via coarse-mesh numerical fluxes, etc.---uses the operator-coarsening ideas of ref.~\cite{dgmg}. In this technique, the discrete gradient and penalty operators underlying \eqref{eq:blockform} are coarsened solely based on the interpolation operator hierarchy, using simple block-sparse linear algebra; these methods are further described in \iffalse the Supplementary Material \fi \cref{sec:op}. Returning to \cref{algo:vcycle}, note that the V-cycle computes coarse-grid corrections for both the velocity and pressure, simultaneously, and there is no need to strictly enforce the divergence constraint on any level of the hierarchy. Regarding the bottom solver, in this work a direct solver using a symmetric indefinite factorisation of ${\mathcal A}_h$ on the coarsest mesh is used, together with an appropriate treatment of its associated trivial kernel.\footnote{With periodic boundary conditions, the Stokes problem has a trivial kernel of dimension $d+1$, spanned by constant velocity and pressure fields; with velocity Dirichlet boundary conditions, the kernel is one-dimensional, spanned by constant pressure fields; with stress boundary conditions in viscous-stress form, the kernel is spanned by constant velocity fields as well as less trivial velocity modes such as, e.g., the velocity field $(x,y) \mapsto (-y,x)$ in 2D. The bottom solver robustly treats these modes through a simple least squares approach which (pre)computes the symmetric eigendecomposition of $\mathcal A_h$, ``snapping'' any nearly-zero eigenvalues to exactly zero.}

As is typical, applying more and more pre- and post-smoothing steps increases the convergence rate of the multigrid solver, but at greater computational cost. According to a variety of numerical experiments, a general observation made in this work is that a V-cycle with three pre- and post-smoothing steps is a good all-rounder, based on the metric of fastest computation time in reducing solution error by a given factor. On occasion, four pre-smoothing steps and two post-smoothing steps, or vice versa, performs marginally better, but on a problem-specific basis. Naturally, the optimal choice of multigrid design parameters is implementation- and problem-dependent, influenced by a wide variety of aspects, e.g., the relative computational costs of interpolation, restriction, and relaxation operators, or computing hardware characteristics, such as shared memory or distributed memory architectures and their associated memory communication costs. Further comments on V-cycle design, or counterparts such as W-cycles, are provided in the concluding remarks.

To complete the description of the multigrid method, we note that although the V-cycle can be used as a standalone iterative solver, solver efficiency can be further accelerated by using it as a preconditioner of a Krylov method \cite{BenziGolubLiesen2005}. In this work, we have used  a single V-cycle (with $\nu_1=\nu_2=3$ pre- and post-smoothing steps) as a left-preconditioner of the GMRES method. Specifically, applying the V-cycle to an initial guess of zero on the fine mesh results in a linear operator, denoted in the remainder of this paper as $V$; the preconditioned system is then $V \mathcal A$. For simplicity, we do not consider restarted variants of GMRES here, in part because experiments show that convergence is generally attained in as few as 5-15 steps for a ten-fold reduction in the order of magnitude of the residual. Moreover, convergence behaviour is generally smooth during the iterations, such that the residual reduces in norm by a nearly constant factor each iteration of the GMRES method.

\section{Influence of Pressure Penalty Stabilisation on Multigrid Efficiency}
\label{sec:tauinfluence}

As outlined, the Stokes multigrid solver consists of a standard V-cycle, using a block Gauss-Seidel relaxation method in which each block corresponds to the collective set of degrees of freedom, of both velocity and pressure, on each mesh element. Key to rapid multigrid convergence is an apt choice of the user-defined pressure penalty stabilisation prefactor parameter $\tau$ in \eqref{eq:taueqn}. In general terms, if $\tau$ is below some positive threshold, the V-cycle fails to converge; above this threshold, there is a range of values for which convergence rates can match that of fast geometric multigrid methods for scalar Poisson problems; and, beyond this range, multigrid efficiency will degrade.

\begin{figure
\centering
\includegraphics[scale=0.91]{tau-spectral}
\caption{Spectral properties of the multigrid-preconditioned system $V{\mathcal A}$ as a function of pressure penalty stabilisation parameter; note the quasi-logarithmic axis. Results correspond to a two-dimensional Stokes problem in standard form with $p = 2$; similar characteristics are obtained with other $p$, or in 3D, or with the viscous-stress form of the Stokes equations.}
\label{fig:tauspectral}
\end{figure}

To illustrate this behaviour, as well as the suitability of the V-cycle as a preconditioner for the Stokes system ${\mathcal A}$, we examine the spectral properties of the preconditioned system $V {\mathcal A}$. In this particular example, we consider the standard form of the Stokes equations ($\gamma = 0$), with periodic boundary conditions on a $128 \times 128$ Cartesian grid using $p = 2$ biquadratic elements. As a function of $\tau$, \cref{fig:tauspectral} plots three quantities concerning the spectrum\footnote{The spectral analysis deliberately excludes the zero eigenvalues associated with the trivial kernel of the Stokes operator. In addition, regarding the results plotted in \cref{fig:tauspectral}, the spectral extremes have been estimated via the GMRES method, through computation of the eigenvalues of the upper Hessenberg matrix of the corresponding Arnoldi iteration. Although only an approximation to the true spectrum of $V{\mathcal A}$, estimation via GMRES is significantly more efficient than forming and computing the spectrum of the dense matrix $V{\mathcal A}$; furthermore, the computed quantities are found in practice to be sufficiently accurate for the present purpose.} of $V{\mathcal A}$: (i) the real part of the right-most eigenvalue, $\max_i \textsf{Re}\,\lambda_i$; (ii) the real part of the left-most eigenvalue, $\min_i \textsf{Re}\,\lambda_i$; and (iii) the greatest imaginary part, $\max_i \textsf{Im}\,\lambda_i$. Ideally, the eigenvalues of $V{\mathcal A}$ should be clustered around $1$, and we observe this is the case when $\tau \approx 0.1$; furthermore, near this value, the eigenvalues are nearly real. However, if $\tau$ is too small, then the left-most eigenvalue of $V{\mathcal A}$ crosses the imaginary axis; in \cref{fig:tauspectral} this occurs when $\tau \lessapprox 10^{-2}$. This represents a breakdown of the V-cycle, as then the (non-trivial) eigenvalues of $V{\mathcal A}$ cease to be bounded away from zero. (In fact, numerical experiments examining the efficacy of the V-cycle as a standalone iterative method show that the V-cycle ceases to have spectral radius less than one at this same point.) Meanwhile, if $\tau$ is too large, e.g., $\tau \gtrapprox 1$, then the eigenvalues of $V{\mathcal A}$ begin to diverge away from $1$.

\begin{figure
\centering
\includegraphics[scale=0.91]{tau-rho}
\caption{Multigrid efficiency as a function of pressure penalty stabilisation parameter. Here, the convergence rate $\rho$, defined by eq.~\ref{eq:rho}, measures the average reduction factor per iteration in the residual of the V-cycle preconditioned GMRES method.}
\label{fig:taurho}
\end{figure}

In this work, the primary metric used to assess multigrid efficiency is the convergence rate $\rho$ of the multigrid-preconditioned GMRES method. This metric correlates with the spectral characteristics described above, and also to the performance of the V-cycle as a standalone iterative solver. Here, $\rho$ is defined as the average residual reduction factor per iteration of the (left) preconditioned GMRES method,
\begin{equation} \label{eq:rho} \rho = \exp \Bigl( \frac1n \log \frac{\|V{\mathcal A} x_n - V b\|_2}{\|V{\mathcal A} x_0 - V b\|_2} \Bigr), \end{equation}
where $n$ is the number of iterations required to reduce the residual by a factor of $10^8$ from its starting value. In particular, a right hand side of $b = 0$ is used, with initial guess $x_0$ given by a randomly generated $(d+1)$-dimensional vector field. With high probability, the randomly generated field contains modes which are damped slowest by the multigrid method, and thus $\rho$ represents a typical ``worst case'' convergence rate. \cref{fig:taurho} shows the convergence rate $\rho$ as a function of $\tau$ for the same example considered in \cref{fig:tauspectral}. Optimal convergence is attained when $\tau \approx 0.1$, precisely when the spectrum of $V{\mathcal A}$ is tightly clustered around 1; meanwhile, when $\tau \lessapprox 10^{-2}$, $\rho$ is approximately one, representing the fact that GMRES is unable to effectively reduce the residual owing to a breakdown of the preconditioner.

Similar behaviour to that seen in \cref{fig:tauspectral} and \cref{fig:taurho} is observed for other choices of polynomial degree $p$, for different grid sizes, for the viscous-stress form of the Stokes equations, and in 3D as well as 2D. In all cases, the convergence rate $\rho$ exhibits a well-defined valley as a function of $\tau$. At present, a formula for the corresponding optimal value of $\tau$ is not known. In this work, a simple one-dimensional parameter sweep was used to find the optimal value of $\tau$ on successive grid sizes $n \times n \,(\,\times n)$ for $n = 4,8,16,\ldots$ up to $n = 256$ in 2D and $n = 128$ in 3D; experiments indicated that the optimum essentially converges around $n = 64$ or $128$, beyond which $\argmin_{\tau} \rho$ is relatively insensitive to the grid size. \cref{tab:optimaltau} contains the results of the search for a variety of $p$ in 2D and 3D. In addition to the optimal value, a ``window'' of acceptable $\tau$ may also be computed: one can search for all $\tau$ such that $\rho(\tau) \leq (\min \rho)^{1 - \epsilon}$, where $0 \leq \epsilon < 1$ represents a user-defined threshold for which the number of multigrid iterations increases by a factor of $1/(1-\epsilon)$ above the optimal minimum.
For example, with $\epsilon = 1/9$, the corresponding range of $\tau$ values yields at most 12.5\% more iterations than optimal; \iffalse Table SM1 in the Supplementary Material \fi \cref{tab:tauwindow} contains the corresponding ranges, and shows that, even if $\tau$ is not chosen exactly at the optimum, there is a relatively wide range of values for $\tau$ that will nevertheless attain good multigrid efficiency. Finally, numerical experiments indicate that the pressure penalty parameter has very little influence on the velocity or pressure discretisation error (see, for example, the demonstration given in \cref{app:tauaccuracy}); this is ideal, as it allows us to concentrate mainly on the impact of $\tau$ on multigrid performance. Further comments concerning the selection of $\tau$ are given in the concluding remarks.

\begin{table
\centering%
\sffamily\footnotesize%
\caption{Optimal values of pressure penalty stabilisation parameter, attaining minimal multigrid iteration count.}
\begin{tabular}{lcccccc}
\multirow{2}{*}{Stokes form} & \multirow{2}{*}{$d$} & \multicolumn{5}{c}{Polynomial degree $p$} \\
 & & 1 & 2 & 3 & 4 & 5 \\
\midrule
\multirow{2}{*}{Standard} & 2D & 0.19 & 0.10 & 0.086 & 0.019 & 0.031 \\
& 3D & 0.12 & 0.088 & 0.084 & -- & -- \\
\midrule
\multirow{2}{*}{Viscous stress} & 2D & 0.14 & 0.046 & 0.034 & 0.0095 & 0.011 \\
& 3D & 0.12 & 0.039 & 0.040 & -- & -- \\
\bottomrule
\end{tabular}
\label{tab:optimaltau}
\end{table}

\begin{table
\centering%
\caption{Range of pressure penalty stabilisation parameters for which the number of multigrid iterations is at most 12.5\% more than optimal.}%
\sffamily%
\footnotesize%
\newcommand{1.5in}{0.83em}%
\begin{tabular}{l@{\hspace{1.5in}}c@{\hspace{1.5in}}c@{\hspace{1.5in}}c@{\hspace{1.5in}}c@{\hspace{1.5in}}c@{\hspace{1.5in}}c}
\multirow{2}{*}{Stokes form} & \multirow{2}{*}{$d$} & \multicolumn{5}{c}{Polynomial degree $p$} \\
 & & 1 & 2 & 3 & 4 & 5 \\
\midrule
\multirow{2}{*}{Standard} & 2D & {(0.15, 0.28)} & {(0.082, 0.12)} & {(0.067, 0.11)} & {(0.013, 0.029)} & {(0.021, 0.041)} \\
& 3D & {(0.061, 0.30)} & {(0.070, 0.11)} & {(0.064, 0.12)} & -- & -- \\
\midrule
\multirow{2}{*}{Viscous stress} & 2D & {(0.091, 0.18)} & {(0.040, 0.056)} & {(0.027, 0.043)} & {(0.0058, 0.021)} & {(0.0072, 0.020)} \\
& 3D & {(0.025, 0.19)} & {(0.031, 0.066)} & {(0.021, 0.059)} & -- & -- \\
\bottomrule
\end{tabular}
\label{tab:tauwindow}
\end{table}

\section{Multigrid Efficiency for the Time-Independent Stokes Equations}
\label{sec:timeindep}

In the next two sections, multigrid performance is examined for a variety of Stokes problems, in both standard and viscous-stress form, and with different types of boundary conditions. We also consider test cases with variable viscosity, multi-phase problems exhibiting large discontinuities in $\mu$ across an embedded interface, and curved geometry problems which use semi-unstructured meshes. Our primary focus is on demonstrating effective multigrid performance under the action of mesh refinement. The order of accuracy in the velocity and pressure, in the $L^2$ and maximum error norms, is also measured and reported. In two dimensions, grid sizes typically range from $4 \times 4$ up to $1024 \times 1024$, with polynomial degrees $p = 1, 2, \ldots, 5$; in three dimensions, owing to limited computing resources, only $p = 1, 2,$ and $3$ is considered (i.e., trilinear, triquadratic, and tricubic polynomials) on grid sizes up to $128 \times 128 \times 128$ (the largest of these problems has half a billion degrees of freedom and requires 1TB of memory to store just the block-diagonal component of the block-sparse matrix ${\mathcal A}$). In the remainder of this paper, for every test case, the pressure penalty stabilisation parameter is chosen equal to the values reported in \cref{tab:optimaltau}. Multigrid convergence rates are assessed using the average reduction factor in residual per iteration of the GMRES method, i.e., using \eqref{eq:rho}, on a test problem with right hand side $b = 0$, initial guess defined by a randomly generated $(d+1)$-dimensional vector field, over as many iterations as necessary to reduce the initial residual by a factor of $10^8$.

\subsection{Periodic boundary conditions}
\label{prob:periodic}

\begin{figure
\centering%
\includegraphics[scale=0.91]{periodic-stokes-x}%
\caption{Measured multigrid convergence rates when solving the Stokes equations in standard form, with $\mu = 1$ and periodic boundary conditions.}%
\label{fig:periodic-stokes}%
\end{figure}

We begin with perhaps the simplest Stokes problem, i.e., the Stokes equations in standard form, with $\mu = 1$ and periodic boundary conditions, on the unit square/cube domain $\Omega = (0,1)^d$. \Cref{fig:periodic-stokes} plots the measured convergence rate $\rho$ as a function of grid size $n \times n\,(\times\, n)$. With the exception of $p = 1$, ideal convergence rates are attained. For example, $\rho \approx 0.05$ corresponds to needing only seven iterations to reduce the residual by a factor of $10^9$. When $p = 1$, however, slower convergence is seen---compared to the higher-degree cases, it appears that the an element-wise block Gauss-Seidel relaxation method is less effective for a bilinear and trilinear LDG discretisation of the Stokes equations. This behaviour is consistently observed across all of the presented tests; see also the next two examples. Regarding the order of accuracy, numerical experiments show that both velocity and pressure attain order $p + 1$ accuracy, in both the $L^2$ and maximum error norms, for all $p$ considered, in 2D and 3D; see \iffalse Fig.~SM3 in the Supplementary Material. \fi \cref{sec:gridconv}.

\subsection{Dirichlet boundary conditions}
\label{prob:dirichlet}

\begin{figure
\centering%
\includegraphics[scale=0.91]{dirichlet-stokes-x}%
\caption{Measured multigrid convergence rates when solving the Stokes equations in standard form, with $\mu = 1$ and velocity Dirichlet boundary conditions.}%
\label{fig:dirichlet-stokes}%
\end{figure}

\begin{figure}
\centering
\newcommand{1.5in}{1.5in}
\begin{tabular}{ccc}
{\setlength{\fboxsep}{0pt}\setlength{\fboxrule}{0.5pt}\fbox{\includegraphics[width=1.5in]{uerr}}}
 & {\setlength{\fboxsep}{0pt}\setlength{\fboxrule}{0.5pt}\fbox{\includegraphics[width=1.5in]{verr}}}
 & {\setlength{\fboxsep}{0pt}\setlength{\fboxrule}{0.5pt}\fbox{\includegraphics[width=1.5in]{perr}}} \\[1em]
\makebox[2cm][c]{{\begingroup%
\setlength{\unitlength}{2cm}%
\def\ticwidth{0.03}%
\def\ticsep{-0.02}%
\setlength\fboxsep{0pt}%
\begin{picture}(1,0.1)%
	\put(0,\ticwidth){\fbox{\includegraphics[width=0.995\unitlength,height=0.05\unitlength]{bluegreen}}}%
	\put(0.0035,0){\line(0,1){\ticwidth}}%
	\put(0.25,0){\line(0,1){\ticwidth}}%
	\put(0.5,0){\line(0,1){\ticwidth}}%
	\put(0.75,0){\line(0,1){\ticwidth}}%
	\put(1.0055,0){\line(0,1){\ticwidth}}%
	\put(0.0,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$-0.000015\hphantom{-}$}}}}%
	\put(0.5,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0$}}}}%
	\put(1,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0.000015$}}}}%
	\put(0.5,0.23){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.75}{$u-u_h$}}}}%
\end{picture}%
\global\let\ticwidth\undefined%
\global\let\ticsep\undefined%
\endgroup}}
&
\makebox[2cm][c]{{\begingroup%
\setlength{\unitlength}{2cm}%
\def\ticwidth{0.03}%
\def\ticsep{-0.02}%
\setlength\fboxsep{0pt}%
\begin{picture}(1,0.1)%
	\put(0,\ticwidth){\fbox{\includegraphics[width=0.995\unitlength,height=0.05\unitlength]{bluegreen}}}%
	\put(0.0035,0){\line(0,1){\ticwidth}}%
	\put(0.25,0){\line(0,1){\ticwidth}}%
	\put(0.5,0){\line(0,1){\ticwidth}}%
	\put(0.75,0){\line(0,1){\ticwidth}}%
	\put(1.0055,0){\line(0,1){\ticwidth}}%
	\put(0.0,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$-0.000015\hphantom{-}$}}}}%
	\put(0.5,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0$}}}}%
	\put(1,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0.000015$}}}}%
	\put(0.5,0.23){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.75}{$v-v_h$}}}}%
\end{picture}%
\global\let\ticwidth\undefined%
\global\let\ticsep\undefined%
\endgroup}}
&
\makebox[2cm][c]{{\begingroup%
\setlength{\unitlength}{2cm}%
\def\ticwidth{0.03}%
\def\ticsep{-0.02}%
\setlength\fboxsep{0pt}%
\begin{picture}(1,0.1)%
	\put(0,\ticwidth){\fbox{\includegraphics[width=0.995\unitlength,height=0.05\unitlength]{greenblue}}}%
	\put(0.0035,0){\line(0,1){\ticwidth}}%
	\put(0.25,0){\line(0,1){\ticwidth}}%
	\put(0.5,0){\line(0,1){\ticwidth}}%
	\put(0.75,0){\line(0,1){\ticwidth}}%
	\put(1.0055,0){\line(0,1){\ticwidth}}%
	\put(0.0,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$-0.004\hphantom{-}$}}}}%
	\put(0.5,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0$}}}}%
	\put(1,\ticsep){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.6}{$0.004$}}}}%
	\put(0.5,0.23){\makebox(0,0)[ct]{\rotatebox{0}{\scalebox{0.75}{$p-p_h$}}}}%
\end{picture}%
\global\let\ticwidth\undefined%
\global\let\ticsep\undefined%
\endgroup}}
\end{tabular}
\caption{Illustration of the discrete error for the test case considered in \S\ref{prob:dirichlet}, corresponding to a single-phase Stokes problem in standard form with velocity Dirichlet boundary conditions. The error in velocity $\mathbf u = (u,v)$ and pressure $p$ is shown in the case of a $16 \times 16$ Cartesian mesh, for $p = 3$ bicubic polynomials. Note the numerical boundary layer in pressure, which according to grid convergence analyses, does not impact the maximum norm optimal order accuracy of the velocity field.}
\label{fig:2dsquare}
\end{figure}

The next test problem is identical to the previous, but with velocity Dirichlet boundary conditions imposed on $\partial \Omega$. \Cref{fig:dirichlet-stokes} plots the measured multigrid convergence rates and shows that, with the exception of $p = 1$, excellent multigrid performance is attained, similar to that of the periodic case in \cref{fig:periodic-stokes}. When $p = 1$, we observe a stronger failure of multigrid efficiency, with $\rho$ diverging toward one as the mesh is refined. Another difference compared to the periodic case concerns the order of accuracy: numerical experiments show that the velocity attains order $p + 1$, in both the $L^2$ and maximum error norms, whereas the pressure field attains order $p+\tfrac12$ in the $L^2$ norm, and order $p$ in the maximum norm; see \iffalse Fig.~SM4 in the Supplementary Material.\fi \cref{sec:gridconv}. The order reduction in the computed pressure field, as compared to the order $p + 1$ observed in the periodic case, is due to a numerical boundary layer \iffalse (see, e.g., Fig.~SM2 in the SM); \fi (see, e.g., \cref{fig:2dsquare}); it is important to note, however, that this numerical boundary layer does not impact the optimal order accuracy of the computed velocity. See also the discussion of Cockburn \textit{et al} \cite{CockburnKanschatSchotzauSchwab2002}, wherein a priori estimates shows that order $p$ accuracy in pressure is to be expected for this LDG discretisation of the Stokes equations.

\subsection{Stress boundary conditions}
\label{prob:stress}

\begin{figure
\centering%
\includegraphics[scale=0.91]{neumann-stress-x}%
\caption{Measured multigrid convergence rates when solving the Stokes equations in viscous-stress form, with $\mu = 1$ and stress boundary conditions.}%
\label{fig:neumann-stress}%
\end{figure}

We next test performance of the Stokes multigrid solver when stress boundary conditions are imposed, in which case the pertinent form of the Stokes equations is the viscous-stress form ($\gamma = 1$). Similar to the previous test problems, eq.~\eqref{eq:stokes2} with $\mu = 1$ is solved on a $n \times n\,(\times\,n)$ mesh of a unit square/cube domain $\Omega = (0,1)^d$, with stress boundary conditions ${\bm \sigma} \cdot \mathbf n = {\mathbf h}_{\partial}$ imposed on $\partial \Omega$. \Cref{fig:neumann-stress} plots the measured multigrid convergence rates. Compared to the previous two test problems (which employed the standard form of the Stokes equations), we observe a mild increase in $\rho$ for the cases in which $p > 1$, most visible in 2D. The slight increase in $\rho$ is mainly attributed to the viscous-stress form of the Stokes equations, and not solely to the imposition of stress boundary conditions. (Indeed, if one imposes Neumann-like boundary conditions for the Stokes equations in standard form, convergence rates similar to those in \cref{fig:periodic-stokes} and \cref{fig:dirichlet-stokes} are obtained.) In \cref{fig:neumann-stress}, for $p = 1$, we once again see less than ideal multigrid efficiency, though with marginal improvements compared to the case of velocity Dirichlet boundary conditions. Since the case of $p = 1$ is generally not of significant practical interest in the context of high-order accurate DG methods, we will focus on degrees $p > 1$ in the remainder of the presented results. Meanwhile, for the current test problem, numerical experiments show that the velocity field attains order $p+1$ accuracy in the maximum norm, while the pressure field attains order $p + \tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm, see \cref{sec:gridconv}; as a general rule, our results indicate that whenever boundary conditions (or interfacial jump conditions) are imposed, the pressure field loses one order of accuracy near the boundary (or interface) owing to a numerical boundary layer, but this never affects the optimal order accuracy of the computed velocity field.

\subsection{Variable viscosity}
\label{prob:variableviscosity}

\begin{figure
\centering%
\includegraphics[scale=0.91]{variable-mu-stress-x}%
\caption{Measured multigrid convergence rates when solving the Stokes equations in viscous-stress form with a non-constant viscosity function $\mu: \Omega \to {\mathbb R}^+$ given by \eqref{eq:mu}, and stress boundary conditions.}%
\label{fig:variable-mu-stress}%
\end{figure}

In the next example, we consider the possibility of a non-constant viscosity function $\mu : \Omega \to {\mathbb R}^+$, varying throughout the domain. This problem serves three main purposes: (i) to test the application of the local inverse scaling by $\mu$ of the pressure penalty stabilisation parameter suggested in \eqref{eq:taueqn}; (ii) to demonstrate whether or not multigrid efficiency is impacted by variable ellipticity coefficient; and (iii) to examine the order of accuracy of the discrete solution in the variable coefficient case. Specifically, we consider the Stokes equations in viscous-stress form on the unit square/cube domain $\Omega = (0,1)^d$, where $\mu : \Omega \to {\mathbb R}^+$ is given by
\begin{equation} \label{eq:mu} \mu = \begin{cases} 1 + \tfrac12 \sin 4\pi x \sin 4 \pi y & \text{in 2D,} \\ 1 + \tfrac12 \sin 4\pi x \sin 4\pi y \sin 4 \pi z & \text{in 3D,} \end{cases} \end{equation}
together with stress boundary conditions on $\partial\Omega$. Measured multigrid convergence rates, plotted in \cref{fig:variable-mu-stress}, show very similar behaviour to the constant-viscosity test problem of the previous example (\cref{fig:neumann-stress}). Meanwhile, numerical experiments examining the order of accuracy (see \cref{sec:gridconv}) confirm that the velocity field attains order $p+1$ in the maximum norm, while pressure attains order $p+\tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm.

\subsection{Curved domain geometry}
\label{prob:sphere}

\begin{figure
\centering\sffamily\footnotesize%
\begin{minipage}{0.3\textwidth}\centering%
\includegraphics[width=\textwidth]{spheremesh}\\[0.5em]
(a) Circular domain
\end{minipage}
\qquad\qquad
\begin{minipage}{0.3\textwidth}\centering%
\includegraphics[width=\textwidth]{circlemesh}\\[0.5em]
(b) Circular interface
\end{minipage}%
\caption{Examples of implicitly defined meshes generated with a background $16 \times 16$ Cartesian grid. (a) Circular domain, used in the single-phase Stokes problem examined in \cref{fig:sphere}. (b) Square domain with an embedded circular interface, used in the two-phase time-dependent Stokes problem considered later in this article, in \cref{fig:water-bubble}.}%
\label{fig:circle}%
\end{figure}

In the last single-phase example of this section, we consider a Stokes problem in a curved domain. Here, and in all others involving curved geometry, we make use of implicitly defined meshes, building on the implicit mesh DG framework developed in prior work \cite{ImplicitMeshPartOne,ImplicitMeshPartTwo}. Briefly, an implicitly defined mesh uses one or more level set functions, describing the domain geometry and/or embedded interfaces, to cut through the cells of a background quadtree or octree; tiny cut cells are then merged with neighbouring cells to create a mesh such that elements adjacent to curved geometry have their shape defined implicitly by the level set functions. In particular, the resulting mesh is interface- and boundary-conforming, thereby sharply representing the implicitly-defined geometry. To use an implicitly defined mesh with a DG method, the main task is to compute quadrature schemes for the elements and faces whose geometry is implicitly defined; these are computed using the high-order accurate algorithms detailed in refs.~\cite{HighorderImplicitQuad,algoim}, and then used when computing mass matrices, discrete gradient operators, $L^2$ projections, etc. For details on the implicit mesh DG framework, see refs.~\cite{ImplicitMeshPartOne,ImplicitMeshPartTwo}, and for illustrations demonstrating the associated (implicitly-formed) mesh hierarchy underlying the multigrid method, see ref.~\cite{dgmg}.

\begin{figure
\centering%
\includegraphics[scale=0.91]{sphere-x}%
\caption{Measured multigrid convergence rates when solving the Stokes equations in standard form, with $\mu = 1$ in a unit diameter spherical domain using implicitly defined meshes together with velocity Dirichlet boundary conditions. Here, $n$ denotes the number of cells in the background uniform Cartesian $n\times n\,(\,\times n)$ grid used to build the corresponding implicitly defined mesh; see \cref{fig:circle}(a).}%
\label{fig:sphere}%
\end{figure}

In the current test problem, we consider solving the Stokes equations, in standard form, on a unit diameter circle or sphere, $\Omega = \{ x \in {\mathbb R}^d : |x| < \tfrac12 \}$, with $\mu = 1$ and velocity Dirichlet boundary conditions. An example of the corresponding implicitly defined mesh generated by a background $16 \times 16$ Cartesian grid is shown in \cref{fig:circle}(a). Measured multigrid convergence rates for this test case are shown in \cref{fig:sphere} for polynomial degrees $p = 2$ and $p = 3$. Overall, efficient multigrid convergence rates are witnessed, however once $n \geq 256$ we observe that $\rho$ increases in value compared to previous test cases; we attribute this to the worsening conditioning of the discrete Stokes operator. In particular, inspection of numerical results shows that the conditioning of ${\mathcal A}$ exceeds $10^8$ in these cases. Since the V-cycle is built by coarsening the components of the fine mesh operator ${\mathcal A}$, and since $\rho$ is measured according to how many GMRES iterations it takes to reduce the (preconditioned) residual by a factor $10^8$, it follows that $\rho$ may be impacted by this conditioning.
Regarding the grid convergence analysis (see \cref{sec:gridconv}), numerical results once again confirm that the velocity field attains order $p+1$ in the maximum norm, whereas the pressure attains order $p+\tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm.

\subsection{Multi-phase Stokes equations with interfacial jump conditions}
\label{prob:square}

In this last example for the steady-state Stokes equations, we consider a challenging situation in which the viscosity coefficient $\mu$ exhibits jumps several orders in magnitude across an embedded interface. In particular, let $\Omega = (0,1)^d$ be divided into an interior square or cubic phase, $\Omega_1 = (\tfrac14,\tfrac34)^d$, together with an exterior phase, $\Omega_2 = \Omega \setminus \overline{\Omega_1}$. Four different ratios of viscosity jump across the interface are considered: in every case, the exterior phase has unit viscosity, $\mu_2 = 1$, whereas the interior phase will have one of four values, $\mu_1 \in \{10^{-6}, 10^{-3}, 10^3, 10^6\}$. We consider the multi-phase viscous-stress form of the Stokes equations, \eqref{eq:govern1}--\eqref{eq:govern3} with $\gamma = 1$, for which interfacial jump conditions in velocity and the stress tensor are imposed on $\Gamma = \partial \Omega_1 \cup \partial \Omega_2$.
To tackle this problem, we apply two distinct but complementary strategies: viscosity-upwinded weighted fluxes on interfacial mesh faces, and a diagonal scaling to improve the conditioning of the Stokes operator; these are discussed next.

Prior work on LDG methods for elliptic interface problems \cite{fluxx} shows that, in order to obtain ideal multigrid efficiency and solution accuracy, one should apply a biasing strategy for the LDG numerical fluxes on interfacial faces, wherein $\hat{\mathbf u}$ and $\hat{\bm{\sigma}}$ (see \iffalse the Supplementary Material \fi \cref{sec:ldg}) are biased toward one phase or the other, depending on the local viscosity coefficients. Here, we build on these ideas and extend it the multi-phase Stokes case.

The strategy of biasing can be intuitively motivated as follows. Suppose that the interior phase has a vastly smaller viscosity than the exterior phase, i.e., $\mu_1 = \epsilon$ with $0 < \epsilon \ll 1$. Suppose also that the velocity $\mathbf u$ and its gradient is unit order in magnitude near the interface. From a physical standpoint, the pressure scales as $p \sim \mu U/L$ (where $U$ is a typical velocity scale and $L$ is a typical length scale), and so, suppose also that pressure and its gradient is unit order magnitude in the exterior phase, and scales with $\epsilon$ in the interior phase. Then, in the limit of vanishingly small $\epsilon$, the stress jump condition in \eqref{eq:govern2} approximately reduces to phase $\Omega_2$ having a stress-like boundary condition, $\mu_2 (\nabla \mathbf u + \nabla \mathbf u^T) \mathbf n - p \mathbf n \approx {\mathbf h}_{12}$ on $\Gamma$. Except for the modes associated with the trivial kernel of the corresponding Stokes operator, this is enough to determine the solution $(\mathbf u_2,p_2)$ in $\Omega_2$ and thus the phase $\Omega_2$ (nearly) decouples from phase $\Omega_1$. Once $(\mathbf u_2,p_2)$ is found, the remaining jump condition in \eqref{eq:govern2} essentially reduces to a Dirichlet boundary condition for the Stokes problem in phase $\Omega_1$, i.e., $\mathbf u_1|_\Gamma \approx  {\mathbf u}_2|_\Gamma - {\mathbf g}_{12}$, which is enough to determine $(\mathbf u_1,p_1)$ up to the modes associated with its kernel. Therefore, in the limit of vanishingly small $\epsilon$, the highly viscous phase essentially ``sees'' a stress boundary condition on $\Gamma$ (whose data is nearly independent of the other phase), and the nearly inviscid phase ``sees'' a Dirichlet boundary condition on $\Gamma$ (whose data depends on the solution across the interface). It follows, therefore, that it may be advantageous to bias the numerical fluxes of the LDG formulation in the same way, to reflect the physical nature of the interfacial jump conditions in \eqref{eq:govern2}.

For elliptic interface problems, a common approach is to bias the fluxes according to a kind of harmonic weighting; in ref.~\cite{fluxx}, a stronger kind of biasing is advocated, denoted viscosity-upwinded weighting. Here, we adopt the same strategy and refer the reader to ref.~\cite{fluxx} for discussion, much of which is directly analogous to the multi-phase Stokes case, and to \iffalse the Supplementary Material \fi \cref{sec:ldg} with details on how the multi-valued numerical flux functions $\hat\mathbf u$ and $\hat\bm{\sigma}$ are chosen for interfacial mesh faces.

In addition to the application of viscosity-upwinded numerical fluxes, one other numerical technique is employed to improve the conditioning of the multi-phase Stokes equations. For a single-phase, constant-coefficient Stokes problem with viscosity $\mu$, i.e., $\{-\mu \nabla^2 \mathbf u + \nabla p = \mathbf f, -\nabla \cdot \mathbf u = f\}$, the largest positive eigenvalue of the discretised Stokes operator scales as $\mu/h^2$, whereas all negative eigenvalues scale inversely-proportional to $\mu$ (and independently of $h$). Thus, unlike a Poisson problem, whose conditioning is independent of an arbitrary multiple of its ellipticity coefficient, the conditioning of the Stokes problem worsens quadratically in $\mu$ as $\mu$ is made arbitrarily large. In essence, the two operators of the momentum equations, $\mu \nabla^2 \mathbf u$ and $\nabla p$, are not on equal footing if $\mu$ is large and $\mathbf u$ and $p$ are treated as independent variables. However, as noted earlier, from a physical point of view, the pressure scales as $p \sim \mu U/L$, and so the magnitude of $p$ depends on $\mu$. This apparent ill-conditioning of the Stokes operator can be easily remedied by rescaling one of the variables $\mathbf u$ or $p$ by $\mu$ or $1/\mu$, respectively, effectively recasting the Stokes equations into a unit-viscosity form, $-\nabla^2 \tilde{\mathbf u} + \nabla \tilde p = \tilde{\mathbf f}$.

This issue of ill-conditioning is exaggerated for the multi-phase Stokes equations, in which the largest positive eigenvalue scales as $\max_i \mu_i /h^2$. However, as motivated by the single-phase case, a simple rescaling of the solution variables $\mathbf u$ and $p$ can be used to mitigate the issue (see also refs.~\cite{FuruichiMayTackley2011,AksoyluUnlu2014,BorzacchielloLericheBlottiereGuillet2017}). In this work, we achieve this through a diagonal pre- and post-scaling of the Stokes operator, i.e., replace the Stokes operator in \eqref{eq:blockform} with
\[ \begin{pmatrix} \alpha & 0 \\ 0 & \beta \end{pmatrix} \begin{pmatrix} A & M {\mathcal G} \\ {\mathcal G}^\mathsf{T} M & -E \end{pmatrix} \begin{pmatrix} \alpha & 0 \\ 0 & \beta \end{pmatrix} \]
where $\alpha$ and $\beta$ are diagonal matrices whose entries equal $1/\sqrt{\mu}$ and $\sqrt{\mu}$, respectively. (More precisely, for every velocity field $\mathbf u$ and pressure field $p$, we have that $\smash{(\alpha \mathbf u)|_E = \mu_E^{-1/2} {\mathbf u}|_E}$ and $\smash{(\beta p)|_E = \mu_E^{1/2} p|_E}$ for every element $E \in {\mathcal E}$, where $\mu_E$ is the viscosity on element $E$.) This pre- and post-scaling is built into the multigrid schemes by replacing the Stokes operator with the scaled version on all levels of the hierarchy. Note that it is straightforward to adapt the original linear system to the scaled approach: instead of solving ${\mathcal A}x = b$, one instead solves for $\tilde {\mathcal A} \tilde x = \tilde b$, where $\tilde {\mathcal A}$ is the scaled Stokes operator and $\tilde b = \operatorname{diag}(\alpha,\beta) b$, whereupon solving for $\tilde x$ gives the original unscaled solution as $x = \operatorname{diag}(\alpha,\beta) \tilde x$.

\begin{figure}[t
\centering\sffamily\footnotesize%
(a) $\tfrac{\mu_1}{\mu_2} = 10^{-3}$\\[0em]%
\includegraphics[scale=0.91]{square-3-x}\\[-1em]%
(b) $\tfrac{\mu_1}{\mu_2} = 10^{+3}$\\[-1em]%
\includegraphics[scale=0.91]{square+3-x}\\[-1em]%
(c) $\tfrac{\mu_1}{\mu_2} = 10^{-6}$\\[-1em]%
\includegraphics[scale=0.91]{square-6-x}\\[-1em]%
(d) $\tfrac{\mu_1}{\mu_2} = 10^{+6}$\\[-1em]%
\includegraphics[scale=0.91]{square+6-x}%
\caption{Measured multigrid convergence rates when solving the multi-phase Stokes equations in viscous-stress form, wherein $\Omega_1 = (\tfrac14,\tfrac34)^d$ has viscosity $\mu_1$ and $\Omega_2 = (0,1)^d \setminus \overline{\Omega_1}$ has viscosity $\mu_2$, with the viscosity ratio as indicated, and periodic boundary conditions. In each case, 2D results are plotted in the left column, and 3D results in the right column.}%
\label{fig:square}%
\end{figure}

With the application of viscosity-upwinded numerical fluxes and the simple diagonal scaling to improve conditioning, \cref{fig:square} plots the measured multigrid convergence rates when solving the multi-phase Stokes equation in viscous-stress form with periodic boundary conditions, for the four different viscosity ratios considered, $\tfrac{\mu_1}{\mu_2} \in \{10^{-6}, 10^{-3}, 10^3, 10^6\}$. Compared to previous test problems, we observe marginally slower convergence rates, representative of the challenging Stokes problems at hand. Nevertheless, good convergence rates are obtained across significant viscosity ratios, i.e., typically 7--10 iterations of multigrid-preconditioned GMRES for a $10^8$ reduction in the residual norm.  Numerical experiments examining the order of accuracy (see \iffalse the Supplementary Material \fi \cref{sec:gridconv}) show that the velocity attains order $p+1$ in the maximum norm, while pressure attains order $p+\tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm. As in other test problems, a numerical boundary layer exists in the pressure field: in the present case, the layer is adjacent to the interface, but, as in other test problems, does not impact the optimal order accuracy of the velocity.

\section{Multigrid Efficiency for the Time-Dependent Stokes Equations}

So far, we have focused on the efficacy of the multigrid solver when applied to the steady-state Stokes equations. The time-dependent Stokes equations, however, may pose additional challenges, owing to the competing effects of the temporal derivative and viscous stress operator \cite{BenziGolubLiesen2005}. To illustrate, consider the time-dependent Stokes equations in the following form:
\begin{equation} \label{eq:timedep1} \left. \begin{aligned} \frac{\rho_i}{\delta} \mathbf u -\nabla \cdot \bigl(\mu_i (\nabla \mathbf u + \gamma \nabla \mathbf u^\mathsf{T}) \bigr) + \nabla p &= \mathbf f \\ -\nabla \cdot \mathbf \mathbf u &= f \end{aligned} \right\} \text{ in $\Omega_i$,} \end{equation}
where $\rho_i$ is a phase-dependent density and $\delta$ is a parameter proportional to the time step $\Delta t$ of a temporal integration method. Its LDG discretisation is a straightforward modification to the corresponding steady-state Stokes equations, and reads
\begin{equation} \label{eq:timedep2} \begin{pmatrix} \frac{1}{\delta} M_\rho + A & M {\mathcal G} \\ {\mathcal G}^\mathsf{T} M & -E \end{pmatrix} \begin{pmatrix} \mathbf u_h \\ p_h \end{pmatrix} = \begin{pmatrix} {\mathbf b}_{\mathbf u} \\ b_p \end{pmatrix}, \end{equation}
where $M_\rho$ is a $\rho$-weighted mass matrix, with all other operators unchanged.

When the viscous operator dominates (i.e., a small effective Reynolds number such that $\mu \delta/\rho$ is sufficiently large, depending on mesh resolution), the dominant operator is the steady-state Stokes equations, with a small $\rho/\delta$-weighted identity shift added to the viscous operator; in this case, one may expect a good steady-state Stokes solver to be effective. On the other hand, in the case when viscous effects are weak (i.e., a large effective Reynolds number such that $\mu \delta/\rho$ is small or the mesh is unable to resolve viscous effects), then \eqref{eq:timedep1} essentially reduces to a Helmholtz-Hodge projection problem (having strong connections to Chorin's projection method for solving the incompressible Navier-Stokes equations \cite{Chorin,GuermondMinevShen,sdpc}). In the latter case, a solver designed specifically for the steady-state Stokes equations may deteriorate.

In this work, a simple strategy is employed to automatically account for these two competing effects, resulting in fast multigrid solvers across a full range of Reynolds numbers. In essence, the strategy chooses the pressure penalty stabilisation parameter $\tau_p$ in \eqref{eq:eop} according to the expected scaling of the maximal eigenvalue of the operator $\tfrac{\rho}{\delta} {\mathbb I} - \nabla \cdot (\mu(\nabla + \gamma \nabla^\mathsf{T}))$.
In the steady-state case, e.g., $\rho = 0$, the penalty parameter should scale as $\tau_p \sim h/\mu$. However, in the degenerate time-dependent Stokes case, in which $\mu = 0$, \eqref{eq:timedep2} reduces to an LDG method, written in flux form, for computing the solution of a Poisson problem with operator $\smash{\nabla \cdot (\tfrac{\delta}{\rho} \nabla p)}$; the appropriate scaling of the pressure penalty stabilisation parameter in this case is then $\tau_p \sim \delta/(\rho h)$ \cite{ldg,fluxx}. Both of these scaling statements may be summarised as follows: the penalty parameter should scale such that $\tau \sim (h \Lambda)^{-1}$, where $\Lambda$ is the maximal eigenvalue of either the discretised operator $-\nabla \cdot (\mu (\nabla + \gamma \nabla^\mathsf{T}))$ (in the former case) or $\tfrac{\rho}{\delta} {\mathbb I}$ (in the latter case). This leads to a simple, but effective, idea to treat the general case: let $\tau \sim (h \Lambda)^{-1}$ where $\Lambda$ is the sum of the expected scalings of the maximal eigenvalues of the two operators. Using this idea, the definition of the pressure penalty parameter in \eqref{eq:eop}--\eqref{eq:taueqn} is replaced with the following relation:
\begin{equation} \label{eq:taugeneral} \tau_p = \Bigl( \frac{h \rho}{\tau_0 \delta} + \frac{\mu}{\tau h} \Bigr)^{-1}. \end{equation}
Note that \eqref{eq:taugeneral} reduces to the correct penalty parameter for the steady-state Stokes case when $\rho = 0$, i.e., $\tau_p = \tau h/\mu$, where $\tau$ is the multigrid-optimal parameter given in \cref{tab:optimaltau}; in the other extreme, when $\mu = 0$ (or $\delta$ is vanishingly small), \eqref{eq:taugeneral} reduces to $\tau = \tau_0 \delta / (h \rho)$, appropriate for a scalar Poisson problem with ellipticity coefficient $\tfrac{\delta}{\rho}$, where $\tau_0$ is a user-defined constant prefactor. In the general situation, i.e., between these two extremes, it is important to note the scaling of \eqref{eq:taugeneral} may change across the multigrid hierarchy---on a very fine mesh, $h$ is small so that the second term may dominate; on a coarse mesh, however, $h$ is large and the first term may dominate. This should be taken into account when building the multigrid method---in a purely-geometric approach, wherein $\mathcal A_h$ is explicitly built on every level, \eqref{eq:taugeneral} can be utilised directly; for the operator coarsening strategy,
a simple modification to existing schemes is discussed in \iffalse the Supplementary Material. \fi \cref{sec:op}.

The penalisation choice in \eqref{eq:taugeneral} blends across the two extremes of the time-dependent Stokes equations: a steady-state Stokes problem at one extreme, and a pure Poisson problem (written in flux form) at the other. Across a full range of Reynolds numbers, convergence results confirm optimal order of accuracy and that this choice of $\tau$ neither saturates the discretisation error nor underpenalises. As in the case of steady-state Stokes problems, we found that a V-cycle with three pre- and post-smoothing steps (used throughout this work) resulted in approximately optimal convergence speed, independent of the effective Reynolds number; occasionally a different combination of pre- and post-smoothing steps (at least two and at most four) can be slightly faster, but as mentioned earlier, the precise optimal value depends on a host of implementation, hardware, and problem-dependent factors.

\subsection{\hspace{-0.04em}Single-phase time-dependent Stokes equations}
\label{prob:timedep1}

\begin{figure
\centering\sffamily\footnotesize%
(a) $\mu = 10^{-2}$ corresponding to $\textsf{Re} \approx 100$\\[0.2em]
\includegraphics[scale=0.91]{time-single-phase-2-x}\\[1em]%
(b) $\mu = 10^{-4}$ corresponding to $\textsf{Re} \approx 10,\!000$\\[0.2em]
\includegraphics[scale=0.91]{time-single-phase-4-x}%
\caption{Measured multigrid convergence rates when solving the time-dependent, single-phase Stokes equations in standard form with $\mu$ as indicated and $\delta = 0.1 h$, together with Neumann boundary conditions.}%
\label{fig:time-single-phase-24}%
\end{figure}

To demonstrate multigrid performance on time-dependent single-phase Stokes problems, we consider two effective Reynolds numbers $\textsf{Re} = \rho U L/\mu$, one in which $\textsf{Re} \approx 100$, representing a viscous-dominated case (but where the time-derivative operator nevertheless influences performance characteristics) and the other with $\textsf{Re} \approx 10,\!000$ (wherein the time-derivative operator definitively dominates). In both cases, the velocity and length scales are unitary, $U = 1$, $L = 1$; density is set to $\rho = 1$, while viscosity satisfies $\mu = 10^{-2}$ in the former case and $\mu = 10^{-4}$ in the latter. In addition, we set $\delta = 0.1 h$, where $h$ is the element size on the finest-level mesh, representing a typical scenario of applying the time-dependent Stokes equations in a temporal integration method with CFL about $0.1$. As before, multigrid efficiency is quantified through the average convergence rate $\rho$; \cref{fig:time-single-phase-24} plots the results (note the magnified vertical axis).\footnote{In the results of this section, we have reincluded the case of $p = 1$, mainly to serve as a point of interest: specifically, for steady-state Stokes problems, we noted in \cref{sec:timeindep} that sub-optimal multigrid efficiency may occur when $p = 1$; for unsteady Stokes problems, and depending on the effective Reynolds number, ideal multigrid convergence can be restored in this case, as seen in \cref{fig:time-single-phase-24,fig:water-bubble}.} In \cref{fig:time-single-phase-24}(a), we observe an upwards trend in $\rho$ as the mesh is refined; this corresponds to the fact that as the mesh is refined, eventually the viscous operator will definitively dominate and multigrid convergence rates similar to the steady-state Stokes equations will be attained. In the weakly viscous case with $\textsf{Re} \approx 10,\!000$, \cref{fig:time-single-phase-24}(b) shows exceptionally fast multigrid convergence rates, with $\rho \approx 0.01$, corresponding to needing only 4 GMRES iterations to achieve a factor $10^8$ reduction in the initial residual. In the context of solving the incompressible Navier-Stokes equations, these results suggest that a fast non-stationary Stokes solver may outperform a fast projection method solver; further remarks on this topic are given in the conclusions. Meanwhile, grid convergence experiments examining the order of accuracy (see \cref{sec:gridconv}) show a departure from the typical results seen elsewhere in this work. In the weakly-viscous case with $\textsf{Re} \approx 10,\!000$, the velocity is order $p+1$ in the maximum norm, and so is the pressure, despite the presence of boundary conditions. This apparent ``super convergence'' in pressure is attributed to the property that, for very large Reynolds numbers, the time-dependent Stokes equations nearly reduce to a Helmholtz-Hodge projection, where one may naturally expect to attain optimal order accuracy in the pressure field, see, e.g., ref.~\cite{ImplicitMeshPartOne}. On the other hand, in the viscous-dominated case with $\textsf{Re} \approx 100$, the pressure reduces to order $p+\tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm, while the velocity maintains order $p+1$ in the maximum norm. This is perhaps expected, based on the order of accuracy results reported elsewhere in this work, together with the property that the viscous-dominated case essentially represents a mildly perturbed stationary Stokes problem.

\subsection{Multi-phase time-dependent Stokes equations}
\label{prob:timedep2}

In our last two examples, we combine the challenging aspects of a multi-phase Stokes problem together with non-stationary effects and consider a problem in which both the density $\rho$ and viscosity $\mu$ have discontinuities several orders in magnitude across an embedded interface. The specific parameters considered correspond to two scenarios: a water bubble surrounded by gas, and a gas bubble surrounded by water. Specifically, $\rho_{\text{water}} = 1$, $\rho_{\text{gas}} = 0.001$, $\mu_{\text{water}} = 1$, and $\mu_{\text{gas}} = 0.0002$ (approximately accurate values for water and air at ambient temperature in CGS units). The radius of the bubble is $0.3$, centred in a unit square/cubic domain, $\Omega = (0,1)^d$; equations \eqref{eq:timedep1} with $\gamma = 1$ and \eqref{eq:govern2}--\eqref{eq:govern3} are solved with velocity Dirichlet boundary conditions and stress jump conditions across the gas-water interface. As in the test case on curved domain geometry, a level set function describing the interface geometry is used together with an implicit mesh DG framework to create a semi-unstructured, interface-conforming mesh; an example is shown in \cref{fig:circle}(b). As before, the time step is set equal to $\delta = 0.1 h$, where $h$ is the typical element size on the finest mesh, representing the application of a time stepping method with CFL number about $0.1$.

\begin{figure
\centering\sffamily\footnotesize%
(a) Water bubble surrounded by gas\\[0.2em]
\includegraphics[scale=0.91]{water-bubble-x}\\[1em]%
(b) Gas bubble surrounded by water\\[0.2em]
\includegraphics[scale=0.91]{air-bubble-x}%
\caption{Measured multigrid convergence rates when solving the time-dependent multi-phase Stokes equations in viscous-stress form, with $\rho_{\text{water}} = 1$, $\rho_{\text{gas}} = 0.001$, $\mu_{\text{water}} = 1$, $\mu_{\text{gas}} = 0.0002$, and $\delta = 0.1 h$, together with velocity Dirichlet boundary conditions. Here, $n$ denotes the number of cells in the background uniform Cartesian $n\times n\,(\,\times n)$ grid used to build the corresponding implicitly defined mesh; see \cref{fig:circle}(b).}%
\label{fig:water-bubble}%
\end{figure}

This example combines three distinct but complementary strategies developed in prior examples: (i) viscosity-upwinded numerical fluxes, to robustly and accurately handle the large jump in viscosity across the interface; (ii) diagonal pre- and post-scaling of the Stokes operator, to remove the unnecessary ill-conditioning caused by viscosity coefficients differing several orders in magnitude; and (iii) a pressure penalty stabilisation parameter controlled by \eqref{eq:taugeneral}, to automatically adjust penalisation behaviour between the two extremes of the time-dependent Stokes equations. As such, this example serves to demonstrate a variety of subtleties in efficiently solving the multi-phase time-dependent Stokes equations, but as shown in the results plotted in \cref{fig:water-bubble}, combined together, one can attain highly efficient multigrid solvers. For the particular scales chosen in this problem, the water-bubble and gas-bubble problems correspond to a unit-order Reynolds number. As such, grid convergence analyses confirms the expectation that the velocity attains order $p+1$ in the maximum norm, and pressure attains order $p+\tfrac12$ in the $L^2$ norm and order $p$ in the maximum norm.

\section{Concluding Remarks}

In this paper, we devised efficient geometric multigrid solvers for the Stokes equations discretised by local discontinuous Galerkin (LDG) methods. The approach follows standard geometric multigrid concepts, utilising a V-cycle and a simple block Gauss-Seidel relaxation method free from under-relaxation parameters. With a suitably-chosen pressure penalty stabilisation parameter, a wide array of tests showed that the Stokes multigrid solver can match the speed of classical geometric multigrid methods for Poisson problems \cite{dgmg,fluxx}. For example, typical convergence rates ranged from $\rho \approx 0.05$ to $0.1$, corresponding to needing about six to eight iterations for a $10^8$ reduction of the residual; fewer iterations are required for unsteady time-dependent Stokes equations.
To implement the multigrid algorithm, one possibility is to explicitly build the mesh and associated LDG operators on every level of the hierarchy; an alternative method, not requiring the explicit formation of coarse meshes, can be implemented based on the operator coarsening algorithms detailed in \iffalse the Supplementary Material.\fi \cref{sec:op}. In addition, we also extended the LDG methods of ref.~\cite{CockburnKanschatSchotzauSchwab2002} to variable-viscosity and multi-phase problems exhibiting interfacial stress jump conditions; across all test problems, grid convergence analyses demonstrated order $p+1$ accuracy in the maximum norm for the computed velocity field, and at least order $p$ accuracy in the maximum norm for pressure.

Key findings of this work include the following aspects:
\begin{itemize}
\item Multigrid efficiency depends on an appropriate choice of the pressure penalty stabilisation parameter $\tau_p$. In general, for the steady-state Stokes equations, $\tau_p$ scales linearly with the local element size $h$ and inversely-proportional to the local viscosity, leaving the end-user to define the prefactor $\tau$ in the formula $\tau_p = \tau\, h/\mu$. \Cref{tab:optimaltau} provides values of $\tau$ resulting in optimal multigrid convergence rates, showing that $\tau$ depends on the polynomial degree of the DG space, the spatial dimension, and whether the standard form or viscous-stress form of the Stokes equations are being solved. \iffalse Table SM1 in the Supplementary Material \fi \Cref{tab:tauwindow} shows that an associated range of values exist for $\tau$ within which multigrid performance is very close to optimal. (See also Kanschat~\cite{Kanschat2005}, in which the performance of an approximate Schur form block preconditioner for LDG was also noted to depend crucially on the pressure penalty parameter, although the scaling of $\tau$ determined therein is very different to the results of the present work.)

\item In the case of the time-dependent Stokes equations, \eqref{eq:taugeneral} implements an effective strategy to automatically adjust the penalty parameter based on the competing effects of the viscous and time-derivative operators.
\item For multi-phase problems, in which the viscosity may exhibit jumps several orders in magnitude across an embedded interface, one should bias the numerical fluxes (here implemented with viscosity-upwinded weighting) to reflect the physical nature of the interfacial stress and velocity jump conditions. In addition, it can be advantageous to employ a diagonal preconditioning of the Stokes operator to correct the potential adverse affects large viscosities may have on the spectral characteristics of the saddle point problem.
\end{itemize}

A formula for the optimal pressure penalty parameter, $\tau$, is not known at this time. In this work, a simple one-dimensional parameter sweep was used to find $\tau$ for particular choices of $p$, but a theoretical result would be ideal. Here, it may be possible to use analytical tools such as convergence criteria for block Gauss-Seidel methods \cite{ElsnerMhermann1991}, local mode analyses \cite{Sivaloganathan1991,GasparNotayOosterleeRodrigo2014} or local Fourier analyses \cite{HeMacLachlan2017,FarrellHeMacLachlan2019} to help determine a formula. On a similar note, it may also be possible to derive analytical proofs of the convergence of the overall multigrid method; see also, e.g., \cite{SchoberlZulehner2003,Manservisi2006,DrzisgaJohnRudeWohlmuthZulehner2018,KanschatMao2015,HongKrausXuZikatanov2016}.

Other areas of study include the following. In this work, we mainly considered structured meshes (such as Cartesian grids) and semi-unstructured, non-conforming implicitly defined meshes resulting from cell merging procedures. Efficacy of the multigrid Stokes solver on fully unstructured meshes is also worthy of examination--for highly anisotropic meshes, one may need to group elements into clusters for the block Gauss-Seidel method to be effective \cite{ThompsonFerziger1989,OosterleeWesseling1993,WobkerTurek2009}. Another possibility is to accelerate the Gauss-Seidel method through low-degree Chebyshev iterations, as was noted by Farrel \textit{et al} \cite{FarrellHeMacLachlan2019} for Vanka-type smoothers. Meanwhile, a wide variety of work on developing multigrid methods for the Stokes equations report that W-cycles can be more effective than V-cycles, or even variable V-cycles which change the smoothing counts between levels of the hierarchy; see, e.g., \cite{GasparNotayOosterleeRodrigo2014,DrzisgaJohnRudeWohlmuthZulehner2018,FarrellHeMacLachlan2019,KanschatMao2015,AdlerBensonMacLachlan2016}. Although a V-cycle with fixed pre-smoothing and post-smoothing steps was found highly effective in this work, these alternative strategies could prove useful for different applications. In the presented results, we also mentioned that the block Gauss-Seidel method is less effective when $p = 1$, i.e., for a bilinear or trilinear LDG discretisation of the Stokes equations; although $p = 1$ is rarely of interest for DG methods, one possibility here is to use a $p > 1$ method as a preconditioner for the $p = 1$ system, or perhaps cluster elements into larger blocks for the Gauss-Seidel method. Meanwhile, the results of this paper could also be used to inform the design of algebraic multigrid methods \cite{Notay2017}.

Finally, we remark that the measured multigrid convergence rates for the time-dependent Stokes equations are particularly encouraging, indicating a strong potential for developing fast solvers for the general incompressible Navier-Stokes equations. In particular, preliminary work indicates that a fast Stokes solver could outperform the well-known and widely-applied projection method of Chorin \cite{Chorin,GuermondMinevShen}, and could be integrated into arbitrary-order time-stepping methods for Navier-Stokes \cite{sdpc}; this will be further investigated in future work.

\subsection*{Acknowledgements} This research was supported by the Applied Mathematics Program of the U.S.~DOE Office of Advanced Scientific Computing Research under contract number DE-AC02-05CH11231. Some computations used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S.~DOE Office of Science User Facility operated under Contract No.~DE-AC02-05CH11231.

\clearpage

redpajama_set_name: RedPajamaArXiv


text:
1. Supervision:
The Hirer will, during the period of the hire, be accountable as the 'Responsible Person' for supervision of the premises, the fabric and the contents, their care, safety from damage however slight, and behaviour of all persons using the premises, whatever their capacity.
2. Use of premises:
a) The Hirer shall not use the premises for any purpose other than that described in the hiring agreement and shall not sub-hire or use the premises or allow the premises to be used for any unlawful purpose or in any unlawful way, nor do anything or bring onto the premises anything which may endanger the same or render invalid any insurance policies in respect therof.
b) No fireworks are permitted in the vicinity of the Hall, the grassed area and Car Park. Express permission is required for the erection of marquees, tents or other structures at the rear of the Hall. This also applies to the use of BBQs or other cooking arrangements.
c) Use of laser beams or artificial smoke is not permitted inside the building.
d) Where children are present in the letting and have access to the external area to the rear of the Hall via the Emergency Door, users must ensure that the gates are secured and children do not leave this area.
3. Care of premises and property:
No application will be accepted from any person under the age of 18. The Hirer will be responsible for any loss, theft or damage occasioned to the premises or equipment during the term of hiring and any such loss will be charged on the Hirer.
4. Liability:
The Trustees will not accept liability for any loss or injury sustained by the hirer, his/her guests or other persons attending the function for which the premises have been hired unless and only if it is attributable to the Trustees' negligence or that of their representatives.
5. Access by Trustees:
The Hirer shall on demand afford access to any of the Trustees or their agents and conform to any reasonable instructions that they may give.
6. Consumption of alcohol:
The premises are licensed for the sale of alcohol for hirers who are Crawley village residents. This is subject to the Trustees' discretion and agreement.
For non resident hirers who wish alcohol to be sold on the premises, a Temporary Events Notice must be applied for from the Winchester City Council Licensing authority. However these are limited in number and any such application must first have the approval of the Trustees.
In order to comply with the Licensing Act 2003 the hirer, or person nominated by the hirer, will take full responsibility for supervising the sale and consumption of alcohol on the premises and for ensuring that no alcohol is sold or consumed by a person under the age of 18 years.
The Crawley Village Hall Management Committee (ie the Trustees) have adopted the "Challenge 21" scheme. This requires the hirer to challenge any person that appears under the age of 21 for photographic ID prior to any sale of alcohol. If the person cannot prove they are over 18 years of age, they must be refused the sale of alcohol.
7. Music licences:
The Hall does not hold either a PPL licence for the playing of sound recordings where the music is still in copyright, or a PRS for Music licence for live performance of music under copyright. The Trustees' understanding is that such licences are not required for private functions (eg weddings, birthdays, christenings) but are required where the Hall is hired for public or commercial functions where music under copyright will be played and/or performed. For such functions, the Trustees will ask for confirmation from the hirers of the Hall that either: a) no copyrighted music will be played or performed; or b) that the relevant licence(s) are held by the hirers of the Hall.
8. Removal of property:
Property belonging to the Trustees must not be removed from the premises under any circumstances. No equipment is to be taken outside the building except by specific approval of the Trustees or their agents.
9. Maximum numbers:
The Hall is licensed in respect of Music, Singing, Dancing, Theatrical Productions and Films. The maximum number allowed in the Main Hall at any one time is 200, and in the Finch (Meeting) Room 40.
10. Smoking:
Smoking is not permitted within any part of the building.
11. Safety precautions:
All doors, gangways, exit notices and exits must be kept clear at all times. Fire exit doors must not be locked while the Hall is in use. Fire precautions must be observed and the Hirer must undertake to know the location of all Fire Extinguishers and Exit Doors. (See also Item 19 of the Terms of Letting regarding Access).
12. Safety of electrical equipment:
It is a requirement that the Hirer holds a valid Electrical Test Certificate for each electrical appliance or cable brought into the premises during the hire period.
13. Cancellation of booking:
If the Hirer cancels a booking after payment has been made, then 50% of the Total Hire Charge Less Discount as indicated on the Hire invoice is non-refundable. Should the cancellation be made less than fourteen days before the date of the first booking indicated on the invoice, then 100% of the Total Hire Charge Less Discount will become non-refundable.
14. Car parking:
Car parking is limited to a small number of spaces indicated in an area opposite the main ramped entrance to the Hall. There are also two car parking spaces available for the disabled immediately adjacent to the ramp. Parking is permitted on the Main Street but care must be taken not to block any entrances to residents' drives. (No cars must be parked FOR ANY REASON on the access road to the Cottages located at the rear of the Hall. This restriction also applies to the parking areas and circulation space in front of the Cottages. The Hirer must acquaint guests arriving by car of these restrictions).
15. Animals:
The Hirer shall ensure that no animals (including birds and reptiles) except Guide Dogs are brought into the building, other than with special approval of the Trustees or their agents. At no time shall animals be allowed to enter the kitchen.
16. Health & hygiene:
The Hirer shall, if preparing, serving or selling food, observe all relevant food, health and hygiene legislation and associated regulations.
17. Accidents & dangerous occurences:
The Hirer must report all accidents involving injury to the public to the Trustees or their agents as soon as possible. This also applies to the failure of equipment belonging to the Hall or brought in by the Hirer.
18. Compliance with legislation & licenses:
The Hirer shall, as the 'Responsible Person', ensure that any activities for children under eight years of age comply with the provisions of the Children's Acts and that only fit and proper persons have access to the children. Equally that all conditions and restrictions of theatre, film and other entertainment licences issued by Winchester City Council are adhered to.
19. Noise & nuisance:
The Hirer shall ensure that the minimum of noise is made during the period of the hire, including arrival and departure. (The Hall is situated in a quiet residential area of the Village and the utmost care shall be taken to avoid undue and unnecessary disturbance to the nearby residents). MUSIC must end by 11.15pm on Saturday and 10.45pm on other days. In no circumstances are fireworks to be set off in the vicinity of the Hall, grassed area or Car Park.
20. Access to the Hall:
The main ramped access to the Hall must be used at all times. The Access Code will be advised to the Hirer when the hire charge is paid. (The old Access Door from the Main Street will not normally be used as an entrance to the Hall. However, as an identified FIRE EXIT, when the Main Hall is in use, this door will still be used as an Emergency Exit and must be kept as such in readiness).
21. Block bookings:
Special rates are available to organizations that wish to hire the Hall on a regular basis. Special rates will be negotiated with such organizations that are judged to fall within this category.
22. Charges for hire of the Hall:
A schedule of standard charges is available that provides for varying rates for hire of different parts of the Hall. (In most cases of hire, kitchen access will be freely available). There are also different rates for week-end and week-day hiring. Discounts are available for the hire by Village Residents for essentially domestic and family occasions and for regular users. Whole-day hire of the Hall will also attract a reduction from the standard rates.
23. Kitchen facilities:
The kitchen is modern and well equipped with two good size electric ovens, an associated five ring gas hob and extractor fan. A hot water boiler and kettles are also available. There is a refrigerator, a dish-washer and hot water on tap. See Item 23 for details of other equipment. The kitchen facilities are available to most hirers but it is essential that the various items of equipment are used with great care, particularly the use of the ovens and gas hob, The 'Responsible Person' must ensure that users, including outside caterers, are aware of their responsibilities in this respect. All equipment must be turned off before finally vacating the Hall.
24. The Village Hall inventory:
The Hall is equipped with an ample number of tables and chairs which are stored off the Main Hall. Also available in the kitchen are every-day items of crockery (cups, saucers and side plates suitable for serving coffee, tea, sandwiches and cakes). These, together with a restricted number of goblets and tumblers, are readily available in the kitchen cupboards.
Additionally, and available by arrangement with the Trustees, are up to 100 place settings (dinner, pudding and side plates together with associated cutlery). For these to be available the Trustees will need to be informed of the number of place settings required. Every care must be taken to ensure that all the items made available are returned . These will be checked and any items which are missing must be paid for.
25. Tidying the Hall after use:
Every care is taken by the Trustees to ensure that the Hall is in extremely good order prior to a letting and that all items of equipment, including the toilet accommodation (male, female and disabled), are maintained to a very high standard. The Hirer is requested to ensure that the condition of the Hall as far as possible is left as found on arrival and that all floors are swept, paying particular attention to the state of entrance hall when adverse weather conditions occur. All rubbish should be placed in the black plastic bags provided (in the kitchen) and taken away by the Hirer. A deposit of £50 is to be paid with each booking for Wedding receptions (or similar events). This will be refunded to the Hirer if the Hall is left clean and tidy after the event.

redpajama_set_name: RedPajamaCommonCrawl