# Mass-lumping discretization and solvers
for distributed elliptic optimal control problems

Ulrich Langer111Institute of Numerical Mathematics, Johannes Kepler University
Linz, and Johann Radon Institute for Computational and Applied Mathematics of
the Austrian Academy of Sciences, Altenberger Straße 69, 4040 Linz, Austria,
Email<EMAIL_ADDRESS>Richard Löscher222Institut für Angewandte
Mathematik, Technische Universität Graz, Steyrergasse 30, 8010 Graz, Austria,
Email<EMAIL_ADDRESS>Olaf Steinbach333Institut für Angewandte
Mathematik, Technische Universität Graz, Steyrergasse 30, 8010 Graz, Austria,
Email<EMAIL_ADDRESS>Huidong Yang444Faculty of Mathematics,
University of Vienna, Oskar–Morgenstern–Platz 1, 1090 Wien, Austria, and
Christian Doppler Laboratory for Mathematical Modeling and Simulation of Next
Generations of Ultrasound Devices (MaMSi), Oskar–Morgenstern–Platz 1, 1090
Wien, Austria, Email<EMAIL_ADDRESS>
###### Abstract

The purpose of this paper is to investigate the effects of the use of mass-
lumping in the finite element discretization of the reduced first-order
optimality system arising from a standard tracking-type, distributed elliptic
optimal control problem with $L_{2}$ regularization. We show that mass-lumping
will not affect the $L_{2}$ error between the desired state and the computed
state, but will lead to a Schur-complement system that allows for a fast
matrix-by-vector multiplication. We show that the use of the Schur-Complement
Preconditioned Conjugate Gradient method in a nested iteration setting leads
to an asymptotically optimal solver with respect to the complexity.
Keywords: Elliptic optimal control problems, $L_{2}$ regularization, finite
element discretization, mass lumping, preconditioned conjugate gradient
method, nested iteration.

##  1 Introduction

We consider the following tracking-type, distributed elliptic optimal control
problem with standard $L_{2}$ regularization: find the state $y_{\varrho}\in
Y=H_{0}^{1}(\Omega)$ and the control $u_{\varrho}\in U=L_{2}(\Omega)$
minimizing the cost functional

$J(y_{\varrho},u_{\varrho}):=\frac{1}{2}\,\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}^{2}+\frac{\varrho}{2}\,\|u_{\varrho}\|_{L_{2}(\Omega)}^{2},$
(1)

subject to (s.t.) the elliptic boundary value model problem

$-\Delta y_{\varrho}=u_{\varrho}\text{ in }\Omega,\quad y_{\varrho}=0\text{ on
}\partial\Omega,$ (2)

for some given desired state (target) $y_{d}\in L_{2}(\Omega)$, and some
regularization parameter $\varrho>0$, where $\Omega\subset\mathbb{R}^{d}$,
$d=1,2,3$, is a bounded Lipschitz domain with the boundary $\partial\Omega$.
We here use the standard notations for Lebesgue and Sobolev spaces. Since the
state equation (2) has a unique solution $y_{\varrho}\in Y$ for every given
control $u_{\varrho}\in U$, the optimal control problem (1)-(2) has a unique
solution $(y_{\varrho},u_{\varrho})\in Y\times U$ too; see, e.g., [29], [19],
or [41]. Moreover, the state $y_{\varrho}$ obviously belongs to
$H^{\Delta}(\Omega)=\\{y\in H_{0}^{1}(\Omega):\Delta y\in L_{2}(\Omega)\\}$,
and the solution operator $S$ mapping $u_{\varrho}$ to $y_{\varrho}$ (control-
to-state map) is an isomorphism between $L_{2}(\Omega)$ and
$H^{\Delta}(\Omega)$.

The finite element (fe) discretization of the reduced (after elimination of
the control $u_{\varrho}$) optimality system, which defines the solution to
the optimal control problem (1)-(2), leads to the solution of a large-scale
symmetric, but indefinite linear system of algebraic equations for defining
the fe nodal adjoint state vector $\mathbf{p}_{h}\in\mathbb{R}^{n_{h}}$ and
the fe nodal state vector $\mathbf{y}_{h}\in\mathbb{R}^{n_{h}}$ such that

$\begin{pmatrix}\varrho^{-1}M_{h}&K_{h}\\\
K_{h}&-M_{h}\end{pmatrix}\begin{pmatrix}\mathbf{p}_{h}\\\
\mathbf{y}_{h}\end{pmatrix}=\begin{pmatrix}\mathbf{0}_{h}\\\
-\mathbf{y}_{dh}\end{pmatrix},$ (3)

where the stiffness matrix $K_{h}$ and the mass matrix $M_{h}$ are symmetric
and positive definite (spd), $\mathbf{y}_{dh}\in\mathbb{R}^{n_{h}}$ is nothing
but the fe load vector representing the desired state $y_{d}$, and $h$ denotes
a suitable discretization parameter. For fixed $\varrho$, discretization error
estimates can be found, e.g., in [19]. There is a huge number of publications
on efficient preconditioned iterative solvers for symmetric, but indefinite
systems in general; see, e.g., the unified approach proposed in [42], the
survey paper [8], the review article [31], the books [15] and [6], the more
recent papers [1, 3, 4, 33], and the literature cited therein. Special
iterative solvers for discrete optimality systems such as (3) should be not
only robust with respect to (wrt) the mesh refinement quantified by the
discretization parameter $h$ but also wrt the regularization parameter
$\varrho$ that can be quite small depending on the cost that we are willing to
pay. Such kind of $h$ and $\varrho$ robust preconditioned iterative methods
have been proposed and investigated in [1, 5, 35, 37, 43]; see also [2, 14,
34, 36, 40], for handling control and state constraints, and the references
therein. Alternatively, we can use all-at-once multigrid methods to solve
saddle-point problems such as (3) efficiently; see, e.g., [38] and the review
paper [10].

In this paper, we are interested in the case $\varrho=h^{4}$ leading to
asymptotically optimal balanced estimates of the $L_{2}$-error between the
desired state $y_{d}$ and the computed finite element state $y_{\varrho h}$
that is related to $\mathbf{y}_{h}$ by the fe isomorphism; see [23].
Asymptotically optimal preconditioned iterative solvers for the saddle-point
system (3) were proposed in [23] and [22] for constant and variable $L_{2}$
regularizations, respectively. More precisely, it turns out that very cheap
preconditioners for the MINRES and BP-CG can be constructed on the basis of
simple diagonal approximations of the mass matrix $M_{h}$. Of course, we can
further reduce the saddle-point system (3) to the Schur-Complement (SC) system

$(\varrho K_{h}M_{h}^{-1}K_{h}+M_{h})\mathbf{y}_{h}=\mathbf{y}_{dh}$ (4)

by eliminating the adjoint state $\mathbf{p}_{h}$. The system matrix is spd,
and we would like to solve this system by means of the Preconditioned
Conjugate Gradient (PCG) method. Although we can use very cheap diagonal
matrices $D_{h}$ such as $\mbox{diag}(M_{h})$ or the lumped mass matrix
$\mbox{lump}(M_{h})$ as preconditioners that are spectrally equivalent to the
Schur complement $\varrho K_{h}M_{h}^{-1}K_{h}+M_{h}$ for $\varrho=h^{4}$ [22,
23], we cannot simply replace the mass matrix $M_{h}^{-1}$ by the lumped mass
matrix $(\mbox{lump}(M_{h}))^{-1}$ in the Schur complement without a precise
analysis of the impact of this replacement to the discretization error. In
Section 2, we just provide this analysis, and show that, in the case of
continuous, piecewise linear (Courant’s) finite element spaces, the
discretization error is not affected at all. This theoretical result is
supported by our numerical results presented in Section 4. Now we have to
solve the mass-lumped SC system

$(\varrho
K_{h}(\mbox{lump}(M_{h}))^{-1}K_{h}+M_{h})\hat{\mathbf{y}}_{h}=\mathbf{y}_{dh}$
(5)

instead of the original SC system (4). Using the diagonal preconditioner
$D_{h}=\mbox{lump}(M_{h})$, we can now solve (5) in asymptotically optimal
complexity for some fixed relative accuracy. In Section 3, we show how we can
use this SC-PCG in a nested iteration setting in order to compute a fe
approximation to the desired state $y_{d}$, which differs from $y_{d}$ in the
$L_{2}$-norm in the order of the discretization error, with asymptotically
optimal complexity $O(n_{h})$. These theoretical results are again
quantitatively illustrated by numerical experiments in Section 4.

##  2 Mass-Lumping and Error Analysis

The first-order optimality system, derived from (1)–(2), is given by the
equations

$-\Delta y_{\varrho}=u_{\varrho},\;-\Delta
p_{\varrho}=y_{\varrho}-y_{d},\;\mbox{and}\;p_{\varrho}+\varrho
u_{\varrho}=0\quad\mbox{in}\;\Omega,$ (6)

with the boundary conditions

$y_{\varrho}=0\;\mbox{and}\;p_{\varrho}=0\quad\mbox{on}\;\partial\Omega.$ (7)

Eliminating the control $u_{\varrho}$, we arrive at the reduced first-order
optimality system, the variational form of which reads as follows: find
$(y_{\varrho},p_{\varrho})\in H_{0}^{1}(\Omega)\times H_{0}^{1}(\Omega)$ such
that

$\displaystyle\frac{1}{\varrho}\left\langle
p_{\varrho},q\right\rangle_{L_{2}(\Omega)}+\left\langle\nabla
y_{\varrho},\nabla q\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle 0,\;\;\forall q\in H_{0}^{1}(\Omega),$ (8)
$\displaystyle-\left\langle\nabla p_{\varrho},\nabla
v\right\rangle_{L_{2}(\Omega)}+\left\langle
y_{\varrho},v\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle\left\langle y_{d},v\right\rangle_{L_{2}(\Omega)},\;\;\forall
v\in H_{0}^{1}(\Omega).$ (9)

Introducing the variable
$\tilde{p}_{\varrho}=\frac{1}{\sqrt{\varrho}}p_{\varrho}$, we further derive
the scaled system

$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle\tilde{p}_{\varrho},q\right\rangle_{L_{2}(\Omega)}+\left\langle\nabla
y_{\varrho},\nabla q\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle 0,\;\;\forall q\in H_{0}^{1}(\Omega),$ (10)
$\displaystyle-\left\langle\nabla\tilde{p}_{\varrho},\nabla
v\right\rangle_{L_{2}(\Omega)}+\frac{1}{\sqrt{\varrho}}\left\langle
y_{\varrho},v\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle
y_{d},v\right\rangle_{L_{2}(\Omega)},\;\;\forall v\in H_{0}^{1}(\Omega).$ (11)

For simplicity, we assume from now on that $\Omega\subset\mathbb{R}^{d}$ is
polygonally ($d=2$) or polyhedrally ($d=3$) bounded. Let
$\mathcal{T}_{h}=\\{\tau_{e}\\}_{e=1}^{N_{h}}$ be an admissible, globally
quasi-uniform and shape-regular decomposition of $\Omega$ into simplicies
$\tau_{e}$, with the mesh-size $h_{e}=|\tau_{e}|^{1/d}$, such that
$\overline{\Omega}=\bigcup_{e=1}^{N_{h}}\overline{\tau}_{e}$. Let
$S_{h}^{1}(\mathcal{T}_{h})=\text{span}\\{\varphi_{j}^{h}\\}_{j=1}^{\overline{n}_{h}}$
denote the space of piecewise linear, globally continuous functions spanned by
the Lagrange basis functions $\varphi_{j}^{h}$ (hat functions), which fulfil
the equations

$\displaystyle\sum_{j=1}^{\overline{n}_{h}}\varphi_{j}^{h}(x)=1\;\;\forall
x\in\Omega,\quad\text{and}\quad\varphi_{j}^{h}(x_{i})=\delta_{i,j}\text{ for
each node }x_{i},\,i=1,\ldots,\overline{n}_{h}.$ (12)

Further, we define $V_{h}:=S_{h}^{1}(\mathcal{T}_{h})\cap
H_{0}^{1}(\Omega)=\text{span}\\{\varphi_{j}^{h}\\}_{j=1}^{n_{h}}$, where we
assume that the ordering of the basis functions is such that the indices
$j=1,\ldots,n_{h}$ correspond to vertices $x_{j}\in\Omega$ and
$j=n_{h}+1,\ldots,\overline{n}_{h}$ corresponds to the vertices on the
boundary, $x_{j}\in\partial\Omega$. We refer to the books [11, 16, 39] for
more details on standard finite element discretizations of elliptic PDEs.

A conforming discretization of (10)-(11) is then to find $(y_{\varrho
h},\tilde{p}_{\varrho h})\in V_{h}\times V_{h}$ such that

$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle\tilde{p}_{\varrho
h},q_{h}\right\rangle_{L_{2}(\Omega)}+\left\langle\nabla y_{\varrho h},\nabla
q_{h}\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$ $\displaystyle
0,\;\;\forall q_{h}\in V_{h},$ (13)
$\displaystyle-\left\langle\nabla\tilde{p}_{\varrho h},\nabla
v_{h}\right\rangle_{L_{2}(\Omega)}+\frac{1}{\sqrt{\varrho}}\left\langle
y_{\varrho h},v_{h}\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle
y_{d},v_{h}\right\rangle_{L_{2}(\Omega)},\;\;\forall v_{h}\in V_{h}.$ (14)

In [23], we were able to show the following result for the $L_{2}$ error
between the desired state $y_{d}$ and the computed finite element state
$y_{\varrho h}$.

######  Theorem 1 ([23, Corollary 1]).

Let $(y_{\varrho h},\tilde{p}_{\varrho h})\in V_{h}\times V_{h}$ be the unique
solution of the coupled finite element variational formulation (13)-(14). Let
$y_{d}\in H_{0}^{s}(\Omega)$ for $s\in[0,1]$ or $y_{d}\in H^{s}(\Omega)\cap
H_{0}^{1}(\Omega)$ for $s\in(1,2]$. Then

$\displaystyle\|y_{\varrho h}-y_{d}\|_{L_{2}(\Omega)}\leq
c\,h^{s}\,\|y_{d}\|_{H^{s}(\Omega)},$

provided that $\varrho=h^{4}$.

We recall that $v_{h}\in V_{h}$ can be represented in the form

$\displaystyle v_{h}(x)=\sum_{i=1}^{n_{h}}v_{i}\varphi_{i}^{h}(x),$ (15)

where $v_{i}=v_{h}(x_{i})$. Thus, we can associate each finite element
function with its coefficient vector via the finite element isomorphism
$v_{h}\leftrightarrow\mathbf{v}$, where $\mathbf{v}[i]=v_{i}.$ With this, the
matrix system corresponding to the fe scheme (13)-(14) can be written in the
form: find
$(\mathbf{y}_{h},\mathbf{\tilde{p}}_{h})\in\mathbb{R}^{n_{h}}\times\mathbb{R}^{n_{h}}$
such that

$\displaystyle\begin{pmatrix}\frac{1}{\sqrt{\varrho}}M_{h}&K_{h}\\\
-K_{h}^{\top}&\frac{1}{\sqrt{\varrho}}M_{h}\end{pmatrix}\begin{pmatrix}\mathbf{\tilde{p}}_{h}\\\
\mathbf{y}_{h}\end{pmatrix}=\begin{pmatrix}\mathbf{0}\\\
\frac{1}{\sqrt{\varrho}}\mathbf{y}_{dh}\end{pmatrix},$ (16)

where $M_{h}$ and $K_{h}$ denote the mass resp. stiffness matrix with the
entries

$\displaystyle
M_{h}[i,j]=\int_{\Omega}\varphi_{j}^{h}(x)\varphi_{i}^{h}(x)\,dx\quad\text{and}\quad
K_{h}[i,j]=\int_{\Omega}\nabla\varphi_{j}^{h}(x)\cdot\nabla\varphi_{i}^{h}(x)\,dx,$

and the load vector

$\displaystyle\mathbf{y}_{dh}[i]=\int_{\Omega}y_{d}(x)\varphi_{i}^{h}(x)\,dx.$

We note that the system (16) is equivalent to (4). Moreover, when eliminating
$\mathbf{\tilde{p}}_{h}$ resp. $\mathbf{p}_{h}$, we arrive at the same Schur-
complement system (4).

As already mentioned in the introductionary Section 1, we would like to
replace the inverse of the mass matrix $M_{h}$ in the spd Schur complement
$\varrho K_{h}M_{h}^{-1}K_{h}+M_{h}$ by the inverse of the lumped mass matrix
$\text{lump}(M_{h})$ that is diagonal. The entries of the lumped mass matrix
$\text{lump}(M_{h})$ are given as

$\displaystyle\text{lump}(M_{h})[i,j]=\delta_{i,j}\sum_{k=1}^{\overline{n}_{h}}\widetilde{M}_{h}[i,k],\quad
i,j=1,\ldots,n_{h},$ (17)

where
$\widetilde{M}_{h}\in\mathbb{R}^{\overline{n}_{h}}\times\mathbb{R}^{\overline{n}_{h}}$
denotes the mass matrix on $S_{h}^{1}(\mathcal{T}_{h})$ with entries

$\widetilde{M}_{h}[i,j]:=\int_{\Omega}\varphi_{j}^{h}(x)\varphi_{i}^{h}(x)\,dx,\quad
i,j=1,\ldots,\overline{n}_{h}.$

The Schur complement system is then given by (5) with the Schur complement
$S_{h}=\varrho K_{h}(\text{lump}(M_{h}))^{-1}K_{h}+M_{h}$ as system matrix.
Unique solvability of the discrete system follows immediately, since
$M_{h}=M_{h}^{\top}>0$ is symmetric and positive definite (spd) and
$K_{h}(\text{lump}(M_{h}))^{-1}K_{h}>0$ is spd, as $\text{lump}(M_{h})$ and
$K_{h}$ are spd.

Now the aim is to show an equivalent result to Theorem 1, when using the
lumped mass matrix. We will exploit ideas from [7]. The discrete variational
formulation for the lumped case reads as follows: find $(\hat{y}_{\varrho
h},\hat{p}_{\varrho h})\in V_{h}\times V_{h}$ such that

$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle\hat{p}_{\varrho
h},q_{h}\right\rangle_{h}+\left\langle\nabla\hat{y}_{\varrho h},\nabla
q_{h}\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$ $\displaystyle
0,\;\;\forall q_{h}\in V_{h},$ (18)
$\displaystyle-\left\langle\nabla\hat{p}_{\varrho h},\nabla
v_{h}\right\rangle_{L_{2}(\Omega)}+\frac{1}{\sqrt{\varrho}}\left\langle\hat{y}_{\varrho
h},v_{h}\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\varrho}}\left\langle
y_{d},v_{h}\right\rangle_{L_{2}(\Omega)},\;\;\forall v_{h}\in V_{h},$ (19)

where $\left\langle
p_{h},q_{h}\right\rangle_{h}=\mathbf{q}_{h}^{\top}\text{lump}(M_{h})\mathbf{p}_{h}$
denotes the underintegrated inner product on $L_{2}(\Omega)$ that is nothing
but the realization of the mass lumping.

######  Lemma 1.

For $p_{h},q_{h}\in V_{h}$ the realization of the lumped mass matrix admits
the representation

$\displaystyle\left\langle
p_{h},q_{h}\right\rangle_{h}=\int_{\Omega}I_{h}^{1}(p_{h}q_{h})\,dx,$

where $I_{h}^{1}:\mathcal{C}(\overline{\Omega})\to V_{h}$ denotes the
interpolation operator, given as

$\displaystyle
I_{h}^{1}v(x)=\sum_{i=1}^{n_{i}}v_{i}\varphi_{i}^{h}(x),\,x\in\overline{\Omega},$

where $v_{i}=v(x_{i})$, $v\in\mathcal{C}(\overline{\Omega})$. Furthermore, it
holds that

$\displaystyle\frac{1}{d+2}\,\|p_{h}\|_{h}^{2}\leq\|p_{h}\|_{L_{2}(\Omega)}^{2}\leq\|p_{h}\|_{h}^{2}:=\left\langle
p_{h},p_{h}\right\rangle_{h}\text{ for all }p_{h}\in V_{h}.$ (20)

###### Proof.

With the representation (15), we have the coefficient vectors
$\mathbf{p}_{h}\leftrightarrow p_{h}$ and $\mathbf{q}_{h}\leftrightarrow
q_{h}$. Then we compute with (17), using (12),

$\displaystyle\left\langle p_{h},q_{h}\right\rangle_{h}$
$\displaystyle=\mathbf{q}_{h}^{\top}\text{lump}(M_{h})\mathbf{p}_{h}=\sum_{i=1}^{n_{h}}\sum_{j=1}^{n_{h}}p_{i}q_{j}\delta_{i,j}\sum_{k=1}^{\overline{n}_{h}}\widetilde{M}_{h}[i,k]$
$\displaystyle=\sum_{i=1}^{n_{h}}p_{i}q_{i}\sum_{k=1}^{\overline{n}_{h}}\widetilde{M}_{h}[i,k]=\sum_{i=1}^{n_{h}}p_{i}q_{i}\sum_{k=1}^{\overline{n}_{h}}\int_{\Omega}\varphi_{i}^{h}(x)\varphi_{k}^{h}(x)\,dx$
$\displaystyle=\int_{\Omega}\underbrace{\sum_{i=1}^{n_{h}}p_{i}q_{i}\varphi_{i}^{h}(x)}_{I_{h}^{1}(p_{h}q_{h})}\underbrace{\sum_{k=1}^{\overline{n}_{h}}\varphi_{k}^{h}(x)}_{=1}\,dx.$

The estimate (20) follows from, e.g., [39, Lemma 9.4],

$\frac{|\tau_{e}|}{(d+1)(d+2)}\sum\limits_{x_{i}\in\overline{\tau}_{e}}p_{i}^{2}\leq\|p_{h}\|^{2}_{L_{2}(\tau_{e})}\leq\frac{|\tau_{e}|}{d+1}\sum\limits_{x_{i}\in\overline{\tau}_{e}}p_{i}^{2},$

and

$\int_{\tau_{e}}I_{h}^{1}(p_{h}^{2})(x)\,dx=\sum\limits_{x_{i}\in\overline{\tau}_{e}}p_{i}^{2}\int_{\tau_{e}}\varphi_{i}^{h}(x)\,dx=\frac{|\tau_{e}|}{d+1}\sum\limits_{x_{i}\in\overline{\tau}_{e}}p_{i}^{2},$

when summing up over all elements $\tau_{e}$.

∎

With this representation, we can compute the consistency error.

######  Lemma 2.

Let $h=\max_{e=1,\ldots,N_{h}}h_{e}$. Then, for $p_{h},q_{h}\in V_{h}$, it
holds

$\displaystyle\left|\left\langle
p_{h},q_{h}\right\rangle_{L_{2}(\Omega)}-\left\langle
p_{h},q_{h}\right\rangle_{h}\right|$
$\displaystyle=\left|\int_{\Omega}\Big{[}p_{h}(x)q_{h}(x)-I_{h}^{1}(p_{h}q_{h})(x)\Big{]}\,dx\right|$
$\displaystyle\leq c\,h^{2}\Big{(}\varepsilon^{2}\,\|\nabla
p_{h}\|_{L_{2}(\Omega)}^{2}+\frac{1}{\varepsilon^{2}}\|\nabla
q_{h}\|_{L_{2}(\Omega)}^{2}\Big{)},$

for any $\varepsilon>0$.

###### Proof.

The first representation follows from Lemma 1. Let $\tau_{e}$ be a simplicial
finite element with the nodes $x_{e_{i}}$, $i=1,\ldots,d+1$. The associated
nodal values of a piecewise linear finite element function $p_{h}$ are the
coefficients $p_{e_{i}}$, $i=1,\ldots,d+1$. In particular, for $d=1$ and
$x\in\tau_{e}$, we then compute

$\displaystyle\int_{\tau_{e}}\Big{[}p_{h}(x)q_{h}(x)-I_{h}^{1}(p_{h}q_{h})(x)\Big{]}\,dx$
$\displaystyle=\,\int_{x_{e_{1}}}^{x_{e_{2}}}\left(\Big{[}p_{e_{1}}+\frac{x-x_{e_{1}}}{h_{e}}(p_{e_{2}}-p_{e_{1}})\Big{]}\Big{[}q_{e_{1}}+\frac{x-x_{e_{1}}}{h_{e}}(q_{e_{2}}-q_{e_{1}})\Big{]}\right.$
$\displaystyle\hskip
113.81102pt\left.-\Big{[}p_{e_{1}}q_{e_{1}}+\frac{x-x_{e_{1}}}{h_{e}}(p_{e_{2}}q_{e_{2}}-p_{e_{1}}q_{e_{1}})\Big{]}\right)dx$
$\displaystyle=\frac{1}{6}\,h_{e}\,(p_{e_{2}}-p_{e_{1}})\,(q_{e_{2}}-q_{e_{1}})\,=\,\frac{1}{6}\,h_{e}^{2}\,\int_{x_{e_{1}}}^{x_{e_{2}}}\frac{p_{e_{2}}-p_{e_{1}}}{h_{e}}\,\frac{q_{e_{2}}-q_{e_{1}}}{h_{e}}\,dx$
$\displaystyle=\frac{1}{6}\,h_{e}^{2}\,\int_{x_{e_{1}}}^{x_{e_{2}}}p_{h}^{\prime}(x)\,q_{h}^{\prime}(x)\,dx\,\leq\,\frac{1}{6}\,h_{e}^{2}\,\|\nabla_{x}p_{h}\|_{L^{2}(\tau_{e})}\|\nabla_{x}q_{h}\|_{L^{2}(\tau_{e})}\,.$

For $d=2$ and $x\in\tau_{e}$, we introduce the representation
$x=x_{e_{1}}+J_{e}\eta$ with respect to the reference element
$\tau=\\{\eta\in{\mathbb{R}}^{2}:\eta_{1}\in(0,1),\eta_{2}\in(0,1-\eta_{1})\\}$
and we write $p_{h}(x)=p_{h}(x_{e_{1}}+J_{e}\eta)=\widetilde{p}_{h}(\eta)$,
$\eta\in\tau$. Similar as in the case $d=1$ we then compute, using
$\mbox{det}\,J_{e}=2\,|\tau_{e}|$,

$\displaystyle\int_{\tau_{e}}\Big{[}p_{h}(x)q_{h}(x)-I_{h}^{1}(p_{h}q_{h})(x)\Big{]}\,dx\,=\,\int_{\tau}\Big{[}\widetilde{p}_{h}(\eta)\widetilde{q}_{h}(\eta)-I_{h}^{1}(\widetilde{p}_{h}\widetilde{q}_{h})(\eta)\Big{]}\,\mbox{det}J_{e}\,d\eta$
$\displaystyle=\frac{|\tau_{e}|}{12}\Big{[}(p_{0}-p_{2})(q_{2}-q_{0})+(p_{1}-p_{0})(q_{0}-q_{1})+(p_{1}-p_{2})(q_{2}-q_{1})\Big{]}$
$\displaystyle\leq\frac{|\tau_{e}|}{12}\Big{[}(p_{0}-p_{2})^{2}+(p_{1}-p_{0})^{2}+(p_{1}-p_{2})^{2}\Big{]}^{1/2}$
$\displaystyle\hskip
85.35826pt\cdot\Big{[}(q_{2}-q_{0})^{2}+(q_{0}-q_{1})^{2}+(q_{2}-q_{1})^{2}\Big{]}^{1/2}\,.$

With

$(p_{1}-p_{2})^{2}=(p_{1}-p_{0}+p_{0}-p_{2})^{2}\leq
2\,(p_{1}-p_{0})^{2}+2\,(p_{0}-p_{2})^{2},$

we further have, e.g., [39, Lemma 9.1],

$\displaystyle\int_{\tau_{e}}\Big{[}p_{h}(x)q_{h}(x)-I_{h}^{1}(p_{h}q_{h})(x)\Big{]}\,dx$
$\displaystyle\hskip
28.45274pt\leq\,\frac{|\tau_{e}|}{4}\Big{[}(p_{0}-p_{2})^{2}+(p_{1}-p_{0})^{2}\Big{]}^{1/2}\Big{[}(q_{2}-q_{0})^{2}+(q_{0}-q_{1})^{2}\Big{]}^{1/2}$
$\displaystyle\hskip
28.45274pt=\,\frac{|\tau_{e}|}{4}\left[2\int_{\tau}|\nabla_{\eta}\widetilde{p}_{h}|^{2}d\eta\right]^{1/2}\left[2\int_{\tau}|\nabla_{\eta}\widetilde{q}_{h}|^{2}d\eta\right]^{1/2}$
$\displaystyle\hskip
28.45274pt=\,\frac{|\tau_{e}|}{2}\,\|\nabla_{\eta}\widetilde{p}_{h}\|_{L^{2}(\tau)}\|\nabla_{\eta}\widetilde{q}_{h}\|_{L^{2}(\tau)}\,\leq\,c\,h_{e}^{2}\,\|\nabla_{x}p_{h}\|_{L^{2}(\tau_{e})}\|\nabla_{x}q_{h}\|_{L^{2}(\tau_{e})}\,.$

For $d=3$, we proceed in the same way. Now the reference element is given by
$\tau=\\{\eta\in{\mathbb{R}}^{3}:\eta_{1}\in(0,1),\eta_{2}\in(0,1-\eta_{1}),\eta_{3}\in(0,1-\eta_{1}-\eta_{2})\\}$,
and $\mbox{det}\,J_{e}=6\,|\tau_{e}|$. Then,

$\displaystyle\int_{\tau_{e}}\Big{[}p_{h}(x)q_{h}(x)-I_{h}^{1}(p_{h}q_{h})(x)\Big{]}\,dx$
$\displaystyle=\,6\,|\tau_{e}|\,\int_{0}^{1}\int_{0}^{1-\eta_{1}}\int_{0}^{1-\eta_{1}-\eta_{2}}\Big{[}\widetilde{p}_{h}(\eta)\widetilde{q}_{h}(\eta)-I_{h}^{1}(\widetilde{p}_{h}\widetilde{q}_{h})(\eta)\Big{]}d\eta_{3}d\eta_{2}d\eta_{1}$
$\displaystyle=\frac{|\tau_{e}|}{20}\Big{[}(p_{0}-p_{1})(q_{1}-q_{0})+(p_{0}-p_{2})(q_{2}-q_{0})+(p_{0}-p_{3})(q_{3}-q_{0})$
$\displaystyle\hskip
85.35826pt+(p_{1}-p_{2})(q_{2}-q_{1})+(p_{1}-p_{3})(q_{3}-q_{1})+(p_{2}-p_{3})(q_{3}-q_{2})\Big{]}$
$\displaystyle\leq\frac{|\tau_{e}|}{20}\Big{[}(p_{0}-p_{1})^{2}+(p_{0}-p_{2})^{2}+(p_{0}-p_{3})^{2}+(p_{1}-p_{2})^{2}+(p_{1}-p_{3})^{2}+(p_{2}-p_{3})^{2}\Big{]}^{1/2}$
$\displaystyle\hskip
28.45274pt\cdot\Big{[}(q_{1}-q_{0})^{2}+(q_{2}-q_{0})^{2}+(q_{3}-q_{0})^{2}+(q_{2}-q_{1})^{2}+(q_{3}-q_{1})^{2}+(q_{3}-q_{2})^{2}\Big{]}^{1/2}$
$\displaystyle\leq\frac{|\tau_{e}|}{4}\Big{[}(p_{0}-p_{1})^{2}+(p_{0}-p_{2})^{2}+(p_{0}-p_{3})^{2}\Big{]}^{1/2}\Big{[}(q_{1}-q_{0})^{2}+(q_{2}-q_{0})^{2}+(q_{3}-q_{0})^{2}\Big{]}^{1/2}$
$\displaystyle=\frac{|\tau_{e}|}{4}\left[6\int_{\tau}|\nabla_{\eta}\widetilde{p}_{h}|^{2}d\eta\right]^{1/2}\left[6\int_{\tau}|\nabla_{\eta}|\widetilde{q}_{h}|^{2}d\eta\right]^{1/2}$
$\displaystyle=\frac{3}{2}\,|\tau_{e}|\,\|\nabla_{\eta}\widetilde{p}_{h}\|_{L^{2}(\tau)}\|\nabla_{\eta}\widetilde{q}_{h}\|_{L^{2}(\tau)}\,\leq\,c\,h_{e}^{2}\,\|\nabla_{x}p_{h}\|_{L^{2}(\tau_{e})}\|\nabla_{x}q_{h}\|_{L^{2}(\tau_{e})}\,.$

Hence, using Young’s inequality,

$\|\nabla_{x}p_{h}\|_{L_{2}(\tau_{e})}\|\nabla_{x}q_{h}\|_{L_{2}(\tau_{e})}\leq\frac{1}{2}\,\left(\varepsilon^{2}\,\|\nabla_{x}p_{h}\|^{2}_{L_{2}(\tau_{e})}+\frac{1}{\varepsilon^{2}}\,\|\nabla_{x}q_{h}\|^{2}_{L_{2}(\tau_{e})}\right),$

and summing up over all elements $\tau_{e}$, this gives the desired estimate.
∎

We need one more preliminary result, before we can state the main theorem.

######  Lemma 3.

Let $(y_{\varrho},\tilde{p}_{\varrho})\in H_{0}^{1}(\Omega)\times
H_{0}^{1}(\Omega)$ be the unique solution of the reduced optimality system
(10) and (11). Then there holds the regularization error estimate

$\displaystyle\|y_{\varrho}-y_{d}\|_{H^{-1}(\Omega)}\leq
c\,\sqrt{\varrho}\,|y_{d}|_{H^{1}(\Omega)}\quad\text{for }y_{d}\in
H_{0}^{1}(\Omega).$

Additionally, for $y_{d}\in H^{\Delta}(\Omega)\cap H_{0}^{1}(\Omega)$, there
holds

$\displaystyle\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}\leq\sqrt{\varrho}\,\|\Delta
y_{d}\|_{L_{2}(\Omega)}\quad\text{and}\quad\|\Delta
y_{\varrho}\|_{L_{2}(\Omega)}\leq\|\Delta y_{d}\|_{L_{2}(\Omega)}.$

###### Proof.

The first estimate is given in [32, Theorem 4.1, (4.7)]. The second and third
estimate can be found in [23, Lemma 1, (2.5)] and in the proof of this lemma.
But for clarity, we will recall the proof. We note that, by the optimality
system, we have the equations $-\Delta y_{\varrho}=u_{\varrho}$, and
$p_{\varrho}=-\varrho u_{\varrho}$. First, assuming the regularity
$y_{\varrho},y_{d}\in H^{\Delta}(\Omega)\cap H_{0}^{1}(\Omega)$, using (8) and
(9), and integration by parts, we obtain

$\displaystyle\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}^{2}$
$\displaystyle=\left\langle
y_{\varrho}-y_{d},y_{\varrho}-y_{d}\right\rangle_{L_{2}(\Omega)}=\left\langle\nabla
p_{\varrho},\nabla(y_{\varrho}-y_{d})\right\rangle_{L_{2}(\Omega)}$
$\displaystyle=\left\langle p_{\varrho},-\Delta
y_{\varrho}\right\rangle_{L_{2}(\Omega)}+\left\langle p_{\varrho},\Delta
y_{d}\right\rangle_{L_{2}(\Omega)}$ $\displaystyle=-\varrho\,\left\langle
u_{\varrho},-\Delta
y_{\varrho}\right\rangle_{L_{2}(\Omega)}-\varrho\,\left\langle
u_{\varrho},\Delta y_{d}\right\rangle_{L_{2}(\Omega)}$
$\displaystyle=-\varrho\,\|\Delta
y_{\varrho}\|_{L_{2}(\Omega)}^{2}+\varrho\,\left\langle\Delta
y_{\varrho},\Delta y_{d}\right\rangle_{L_{2}(\Omega)}.$

From this we conclude

$\displaystyle\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}^{2}+\varrho\,\|\Delta
y_{\varrho}\|_{L_{2}(\Omega)}^{2}\leq\varrho\,\|\Delta
y_{\varrho}\|_{L_{2}(\Omega)}\|\Delta y_{d}\|_{L_{2}(\Omega)},$

and further

$\displaystyle\|\Delta y_{\varrho}\|_{L_{2}(\Omega)}\leq\|\Delta
y_{d}\|_{L_{2}(\Omega)}\quad\text{and}\quad\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}\leq\sqrt{\varrho}\,\|\Delta
y_{d}\|_{L_{2}(\Omega)}.$

From the first estimate we conclude that $y_{d}\in H^{\Delta}(\Omega)\cap
H_{0}^{1}(\Omega)$ implies $y_{\varrho}\in H^{\Delta}(\Omega)\cap
H_{0}^{1}(\Omega)$. This confirms that it is sufficient to require the
regularity on $y_{d}$ only. ∎

The main statement of this paper is formulated in the following theorem.

######  Theorem 2.

Let $(\hat{y}_{\varrho h},\hat{p}_{\varrho h})\in V_{h}\times V_{h}$ be the
unique solution of the variational formulation (18) and (19). Assume that
$\mathcal{T}_{h}$ is globally quasi-uniform such that a global inverse
inequality holds true. Further, choose $\varrho=h^{4}$. Then,

$\displaystyle\|\hat{y}_{\varrho
h}-y_{d}\|_{L_{2}(\Omega)}\leq\begin{cases}c\,h\,\|y_{d}\|_{H^{1}_{0}(\Omega)},&\text{if
}y_{d}\in H_{0}^{1}(\Omega),\\\\[2.84526pt]
c\,h^{2}\,\|y_{d}\|_{H^{2}(\Omega)},&\text{if }y_{d}\in H^{2}(\Omega)\cap
H_{0}^{1}(\Omega)\text{ and }\Omega\text{ is convex.}\end{cases}$

###### Proof.

Let $(y_{\varrho h},\tilde{p}_{\varrho h})\in V_{h}\times V_{h}$ be the unique
solution of (13) and (14). By the triangle inequality, we get that

$\displaystyle\|\hat{y}_{\varrho
h}-y_{d}\|_{L_{2}(\Omega)}\leq\|\hat{y}_{\varrho h}-y_{\varrho
h}\|_{L_{2}(\Omega)}+\|y_{\varrho h}-y_{d}\|_{L_{2}(\Omega)}.$

By Theorem 1, the second term fulfils the estimate. Thus, it is sufficient to
bound the first term. Therefore, subtracting the variational formulation (13)
and (14) from (18) and (19) with $v_{h}=\hat{y}_{\varrho h}-y_{\varrho h}$ and
$q_{h}=\hat{p}_{\varrho h}-\tilde{p}_{\varrho h}$, we obtain the equalities

$\displaystyle\frac{1}{\sqrt{\varrho}}\,\|\hat{y}_{\varrho h}-$ $\displaystyle
y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}\,=\,\frac{1}{\sqrt{\varrho}}\,\left\langle\hat{y}_{\varrho
h}-y_{\varrho h},\hat{y}_{\varrho h}-y_{\varrho
h}\right\rangle_{L_{2}(\Omega)}$
$\displaystyle=\left\langle\nabla(\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}),\nabla(\hat{y}_{\varrho h}-y_{\varrho h})\right\rangle_{L_{2}(\Omega)}$
$\displaystyle=\frac{1}{\sqrt{\varrho}}\,\left(\left\langle\tilde{p}_{\varrho
h},\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\right\rangle_{L_{2}(\Omega)}-\left\langle\hat{p}_{\varrho
h},\hat{p}_{\varrho h}-\tilde{p}_{\varrho h}\right\rangle_{h}\right)$
$\displaystyle=\frac{1}{\sqrt{\varrho}}\,\left(\left\langle\tilde{p}_{\varrho
h},\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\right\rangle_{L_{2}(\Omega)}-\|\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\|_{h}^{2}-\left\langle\tilde{p}_{\varrho h},\hat{p}_{\varrho
h}-\tilde{p}_{\varrho h}\right\rangle_{h}\right).$

Multiplying by $\sqrt{\varrho}$ and using Lemma 1, we further get

$\displaystyle\|\hat{y}_{\varrho h}-y_{\varrho h}\|_{L_{2}(\Omega)}^{2}$
$\displaystyle+\|\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}\leq\|\hat{y}_{\varrho h}-y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\|_{h}^{2}$ $\displaystyle=\left\langle\tilde{p}_{\varrho
h},\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\right\rangle_{L_{2}(\Omega)}-\left\langle\tilde{p}_{\varrho
h},\hat{p}_{\varrho h}-\tilde{p}_{\varrho h}\right\rangle_{h}$
$\displaystyle=\int_{\Omega}\left[\tilde{p}_{\varrho h}(x)(\hat{p}_{\varrho
h}(x)-\tilde{p}_{\varrho h}(x))-I_{h}^{1}(\tilde{p}_{\varrho
h}(\hat{p}_{\varrho h}-\tilde{p}_{\varrho h})(x)\right]\,dx.$

With Lemma 2, choosing $p_{h}=\tilde{p}_{\varrho h}$ and
$q_{h}=\hat{p}_{\varrho h}-\tilde{p}_{\varrho h}$, we estimate, for some
$\varepsilon>0$ to be specified,

$\displaystyle\|\hat{y}_{\varrho h}-y_{\varrho h}\|_{L_{2}(\Omega)}^{2}$
$\displaystyle+\|\hat{p}_{\varrho h}-\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}$ $\displaystyle\leq
c\,h^{2}\,\left(\varepsilon^{2}\,\|\nabla\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\frac{1}{\varepsilon^{2}}\,\|\nabla(\hat{p}_{\varrho
h}-\tilde{p}_{\varrho h})\|_{L_{2}(\Omega)}^{2}\right).$

Using an inverse inequality, we estimate the second term by

$\displaystyle\frac{ch^{2}}{\varepsilon^{2}}\,\|\nabla(\hat{p}_{\varrho
h}-\tilde{p}_{\varrho
h})\|_{L_{2}(\Omega)}^{2}\leq\frac{cc_{I}}{\varepsilon^{2}}\,\|\hat{p}_{\varrho
h}-\tilde{p}_{\varrho h}\|_{L_{2}(\Omega)}^{2}=\|\hat{p}_{\varrho
h}-\tilde{p}_{\varrho h}\|_{L_{2}(\Omega)}^{2},$

when choosing $\varepsilon=\sqrt{cc_{I}}$. Thus, it holds

$\displaystyle\|\hat{y}_{\varrho h}-y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}\leq\tilde{c}\,h^{2}\,\|\nabla\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}.$ (21)

Now it is sufficient to bound $\|\nabla\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}$ suitably. Let $(y_{\varrho},\tilde{p}_{\varrho})\in
H_{0}^{1}(\Omega)\times H_{0}^{1}(\Omega)$ be the unique solution of the
coupled variational formulation (8) and (9). Using the triangle inequality and
the trivial inequality $(a+b)^{2}\leq 2(a^{2}+b^{2})$, we get

$\displaystyle\|\nabla\tilde{p}_{\varrho h}\|_{L_{2}(\Omega)}^{2}\leq
2\,\Big{(}\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}^{2}+\|\nabla(\tilde{p}_{\varrho
h}-\tilde{p}_{\varrho})\|_{L_{2}(\Omega)}^{2}\Big{)}.$ (22)

For $(y_{\varrho},\widetilde{p}_{\varrho})\in H_{0}^{1}(\Omega)\times
H_{0}^{1}(\Omega)$ and $(y_{\varrho h},\widetilde{p}_{\varrho h})\in
V_{h}\times V_{h}$ as solutions of (8)-(9) and (13)-(14), respectively, we can
show, as in the proof of [23, Theorem 1], using an inverse inequality and
$\varrho=h^{4}$, that Cea’s Lemma

$\displaystyle h^{-2}\,\|y_{\varrho}-y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(y_{\varrho}-y_{\varrho
h})\|_{L_{2}(\Omega)}^{2}$ $\displaystyle\hskip
28.45274pt+\,h^{-2}\,\|\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h})\|_{L_{2}(\Omega)}^{2}$ $\displaystyle\hskip
56.9055pt\leq\,c\,\Big{[}h^{-2}\,\|y_{\varrho}-v_{h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(y_{\varrho}-v_{h})\|_{L_{2}(\Omega)}^{2}$
$\displaystyle\hskip
85.35826pt+\,h^{-2}\,\|\tilde{p}_{\varrho}-q_{h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(\tilde{p}_{\varrho}-q_{h})\|_{L_{2}(\Omega)}^{2}\Big{]}$

holds true for all $(v_{h},q_{h})\in V_{h}\times V_{h}$. Further, using best
approximation results, we get, for $y_{d}\in H_{0}^{1}(\Omega)$, that

$\displaystyle h^{-2}\,\|y_{\varrho}-y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(y_{\varrho}-y_{\varrho
h})\|_{L_{2}(\Omega)}^{2}$ $\displaystyle\hskip
28.45274pt+\,h^{-2}\,\|\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h})\|_{L_{2}(\Omega)}^{2}\,\leq\,c\,|y_{d}|_{H^{1}(\Omega)}^{2},$

and, for $y_{d}\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$, that

$\displaystyle h^{-2}\,\|y_{\varrho}-y_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(y_{\varrho}-y_{\varrho
h})\|_{L_{2}(\Omega)}^{2}$ $\displaystyle\hskip
28.45274pt+\,h^{-2}\,\|\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}^{2}+\|\nabla(\tilde{p}_{\varrho}-\tilde{p}_{\varrho
h})\|_{L_{2}(\Omega)}^{2}\,\leq\,c\,h^{2}\,|y_{d}|_{H^{2}(\Omega)}^{2}.$

From this we immediately conclude that

$\displaystyle\|\nabla(\tilde{p}_{\varrho
h}-\tilde{p}_{\varrho})\|_{L_{2}(\Omega)}^{2}\leq\begin{cases}c\,|y_{d}|_{H^{1}(\Omega)}^{2},&\text{for
}y_{d}\in H^{1}_{0}(\Omega),\\\\[2.84526pt]
c\,h^{2}\,|y_{d}|_{H^{2}(\Omega)}^{2},&\text{for }y_{d}\in
H_{0}^{1}(\Omega)\cap H^{2}(\Omega).\end{cases}$ (23)

For the first term, we use the first inequality of Lemma 3 for $y_{d}\in
H_{0}^{1}(\Omega)$ to estimate

$\displaystyle\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}^{2}$
$\displaystyle=\left\langle\nabla\tilde{p}_{\varrho},\nabla\tilde{p}_{\varrho}\right\rangle_{L_{2}(\Omega)}=\frac{1}{\sqrt{\varrho}}\,\left\langle
y_{\varrho}-y_{d},\tilde{p}_{\varrho}\right\rangle_{L_{2}(\Omega)}$
$\displaystyle\leq\frac{1}{\sqrt{\varrho}}\,\|y_{\varrho}-y_{d}\|_{H^{-1}(\Omega)}\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}\,\leq\,c\,|y_{d}|_{H^{1}(\Omega)}\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)},$

from which we conclude

$\displaystyle\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}\leq
c\,|y_{d}|_{H^{1}(\Omega)}.$ (24)

For a convex domain $\Omega$, we have $H^{\Delta}(\Omega)=H^{2}(\Omega)$. Thus
$y_{d}\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$, and we get with the same
reasoning, using the second inequality of Lemma 3,

$\displaystyle\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}^{2}$
$\displaystyle=\frac{1}{\sqrt{\varrho}}\,\left\langle
y_{\varrho}-y_{d},\tilde{p}_{\varrho}\right\rangle_{L_{2}(\Omega)}$
$\displaystyle\leq\frac{1}{\sqrt{\varrho}}\,\|y_{\varrho}-y_{d}\|_{L_{2}(\Omega)}\|\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}\,\leq\,|y_{d}|_{H^{2}(\Omega)}\|\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}.$

Now, we recall that, from the optimality system (6)-(9), we have

$\displaystyle\tilde{p}_{\varrho}=\frac{1}{\sqrt{\varrho}}\,p_{\varrho}=-\sqrt{\varrho}\,u_{\varrho}=\sqrt{\varrho}\,\Delta
y_{\varrho}\quad\text{in }\Omega,$

and thus, with Lemma 3, we obtain

$\displaystyle\|\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}=\sqrt{\varrho}\,\|\Delta
y_{\varrho}\|_{L_{2}(\Omega)}\leq\sqrt{\varrho}\,\|\Delta
y_{d}\|_{L_{2}(\Omega)}\leq\sqrt{\varrho}\,|y_{d}|_{H^{2}(\Omega)}.$

Thus, for $\varrho=h^{4}$ and $y_{d}\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$,
we arrive at the estimate

$\displaystyle\|\nabla\tilde{p}_{\varrho}\|_{L_{2}(\Omega)}\leq\,h\,|y_{d}|_{H^{2}(\Omega)}.$
(25)

Now, combining (21) with (22), (23),(24), and (25), this gives

$\displaystyle\|\hat{y}_{\varrho h}-y_{\varrho h}\|_{L_{2}(\Omega)}\leq
h\,\|\nabla\tilde{p}_{\varrho
h}\|_{L_{2}(\Omega)}\leq\begin{cases}c\,h\,|y_{d}|_{H^{1}(\Omega)},&\text{for
}y_{d}\in H_{0}^{1}(\Omega),\\\\[2.84526pt]
c\,h^{2}\,|y_{d}|_{H^{2}(\Omega)},&\text{for }y_{d}\in H^{2}(\Omega)\cap
H_{0}^{1}(\Omega).\end{cases}$

∎

######  Lemma 4.

Let $(\hat{y}_{\varrho h},\hat{p}_{\varrho h})\in V_{h}\times V_{h}$ be the
unique solution of the coupled variational formulation (18) and (19). Then,
for $y_{d}\in L_{2}(\Omega)$, we get the error estimate

$\displaystyle\|\hat{y}_{\varrho
h}-y_{d}\|_{L_{2}(\Omega)}\leq\|y_{d}\|_{L_{2}(\Omega)}.$

###### Proof.

Choosing $q_{h}=\hat{p}_{\varrho h}$ and $v_{h}=\hat{y}_{\varrho h}$ in (18)
and (19), summing up the equations, and multiplying with $\sqrt{\varrho}$,
this gives

$\displaystyle\left\langle\hat{p}_{\varrho
h},\hat{p}_{\varrho}\right\rangle_{h}+\left\langle\hat{y}_{\varrho
h},\hat{y}_{\varrho h}\right\rangle_{L_{2}(\Omega)}=\left\langle
y_{d},\hat{y}_{\varrho h}\right\rangle_{L_{2}(\Omega)}.$

Rewriting this equality gives

$\displaystyle\left\langle\hat{p}_{\varrho h},\hat{p}_{\varrho
h}\right\rangle_{h}+\left\langle\hat{y}_{\varrho h}-y_{d},\hat{y}_{\varrho
h}-y_{d}\right\rangle_{L_{2}(\Omega)}=\left\langle y_{d}-\hat{y}_{\varrho
h},y_{d}\right\rangle_{L_{2}(\Omega)},$

which yields the desired estimate. ∎

######  Theorem 3.

Let $(\hat{y}_{\varrho h},\hat{p}_{\varrho h})\in V_{h}\times V_{h}$ be the
unique solution of the coupled variational formulation (18)-(19). For
$y_{d}\in H_{0}^{s}(\Omega)$, $s\in[0,1]$, and $y_{d}\in H^{s}(\Omega)\cap
H_{0}^{1}(\Omega)$, $s\in(1,2]$, there holds the error estimate

$\displaystyle\|\hat{y}_{\varrho
h}-y_{d}\|_{L_{2}(\Omega)}\,\leq\,c\,h^{s}\,\|y_{d}\|_{H^{s}(\Omega)}.$ (26)

###### Proof.

This is a direct consequence of Theorem 2 and Lemma 4, together with a space
interpolation argument. ∎

##  3 Nested PCG Iteration

Finally, we have to solve the spd mass-lumped Schur-complement system (5) that
we now write in the compact from: find $\mathbf{y}_{h}\in\mathbb{R}^{n_{h}}$
such that

$S_{h}\mathbf{y}_{h}=\mathbf{y}_{dh},$ (27)

where $S_{h}=\varrho K_{h}D_{h}^{-1}K_{h}+M_{h}$, and $D_{h}$ is the lumped
mass matrix $\text{lump}(M_{h})$. For simplicity, we omit the hat over
$\mathbf{y}_{h}$ in (27) and throughout this section. The fast solution of the
symmetric and indefinite system (3) with the original mass matrix $M_{h}$
instead of $D_{h}$ was studied in [23]. Since the matrix $D_{h}$ is diagonal,
the matrix-by-vector multiplication $S_{h}*\mathbf{y}_{h}$ can now be
performed efficiently. Therefore, we can use the PCG method for solving (27).
Moreover, it turns out that $D_{h}$ can also serve as preconditioner in the
case $\varrho=h^{4}$ that leads to the optimally balanced estimate of
$\|{\hat{y}}_{\varrho h}-y_{d}\|_{L_{2}(\Omega)}$ as was shown in Section 2;
see Theorem 2 and Theorem 3. First of all, we can easily show that $D_{h}$ is
spectrally equivalent to $M_{h}$, i.e., there exist positive, $h$-independent
constants $\underline{c}_{\text{\tiny MD}}$ and $\overline{c}_{\text{\tiny
MD}}$ such that

$\underline{c}_{\text{\tiny MD}}D_{h}\leq M_{h}\leq\overline{c}_{\text{\tiny
MD}}D_{h},$ (28)

where $\underline{c}_{\text{\tiny MD}}=\lambda_{\text{\tiny
min}}=\lambda_{\text{\tiny min}}(D_{\tau}^{-1}M_{\tau})=1/(d+2)$ and
$\overline{c}_{\text{\tiny MD}}=\lambda_{\text{\tiny
max}}=\lambda_{\text{\tiny max}}(D_{\tau}^{-1}M_{\tau})=1$ are the minimal
eigenvalue and maximal eigenvalue of the small generalized eigenvalue problem
$M_{\tau}\mathbf{v}_{\tau}=\lambda D_{\tau}\mathbf{v}_{\tau}$ in
$\mathbb{R}^{d+1}$, respectively. $M_{\tau}$ and $D_{\tau}$ denote resp. the
mass matrix and the lumped mass matrix corresponding to the reference element
(unit simplex) $\tau$ to which every element $\tau_{e}$ from $\mathcal{T}_{h}$
is mapped by an affine-linear mapping $x=x_{e_{1}}+J_{e}\eta$. We note that
the spectral equivalence inequalities (28) are nothing but the algebraic
version of the inequalities (20) where the constants were already explicitly
computed.

In order to estimate the Schur complement $S_{h}=\varrho
K_{h}D_{h}^{-1}K_{h}+M_{h}$ by the mass matrix $M_{h}$ from below and above in
the spectral sense, it is obviously enough to estimate $\varrho
K_{h}D_{h}^{-1}K_{h}$ from above by $M_{h}$. Using the spectral equivalence
inequalities (28), Cauchy’s inequality, local inverse inequalities, and
$\varrho=h^{4}$, we get

$\displaystyle(\varrho K_{h}D_{h}^{-1}K_{h}\mathbf{v}_{h},\mathbf{v}_{h})$
$\displaystyle=$
$\displaystyle\varrho\,(D_{h}^{-1}K_{h}\mathbf{v}_{h},K_{h}\mathbf{v}_{h})\,\leq\,\overline{c}_{\text{\tiny
MD}}\,\varrho\,(M_{h}^{-1}K_{h}\mathbf{v}_{h},K_{h}\mathbf{v}_{h})$ (29)
$\displaystyle=$ $\displaystyle\overline{c}_{\text{\tiny
MD}}(K_{h}(\varrho^{-1}M_{h})^{-1}K_{h}\mathbf{v}_{h},\mathbf{v}_{h})$
$\displaystyle=$ $\displaystyle\overline{c}_{\text{\tiny
MD}}\sup_{\mathbf{q}_{h}\in\mathbb{R}^{n_{h}}}\frac{(K_{h}\mathbf{v}_{h},\mathbf{q}_{h})^{2}}{(\varrho^{-1}M_{h}\mathbf{q}_{h},\mathbf{q}_{h})}$
$\displaystyle=$ $\displaystyle\overline{c}_{\text{\tiny MD}}\sup_{q_{h}\in
V_{h}}\frac{\displaystyle\left[\int_{\Omega}\varrho^{1/4}\nabla
v_{h}\cdot\varrho^{-1/4}\nabla
q_{h}dx\right]^{2}}{\displaystyle\int_{\Omega}\varrho^{-1}[q_{h}(x)]^{2}dx}$
$\displaystyle\leq$ $\displaystyle\overline{c}_{\text{\tiny MD}}\sup_{q_{h}\in
V_{h}}\frac{\|\varrho^{1/4}\nabla
v_{h}\|_{L_{2}(\Omega)}^{2}\|\varrho^{-1/4}\nabla
q_{h}\|_{L_{2}(\Omega)}^{2}}{\displaystyle\int_{\Omega}\varrho^{-1}[q_{h}(x)]^{2}dx}$
$\displaystyle=$ $\displaystyle\overline{c}_{\text{\tiny
MD}}\,\|\varrho^{1/4}\nabla v_{h}\|_{L_{2}(\Omega)}^{2}\sup_{q_{h}\in
V_{h}}\frac{\displaystyle\sum\limits_{\tau_{e}\in\mathcal{T}_{h}}h^{-2}\int_{\tau_{e}}|\nabla
q_{h}|^{2}dx}{\displaystyle\int_{\Omega}\varrho^{-1}[q_{h}(x)]^{2}dx}$
$\displaystyle\leq$ $\displaystyle\overline{c}_{\text{\tiny
MD}}\,\|\varrho^{1/4}\nabla v_{h}\|_{L_{2}(\Omega)}^{2}\sup_{q_{h}\in
V_{h}}\frac{\displaystyle\sum\limits_{\tau_{e}\in\mathcal{T}_{h}}h^{-4}\,c_{\text{\tiny
inv}}^{2}\int_{\tau_{e}}(q_{h})^{2}\,dx}{\displaystyle\sum_{\tau_{e}\in\mathcal{T}_{h}}h^{-4}\int_{\tau_{e}}(q_{h})^{2}dx}$
$\displaystyle=$ $\displaystyle c_{\text{\tiny
inv}}^{2}\,\overline{c}_{\text{\tiny MD}}\,\|\varrho^{1/4}\nabla
v_{h}\|_{L_{2}(\Omega)}^{2}$ $\displaystyle=$ $\displaystyle c_{\text{\tiny
inv}}^{2}\,\overline{c}_{\text{\tiny
MD}}\,\sum\limits_{\tau_{e}\in\mathcal{T}_{h}}h^{2}\int_{\tau_{e}}|\nabla
v_{h}|^{2}dx$ $\displaystyle\leq$ $\displaystyle c_{\text{\tiny
inv}}^{4}\,\overline{c}_{\text{\tiny
MD}}\,\sum\limits_{\tau_{e}\in\mathcal{T}_{h}}\int_{\tau_{e}}(v_{h})^{2}dx$
$\displaystyle=$ $\displaystyle c_{\text{\tiny
inv}}^{4}\,\overline{c}_{\text{\tiny
MD}}\,(M_{h}\mathbf{v}_{h},\mathbf{v}_{h}),\;\;\forall\mathbf{v}_{h}\in
V_{h},$

where $c_{\text{\tiny inv}}$ is the universal positive constant in the local
inverse inequalities

$\|\nabla w_{h}\|_{L_{2}(\tau_{e})}\leq c_{\text{\tiny
inv}}\,h_{e}^{-1}\,\|w_{h}\|_{L_{2}(\tau_{e})}\quad\forall w_{h}\in
V_{h},\;\forall\tau_{e}\in\mathcal{T}_{h}.$ (30)

Here the local mesh size $h_{e}$ can be replaced by the global mesh size $h$
since we assumed quasi-uniform and shape-regular mesh $\mathcal{T}_{h}$. The
local inverse inequalities (30) can again be proved by mapping $\tau_{e}$ to
the unit simplex $\tau$. In this way the constant $c_{\text{\tiny inv}}$ can
even be computed explicitly in dependence of the mesh characteristics [12].
Therefore, we have just proved the spectral equivalence theorem that is
fundamental for the efficient solution of the spd mass-lumped Schur-complement
system (27) by means of PCG iteration.

######  Theorem 4.

Let us assume that the mesh $\mathcal{T}_{h}$ is globally quasi-uniform with
the global mesh-size $h$, and $\varrho=h^{4}$. Then the spectral equivalence
inequalities

$\underline{c}_{\text{\tiny SD}}D_{h}\leq\underline{c}_{\text{\tiny
SM}}M_{h}\leq S_{h}=\varrho
K_{h}D_{h}^{-1}K_{h}+M_{h}\leq\overline{c}_{\text{\tiny
SM}}M_{h}\leq\overline{c}_{\text{\tiny SD}}D_{h},$ (31)

hold with the spectral equivalence constants

$\underline{c}_{\text{\tiny SM}}=1,\quad\overline{c}_{\text{\tiny
SM}}=c_{\text{\tiny inv}}^{4}\overline{c}_{\text{\tiny
MD}}+1,\quad\underline{c}_{\text{\tiny SD}}=\underline{c}_{\text{\tiny
MD}}=\lambda_{\text{\tiny min}}=\lambda_{\text{\tiny
min}}(D_{T}^{-1}M_{T})=\frac{1}{d+2},$

and $\overline{c}_{\text{\tiny SD}}=\overline{c}_{\text{\tiny
MD}}^{2}c_{\text{\tiny inv}}^{4}+\overline{c}_{\text{\tiny MD}}$, where
$\overline{c}_{\text{\tiny MD}}=\lambda_{\text{\tiny
max}}=\lambda_{\text{\tiny max}}(D_{T}^{-1}M_{T})=1$.

###### Proof.

The spectral equivalence inequalities (31) immediately follow from the
inequalities (28), and (29). ∎

######  Remark 1.

The spectral estimate (29) can also be proved by Fourier analysis when one
expands the vectors $\mathbf{v}_{h}$ into the orthonormal eigenvector basis
corresponding to the eigenvalue problem $K_{h}\mathbf{e}_{h}=\lambda
D_{h}\mathbf{e}_{h}$ as it was done in [23] for $D_{h}=M_{h}$. In [22], we
provide a rigorous analysis of the variable $L_{2}$ regularization with a
technique that is different from the technique used for proving (29). We note
that the latter technique can be used to analyse the case of constant and
variable energy regularizations for state equations leading to non-symmetric
fe stiffness matrices $K_{h}$ such as convection-diffusion problems as well as
parabolic and hyperbolic problems when using space-time finite element
discretizations; see [27, 30].

Now we can efficiently solve the mass-lumped Schur-complement system (5)
respectively (27) by means of the PCG methods because, thanks to mass lumping,
the matrix-vector multiplication $S_{h}*\mathbf{y}_{h}^{k}$ can be performed
in asymptotically optimal complexity $O(n_{h})$, and, at the same time, the
lumped mass matrix $D_{h}=\text{lump}(M_{h})$ is a perfect preconditioner.
More precisely, let $\mathbf{y}_{h}^{k}\in\mathbb{R}^{n_{h}}$ be the $k$th PCG
iterate. Due to the spectral equivalence inequalities (28) and (31), and the
well-known convergence rate estimate for the PCG method (see, e.g., [39,
Chapter 13]), we can estimate the $L_{2}$ error $\|\hat{y}_{\varrho
h}-y_{\varrho h}^{k}\|_{L_{2}(\Omega)}$ between the fe functions
$\hat{y}_{\varrho h}(x)=\sum_{i=1}^{n_{h}}y_{i}\varphi_{i}^{h}(x)\in V_{h}$
and $y_{\varrho h}^{k}(x)=\sum_{i=1}^{n_{h}}y_{i}^{k}\varphi_{i}^{h}(x)\in
V_{h}$ corresponding to the solution
$\mathbf{y}_{h}=(y_{i})_{i=1,\ldots,n_{h}}\in\mathbb{R}^{n_{h}}$ of the Schur
complement system (27) and the $k$-th PCG iterate
$\mathbf{y}_{h}^{k}=(y_{i}^{k})_{i=1,\ldots,n_{h}}\in\mathbb{R}^{n_{h}}$,
respectively, as follows:

$\displaystyle\|\hat{y}_{\varrho h}-y_{\varrho h}^{k}\|_{L_{2}(\Omega)}$
$\displaystyle=$
$\displaystyle\|\mathbf{y}_{h}-\mathbf{y}_{h}^{k}\|_{\mathbf{M}_{h}}:=(\mathbf{M}_{h}(\mathbf{y}_{h}-\mathbf{y}_{h}^{k}),\mathbf{y}_{h}-\mathbf{y}_{h}^{k})^{1/2}$
(32) $\displaystyle\leq$
$\displaystyle(\mathbf{S}_{h}(\mathbf{y}_{h}-\mathbf{y}_{h}^{k}),\mathbf{y}_{h}-\mathbf{y}_{h}^{k})^{1/2}$
$\displaystyle=$
$\displaystyle\|\mathbf{y}_{h}-\mathbf{y}_{h}^{k}\|_{\mathbf{S}_{h}}\leq\,2\,q^{k}\,\|\mathbf{y}_{h}-\mathbf{y}_{h}^{0}\|_{\mathbf{S}_{h}}$
$\displaystyle\leq$ $\displaystyle 2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,q^{k}\,\|\mathbf{y}_{h}-\mathbf{y}_{h}^{0}\|_{\mathbf{M}_{h}}=2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,q^{k}\,\|\hat{y}_{\varrho h}-y_{\varrho h}^{0}\|_{L_{2}(\Omega)},$

where
$q=(\sqrt{\text{cond}_{2}(\mathbf{D}_{h}^{-1}\mathbf{S}_{h})}-1)/(\sqrt{\text{cond}_{2}(\mathbf{D}_{h}^{-1}\mathbf{S}_{h})}+1)<1$,
and
$\text{cond}_{2}(\mathbf{D}_{h}^{-1}\mathbf{S}_{h})=\lambda_{\text{max}}(\mathbf{D}_{h}^{-1}\mathbf{S}_{h})/\lambda_{\text{min}}(\mathbf{D}_{h}^{-1}\mathbf{S}_{h})$
denotes the spectral condition number of $\mathbf{D}_{h}^{-1}\mathbf{S}_{h}$
that can be bounded by the constant

$\frac{\overline{c}_{\text{\tiny SD}}}{\underline{c}_{\text{\tiny
SD}}}=\frac{\overline{c}_{\text{\tiny MD}}^{2}c_{\text{\tiny
inv}}^{4}+\overline{c}_{\text{\tiny MD}}}{\underline{c}_{\text{\tiny
MD}}}=\frac{\lambda_{\text{\tiny max}}(D_{T}^{-1}M_{T})^{2}c_{\text{\tiny
inv}}^{4}+\lambda_{\text{\tiny max}}(D_{T}^{-1}M_{T})}{\lambda_{\text{\tiny
min}}(D_{T}^{-1}M_{T})}=(d+2)(c_{\text{\tiny inv}}^{4}+1)$

that is independent of $h$. Using the triangle inequality, the $L_{2}$-norm
discretization error estimate (26) from Theorem 3, the $L_{2}$-norm iteration
error estimate (32), the inequality

$\|\hat{y}_{\varrho h}\|_{L_{2}(\Omega)}\leq\|y_{d}\|_{L_{2}(\Omega)}$ (33)

that follow from (18)-(19) when we choose the test functions
$q_{h}=\hat{p}_{\varrho h}$ and $v_{h}=\hat{y}_{\varrho h}$, we finally arrive
at $L_{2}$-norm estimate between desired state $y_{d}$ and the $k$th PCG
iterate $y_{\varrho h}^{k}$ computed by the PCG method:

$\displaystyle\|y_{d}-y_{\varrho h}^{k}\|_{L_{2}(\Omega)}$ $\displaystyle\leq$
$\displaystyle\|y_{d}-\hat{y}_{\varrho h}\|_{L_{2}(\Omega)}+\|\hat{y}_{\varrho
h}-y_{\varrho h}^{k}\|_{L_{2}(\Omega)}$ (34) $\displaystyle\leq$
$\displaystyle
c\,h^{s}\,\|y_{d}\|_{H^{s}(\Omega)}+2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,q^{k}\,\|\hat{y}_{\varrho h}-y_{\varrho h}^{0}\|_{L_{2}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
h^{s}\Big{(}c\,\|y_{d}\|_{H^{s}(\Omega)}+2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,\|\hat{y}_{\varrho h}-y_{\varrho h}^{0}\|_{L_{2}(\Omega)}\Big{)}$
$\displaystyle\leq$ $\displaystyle
h^{s}\Big{(}c\,\|y_{d}\|_{H^{s}(\Omega)}+2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,\|y_{d}\|_{L_{2}(\Omega)}\Big{)}=h^{s}c(y_{d})$

provided that $q^{k}\leq h^{s}$ and that the initial guess $y_{\varrho h}^{0}$
is chosen to be zero. Therefore, $k\geq\ln h^{-s}/\ln q^{-1}$ ensures that the
PCG computes an approximation $y_{\varrho h}^{k}$ to the desired state $y_{d}$
that differs from $y_{d}$ in the same order $O(h^{s})$ as the discretization
error $\|y_{d}-\hat{y}_{\varrho h}\|_{L_{2}(\Omega)}$ in the $L_{2}$ norm.
Moreover, this can be done with $O(n_{h}\ln h^{-1})=O(h^{-d}\ln h^{-1})$
arithmetical operations, i.e. the complexity is asymptotically optimal up to
the logarithmical factor $\ln h^{-1}$.

This logarithmical factor can be avoided in a nested iteration setting on a
sequence of refined (nested) meshes. Indeed, let us consider a sequence of
uniformly refined meshes $\mathcal{T}_{\ell}=\mathcal{T}_{h_{\ell}}$ with the
mesh size $h_{\ell}$ and the optimally balanced regularization parameter
$\varrho_{\ell}=h_{\ell}^{4}$, $\ell=1,\ldots,L$, where
$h_{\ell}=h_{\ell-1}/2$, $\ell=2,\ldots,L$. Thus, the coarsest mesh
corresponds to the subindex $1$, whereas the finest mesh is related to $L$. On
every mesh $\mathcal{T}_{\ell}$, $\ell=1,\ldots,L$, we have to solve the mass-
lumped Schur-complement system (27) that we now write in form: find
$\mathbf{y}_{\ell}=\mathbf{y}_{h_{\ell}}\in\mathbb{R}^{n_{\ell}}=\mathbb{R}^{n_{h_{\ell}}}$
such that

$S_{\ell}\mathbf{y}_{\ell}=\mathbf{y}_{d\ell}$ (35)

where $S_{\ell}=\varrho_{\ell}K_{\ell}D_{\ell}^{-1}K_{\ell}+M_{\ell}$,
$K_{\ell}=K_{h_{\ell}}$. $D_{\ell}=D_{h_{\ell}}$, $M_{\ell}=M_{h_{\ell}}$,
$\mathbf{y}_{d\ell}=\mathbf{y}_{dh_{\ell}}$, and
$\varrho_{\ell}=h_{\ell}^{4}$.

Now the nested iteration algorithm works as follows. First we solve the
coarse-mesh problem (35), $\ell=1$, sufficiently accurate. More precisely, we
compute an iterate $\mathbf{y}_{1}^{k_{1}}\in\mathbb{R}^{n_{1}}$ corresponding
to the fe function $y_{1}^{k_{1}}=y_{\varrho_{1}h_{1}}^{k_{1}}\in
V_{1}=V_{h_{1}}$ (short: $\mathbf{y}_{1}^{k_{1}}\leftrightarrow
y_{1}^{k_{1}}$) such that

$\|y_{d}-y_{1}^{k_{1}}\|_{L_{2}(\Omega)}\leq h_{1}^{s}c(y_{d}).$ (36)

Due to (34), this can be done with $k_{1}\geq\ln h_{1}^{-s}/\ln q^{-1}$ PCG
iterations starting with $y_{1}^{0}=0$. Now, let us assume that, on level
$\ell-1\in\\{1,\ldots,L-1\\}$, the iterate $y_{\ell-1}^{k_{\ell-1}}\in
V_{\ell-1}$ fulfills the estimate

$\|y_{d}-y_{\ell-1}^{k_{\ell-1}}\|_{L_{2}(\Omega)}\leq
h_{\ell-1}^{s}c(y_{d}).$ (37)

Let $y_{\ell}^{0}=I_{\ell-1}^{\ell}y_{\ell-1}^{k_{\ell-1}}$ be the affine-
linear interpolate of $y_{\ell-1}^{k_{\ell-1}}$. We note that
$y_{\ell}^{0}=I_{\ell-1}^{\ell}y_{\ell-1}^{k_{\ell-1}}=y_{\ell-1}^{k_{\ell-1}}\in
V_{\ell-1}\subset V_{\ell}$ since the meshes are nested. Then we get the
estimate

$\displaystyle\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\|y_{d}-\hat{y}_{\ell}\|_{L_{2}(\Omega)}+\|\hat{y}_{\ell}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$
(38) $\displaystyle\leq$ $\displaystyle
ch_{\ell}^{s}\|y_{d}\|_{H^{s}(\Omega)}+2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,q^{k_{\ell}}\,\|\hat{y}_{\ell}-y_{\ell}^{0}\|_{L_{2}(\Omega)},$

where $\hat{y}_{\ell}=\hat{y}_{\varrho_{\ell}h_{\ell}}\in V_{\ell}$ is the
exact state solution of the finite element scheme (18)-(19) corresponding to
the solution $\mathbf{y}_{\ell}\in\mathbb{R}^{n_{\ell}}$ of the mass-lumped
Schur-complement system (35). Now using
$y_{\ell}^{0}=I_{\ell-1}^{\ell}y_{\ell-1}^{k_{\ell-1}}$, the triangle
inequality, Theorem 3, and estimate (37), we can continue to estimate the last
term in (38) as follows:

$\displaystyle\|\hat{y}_{\ell}-y_{\ell}^{0}\|_{L_{2}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\|\hat{y}_{\ell}-y_{d}\|_{L_{2}(\Omega)}+\|y_{d}-I_{\ell-1}^{\ell}y_{\ell-1}^{k_{l-1}}\|_{L_{2}(\Omega)}$
(39) $\displaystyle\leq$
$\displaystyle\|\hat{y}_{\ell}-y_{d}\|_{L_{2}(\Omega)}+\|y_{d}-y_{\ell-1}^{k_{l-1}}\|_{L_{2}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
ch_{\ell}^{s}\|y_{d}\|_{H^{s}(\Omega)}+h_{\ell-1}^{s}c(y_{d})$
$\displaystyle\leq$ $\displaystyle
h_{\ell}^{s}(c\|y_{d}\|_{H^{s}(\Omega)}+2^{s}c(y_{d}))$

Inserting (39) into (38), we get

$\displaystyle\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
h_{\ell}^{s}\left[c\|y_{d}\|_{H^{s}(\Omega)}+2\,\overline{c}_{\text{\tiny
SM}}^{1/2}\,q^{k_{\ell}}(c\|y_{d}\|_{H^{s}(\Omega)}+2^{s}c(y_{d}))\right]$
(40) $\displaystyle\leq$ $\displaystyle h_{\ell}^{s}c(y_{d})$

provided that
$q^{k_{\ell}}(c\|y_{d}\|_{H^{s}(\Omega)}+2^{s}c(y_{d}))\leq\|y_{d}\|_{L_{2}(\Omega)}$.
The latter inequality is ensured if we performed not more than

$k_{\ell}=k_{*}\geq\ln(q(y_{d})^{-1})/\ln q^{-1}$ (41)

nested iterations, where
$q(y_{d})=\|y_{d}\|_{L_{2}(\Omega)}/((1+2^{s})c\|y_{d}\|_{H^{s}(\Omega)}+2^{1+s}\|y_{d}\|_{L_{2}(\Omega)})<1.$
Here we exclude the trivial case that $y_{d}=0$. Therefore, we have proved the
following nested iteration theorem by induction.

######  Theorem 5.

If the coarse mesh problem on level $l=1$ is solved by $k_{1}$ PCG iterations
with the initial guess $\mathbf{y}_{1}^{0}=\mathbf{0}_{1}$ such that (38)
holds, and if $k_{*}$ nested PCG iterations are used on all nested levels
$\ell=2,\ldots,L$, i.e. $k_{2}=\ldots=k_{L}=k_{*}$ defined by (41), then the
last iterate $y_{L}^{k_{L}}\leftrightarrow\mathbf{y}_{L}^{k_{L}}$ on the
finest level $\ell=L$ differs from given desired state $y_{d}$ in the order of
the discretization error $O(h_{L}^{s})$ with respect to the $L_{2}(\Omega)$
norm. More precisely, we get the estimate

$\|y_{d}-y_{L}^{k_{L}}\|_{L_{2}(\Omega)}\leq h_{L}^{s}c(y_{d}).$ (42)

The computation of $y_{L}^{k_{L}}\leftrightarrow\mathbf{y}_{L}^{k_{L}}$
requires not more than $O(n_{L})=O(h_{L}^{-d})$ arithmetical operations and
memory, i.e. the nested iteration procedure proposed is asymptotically
optimal.

###### Proof.

The proof of estimate (42) follows from above by induction. The complexity
analysis is based on simple use of the geometric series. ∎

We stop the nested iteration process as soon as we arrive at some desired
relative accuracy $\varepsilon\in(0,1)$ such that

$\|y_{d}-y_{L}^{k_{L}}\|_{L_{2}(\Omega)}\leq\varepsilon\,\|y_{d}\|_{L_{2}(\Omega)}.$
(43)

The apriori estimate (42) immediately yields that estimate (43) is guaranteed
when
$ch_{L}^{s}\|y_{d}\|_{H^{s}(\Omega)}\leq\varepsilon\,\|y_{d}\|_{L_{2}(\Omega)}$,
but in practice we directly check (43) because all quantities are computable.

We will also stop the nested iteration if the cost for the control
$u_{\ell}^{k_{\ell}}$ becomes too large, where
$u_{\ell}^{k_{\ell}}\leftrightarrow\mathbf{u}_{\ell}^{k_{\ell}}$ is computed
from the fe state equation

$\mathbf{u}_{\ell}^{k_{\ell}}=-\varrho_{\ell}^{-1}\mathbf{p}_{\ell}^{k_{\ell}}=D_{\ell}^{-1}K_{\ell}\mathbf{y}_{\ell}^{k_{\ell}}$

More precisely, let $c_{u}>0$ be a given threshold for the control cost that
we are willing to pay. Then we stop the nested iteration if

$\|u_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}=(M_{\ell}\mathbf{u}_{\ell}^{k_{\ell}},\mathbf{u}_{\ell}^{k_{\ell}})\leq\overline{c}_{\text{\tiny
MD}}(D_{\ell}\mathbf{u}_{\ell}^{k_{\ell}},\mathbf{u}_{\ell}^{k_{\ell}})=\overline{c}_{\text{\tiny
MD}}(K_{\ell}\mathbf{y}_{\ell}^{k_{\ell}},\mathbf{y}_{\ell}^{k_{\ell}})\leq
c_{u},$

but $\|u_{\ell+1}^{k_{\ell+1}}\|_{L_{2}(\Omega)}>c_{u}$, where
$\overline{c}_{\text{\tiny MD}}=1$. Then we set $L=\ell$.

Now we may proceed with cascadic nested iteration freezing the cost
(regularization) parameter $\varrho_{L}=h_{L}^{4}$ and refining the mesh only
, i.e.

$\varrho_{\ell+1}=\varrho_{L}=h_{L}^{4}=\text{const.}\;\;\mbox{and}\;\;h_{\ell+1}=h_{\ell}/2\;\mbox{for}\;\ell=L,\ldots,L+J-1,$
(44)

in order to improve the approximation of the control. We note that this
further mesh refinement will not improve the approximation to the desired
state $y_{d}$ since this error is defined by the frozen cost parameter
$\varrho_{L}$. Since we only use a few additional levels for the improvement
of the control, we can proceed with the PCG preconditioned by $D_{\ell}$ as
before as nested iteration, but replacing $\overline{c}_{\text{\tiny SD}}$ by
$\overline{c}_{\text{\tiny SD},\ell}=\overline{c}_{\text{\tiny
SD}}2^{4(\ell-L)}$ for $\ell=L+1,\ldots,L+J$. If we want to add many levels,
i.e. $J>>1$, then we may use some cascadic full multigrid Schur complement
iteration using level $L$ as coarse mesh and
$y_{L+1}^{0}=I_{L}^{L+1}y_{L}^{k_{L}}=y_{L}^{k_{L}}\in V_{L}\subset V_{L+1}$
as initial guess; see [18, 11] for $L_{2}$ convergent multigrid methods.

##  4 Numerical Results

In our numerical experiments, we consider the discontinuous desired state

$y_{d}=\begin{cases}1&\textup{ in }\;(0.25,0.75)^{3},\\\ 0&\textup{ in
}\;\overline{\Omega}\setminus(0.25,0.75)^{3},\end{cases}$

in the computational domain $\Omega=(0,1)^{3}\subset\mathbb{R}^{d=3}$. This
discontinuous desired state $y_{d}$ does not belong to $Y=H^{1}_{0}(\Omega)$,
and has a rather low Sobolev regularity. More precisely, $y_{d}\in
H^{1/2-\varepsilon}(\Omega)$ for any $\varepsilon>0$. This discontinuous
target has been utilized in the work [23, 26, 32] in both the cases of $L_{2}$
and energy ($H^{-1}$) regularization for distributed elliptic optimal control
problems. So, we can easily compare the numerical results presented below for
the mass-lumping discretization of the control term in the reduced optimality
system with those of the $L_{2}$ regularization without mass lumping and the
$H^{-1}$ regularization.

We decompose the domain $\Omega=(0,1)^{3}$ into uniformly refined tetrahedral
elements $\tau_{e}$, and start with an initial mesh that contains $384$
tetrahedral elements and $125$ vertices, leading to the mesh size $h=2^{-2}$.
From such a mesh, we make successive refinements on the levels $\ell=1,...,8$.
On the finest level $\ell=L=8$, we have $135,005,697$ vertices,
$h=2^{-9}=1.9531$e$-3$, and $\varrho=h^{4}=2^{-36}=1.4552$e$-11$. Further, we
run tests on the adaptively refined meshes, in which we have employed the
standard red-green refinement of tetrahedral elements, and we have chosen the
locally varying regularization parameter $\varrho_{\tau}=h_{e}^{4}$ on each
tetrahedral element $\tau_{e}$. The adaptive procedure is simply based on the
localization of the error $\|y_{d}-{\tilde{y}}_{\ell}\|_{L_{2}(\Omega)}$ that
is explicitly computable for any known fe approximation ${\tilde{y}}_{\ell}$
to the given desired state $y_{d}$; see [21] for a detailed description.

As described in Section 3, thanks to the replacement of the mass matrix
$M_{h}$ by its diagonal approximation $D_{h}=\text{lump}(M_{h})$, we can
efficiently solve the spd mass-lumped Schur-complement system (27) by means of
the PCG preconditioned by $D_{h}$. We first use the initial guess
$\mathbf{y}_{\ell}^{0}=\mathbf{0}$, and terminate the iteration as soon as the
preconditioned residual is reduced by a factor $10^{6}$. The number of PCG
iterations (Its)and the computational time (Time) in seconds (s) are provided
in Table 1 for both uniform and adaptive refinements.

Therein, we observe the robustness of our proposed preconditioner for (27)
with respect to both the mesh size and local adaptivity under the choice of
$\varrho_{\tau}=h_{e}^{4}$. We only see slightly more iterations for the
adaptive refinements in comparison to the uniform refinements.

$\ell$ | Adaptive | Uniform
---|---|---
#Dofs | error | Its (Time) | #Dofs | error | Its (Time)
$1$ | $125$ |  $3.26$e$-1$ |  $10$ ($6.3$e$-4$) | $125$ |  $3.26$e$-1$ |  $10$ ($6.4$e$-4$)
$2$ | $223$ |  $2.35$e$-1$ |  $62$ ($6.5$e$-3$) | $729$ |  $2.25$e$-1$ |  $55$ ($1.8$e$-2$)
$3$ | $1,044$ |  $1.86$e$-1$ |  $106$ ($5.5$e$-2$) | $4,913$ |  $1.59$e$-1$ |  $79$ ($1.9$e$-1$)
$4$ | $4,548$ |  $1.32$e$-1$ |  $123$ ($2.8$e$-1$) | $35,937$ |  $1.12$e$-1$ |  $85$ ($1.7$e$-0$)
$5$ | $10,524$ |  $1.05$e$-1$ |  $116$ ($6.3$e$-1$) | $274,625$ |  $7.96$e$-2$ |  $81$ ($2.3$e$+1$)
$6$ | $25,807$ |  $8.35$e$-2$ |  $113$ ($1.6$e$-0$) | $2,146,689$ |  $5.62$e$-2$ |  $74$ ($1.8$e$+2$)
$7$ | $91,520$ |  $6.03$e$-2$ |  $100$ ($5.7$e$-0$) | $16,974,593$ |  $3.97$e$-2$ |  $68$ ($1.4$e$+3$)
$8$ | $118,334$ |  $5.62$e$-2$ |  $102$ ($7.7$e$-0$) | $135,005,697$ |  $2.81$e$-2$ |  $66$ ($1.3$e$+4$)
$9$ | $432,195$ |  $4.08$e$-2$ |  $93$ ($3.5$e$+1$) |  |  |
${10}$ | $473,638$ |  $3.97$e$-2$ |  $95$ ($6.3$e$+1$) |  |  |
${11}$ | $1,843,740$ |  $2.84$e$-2$ |  $91$ ($2.5$e$+2$) |  |  |
${12}$ | $1,937,983$ |  $2.79$e$-2$ |  $92$ ($3.2$e$+2$) |  |  |
${13}$ | $7,681,306$ |  $1.99$e$-2$ |  $91$ ($6.3$e$+2$) |  |  |
${14}$ | $7,922,574$ |  $1.96$e$-2$ |  $93$ ($1.0$e$+3$) |  |  |
${15}$ | $31,496,575$ |  $1.39$e$-2$ |  $83$ ($3.6$e$+3$) |  |  |
${16}$ | $32,000,845$ |  $1.38$e$-2$ |  $84$ ($5.4$e$+3$) |  |  |
${17}$ | $127,607,911$ |  $9.84$e$-3$ |  $68$ ($1.6$e$+4$) |  |  |
Table 1: Comparison of the PCG iterations (Its) and computational time (Time)
in seconds (indicated in the parentheses) for solving (27) for both adaptive
and uniform refinements using the non-nested iterations, where error =
$\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$.

As shown in the theoretical part, solving the Schur complement equation with
the lumped mass does not deteriorate the convergence of our finite element
approximation. This is confirmed in our numerical experiments. The comparison
of convergence on both uniform and adaptive refinements is given in Figure 1.
We observe the convergence rate $h^{0.5}$ for the uniform refinement as
predicted by Theorem3, and a much better convergence rate $h^{0.75}$ for the
adaptive refinements; see [21] for the case of variable energy regularization.
There one can also find an explanation of the convergence rate that can be
achieved via this adaptive procedure.

Figure 1: Comparison of the convergence history obtained from non-nested
iterations for uniform and adaptive refinements, where
$\|e\|=\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$.

In order to further reduce the computational cost, we utilize nested PCG
iterations as described in Section 3. Here, on the coarsest level $\ell=1$, we
run the PCG iterations until the relative preconditioned residual reaches
$10^{-6}$. On the refined levels $\ell=2,3,...$, we have utilized an adaptive
tolerance

$\alpha\left[n_{\ell}/n_{\ell-1}\right]^{-\frac{\beta}{3}},\;\ell=2,3,\dots,$
(45)

for the relative preconditioned residual, with $\alpha$ being a scaling
factor, $\beta=0.5$ and $0.75$ for the uniform and adaptive refinement,
respectively, and $n_{\ell}$ the number of degrees of freedom (#Dofs) on the
mesh level $\ell=1,2,\ldots$ . The solution on the level $\ell-1$ is used as
an initial guess for the PCG iteration on the next finer level $\ell$. The
reduced number of nested iterations (Its) on both uniform and adaptive
refinements is given in Table 2, where we have chosen $\alpha=0.5$ and
$\alpha=1$ for the adaptive and uniform refinement, respectively. From this,
we easily see much fewer iteration numbers and significantly less
computational time in seconds in comparison with the case of non-nested
iterations as shown in Table 1, without loss of accuracy of the numerical
approximations; see Figure 2 for a comparison of convergence history for both
uniform and adaptive refinements using the nested iterations.

$\ell$ | Adaptive | Uniform
---|---|---
#Dofs | error | Its (Time) | #Dofs | error | Its (Time)
$1$ | $125$ |  $3.26$e$-1$ |  $10$ ($6.3$e$-4$) | $125$ |  $3.26$e$-1$ |  $10$ ($6.3$e$-4$)
$2$ | $223$ |  $3.30$e$-1$ |  $1$ ($2.6$e$-4$) | $729$ |  $2.27$e$-1$ |  $8$ ($2.9$e$-3$)
$3$ | $1,067$ |  $1.84$e$-1$ |  $19$ ($1.1$e$-2$) | $4,913$ |  $2.25$e$-1$ |  $1$ ($4.6$e$-3$)
$4$ | $4,705$ |  $1,28$e$-1$ |  $13$ ($3.3$e$-2$) | $35,937$ |  $1.08$e$-1$ |  $9$ ($1.9$e$-1$)
$5$ | $15,368$ |  $1.00$e$-1$ |  $17$ ($1.4$e$-1$) | $274,625$ |  $8.22$e$-2$ |  $8$ ($1.5$e$-0$)
$6$ | $30,996$ |  $8.45$e$-2$ |  $17$ ($4.0$e$-1$) | $2,146,689$ |  $5.60$e$-2$ |  $9$ ($1.4$e$+1$)
$7$ | $94,176$ |  $6.30$e$-2$ |  $19$ ($1.3$e$-0$) | $16,974,593$ |  $3.98$e$-2$ |  $9$ ($2.1$e$+2$)
$8$ | $129,760$ |  $5.68$e$-2$ |  $18$ ($1.7$e$-0$) | $135,005,697$ |  $2.81$e$-2$ |  $9$ ($2.2$e$+3$)
$9$ | $440,572$ |  $4.18$e$-2$ |  $17$ ($1.2$e$+1$) |  |  |
${10}$ | $488,124$ |  $4.03$e$-2$ |  $17$ ($1.3$e$+1$) |  |  |
${11}$ | $1,860,339$ |  $2.90$e$-2$ |  $18$ ($6.1$e$+1$) |  |  |
${12}$ | $1,958,388$ |  $2.85$e$-2$ |  $16$ ($5.9$e$+1$) |  |  |
${13}$ | $7,254,384$ |  $2.06$e$-2$ |  $18$ ($2.6$e$+2$) |  |  |
${14}$ | $7,408,106$ |  $2.04$e$-2$ |  $16$ ($2.1$e$+2$) |  |  |
${15}$ | $29,094,073$ |  $1.47$e$-2$ |  $17$ ($6.9$e$+2$) |  |  |
${16}$ | $29,682,531$ |  $1.44$e$-2$ |  $16$ ($7.6$e$+2$) |  |  |
${17}$ | $116,229,104$ |  $1.04$e$-2$ |  $16$ ($3.7$e$+3$) |  |  |
Table 2: Comparison of the PCG iterations (Its) and computational time (Time)
in seconds (indicated in the parentheses) for solving (27) for both adaptive
($\alpha=0.5$, $\beta=0.75$) and uniform ($\alpha=1$, $\beta=0.5$) refinements
using the nested iteration approach, where error =
$\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$ Figure 2:  Comparison of the
convergence history obtained from nested iterations for uniform and adaptive
refinements, where $\|e\|=\|y_{d}-y_{\ell}^{k_{\ell}}\|_{L_{2}(\Omega)}$.

Another approach to reduce the computational time, especially, in the case of
uniform refinement is the parallelization of the PCG solver. The
parallelization of the conjugate gradient algorithm is now a standard
procedure [13]. The crucial point is always the preconditioner. It is clear
that the parallelization of a diagonal preconditioner is much easier than the
parallelization of a multigrid preconditioner. More precisely, the
parallelization of a diagonal preconditioner such as $D_{h}$ is trivial. For
parallel performance studies, we have utilized the open source
MFEM555https://mfem.org/. We observe from the diagonals of Table 3, e.g. from
level $7$ with $16$ cores to level $8$ with $512$ cores (always factor $8$),
almost constant time, i.e. a good weak scaling behavior, whereas the
horizontal lines show an almost perfect strong scaling. The latter one is also
illustrated in Figure 3 for $\ell=7$ and $\ell=8$ corresponding to
$16,974,593$ and $135,005,697$ Dofs, respectively. The largest problem with
$135,005,697$ Dofs can be solved in $6.4$ seconds using $512$ cores. Similar
scaling behaviors are also observed for the nested iterations approach; see
Table 4 and Figure 4. The computational time in seconds (s) using nested
iterations is further reduced by a factor of about $7$ in comparison with the
non-nested iterations. Using $512$ cores, the largest problem with
$135,005,697$ Dofs is solved in $1$ second. Finally, we made some performance
tests for the adaptive refinement using the nested iteration setting. The
results are given in Table 5. We observe relatively good scaling in this case
as well. Here, we have used the non-conforming simplicial complex and load
balance from the open source MFEM.

We note that we used different computers and different codes for the single-
core and parallel computations. More precisely, we used the shared-memory
computer MACH2666https://www3.risc.jku.at/projects/mach2/, that provides a big
memory, and the distributed-memory computer
RADON1777https://www.oeaw.ac.at/ricam/hpc for the single-core and parallel
computations, respectively.

$\ell$ | #Cores
---|---
16 | 32 | 64 | 128 | 256 | 512
4 |  $85$ ($4.0$e$-2$) | - | - | - | - | -
5 |  $84$ ($3.4$e$-1$) |  $84$ ($1.6$e$-1$) |  $84$ ($6.4$e$-2$) | - | - | -
6 |  $82$ ($2.9$e$-0$) |  $82$ ($1.4$e$-0$) |  $82$ ($7.3$e$-1$) |  $82$ ($3.7$e$-1$) |  $82$ ($1.8$e$-1$) |  $82$ ($8.0$e$-2$)
7 |  $80$ ($2.5$e$+1$) |  $80$ ($1.2$e$+1$) |  $80$ ($6.3$e$-0$) |  $80$ ($3.0$e$-0$) |  $80$ ($1.5$e$-0$) |  $80$ ($8.3$e$-1$)
8 | - | - |  $77$ ($5.1$e$+1$) |  $77$ ($2.5$e$+1$) |  $77$ ($1.3$e$+1$) |  $77$ ($6.4$e$-0$)
Table 3: Parallel performance on a distributed computer system for uniform refinement and non-nested iterations.  $\ell$ | #Cores
---|---
16 | 32 | 64 | 128 | 256 | 512
4 |  $9$ ($5.8$e$-3$) | - | - | - | - | -
5 |  $9$ ($5.0$e$-2$) |  $9$ ($2.7$e$-2$) |  $9$ ($1.1$e$-2$) | - | - | -
6 |  $8$ ($3.5$e$-1$) |  $8$ ($1.8$e$-1$) |  $8$ ($9.4$e$-2$) |  $8$ ($5.2$e$-2$) |  $8$ ($2.8$e$-2$) |  $8$ ($1.2$e$-2$)
7 |  $9$ ($3.1$e$-0$) |  $9$ ($1.6$e$-0$) |  $9$ ($8.1$e$-1$) |  $9$ ($4.2$e$-1$) |  $9$ ($2.2$e$-1$) |  $9$ ($1.2$e$-1$)
8 | - | - |  $11$ ($7.9$e$-0$) |  $11$ ($4.0$e$-0$) |  $11$ ($2.0$e$-0$) |  $11$ ($1.0$e$-0$)
Table 4: Parallel performance on a distributed computer system for uniform refinement and nested iterations. Figure 3: Strong scalability and computational time in seconds (s) with respect to the number of cores for uniform refinement and non-nested iterations Figure 4: Strong scalability and computational time in seconds (s) with respect to the number of cores for uniform refinement and nested iterations #Dofs | #Cores
---|---
16 | 32 | 64 | 128 | 256
$2.76154$e+$6$ |  $16$ ($1.0$e$-0$) |  $16$ ($5.3$e$-1$) |  16 ($2.8$e$-1$) |  16 ($1.6$e$-1$) |  16 ($1.0$e$-1$)
$1.06728$e+$7$ | - |  17 ($2.3$e$-0$) |  17 ($1.2$e$-0$) |  16 ($6.2$e$-1$) |  17 ($3.3$e$-1$)
Table 5: Parallel performance on a distributed computer system for adaptive
refinement and nested iterations.

##  5 Conclusions and Outlook

We provide a rigorous analysis of the discretization error
$\|y_{d}-{\hat{y}}_{\varrho h}\|_{L_{2}(\Omega)}$ when replacing the mass
matrix $M_{h}$ arising from the regularization term in the reduced optimality
system by its lumped version $D_{h}=\text{lump}(M_{h})$. It turns out that the
asymptotic behavior of the error is not affected by mass lumping when using
affine-linear finite elements. More precisely, we again get the upper bound
$ch^{s}\|y_{d}\|_{H^{s}(\Omega)}$, $s\in[0,2]$, for the choice $\varrho=h^{4}$
that provides the optimal balance between the regularization parameter
$\varrho$ and the mesh-size $h$. Moreover, this replacement of $M_{h}$ by
$D_{h}$ opens the way to reduce the discrete reduced optimal optimality system
further to a spd Schur complement problem that can efficiently be solved by
PCG since now the matrix-by-vector multiplication is cheap and, surprisingly,
$D_{h}$ is a diagonal preconditioner that is spectrally equivalent to the
Schur complement $S_{h}$. This PCG can efficiently be parallelized as the
numerical results show. These findings provide the perfect ingredients for a
nested PCG iteration producing iterates $y_{\ell}^{k_{\ell}}$ that differ from
the desired state $y_{d}$ in the order $O(h_{\ell}^{s})$ of the discretization
error in asymptotically optimal complexity $O(h_{\ell}^{-d})$. The nested
iteration process will be stopped when some relative accuracy
$\varepsilon\in(0,1)$ of the error is reached, or the cost we are willing to
pay in terms of the control energy density $\|u_{L}\|^{2}_{L_{2}(\Omega)}$
becomes too large. In this case, we can freeze the regularization (cost)
parameter $\varrho_{L}=h_{L}^{4}$, and continue the nested iteration process
with mesh refinement only in order to improve to the approximation of the
control.

We provide not only numerical results for the case of uniform refinement that
nicely demonstrated the theoretical predictions but also for adaptive
refinement when using variable regularization. The numerical results show that
this adaptive approach works well, but a rigorous numerical analysis is still
missing. Further investigation comprises this analysis, and the generalization
to larger classes of PDEs like elliptic diffusion-convection-reaction,
parabolic and hyperbolic state equations. We refer to [25, 24, 27] and [30]
when using space-time fe discretization for parabolic and hyperbolic initial-
boundary value problems, respectively. Another future research topic are the
consideration of control and state (box) constraints in the framework
discussed here; see [17] for first results. Finally, we mention that singular-
perturbed problems as discussed here also appear in fluid mechanics where they
are known as (discrete) differential filter that provide approximate
deconvolution models of turbulence [9, 28, 20].

## References

  * [1] O. Axelsson and J. Karátson.  Superior properties of the PRESB preconditioner for operators on two-by-two block form with square blocks.  Numer. Math, 146(2):335–368, 2020.
  * [2] O. Axelsson, M. Neytcheva, and A. Ström.  An efficient preconditioning method for state box-constrained optimal control problems.  J. Numer. Math., 26(4):185–207, 2018.
  * [3] Z.-Z. Bai.  Regularized HSS iteration methods for stabilized saddle-point problems.  IMA J. Numer. Anal., 39(4):1888–1923, 2019.
  * [4] Z.-Z. Bai and M. Benzi.  Regularized HSS iteration methods for saddle-point linear systems.  BIT Numer. Math., 57(2):287–311, 2017.
  * [5] Z.-Z. Bai, M. Benzi, F. Chen, and Z.-Q. Wang.  Preconditioned MHSS iteration methods for a class of block two-by-two linear systems with applications to distributed control problems.  IMA J. Numer. Anal., 33(1):343–369, 2013.
  * [6] Z.-Z. Bai and J.-Y. Pan.  Matrix Analysis and Computations.  SIAM, 2021.
  * [7] R. Becker and P. Hansbo.  A simple pressure stabilization method for the Stokes equation.  Comm. Numer. Methods Engrg., 24(11):1421–1430, 2008.
  * [8] M. Benzi, G. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * [9] L. Berselli, T. Iliescu, and W. Layton.  Mathematics of large eddy simulation of turbulent flows.  Scientific Computation. Springer-Verlag, Berlin, 2006.
  * [10] A. Borzi and V. Schulz.  Multigrid methods for PDE optimization.  SIAM Review, 51(2):361–395, 2009.
  * [11] D. Braess.  Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics.  Cambridge University Press, Cambridge, 2007.
  * [12] S. Chen and J. Zhao.  Estimations of the constants in inverse inequalities for finite element functions.  J. Comput. Math., 31(5):522–531, 2013.
  * [13] C. Douglas, G. Haase, and U. Langer.  A Tutorial on Elliptic PDE Solvers and Their Parallelization.  Software, Environments, and Tools,. SIAM, Philadelphia, 2003.
  * [14] I. Dravins and M. Neytcheva.  On the Numerical Solution of State- and Control-constrained Optimal Control Problems.  Department of Information Technology, Uppsala Universitet, 2021.
  * [15] H. C. Elman, D. J. Silvester, and A. J. Wathen.  Finite elements and fast iterative solvers: With applications in incompressible fluid dynamics.  Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2005.
  * [16] A. Ern and J.-L. Guermond.  Theory and Practice of Finite Elements.  Springer-Verlag, New York, 2004.
  * [17] P. Gangl, R. Löscher, and O. Steinbach.  Regularization and finite element error estimates for distributed optimal control problems with energy regularization and state or control constraints.  Technical report, TU Graz, 2023.  in preparation.
  * [18] W. Hackbusch.  Multi-Grid Methods and Applications.  Springer Verlag, 1985.
  * [19] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE Constraints, volume 23.  Springer-Verlag, Berlin, 2009.
  * [20] V. John.  Finite Element Methods for Incompressible Flow Problems, volume 51 of Springer Series in Computational Mathematics.  Springer, 2016.
  * [21] U. Langer, R. Löscher, O. Steinbach, and H. Yang.  An adaptive finite element method for distributed elliptic optimal control problems with variable energy regularization.  Technical Report arXiv:2209.08811, arXiv.org, 2022.
  * [22] U. Langer, R. Löscher, O. Steinbach, and H. Yang.  Robust iterative solvers for algebraic systems arising from elliptic optimal control problems.  Berichte aus dem Institut für Angewandte Mathematik 2023/2, Technische Universität Graz, Institut für Angewandte Mathematik, February 2023\.  submitted to LSSC 2023 Proceedings.
  * [23] U. Langer, R. Löscher, O. Steinbach, and H. Yang.  Robust finite element discretization and solvers for distributed elliptic optimal control problems.  Comput. Meth. Appl. Math., 2023.
  * [24] U. Langer, O. Steinbach, F. Tröltzsch, and H. Yang.  Space-time finite element discretization of parabolic optimal control problems with energy regularization.  SIAM J. Numer. Anal., 59:675–695, 2021.
  * [25] U. Langer, O. Steinbach, F. Tröltzsch, and H. Yang.  Unstructured space-time finite element methods for optimal control of parabolic equations.  SIAM J. Sci. Comput., 43:A744–A771, 2021.
  * [26] U. Langer, O. Steinbach, and H. Yang.  Robust discretization and solvers for elliptic optimal control problems with energy regularization.  Comput. Meth. Appl. Math., 22:97–111, 2022.
  * [27] U. Langer, O. Steinbach, and H. Yang.  Robust space-time finite element error estimates for parabolic distributed optimal control problems with energy regularization.  Technical Report arXiv:2206.06455, arXiv.org, 2022.
  * [28] W. Layton and L. Rebholz.  Approximate deconvolution models of turbulence: Analysis, phenomenology and numerical analysis, volume 2042 of Lecture Notes in Mathematics.  Springer, Heidelberg, 2012.
  * [29] J. L. Lions.  Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles.  Dunod Gauthier-Villars, Paris, 1968.
  * [30] R. Löscher and O. Steinbach.  Space-time finite element methods for distributed optimal control of the wave equation.  Technical Report arXiv:2211.02562, arXiv.org, 2022.
  * [31] K.-A. Mardal and R. Winther.  Preconditioning discretizations of systems of partial differential equations.  Numer. Linear Algebra Appl., 18(1):1–40, 2011.
  * [32] M. Neumüller and O. Steinbach.  Regularization error estimates for distributed control problems in energy spaces.  Math. Methods Appl. Sci., 44(5):4176–4191, 2021.
  * [33] Y. Notay.  Convergence of some iterative methods for symmetric saddle point linear systems.  SIAM J. Matrix Anal. Appl., 40(1):122–146, 2019.
  * [34] J. Pearson, M. Stoll, and A. Wathen.  Preconditioners for state-constrained optimal control problems with moreau-yosida penalty function.  Numer. Linear Algebra Appl., 21(1):81–97, 2014.
  * [35] J. Pearson and A. Wathen.  A new approximation of the Schur complement in preconditioners for PDE-constrained optimization.  Numer. Linear Algebra Appl., 12(5):816–829, 2012.
  * [36] A. Schiela and S. Ulbrich.  Operator preconditioning for a class of inequality constrained optimal control problems.  SIAM J. Optim., 24(1):435–466, 2014.
  * [37] J. Schöberl and W. Zulehner.  Symmetric indefinite preconditioners for saddle point problems with applications to PDE-constrained optimization problems.  SIAM J. Matrix Anal. Appl., 29:752–773, 2007.
  * [38] V. Schulz and G. Wittum.  Transforming smoothers for pde constrained optimization problems.  Comput. Visual. Sci., 11:207–219, 2008.
  * [39] O. Steinbach.  Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements.  Springer, New York, 2008.
  * [40] M. Stoll and A. Wathen.  Preconditioning for partial differential equation constrained optimization with control constraints.  Numer. Linear Algebra Appl., 19:53–71, 2012.
  * [41] F. Tröltzsch.  Optimal control of partial differential equations: Theory, methods and applications, volume 112 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, Rhode Island, 2010.
  * [42] W. Zulehner.  Analysis of iterative methods for saddle point problems: a unified approach.  Math. Comp., 71(238):479–505, 2002.
  * [43] W. Zulehner.  Nonstandard norms and robust estimates for saddle point problems.  SIAM J. Matrix Anal. Appl., 32(2):536–560, 2011.