# Linearization Errors in Discrete Goal-Oriented Error Estimation

Brian N. Granzow  D. Thomas Seidl  Stephen D. Bond  Sandia National
Laboratories
P.O. Box 5800
Albuquerque, NM 87185-1321

###### Abstract

This paper is concerned with goal-oriented _a posteriori_ error estimation for
nonlinear functionals in the context of nonlinear variational problems solved
with continuous Galerkin finite element discretizations. A two-level, or
discrete, adjoint-based approach for error estimation is considered. The
traditional method to derive an error estimate in this context requires
linearizing both the nonlinear variational form and the nonlinear functional
of interest which introduces linearization errors into the error estimate. In
this paper, we investigate these linearization errors. In particular, we
develop a novel discrete goal-oriented error estimate that accounts for
traditionally neglected nonlinear terms at the expense of greater
computational cost. We demonstrate how this error estimate can be used to
drive mesh adaptivity. We show that accounting for linearization errors in the
error estimate can improve its effectivity for several nonlinear model
problems and quantities of interest. We also demonstrate that an adaptive
strategy based on the newly proposed estimate can lead to more accurate
approximations of the nonlinear functional with fewer degrees of freedom when
compared to uniform refinement and traditional adjoint-based approaches.

###### keywords:

adjoint, goal-oriented, a posteriori error estimation, linearization error,
adaptive mesh, finite element

††journal: arXiv

##  1 Introduction

Finite element simulations have become ubiquitous in engineering practice. It
is natural to question the accuracy of a given finite element solution and, in
particular, the accuracy of post-processed output quantities upon which
critical design decisions may be made. For instance, an aerospace engineer may
be interested in accurately assessing the drag over an airfoil design. _A
posteriori_ error estimation provides a mechanism to approximate
discretization errors in this context, and, in particular, goal-oriented _a
posteriori_ error estimation can approximate the discretization error in a
chosen functional quantity of interest (QoI). Localized information obtained
from a goal-oriented error estimate can then be used to control the QoI
discretization error by adapting the finite element mesh via local mesh
modifications [23, 3].

Goal-oriented _a posteriori_ error estimation has been successfully applied in
a wide variety of contexts over the past thirty years [26, 27, 32, 34, 9, 16,
2, 20, 19]. Presently, we consider _two-level_ or _discrete_ approaches for
goal-oriented error estimation. Modern finite element software packages [7, 6,
21, 5] significantly ease the burden of implementing a two-level adjoint-based
error estimation scheme which may, in part, explain why such schemes have been
applied in a variety of contexts recently [4, 1, 10, 13, 15, 18, 31].
Presently, we have used the software PUMI [21] to implement parallel, adaptive
finite element simulations driven by goal-oriented error estimation.

Inherent in the traditional derivation of a discrete goal-oriented error
estimate are linearizations of the PDE residual and of the QoI. Once
linearized, the higher-order terms are subsequently disregarded in typical
discrete goal-oriented error estimate derivations [16]. We refer to these
neglected higher-order terms as _linearization errors_. The purpose of this
paper is to investigate the effect of including these linearization errors in
a goal-oriented error estimate and the effect that localizing these
linearization errors has on mesh adaptation.

We briefly review the relevant portions of the literature that concern
linearization errors in the context of goal-oriented error estimation.
Previous authors have included the adjoint residual in the derivation of goal-
oriented error estimates to obtain an error estimate that is _locally_ third-
order, but _globally_ second-order [32, 9, 22]. Additionally, Venditti and
Darmofal [32] quantified the error due to nonlinear effects when using a two-
level approach, but did not include these effects in their error estimate.
Şimşek et al. [30] explicitly accounted for the residual linearization error
in the derivation of a two-level goal-oriented error estimate for various
time-dependent PDEs, but did not consider nonlinear QoIs and did not consider
the _localization_ of higher-order effects to drive mesh adaptivity. Finally,
more recently Endtmayer et al. [14] and Dolejsi and Congreve [12] developed
methods to incorporate linearization errors (as well as errors that arise from
iteratively solving nonlinear algebraic equations) in a variational setting
that depend on the strong-form of the adjoint residual.

The contributions of this paper can be stated as follows. First, by
considering linearization errors, we derive a novel two-level adjoint-based
error estimate that exactly represents the QoI discretization error between
two function spaces and is independent of any form of the adjoint residual.
Second, we prove that the solution of a requisite nonlinear scalar problem for
this newly proposed estimate can be known _a priori_ when quadratic QoIs are
considered. Third, we outline a convenient adjoint verification procedure that
makes use of linearization errors and we highlight how this procedure can be
applied to existing adjoint-based error estimation codes. We then present an
adaptive mesh scheme that utilizes the newly proposed error estimate. Last, we
demonstrate that the novel error estimate can lead to more accurate error
estimates and more optimal meshes in certain scenarios when compared to a
traditional adjoint-weighted residual error estimate.

The remainder of this paper is structured as follows. First, we provide a
review of discrete goal-oriented error estimation in a general setting and
motivate our interest in studying linearization errors in this context. Next,
we introduce and derive a novel discrete adjoint-based error estimate that is
free of linearization errors. The evaluation of this estimate requires the
solution of an auxiliary nonlinear scalar problem. We then provide the
solution to this nonlinear scalar problem when quadratic QoIs are considered.
Next, we highlight that linearization errors can be used to perform adjoint
verification in a traditional adjoint-weighted residual error estimation
context. We then discuss the localization of the newly proposed error estimate
and how this localization can be used to adaptively modify meshes to reduce
QoI discretization errors. Next, we investigate a nonlinear model problem with
multiple exemplar QoIs to verify and demonstrate the effectiveness of the
newly proposed estimate. We conclude by summarizing our findings and proposing
avenues for future work.

##  2 Review of Discrete Goal-Oriented Error Estimation

###  2.1 Two-Level Strategy

Let $\Omega\subset\mathbb{R}^{n_{\text{dim}}}$ denote an open bounded domain
in $n_{\text{dim}}$ spatial dimensions with a Lipschitz continuous boundary
$\Gamma$. Consider the stationary nonlinear variational problem: find
$u\in\mathcal{V}$ such that

$\mathcal{N}(u;w)=\mathcal{L}(w),\quad\forall\,w\in\mathcal{V}.$ (1)

Here, $\mathcal{V}$ denotes an appropriate Sobolev space for the underlying
PDE, $\mathcal{N}:\mathcal{V}\times\mathcal{V}\to\mathbb{R}$ denotes a
semilinear form, nonlinear in its first argument and linear in its second,
$\mathcal{L}:\mathcal{V}\to\mathbb{R}$ denotes a linear functional, $u$
denotes the exact weak solution to the chosen PDE, and $w$ denotes an
arbitrary weighting function. Let
$\cup_{e=1}^{n_{el}}\overline{\Omega^{e}}=\overline{\Omega}$ denote a
partitioning of the domain $\Omega$ into $n_{el}$ non-overlapping elements,
such that $\Omega^{i}\cap\Omega^{j}=\varnothing$. Let
$\mathcal{V}^{H}\subset\mathcal{V}$ denote a finite-dimensional space. A
finite element formulation of the variational problem (1) can be stated: find
$u^{H}\in\mathcal{V}^{H}$ such that

$\mathcal{N}(u^{H};w^{H})=\mathcal{L}(w^{H}),\quad\forall\,w^{H}\in\mathcal{V}^{H},$
(2)

which results in a system of $N$ nonlinear algebraic equations
$\boldsymbol{R}^{H}:\mathbb{R}^{N}\to\mathbb{R}^{N}$

$\boldsymbol{R}^{H}(\boldsymbol{u}^{H})=\boldsymbol{0},$ (3)

where $\boldsymbol{u}^{H}\in\mathbb{R}^{N}$ denotes the vector of nodal
solution coefficients. The problem posed by equation (3) is referred to as the
_primal problem_. The primal problem (3) is solved via Newton’s method, where
iterations of the steps:

$\displaystyle\left[\frac{\partial\boldsymbol{R}^{H}}{\partial\boldsymbol{u}^{H}}\biggr{|}_{\boldsymbol{u}^{H}_{i}}\right]\,\delta\boldsymbol{u}^{H}_{i}=-\boldsymbol{R}^{H}(\boldsymbol{u}^{H}_{i}),$
(4)
$\displaystyle\boldsymbol{u}^{H}_{i+1}=\boldsymbol{u}^{H}_{i}+\delta\boldsymbol{u}^{H}_{i},$

are performed until the convergence criterion
$\|\boldsymbol{R}^{H}(\boldsymbol{u}^{H})\|_{2}<\text{TOL}$ is satisfied.
Here, TOL is some user-specified convergence tolerance,
$\boldsymbol{u}^{H}_{i}$ denotes the solution vector at the $i$th Newton
iteration, and $\delta\boldsymbol{u}^{H}_{i}$ denotes the incremental solution
update at the $i$th iteration.

Let $\mathcal{J}(u):\mathcal{V}\to\mathbb{R}$ denote a functional that is of
particular interest to a given simulation, referred to as the _quantity of
interest_. The aim of goal-oriented error estimation is to approximate the
discretization error in the QoI,
$\mathcal{E}:=\mathcal{J}(u)-\mathcal{J}(u^{H})$. To this end, we adopt a two-
level error estimation strategy [32, 33, 34, 16], which is sometimes referred
to as a _discrete_ approach for goal-oriented error estimation. For this
strategy, a second finite-dimensional space $\mathcal{V}^{h}$ is chosen, such
that $\mathcal{V}^{H}\subset\mathcal{V}^{h}\subset\mathcal{V}$. In this
context, $\mathcal{V}^{H}$ and $\mathcal{V}^{h}$ are referred to the _coarse
space_ and _fine space_ , respectively.

The finite element problem (2) posed on the fine space $\mathcal{V}^{h}$ leads
to a system of $n$ nonlinear algebraic equations
$\boldsymbol{R}^{h}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, where $n>N$, and

$\boldsymbol{R}^{h}(\boldsymbol{u}^{h})=\boldsymbol{0}.$ (5)

Here $\boldsymbol{u}^{h}$ denotes the solution vector of nodal coefficients
for the fine problem. Additionally, the QoI can be evaluated using the coarse
and fine space discretizations, which we denote as
$\mathcal{J}^{H}:\mathbb{R}^{N}\to\mathbb{R}$ and
$\mathcal{J}^{h}:\mathbb{R}^{n}\to\mathbb{R}$, respectively. Let
$\boldsymbol{I}^{H}_{h}:\mathcal{V}^{H}\to\mathcal{V}^{h}$ denote the
prolongation operator from the coarse space to the fine space defined by
interpolation and let $(\cdot)^{H}_{h}:=\boldsymbol{I}^{H}_{h}(\cdot)^{H}$
denote the prolongation of a coarse space vector $(\cdot)^{H}$ onto the fine
space $\mathcal{V}^{h}$. Because the exact error $\mathcal{E}$ is generally
unknowable, two-level approaches for goal-oriented error estimation introduce
computable approximations for the QoI discretization error $\mathcal{E}^{h}$
between the coarse and fine space, where

$\displaystyle\mathcal{E}^{h}$
$\displaystyle:=\mathcal{J}^{h}(\boldsymbol{u}^{h})-\mathcal{J}^{H}(\boldsymbol{u}^{H})$
(6)
$\displaystyle=\mathcal{J}^{h}(\boldsymbol{u}^{h})-\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h}).$

Here the error $\mathcal{E}^{h}$ serves as a proxy for the true discretization
error $\mathcal{E}$.

###  2.2 Choices for the Fine Space

There exist several common choices for the fine space $\mathcal{V}^{h}$. Of
primary importance in making this choice is that $\mathcal{E}^{h}$ is a _good_
proxy for the true QoI discretization error $\mathcal{E}$. With this in mind,
the fine space could be obtained by uniformly refining the mesh used for the
primal problem. We refer to this approach as $h$-enrichment. Alternatively,
the fine space could be defined by using a finite element basis with an
increased polynomial order when compared to the coarse space basis. We refer
to this approach as $p$-enrichment. Lastly, a combination of both $h$\- and
$p$-enrichment could be used to define the fine space, which we refer to as
$hp$-enrichment. As a note, more esoteric approaches can be chosen for the
fine space that we will not consider in the present context [19, 28].

The choice of the fine space affects the asymptotic behavior of the error
estimate $\mathcal{E}^{h}$. For example, consider a QoI that converges at a
rate of $k$ on the coarse space and at a rate of $k+l$ on the fine space, such
that
$\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h})=c|\Omega^{e}_{H}|^{k}$
and
$\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{h})=c|\Omega^{e}_{h}|^{k+l}$.
Here, $|\Omega^{e}_{H}|:=\text{meas}(\Omega^{e}_{H})$ and
$|\Omega^{e}_{h}|:=\text{meas}(\Omega^{e}_{h})$ denote the characteristic mesh
size on the coarse space and fine space, respectively. The ratio of the error
in the QoI between the two spaces $\mathcal{E}^{h}$ and the true
discretization error $\mathcal{E}$ can then be written as

$\displaystyle\frac{\mathcal{E}^{h}}{\mathcal{E}}$
$\displaystyle=\frac{\left(\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h})\right)-\left(\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{h})\right)}{\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h})}$
(7)
$\displaystyle=\frac{c\,|\Omega^{e}_{H}|^{k}-c\,|\Omega^{e}_{h}|^{k+l}}{c\,|\Omega^{e}_{H}|^{k}}$
$\displaystyle=1-\left[\frac{|\Omega^{e}_{h}|}{|\Omega^{e}_{H}|}\right]^{k}|\Omega^{e}_{h}|^{l}.$

When $p$-enrichment is utilized, $l>0$ and the characteristic mesh size on
both spaces is equal, with $|\Omega^{e}_{H}|=|\Omega^{e}_{h}|$. This results
in the asymptotic behavior $\mathcal{E}^{h}/\mathcal{E}\to 1$ as
$|\Omega^{e}_{H}|\to 0$. When uniform $h$-enrichment is utilized, the ratio of
characteristic mesh sizes is given as $|\Omega^{e}_{h}|/|\Omega^{e}_{H}|=1/2$
and there is no increase in the convergence rate of the QoI, with $l=0$. This
results in the asymptotic behavior $\mathcal{E}^{h}/\mathcal{E}\to
1-(1/2)^{k}$ as $|\Omega^{e}_{H}|\to 0$. Presently, we choose the fine space
to be $p$-enrichment, though this is not a requirement for the following
discussions.

###  2.3 Traditional Adjoint-Weighted Residual Estimate

Given a choice for the fine space and the notation introduced in Section 2.1,
the residual equations on the fine space can be expanded about the prolonged
coarse solution $\boldsymbol{u}^{H}_{h}$ as

$\cancelto{0}{\boldsymbol{R}^{h}(\boldsymbol{u}^{h})}=\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})+\left[\frac{\partial\boldsymbol{R}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]\boldsymbol{e}^{h}+\boldsymbol{E}^{\mathcal{R}}_{L}.$
(8)

Here the left-hand side vanishes by virtue of the governing equations on the
fine space (5),
$\boldsymbol{e}^{h}:=\left(\boldsymbol{u}^{h}-\boldsymbol{u}^{H}_{h}\right)$
denotes the solution discretization error between the coarse and fine space,
and $\boldsymbol{E}^{\mathcal{R}}_{L}$ denotes higher-order terms in the
residual Taylor expansion that we will refer to as the _residual linearization
error_. Similarly, the QoI can be expanded about the prolonged coarse solution
in a Taylor series as

$\mathcal{J}^{h}(\boldsymbol{u}^{h})=\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h})+\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]\boldsymbol{e}^{h}+\mathcal{E}^{\mathcal{J}}_{L},$
(9)

where $\mathcal{E}^{\mathcal{J}}_{L}$ denotes higher-order terms that we will
refer to as the _QoI linearization error_. Neglecting the residual
linearization error in the expansion (8), the discretization error can be
approximated linearly as

$\boldsymbol{e}^{h}\approx\boldsymbol{e}^{h}_{L}:=-\left[\frac{\partial\boldsymbol{R}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]^{-1}\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h}).$
(10)

Let $\boldsymbol{z}^{h}$ denote the solution to the so-called _adjoint
problem_ , given as

$\left[\frac{\partial\boldsymbol{R}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]^{T}\boldsymbol{z}^{h}=\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]^{T}.$
(11)

Substituting the linearized approximation for the discretization error (10)
into the QoI Taylor expansion (9), neglecting the QoI linearization error, and
utilizing the adjoint solution (11) yields the well-known adjoint-weighted
residual error estimate:

$\mathcal{E}^{h}\approx\eta_{1}:=-\boldsymbol{z}^{h}\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{h}_{H}).$
(12)

The introduction of the adjoint problem is motivated by two considerations.
First, it replaces a potentially cost prohibitive nonlinear primal solve on
the fine space with a cheaper linear solve on the fine space. Depending on the
nonlinearity of the primal PDE and the choice of fine space, the cost of
solving the adjoint problem (11) on the fine space can be comparable or
cheaper than the cost of solving the primal problem on the coarse space (3).
Second, and perhaps most important, the adjoint solution itself provides
information about the spatial distribution of contributions to the QoI
discretization error. Concretely, this is because $\boldsymbol{z}^{h}$
represents the linear sensitivity of the QoI with respect to perturbations in
the PDE residual.

For example, consider a point-wise QoI inside of the computational domain. We
could naïvely solve for the linearized discretization error (10) for the
purpose of replacing a nonlinear solve on the fine space with a linear solve
on the fine space. We could then use the linearized discretization error
$\boldsymbol{e}^{h}_{L}$ to approximate the fine space solution as
$\boldsymbol{u}^{h}\approx\boldsymbol{u}^{H}_{h}+\boldsymbol{e}^{h}_{L}$ and,
in turn, approximate the QoI error $\mathcal{E}^{h}$. In fact, this process
would yield an identical error estimate to $\eta_{1}$. However, in doing so,
we would obtain no information about the spatial distribution of the QoI error
in the domain except at the point of interest and would thus have limited
tools to control the discretization error in the QoI via mesh adaptivity.

###  2.4 Motivation For Studying Linearization Errors

Much has been said about the error estimate $\eta_{1}$ [16]. It is, however,
less common to consider the effect of higher-order information when deriving
goal-oriented error estimates [30, 14, 12], particularly in a discrete adjoint
setting. This is the primary aim of the present work. Consider, for example,
the function $f(x)=x^{2}$, whose Taylor expansion about the point $0$ can be
expressed exactly as
$f(x)=f(0)+f^{\prime}(0)x+\frac{1}{2}f^{\prime\prime}(0)x^{2}$. At any
evaluation point $a$, the entire Taylor series expansion is encoded in the
second derivative term, with the linearized approximation $f(a)-f(0)\approx
f^{\prime}(0)a$ providing no meaningful information about the difference
$f(a)-f(0)$ as it simply evaluates to $0$. Unfortunately, this is exactly the
approximation that is made in the derivation of the traditional adjoint-
weighted residual estimate. If instead we include higher-order information in
the function’s Taylor expansion, where
$f(a)-f(0)=f^{\prime}(0)a+f^{\prime\prime}(0)a^{2}$, we recognize that we will
recover the exact difference $f(0)-f(a)$.

Further, let $u^{H}$ denote the function corresponding to the nodal
coefficient vector $\boldsymbol{u}^{H}$ and let $e:=u-u^{H}$ denote the exact
discretization error. Consider the quadratic functional
$\mathcal{J}(u)=\int_{\Omega}\nabla u\cdot\nabla u\;\text{d}\Omega$. The
discretization error in this functional can be exactly represented as

$\mathcal{J}(u)-\mathcal{J}(u^{H})=2\int_{\Omega}\nabla u^{H}\cdot\nabla
e\;\text{d}\Omega+\int_{\Omega}\nabla e\cdot\nabla e\;\text{d}\Omega.$ (13)

Here, we recognize that the first integral on the right-hand side of the error
expression (13) corresponds to the linearization used for the QoI in the
traditional adjoint-weighted residual approach and the second integral
corresponds to the linearization error introduced in the estimate. Naturally,
for a convergent discretization $e\to 0$ as $H\to 0$, and one would expect the
significance of the linearization error to quickly diminish as $H\to 0$ since
it is $\mathcal{O}(e^{2})$. However, the integrand in the linearization error
term is strictly positive, while the first integral has the potential to tend
to zero quickly relative to the linearization error due to subtractive
cancellation. In this scenario, if the linearization error is neglected then
there could be a significant under-prediction of the actual error.

These are our primary motivations for studying linearization errors and
including higher-order information in the derivation of a goal-oriented error
estimate. While the above scenarios may seem pathological, we provide concrete
examples where the traditional error estimate $\eta_{1}$ performs poorly for
nonlinear QoIs in Section 6. While we are presently concerned with extending
discrete adjoint-weighted error estimates to include nonlinear effects, it is
also worth noting that the reliability of the traditional dual weighted
residual error [9] estimate, which is derived in a continuous variational
setting, has also been demonstrated to significantly under-predict errors [25]
in specific contexts.

##  3 An Estimate Free of Linearization Errors

###  3.1 Derivation

We can entirely account for the linearization errors
$\boldsymbol{E}^{\mathcal{R}}_{L}$ and $\mathcal{E}^{\mathcal{J}}_{L}$ in the
derivation of a two-level adjoint-based error estimate. From the mean-value
theorem, there exists a vector $\boldsymbol{u}^{*}$ where the QoI
linearization error $\mathcal{E}^{\mathcal{J}}_{L}$ vanishes, such that

$\mathcal{J}^{h}(\boldsymbol{u}^{h})=\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h})+\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}}\right]\boldsymbol{e}^{h}.$
(14)

Here $\boldsymbol{u}^{*}$ denotes a point on a linear path between the coarse-
space solution $\boldsymbol{u}^{H}_{h}$ and the fine-space solution
$\boldsymbol{u}^{h}$, written as

$\boldsymbol{u}^{*}(\theta)=\boldsymbol{u}^{H}_{h}+\theta\boldsymbol{e}^{h},\quad\theta\in[0,1].$
(15)

Finding the appropriate value of $\theta$, and thus $\boldsymbol{u}^{*}$,
amounts to solving the nonlinear scalar equation

$f(\theta):=\mathcal{E}^{h}-\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}(\theta)}\right]\boldsymbol{e}^{h}=0.$
(16)

We solve this problem with Newton’s method, where the iterations

$\theta_{i+1}=\theta_{i}-\frac{f(\theta_{i})}{f^{\prime}(\theta_{i})},$ (17)

are performed until $|f(\theta)|<\text{TOL}$ is satisfied. Here, TOL denotes
some user-specified convergence tolerance and $f^{\prime}$ denotes the
derivative of the function $f$ with respect to $\theta$, which can be
expressed using the chain rule as

$f^{\prime}(\theta)=-\left[\boldsymbol{e}^{h}\right]^{T}\left[\frac{\partial^{2}\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}(\theta)}\right]\boldsymbol{e}^{h}.$
(18)

Here, the cost of solving the nonlinear scalar problem (16) involves only
assembly and does not require matrix inversion.

Unfortunately, no direct analogue to the mean-value theorem exists for vector-
valued functions [24]. Thus, the residual linearization error
$\boldsymbol{E}^{\mathcal{R}}_{L}$ must be accounted for in a different
manner. Given the vector $\boldsymbol{u}^{*}$, we introduce the _modified
adjoint problem_ :

$\left[\frac{\partial\boldsymbol{R}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]^{T}\boldsymbol{z}^{*}=\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}}\right]^{T}.$
(19)

Solving for the discretization error $\boldsymbol{e}^{h}$ in the residual
Taylor expansion (8), substituting the resulting expression for
$\boldsymbol{e}^{h}$ into the modified QoI Taylor expansion (14), and using
the solution $\boldsymbol{z}^{*}$ to the modified adjoint (19) problem yields
the exact error representation

$\mathcal{E}^{h}=-\boldsymbol{z}^{*}\cdot\left(\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})+\boldsymbol{E}^{\mathcal{R}}_{L}\right).$
(20)

For the purposes of mesh adaptivity, it is desirable to express the error only
in terms of the residual scaled by a weighting vector. To this end, we
endeavor to find a modified vector $\boldsymbol{z}^{**}$ that satisfies

$\boldsymbol{z}^{**}\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})=\boldsymbol{z}^{*}\cdot\left(\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})+\boldsymbol{E}^{\mathcal{R}}_{L}\right).$
(21)

The exact least-squares solution to the above equality is

$\boldsymbol{z}^{**}=\boldsymbol{z}^{*}+\frac{\boldsymbol{z}^{*}\cdot\boldsymbol{E}^{\mathcal{R}}_{L}}{\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})}\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h}),$
(22)

and the QoI discretization between the coarse and fine space can be exactly
expressed as

$\mathcal{E}^{h}=\eta_{2}:=-\boldsymbol{z}^{**}\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h}).$
(23)

Unlike the traditional adjoint-weighted residual error estimate $\eta_{1}$,
the estimate $\eta_{2}$ accounts for the contributions to the error that arise
from the residual and QoI linearization errors,
$\boldsymbol{E}^{\mathcal{R}}_{L}$ and $\mathcal{E}^{\mathcal{J}}_{L}$,
respectively. Previous techniques that include higher-order effects in goal-
oriented error estimates have been derived in a continuous adjoint setting and
include the strong form of the adjoint residual in their evaluation [9, 30,
14, 12], a term which can be difficult to compute or even derive [28]. In the
present discrete adjoint context, the estimate $\eta_{2}$ includes higher-
order terms without requiring an evaluation of the adjoint residual.

The evaluation of $\eta_{2}$ is not without practical downsides, however, as
it requires the solution of the nonlinear primal problem on the fine space (5)
to evaluate the term $\boldsymbol{e}^{h}$, which is present in the nonlinear
scalar equation (16) and in the evaluation of the residual linearization error
(8). This is a computationally expensive proposition. Similarly to Nochetto et
al. [25], we propose that this method can be used to safeguard traditional
adjoint-weighted residual methods at coarse mesh resolutions, where the
solution of the fine space nonlinear problem remains tractable. At such mesh
resolutions, the traditional adjoint-weighted residual estimate $\eta_{1}$ may
severely under-predict the QoI discretization error and cause an iterative
adaptive algorithm to terminate prematurely. We demonstrate the under-
prediction behavior of $\eta_{1}$ in Section 6 and have provided a potential
explanation for this behavior for nonlinear QoIs in Section 2.4.

In fact, it is natural to question why it is even necessary to evaluate the
term $\boldsymbol{z}^{**}$ when the QoI discretization error $\mathcal{E}^{h}$
can be found exactly when the coarse space solution $\boldsymbol{u}^{H}$ and
fine space solution $\boldsymbol{u}^{h}$ are known. The answer, analogously to
the traditional adjoint-weighted residual approach, is that the term
$\boldsymbol{z}^{**}$ provides _local_ information about the spatial
distribution of the QoI discretization error which can be used for the
purposes of mesh adaptivity. In particular, the term $\boldsymbol{z}^{**}$
also provides local information about the residual and QoI linearization
errors $\boldsymbol{E}^{\mathcal{R}}_{L}$ and $\mathcal{E}^{\mathcal{J}}_{L}$,
respectively, whereas the traditional adjoint solution $\boldsymbol{z}^{h}$
does not. As an avenue for future work, one could attempt to approximate the
term $\boldsymbol{e}^{h}$ with patch recovery techniques [36, 35, 29, 17] to
reduce the computational cost of the newly proposed estimate.

Lastly, we remark that while the estimate $\eta_{2}$ exactly represents the
QoI discretization error $\mathcal{E}^{h}$ between the coarse and the fine
spaces, there is a remaining error term
$\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{h})$, such that the exact QoI
discretization is given as
$\mathcal{E}=\eta_{2}+\mathcal{J}(u)-\mathcal{J}^{h}(\boldsymbol{u}^{h})$.
Naturally, this remainder term is, in general, unknowable. However, as another
avenue for future work, one could envision using traditional functional
analysis techniques found in the goal-oriented error estimation literature to
bound this remaining term and combine its effects with the estimate
$\eta_{2}$.

###  3.2 Quadratic Quantities of Interest

######  Proposition 1.

If $\mathcal{J}(u)$ is a quadratic functional, then the solution to the
problem (16) is $\theta=\nicefrac{{1}}{{2}}$.

###### Proof.

Let $\mathcal{J}(u)$ denote a quadratic functional. We begin by equating the
expression for $\mathcal{E}^{h}$ obtained from the Taylor expansion (14) with
the exact Taylor expansion for a quadratic functional expanded about the
prolonged coarse solution, yielding

$\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}}\right]\boldsymbol{e}^{h}=\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]\boldsymbol{e}^{h}+\frac{1}{2}\left[\boldsymbol{e}^{h}\right]^{T}\left[\frac{\partial^{2}\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\right]\boldsymbol{e}^{h}.$
(24)

Because the QoI is quadratic, its first derivative is a linear operator acting
on the solution state and its second derivative is constant, independent of
the solution state. The above equality must hold for arbitrary discretization
errors $\boldsymbol{e}^{h}$, with

$\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}}\right]=\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{H}_{h}}\right]+\frac{1}{2}\left[\boldsymbol{e}^{h}\right]^{T}\left[\frac{\partial^{2}\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\right].$
(25)

Due to the linearity of the QoI’s first derivative, we can write

$\left[\frac{\partial\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}}\biggr{|}_{\boldsymbol{u}^{*}-\boldsymbol{u}^{H}_{h}}\right]=\frac{1}{2}\left[\boldsymbol{e}^{h}\right]^{T}\left[\frac{\partial^{2}\mathcal{J}^{h}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\right],$
(26)

and because the QoI is quadratic we can write

$\left(\boldsymbol{u}^{*}-\boldsymbol{u}^{H}_{h}\right)^{T}\left[\frac{\partial^{2}\mathcal{J}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\right]=\frac{1}{2}\left[\boldsymbol{e}^{h}\right]^{T}\left[\frac{\partial^{2}\mathcal{J}}{\partial\boldsymbol{u}^{h}\partial\boldsymbol{u}^{h}}\right].$
(27)

Using the definition of $\boldsymbol{e}^{h}$ and because the above expansions
hold for arbitrary quadratic QoIs, we obtain

$\boldsymbol{u}^{*}=\frac{1}{2}\left(\boldsymbol{u}^{H}_{h}+\boldsymbol{u}^{h}\right).$
(28)

From the definition of $\boldsymbol{u}^{*}$, the result
$\theta=\nicefrac{{1}}{{2}}$ follows immediately. ∎

This result allows us to entirely circumvent a nonlinear solve when a
quadratic QoI is chosen, and motivates an initial guess of
$\theta_{0}=\nicefrac{{1}}{{2}}$ for Newton’s method when solving (16).

##  4 Adjoint Verification

One surprising benefit of studying linearization errors manifests itself in
the realm of verification. For complex nonlinear physics, it is natural to
question whether the implementation of the adjoint problem (11) leads to the
correct adjoint solution $\boldsymbol{z}^{h}$. Answering this question with
traditional verification techniques may be a burdensome and non-trivial task.
It turns out, however, that exactly solving for the residual linearization
error $\boldsymbol{E}^{\mathcal{R}}_{L}$ provides a convenient and
straightforward mechanism for adjoint verification.

This verification procedure can be described as follows. Assume that we have
the machinery to compute every term on the right-hand side of equation (8), so
that we can solve for the residual linearization error
$\boldsymbol{E}^{R}_{L}$ exactly. Additionally, assume that we have the
machinery to solve the adjoint problem (11) on the fine space for
$\boldsymbol{z}^{h}$. Choose $\mathcal{J}(u)$ to be a _linear_ functional,
such that $\mathcal{E}^{\mathcal{J}}_{L}=0$. We can exactly express the error
$\mathcal{E}^{h}$ in this context as a function of the adjoint solution
$\boldsymbol{z}^{h}$:

$\displaystyle\mathcal{E}^{h}$
$\displaystyle:=\mathcal{J}^{h}(\boldsymbol{u}^{h})-\mathcal{J}^{h}(\boldsymbol{u}^{H}_{h}),$
(29)
$\displaystyle=-\boldsymbol{z}^{h}\cdot\left(\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})+\boldsymbol{E}^{\mathcal{R}}_{L}\right),$
$\displaystyle=\eta_{1}+\eta^{\mathcal{R}}_{L},$

where
$\eta^{\mathcal{R}}_{L}:=-\boldsymbol{z}^{h}\cdot\boldsymbol{E}^{\mathcal{R}}_{L}$.
As a verification check, we ensure that the effectivity $\mathcal{I}^{v}$,
defined as

$\mathcal{I}^{v}:=\frac{\eta_{1}+\eta^{\mathcal{R}}_{L}}{\mathcal{E}^{h}}$
(30)

evaluates to 1. Here we make several remarks. First, the evaluation of
$\boldsymbol{E}^{\mathcal{R}}_{L}$ from the Taylor expansion (8) involves only
data obtained from the primal problem. Thus, if one has confidence in the
correctness of the primal problem implementation, then the linearization error
is accurately obtained through linear algebra. Second, the denominator of the
effectivity (30) contains only data from the primal problem and is independent
of the adjoint solution. This implies that $\mathcal{I}^{v}=1$ is a true
verification check in the sense that it relates the adjoint solution to a
computable value that does not involve the adjoint implementation.

For codes that were not designed with two-level error estimation in mind but
that still have an adjoint capability, the above verification procedure can be
performed using $h$-enrichment for the fine space. In this context, the primal
problem can readily be solved on the coarse and fine mesh. The only source of
additional development effort stems from the prolongation of the coarse
solution to the fine mesh and the evaluation of the terms in equation (8).
Depending on software design, the latter step may not be trivial, but seems
well worth the effort for the verification benefit it provides.

##  5 Mesh Adaptivity

###  5.1 Error Localization

To control discretization errors, it is necessary to _localize_ the error
estimate into contributions to the error at the mesh entity level. These
contributions are often called _correction indicators_. Both estimates,
$\eta_{1}$ and $\eta_{2}$, take the form of the PDE residual weighted by an
adjoint vector. Because the original variational problem (2) is posed with a
Galerkin finite element method, Galerkin orthogonality allows us to write
equivalent expressions for the estimates $\eta_{1}$ and $\eta_{2}$ as

$\eta_{1}=-\left(\boldsymbol{z}^{h}-\boldsymbol{z}^{h}_{H}\right)\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h}),$
(31)

and

$\eta_{2}=-\left(\boldsymbol{z}^{**}-\boldsymbol{z}^{**}_{H}\right)\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h}),$
(32)

respectively. Here $\boldsymbol{z}^{h}_{H}$ and $\boldsymbol{z}^{**}_{H}$
denote the restriction via interpolation of the adjoint solutions
$\boldsymbol{z}^{h}$ and $\boldsymbol{z}^{**}$ onto the coarse space
$\mathcal{V}^{H}$, respectively. The subtraction of the coarse space
representation of the adjoint solution is common in the literature and is
done, in part, so that the residual is weighted by a term that indicates where
the adjoint is _poorly approximated_ rather than by a term that indicates
where the adjoint is large in value [22].

A common localization approach in finite volume and discontinuous Galerkin
methods considers decomposing the adjoint-weighted residual into a sum of
absolute values of element-level adjoint-weighted residuals, for instance as
$\eta_{1}\leq\sum_{e=1}^{n_{el}}\left|\left(-(\boldsymbol{z}^{h}-\boldsymbol{z}^{h}_{H})\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})\right)\bigr{|}_{\Omega_{e}}\right|$,
where the absolute value term at each element is used as the correction
indicator. This approach, however, may not be optimal for continuous Galerkin
finite element methods as it does not account for inter-element cancellation,
where the sum of the resultant indicators may vastly overestimate the error
and, in turn, lead to a suboptimal adaptive strategy [16].

Presently, we make use of a variational technique for localization developed
by Richter and Wick [28], which expresses the estimates (31) and (32)
analogously in variational form and achieves localization by inserting a
partition of unity $\sum_{i}\phi^{i}=1$ into the weighting function slot of
the variational residual $\mathcal{R}$, to be defined. We take this partition
of unity to be linear Lagrange basis functions which results in the error
being localized to the $n_{vtx}$ mesh vertices. Concretely, this results in
localized error contributions $\eta^{i}_{1}$ and $\eta^{i}_{2}$,
$i=1,2,\dots,n_{vtx}$ at mesh vertices of the form

$\eta_{1}=\sum_{i=1}^{n_{vtx}}\underbrace{-\mathcal{R}(u^{H};(z^{h}-z^{h}_{H})\phi^{i})}_{\eta^{i}_{1}}$
(33)

and

$\eta_{2}=\sum_{i=1}^{n_{vtx}}\underbrace{-\mathcal{R}(u^{H};(z^{**}-z^{**}_{H})\phi^{i})}_{\eta^{i}_{2}},$
(34)

where $\mathcal{R}:\mathcal{V}\times\mathcal{V}\to\mathbb{R}$ denotes the
variational residual $\mathcal{R}(u,w):=\mathcal{L}(w)-\mathcal{N}(u;w)$ and
$u^{H},z^{h},z^{h}_{H},z^{**},z^{**}_{H}$ denote the function counterparts of
the nodal coefficient vectors
$\boldsymbol{u}^{H},\boldsymbol{z}^{h},\boldsymbol{z}^{h}_{H},\boldsymbol{z}^{**}$,
and $\boldsymbol{z}^{**}_{H}$, respectively, as determined by the function
spaces $\mathcal{V}^{H}$ and $\mathcal{V}^{h}$. As a final step, to compute an
element level correction indicator $\eta^{e}$, we interpolate the value of the
vertex-based error contribution $\eta^{i}$ to element centers and take its
absolute value.

Here we remark that a primary motivation in the localization of the estimate
$\eta_{2}$ is to additionally localize the effects of the linearization errors
$\mathcal{E}^{\mathcal{J}}_{L}$ and $\boldsymbol{E}^{\mathcal{R}}_{L}$ during
the mesh adaptation process. That is, the localized indicator $\eta^{i}_{2}$
will also indicate regions in the domain where the linearization errors are
large. By accurately resolving the mesh in such regions, we endeavor to obtain
more optimal meshes than those obtained from the use of the traditional
adjoint-weighted residual estimate $\eta_{1}$.

###  5.2 Mesh Size Field

Once element-level correction indicators $\eta^{e}$, $e=1,2,\dots,n_{el}$ have
been computed, they must be used to modify the finite element mesh. Naturally,
we would like to refine regions of the mesh that the indicators $\eta^{e}$
suggest contribute strongly to the QoI discretization error and, if
applicable, coarsen regions of the mesh that do not greatly affect the QoI
discretization error. To this end, we use isotropic _conformal_ mesh
adaptation driven by a _mesh size field_ , which produces an adapted mesh
without hanging nodes. This process achieves refinement and coarsening via a
series of edge splits, swaps, and collapses [23, 3]. The mesh size field
indicates the desired characteristic element size in the adapted mesh. We
choose a mesh size field that attempts to equidistribute the error in the
resultant mesh given a target $T$ number of elements, as described by Bousetta
et al. [11]. Here, the size field is found by assuming that the error is
uniformly distributed in the resultant adapted mesh and comparing the ratio of
the element contribution to the error in the current mesh to the desired
adapted mesh as related by expected convergence rates. This leads to an
optimization problem that can be solved analytically for the new mesh size.
Let $p$ denote the polynomial interpolant order of the function space
$\mathcal{V}^{H}$ in $d$ spatial dimensions. Let
$|\Omega^{e}_{H}|:=\text{meas}(\Omega^{e})$ denote the characteristic mesh
size of element $\Omega^{e}$. The desired new element size
$|\Omega^{e}_{H}|^{\text{new}}$ can then be expressed as

$\left|\Omega^{e}_{H}\right|^{\text{new}}=\left(T^{-1}\sum_{e=1}^{n_{el}}\left(\eta^{e}\right)^{\frac{2d}{2p+d}}\right)^{\frac{1}{d}}\left(\eta^{e}\right)^{\frac{-2}{2p+d}}\left|\Omega^{e}_{H}\right|.$
(35)

As a final consideration we clamp the value of the new element size to be no
less than half of the current element size and no greater than twice the
current element size, such that
$\nicefrac{{1}}{{2}}|\Omega^{e}_{H}|\leq|\Omega^{e}_{H}|^{\text{new}}\leq
2|\Omega^{e}_{H}|$.

##  6 Results

###  6.1 A Nonlinear Poisson’s Problem

As a model problem, we consider a nonlinear Poisson’s equation, written as

$\begin{cases}\begin{aligned} -\nabla\cdot[(1+\alpha u^{2})\nabla
u]&=f,&&\text{in}\;\Omega,\\\ u&=0,&&\text{on}\;\Gamma.\\\
\end{aligned}\end{cases}$ (36)

where $\alpha\in\mathbb{R}_{\geq 0}$ denotes a scaling parameter that controls
the degree of nonlinearity. The forcing function $f$ is assumed to be
sufficiently smooth, such that $f\in L^{2}(\Omega)$. Let $\mathcal{V}$ denote
the Sobolev space $H^{1}(\Omega)$ whose trace vanishes on the boundary
$\Gamma$, such that $\mathcal{V}:=\left\\{v\in
H^{1}(\Omega):v=0\;\text{on}\;\Gamma\right\\}$. Placing the model problem in
weak form yields the semilinear form

$\mathcal{N}(u;w):=\int_{\Omega}(1+\alpha u^{2})\nabla u\cdot\nabla
w\,\text{d}\Omega,$ (37)

and the linear functional

$\mathcal{L}(w):=\int_{\Omega}fw\,\text{d}\Omega.$ (38)

Let $\mathcal{V}^{H}$ and $\mathcal{V}^{h}$ denote finite dimensional function
spaces defined as

$\mathcal{V}^{H}:=\left\\{u^{H}:u^{H}\in\mathcal{V},\,u^{H}|_{\Omega^{e}}\in\mathbb{P}^{1}(\Omega)^{n_{\text{dim}}}\right\\},$
(39)

and

$\mathcal{V}^{h}:=\left\\{u^{h}:u^{h}\in\mathcal{V},\,u^{h}|_{\Omega^{e}}\in\mathbb{P}^{2}(\Omega)^{n_{\text{dim}}}\right\\},$
(40)

respectively, where $\mathbb{P}^{p}$ denotes the space of piecewise
polynomials of order $p$ over the elements $\Omega^{e}$, $e=1,2,\dots,n_{el}$.

For all results below, we consider a square geometry with a square hole, such
that
$\Omega:=\\{\boldsymbol{x}:\boldsymbol{x}\in(-1,1)\times(-1,1)\setminus[-\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}]\times[-\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}]\\}$,
discretized with triangular elements. To eliminate or drastically reduce
additional sources of numerical errors, all results use the following
operations: requisite integral evaluations are computed using a $6^{th}$
order, 12 point Gaussian quadrature rule over triangles, nonlinear systems are
solved using an absolute tolerance of $\text{TOL}=10^{-10}$, linear systems
are solved with GMRES using an absolute tolerance of $10^{-12}$, and
derivative quantities (the residual Jacobian, the residual Jacobian adjoint,
the QoI gradient, and the QoI Hessian) are computed with automatic
differentiation.

We consider the following four quantities of interest:

$\displaystyle\mathcal{J}_{1}(u)$
$\displaystyle:=\int_{\Omega}u\,\text{d}\Omega,$ (41)
$\displaystyle\mathcal{J}_{2}(u)$
$\displaystyle:=\int_{\Omega_{s}}u^{3}\,\text{d}\Omega,$
$\displaystyle\mathcal{J}_{3}(u)$ $\displaystyle:=\int_{\Omega_{s}}\nabla
u\cdot\nabla u\,\text{d}\Omega,$ $\displaystyle\mathcal{J}_{4}(u)$
$\displaystyle:=\int_{\Omega_{s}}\sqrt{\nabla u\cdot\nabla
u}\,\text{d}\Omega.$

Here,
$\Omega_{s}:=\\{\boldsymbol{x}:\boldsymbol{x}\in(0,1)\times(-1,0)\setminus[0,\nicefrac{{1}}{{2}}]\times[-\nicefrac{{1}}{{2}},0]\\}$
denotes a subdomain of the overall domain $\Omega$. Figure 1 illustrates the
domain $\Omega$, its initial mesh, and the definition of the subdomain
$\Omega_{s}$.

Figure 1: The domain $\Omega$, its initial mesh, and the subdomain
$\Omega_{s}$ (red outline).

####  6.1.1 A Manufactured Solution

We first consider a manufactured solution on the domain $\Omega$ of the form:

$u(x,y)=\sin(2\pi x)\sin(2\pi y)\exp(\nicefrac{{5}}{{2}}(x+y)),$ (42)

which satisfies the homogeneous Dirichlet boundary conditions. The right-hand
side forcing function $f$ is found by substituting the exact expression for
$u$ into the original governing equations (36). The solution and its gradient
components are shown in Figure 2. The numerical values of the QoIs
$\mathcal{J}_{i},\;i=1,2,3,4,$ are given in A.

Figure 2: The manufactured solution (left), the $x$-component of the
manufactured solution gradient (center), and the $y$-component of the
manufactured solution gradient (right) with $20$ linearly-spaced contour
lines.

####  6.1.2 Adjoint Verification

As an initial test with the manufactured solution, we consider the mesh with
$192$ elements shown in Figure 1. We then consider the following values of
$\alpha\in\left\\{0,10^{-4},10^{-3},10^{-2},10^{-1}\right\\}$, where
$\alpha=0$ corresponds to a linear PDE and increasing values of $\alpha$
correspond to a greater degree of nonlinearity in the PDE. We solve the primal
problem on the coarse and fine spaces, $\mathcal{V}^{H}$ and
$\mathcal{V}^{h}$. For context, Newton’s method takes $0,4,7,11,17$ iterations
to solve the primal problem on both spaces for the corresponding increasing
values of $\alpha$. Using the coarse and fine space primal solutions, we
compute the quantities $\mathcal{E}^{h}$, $\eta_{1}$ and
$\eta^{\mathcal{R}}_{L}$, as required for the verification procedure outlined
in Section 4. We consider the linear quantity of interest $\mathcal{J}_{1}(u)$
and demonstrate that the adjoint solution $\boldsymbol{z}^{h}$ passes the
verification check, as shown in the right-most column of Table 1.
Additionally, the second column of this table illustrates that the norm of the
residual linearization error $\boldsymbol{E}^{\mathcal{R}}_{L}$ is increasing
with increasing values of $\alpha$. Lastly, the sixth column of these tables
illustrates that the traditional error estimate $\eta_{1}$ does not entirely
capture the QoI discretization error $\mathcal{E}^{h}$ between the fine and
coarse spaces due to the residual linearization error
$\boldsymbol{E}^{\mathcal{R}}_{L}$. However, correcting the traditional
estimate with the term $\eta^{\mathcal{R}}_{L}$ that corresponds to residual
linearization errors entirely recovers the error $\mathcal{E}^{h}$.

$\alpha$ | $||\boldsymbol{E}^{\mathcal{R}}_{L}||_{2}$ | $\mathcal{E}^{h}$ | $\eta_{1}$ | $\eta^{\mathcal{R}}_{L}$ | $\nicefrac{{\eta_{1}}}{{\mathcal{E}^{h}}}$ | $\nicefrac{{(\eta_{1}+\eta^{\mathcal{R}}_{L})}}{{\mathcal{E}^{h}}}$
---|---|---|---|---|---|---
$0$ | 1.8785e$-$13 | 3.2696e$-$01 | 3.2696e$-$01 | 1.6688e$-$15 | 1.0000e+00 | 1.0000e+00
$10^{-4}$ | 1.7321e+00 | 2.9114e$-$01 | 3.0420e$-$01 | -1.3061e$-$02 | 1.0449e+00 | 1.0000e+00
$10^{-3}$ | 1.2996e+01 | 1.3244e$-$01 | 1.9139e$-$01 | -5.8949e$-$02 | 1.4451e+00 | 1.0000e+00
$10^{-2}$ | 8.5432e+01 | 7.6441e$-$04 | 1.1158e$-$01 | -1.1082e$-$01 | 1.4597e+02 | 1.0000e+00
$10^{-1}$ | 7.7580e+02 | 9.4845e$-$02 | 2.6463e$-$01 | -1.6978e$-$01 | 2.7901e+00 | 1.0000e+00
Table 1: Adjoint verification data for the linear QoI $\mathcal{J}_{1}(u)$ for
the manufactured solution.

####  6.1.3 Modified Adjoint Verification

As a second test with the manufactured solution, we consider the mesh with
$192$ elements shown in Figure 1, choose $\alpha=10^{-2}$, and consider the
nonlinear QoIs $\mathcal{J}_{2}(u),\mathcal{J}_{3}(u),$ and
$\mathcal{J}_{4}(u)$, so the residual and QoI linearization errors,
$\boldsymbol{E}^{\mathcal{R}}_{L}$ and $\mathcal{E}^{\mathcal{J}}_{L}$,
respectively, are nonzero. As a verification check, we ensure the estimate
$\eta_{2}$ exactly recovers the QoI discretization error $\mathcal{E}^{h}$
between the coarse and fine space, such that
$\nicefrac{{\eta_{2}}}{{\mathcal{E}^{h}}}=1$. This provides assurance that the
value $\boldsymbol{z}^{**}$ is computed correctly. Table 2 shows results
obtained from this verification check. In particular, the right-most column
illustrates that the traditional error estimate $\eta_{1}$ does not exactly
recover the error $\mathcal{E}^{h}$ in the presence of residual and
linearization errors. In contrast, the right-most column of the table
demonstrates that the newly proposed estimate $\eta_{2}$ does indeed exactly
recover the error $\mathcal{E}^{h}$, thus completing the verification check.
Figure 3 illustrates the adjoint solutions $z:=\lim_{H\to
0}\boldsymbol{z}^{H}$ for the QoIs $\mathcal{J}_{2}(u),\mathcal{J}_{3}(u)$ and
$\mathcal{J}_{4}(u)$. Lastly, Figure 4 illustrates the form of the function
$f(\theta)$, as described by equation (16), for the QoI $\mathcal{J}_{2}(u)$.
Here we remark that even though the functional form of $f(\theta)$ is
quadratic in $\theta$, it remains very close to linear behavior in the domain
$\theta\in[0,1]$, and its endpoints are very nearly equal and opposite. As a
result, our initial guess of $\theta_{0}=\nicefrac{{1}}{{2}}$ for Newton’s
method seems well justified in this instance.

| $\mathcal{E}^{h}$ | $\eta_{1}$ | $\eta_{2}$ | $\nicefrac{{\eta_{1}}}{{\mathcal{E}^{h}}}$ | $\nicefrac{{\eta_{2}}}{{\mathcal{E}^{h}}}$
---|---|---|---|---|---
$\mathcal{J}_{2}(u)$ | 1.1417e+00 | 8.7186e$-$01 | 1.1417e+00 | 7.6362e$-$01 | 1.0000e+00
$\mathcal{J}_{3}(u)$ | 1.5043e+01 | -1.8944e+00 | 1.5043e+01 | -1.2593e$-$01 | 1.0000e+00
$\mathcal{J}_{4}(u)$ | 3.3263e$-$01 | -9.4750e$-$02 | 3.3263e$-$01 | -2.8485e$-$01 | 1.0000e+00
Table 2: Modified adjoint verification data for the nonlinear QoIs
$\mathcal{J}_{2}(u),\mathcal{J}_{3}(u),\mathcal{J}_{4}(u)$ for the
manufactured solution when $\alpha=10^{-2}$.

Figure 3: The adjoint solution $z$ corresponding to the model problem with a
manufactured solution when $\alpha=10^{-2}$ for the nonlinear QoIs
$\mathcal{J}_{2}(u)$ (left), $\mathcal{J}_{3}(u)$ (center), and the
$\mathcal{J}_{4}(u)$ (right) with 20 linearly spaced contour curves. Figure 4:
The function $f(\theta)$ plotted at $10$ evenly spaced points in $[0,1]$ for
the QoI $\mathcal{J}_{2}(u)$ for the manufactured solution when
$\alpha=10^{-2}$.

####  6.1.4 Asymptotic Behavior

In this section, we investigate the behavior of the nonlinear functional
quantities and the estimates $\eta_{1}$ and $\eta_{2}$ as $|\Omega^{e}_{H}|\to
0$. For all examples, we begin with an initial mesh with $192$ elements, as
shown in Figure 1 and consider a sequence of uniformly refined meshes, where
the number of elements is $n_{el}=192,768,3072,12288$ on this sequence of
meshes. For each QoI, at each mesh instance we solve the primal problems (3)
and (5), solve the traditional adjoint equation (11), solve the modified
adjoint equation (19), compute the error estimate $\eta_{1}$ from (12), and
compute the error estimate $\eta_{2}$ from (23), where the scaling parameter
$\alpha$ is chosen to be $\alpha=10^{-2}$.

Using the error estimates $\eta_{1}$ and $\eta_{2}$, we compute four
effectivity indices:
$\nicefrac{{\eta_{1}}}{{\mathcal{E}}},\nicefrac{{\eta_{2}}}{{\mathcal{E}}},\nicefrac{{\eta_{1}}}{{\mathcal{E}^{h}}}$,
and $\nicefrac{{\eta_{2}}}{{\mathcal{E}^{h}}}$. The first two effectivities
compare the error estimates to the exact QoI discretization error
$\mathcal{E}$, while the second two effectivities compare the error estimates
to the QoI discretization error between the coarse and fine spaces,
$\mathcal{E}^{h}$. As discussed in the previous section, the effectivity
$\nicefrac{{\eta_{2}}}{{\mathcal{E}^{h}}}$ should evaluate identically to $1$.
In contrast, the effectivity $\nicefrac{{\eta_{1}}}{{\mathcal{E}^{h}}}$
illustrates the effect that neglecting the linearization errors
$\boldsymbol{E}^{\mathcal{R}}_{L}$ and $\mathcal{E}^{\mathcal{J}}_{L}$ has on
recovering the two-space QoI error $\mathcal{E}^{h}$. If these linearization
errors become small as $|\Omega^{e}_{H}|\to 0$, as one would hope, then the
effectivity $\nicefrac{{\eta_{1}}}{{\mathcal{E}^{h}}}$ should approach $1$.
Lastly, if the error estimates are _effective_ , then the effectivity indices
with respect to the exact error $\mathcal{E}$ will tend towards $1$ as
$|\Omega^{e}_{H}|\to 0$.

In addition to computing the effectivity indices described above, we also
compute _corrected_ functional approximations using the error estimates. These
corrected functional values take the form
$\mathcal{J}^{H}(\boldsymbol{u}^{H})+\eta_{1}$ for the traditional error
estimate $\eta_{1}$ and $\mathcal{J}^{H}(\boldsymbol{u}^{H})+\eta_{2}$ for the
newly proposed estimate $\eta_{2}$. The corrected functional values should
provide a more accurate representation of the QoI, which we measure by
computing their errors with respect to the exact solution
$\mathcal{J}(u)-(\mathcal{J}^{H}(\boldsymbol{u}^{H})+\eta_{1})$ and
$\mathcal{J}(u)-(\mathcal{J}^{H}(\boldsymbol{u}^{H})+\eta_{2})$. If the
correction improves the accuracy of the QoI approximation, then these computed
errors will converge at a faster rate than the exact QoI discretization error,
$\mathcal{J}(u)-\mathcal{J}^{H}(\boldsymbol{u}^{H})$. Figures 5, 6, and 7
illustrate the asymptotic behaviors of the error estimates $\eta_{1}$ and
$\eta_{2}$ for the QoIs $\mathcal{J}_{2}(u),\mathcal{J}_{3}(u)$, and
$\mathcal{J}_{4}(u)$, respectively. For each figure, the plot on the left
illustrates the behavior of the effectivity indices and the plot on the right
illustrates the behavior of computed error metrics as $|\Omega^{e}_{H}|\to 0$.

Figure 5: Asymptotic behavior of the estimates $\eta_{1}$ and $\eta_{2}$ for
the QoI $\mathcal{J}_{2}(u)$ for the manufactured solution when
$\alpha=10^{-2}$ on a sequence of uniformly refined meshes.

Figure 5 illustrates that both estimates $\eta_{1}$ and $\eta_{2}$ are
effective as the mesh size goes to zero for the QoI $\mathcal{J}_{2}(u)$ for
the chosen problem of interest, with $\nicefrac{{\eta_{1}}}{{\mathcal{E}}}\to
1$ and $\nicefrac{{\eta_{2}}}{{\mathcal{E}}}\to 1$ as $|\Omega^{e}_{H}|\to 0$.
However, the plot on the left of Figure 5 illustrates that the newly proposed
estimate $\eta_{2}$ is more effective than that of the traditional estimate
$\eta_{1}$ at coarse mesh resolutions. Additionally, the plot on the left of
Figure 5 demonstrates that both estimates result in a corrected functional
evaluation that converges at a faster rate than that of
$\mathcal{J}_{2}^{H}(\boldsymbol{U}^{H})$. In fact, the rate of this corrected
functional is nearly identical for the two estimates, but the newly proposed
estimate $\eta_{2}$ results in a constant shift towards a lower error and, as
a result, a more accurate approximation of the true functional value
$\mathcal{J}_{2}(u)$.

Figure 6: Asymptotic behavior of the estimates $\eta_{1}$ and $\eta_{2}$ for
the QoI $\mathcal{J}_{3}(u)$ for the manufactured solution when
$\alpha=10^{-2}$ on a sequence of uniformly refined meshes.

Figure 7: Asymptotic behavior of the estimates $\eta_{1}$ and $\eta_{2}$ for
the QoI $\mathcal{J}_{4}(u)$ for the manufactured solution when
$\alpha=10^{-2}$ on a sequence of uniformly refined meshes.

However, Figures 6 and 7 surprisingly illustrate that only the newly proposed
estimate $\eta_{2}$ is effective as the mesh size goes to zero for the QoIs
$\mathcal{J}_{3}(u)$ and $\mathcal{J}_{4}(u)$ for the chosen problem. In fact,
the error approximated by the traditional adjoint-weighted residual estimate
$\eta_{1}$ incurs a sign error. As a result, the corrected functional using
$\eta_{1}$ is actually _less accurate_ than that of the coarse-space
functional evaluation $\mathcal{J}^{H}(\boldsymbol{u}^{H})$, as seen by the
right-hand plots of Figures 6 and 7. In contrast, the newly proposed estimate
$\eta_{2}$ still results in a corrected functional evaluation that converges
at a faster rate than the coarse-space functional evaluation for the QoIs
$\mathcal{J}_{3}(u)$ and $\mathcal{J}_{4}(u)$.

$n_{el}$ | $||\boldsymbol{E}^{\mathcal{R}}_{L}||_{2}$
---|---
192 | 8.54316e+01
768 | 1.76089e+01
3072 | 3.05150e+00
12287 | 4.26147e$-$01
Table 3: Residual linearization errors for the manufactured solution with
$\alpha=10^{-2}$ at different mesh resolutions.

This leads to an insight of the present work: _the effectiveness of the
traditional adjoint-weighted residual can be QoI dependent_. That is to say,
it appears that the QoI linearization error $\mathcal{E}^{\mathcal{J}}_{L}$
can have a large impact on the effectiveness of the estimate $\eta_{1}$. We
speculate that this effect is typically unobserved in common goal-oriented
error estimation settings, where QoIs such as lift or drag over an airfoil are
_linear_ functionals, such that $\mathcal{E}^{\mathcal{J}}_{L}=0$. Further,
for our specific PDE of interest, the stark contrast in the asymptotic
behavior of the estimate $\eta_{1}$ for the QoI $\mathcal{J}_{2}(u)$ and the
QoIs $\mathcal{J}_{3}(u)$ and $\mathcal{J}_{4}(u)$ indicates that the residual
linearization error $\boldsymbol{E}^{\mathcal{R}}_{L}$ is not the culprit for
this loss of effectivity, since this term is independent of the chosen QoI.
That is to say, as the mesh size decreases,
$\boldsymbol{E}^{\mathcal{R}}_{L}\to\boldsymbol{0}$, which is confirmed in
Table 3.

####  6.1.5 Mesh Adaptivity

In this section, we investigate the behavior of the nonlinear functional
quantities and the use of the localized error estimates $\eta^{i}_{1}$ and
$\eta^{i}_{2}$, $i=1,2,\dots,n_{vtx}$, in the adaptive procedure described in
Section 5. For all examples, we begin with an initial mesh with $192$
elements, as shown in Figure 1 and perform five iterations of the process:

$\text{solve primal PDE}\rightarrow\text{solve adjoint
PDE}\rightarrow\text{estimate error}\rightarrow\text{adapt mesh}$

where either the traditional adjoint-weighted residual error estimate
$\eta_{1}$ or the newly proposed error estimate $\eta_{2}$ is used to drive
mesh modification. During each adaptive iteration, the mesh size field is
specified according to equation (35) so that the desired number of elements
$T$ is twice the number of elements in the previous mesh, and we evaluate the
QoI $\mathcal{J}^{H}(\boldsymbol{u}^{H})$ on the coarse space and measure its
error $\mathcal{E}$ with respect to the exact QoI value $\mathcal{J}(u)$.
Additionally, we compute the effectivity of both error estimates,
$\nicefrac{{\eta_{1}}}{{\mathcal{E}}}$ and
$\nicefrac{{\eta_{2}}}{{\mathcal{E}}}$, respectively. Figures 8, 9, and 10
illustrate the behavior of the adaptive schemes when $\eta_{1}$ and $\eta_{2}$
are chosen to drive mesh adaptivity for the QoIs
$\mathcal{J}_{2}(u),\mathcal{J}_{3}(u),$ and $\mathcal{J}_{4}(u)$,
respectively. The plot on the left of these figures shows the effectivities of
the estimates $\eta_{1}$ and $\eta_{2}$ as the mesh is adapted. The plot on
the right demonstrates the convergence behavior of the error $\mathcal{E}$ as
the mesh is adapted. Additionally, the plot on the right of these figures
compares the convergence behavior of the adaptive schemes to that of uniform
refinement.

Figure 8: Behavior of the adaptive scheme using the estimates $\eta_{1}$ and
$\eta_{2}$ for the QoI $\mathcal{J}_{2}(u)$ for the manufactured solution when
$\alpha=10^{-2}$.

Figure 8 illustrates that both estimates $\eta_{1}$ and $\eta_{2}$ are
effective as the mesh is adapted according to their localized indicators (33)
and (34), respectively, where the newly proposed estimate $\eta_{2}$ is
slightly more accurate at coarser mesh resolutions. Despite this fact,
however, the right-hand side of Figure 8 illustrates that the localizations of
each estimate lead to adaptive schemes that are very similar in terms of error
convergence, while adapting based on the newly proposed error estimate
$\eta_{2}$ only marginally improves the accuracy of the QoI when contrasted to
adapting based on the traditional adjoint-weighted residual error estimate
$\eta_{1}$. Both adaptive schemes, though, considerably outperform uniform
refinement in terms of reducing the error at a fixed number of degrees of
freedom in the problem.

Figure 9: Behavior of the adaptive scheme using the estimates $\eta_{1}$ and
$\eta_{2}$ for the QoI $\mathcal{J}_{3}(u)$ for the manufactured solution when
$\alpha=10^{-2}$.

Figure 10: Behavior of the adaptive scheme using the estimates $\eta_{1}$ and
$\eta_{2}$ for the QoI $\mathcal{J}_{4}(u)$ for the manufactured solution when
$\alpha=10^{-2}$.

In contrast, the left-hand plots of Figures 9 and 10 highlight that only the
newly proposed estimate $\eta_{2}$ is effective for the QoIs
$\mathcal{J}_{3}(U)$ and $\mathcal{J}_{4}(u)$ as adaptive iterations are
performed. Again, the traditional adjoint-weighted residual estimate
$\eta_{1}$ greatly under-predicts the QoI discretization error for these two
quantities of interest, which results in an effectivity
$\nicefrac{{\eta_{1}}}{{\mathcal{E}}}$ close to $0$. The right-hand plot of
Figure 9 also illustrates that driving mesh adaptivity with the estimate
$\eta_{2}$ leads to a significantly more accurate adaptive scheme for the QoI
$\mathcal{J}_{3}(u)$ when compared to driving adaptivity with the estimate
$\eta_{1}$ or when compared to uniform refinement. The right-hand plot of
Figure 10 tells a similar story for the QoI $\mathcal{J}_{4}(u)$, however, the
improvement of the adaptive scheme using the estimate $\eta_{2}$ over the
estimate $\eta_{1}$ is not as pronounced for this QoI.

This leads to another insight of the present work: _incorporating localized
information about the linearization errors $\boldsymbol{E}^{\mathcal{R}}_{L}$
and $\mathcal{E}^{\mathcal{J}}_{L}$ into an adaptive procedure can lead to
more optimal meshes in certain scenarios_. For our particular chosen PDE, the
fact that $\eta_{1}$ and $\eta_{2}$ lead to adaptive schemes that provide
quite similar behavior for the QoI $\mathcal{J}_{2}(u)$ but not for
$\mathcal{J}_{3}(u)$ and $\mathcal{J}_{4}(u)$ suggests that localizing the QoI
linearization $\mathcal{E}^{\mathcal{J}}_{L}$ may provide more benefit than
localizing the residual linearization error
$\boldsymbol{E}^{\mathcal{R}}_{L}$. We have presently made no effort to
decouple these effects, however, as an avenue for future work one could
investigate the approximation
$\mathcal{E}^{h}\approx-\boldsymbol{z}^{*}\cdot\boldsymbol{R}^{h}(\boldsymbol{u}^{H}_{h})$,
which does not account for the residual linearization error
$\boldsymbol{E}^{\mathcal{R}}_{L}$. In general, the interaction between the
two linearization error terms will depend on the structure of the primal PDE
and the chosen QoI.

####  6.1.6 A Problem with Gradient Singularities

In this section, we consider the model problem (36) over the domain $\Omega$
shown in Figure 1 with a forcing function $f$ chosen to be $f=100$ and the
scaling parameter $\alpha=10^{-2}$. For context, Newton’s method takes $4$
iterations for the solution of the primal problem on the coarse space and the
fine space for all results shown below. The choice of forcing function $f$ and
domain $\Omega$ leads to gradient singularities in the solution $u$ at each
interior corner of the domain $\Omega$. In a general sense, without accurately
resolving these singularities, so-called “pollution” error [8] will affect the
accuracy of the finite element approximation to the QoI throughout the domain.
Figure 11 illustrates the solution $u$ and its gradient magnitude $|\\!|\nabla
u|\\!|$ of the model problem with the chosen constant forcing function.

Figure 11: The solution with gradient singularities (left), the $L_{2}$ norm
of the solution gradient (right).

We investigate the behavior of all chosen functional quantities
$\mathcal{J}_{i}(u),i=1,2,3,4,$ and the use of the localized error estimates
$\eta^{i}_{1}$ and $\eta^{i}_{2},i=1,2,\dots,n_{vtx}$, in the adaptive process
described in Section 5. In general, the exact values of the QoIs are not
quantifiable analytically. To approximate these exact values, the primal
problem is solved on a mesh with about one half million degrees of freedom
using the fine space $\mathcal{V}^{h}$. This mesh is finer at every spatial
location than the final adapted meshes in the following results, and the
reference values computed by this process for each QoI are shown in B.

Figure 12: Convergence behavior for the example problem with gradient
singularities using uniform refinement and adaptivity driven by the estimates
$\eta_{1}$ and $\eta_{2}$ for the $\mathcal{J}_{1}(u)$ (left) and
$\mathcal{J}_{2}(u)$ (right).

For all examples, we begin with an initial mesh with $192$ elements, as shown
in Figure 1 and perform five iterations of the process:

$\text{solve primal PDE}\rightarrow\text{solve adjoint
PDE}\rightarrow\text{estimate error}\rightarrow\text{adapt mesh}$

where either the traditional adjoint-weighted residual error estimate
$\eta_{1}$ or the newly proposed estimate $\eta_{2}$ is used to drive mesh
modification. At each adaptive iteration, the mesh size field is specified
according to equation (35) so that the target number of elements $T$ in the
resultant mesh is twice the number of elements in the current mesh. At each
mesh instance, we evaluate the chosen QoI
$\mathcal{J}^{H}(\boldsymbol{u}^{H})$ on the coarse space and measure its
error $\mathcal{E}$ with respect to the reference value that approximates its
exact value $\mathcal{J}(u)$.

Figure 13: Convergence behavior for the example problem with gradient
singularities using uniform refinement and adaptivity driven by the estimates
$\eta_{1}$ and $\eta_{2}$ for the $\mathcal{J}_{3}(u)$ (left) and the
$\mathcal{J}_{4}(u)$ (right).

Figure 12 illustrates the convergence behavior of the QoI discretization error
when the mesh is uniformly refined and when $\eta_{1}$ and $\eta_{2}$ are
chosen to drive mesh adaptivity for $\mathcal{J}_{1}(u)$ and
$\mathcal{J}_{2}(u)$. Figure 13 similarly illustrates this behavior for
$\mathcal{J}_{3}(u)$ and $\mathcal{J}_{4}(u)$. Interestingly, Figure 12 shows
that both estimates, $\eta_{1}$ and $\eta_{2}$, result in adaptive convergence
histories with identical errors. This is because both schemes lead to
identical meshes at each adaptive iteration, as illustrated in Figure 14,
where the resultant mesh of the third, fourth, and fifth adaptive iteration is
shown. In addition, the plot on the left of Figure 12 illustrates that both
estimates, $\eta_{1}$ and $\eta_{2}$, only begin to provide improved
convergence behavior at finer mesh resolutions, whereas the plot on the right
of Figure 12 demonstrates that both estimates improve the QoI convergence
behavior, where the error is reduced by over half an order of magnitude at the
final mesh resolution of around $10,000$ degrees of freedom. This may, in
part, be explained by the fact that $\mathcal{J}_{1}(u)$ is a _global_ QoI,
while $\mathcal{J}_{2}(u)$ is restricted to only a localized portion of the
domain.

Figure 14: Sequence of meshes obtained for the example problem with gradient
singularities for $\mathcal{J}_{2}(u)$ using either the estimate $\eta_{1}$ or
$\eta_{2}$ at the third (left), fourth (center), and fifth (right) adaptive
iterations.

Figure 13 highlights that, for the QoIs $\mathcal{J}_{3}(u)$ and
$\mathcal{J}_{4}(u)$, the error convergence behavior using either the estimate
$\eta_{1}$ or $\eta_{2}$ is improved when compared to uniform refinement, but
that using the estimate $\eta_{2}$ provides a further reduction in error when
compared to $\eta_{1}$. This result is consistent with the behavior seen in
the previous manufactured solution section. Additionally, Figures 15 and 16
compare the final meshes obtained from the iterative adaptive process when
considering $\eta_{1}$ or $\eta_{2}$ for the QoIs $\mathcal{J}_{3}(u)$ and
$\mathcal{J}_{4}(u)$, respectively. Here, we make several remarks. First, for
both QoIs, both estimates heavily refine the gradient singularity nearest to
the local QoI subdomain, refine to a lesser degree the next two nearest
gradient singularities, and do not resolve at all the furthest gradient
singularity. Second, comparing the image on the right to the image on the left
of these figures highlights that both estimates lead to qualitatively distinct
meshes. This can, in part, be explained by the fact that the newly proposed
estimate $\eta_{2}$ includes linearization errors in its localization, while
the traditional adjoint-weighted residual $\eta_{1}$ does not.

Figure 15: Final adapted meshes obtained for the example problem with gradient
singularities for $\mathcal{J}_{3}(u)$ using the estimate $\eta_{1}$ (left)
and the estimate $\eta_{2}$ (right).

Figure 16: Final adapted meshes obtained for the example problem with gradient
singularities for the QoI $\mathcal{J}_{4}(u)$ using the estimate $\eta_{1}$
(left) and the estimate $\eta_{2}$ (right).

###  6.2 Finite Deformation Elasticity

In this section, we consider finite deformation elasticity with a neo-Hookean
material model in three spatial dimensions. Let the domain boundary $\Gamma$
be decomposed into $\Gamma_{G}\subseteq\Gamma$ and
$\Gamma_{H}=\Gamma\setminus\Gamma_{G}$. Let
$\boldsymbol{X}\in\Omega\subset\mathbb{R}^{3}$ denote a point in a reference
configuration that transforms to a point $\boldsymbol{x}\in\Omega_{t}$ in a
deformed configuration $\Omega_{t}\subset\mathbb{R}^{3}$ after the reference
domain $\Omega$ undergoes some deformation. Let
$\boldsymbol{u}=\boldsymbol{x}-\boldsymbol{X}$ denote the displacement. Let
$\boldsymbol{F}:=\boldsymbol{I}+\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{X}}$
denote the deformation gradient, where $\boldsymbol{I}$ denotes the second-
order identity tensor. Let $J:=\det(\boldsymbol{F})$ denote the determinant of
the deformation gradient. In this context, the balance of linear momentum
without inertial and body forces can be written as

$\begin{cases}\begin{aligned}
-\nabla\cdot\boldsymbol{P}&=\boldsymbol{0},\quad&&\boldsymbol{X}\in\Omega,\\\
\boldsymbol{u}&=\boldsymbol{G},&&\boldsymbol{X}\in\Gamma_{G},\\\
\boldsymbol{P}\cdot\boldsymbol{N}&=\boldsymbol{0},&&\boldsymbol{X}\in\Gamma_{H}.\end{aligned}\end{cases}$
(43)

Here, $\boldsymbol{G}$ denotes an externally applied displacement,
$\boldsymbol{N}$ denotes the unit outward normal to the boundary $\Gamma_{H}$,
$\boldsymbol{P}:=J\boldsymbol{\sigma}\boldsymbol{F}^{-T}$ denotes the first
Piola-Kirchhoff stress tensor, and $\boldsymbol{\sigma}$ denotes the Cauchy
stress tensor. We consider a neo-Hookean constitutive model determined by the
relationship

$\boldsymbol{\sigma}=\mu
J^{-\nicefrac{{5}}{{3}}}\text{dev}(\boldsymbol{F}\boldsymbol{F}^{T})+\frac{\kappa}{2}(J-1/J)\boldsymbol{I},$
(44)

where $\kappa$ denotes the bulk modulus, $\mu$ denotes the shear modulus, and
$\text{dev}(\cdot)$ denotes the deviatoric portion of a second-order tensor.
Placing the problem (43) in weak form yields the semilinear form

$\mathcal{N}(\boldsymbol{u};\boldsymbol{w}):=\int_{\Omega}\boldsymbol{P}(\boldsymbol{u}):\nabla\boldsymbol{w}\,\text{d}\Omega,$
(45)

and the linear functional $\mathcal{L}(\boldsymbol{w})=0$.

Figure 17: The domain $\Omega$ and its initial mesh for the nonlinear
elasticity example (left), the deformed domain $\Omega_{t}$ after loading
(exaggerated by a factor of 5) with the von Mises stress $\sigma_{vm}$ (MPa)
plotted (center), and the norm of the adjoint solution vector $\boldsymbol{z}$
(right).

As a model problem, we consider the domain $\Omega$ shown in Figure 17 which
corresponds to a cube of dimensions $5\text{mm}\times 5\text{mm}\times
5\text{mm}$ with a cube of size $2.5\text{mm}\times 2.5\text{mm}\times
2.5\text{mm}$ missing from its uppermost corner. Displacements are set to
zero, $\boldsymbol{u}=\boldsymbol{0}$, on the minimal $y$-face of the geometry
and the problem is driven by a prescribed $y$-displacement of
$u_{y}=0.1\text{mm}$, which corresponds to a $2\%$ strain in the
$y$-direction. We consider a material with elastic modulus $E=192.7$GPa and
Poisson’s ratio $\nu=0.27$, from which the bulk and shear moduli can be
determined by the relationships $\kappa=\frac{E}{3(1-2\nu)}$ and
$\mu=\frac{E}{2(1+\nu)}$, respectively.

We again consider solving the finite element problem corresponding to the
governing equations on a coarse space with piecewise linear Lagrange basis
functions and on a fine space with piecewise quadratic Lagrange basis
functions. As a quantity of interest, we consider the von Mises stress
$\sigma_{vm}$ integrated over the domain

$\mathcal{J}_{vm}(\boldsymbol{u}):=\int_{\Omega}\,\sigma_{vm}\;\text{d}\Omega,$
(46)

where $\sigma_{vm}:=\sqrt{\frac{3}{2}s_{ij}s_{ij}}$ (with summation on
repeated indices implied), and $\boldsymbol{s}$ denotes the deviatoric portion
of the Cauchy stress tensor, as defined by $\boldsymbol{s}:=\mu
J^{-\nicefrac{{5}}{{3}}}\text{dev}(\boldsymbol{F}\boldsymbol{F}^{T})$.

Figure 18: Behavior of the adaptive scheme using the estimates $\eta_{1}$ and
$\eta_{2}$ for the QoI $\mathcal{J}_{vm}(u)$ for the finite deformation
elasticity example.

We consider the initial mesh of the domain $\Omega$ shown in the left-most
image of Figure 17 and perform $10$ iterations of the process:

$\text{solve primal PDE}\rightarrow\text{solve adjoint
PDE}\rightarrow\text{estimate error}\rightarrow\text{adapt mesh}$

where either the traditional adjoint-weighted residual $\eta_{1}$ estimate or
the newly proposed estimate $\eta_{2}$ is used to drive mesh adaptivity. At
each adaptive iteration, the mesh size field is set according to equation (35)
so that the target number of elements $T$ in the resultant mesh is twice the
number of elements in the current mesh. At each mesh instance, we evaluate the
QoI $\mathcal{J}^{H}_{e}(\boldsymbol{u}^{H})$ on the space and measure its
error with respect to a reference approximation for the exact value of the
QoI, compute an approximation to the exact error $\mathcal{E}$ using this
reference value, and compute the effectivity indices $\eta_{1}/\mathcal{E}$
and $\eta_{2}/\mathcal{E}$. To find the reference value, the primal problem is
solved on a mesh with $17.3$ million degrees of freedom using the fine space
$\mathcal{V}^{h}$, where this mesh is finer at every spatial location than the
final adapted meshes used in the results below, which contain about $300,000$
degrees of freedom. The reference value is found to be
$\mathcal{J}_{vm}(\boldsymbol{u})\approx 3.731860\text{GPa}$.

Figure 18 illustrates the behavior of the error estimates $\eta_{1}$ and
$\eta_{2}$ as adaptive iterations are performed. The left-hand side of Figure
18 illustrates illustrates that the newly proposed estimate $\eta_{2}$ is more
effective for the QoI $\mathcal{J}_{vm}(\boldsymbol{u})$ than the traditional
adjoint-weighted residual estimate $\eta_{1}$. Note that for this problem the
error in the QoI is negative, which is why the effectivity index for
$\eta_{1}$ is consistently greater than one. This is in contrast to the
examples in the previous section. The right-hand side of Figure 18 highlights
once again that the localizations of each estimate lead to adaptive schemes
that are different in terms of convergence behavior. Specifically, driving
mesh adaptivity with the estimate $\eta_{2}$ leads to a significantly more
accurate adaptive scheme for the QoI $\mathcal{J}_{vm}(\boldsymbol{u})$ when
compared to driving adaptivity with the estimate $\eta_{1}$ or when compared
to uniform refinement.

Figure 19: Meshes using the $\eta_{2}$ adative scheme: the adapted mesh at the
fifth adaptive iteration (left), the adapted mesh at the final adaptive
iteration (center), a close-up and cutaway of the final adaptive mesh (right).

##  7 Conclusions

In this paper, we have considered goal-oriented _a posteriori_ error
estimation for nonlinear QoIs when solving nonlinear PDEs with continuous
Galerkin finite element methods. As a starting point, we have reviewed the
derivation of a traditional two-level, or discrete, adjoint-weighted residual
approach to approximate QoI discretization errors. Inherent in this derivation
is the linearization of both the nonlinear PDE residual and the nonlinear
functional of interest, where higher-order terms are subsequently disregarded.
We refer to these higher-order terms as linearization errors, as they
introduce an error into the error estimate. By considering these linearization
errors, we have derived a novel two-level adjoint-based error estimate that
exactly represents the QoI discretization error between two function spaces at
the expense of increased computational complexity.

We have motivated our interest in this novel error estimate with several _a
priori_ justifications and, in particular, have reasoned that the traditional
adjoint-weighted residual error estimate can under-predict the QoI
discretization error for nonlinear QoIs. We have demonstrated that this under-
prediction can in fact occur for two example problems with a nonlinear
Poisson’s equation in two spatial dimensions and that this behavior can be QoI
dependent. For these problems, we have demonstrated that the newly proposed
error estimate does not suffer from this under-prediction and is effective for
each investigated QoI. Additionally, we have demonstrated that the newly
proposed error estimate is more effective than the traditional adjoint-
weighted residual estimate for a finite deformation elasticity example in
three spatial dimensions.

The solution of a nonlinear scalar problem is required for the newly proposed
error estimate. We have provided the analytic solution to this problem when
quadratic QoIs are considered and have outlined how Newton’s method can be
used to solve the general nonlinear scalar problem. We have demonstrated
heuristically that the analytic solution for quadratic QoIs is a good initial
guess for Newton’s method when more general nonlinear QoIs are considered.

We have shown how the newly proposed error estimate can be localized and
incorporated into an iterative adaptive mesh algorithm. We have demonstrated
that the newly proposed estimate can lead to more accurate approximations of a
nonlinear QoI with fewer degrees of freedom when compared to uniform
refinement and traditional adjoint-based approaches by considering a nonlinear
Poisson’s example problem and a finite deformation elasticity example problem.

Orthogonally, we have demonstrated that the evaluation of the PDE residual
linearization error can lead to a convenient verification check for the
adjoint solution used in traditional adjoint-weighted residual error
estimates.

As an avenue for future work, we have suggested incorporating known recovery
techniques into the novel estimate to reduce its computational complexity. One
could also investigate the effect that the choice of coarse and fine spaces
has on the newly proposed error estimate, as we have only considered piecweise
linear and quadratic Lagrange basis functions, respectively, in the present
work.

##  8 Acknowledgments

Supported by the Advanced Simulation and Computing program at Sandia National
Laboratories, a multimission laboratory managed and operated by National
Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary
of Honeywell International Inc. for the U.S. Department of Energy’s National
Nuclear Security Administration under contract DE-NA0003525. This article
describes objective technical results and analysis. Any subjective views or
opinions that might be expressed in the article do not necessarily represent
the views of the U.S. Department of Energy or the United States Government.

## References

  * [1] Behzad R Ahrabi, W Kyle Anderson, and James C Newman III.  An adjoint-based hp-adaptive stabilized finite-element method with shock capturing for turbulent flows.  Computer Methods in Applied Mechanics and Engineering, 318:1030–1065, 2017.
  * [2] Burak Aksoylu, Stephen D. Bond, Eric C. Cyr, and Michael Holst.  Goal-oriented adaptivity and multilevel preconditioning for the Poisson-Boltzmann equation.  Journal of Scientific Computing, 52:202–225, 2012.
  * [3] Frédéric Alauzet, Xiangrong Li, E. Seegyoung Seol, and Mark S. Shephard.  Parallel anisotropic 3D mesh adaptation by mesh modification.  Engineering with Computers, 21(3):247–258, 2006.
  * [4] Zaib Ali, Paul G Tucker, and Shahrokh Shahpar.  Optimal mesh topology generation for cfd.  Computer Methods in Applied Mechanics and Engineering, 317:431–457, 2017.
  * [5] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E Rognes, and Garth N Wells.  The fenics project version 1.5.  Archive of numerical software, 3(100), 2015.
  * [6] R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J. Cerveny, V. Dobrev, Y. Dudouit, A. Fisher, Tz. Kolev, W. Pazner, M. Stowell, V. Tomov, I. Akkerman, J. Dahm, D. Medina, and S. Zampini.  MFEM: A modular finite element methods library.  Computers & Mathematics with Applications, 81:42–74, 2021.
  * [7] Daniel Arndt, Wolfgang Bangerth, Marco Feder, Marc Fehling, Rene Gassmöller, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Peter Munch, Jean-Paul Pelteret, Simon Sticko, Bruno Turcksin, and David Wells.  The deal.II library, version 9.4.  Journal of Numerical Mathematics, 30(3):231–246, 2022.
  * [8] I. Babuška, T. Strouboulis, A. Mathur, and C.S. Upadhyay.  Pollution-error in the h-version of the finite-element method and the local quality of a-posteriori error estimators.  Finite Elements in Analysis and Design, 17(4):273–321, 1994.
  * [9] Roland Becker and Rolf Rannacher.  An optimal control approach to _a posteriori_ error estimation in finite element methods.  Acta Numerica, 10:1–102, 2001.
  * [10] Patrick J Blonigan and Eric J Parish.  Evaluation of dual-weighted residual and machine learning error estimation for projection-based reduced-order models of steady partial differential equations.  Computer Methods in Applied Mechanics and Engineering, 409:115988, 2023.
  * [11] Ramzy Boussetta, Thierry Coupez, and Lionel Fourment.  Adaptive remeshing based on a posteriori error estimation for forging simulation.  Computer Methods in Applied Mechanics and Engineering, 195(48-49):6626–6645, 2006.
  * [12] Vit Dolejší and Scott Congreve.  Goal-oriented error analysis of iterative Galerkin discretizations for nonlinear problems including linearization and algebraic errors, 2023.
  * [13] Geoff Donoghue and Masayuki Yano.  Spatio-stochastic adaptive discontinuous galerkin methods.  Computer Methods in Applied Mechanics and Engineering, 374:113570, 2021.
  * [14] Bernhard Endtmayer, Ulrich Langer, and Thomas Wick.  Two-side a posteriori error estimates for the dual-weighted residual method.  SIAM Journal on Scientific Computing, 42(1):A371–A394, 2020.
  * [15] Krzysztof J Fidkowski.  Output-based error estimation and mesh adaptation for unsteady turbulent flow simulations.  Computer Methods in Applied Mechanics and Engineering, 399:115322, 2022.
  * [16] Krzysztof J. Fidkowski and David L. Darmofal.  Review of output-based error estimation and mesh adaptation in computational fluid dynamics.  AIAA Journal, 49(4):673–694, 2011.
  * [17] O.A. González-Estrada, E. Nadal, J.J. Ródenas, P. Kerfriden, S.P.A. Bordas, and F.J. Fuenmayor.  Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery.  Computational Mechanics, 53(5):957–976, 2014.
  * [18] Brian N. Granzow, Assad A. Oberai, and Mark S. Shephard.  Adjoint-based error estimation and mesh adaptation for stabilized finite deformation elasticity.  Computer Methods in Applied Mechanics and Engineering, 337:263–280, 2018.
  * [19] Brian N. Granzow, Mark S. Shephard, and Assad A. Oberai.  Output-based error estimation and mesh adaptation for variational multiscale methods.  Computer Methods in Applied Mechanics and Engineering, 322:441–459, 2017.
  * [20] Michael Holst, Sara Pollock, and Yunrong Zhu.  Convergence of goal-oriented adaptive finite element methods for semilinear problems.  Computing and Visualization in Science, 17:43–63, 2015.
  * [21] Daniel A Ibanez, E Seegyoung Seol, Cameron W Smith, and Mark S Shephard.  Pumi: Parallel unstructured mesh infrastructure.  ACM Transactions on Mathematical Software (TOMS), 42(3):1–28, 2016\.
  * [22] Steven M. Kast.  An introduction to adjoints and output error estimation in computational fluid dynamics, 2017.
  * [23] Xiangrong Li, Mark S. Shephard, and Mark W. Beall.  3D anisotropic mesh adaptation by mesh modification.  Computer Methods in Applied Mechanics and Engineering, 194(48-49):4915–4950, 2005.
  * [24] Jorge Nocedal and Stephen J. Wright.  Numerical optimization.  Springer, 1999.
  * [25] Ricardo H Nochetto, Andreas Veeser, and Marco Verani.  A safeguarded dual weighted residual method.  IMA journal of Numerical Analysis, 29(1):126–140, 2009.
  * [26] J. Peraire and A.T. Patera.  Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement.  Studies in Applied Mechanics, 47:199–216, 1998.
  * [27] Serge Prudhomme and J. Tinsley Oden.  On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors.  Computer Methods in Applied Mechanics and Engineering, 176(1-4):313–331, 1999.
  * [28] Thomas Richter and Thomas Wick.  Variational localizations of the dual weighted residual estimator.  Journal of Computational and Applied Mathematics, 279:192–208, 2015\.
  * [29] J.J. Ródenas, M. Tur, F.J. Fuenmayor, and A. Vercher.  Improvement of the superconvergent patch recovery technique by the use of constraint equations: the SPR-C technique.  International Journal for Numerical Methods in Engineering, 70(6):705–727, 2007.
  * [30] G. Şimşek, X. Wu, K.G. van der Zee, and E.H. van Brummelen.  Duality-based two-level error estimation for time-dependent PDEs: Application to linear and nonlinear parabolic equations.  Computer Methods in Applied Mechanics and Engineering, 288:83–109, 2015.
  * [31] Michael K Sleeman and Masayuki Yano.  Goal-oriented model reduction for parametrized time-dependent nonlinear partial differential equations.  Computer Methods in Applied Mechanics and Engineering, 388:114206, 2022.
  * [32] David A. Venditti and David L. Darmofal.  Adjoint error estimation and grid adaptation for functional outputs: Application to quasi-one-dimensional flow.  Journal of Computational Physics, 164(1):204–227, 2000.
  * [33] David A. Venditti and David L. Darmofal.  Grid adaptation for functional outputs: Application to two-dimensional inviscid flows.  Journal of Computational Physics, 176(1):40–69, 2002.
  * [34] David A. Venditti and David L. Darmofal.  Anisotropic grid adaptation for functional outputs: Application to two-dimensional viscous flows.  Journal of Computational Physics, 187(1):22–46, 2003.
  * [35] Nils-Erik Wiberg, Fethi Abdulwahab, and Saulius Ziukas.  Enhanced superconvergent patch recovery incorporating equilibrium and boundary conditions.  International Journal for Numerical Methods in Engineering, 37(20):3417–3440, 1994.
  * [36] Olgierd Cecil Zienkiewicz and Jian Zhong Zhu.  The superconvergent patch recovery and _a posteriori_ error estimates. Part 1: The recovery technique.  International Journal for Numerical Methods in Engineering, 33(7):1331–1364, 1992.

##  Appendix A QoI Values for the Manufactured Solution

The QoI values for the manufactured solution and domain defined in Section
6.1.1 are:

$\displaystyle\mathcal{J}_{1}(u)$
$\displaystyle=\frac{64\pi^{2}(e^{\nicefrac{{5}}{{2}}}-1)^{2}(e^{5}+e^{\nicefrac{{5}}{{2}}}+1)}{e^{5}(16\pi^{2}+25)^{2}}\approx
2.57052599061823,$ (47) $\displaystyle\mathcal{J}_{2}(u)$
$\displaystyle=\frac{65536\pi^{6}(e^{15}+e^{\nicefrac{{45}}{{4}}}+e^{\nicefrac{{15}}{{4}}}+1)}{9e^{\nicefrac{{15}}{{2}}}(256\pi^{4}+4000\pi^{2}+5625)^{2}}\approx
2.64090163593838,$ $\displaystyle\mathcal{J}_{3}(u)$
$\displaystyle=\frac{32\pi^{4}(e^{\nicefrac{{5}}{{2}}}-1)^{2}(e^{5}+e^{\nicefrac{{5}}{{2}}}+1)}{25e^{5}(16\pi^{2}+25)}\approx
9.28106693871883\times 10^{1},$ $\displaystyle\mathcal{J}_{4}(u)$
$\displaystyle\approx 5.67945022.$

Here $\mathcal{J}_{4}(u)$ was approximated numerically using Mathematica.

##  Appendix B QoI Values for the Solution with Gradient Singularities

The QoI values for the problem definition defined in Section 6.1.6 are:

$\displaystyle\mathcal{J}_{1}(u)$ $\displaystyle\approx 6.540644835,$ (48)
$\displaystyle\mathcal{J}_{2}(u)$ $\displaystyle\approx 1.238067612\times
10^{1},$ $\displaystyle\mathcal{J}_{3}(u)$ $\displaystyle\approx
1.596007278\times 10^{2},$ $\displaystyle\mathcal{J}_{4}(u)$
$\displaystyle\approx 9.391778787.$