We present the full approximation scheme constraint decomposition (FASCD) multilevel method for solving variational inequalities (VIs).  FASCD is a common extension of both the full approximation scheme (FAS) multigrid technique for nonlinear partial differential equations, due to A. Brandt, and the constraint decomposition (CD) method introduced by X.-C. Tai for VIs arising in optimization.  We extend the CD idea by exploiting the telescoping nature of certain function space subset decompositions arising from multilevel mesh hierarchies.  When a reduced-space (active set) Newton method is applied as a smoother, with work proportional to the number of unknowns on a given mesh level, FASCD V-cycles exhibit nearly mesh-independent convergence rates, and full multigrid cycles are optimal solvers.  The example problems include differential operators which are symmetric linear, nonsymmetric linear, and nonlinear, in unilateral and bilateral VI problems.

multigrid, variational inequalities, full approximation scheme, constraint decomposition

65K15, 35M86, 90C33

§ INTRODUCTION

The constraint decomposition (CD) methods of Tai [32] are designed to solve variational inequality (VI) problems arising as the Karush–Kuhn–Tucker optimality conditions for the minimization of a convex functional over a convex set.  In a finite element (FE) context, these methods converge for multilevel and domain decompositions.  The multilevel CD method has almost grid-independent error bounds, and near-optimal complexity, for elliptic, linear obstacle problems <cit.>, <cit.>.

Here we extend Tai's multilevel CD method by replacing the coarser-level corrections with a full approximation scheme (FAS; [5, 8]) approach.  Specifically, our full approximation scheme constraint decomposition (FASCD) method enjoys the following advantages:

* The FAS coarse corrections are adapted for general VI problems. The method requires neither an objective function nor linearity of the residual.

* Box constraints, i.e. both upper and lower obstacles, are allowed.

* The method is strictly admissible, and so can be applied even if the residual operator is defined only on those functions which satisfy the constraints.

* V-cycles (Algorithm <ref>) and full multigrid (FMG; [34]) cycles (Algorithm <ref>) are proposed.

* We identify a “telescoping” subset property which can be exploited on the up-smoothing side of a V-cycle for greater performance.

* The method is (nonlinear) smoother-agnostic.  When a (mesh-dependent) reduced-space Newton method [3] is used as relaxation, we demonstrate that application of FMG is an optimal solver for the examples in Section <ref>.

Prior multilevel algorithms for VIs have limitations compared to the FASCD framework.  The projected full approximation scheme (PFAS) algorithm of Brandt & Cryer [6] applies to linear complementarity problems and it uses projected (truncated) Gauss–Seidel smoothers.  While FASCD reduces to general, nonlinear FAS for PDEs when the inequality constraints are removed, PFAS reduces to linear multigrid.  Another work which combines FAS and VIs is [24], but the problems solved, and the approach, are completely different.  For an optimal control problem with an elliptic VI constraint, [24] applies FAS to a system of PDEs derived from penalizing the VI and the control constraints; their V- and W-cycles are not altered in any manner to maintain admissibility.  On the other hand, the multilevel CD iteration of Tai [32] is formulated in terms of an objective function, so as stated it only applies to optimization-type problems, and furthermore each of its coarser-level problems refers to the finest-level discretization.  The latter limitation is removed by Gräser & Kornhuber's <cit.> reformulation using defect constraints and monotone injection operators; see Section <ref> below.  The monotone multigrid methods of Kornhuber and co-authors [22, 29] appear to be formulated only with Gauss–Seidel type smoothers, and may have been applied only to optimization-type VI problems <cit.>.  Furthermore, unlike the fastest variants of monotone multigrid, no zeroing of the basis functions associated with the active set (truncation) is employed in FASCD, and this simplifies the implementation.  Nonetheless FASCD coarse grid corrections can update the estimate of the active set on the fine grid, and indeed the numerical evidence indicates that FASCD coarse corrections are effective even when the coarse meshes are far too coarse to resolve the active set.  On a problem with a complicated coincidence set, for which the truncated Newton/monotone multigrid methods were apparently designed <cit.>, the performance of FASCD is at least as good (Example <ref>).  A limitation of FASCD is that, as formulated, it only solves bound-constrained VI problems.

In our experiments, FASCD V-cycles exhibit essentially the same mesh-independent convergence rates for VI problems that FAS multigrid V-cycles exhibit for the corresponding unconstrained PDE problems.  The FMG results in Section <ref> show near-textbook multigrid efficiency, yielding solutions within discretization error in work comparable to a few smoother sweeps on the finest-level mesh [7].  Specifically, we see this performance for classical (unilateral, Laplacian) obstacle problems, $p$-Laplacian obstacle problems, and box-constrained advection-diffusion equations.  Furthermore, V-cycle results are largely insensitive to the geometric complexity of the coincidence set and free boundary.

All iterates in FASCD, on all mesh levels, are admissible, and thus the nonlinear operator in the VI needs only be defined for admissible states.  For example, we solve an ice sheet geometry problem (Examples <ref> and <ref>) in which the ice flow formulas are only meaningful for surface elevation iterates that do not penetrate the bedrock.  (Non-admissible methods, including semi-smooth methods [3], would require unnatural modifications of the operator formula.)  FASCD solves this low-regularity ice sheet problem, wherein the solution has unbounded gradient at the free boundary, by a small, mesh-independent number of FMG iterations (Example <ref>). The same Example investigates the parallel weak scaling [11] of our implementation.

The paper is organized as follows.  In Section <ref> we recall the theory of coercive and monotone VIs, and describe several motivating examples.  Section <ref> reviews CDs and the iterations they induce.  Section <ref> sets up multilevel FE hierarchies, and then Section <ref> extends the CD method via defect constraints generated by monotone injection operators [22].  Here we propose an apparently new understanding of telescoping sets for up-smoothing, based on an “incomplete” CD and an associated iteration.  Section <ref> then states the V-cycle, with Section <ref> addressing convergence criteria, FMG cycles, and our smoother choice.  Section <ref> gives numerical results from our open source Python/Firedrake [31] implementation.

§ VARIATIONAL INEQUALITIES

Let $\cV$ be a real reflexive Banach space with norm $\|\cdot\|$ and topological dual space $\cV'$.  Denote the dual pairing of $\phi \in \cV'$ and $v\in\cV$ by $\ip{\phi}{v} = \phi(v)$, and define the (Banach space) norm on $\cV'$ by $\|\phi\|' = \sup_{\|v\|=1} |\!\ip{\phi}{v}\!|$.  Let $\cK \subset \cV$ be a nonempty, closed, and convex subset, the constraint set; elements of $\cK$ are said to be admissible.  For a continuous, but generally nonlinear, operator $f:\cK \to \cV'$ and a linear source functional $\ell\in \cV'$ we consider the following variational inequality (VI) problem: find $u^\star\in \cK$ such that

\begin{equation}
\ip{f(u^\star)}{v-u^\star} \ge \ip{\ell}{v-u^\star} \quad \text{for all } v\in \cK. \label{eq:vi}
\end{equation}

Because $f$ is (generally) nonlinear, the source term $\ell$ is not strictly needed in stating this class of problems—by redefining $f$ we may take $\ell=0$—but introducing $\ell$ clarifies the algorithms which follow.

VI (<ref>) generalizes nonlinear systems of equations $f(u^\star)=\ell$ to problems where $u^\star$ is also constrained to be in $\cK$.  Pretending the dual pairing is an inner product, by (<ref>) the angle between $f(u^\star)-\ell$ and any arbitrary vector $v-u^\star$ pointing from $u^\star$ into $\cK$ is at most $90^\circ$; the dual vector $f(u^\star)-\ell$ must point into $\cK$.  If $u^\star$ is in the interior of the constraint set then (<ref>) implies $f(u^\star)=\ell$.
Variational inequalities may generally permit multiple solutions [19].

The following definitions are standard [28].

An operator $f:\cK \to \cV'$ is monotone if

\begin{equation}
\ip{f(u)-f(v)}{u-v} \ge 0 \qquad \text{for all } u,v \in \cK, \label{eq:monotone}
\end{equation}

strictly monotone if equality in (<ref>) implies $u=v$, and coercive if there exists $w \in \cK$ so that

\begin{equation}
\frac{\ip{f(u)-f(w)}{u-w}}{\|u-w\|} \to +\infty \qquad \text{as } \|u\|\to +\infty. \label{eq:coercive}
\end{equation}

It is well-known that if $f:\cK \to \cV'$ is continuous, monotone, and coercive then VI (<ref>) has a solution <cit.>, and that the solution is unique when $f$ is strictly monotone.  As is standard in the calculus of variations [17], coercivity permits a compactness argument for unbounded sets $\cK$; recall that the closed and bounded subsets of a reflexive Banach space are weakly compact.  Some VIs solved in this paper satisfy a stronger inequality than (<ref>).

Let $p>1$.  The map $f:\cK \to \cV'$ is $p$-coercive if there exists $\kappa>0$ such that

\begin{equation}
\ip{f(u)-f(v)}{u-v} \ge \kappa \|u-v\|^p \qquad \text{for all } u,v \in \cK. \label{eq:pcoercive}
\end{equation}

It is easy to see that if $f$ is $p$-coercive then it is monotone, strictly monotone, and coercive.  Thus if $f:\cK \to \cV'$ is continuous and $p$-coercive then there exists a unique $u^\star\in \cK$ solving VI (<ref>).

Let $\Omega \subset \RR^d$ denote a bounded, open set with piecewise-smooth boundary.  Sobolev spaces [17] will be denoted by $W^{k,p}(\Omega)$, for integer $k$ and $1\le p \le \infty$, with $W^{1,p}_0(\Omega)$ denoting the subspace with zero trace.  The following example includes the classical and $p$-Laplacian obstacle problems [14].

Let $1<p<\infty$.  For $u,v \in \cV = W^{1,p}_0(\Omega)$, define $f:\cV \to \cV'$ by

\begin{equation}
\ip{f(u)}{v} = \int_\Omega |\grad u|^{p-2} \grad u \cdot \grad v\dx, \label{eq:plaplacian}
\end{equation}

a continuous map <cit.>.  It follows from inequalities for $|\cdot|^p$ in $\RR^d$, and the Poincaré inequality, that $f$ is $p$-coercive [14]; <cit.>.

For monotone operators $f$, VI (<ref>) generalizes the problem of minimizing a convex function over $\cK$.  Suppose $F:\cK \to \RR$ is lower semi-continuous and (Gâteaux) differentiable with continuous derivative $F':\cK \to \cV'$.  Then $F$ is convex if and only if $F'$ is monotone <cit.>.  Furthermore, Proposition II.2.1 in [15] shows that if $F$ is convex then (<ref>) holds for $f=F'$ if and only if

\begin{equation}
u^\star = \argmin_{v\in\cK} F(v) - \ip{\ell}{v}. \label{eq:minimization}
\end{equation}

The CD methods of Tai [32] address problem (<ref>); the analysis in [32] supposes that an objective $F$ exists with $F'$ 2-coercive.  For Example <ref> we may define

\begin{equation}
F(v) = \frac{1}{p} \int_\Omega |\grad v|^p\dx. \label{eq:plaplacianobjective}
\end{equation}

Then $F$ is a convex functional and $F'(v) = f(v)$ for $f$ given in (<ref>).  For any closed and convex $\cK\subset \cV$, VI problem (<ref>) is then equivalent to optimization problem (<ref>).

Next we give two examples which are not of optimization type (<ref>).  The first is an advection-diffusion problem; it requires the following Lemma.

[16]   Suppose $\bX :\Omega \to \RR^d$ is a bounded and boundedly-differentiable vector field on $\Omega$ with $\Div \bX=0$.  For $u,v \in W^{1,2}(\Omega)$ let $b(u,v) = \int_\Omega (\bX \cdot \grad u) v\dx$.  Then $b(u,u) = \frac{1}{2} \int_{\partial \Omega} u^2 \bX\cdot \bn\dx$ where $\bn$ is the outward normal on $\partial \Omega$.

Integration by parts gives $b(u,v) = - b(v,u) + \int_{\partial \Omega} uv \bX\cdot \bn\dx$.

Suppose $\partial\Omega$ is partitioned into Dirichlet and Neumann portions, $\partial\Omega = \partial_D\Omega \cup \partial_N\Omega$, with $\partial_D\Omega$ of positive Hausdorff measure.  Let $\cV \subset W^{1,2}(\Omega)$ be the space of functions which have zero trace on $\partial_D\Omega$.  Consider a velocity field $\bX$ on $\Omega$ satisfying the conditions of Lemma <ref>, and assume that the flow is outward on the Neumann boundary: $\bX \cdot \bn \ge 0$ on $\partial_N\Omega$.  For $u,v \in \cV$, $\eps>0$, and $\phi \in \cV'$, define

\begin{equation}
\ip{f(u)}{v} = \eps \left(\grad u, \grad v\right)_{L^2(\Omega)} + b(u,v) - \ip{\phi}{v}. \label{eq:advectiondiffusion}
\end{equation}

With this operator, consider VI (<ref>) with $\ell = 0$ for any closed and convex $\cK \subset \cV$. The interior condition in strong form is the following linear advection-diffusion equation

\begin{equation}
-\eps \grad^2 u + \bX\cdot \grad u = \phi.
\label{eq:advectiondiffusionstrong}
\end{equation}

On the other hand, it is easy to see that $|\!\ip{f(u)}{v}\!| \le \left((\eps + \|\bX\|_\infty) \|u\| + \|\phi\|'\right) \|v\|$, so that $f:\cK \to \cV'$ is continuous.  Lemma <ref> says that the bilinear form $b(u,v)$ is skew-symmetric up to a term which is nonnegative by the outward flow assumption.  From the Poincaré inequality,

\begin{align*}
\ip{f(u)-f(v)}{u-v} &= \eps \int_\Omega |\grad u - \grad v|^2\dx + b(u-v,u-v) \\
&= \eps \int_\Omega |\grad u - \grad v|^2\dx + \frac{1}{2} \int_{\partial_N\Omega} (u-v)^2 \bX\cdot\bn \ge \eps C \|u-v\|^2.
\end{align*}

Thus $f$ is 2-coercive, and so the associated VI problem is well-posed. If $\bX \ne 0$ then $f\ne F'$ for any objective $F$, because $\ip{f(u)}{v}$ does not possess the symmetry of a gradient: for $u,v$ which are zero on $\partial \Omega$, note that $\ip{f(u)}{v} - \ip{f(v)}{u} = -2 b(u,v)$.  References [10, 13] consider similar advection-diffusion VI problems over $\cK = \{v\ge 0\}$.

The next VI example is a model for steady ice sheet geometry in a given climate.  It is fully nonlinear and not of optimization type.  Related VI problems arise for any fluid layer subject to processes which can remove fluid mass [10].

Let $\Omega \subset \RR^2$ be a fixed region of land, and suppose $b \in C^1(\Omega)$ is a given bedrock elevation.  Let $a \in L^2(\Omega)$ denote a given “surface mass balance” function, the annually-averaged rate of ice accumulation (snow) minus melt and runoff.  Then the surface elevation $s\in C(\Omega)$ of a steady-state, isothermal, and non-sliding ice sheet (glacier) under the shallow ice approximation (SIA) [23] satisfies the following VI problem: find $s \in \mathcal{K} = \{s\ge b\}$ such that

\begin{equation}
\int_\Omega \Gamma (s-b)^{n+2} |\grad s|^{n-1} \grad s \cdot \grad (v-s) \dx \ge \int_\Omega a (v-s)\dx \quad \text{ for all } v \in \mathcal{K}, \label{eq:siavi}
\end{equation}

where $n>1$ and $\Gamma>0$ are physical constants for ice; $n=3$ is a typical value [23].  The doubly-nonlinear and doubly-degenerate operator here is

\begin{equation}
\ip{f(s)}{v} = \int_\Omega \Gamma (s-b)^{n+2} |\grad s|^{n-1} \grad s \cdot \grad v\dx. \label{eq:siaoperator}
\end{equation}

If the bed is flat ($b=0$) then a power transformation converts (<ref>) into form (<ref>) with $p=n+1$.  Less is understood for general beds $b$, but existence holds for the VI [25], and a power of the solution lives in a Sobolev space: $(s-b)^{2p/(p-1)} \in W^{1,p}(\Omega)$.  This theory, and also observations of ice sheets, shows that $|\grad s|$ is usually unbounded as one approaches the free boundary from the icy ($s>b$) side.

The numerical results in Section <ref> will apply FASCD Algorithms <ref> and <ref> to the above Examples.

§ CONSTRAINT DECOMPOSITION (CD)

Suppose that there are $m<\infty$ closed linear subspaces $\cV_i \subset \cV$, so that the sum

\begin{equation}
\cV = \sum_{i=0}^{m-1} \cV_i \label{eq:subspacedecomp}
\end{equation}

holds in the sense that if $w \in \cV$ then there exist $w_i \in \cV_i$ so that $w = \sum_i w_i$.  We say (<ref>) is a space decomposition of $\cV$ [36].  For a closed, convex subset $\cK \subset \cV$, suppose further that $\cK_i \subset \cV_i$ are nonempty, closed, and convex subsets such that

\begin{equation}
\cK = \sum_{i=0}^{m-1} \cK_i. \label{eq:constraintdecomp}
\end{equation}

The sum in (<ref>) must hold in two senses [33]: (i) if $w \in \cK$ then there exist $w_i \in \cK_i$ so that $w = \sum_i w_i$, and (ii) if $z_i \in \cK_i$ for each $i$ then $\sum_i z_i \in \cK$.  Note that neither decomposition (<ref>) or (<ref>) is required to be unique, and also notice that sense (ii) is automatic for (<ref>) because the $\cV_i$ are subspaces.  Finally, suppose there are bounded, generally nonlinear, decomposition operators $\Pi_i : \cK \to \cK_i$ such that if $v \in \cK$ then

\begin{equation}
v = \sum_{i=0}^{m-1} \Pi_i v,  \label{eq:constraintrestrictionsum}
\end{equation}

with the decomposition operators satisfying the stability property

\begin{equation} \label{eq:decompositionstability}
\left(\sum_{i=0}^{m-1} \| \Pi_i u - \Pi_i v \|^2\right)^{1/2} \lesssim \|u - v\|.
\end{equation}

Clearly (<ref>) implies sense (i) for (<ref>).  A constraint decomposition (CD) of $\cK$ is a choice of $\cV_i,\cK_i,\Pi_i$ satisfying (<ref>)–(<ref>) [32].  Figure <ref> suggests how a low-dimensional example might look; note $\cK_i \not\subset \cK$ in this and many other cases.

A constraint decomposition (CD) for a unilateral obstacle problem on a two-point space $\Omega=\{x_1,x_2\}$, with $\mathcal{V}=\left\{v \,:\, \Omega \to \RR\right\}$ and $\mathcal{K}=\{v\ge \psi\}$.

We will construct a multilevel FE CD discretization in Sections <ref>–<ref>.  However, the CD concept can be applied at the level of the continuum problem.  For example, in a unilateral obstacle problem ($\cK = \{v \ge \psi\}$) one can construct an overlapping domain decomposition of $\cV$ by using a smooth, finite partition of unity $\{\varphi_i\}$, subspaces $\cV_i = \{w \in \cV : w = 0 \text{ if } \varphi_i = 0\}$, constraint sets $\cK_i = \{v \in \cV_i : v \ge \varphi_i \psi\}$, and decomposition operators $\Pi_i(v) = \varphi_i v$.  One may also construct a CD which is a disjoint frequency decomposition of a Hilbert space of periodic functions on a cube $\Omega=(0,1)^d$.  The multilevel CD, over a finite-dimensional FE discretization, will, in a vague sense, approximate such a frequency decomposition.

Associated to any CD are certain iterative methods [32, 36].  The following multiplicative/successive cd-mult algorithm, which generalizes the Gauss–Seidel iteration, replaces $u \in \cK$ with a new iterate $w\in\cK$ which should be closer to $u^\star \in \cK$ solving (<ref>).  The idea is to solve smaller VI problems over each subset $\cK_i$.  Each subset solution corrects the global iterate immediately, so the ordering of the $\cK_i$ is important.

for $i = 0,\dots,m-1$:
find $w_i\in \cK_i$ such that
$\displaystyle \Big<f\Big(\sum_{j<i} w_j + w_i + \sum_{j>i} \Pi_j u\Big),\, v_i - w_i\Big> \ge \ip{\ell}{v_i - w_i} \textnormal{ for all } v_i \in \cK_i$
return $w=\sum_i w_i\in\cK$

The sum “$\sum_{j>i} \Pi_j u$” should be read as “those parts of $u$ which need improvement”.  Tai [32] also proposes an additive/parallel version of the algorithm, but it is not used here.

The VI problems in cd-mult reduce to constrained line searches when the $\cV_i$ are one-dimensional.  If $e_i = w_i - \Pi_i u \in \cV_i$ denotes the correction from $\cK_i$, then the new iterate $w = u + \sum_i e_i$ is computed by adding a correction from each subspace $\cV_i$, but so that $w \in \cK$; the corrections preserve admissibility.  Furthermore, if VI problem (<ref>) corresponds to convex optimization (<ref>) then cd-mult generates monotonically-decreasing objective values.

Suppose $f=F'$ for $F:\cK\to\RR$ convex and differentiable, and let $u\in\cK$.  The output $w \in \cK$ from cd-mult satisfies

\begin{equation}
F(w) - \ip{\ell}{w} \le F(u) - \ip{\ell}{u}.  \label{eq:objectivemonotone}
\end{equation}

cd-mult is meaningful whether or not we are in the optimization case $f=F'$, and even if $f$ is nonlinear, non-local, or defined only on $\cK$, but its practical and efficient FE implementation for such general operators $f$ seems not to have been addressed in the literature.  In particular, references [22, 32] only apply the multilevel CD method to the classical obstacle problem, wherein $f=F'$ is linear, local, and defined on all of $\mathcal{V}$.  In Section <ref> we will extend multilevel CD iterations to general nonlinear operators by applying the full approximation scheme (FAS) idea of Brandt [5], and then illustrate its application to general operators in Examples <ref>–<ref>.

From now on we restrict to the case of box constraints [3, 21], equivalently bilateral obstacle problems, over $\mathcal{V}=W^{1,p}(\Omega)$.  The lower obstacle $\underline{\gamma}$ and upper obstacle $\overline{\gamma}$ are given measurable functions, defined on the closure $\overline{\Omega}$, with values in the extended real line $\tilde\RR = [-\infty,+\infty]$, satisfying the following properties: $\underline{\gamma} < +\infty$, $\overline{\gamma}>-\infty$, and $\underline{\gamma} \le \overline{\gamma}$.  Also suppose $\partial\Omega$ is split into disjoint sets $\partial\Omega = \partial_D \Omega \cup \partial_N \Omega$, with $\partial_D \Omega$ of positive Hausdorff measure, and $\partial_N \Omega$ possibly empty.  Dirichlet data is given by $g_D:\partial_D \Omega \to \RR$, and we assume $\underline{\gamma} \le g_D \le \overline{\gamma}$ on $\partial_D \Omega$.  Now define

\begin{equation}
\cK = \left\{w\,:\,\underline{\gamma} \le w \le \overline{\gamma} \, \text{ and }\, w\big|_{\partial_D \Omega} = g_D\right\} \subset \cV =W^{1,p}(\Omega), \label{eq:originalconstraintset}
\end{equation}

the closed and convex constraint set; this treats inequality constraints and Dirichlet conditions in a unified manner.  Dirichlet values are in the trace sense [17], and $g_D$ is assumed to be as regular as needed.  Any Neumann conditions over $\partial_N \Omega$ are imposed weakly by the operator $f$, via an integral over $\partial_N\Omega$.  We will numerically solve VI problems (<ref>) over constraint sets (<ref>).

§ MULTILEVEL FINITE ELEMENTS

We seek to rapidly compute the approximate solution of box-constrained VI problems (<ref>) and (<ref>) using finite elements (FE).  In this Section a nested multilevel $P_1$ or $Q_1$ element [16] approximation to (<ref>) will be constructed by refinement of a coarse mesh.  The mesh-refinement hierarchy permits the construction of a multilevel CD of the finest-level constraint set.

The mesh hierarchy and FE spaces can be constructed by standard nested refinement.  Let $\Omega \subset \RR^d$ be an open bounded polygon and suppose $\mathcal{T}^0$ is the coarsest mesh, a finite set of non-overlapping, nondegenerate cells with union equal to $\overline{\Omega}$ and maximum cell diameter (mesh size) $h_0>0$.  We assume $\mathcal{T}^0$ is conforming, with no hanging nodes.  For $j=1,\dots,J$, let $\mathcal{T}^j$ be the uniform refinement of $\mathcal{T}^{j-1}$; each $\mathcal{T}^{j-1}$ cell becomes $2^d$ cells in $\mathcal{T}^j$ with mesh size $h_j = 2^{-j} h_0$.  We call $\mathcal{T}^J$ the finest-level mesh and $J+1$ the number of levels.  Let $m_j$ be the number of vertices in $\mathcal{T}^j$, the degrees of freedom at the $j$th level.

Let $\mathcal{V}^j$ be nested FE spaces over each mesh $\mathcal{T}^j$, with arbitrary boundary values:

\begin{equation}
\mathcal{V}^0 \subset \mathcal{V}^1 \subset \dots \subset \mathcal{V}^J \subset \mathcal{V}=W^{1,p}(\Omega).  \label{eq:fe:nestedspaces}
\end{equation}

As is common in the discretization of variational inequalities, we employ lowest-order continuous Lagrange elements $P_1$ or $Q_1$ on simplices and tensor-product cells respectively[Keith & Surowiec have recently proposed a scheme for solving VIs using high-order finite elements [27].]. This choice is convenient because it enjoys nodal monotonicity:

\begin{equation}
\bw \ge \bz \quad \implies \quad w \ge z \label{eq:nodalmonotonicity}
\end{equation}

for all $w,z \in \mathcal{V}^j$, where $\bw,\bz \in \RR^{m_j}$ denote vectors of nodal values. In other words, to enforce an inequality constraint, with $P_1$ or $Q_1$ we only enforce it pointwise on the degrees of freedom. For $k\ge 2$ there are functions in the standard bases of $P_k$ and $Q_k$ which take on negative values <cit.>, violating (<ref>).

To construct box constraints on all meshes in the hierarchy, let $\tilde{\mathcal{V}}^j$ denote the set of functions on the vertices of $\mathcal{T}^j$ with values in $\tilde{\RR} = [-\infty,+\infty]$.  The finest-level lower and upper obstacles $\underline{\gamma}^J, \overline{\gamma}^J \in \tilde{\mathcal{V}}^J$ are constructed by suitable interpolation of $\underline{\gamma}$, $\overline{\gamma}$, e.g. by applying the Clément interpolation operator [12].  Again the obstacles must satisfy $\underline{\gamma}^J < +\infty$, $-\infty < \overline{\gamma}^J$, and $\underline{\gamma}^J \le \overline{\gamma}^J$.  On the finest-level mesh the discrete boundary data $g_D^J \in \mathcal{V}^J$ is assumed to satisfy $g_D^J = g_D$ on $\partial_D \Omega$—this restricts $g_D$ in (<ref>) to be piecewise-linear—and $\underline{\gamma}^J \le g_D^J \le \overline{\gamma}^J$ on those vertices of $\mathcal{T}^J$ which lie in $\partial_D \Omega$.  Then we define

\begin{equation}
\mathcal{K}^J = \big\{w\,:\,\underline{\gamma}^J \le w \le \overline{\gamma}^J \text{ on vertices of } \mathcal{T}^J, \, \text{ and } \, w|_{\partial_D\Omega} = g_D^J|_{\partial_D\Omega}\big\} \subset \mathcal{V}^J. \label{eq:fe:fineconstraintset}
\end{equation}

This nonempty, closed, and convex subset is the finest-level constraint set.

Let $f^J:\mathcal{K}^J \to (\mathcal{V}^J)'$ be a discretization of $f$ in (<ref>).  Note that the functional $\ell^J=\ell$ is defined and continuous over $\mathcal{V}^J$.  We seek the FE solution $u^J \in \mathcal{K}^J$ of the finest-level, finite-dimensional, nonlinear VI

\begin{equation}
\ip{f^J(u^J)}{v-u^J} \ge \ip{\ell^J}{v-u^J} \qquad \text{for all } v\in \cK^J. \label{eq:fe:vi}
\end{equation}

If $f^J$ is continuous and $p$-coercive on $\mathcal{K}^J$, in the sense that (<ref>) holds for all $u,v \in \mathcal{K}^J$, then (<ref>) is well-posed.

While (<ref>) approximates (<ref>), we warn the reader that it is usually not conforming.  That is, even if quadrature and other implementation details for $f^J$ are ignored, and despite the vector space inclusion $\mathcal{V}^J \subset \mathcal{V}$, from definitions (<ref>) and (<ref>) we do not expect $\mathcal{K}^J \subset \mathcal{K}$ in general.  Even for continuous obstacles $\underline{\gamma}, \overline{\gamma}$, admissibility with respect to the FE obstacles $\underline{\gamma}^J, \overline{\gamma}^J$ will generally not imply $\mathcal{K}$-admissibility.  For a practical example, nodal elevation data in a glacier application misses the topographic detail in-between the vertices [9].

The multilevel algorithm defined in Section <ref> requires five transfer operators between meshes.  Of these, three are employed in standard FAS [34]; Table <ref> will help the reader keep track.  Canonical prolongation $P:\mathcal{V}^{j-1}\to\mathcal{V}^j$ implements the vector space inclusion.  Canonical (dual) restriction $R:(\mathcal{V}^j)'\to(\mathcal{V}^{j-1})'$ acts as an identity in the sense that if $\sigma \in (\mathcal{V}^j)'$ then $\ip{R\sigma}{z} = \ip{\sigma}{z}$ for all $z \in \mathcal{V}^{j-1} \subset \mathcal{V}^j$.  These canonical transfers generally involve a change of FE representation, and in particular $R$ is “full-weighting” [34].  To define injection $\iR:\mathcal{V}^j\to\mathcal{V}^{j-1}$, note that the vertices of $\mathcal{T}^{j-1}$ are vertices of $\mathcal{T}^j$ and so for $z\in\mathcal{V}^j$ define $\iR z$ as the element of $\mathcal{V}^{j-1}$ with the same (coarse) nodal values.  Injection is monotonic in the sense that if $w \le z$ then $\iR w \le \iR z$.

name      mapping      linear?
canonical prolongation            $P:\mathcal{V}^{j-1}\to\mathcal{V}^j$

canonical (dual) restriction      $R:(\mathcal{V}^j)'\to(\mathcal{V}^{j-1})'$

(nodal) injection                 $\iR:\mathcal{V}^j\to\mathcal{V}^{j-1}$

maximum injection               $\maxR:\tilde{\mathcal{V}}^j\to\tilde{\mathcal{V}}^{j-1}$     55

minimum injection               $\minR:\tilde{\mathcal{V}}^j\to\tilde{\mathcal{V}}^{j-1}$     55

Transfer operators used in this paper.

The monotone injections $\maxR,\minR$ [32, 22], which will only be applied to constraints in $\tilde{\mathcal{V}}^j$, are nonlinear maps which maximize and minimize, respectively, nodal values over the open supports of coarser-mesh basis functions (i.e. over the stars of the coarse vertices).  To make the construction precise, let $x_q^j$ denote the vertices of mesh $\mathcal{T}^j$, with corresponding basis functions $\psi_q^j$.  For an $\tilde\RR$-valued nodal function $z\in\tilde{\mathcal{V}}^j$ we define functions in $\tilde{\mathcal{V}}^{j-1}$:

\begin{align}
(\maxR z)(x_p^{j-1}) &= \max \{z(x_q^j) \,:\, \psi_p^{j-1}(x_q^j) > 0\}, \label{eq:fe:monotoneinjections} \\
(\minR z)(x_p^{j-1}) &= \min \{z(x_q^j) \,:\, \psi_p^{j-1}(x_q^j) > 0\}. \notag
\end{align}

These operators are constructed to enforce the property

\begin{equation}
\minR z \le z \le \maxR z,  \label{eq:fe:monotoneproperty}
\end{equation}

which does not hold for the ordinary injection operator $\iR$. An illustration of the construction on nested meshes is given in Figure <ref>; implementing these operators on non-nested meshes could be efficiently done using a supermesh [20]. In the next Section we will use $\maxR,\minR$ to construct “level defect constraints” such that the zero function is always admissible.  This construction depends on additional ordering properties which follow from definition (<ref>):

\begin{equation}
z\ge 0 \implies \minR z \ge 0 \quad \text{ and } \quad z \le 0 \implies \maxR z \le 0. \label{eq:fe:monotoneadditional}
\end{equation}

The monotone injection operators $\maxR$ and $\minR$ assign to a coarse degree of freedom (blue) the maximum and minimum respectively of the fine function over the patch of coarse cells sharing its vertex. The coarse mesh is in black and the fine in grey.

§ MULTILEVEL CDS FROM LEVEL DEFECT CONSTRAINTS

A key choice in devising a multilevel scheme for VIs is how to construct bounds for the coarser levels.  Our multilevel CD scheme uses defect constraints [22] derived from an admissible finest-level iterate.  The next Example illustrates why we do not[An exception occurs in FMG; see Section <ref>.] apply injection operators directly to the finest-level box constraints $\underline{\gamma}^J,\overline{\gamma}^J$.

Let $\Omega = (0,L) \subset \RR^1$ and consider a two-level mesh hierarchy ($J=1$).  Let $\underline{\gamma}^1$ be a sawtooth function with minimum value zero, and suppose $\overline{\gamma}^1=C+\underline{\gamma}^1$ with $C>0$ strictly less than the maximum of $\underline{\gamma}^1$.  Figure <ref> shows such functions; here levels $\mathcal{T}^1$ and $\mathcal{T}^0$ have $5$ and $3$ vertices, respectively.  For appropriate boundary conditions the fine-level constraint set $\mathcal{K}^1$ defined by (<ref>) is nonempty.  However, various schemes to coarsen the constraints by direct transfers of $\underline{\gamma}^1,\overline{\gamma}^1$ are doomed to failure.  Suppose we propose

\begin{equation}
\label{eq:fe:badstandard}
\underline{\gamma}^0, \overline{\gamma}^0 \equiv \iR \underline{\gamma}^1, \iR \overline{\gamma}^1. \qquad \text{\emph{[problematic]}}
\end{equation}

Then $\underline{\gamma}^0=0$ and $\overline{\gamma}^0=C$ identically.  In that case the corresponding constraint set $\mathcal{K}^0$ is nonempty, but prolongation $Pw^0$ of an admissible solution of the coarse problem ($w^0\in\mathcal{K}^0$; $\underline{\gamma}^0 \le w^0 \le \overline{\gamma}^0$) is not admissible on the finer level ($Pw^0 \notin \mathcal{K}^1$).  Attempting to fix this by projecting (truncating) into $\mathcal{K}^1$ would introduce high frequencies not present in $w^0$.  On the other hand, if

\begin{equation}
\label{eq:fe:badmonotone}
\underline{\gamma}^0, \overline{\gamma}^0 \equiv \maxR \underline{\gamma}^1, \minR \overline{\gamma}^1 \qquad \text{\emph{[problematic]}}
\end{equation}

then the coarse-level constraint set $\mathcal{K}^0$ is empty because $\underline{\gamma}^0 > \overline{\gamma}^0$ (Figure <ref>, right).

For the finer-level constraints $\underline{\gamma}^1$ (blue) and $\overline{\gamma}^1$ (red) on the left, direct application of $\iR,\minR,\maxR$ to generate coarser-level constraints is problematic.  Possibility (<ref>), which generates an empty coarser-level set $\mathcal{K}^0$, is at right.  Defect constraints (<ref>) avoid this difficulty.

Iterate-dependent defect constraints will satisfy ordering (<ref>) below and bypass the difficulties in the above Example.  Suppose $w^J \in \cK^J$ is an admissible iterate for problem (<ref>).  We define the finest-level defect constraints [22] as functions in $\tilde{\mathcal{V}}^J$:

\begin{equation}
\underline{\chi}^J = \underline{\gamma}^J - w^J, \qquad \overline{\chi}^J = \overline{\gamma}^J - w^J. \label{eq:fe:defectconstraints}
\end{equation}

(If $\underline{\gamma}^J=-\infty$ then $\underline{\chi}^J=-\infty$ by definition, and similarly if $\overline{\gamma}^J=+\infty$ then $\overline{\chi}^J=+\infty$.)  Definition (<ref>) implies that $\underline{\chi}^J \le 0 \le \overline{\chi}^J$ everywhere in $\overline{\Omega}$, including on the Dirichlet boundary $\partial_D\Omega$.    Observe that a corrected iterate $w^J + z^J$ is admissible, $w^J + z^J \in \cK^J$, if and only if the perturbation $z^J \in \mathcal{V}^J$ is zero on $\partial_D\Omega$ and bounded by the defect constraints, $\underline{\chi}^J \le z^J \le \overline{\chi}^J$.

Now descending to coarser levels, $j=J-1,\dots,1,0$, we apply the monotone injections to generate level defect constraints (LDCs) in $\tilde{\mathcal{V}}^j$:[The definition of $\underline{\chi}^j$ matches that in Section 4 of [22]; see (4.22) and Figure 4.1 therein.]

\begin{equation}
\underline{\chi}^{j} = \maxR \underline{\chi}^{j+1}, \qquad \overline{\chi}^{j} = \minR \overline{\chi}^{j+1}. \label{eq:fe:chilevels}
\end{equation}

Observe that by (<ref>) and (<ref>),

\begin{equation}
\underline{\chi}^{J} \le \dots \le \underline{\chi}^0 \le 0 \le \overline{\chi}^0 \le \dots \le \overline{\chi}^J. \label{eq:fe:chiordering}
\end{equation}

The LDC $\underline{\chi}^j$ is the most-negative function from $\tilde{\mathcal{V}}^j$ such that if another element of $\tilde{\mathcal{V}}^j$ is above $\underline{\chi}^j$ then it is also above $\underline{\chi}^{j+1}$; similar comments apply to $\overline{\chi}^{j}$.  Note that for Algorithm <ref> below, any system satisfying (<ref>) will suffice; the particular formulas in (<ref>) are not required.

This construction of multilevel box constraints, likewise the unilateral construction in [22], permits large, yet admissible, coarse corrections, which promotes multilevel solver efficiency.  When coarsening the problem downward in a V-cycle, one cannot predict the upcoming coarser corrections and so one must construct constraint sets from which any sum will be admissible.  With this in mind, and looking forward to telescoping sums (<ref>), we compute the LDC differences [22]:

\begin{equation}
\underline{\phi}^j = \underline{\chi}^j - \underline{\chi}^{j-1}, \qquad \overline{\phi}^j = \overline{\chi}^j - \overline{\chi}^{j-1}.  \label{eq:fe:philevels}
\end{equation}

(Take $\underline{\chi}^{-1}=\overline{\chi}^{-1}=0$ so that $\underline{\phi}^0=\underline{\chi}^0$ and $\overline{\phi}^0=\overline{\chi}^0$.  By definition, if $\underline{\chi}^j(x_p)=-\infty$ then we set $\underline{\phi}^j(x_p)=-\infty$ regardless of the value of $\underline{\chi}^{j-1}(x_p)$, and etc.)  Note that while $\underline{\phi}^{j},\overline{\phi}^{j} \in \tilde{\mathcal{V}}^J$ again bracket zero, they are not ordered as in (<ref>).

The corresponding downward and upward constraint sets, respectively, are closed and convex subsets of the FE spaces $\mathcal{V}^j$:

\begin{align}
\mathcal{D}^j = \left\{v \in \mathcal{V}^j \,:\, \underline{\phi}^j \le v \le \overline{\phi}^j \text{ and } \, v|_{\partial_D\Omega} = 0\right\}, \label{eq:fe:constraintsets} \\
\mathcal{U}^j = \left\{v \in \mathcal{V}^j \,:\, \underline{\chi}^j \le v \le \overline{\chi}^j \text{ and } \, v|_{\partial_D\Omega} = 0\right\}. \notag
\end{align}

Because $\underline{\chi}^j \le \underline{\phi}^j \le 0 \le \overline{\phi}^j \le \overline{\chi}^j$, it follows that $\mathcal{D}^j \subseteq \mathcal{U}^j$.  Note that the finest-level set $\mathcal{U}^J$ is equivalent to the original constraint set $\mathcal{K}^J$: for all $z^J \in \mathcal{V}^J$ , $z^J \in \mathcal{U}^J$ if and only if $w^J+z^J \in \mathcal{K}^J$.

The V-cycle solver in the next Section (Algorithm <ref>) will compute corrections from $\mathcal{D}^j$ and $\mathcal{U}^j$.  For general box constraints the upward direction is, however, based on incomplete CDs as we now explain.  Each upward set $\mathcal{U}^j$ is generally only a superset of a sum of downward sets $\mathcal{D}^j$, though this becomes a true CD satisfying (<ref>)–(<ref>) in unilateral cases.  The following Lemma, apparently new, extends the unilateral CD construction underlying Algorithm 4.2 in [22]; see also equation (59) in [32].

\begin{equation}
\mathcal{U}^j \supseteq \mathcal{D}^j + \mathcal{D}^{j-1} + \dots + \mathcal{D}^0 \label{eq:fe:downwardsuminclusion}
\end{equation}

for $j=0,1,\dots,J$.  In unilateral cases, where either $\underline{\gamma}^J=-\infty$ identically or $\overline{\gamma}^J=+\infty$ identically, inclusion (<ref>) becomes a full CD with $\mathcal{U}^j=\sum_{i=0}^j \mathcal{D}^i$.

From nesting (<ref>), subspace decomposition (<ref>) holds, and definition (<ref>) shows $\mathcal{D}^i \subset \mathcal{U}^i \subset \cV^i$.  Furthermore, telescoping sums hold for any $j=0,1,\dots,J$:

\begin{equation}
\sum_{i=0}^j \underline{\phi}^i = \underline{\chi}^j, \qquad \sum_{i=0}^j \overline{\phi}^i = \overline{\chi}^j.  \label{eq:fe:telescoping}
\end{equation}

Thus if $y^i \in \mathcal{D}^i$ for $0 \le i \le j$ then

\begin{equation}
\underline{\chi}^j = \sum_{i=0}^j \underline{\phi}^i \le \sum_{i=0}^j y^i \le \sum_{i=0}^j \overline{\phi}^i \le \overline{\chi}^j, \label{eq:fe:lemmaordering}
\end{equation}

so (<ref>) holds for any box constraints.

In unilateral cases we can also construct decomposition operators $\Pi_i$ and thereby show (<ref>).  For concreteness suppose $\overline{\gamma}^J=+\infty$, thus $\overline{\chi}^j=+\infty$ and $\overline{\phi}^i = +\infty$ for all levels; the $\underline{\gamma}^J=-\infty$ case is similar.  For $v\in \mathcal{V}^j$ and $0\le i \le j$, let $I_{j\to i}^\ominus\,:\,\tilde{\cV}^j \to \tilde{\cV}^i$ be the operator applying the minimum (monotone) injection $j-i$ times: $I_{j\to i}^\ominus v = \minR(\cdots(\minR v)\cdots) $.  (Set $I_{j\to j}^\ominus v = v$ and $I_{j\to -1}^\ominus=0$ by definition.)  Now define nonlinear decomposition operators $\Pi_i:\mathcal{U}^j \to \mathcal{D}^i$:

\begin{equation}
\Pi_i z^j \coloneqq I_{j\to i}^\ominus(z^j - \underline{\chi}^j) - I_{j\to i-1}^\ominus(z^j - \underline{\chi}^j) + \underline{\phi}^i.  \label{eq:fe:unilateraldecompositionoperator}
\end{equation}

(Compare equation (4.9) in [22].)  Property (<ref>) implies $I_{j\to i}^\ominus(z^j - \underline{\chi}^j) \ge I_{j\to i-1}^\ominus(z^j - \underline{\chi}^j)$, thus $\Pi_i z^j \ge \underline{\phi}^i$, so $\Pi_i z^j \in \mathcal{D}^i$, and furthermore stability (<ref>) holds <cit.>.  On the other hand, the following sum telescopes, leaving only its first term plus the sum in (<ref>):

\begin{align*}
\sum_{i=0}^j \Pi_i z^j &= z^j - \underline{\chi}^j + \sum_{i=0}^j \underline{\phi}^i = z^j.
\end{align*}

This shows (<ref>), thus (<ref>), and that (<ref>) is equality.

One might attempt to convert (<ref>) into a full CD for arbitrary box constraints by constructing different decomposition maps $\Pi_i$, for example by splitting $z^j = \min\{z^j,0\} + \max\{z^j,0\}$ and applying unilateral formulas (<ref>) to $\min\{z^j,0\}$ and $\max\{z^j,0\}$ separately.  The above proof shows that the lower bound $\Pi_i (\min\{z^j,0\}) \ge \underline{\phi}^i$ indeed holds.  However, the following Example shows why $\Pi_i(\min\{z^j,0\})$ may have an arbitrarily large maximum, and thus exceed any finite upper LDC $\overline{\phi}^i$.  In fact, $\Pi_i$ defined above does not generally map all non-positive elements of $\mathcal{U}^J$ into $\mathcal{D}^i$ if a finite upper LDC is present.  The Example does not exclude other strategies for showing full CDs in the general box-constrained case.

Let $a > 0$ be an arbitrary positive number.  Consider a 2-level hierarchy ($J=1$), and suppose the lower defect constraint is the constant function $\underline{\chi}^1=-a$.  Note $\underline{\chi}^0=\maxR \underline{\chi}^1=-a$ also, so $\underline{\phi}^1=0$ and $\underline{\phi}^0=-a$.  Suppose $z^1$ is a “sawtooth” function which takes on value $-a$ on the coarse vertices and value $0$ on the other vertices.  Then, for any $\overline{\chi}^1\ge 0$, we have $\overline{\chi}^1 \ge 0 \ge z^1\ge \underline{\chi}^1$, so $z^1 \in \mathcal{U}^1$.  However, $I_{1\to 0}^\ominus(z^1 - \underline{\chi}^1) = \minR(z^1 - \underline{\chi}^1) = 0$ identically.  From (<ref>), $\Pi_1 z^1$ simplifies to $\Pi_1 z^1 = z^1 - \underline{\chi}^1$.  But then $\Pi_1 z^1$ attains a maximum $a\ge 0$.  This maximum can be chosen to exceed any upper defect obstacle $\overline{\phi}^1$, and $\Pi_1 z^1 \notin \mathcal{D}^1$ would then occur.

Each upward set $\mathcal{U}^j$ can also be decomposed, incompletely in the sense of (<ref>), down to an arbitrary coarser level.  The next Lemma, which does not follow from Lemma <ref>, justifies upward admissibility in Algorithm <ref>.

For any $j=0,1,\dots,J$ and $0\le k\le j$,

\begin{equation}
\mathcal{U}^j \supseteq \mathcal{D}^j + \mathcal{D}^{j-1} + \dots + \mathcal{D}^{k+1} + \mathcal{U}^k \label{eq:fe:upwardsuminclusion}
\end{equation}

In unilateral cases, where either $\underline{\gamma}^J=-\infty$ or $\overline{\gamma}^J=+\infty$, actually $\mathcal{U}^j = \mathcal{U}^k + \sum_{i=k+1}^j \mathcal{D}^i$ is a full CD.

Definition (<ref>) shows inclusion (<ref>).  In unilateral cases we may follow the proof of Lemma <ref> with (<ref>) unchanged for $i=k+1,\dots,j$, but modified in the $i=k$ case: $\Pi_k z^j = I_{j\to k}^\ominus(z^j - \underline{\chi}^j) + \underline{\chi}^k$.  Otherwise the proof goes through as before.

Inclusions (<ref>) and (<ref>) show the admissibility of all iterates and corrections in Algorithm <ref> below.  These inclusions are special to multilevel schemes based on defect constraints, and we have seen that they hold regardless of whether decomposition operators $\Pi_i$ can be constructed.  In fact we can propose a new, incomplete variant of the cd-mult iteration.  For $w^J\in\mathcal{K}^J$, and existing downward perturbations $y^i\in\mathcal{D}^i$, the iteration applies in the special case where inclusion (<ref>) holds:

icd-tele(w^J, {y^i}_i=1^J):
for $j = 0,\dots,J$:
find $z^j\in \mathcal{U}^j$ so that for all $v \in \mathcal{U}^j$,
$\displaystyle \Big<f\Big(w^J + \sum_{i>j} y^i + z^j\Big),v-z^j\Big> \ge \ip{\ell}{v-z^j}$
return $\tilde w^J=w^J + z^J$

Note that the downward perturbations $y^i$, which apparently disappear as icd-tele proceeds, can actually play a practical role in solving each VI problem, namely when constructing an initial iterate for $z^j$; see the next Section.

Algorithm <ref> below incorporates icd-tele for up-smoothing in the larger constraint sets $\mathcal{U}^j$.  The observation that this can be more efficient, illustrated by Example <ref> in particular, seems to be new.  For comparison, $V(1,1)$ cycles in [22] use half-size versions of the smaller sets $\mathcal{D}^j$ both on descent and ascent; see also <cit.>.  However, any convergence analysis extending the arguments of [22, 32] might require full CDs.

§ FASCD V-CYCLE

The V-cycle in this Section extends the multilevel CD schemes of Tai [32] and Gräser & Kornhuber <cit.> by following a nonlinear full approximation scheme (FAS) approach [7]; it reduces to FAS for PDEs when inequality constraints are removed.  Each V-cycle generates a correction to a finest-level iterate $w^J \in \mathcal{K}^J$, via information from each coarser FE subspace $\mathcal{V}^j$, by approximately solving VI problems in constraint sets $\mathcal{D}^j$ and $\mathcal{U}^j$ during the downward and upward stages, respectively.  The corrections are inductively constructed so that iterates remain admissible at every level.  A key idea is that up-smoothing solves less-restricted VI problems than down-smoothing.  Going downward one does not know what coarse corrections are coming, thus one must choose smaller constraint sets ($\mathcal{D}^j$) such that the correction sum is admissible once we return to a given level.  Going upward, however, no further coarse corrections can violate admissibility, so we relax in the larger sets $\mathcal{U}^j$.

Suppose $y^i \in \mathcal{D}^i$, for $i=J,\dots,j+1$, are already-computed corrections during the downward part of the V-cycle (Figure <ref>, left).  By (<ref>) the correction at this point is admissible, namely $y^J + \dots + y^{j+1} \in \mathcal{U}^J$, equivalently $w^J + y^J + \dots + y^{j+1} \in \mathcal{K}^J$.  Relaxation then occurs at the next-coarser level, yielding the next correction $y^j \in \mathcal{D}^j$.  Solving the coarsest-level VI results in the full downward correction $y^J + \dots + y^1 + z^0 \in \mathcal{U}^J$.  Starting upward, $y^1 + z^0 \in \mathcal{U}^1$ is an initial iterate for relaxation on level $j=1$, which then yields $z^1 \in \mathcal{U}^1$ and an admissible correction $y^J + \dots + y^2 + z^1 \in \mathcal{U}^J$, using inclusion (<ref>).  Continuing in this “telescoping” way, up-smoothing in the constraint sets $\mathcal{U}^j$, starting from initial iterate $y^j+z^{j-1}$, finally generates $z^J\in \mathcal{U}^J$ as the V-cycle correction.  The new finest-level iterate is $\tilde{w}^J = w^J + z^J \in \mathcal{K}^J$.

Left: The FASCD V-cycle Algorithm <ref> computes downward corrections $y_j \in \mathcal{D}^j$, but the upward corrections $z_j\in\mathcal{U}^j$ are in larger sets. Right: FMG Algorithm <ref> generates the initial iterate on finer levels by injection and truncation (doubled edges).

We must specify the problems solved at each mesh level.  As in FAS multigrid for nonlinear PDEs [7, 8, 34], both the solution approximation and the residual must be coarsened, and a source functional created.  On mesh level $j$ we assume a re-discretized nonlinear operator $f^j$, with output in $(\cV^j)'$, which is an FE discretization of $f$ in (<ref>).  At the start of each V-cycle we assume that an admissible iterate $w^J \in \mathcal{K}^J$ has been used to generate the LDCs and constraint sets $\mathcal{D}^j$, $\mathcal{U}^j$ at each level (Section <ref>); they do not depend on the coarse-level corrections.

As we descend we define the initial $j$th-level iterates

\begin{equation}
w^j = \iR(w^{j+1} + y^{j+1}),  \label{eq:fe:definew}
\end{equation}

from which we define $\ell^j \in (\cV^j)'$:

\begin{equation}
\ell^j = \begin{cases} \ell^J, & j=J \\
f^j(w^j) + R\left(\ell^{j+1}-f^{j+1}(w^{j+1}+y^{j+1})\right), & j<J. \end{cases} \label{eq:fe:levelsource}
\end{equation}

As is typical of FAS <cit.>, injection $\iR$ coarsens the (primal) iterates, but canonical restriction $R$ coarsens (dual) linear functionals.  Descending through the levels $j=J,J-1,\dots,1$ we solve the following VI problems for $y^j \in \mathcal{D}^j$:

\begin{equation}
\ip{f^j(w^j + y^j)}{v-y^j} \ge \ip{\ell^j}{v-y^j} \qquad \text{for all } v\in \mathcal{D}^j. \label{eq:fe:downvi}
\end{equation}

The initial iterate is $y_{(0)}^j=0$; this is admissible because $\underline{\phi}^j \le 0 \le \overline{\phi}^j$.

On ascent, $j=0,1,\dots,J$, the same VI is solved but for $z^j \in \mathcal{U}^j$:

\begin{equation}
\ip{f^j(w^j + z^j)}{v-z^j} \ge \ip{\ell^j}{v-z^j} \qquad \text{for all } v\in \mathcal{U}^j. \label{eq:fe:upvi}
\end{equation}

However, the nontrivial initial iterate for this problem is found by prolonging $z^{j-1} \in \mathcal{U}^{j-1}$ and then updating the current (i.e. down-smoothed) correction:

\begin{equation}
z_{(0)}^j = y^j + P z^{j-1}.  \label{eq:fe:upwardinitial}
\end{equation}

An exception is that $z_{(0)}^0=0$ on the coarsest level.  Note $z_{(0)}^j$ is admissible by (<ref>).

Formula (<ref>) for $\ell^j$ is a classical FAS construction; compare equation (8.5b) in [7] or equation (5.3.14) in [34].  It has the following explanation, inductively on VI problems (<ref>).  Because the down-smoother was only approximate, problem (<ref>) was not solved exactly by the finer-level correction $y^{j+1} \in \mathcal{D}^{j+1}$.  Suppose we seek an improved correction by adding $y \in \mathcal{D}^{j}$, to give $y^{j+1}+y \in \mathcal{D}^{j+1}+\mathcal{D}^j \subset \cV^{j+1}$, so that

\begin{equation}
\ip{f^{j+1}(w^{j+1}+y^{j+1}+y)}{v-y} \ge \ip{\ell^{j+1}}{v-y} \quad \text{\emph{[notional]}} \label{eq:fe:downvinotional}
\end{equation}

for all $v \in \mathcal{D}^j$.  (Note $v-y = (y^{j+1}+v)-(y^{j+1}+y)$ and that nesting $\cV^j \subset \cV^{j+1}$ is assumed here.)

However, because the down-smoother has made progress, $y$ should be approximate-able by solving a coarser, $j$th-level problem generated by modifying (<ref>) in 4 steps: (i) subtract the computable quantity $f^{j+1}(w^{j+1}+y^{j+1}) \in (\mathcal{V}^{j+1})'$ from both sides, (ii) replace the relaxed residual $\ell^{j+1}-f^{j+1}(w^{j+1}+y^{j+1})$ on the right by its restriction, (iii) replace $w^{j+1}+y^{j+1}$, where it appears on the left, by its injection, and (iv) where it appears on the left, replace $f^{j+1}$ by the coarser rediscretization $f^j$.  These steps yield a problem for the $j$th-level correction $y=y^j \in \mathcal{D}^j$:

\begin{align}
&\ip{f^j\left(\iR(w^{j+1}+y^{j+1})+y^j\right)}{v-y^j} - \ip{f^j\left(\iR(w^{j+1}+y^{j+1})\right)}{v-y^j} \label{eq:fe:downvicluttered} \\
&\qquad \ge \ip{R\left(\ell^{j+1}-f^{j+1}(w^{j+1}+y^{j+1})\right)}{v-y^j}, \notag
\end{align}

for all $v\in \mathcal{D}^{j}$.  Though in cluttered form, this is VI problem (<ref>), once definitions (<ref>) and (<ref>) are applied.

The explanation for (<ref>) and (<ref>) is simpler.  Whether or not problem (<ref>) has been exactly solved by the computed $y^j\in\mathcal{D}^j$, when going upward we can expand problem (<ref>) into the larger set $\mathcal{U}^j \supset \mathcal{D}^j$, and (<ref>) assures admissibly upward from any approximate solution to (<ref>).  Generally $y^j$ is not the solution to (<ref>); e.g. in the optimization case $f=F'$ it would not be optimal in the enlarged constraint set.  Our current best estimate of the solution to (<ref>) is $z^j \approx y^j + Pz^{j-1}\in \mathcal{U}^j$, so this is the initial iterate for relaxation.

These ideas come together in Algorithm <ref>, fascd-vcycle, which updates $w^J$.  The smooth and solve procedures are assumed to do in-place modifications of their final arguments, without modifying the others.

$\underline{\chi}^J, \, \overline{\chi}^J = \underline{\gamma}^J - w^J, \, \overline{\gamma}^J - w^J {\large \strut}$

for $j=J$ downto $j=1$
$\underline{\chi}^{j-1}, \, \overline{\chi}^{j-1} = \maxR \underline{\chi}^j, \, \minR \overline{\chi}^j {\large \strut}$

$\underline{\phi}^j, \, \overline{\phi}^j = \underline{\chi}^j - P\underline{\chi}^{j-1}, \, \overline{\chi}^j - P\overline{\chi}^{j-1} {\large \strut}$

$y^j = 0$

$\text{\pr{smooth}}^{\text{\id{down}}}(\ell^j,\underline{\phi}^j,\overline{\phi}^j,w^j;y^j)$  problem (<ref>) in $\mathcal{D}^j$

$w^{j-1} = \iR(w^j + y^j)$

$\ell^{j-1} = f^{j-1}(w^{j-1}) + R \left(\ell^j - f^j(w^j+y^j)\right)$
$z^0 = 0$

$\text{\pr{solve}}(\ell^0,\underline{\chi}^0,\overline{\chi}^0,w^0;z^0)$  problem (<ref>) in $\mathcal{U}^0$

for $j=1$ to $j=J$
$z^j = y^{j} + P z^{j-1}$

$\text{\pr{smooth}}^{\text{\id{up}}}(\ell^j,\underline{\chi}^j,\overline{\chi}^j,w^j;z^j)$  problem (<ref>) in $\mathcal{U}^j$
return $w^J+z^J$

The FASCD V-cycle, an iteration for solving VI problem (<ref>).

With 𝕀down $=1$ and 𝕀up $=1$ smoother iterations on each level, Algorithm <ref> is a $V(1,1)$ cycle.  The $V(1,0)$ cycle with 𝕀up $=0$ is similar to Algorithm 4.2 in [22], but here it is generalized by applying FAS-type corrections and allowing any smoother.  In unilateral cases one $V(1,1)$ iteration of fascd-vcycle approximates applications of cd-mult (Section <ref>) over CD (<ref>), followed by icd-tele (Section <ref>) over (<ref>), and in this sense $V(1,1)=V(1,0)+V(0,1)$.

Regarding the storage required, if all box constraints are nontrivial then 7 vectors must be allocated on each level: $\underline{\chi}^j$, $\overline{\chi}^j$, $\underline{\phi}^j$, $\overline{\phi}^j$, $w^j$, $\ell^j$, $y^j$; note vectors $z^j$ and $y^j$ may use the same storage.  On the finest level one must also store $\underline{\gamma}^J$ and $\overline{\gamma}^J$.  Therefore, using $m_j=2^d m_{j-1}$ and a geometric series argument, the total storage is at most

\begin{equation}
9 m_J + 7 m_{J-1} + \dots + 7 m_1 + 7 m_0 \le \left(9 + \frac{7}{2^d - 1}\right) m_J,
\end{equation}

namely $16m_J$, $12m_J$, and $10m_J$ in dimensions $d=1,2,3$, respectively.  For comparison, a single-level method needs $4 m_J$ storage (i.e. $\underline{\gamma}^J,\overline{\gamma}^J,\ell^J,w^J$).  In unilateral cases the storage can be reduced at all levels.

§ IMPLEMENTATION

In this section we propose a convergence criterion, construct the full multigrid (FMG) extension of Algorithm <ref>, and summarize the reduce-space Newton method used as the smoother and coarse solver.

Repeated application of Algorithm <ref> to solve finite-dimensional VI problem (<ref>) will not generally reduce the residual $f^J(w^J) - \ell^J$ to zero everywhere.  However, the VI is equivalent to a mixed complementarity problem (MCP) [18], which says, by definition, that certain conditions hold for $u^J$ at nodes $x_p \in \mathcal{T}^J$.  Let $r(u^J)_p = \ip{f^J(u^J)-\ell^J}{\psi_p}$ be the ordinary (pointwise) residual.  If $x_p \notin \partial_D\Omega$ then exactly one of these four conditions must hold:

* $\underline{\gamma}^J(x_p)<u^J(x_p)<\overline{\gamma}^J(x_p)$ and $r(u^J)_p = 0$,

* $\underline{\gamma}^J(x_p)=u^J(x_p)<\overline{\gamma}^J(x_p)$ and $r(u^J)_p \ge 0$,

* $\underline{\gamma}^J(x_p)<u^J(x_p)=\overline{\gamma}^J(x_p)$ and $r(u^J)_p \le 0$, or

* $\underline{\gamma}^J(x_p)=u^J(x_p)=\overline{\gamma}^J(x_p)$.

Otherwise, if $x_p \in \partial_D\Omega$ then

* $u^J(x_p)=g_D^J(x_p)$.

These 5 cases might be labeled inactive, lower active, upper active, pinched (both active), and boundary, respectively.  Where pinched or boundary apply, $r(u^J)_p$ can have any value.

This MCP interpretation allows us to construct expressions whose norms should approach zero.  Define the Fischer-Burmeister function [3, 35] on $a,b\in\RR$:

\begin{equation}
\phi(a,b) = a + b - \sqrt{a^2 + b^2}. \label{eq:phiFB}
\end{equation}

Observe that $\phi$ is semi-smooth, thus continuous, and that $\phi(a,b)=0$ if and only if $a,b$ are complementary, i.e. $a\ge 0$, $b\ge 0$, and $ab=0$.  For any $w^J\in \mathcal{K}^J$ denote $\underline{g}_p = w^J(x_p) - \underline{\gamma}^J(x_p)$ and $\overline{g}_p = \overline{\gamma}^J(x_p) - w^J(x_p)$, the nonnegative gaps at $x_p$.  Now define the semi-smooth residual

\begin{equation}
\rSS(w^J)_p = \begin{cases}
\max\big\{\phi\big(\underline{g}_p, r(w^J)_p\big), \phi\left(\overline{g}_p, -r(w^J)_p\right)\big\}, & x_p \notin \partial_D\Omega, \\
w^J(x_p) - g_D^J(x_p), & x_p \in \partial_D\Omega.
\end{cases} \label{eq:rSS}
\end{equation}

Straightforward modifications apply for unilaterally-constrained or unconstrained degrees of freedom, where $\underline{\gamma}^J(x_p) = -\infty$, $\overline{\gamma}^J(x_p) = +\infty$, or both.  VI problem (<ref>) is exactly solved by $u^J\in\cK^J$ if and only if $\rSS(u^J)=0$ identically.

Suppose that Algorithm <ref> is started with $w^{J,0} \in \mathcal{K}^J$ and generates (fine-level) iterates $w^{J,k} \in \mathcal{K}^J$.  Let $\id{atol}$, $\id{rtol}$, $\id{stol} \ge 0$ be given absolute, relative, and step tolerances, respectively.  We employ

\begin{equation}
\|\rSS(w^{J,k})\|_2 < \id{atol} \quad \text{or} \quad \frac{\|\rSS(w^{J,k})\|_2}{\|\rSS(w^{J,0})\|_2} < \id{rtol} \quad \text{or} \quad \frac{\|\Delta w^{J,k}\|_{L^2}}{\|w^{J,k}\|_{L^2}} < \id{stol}, \label{eq:stoppingcriterion}
\end{equation}

where $\Delta w^{J,k} = w^{J,k} -  w^{J,k-1}$, $\|\cdot\|_2$ is the euclidean norm on $\RR^{m_J}$, and $\|\cdot\|_{L^2}$ is the norm on $L^2(\Omega)$, as the stopping criterion.  Defaults in our Firedrake implementation match those of PETSc SNES [2]: $\id{atol}=10^{-50}$, $\id{rtol}=\id{stol}=10^{-8}$; see Section <ref> for specific settings in Examples.

In Section <ref> we will show that V-cycle iterations alone, starting with some finest-mesh initial iterate $w^{J,0}$ and repeated until (<ref>) holds, serve as an effective solver.  However, an FMG cycle, Algorithm <ref> below and Figure <ref> (right), is even better.  The well-known principle <cit.> is to start on the coarsest mesh and solve the problem from any initial iterate $w^0\in\mathcal{K}^0$; one may generate it by injection and projection/truncation of $w^{J,0}$.  Then an admissible initial iterate on the next-finer mesh is generated via prolongation.  However, in the VI context one must also project (truncate) to the bounds after the prolongation; see lines <ref> and <ref> in Algorithm <ref>.  During this “ramp” stage one approximately solves on each level using a fixed number (𝕀rampv) of V-cycles, each down to the coarsest mesh.  When the ramp is complete, one does full-depth V-cycles to convergence.

for $j=J$ downto $j=1$
$\ell^{j-1} = R(\ell^j)$

$\underline{\gamma}^{j-1}, \, \overline{\gamma}^{j-1} = \iR \underline{\gamma}^{j}, \, \iR \overline{\gamma}^{j}$
choose $w^0 \in \mathcal{K}^0$

for $j=1$ to $j=J-1$
$w^j = \max\left\{\underline{\gamma}^{j},\min\{\overline{\gamma}^{j}, Pw^{j-1}\}\right\} \in \mathcal{K}^j$

$w^J = \max\left\{\underline{\gamma}^{J},\min\{\overline{\gamma}^{J}, Pw^{J-1}\}\right\} \in \mathcal{K}^J$

while not (<ref>)
return $w^J$

The FASCD full multigrid (FMG) cycle for solving VI problem (<ref>).

In contrast to the construction of LDCs (Section <ref>) for V-cycles, during the ramp standard injection is used to generate the constraints $\underline{\gamma}^j,\overline{\gamma}^j$ on each level because a full solution iterate is initially needed, not just a correction.  Essentially, our choice is between (<ref>), where the coarse constraint set is nonempty but then the prolonged solution may not be admissible, or (<ref>), where the constraint set may be empty.  We choose (<ref>), with nonempty sets, but we must project (truncate) pointwise to the bounds, potentially adding high-frequency noise to the initial iterate.

A large space of smoothers could be explored, but in the next Section we will only apply an reduced-space Newton method [2, 3], denoted as RS Newton.  A fixed number of preconditioned Krylov iterations approximately solve the arising linear systems, thus $O(m_j)$ work is done per smoother application.  For the symmetric Laplacian (Example <ref>) we apply conjugate gradients (CG) with incomplete Cholesky (IC) preconditioning, but otherwise generalized minimum residual (GMRES) with incomplete LU (ILU) preconditioning is used.  The coarse solver is also RS Newton, iterated to convergence using the MUMPS sparse direct linear solver [1].

§ RESULTS

For $P_1$ and $Q_1$ elements, FASCD Algorithms <ref> and <ref> were implemented in Python using the Firedrake library [31].  The numerical results in this Section, which follow the Section <ref> examples, focus on the scaling of iterations with mesh resolution.  We see that each V-cycle greatly reduces the residual norm, and that a few FMG iterations solve most of these problems to within discretization error.

Consider the classical unilateral 2D Laplacian obstacle problem, namely Example <ref> with $p=2$ and admissible set $\mathcal{K} = \{v \ge \psi\}$.  We consider two cases, over square domains $\Omega$, with different coincidence (active) sets $\{x\in\Omega \,:\, u(x)=\psi(x)\}$ (Figure <ref>).  The standard “ball” problem has an upper hemisphere obstacle <cit.>, but the “spiral” problem's obstacle is constructed to produce the coincidence set shown <cit.>.  For both cases we used a regular coarsest mesh with $m_0=41$ nodes (64 triangles), and the finest mesh had $m_{10} \approx 3.4 \times 10^7$ nodes.  A single RS Newton iteration with three IC-preconditioned CG iterations was used as the smoother, and the initial iterate was $w=\max\{0,\psi\}$.

Table <ref> shows the resulting iteration counts.  For a first experiment we over-solve, with  $=$  $=$  $= 10^{-12}$ in stopping criterion (<ref>), to facilitate comparison to results in <cit.>.  As expected, RS Newton as a solver shows rapidly-increasing iterations, roughly proportional to resolution.  By contrast, V-cycle iterations increase slowly, perhaps logarithmically with the resolution.  For the spiral problem, <cit.> reports that a $m_J \approx 5.2 \times 10^5$ problem needs approximately 25 V-cycle iterations for the best (truncated;  “TNMG” and “TMMG”) monotone multilevel algorithms.[The stopping criterion and asymptotic convergence rates in [22] are based on norm differences from a last numerical iterate with unknown accuracy.  Our method directly measures satisfaction of the VI, via the residual norm $\|\rSS(w)\|_2$.]  The 12 V-cycle iterations needed here, at comparable resolution (8 levels), represent somewhat more arithmetic because our smoother is RS Newton using CG+IC, rather than the simpler projected Gauss-Seidel iteration.  Resolution-dependent asymptotic convergence rates for V-cycles, defined as $\rho_J=(\|\rSS(w^{J,k})\|_2/\|\rSS(w^{J,0})\|_2)^{1/k}$ where $w^{J,k}$ is the first iterate to satisfy (<ref>), are shown in Figure <ref>.  Both iterations and rates are similar for the ball and spiral problems; the geometric complexity of the coincidence set and free boundary minimally influences FASCD V-cycle performance.  Contrast the rate of about $0.1$ for 8 levels here with about $0.4$ reported by <cit.> for TNMG/TMMG at comparable resolution.  These rates are all excellent; our better rate mostly reflects a stronger smoother.

For FMG we instead used the defaults  $=$  $= 10^{-8}$, sufficient to reach discretization error, and in Table <ref> we report the number of V-cycles after the ramp; a small, constant number solves the problem.  Since each V-cycle does the work of about eight finest-level Krylov iterations, the ball results show nearly textbook multigrid efficiency at the highest resolutions despite the low solution regularity.  More V-cycles are needed to clean-up error remaining after the ramp in the spiral problem.

Coincidence sets for the ball and spiral solutions in Example <ref>.

2*levels     2*$m_J$     3c ball     3c spiral
(lr)3-5 (lr)6-8
RS only     V-cycle     FMG     RS only     V-cycle     FMG

1                 $41$       1      1      1       2      1      1

2                $145$       1      3      2       4      4      2

3                $545$       3      6      3       6      5      3

4     $2.1\times 10^3$       6      7      3       7      7      4

5     $8.3\times 10^3$      13      9      4      11      9      4

6     $3.3\times 10^4$      23     11      5      20     10      5

7     $1.3\times 10^5$      46     11      4      35     11      5

8     $5.3\times 10^5$      90     12      4      67     12      5

9     $2.1\times 10^6$         13      3         13      6

10     $8.4\times 10^6$         13      2         14      6

11     $3.4\times 10^7$         14      2         18      6

Iterations for FASCD V-cycle and FMG algorithms in Example <ref>; see text regarding tolerances.  Runs marked were not attempted.

V-cycle asymptotic convergence rates in Example <ref>.

We again consider Example <ref>, with a unilateral constraint, but now for a nonlinear $p$-Laplacian [14] in one dimension.  Setting $p=1.5$ yields a fast diffusion ($D=|\grad u|^{p-2}\to \infty$ as $|\grad u|\to 0$).  Note that Example <ref> below, by constrast, has a doubly-nonlinear operator with $p=4$ degenerate diffusion.

This example reveals the strong contrast between $V(1,0)$ and $V(0,1)$ results (up, down $=0$ in Algorithm <ref>, respectively).  Set $\Omega=(-3,3)$, $\psi(x) = -0.2|x|$, $\mathcal{K} = \{v \ge \psi\} \subset W^{1,p}(\Omega)$, and $\ip\ell v = \int_\Omega g v$, where $g(x)=1$ for $|x|<1$ and $g(x)=-1$ for $|x|>1$, in VI (<ref>).  The exact solution is easily computed for any $p$; see the codes accompanying this paper.  For this experiment the coarsest mesh has $m_0=7$ nodes and the finest $m_9=3073$.  The default smoother uses three iterations of RS Newton, with LU soluton of the linear systems.  We set  $= 10^{-6}$ and  $=$  $= 10^{-12}$.  Figure <ref> illustrates the major steps of a V-cycle, showing each of the coarse-correction VI subproblems and solutions; the result $\tilde w^2$ is visually close to the exact solution.

Visualization of one FASCD V-cycle iteration for $p$-Laplacian Example <ref>, with $J=3$ and a coarsest mesh of $m_0=5$ vertices.  Solutions (i.e. $w^j$, $y^j$, $z^j$) are lines, obstacles ($\underline{\gamma}^j$, $\underline{\phi}^j$, $\underline{\chi}^j$) are strong dots, and the zero function is dotted.  Vertical scale varies.

Iteration results are shown in Table <ref>.  We see that $V(1,0)$ cycles, Algorithm 4.2 in [22] but with a stronger smoother, often do not converge.  By contrast, $V(0,1)$ cycles, up-smoothing in the larger sets $\mathcal{U}^j$ (Section <ref>), are as good as the default $V(1,1)$ cycles (Algorithm <ref>).  FMG is even more efficient; a single iteration suffices at all higher resolutions.   Note that convergence is at rate $O(h^{2.0})$; compare [14].

levels     $m_J$     V(1,0)     V(0,1)     V(1,1)     FMG     $\|$error$\|_\infty$

1        $7$      1      1      1      1     $2.3\times 10^{-1}$

2       $13$      3      3      2      2     $3.3 \times 10^{-2}$

3       $25$     12      4      2      1     $9.1 \times 10^{-3}$

4       $49$     NC      4      2      1     $3.2 \times 10^{-3}$

5       $97$     NC      3      3      1     $5.5 \times 10^{-4}$

6      $193$     NC      3      3      1     $1.7 \times 10^{-4}$

7      $385$     NC      3      3      1     $4.7 \times 10^{-5}$

8      $769$     NC      4      7      1     $9.5 \times 10^{-6}$

9     $1537$     NC      5      6      1     $3.6 \times 10^{-6}$

10     $3073$     NC     11     19      1     $4.9 \times 10^{-7}$

Iterations for V-cycles and FMG in Example <ref>.  NC $=$ did not converge after 50 iterations.

Firedrake applications can exploit the PETSc [2] library's parallel solvers.  The above Examples, which were run in serial to fix the smoother, also run efficiently in parallel.  The next Examples apply parallel solvers and investigate the weak scaling of FASCD FMG.

Consider the advection-diffusion VI problem in Example <ref> with $\Omega=(-1,1)^d$, $d=2,3$, diffusivity $\eps=0.1$ in (<ref>), and zero Dirichlet boundary conditions.  The 2D problem satisfies the coercivity hypothesis, with a divergence-free velocity field $\bX = (7+5y,-5x)$, but in 3D we set $\bX = (7+5y,-5x,2z)$ so $\grad\cdot\bX \ne 0$.  The discontinuous source term $\phi \in L^\infty(\Omega)$ satisfies $\phi\ge 0$ in half of the domain, with positive circular patches, but values are negative in the other half; details are in the associated code.  While the other Examples in this Section are unilateral, here the solution is box-constrained: $0 \le u \le 1$.  Figure <ref> shows coincidence sets for the 2D problem; the 3D sets are comparable in regularity and general pattern.

Coincidence sets for the 2D results in Example <ref>: $u=0$ (left) and $u=1$ (right).

The 2D solver applied $P_1$ elements on a single core, but in 3D the elements were chosen to be $Q_1$ (hexahedra) on $16$ cores.  A coarse mesh with $m_0=16^d$ vertices minimally-resolved the coincidence set.  The local Péclet number $\text{Pe}=\|\bX(x,y,z)\|_2 L \eps^{-1}$, with $L=2$ the width of $\Omega$, varied from 40 to 260, and the mesh Péclet number $m_J^{-1/d}\,\text{Pe}$ varied further with resolution, so the $m_J\le 121^d$ problems were advection-dominated, and those with $m_J\ge 481^d$ diffusion-dominated.  However, the smoother was the same in all cases: two (fixed) RS Newton iterations, with three preconditioned (additive Schwarz method (ASM) with ILU on each process) GMRES iterations.  Tolerances were set to 𝕀rtol $=10^{-5}$ and 𝕀atol $=$ 𝕀stol $= 10^{-9}$.  FMG iteration results are shown in Table <ref>; a few iterations sufficed up to these memory-limited resolutions.

levels     $m_J$     $d=2$     $d=3$

1        $16^d$     1     1

2        $31^d$     2     2

3        $61^d$     2     2

4       $121^d$     2     2

5       $241^d$     2     3

6       $481^d$     2

7       $961^d$     3

8      $1921^d$     3

9      $3841^d$     3

FMG iterations for 2D ($P_1$, serial) and 3D ($Q_1$, 16 processes) solutions of the advection-diffusion bi-obstacle VI problem in Example <ref>.

Consider the steady shallow ice approximation (SIA) obstacle problem (Example <ref>), in which the surface elevation is constrained to exceed the bedrock elevation, $s\ge b$, on $\Omega=(0,L)^2$ with $L=1800$ km.  The data consist of a radial surface mass balance function $a$ <cit.> and a bumpy, synthetic map $b\in C^1(\Omega)$, with an elevation range of $1800$ m; details are in the associated code.  A numerical solution is shown in Figure <ref>, colored by the solution surface flow speed, the magnitude of $\mathbf{u} = (n+1)(n+2)^{-1} \Gamma (s-b)^{n+1} |\grad s|^{n-1} \grad s$.  Note $|\grad s|$ is unbounded at the free boundary; the solution $s$, while continuous, has very low regularity.

Example <ref> solution over a bumpy bedrock topography (grey), with 100-times vertical exaggeration and surface coloring by flow speed.

We applied FASCD FMG using $Q_1$ elements and a smoother of four (fixed) iterations of RS Newton (Section <ref>), with quadratic backtracking line search.  The arising linear systems were approximately solved by three (fixed) iterations of preconditioned GMRES and an ASM+ILU preconditioner.  The coarse mesh was relatively fine, with $m_0=21^2$ vertices and 400 cells.  Because of the low regularity of the solution we used 𝕀rtol $= 2 \times 10^{-4}$ and 𝕀atol $= 10^{-8}$.  The initial iterate was zero ice: $s=b$.  These choices generated a solver with mesh-independent iterations.

For a weak-scaling study we fixed the number of degrees of freedom per process at $m_J/P=641^2 \approx 4.1 \times 10^5$ and solved with $6, \dots, 11$ levels.[The experiments were conducted on ARCHER2, the UK national supercomputer. Each node has 16 memory channels, so a maximum of 16 cores per node were employed.]  The results in Table <ref> show some degradation in time per cycle with increasing core counts, perhaps due to the use of one-level additive Schwarz in the Krylov iteration.  Nevertheless, the finest mesh result, with $m_J=20481^2=4.1 \times 10^8$ degrees of freedom, required less than three minutes to solve for the steady surface elevation of an ice sheet larger than the current Greenland ice sheet ($1.7$ million $\text{km}^2$) at uniform $87.5$ m resolution.

$P$     levels     $m_J$     iterations     time (s)     time$/$iteration (s)

1     6     $641^2$     3     98     32.7

4     7     $1281^2$     3     98     32.7

16     8     $2561^2$     3     124     41.3

64     9     $5121^2$     3     136     45.3

256     10     $10241^2$     2     95     47.5

1024     11     $20481^2$     2     177     88.5

Iterations and wall clock time for parallel FASCD FMG solutions of Example <ref> on $P$ processes.

§ DISCUSSION AND OUTLOOK

As a multilevel solver for VI problems, the FASCD method is a strategy for generating coarser-level problems which, once solved, make helpful admissible additions to a V-cycle update.  The method is smoother-agnostic, and our implementation within the extensible Firedrake [31] and PETSc [2] library framework will allow easy experimentation with other smoothers (projected Gauss–Seidel, semi-smooth Newton, interior-point, etc.), beyond the active-set Newton method chosen for demonstration here.  For smoothers of Newton-Krylov type one could also add linear geometric or algebraic multigrid preconditioning [34] for the linear step problems, instead of our simpler choice of incomplete factorization preconditioners.  In a different direction, a restructured multilevel CD approach which applies the additive/parallel CD iteration [32] might be considered, in pursuit of a VI algorithm analogous to the BPX multigrid method for PDEs [4].  In any case, the performance observed in Section <ref> is already excellent on the tested examples.

§.§ Code availability
The software used to produce the results is archived at tag  in the repository , using the Firedrake version .

[1]
P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L'Excellent, A fully
asynchronous multifrontal solver using distributed dynamic scheduling, SIAM
Journal on Matrix Analysis and Applications, 23 (2001), pp. 15–41.

[2]
S. Balay and 34 others, PETSc/TAO users manual, Tech. Report
ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023.

[3]
S. Benson and T. Munson, Flexible complementarity solvers for
large-scale applications, Optimization Methods and Software, 21 (2006),
pp. 155–168.

[4]
J. H. Bramble, J. E. Pasciak, and J. Xu, Parallel multilevel
preconditioners, Mathematics of computation, 55 (1990), pp. 1–22.

[5]
A. Brandt, Multi-level adaptive solutions to boundary-value
problems, Math. Comp., 31 (1977), pp. 333–390.

[6]
A. Brandt and C. Cryer, Multigrid algorithms for the solution of
linear complementarity problems arising from free boundary problems, SIAM J.
Sci. Stat. Comput., 4 (1983), pp. 655–684.

[7]
A. Brandt and O. Livne, Multigrid Techniques: 1984 Guide with
Applications to Fluid Dynamics, SIAM Press, Philadelphia, revised ed.,

[8]
P. Brune, M. Knepley, B. Smith, and X. Tu, Composing scalable
nonlinear algebraic solvers, SIAM Rev., 57 (2015), pp. 535–565.

[9]
E. Bueler, Stable finite volume element schemes for the shallow ice
approximation, J. Glaciol., 62 (2016), pp. 230–242.

[10]
E. Bueler, Conservation laws for free-boundary fluid layers, SIAM
J. Appl. Math., 81 (2021), pp. 2007–2032,

[11]
E. Bueler, PETSc for Partial Differential Equations:
Numerical Solutions in C and Python, SIAM Press, Philadelphia, 2021.

[12]
C. Carstensen, Clément interpolation and its role in adaptive
finite element error control, vol. 168 of Operator Theory: Advances and
Applications, Birkhäuser Basel, 2006, ch. 2, pp. 27–43,

[13]
J. Chang and K. B. Nakshatrala, Variational inequality approach to
enforcing the non-negative constraint for advection–diffusion equations,
Comput. Methods Appl. Mech. Eng., 320 (2017), pp. 287–334.

[14]
H. J. Choe and J. L. Lewis, On the obstacle problem for quasilinear
elliptic equations of $p$-Laplacian type, SIAM J. Math. Anal., 22 (1991),
pp. 623–638.

[15]
I. Ekeland and R. Temam, Convex Analysis and Variational
Problems, North-Holland, Amsterdam, 1976.

[16]
H. C. Elman, D. J. Silvester, and A. J. Wathen, Finite Elements
and Fast Iterative Solvers: with Applications in Incompressible
Fluid Dynamics, Oxford University Press, Oxford, UK, 2nd ed., 2014.

[17]
L. C. Evans, Partial Differential Equations, Graduate Studies
in Mathematics, American Mathematical Society, Providence, 2nd ed., 2010.

[18]
F. Facchinei and J.-S. Pang, Finite-Dimensional Variational
Inequalities and Complementarity Problems, vol. 1, Springer, 2003.

[19]
P. E. Farrell, M. Croci, and T. M. Surowiec, Deflation for
semismooth equations, Optimization Methods and Software, 35 (2019),
pp. 1248–1271, <https://doi.org/10.1080/10556788.2019.1613655>.

[20]
P. E. Farrell and J. R. Maddison, Conservative interpolation between
volume meshes by local Galerkin projection, Computer Methods in Applied
Mechanics and Engineering, 200 (2011), pp. 89–100,

[21]
M. C. Ferris and J.-S. Pang, Engineering and economic applications
of complementarity problems, SIAM Rev., 39 (1997), pp. 669–713.

[22]
C. Gräser and R. Kornhuber, Multigrid methods for obstacle
problems, J. Comput. Math., 27 (2009), pp. 1–44.

[23]
R. Greve and H. Blatter, Dynamics of Ice Sheets and Glaciers,
Advances in Geophysical and Environmental Mechanics and Mathematics,
Springer, Berlin, Germany, 2009.

[24]
M. Hintermüller and I. Kopacka, A smooth penalty approach and a
nonlinear multigrid algorithm for elliptic MPECs, Computational
Optimization and Applications, 50 (2011), pp. 111–145.

[25]
G. Jouvet and E. Bueler, Steady, shallow ice sheets as obstacle
problems: well-posedness and finite element approximation, SIAM J. Appl.
Math., 72 (2012), pp. 1292–1314.

[26]
G. Jouvet and C. Gräser, An adaptive Newton multigrid method
for a model of marine ice sheets, J. Comput. Physics, 252 (2013),
pp. 419–437.

[27]
B. Keith and T. M. Surowiec, Proximal Galerkin: a
structure-preserving finite element method for pointwise bound constraints.,
2023, <https://doi.org/10.48550/arXiv.2307.12444>.

[28]
D. Kinderlehrer and G. Stampacchia, An Introduction to
Variational Inequalities and their Applications, Academic Press, New
York, 1980.

[29]
R. Kornhuber, Monotone multigrid methods for elliptic variational
inequalities, Numer. Math., 69 (1994), pp. 167–184.

[30]
I. Peral, Multiplicity of solutions for the $p$-Laplacian.
Notes of the Second International School in Functional Analysis and
Applications to Differential Equations, ICTP-Trieste, 1997.

[31]
F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T.
McRae, G.-T. Bercea, G. R. Markall, and P. H. J. Kelly, Firedrake:
automating the finite element method by composing abstractions, ACM Trans.
Math. Softw., 43 (2016), pp. 24:1–24:27.

[32]
X.-C. Tai, Rate of convergence for some constraint decomposition
methods for nonlinear variational inequalities, Numer. Math., 93 (2003),
pp. 755–786.

[33]
X.-C. Tai and P. Tseng, Convergence rate analysis of an asynchronous
space decomposition method for convex minimization, Math. Comput., 71
(2002), pp. 1105–1135.

[34]
U. Trottenberg, C. W. Oosterlee, and A. Schuller, Multigrid,
Elsevier, Oxford, UK, 2001.

[35]
M. Ulbrich, Semismooth Newton Methods for Variational Inequalities
and Constrained Optimization Problems in Function Spaces, SIAM Press,
Philadelphia, 2011.

[36]
J. Xu, Iterative methods by space decomposition and subspace
correction, SIAM Rev., 34 (1992), pp. 581–613.