April 02, 2024

This work reviews goal-oriented a posteriori error control, adaptivity and solver control for finite element approximations to boundary and initial-boundary value problems for stationary and non-stationary partial differential equations, respectively.
In particular, coupled field problems with different physics may require simultaneously the accurate evaluation of several quantities of interest, which is achieved with multi-goal oriented error control. Sensitivity measures are obtained by solving an
adjoint problem. Error localization is achieved with the help of a partition-of-unity. We also review and extend theoretical results for efficiency and reliability by employing a saturation assumption. The resulting adaptive algorithms allow to balance
discretization and non-linear iteration errors, and are demonstrated for four applications: Poisson’s problem, non-linear elliptic boundary value problems, stationary incompressible Navier-Stokes equations, and regularized parabolic \(p\)-Laplace initial-boundary value problems. Therein, different finite element discretizations in two different software libraries are utilized, which are partially accompanied with open-source implementations on GitHub.

**Keywords:** Goal-oriented error control; multi-goal error control; adjoint problems; dual-weighted residual method; partial differential equations; adaptive finite element methods

This work is devoted to a goal-oriented a posteriori error control of single and multiple quantities of interest (QoI) and Adaptive Finite Element Methods (AFEM) for solving stationary and non-stationary, linear and non-linear partial differential equations (PDEs) and systems of PDEs. AFEM were certainly inspired by the pioneer work [1] of Babuška and Rheinboldt. Further important studies on a posteriori error controlled adaptive finite element methods are [2]–[6] to name a few. The first comprehensive monograph [7] by Verfürth was another AFEM milestone. We refer to the survey papers [8], [9], [10], and to the monographs [11], [12], [13], [14], [15], [16], [17], [18] for an overview of a posteriori error estimates and adaptive finite element techniques.

The governing discretization in this work is the finite element method [19]–[24], but isogeometric analysis could be used in a similar fashion [25]–[27]. We consider stationary settings as well as time-dependent cases, modeled as space-time problems and solved by space-time discretization methods; see, e.g., [28]–[36] and the references therein.

In the following, we review the important steps regarding the development of goal-oriented techniques. Starting with the Acta Numerica paper [8] by Eriksson, Estep, Hansbo and Johnson as a survey on adaptive methods, Becker and Rannacher proposed an automated procedure for self-controlled adaptivity with the dual-weighted residual method (DWR) in their Acta Numerica paper [9]. This was followed by the book [13]. Further studies of the early developments include [37]–[46]. An important step towards time-dependent problems in space-time formulations with full space-time adaptivity was the work by Schmich and Vexler [47]. These previously mentioned studies concentrated on single goal-oriented error estimation. In the year \(2003\), Hartmann and Houston proposed goal-oriented a posteriori error estimation for multiple quantities of interest [48], followed by Hartmann’s paper [49] and shortly later [50]. We started ourselves with multigoal-oriented error estimation with [51], [52].

Conceptional developments on goal-oriented error estimation include a safe-guarded method established in [53] and guaranteed bounds derived in [54]. Abstract analyses and reformulations were presented in [55], [56]. First convergence rates of goal-oriented error estimation were obtained in [57], weighted marking [58], bounding techniques for goal-oriented error estimates for elliptic problems [59], goal oriented flux reconstruction [60], goal-oriented error estimation for the fractional-step-theta scheme [61], a partition-of-unity localization including effectivity estimates for single goal functionals [62], partition-of-unity localization for multiple goal functionals [51], space-time partition-of-unity localization [63], [64], equilibrated flux [65], linearization errors [66]. Theoretical work showing optimality of adaptive algorithms was summarized in four axioms of adaptivity in [67] and dates back to residual-based a posteriori error estimates in [2]. First rigorous convergence results for goal-oriented estimators go back to [57], [58], then [55], [68], [69], and recent findings are [70]–[72]. We also mention our own prior work on such theoretical advancements that include effective estimates with common upper bounds of the dual-weighted residual estimator and the error indicators [62]. Efficiency and reliability estimates for the dual-weighted residual method using a saturation assumption were shown in [73], and smart algorithms switching between solving high-order adjoint problems or using interpolations were derived in [74].

A key step in the development of the dual-weighted residual method was the introduction of a partition-of-unity (PU) localization by Richter and Wick [62]; a prototype open-source implementation based on deal.II [75]–[77] can be found on GitHub^{1}. For the general idea of the PU-FEM, we refer the reader to [78], [79]. The PU-DWR method opened the way for addressing goal-oriented error control in practical applications of
time-dependent, non-linear, coupled PDEs and multiphysics applications in a much more convenient way. Along with the PU-DWR method, an indicator index was introduced, which measures the quality of error indicators. The PU-DWR method was extended to
multigoal-oriented error control in [51]. The extension of PU-DWR to space-time goal-oriented error estimation on fully unstructured simplicial decompositions
of the space-time cylinder was done in [63] and a \(cGdG\) (continuous Galerkin in space, discontinuous
Galerkin in time) discretization was done for parabolic problems in [64], including open-source developments in Zenodo [80], [81] and GitHub^{2},^{3}, and finally, space-time PU-DWR for the incompressible Navier-Stokes equations in [82], again with codes on GitHub^{4}.

Single- and multigoal-oriented error control and adaptivity has found numerous applications. These include Poisson’s problem [44], [46], [51], [83], stationary linear elasticity [45], space-time adaptivity for parabolic equations [47],
space-time elasticity and the wave equation [84], [85], multi-rate discretizations of coupled
parabolic/hyperbolic problems [86], \(p\)-Laplace problems [52], [63], [74], [87], elasto-plasticity [88]–[90] and visco-plasticity [91], incompressible flow [73], [74], [82], [92]–[96], transport proplems [97], transport problems with coupled flow [98]–[100], Boussinesq equations and coupled flow with temperature [101]–[103], reactive flows [104]–[106], uncertain inputs [107], Maxwell’s equations [108], elliptic eigenvalue problems [109], [110], modeling errors [41], [111], [112], applications to polygonal meshes [113], conforming and nonconforming approximations with inexact solvers [114], anisotropic mesh refinement [115], heterogeneous multiscale problems [116]–[120], fluid-structure interaction [121]–[132], phase-field slits and fracture [133]–[135],
sea ice simulations [91], optimal control [136]–[145], obstacle/contact
problems and variational inequalities [134], [146]–[151], finite cell methods [152], [153], balancing discretization and iteration errors [52], [87], [154]–[157], financial mathematics [158], goal-oriented model order reduction [159]–[163], neural network enhanced dual-weighted residual technologies [164]–[166], adaptive multiscale predictive modeling [167], and open-source software developments for goal-oriented error estimation [[168]]^{5}, [[169]]^{6},
[[64]; [170]]^{7},^{8}.

In this work, we undertake an extensive review on single- and multigoal-oriented a posteriori error estimation and adaptivity. Our results include classical proofs for educational purposes. But we also add new findings: first, we establish efficiency and reliability estimates with only one saturation assumption in which the strengthened condition from [73] is no longer necessary. Second, we provide details on non-standard discretizations such as non-consistent and non-conformal methods. Third, the code of the basic PU-DWR method [62] is published open-source on GitHub. Fourth, in extension of [63], the space-time PU-DWR method is investigated for multiple goals with singularities. Certain results are illustrated with the elliptic model problem. Our main theoretical results include estimates and algorithms for balancing discretization and non-linear iteration errors. Our algorithms are illustrated with newly designed numerical tests, not published in the literature elsewhere, that have both research and educational character in order to explain the mechanisms of goal-oriented error control and mesh adaptivity. These results include stationary, non-linear situations, incompressible flow, and a space-time \(p\)-Laplace problem.

The outline of this work is as follows. In Section 2, the notation and abstract setting are introduced. Then, in Section 3, single goal-oriented error control is discussed. Section 4 focuses on error estimation for non-standard discretizations. Next, in Section 5, multigoal-oriented error estimators are explained. In Section 6, three applications are considered, and computationally analyzed. The work is concluded in Section 7, and some current research directions are outlined.

Throughout this work, let \(d\) be the space dimension. We denote by \(Q:=\Omega\times I \subset \mathbb{R}^{d+1}\) the space-time domain (cylinder), where \(\Omega \subset \mathbb{R}^{d}\) is the bounded and Lipschitz spatial domain with the boundary \(\Gamma = \partial \Omega\), and \(I:=(0,T)\) denotes the temporal domain (time interval). Here, \(T>0\) is the end time value. Furthermore, we use the standard notations for Lebesgue, Sobolev, and Bochner spaces like \(L_p(\Omega)\), \(W_p^k(\Omega),\mathring{W}_p^k(\Omega), H^k(\Omega)=W_2^k(\Omega), H^k_0(\Omega)=\mathring{W}_2^k(\Omega)\), \(L_p(0,T;\mathring{W}_p^1(\Omega))\) etc.; see, e.g., [19], [20], [171].

Let \(U\) and \(V\) be reflexive Banach spaces with their dual spaces \(U^*\) and \(V^*\), respectively. We consider the abstract operator equation: Find \(u \in U\) such that \[\label{eq:32General32Modelproblem} \mathcal{A}(u)=0 \qquad \text{in }V^*,\tag{1}\] where \(\mathcal{A}:U \mapsto V^*\) represents a non-linear partial differential operator.

Let us assume \(U_h\) and \(V_h\) are conforming discrete spaces, i.e.we have the properties \(U_h \subseteq U\) and \(V_h \subseteq V\) and \(\operatorname{dim}(U_h) < \infty\) and \(\operatorname{dim}(V_h)< \infty\). With this, we can perform a Galerkin-Petrov discretization of the operator equation 1 yielding the discrete problem: Find \(u_h \in U_h\) such that \[\label{eq:32discrete32primal} \mathcal{A}(u_h)(v_h)=0 \qquad \forall v_h \in V_h.\tag{2}\] Once a basis is chosen, this leads to a non-linear system of finite equations. Afterwards, this system of non-linear equations can be solved with some non-linear solver, e.g. Picard’s iteration or Newton’s method [172].

One possible discretization technique is the Finite Element Method (FEM) [19]–[24]. In this section, we assume that \(\Omega\) is a polygonal Lipschitz domain. As an example for the conforming discrete subspaces from Section 2.1, we provide two possible finite element discretizations.

We can decompose our domain into shape-regular simplicial elements \(K \in \mathcal{T}_h\), where \(\bigcup_{K \in \mathcal{T}_h}\overline{K}=\overline{\Omega}\) and for all \(K,K' \in \mathcal{T}_h:\) \(K \cap K' = \emptyset\) if and only if \(K\not=K'\). Let \(\mathbb{P}_k(\hat{K})\) be the space of polynomials of total degree \(k\) on the reference domain \(\hat{K}\) and \[{ P_k:=\{v_h \in \mathcal{C}(\Omega): v_h|_K = v_h(\varphi_K(\cdot)) \in \mathbb{P}_k(\hat{K})\quad \forall K \in \mathcal{T}_h \} }\] the finite element space of continuous finite elements of degree \(k\) on \(\mathcal{T}_h\), where \(\varphi_K: \hat{K} \rightarrow K\) is a regular mapping of the reference element \(\hat{K}\) onto \(K\), e.g., a affine-linear or isoparametric mapping. For more information, we refer to [19], [20], [24].

Another possible way is to decompose our domain into hypercubal elements \(K \in \mathcal{T}_h\), where \(\bigcup_{K \in \mathcal{T}_h}\overline{K}=\overline{\Omega}\) and for all \(K,K' \in \mathcal{T}_h:\) \(K \cap K' = \emptyset\) if and only if \(K\not=K'\). Let \(\mathbb{Q}_k(\hat{K})\) be the tensor product space of polynomials of degree \(k\) on the reference domain \(\hat{K}\) and \[Q_k:=\{v_h \in \mathcal{C}(\Omega): v_h|_K = v_h(\varphi_K(\cdot)) \in \mathbb{Q}_k(\hat{K})\quad \forall K \in \mathcal{T}_h \}\] the finite element space of continuous finite elements of degree \(k\) on \(\mathcal{T}_h\), where \(\varphi_K: \hat{K} \rightarrow K\) is a regular mapping of the reference element \(\hat{K}\) onto \(K\), e.g., a multi-linear or isoparametric mapping. We refer to [19], [20], [24] for more information.

To facilitate adaptive mesh refinement and to avoid connecting elements, we use the concept of hanging nodes. Elements are allowed to have nodes that lie on the midpoints of the faces or edges of neighboring cells. In our implementations, at most, one hanging node is allowed on each face or edge. In three dimensions, this concept is generalized to subplanes and faces because we must deal with two types of lower manifolds. To enforce global continuity (i.e., global conformity), the degrees of freedom located on the interface between different refinement levels have to satisfy additional constraints. They are determined by interpolation of neighboring degrees of freedom. Therefore, hanging nodes do not carry any degrees of freedom. For more details on this, we refer to [23].

In this section, we first provide some background information and we explain the need for goal-oriented error estimation. Moreover, we also introduce and explain adjoint-based error estimation employing the so-called dual-weighted residual (DWR) method.

First, we address the purpose to construct error estimators and adaptive schemes. Error estimators allow to determine approximately the error between approximations and the (unknown) true solution of a given problem statement. One distinguishes:

a priori estimates (estimated before the actual approximate solution is known);

a posteriori error estimates after the approximate solution has been computed.

Based on these error estimations, the error may be further localized in order to ‘refine’ algorithms. Mostly, mesh refinement in space and/or time is of interest, but also model errors or iteration errors can be controlled and balanced [8], [9], [11], [13], [18]. These localized errors can be used to adaptively steer the algorithm, by enhancing the accuracy of the approximate solution, while keeping the computational cost reasonable.

In anticipation of practical realization and optimality of adaptive procedures, the basic algorithm reads:

**Algorithm 1** (AFEM - adaptive finite elements). *The basic algorithm for AFEM reads:*

**Solve**the differential equation on the current mesh \({\cal T}\);**Estimate**the error via a posteriori error estimation to obtain \(\eta\);**Mark**the elements by localizing the error estimator;**Refine/coarsen**the elements with the highest/lowest error contributions using a certain refinement strategy.

In the 70s, a priori and a posteriori error estimates were derived based on global norms, e.g., the \(L^2\)-norm, the \(H^1\)-norm or even classical \(C^0\)- and \(C^1\)-norms. With this, the entire solution is controlled and resulting adaptive schemes act correspondingly. However, in many applications, only certain parts of the numerical solution are of interest. Regarding the geometry (domains \(Q, \Omega, I\)), such parts can be subdomains \(\tilde{Q}\subset Q\), line evaluations \(\Gamma\subset \bar{Q}\) or even simply points \((x,t)\in Q\) at which solution information (values, derivatives) are evaluated. In case of ordinary differential equations (ODEs) or partial differential equation (PDEs), not all solution components may be of interest simultaneously. Clearly, these restrictions cannot be modeled in terms of global norms, but only in terms of so-called goal functionals that specify certain quantities of interest.

In this work, we denote this quantity of interest by \(J\). Even though we are interested in \(J(u)\), all we can obtain is \(J(u_h)\), where \(u_h\) solves 2 , or \(J(\tilde{u})\), where \(\tilde{u}\) is an approximation of \(u_h\). Therefore, we focus on error control of \(J(u)-J(u_h)\), i.e. goal-oriented error estimation for \(J\). There are many techniques for goal-oriented error estimation as, for instance, presented in the works [55], [56], [71], [72], [92], [115], [122], [156], [157], [173], [174].

In order to realize goal-oriented error control in this work, we use the DWR method, which requires solving an adjoint problem. This adjoint formulation follows from the Lagrange formalism and is given by: Find \(z \in V\) such that \[\label{eq:32General32Adjointproblem} \mathcal{A}'(u)(v,z)=J'(u)(v) \qquad \forall v \in U,\tag{3}\] where \(u\) is the solution of 1 and \(\mathcal{A}'\) and \(J'\) describe the Gâteaux-derivatives of \(\mathcal{A}\) and \(J\), respectively. The discrete adjoint problem is given by: Find \(z_h \in V_h\) such that \[\label{eq:32discrete32adjoint} \mathcal{A}'(u_h)(v_h,z_h)=J'(u_h)(v_h) \qquad \forall v_h \in U_h,\tag{4}\] where \(u_h\) is the solution of 2 . The reasoning and details on how the adjoint problem arises will become clear in the following.

The adjoint problem 3 allows us to represent the error in a specific way as shown in the following theorem.

**Theorem 1** (see [52], [87]). *Let \(\tilde{u}
\in U\) and \(\tilde{z} \in V\) be arbitrary but fixed, and let \(u \in U\) be the solution of the model problem 1 , and \(z \in V\) be the solution of the adjoint problem 3 . If \(\mathcal{A} \in \mathcal{C}^3(U,V^*)\) and \(J \in
\mathcal{C}^3(U,\mathbb{R})\), then \[\begin{align} \label{eq:Error32Representation} \begin{aligned} J(u)-J(\tilde{u})&=
\frac{1}{2}\left(\rho(\tilde{u})(z-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}^{(3)} \end{aligned}
\end{align}\qquad{(1)}\] for arbitrary but fixed \(\tilde{u} \in U\) and \(\tilde{z} \in V\), where \[\begin{align}
\label{eq:Error32Estimator:32primal} \rho(\tilde{u})(\cdot) &:= -\mathcal{A}(\tilde{u})(\cdot), \\ \label{eq:Error32Estimator:32adjoint} \rho^*(\tilde{u},\tilde{z})(\cdot) &:= J'(u)-\mathcal{A}'(\tilde{u})(\cdot,\tilde{z}),
\end{align}\] {#eq: sublabel=eq:eq:Error32Estimator:32primal,eq:eq:Error32Estimator:32adjoint} and the remainder term \[\begin{align}
\label{eq:32Error32Estimator:32Remainderterm} \mathcal{R}^{(3)}:=\frac{1}{2}\int_{0}^{1}[J'''(\tilde{u}+se)(e,e,e) -\mathcal{A}'''(\tilde{u}+se)(e,e,e,\tilde{z}+se^*) -3\mathcal{A}''(\tilde{u}+se)(e,e,e)]s(s-1)\,ds,
\end{align}\qquad{(2)}\] with \(e=u-\tilde{u}\) and \(e^* =z-\tilde{z}\).*

**Proof.* The proof can be found in [52], [87]. However, for completeness of our presentation, we will also include the proof here. Let us define \(x := (u,z) \in X:=U \times V\) and \(\tilde{x}:=(\tilde{u},\tilde{z}) \in X\). Moreover, we define the Lagrange function as \[\mathcal{L}(\hat{x}):= J(\hat{u})-\mathcal{A}(\hat{u})(\hat{z}) \quad \forall (\hat{u},\hat{z})=:\hat{x} \in X.\] Since \(\mathcal{A} \in \mathcal{C}^3(U,V^*)\) and \(J \in \mathcal{C}^3(U,\mathbb{R})\) this Lagrange function is in \(\mathcal{C}^3(X,\mathbb{R})\). Using the fundamental theorem of calculus, we can write the difference \(\mathcal{L}(x)-\mathcal{L}(\tilde{x})\) as \[\mathcal{L}(x)-\mathcal{L}(\tilde{x})=\int_{0}^{1} \mathcal{L}'(\tilde{x}+s(x-\tilde{x}))(x-\tilde{x})\,ds.\] Using the trapezoidal rule \[\int_{0}^{1}f(s)\,ds =\frac{1}{2}(f(0)+f(1))+ \frac{1}{2} \int_{0}^{1}f''(s)s(s-1)\,ds,\] with \(f(s):= \mathcal{L}'(\tilde{x}+s(x-\tilde{x}))(x-\tilde{x})\), we end up with \[\begin{align} \mathcal{L}(x)-\mathcal{L}(\tilde{x}) =& \frac{1}{2}(\mathcal{L}'(x)(x-\tilde{x}) +\mathcal{L}'(\tilde{x})(x-\tilde{x})) + \underbrace{\frac{1}{2} \int_{0}^{1}\frac{d^3\mathcal{L}}{ds^3}(\tilde{x}+s(x-\tilde{x}))(x-\tilde{x}) s(s-1)\,ds}_{=\mathcal{R}^{(3)}}. \end{align}\] By using the definition of \(\mathcal{L}\), we obtain that \[J(u)-J(\tilde{u})= \mathcal{L}(x)-\mathcal{L}(\tilde{x}) +\underbrace{\mathcal{A}(u)(z) }_{=0} + \mathcal{A}(\tilde{u})(\tilde{z}) = \frac{1}{2}(\mathcal{L}'(x)(x-\tilde{x}) +\mathcal{L}'(\tilde{x})(x-\tilde{x})) +\mathcal{A}(\tilde{u})(\tilde{z})+ \mathcal{R}^{(3)}.\] We note that \(\mathcal{L}'(x)(y)=0\) for all \(v \in X\). Therefore, the equation from above can be reduced to \[\begin{align} J(u)-J(\tilde{u})=\mathcal{L}(x)-\mathcal{L}(\tilde{x}) +\underbrace{\mathcal{A}(u)(z) }_{=0} + \mathcal{A}(\tilde{u})(\tilde{z}) =& \frac{1}{2}\mathcal{L}'(\tilde{x})(x-\tilde{x}) +\mathcal{A}(\tilde{u})(\tilde{z})+ \mathcal{R}^{(3)}. \end{align}\] Finally, with \[\begin{align} \mathcal{L}'(\tilde{x})(x-\tilde{x}) = &\underbrace{J'(\tilde{u})(e)-\mathcal{A}'(\tilde{u})(e,\tilde{z})}_{=\rho^*(\tilde{u},\tilde{z})(u-\tilde{u})}\underbrace{-A(\tilde{u})(e^*)}_{=\rho(\tilde{u})(z-\tilde{z})}, \end{align}\] we conclude the statement of the theorem. ◻*

**Remark 1**. *A particular example for \(\tilde{u}\) and \(\tilde{z}\) are (finite element) approximations of \(u_h\) and \(v_h\).*

Theoretically, this error identity already gives us an error estimator of the form \[\begin{align} \label{eq:32theory95Error95Estimator} \begin{aligned} J(u)-J(\tilde{u})&= \eta := \frac{1}{2}\left(\rho(\tilde{u})(z-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}^{(3)}. \end{aligned} \end{align}\tag{5}\] However, this error estimator \(\eta\) is not computable, since neither \(u\) nor \(z\) are known in general.

In this subsection, we derive error estimates in enriched spaces that help us to obtain efficiency and reliability results for the error estimator \(\eta\) using a saturation assumption. First, we introduce such enriched spaces \(U_h^{(2)}\) and \(V_h^{(2)}\) with the properties \(U_h \subseteq U_h^{(2)} \subseteq U\) and \(V_h \subseteq V_h^{(2)} \subseteq V\), respectively. For instance, in case of finite elements, this space can be created by using uniform \(h\)- or uniform \(p\)-refinement.

The enriched spaces can be used to formulate the enriched primal problem, which is given by: Find \(u_h^{(2)} \in U_h^{(2)}\) such that: \[\label{eq:32enriched32primal} \mathcal{A}(u_h^{(2)})(v_h^{(2)})=0 \qquad \forall v_h^{(2)} \in V_h^{(2)}.\tag{6}\] Naturally, this can also be done for the adjoint problem leading to the enriched adjoint problem: Find \(z_h^{(2)} \in V_h^{(2)}\) such that \[\label{eq:32enriched32Adjointproblem} \mathcal{A}'(u_h^{(2)})(v_h^{(2)},z_h^{(2)})=J'(u_h^{(2)})(v_h^{(2)}) \qquad \forall v_h^{(2)} \in U_h^{(2)},\tag{7}\] where \(u_h^{(2)}\) solves the enriched primal problem 6 . If we replace \(u\) and \(z\) in the right hand side of 5 with \(u_h^{(2)}\) and \(z_h^{(2)}\), we obtain the approximation \[\label{eq:32def95eta40241} J(u)-J(\tilde{u})\approx \eta^{(2)}:= \frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}^{(3)(2)},\tag{8}\] where \[\begin{align} \label{eq:32Error32Estimator:32Enriched32Remainderterm} \mathcal{R}^{(3)(2)}:=\frac{1}{2}\int_{0}^{1}&[J'''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)}) -\mathcal{A}'''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)},\tilde{z}+se_h^{*(2)}) \\ &-3\mathcal{A}''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)})]s(s-1)\,ds, \end{align}\tag{9}\] with \(e_h^{(2)}=u_h^{(2)}-\tilde{u}\) and \(e_h^{*(2)} =z_h^{(2)}-\tilde{z}\). Such an enriched approach has already been used, e.g., in [9], [13], [62], [96], [103], [145], [169], [175], [176]. With the next theorem, we obtain some results for the error estimator resulting from the enriched approach. Additionally, the theorem allows us to relax differentiability conditions on the solution variables. Specifically, Fréchet-differentiability is only required on \(U_h^{(2)}\) instead of \(U\). For instance, this is important for the parabolic regularized \(p\)-Laplacian in Section 6.4.

**Theorem 2** (see [63]). *Let \(\mathcal{A}: U \mapsto V^*\) and \(J:U\mapsto \mathbb{R}\). Moreover, let \(\mathcal{A}^{(2)} \in \mathcal{C}^3(U_h^{(2)},V_h^{(2)*})\) and \(J_h \in \mathcal{C}^3(U_h^{(2)},\mathbb{R})\) such
that for all \(v_h^{(2)},\psi_h^{(2)} \in U_h^{(2)}\) and \(\phi_h^{(2)} \in V_h^{(2)}\), the equalities \[\begin{align} \label{eq:32A61Ah}
\mathcal{A}(v_h^{(2)})(\phi_h^{(2)})&=\mathcal{A}^{(2)}(v_h^{(2)})(\phi_h^{(2)}), \\ \label{eq:32A3961Ah39} \mathcal{A}'(v_h^{(2)})(\psi_h^{(2)},\phi_h^{(2)})&=(\mathcal{A}^{(2)})'(v_h^{(2)})(\psi_h^{(2)},\phi_h^{(2)}),\\ \label{eq:32J61Jh}
J(\psi_h^{(2)})&=J_h(\psi_h^{(2)}),\\ \label{eq:32J3961Jh39} J'(\psi_h^{(2)})&=J'_h(\psi_h^{(2)}),
\end{align}\] {#eq: sublabel=eq:eq:32A61Ah,eq:eq:32A3961Ah39,eq:eq:32J61Jh,eq:eq:32J3961Jh39} are fulfilled. Here, \(V_h^{(2)*}\) denotes the dual space of \(V_h^{(2)}\). Furthermore,
let us assume that \(J(u) \in \mathbb{R}\), where \(u\in U\) solves the model problem 1 . If \(J(u) \in
\mathbb{R}\), where \(u\in U\) solves the model problem 1 , \(u_h^{(2)} \in U_h^{(2)}\) solves the enriched primal problem 6 and \(z_h^{(2)} \in V_h^{(2)}\) solves the enriched adjoint problem 7 , then for arbitrary but fixed \(\tilde{u} \in U_h^{(2)}\) and \(\tilde{z} \in V_h^{(2)}\) the error representation formula \[\begin{align}
\label{eq:Error32RepresentationEnriched} \begin{aligned} J(u)-J(\tilde{u})&= J(u)-J(u_h^{(2)})+ \frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) +
\mathcal{R}_h^{(3)}, \end{aligned}
\end{align}\qquad{(3)}\] holds, where \(\rho(\tilde{u})(\cdot) := -\mathcal{A}(\tilde{u})(\cdot)\) and \(\rho^*(\tilde{u},\tilde{z})(\cdot) :=
J'(u)-\mathcal{A}'(\tilde{u})(\cdot,\tilde{z})\).*

**Proof.* A similar proof for this theorem is given in [63]. Since \(u_h^{(2)}\) solves the enriched primal problem 6 , and ?? holds, we get \[\mathcal{A}(u_h^{(2)})(v_h^{(2)})=\mathcal{A}^{(2)}(u_h^{(2)})(v_h^{(2)})=0 \quad \forall v_h^{(2)} \in V_h^{(2)}.\] For \(z_h^{(2)}\) solving the enriched adjoint problem 7 , we conclude in combination with ?? and ?? that \[J'(u_h^{(2)})(v_h^{(2)})-\mathcal{A}'(u_h^{(2)})(v_h^{(2)},z_h^{(2)})= J'_h(u_h^{(2)})(v_h^{(2)})-(\mathcal{A}^{(2)})'(u_h^{(2)})(v_h^{(2)},z_h^{(2)})=0 \quad \forall v_h^{(2)} \in U_h^{(2)}.\] This allows us to apply Theorem 1 with \(u=u_h^{(2)}\), \(z=z_h^{(2)}\), \(\mathcal{A}=\mathcal{A}^{(2)}\) and \(J=J_h\). Therefore, we obtain \[\label{eq:32pure32discrete32result} J_h(u_h^{(2)})-J_h(\tilde{u})= \frac{1}{2}\left(\rho_h(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*_h(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u})\right) -\rho_h (\tilde{u})(\tilde{z}) + \mathcal{R}_h^{(3)},\tag{10}\] where \(\rho_h(\tilde{u})(\cdot) := -\mathcal{A}^{(2)}(\tilde{u})(\cdot)\), \(\rho^*_h(\tilde{u},\tilde{z})(\cdot) := J'_h(u)-(\mathcal{A}^{(2)})'(\tilde{u})(\cdot,\tilde{z})\), and \[\begin{align} \mathcal{R}_h^{(3)}:=\frac{1}{2}\int_{0}^{1}&[J_h'''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)}) -(\mathcal{A}^{(2)})'''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)},\tilde{z}+se_h^{*(2)}) \\ &-3(\mathcal{A}^{(2)})''(\tilde{u}+se_h^{(2)})(e_h^{(2)},e_h^{(2)},e_h^{(2)})]s(s-1)\,ds, \end{align}\] with \(e_h^{(2)}=u_h^{(2)}-\tilde{u}\) and \(e_h^{*(2)} =z_h^{(2)}-\tilde{z}\).*

*Again using ?? , ?? and ?? , we notice that \(\rho_h=\rho\) and \(\rho_h^*=\rho^*\) on the enriched spaces \(U_h^{(2)}\) and \(V_h^{(2)}\). This leads to the representation \[\begin{align} J_h(u_h^{(2)})-J_h(\tilde{u})&=
\frac{1}{2}\left(\rho_h(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*_h(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u})\right) -\rho_h (\tilde{u})(\tilde{z}) + \mathcal{R}_h^{(3)}\\ &=
\frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u})\right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}_h^{(3)}.
\end{align}\] In combination with 10 and ?? , we can show \[\begin{align} J(u)-J(\tilde{u})&=J(u)- J(u_h^{(2)})+ J_h(u_h^{(2)})-J_h(\tilde{u})\\ &=J(u)-
J(u_h^{(2)})+\frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u})\right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}_h^{(3)},
\end{align}\] which completes the proof of the theorem. ◻*

**Remark 2**. *With Theorem 2, we know that it is sufficient to fulfill the differentiability conditions \(\mathcal{A}^{(2)} \in \mathcal{C}^3(U_h^{(2)},V_h^{(2)*})\) and \(J_h \in \mathcal{C}^3(U_h^{(2)},\mathbb{R})\) instead of \(\mathcal{A} \in
\mathcal{C}^3(U,V^*)\) and \(J \in \mathcal{C}^3(U,\mathbb{R})\). For instance, point evaluations are well defined on \(U_h^{(2)}=Q_2\), but not in general on \(U\). Additionally, the existence of a solution of the adjoint problem 1 is not required anymore. Instead the existence of the solution \(z_h^{(2)}\in
V_h^{(2)}\) of the enriched adjoint problem 7 is mandatory.*

**Assumption 1** (Saturation Assumption). *Let \(u_h^{(2)} \in U_h^{(2)}\) be the solution of the enriched primal problem 6 . Then there exist \(b_0\in (0,1)\) and \(b_h\in (0,b_0)\) such that \[|J(u)-J(u_h^{(2)})| \leq b_h \, |J(u)-J(\tilde{u})|,\] holds.*

Unfortunately, we are not aware of a general technique to verify the saturation assumption for goal functionals. However, it is a very common assumption in hierarchical based error estimation; see, e.g., [7], [177]–[179]. In the works [73], [178], it was shown that the saturation assumption can fail. However, for specific functionals and PDEs, there are proofs for the saturation assumption; see, e.g., [176], [180]–[187].

**Definition 1** (Efficient and Reliable). *We say an error estimator \(\eta\) is efficient with respect to \(J\), if there exist a constant \(\underline{c} \in \mathbb{R}\) with \(\underline{c}>0\) such that \[\label{eq:32efficient}
\underline{c}|\eta| \leq |J(u)-J(\tilde{u})|.\qquad{(4)}\] We say an error estimator \(\eta\) is reliable with respect to \(J\), if there exist a constant
\(\overline{c} \in \mathbb{R}\) with \(\overline{c}>0\) such that \[\label{eq:32reliable}
\overline{c}|\eta| \geq |J(u)-J(\tilde{u})|.\qquad{(5)}\] *

**Theorem 3**. *Let Assumption 1 be satisfied. Additionally let all assumptions of Theorem 2 be fulfilled. Then the error estimator \(\eta^{(2)}\) defined in 8 , i.e. \[\eta^{(2)}:=\frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}^{(3)(2)},\] is efficient and reliable with the constants
\(\underline{c}=1/(1+b_h)\) and \(\overline{c}=1/(1-b_h)\).*

*Proof.* From Theorem 2, we get that \[\begin{align} \begin{aligned} J(u)-J(\tilde{u})&=
J(u)-J(u_h^{(2)})+ \frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u}) \right) -\rho (\tilde{u})(\tilde{z}) + \mathcal{R}_h^{(3)}, \end{aligned}
\end{align}\] which is equivalent to \[\begin{align} \begin{aligned} J(u)-J(\tilde{u})&= J(u)-J(u_h^{(2)})+ \eta^{(2)}.
\label{help:32show32efficiency32and32reliability} \end{aligned}
\end{align}\tag{11}\] Equation 11 implies that \[\begin{align} \begin{aligned}
|J(u)-J(\tilde{u})|&= |J(u)-J(u_h^{(2)})+ \eta^{(2)}|. \label{help:32show32efficiency32and32reliability32abs} \end{aligned}
\end{align}\tag{12}\] Let us first proof that \(\eta^{(2)}\) is reliable. Since the saturation assumption is valid, we know that \(|J(u)-J(u_h^{(2)})| \leq b_h
|J(u)-J(\tilde{u})|\) holds. Combining this with 12 we obtain \[\begin{align} |J(u)-J(\tilde{u})|&= |J(u)-J(u_h^{(2)})+ \eta^{(2)}| \leq
|J(u)-J(u_h^{(2)})|+ |\eta^{(2)}| \leq b_h |J(u)-J(\tilde{u})|+|\eta^{(2)}|.
\end{align}\] We finally get \[\begin{align} (1-b_h)|J(u)-J(\tilde{u})|\leq |\eta^{(2)}|,
\end{align}\] and consequently \[\begin{align} |J(u)-J(\tilde{u})|\leq \frac{1}{1-b_h}|\eta^{(2)}|.
\end{align}\] Thus, reliability of \(\eta^{(2)}\) follows with the constant \(\overline{c}=1/(1-b_h)\). Now let us prove that \(\eta^{(2)}\) is
efficient as well. We start again with 12 : \[\begin{align} |J(u)-J(\tilde{u})| = |J(u)-J(u_h^{(2)})+ \eta^{(2)}| \geq |\eta^{(2)}|-|J(u)-J(u_h^{(2)})| \geq
|\eta^{(2)}|-b_h |J(u)-J(\tilde{u})|.
\end{align}\] This leads to the estimates \[\begin{align} (1+b_h)|J(u)-J(\tilde{u})|\geq |\eta^{(2)}|,
\end{align}\] and \[\begin{align} |J(u)-J(\tilde{u})|\geq \frac{1}{1+b_h}|\eta^{(2)}|.
\end{align}\] This proves the efficiency of \(\eta^{(2)}\) with the constant \(\underline{c}=1/(1+b_h)\). ◻

If the saturation assumption is fulfilled, the previous theorem shows that the error estimator \(\eta^{(2)}\) defined in 8 is efficient and reliable. Furthermore, this error estimator is also computable and can be employed for measuring discretization and iteration errors. In this subsection, we investigate the different parts of \(\eta^{(2)}\) in more detail. We split \(\eta^{(2)}\) in the following way \[\eta^{(2)}:=\underbrace{\frac{1}{2}\left(\rho(\tilde{u})(z_h^{(2)}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u}) \right)}_{:=\eta_h^{(2)}} \underbrace{-\rho (\tilde{u})(\tilde{z})}_{:=\eta_k} + \underbrace{\mathcal{R}_h^{(3)}}_{:=\eta_{\mathcal{R}}^{(2)}}.\]

This part is of higher order, and usually is neglected in literature. In this work, we will neglect the term as well since it is of higher order regarding the error, and it is connected with high computational cost (depending on the problem). In [73], it was shown that the part indeed is of higher order and can be neglected for \(p\)-Laplace and the Navier-Stokes equations. However, this might not be true in general, as indicated in [66].

The part \[\eta_k:=-\rho (\tilde{u})(\tilde{z}),\] represents the iteration error of the linear or also non-linear solver as presented in [87], [155]. It is the only part of the error estimator which does not depend on the enriched solutions \(u_h^{(2)}\) and \(z_h^{(2)}\). Furthermore, \(\eta_k\) vanishes if \(\tilde{u}\) is the solution of the discrete primal problem 2 . In many algorithms, it is used to stop the non-linear solver. A downside is that, for Newton’s method, the adjoint solution \(\tilde{z}\) is required in every Newton step.

**Theorem 4**. *(Representation of Iteration error; see [73]) Let us assume that \(\tilde{z}\) solves
\[\label{help:32adjoint32disrcrete} \mathcal{A}'(\tilde{u})(\tilde{v},\tilde{z})=J'(\tilde{u})(\tilde{v}) \qquad \forall \tilde{v} \in U_h,\qquad{(6)}\] and let \(\delta \tilde{u}\) solve \[\label{help:32newton32disrcrete} \mathcal{A}'(\tilde{u})(\delta\tilde{u},\tilde{v})=-\mathcal{A}(\tilde{u})(\tilde{v}) \qquad
\forall \tilde{v} \in V_h.\qquad{(7)}\] Then we have the representation \[\rho(\tilde{u},\tilde{z})=J'(\tilde{u})(\delta\tilde{u})\]*

*Proof.* From ?? , ?? and the definition of \(\rho\) in Theorem 3.3 and Theorem 2 it immediately follows that \[\rho(\tilde{u},\tilde{z})=-\mathcal{A}(\tilde{u})(\tilde{z})=\mathcal{A}'(\tilde{u})(\delta
\tilde{u},\tilde{z})=J'(\tilde{u})(\delta\tilde{u}).\] This already concludes the proof. ◻

This identity allows us to use the next Newton update instead of the adjoint solution \(\tilde{z}\), reducing the number of required linear solves from \(2n\) to \(n+1\), where \(n\) is the number of Newton steps. Naturally, this final Newton update \(\delta \tilde{u}\) can be used to update the solution as well.

The part \[\label{eq:32def32primal32and32adjoint32error32estimator} \eta_h^{(2)}:=\frac{1}{2}\bigl(\underbrace{\rho(\tilde{u})(z_h^{(2)}-\tilde{z})}_{:=\eta_{h,p}^{(2)}}+\underbrace{\rho^*(\tilde{u},\tilde{z})(u_h^{(2)}-\tilde{u})}_{:=\eta_{h,a}^{(2)}} \bigr),\tag{13}\] represents the discretization error of the error estimator as discussed in [73], [87], [188]. We would like to mention that, for linear goal functionals \(J\) and affine linear operators \(\mathcal{A}\), the primal part of the discretization error estimator \(\eta_{h,p}\), and the adjoint part \(\eta_{h,a}\), both defined in 13 , coincide. Later on in Subsection 3.10, this part will be localized and used to adapt the mesh.

**Lemma 1**. *Let \(\eta^{(2)}\) be defined as in 8 and let \(\eta_h^{(2)}\) be defined as in 13 . Then \[|\eta_h^{(2)}| - |\eta_k + \eta_{\mathcal{R}}^{(2)}|\leq|\eta^{(2)}|\leq|\eta_h^{(2)}| + |\eta_k + \eta_{\mathcal{R}}^{(2)}|.\]*

*Proof.* We know that \(|\eta^{(2)}|=|\eta_h^{(2)} + \eta_k + \eta_{\mathcal{R}}^{(2)}|\). Therefore, we arrive at the estimates \[|\eta^{(2)}|\leq|\eta_h^{(2)}| + |\eta_k +
\eta_{\mathcal{R}}^{(2)}|,\] and \[|\eta^{(2)}|\geq|\eta_h^{(2)}| - |\eta_k + \eta_{\mathcal{R}}^{(2)}|,\] which concludes the proof. ◻

**Theorem 5** (Almost Efficiency and Reliability). *The two-side discretization error estimate \[\frac{1}{1+b_h}\left(|\eta_h^{(2)}|-|\eta_k + \eta_{\mathcal{R}}^{(2)}|\right)\leq |J(u)-J(\tilde{u})|
\leq\frac{1}{1+b_h}\left(|\eta_h^{(2)}|+|\eta_k + \eta_{\mathcal{R}}^{(2)}|\right)\] is valid provided the saturation assumption [assump:32Saturation] holds. Furthermore, if there exists a \(\alpha_{\eta_h^{(2)}} \in (0,1)\) with \[|\eta_k + \eta_{\mathcal{R}}^{(2)}| \leq
\alpha_{\eta_h^{(2)}} |\eta_h^{(2)}|,\] then the discretization error estimator is efficient and reliable, i. e. \[\underline{c}_h|\eta_h^{(2)}|\leq |J(u)-J(\tilde{u})| \leq\overline{c}_h|\eta_h^{(2)}|,\] with \(\underline{c}_h:=(1-\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}\) and \(\overline{c}_h:=(1+\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}\).*

*Proof.* From Theorem 3, we know that \[\frac{1}{1+b_h}|\eta^{(2)}|\leq
|J(u)-J(\tilde{u})|\leq \frac{1}{1-b_h}|\eta^{(2)}|,\] and, from Lemma 1, we know that \[|\eta_h^{(2)}| - |\eta_k +
\eta_{\mathcal{R}}^{(2)}|\leq|\eta^{(2)}|\leq|\eta_h^{(2)}| + |\eta_k + \eta_{\mathcal{R}}^{(2)}|.\] Combining these two results leads to \[\begin{align}
\frac{1}{1+b_h}\left(|\eta_h^{(2)}|-|\eta_k + \eta_{\mathcal{R}}^{(2)}|\right)&\leq\frac{1}{1+b_h}|\eta^{(2)}|\nonumber\\ &\leq |J(u)-J(\tilde{u})| \label{eq:32thm:32inequality32in32proof}\\ &\leq
\frac{1}{1-b_h}|\eta^{(2)}|\leq\frac{1}{1+b_h}\left(|\eta_h^{(2)}|+|\eta_k + \eta_{\mathcal{R}}^{(2)}|\right).\nonumber
\end{align}\tag{14}\] This is the first statement of the theorem.

From \[|\eta_k + \eta_{\mathcal{R}}^{(2)}| \leq \alpha_{\eta_h^{(2)}} |\eta_h^{(2)}|,\] we can deduce that \[|\eta_h^{(2)}|-|\eta_k + \eta_{\mathcal{R}}^{(2)}| \geq (1-\alpha_{\eta_h^{(2)}})|\eta_h^{(2)}|,\] and \[|\eta_h^{(2)}|+|\eta_k + \eta_{\mathcal{R}}^{(2)}| \leq (1+\alpha_{\eta_h^{(2)}})|\eta_h^{(2)}|.\] Combining this with 14 we obtain \[(1-\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}|\eta_h^{(2)}|\leq |J(u)-J(\tilde{u})| \leq(1+\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}|\eta_h^{(2)}|,\] which provides us the second statement. ◻

This shows that, for the discretization error estimator, we also get efficiency and reliability up to a higher-order term and a term which can be controlled by the accuracy of the non-linear solver. Additionally, if these two terms can be bounded by the inequality \[|\eta_k + \eta_{\mathcal{R}}^{(2)}| \leq \alpha_{\eta_h^{(2)}} |\eta_h^{(2)}|,\] then \(|\eta_h^{(2)}|\) is an efficient and reliable error estimator as well.

In this subsection, we introduce and investigate effectivity indices. These were introduced in [189] in order to measure how well the error estimator approximates the true error. Ideally, the effectivity index approaches \(1\) asymptotically under mesh refinement. However, due to cancellations with contributions with different signs, a (stronger) quality measure for mesh refinement utilizing the triangle inequality resulting into the indicator index was introduced in [62]. As in both quality measures, the true error enters, i.e., the unknown solution \(u\), either academic examples with known \(u\) are taken, or numerical solutions obtained on highly refined meshes (in case the available computational power allows us to do so). Clearly, for complicated non-linear, coupled problems and multiphysics problems, the goal must be to have a previously tested error estimator, which is reliable and efficient, which is then applied to such complicated situations without the need to again measure effectivity and indicator indices.

**Theorem 6** (Bounds on the effectivity index). *We assume that \(|J(u)-J(\tilde{u})|\not=0\). If the assumptions of Theorem 3 are fulfilled, then, for the effectivity index \(I_{\mathrm{eff},+}\) defined by \[\label{def:32Ieff43}
I_{\mathrm{eff},+}:= \frac{\eta^{(2)}}{J(u)-J(\tilde{u})},\qquad{(8)}\] we have the bounds \[\label{eq:32Ieff4332bounds} 1-b_h\leq |I_{\mathrm{eff},+}| \leq
1+b_h.\qquad{(9)}\] If the assumptions of Theorem 5 are fulfilled, then, for the effectivity index
\(I_{\mathrm{eff}}\) defined by \[\label{def:32Ieff} I_{\mathrm{eff}}:= \frac{\eta_h^{(2)}}{J(u)-J(\tilde{u})},\qquad{(10)}\] we have the bounds \[\label{eq:32Ieff32bounds} \frac{1-b_h}{1+\alpha_{\eta_h^{(2)}}}\leq |I_{\mathrm{eff}}| \leq \frac{1+b_h}{1-\alpha_{\eta_h^{(2)}}}.\qquad{(11)}\] *

*Proof.* Here, we follow the ideas in [73], [188]. Theorem 3 provides the result \[\begin{align} \frac{1}{1-b_h}|\eta^{(2)}|\leq|J(u)-J(\tilde{u})|\leq
\frac{1}{1+b_h}|\eta^{(2)}|.
\end{align}\] Now we can divide the inequality from above by \(|J(u)-J(\tilde{u})|\), which leads to \[\begin{align} \frac{1}{1-b_h}\left|\frac{\eta^{(2)}}{J(u)-J(\tilde{u})}\right|\leq 1
\leq \frac{1}{1+b_h}\left|\frac{\eta^{(2)}}{J(u)-J(\tilde{u})}\right|.
\end{align}\] From this, we can easily see the estimates \[1-b_h\leq |I_{\mathrm{eff},+}| \leq 1+b_h.\] The second statement follows from Theorem 5, i.e. \[(1-\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}|\eta_h^{(2)}|\leq |J(u)-J(\tilde{u})|
\leq(1+\alpha_{\eta_h^{(2)}})\frac{1}{1+b_h}|\eta_h^{(2)}|,\] where, by the same argument as above, we get \[\begin{align} \frac{1+\alpha_{\eta_h^{(2)}}}{1-b_h}\left|\frac{\eta^{(2)}}{J(u)-J(\tilde{u})}\right|\leq 1 \leq
\frac{1-\alpha_{\eta_h^{(2)}}}{1+b_h}\left|\frac{\eta^{(2)}}{J(u)-J(\tilde{u})}\right|.
\end{align}\] This is equivalent to \[\frac{1-b_h}{1+\alpha_{\eta_h^{(2)}}}\leq |I_{\mathrm{eff}}| \leq \frac{1+b_h}{1-\alpha_{\eta_h^{(2)}}}.\] ◻

Additionally to the effectivity indices above, we define the primal effectivity indices \(I_{\mathrm{eff},p}\) and adjoint effectivity indices \(I_{\mathrm{eff},a}\) as \[I_{\mathrm{eff},p}:= \frac{\eta_{h,p}^{(2)}}{J(u)-J(\tilde{u})} \qquad \text{and} \qquad I_{\mathrm{eff},a}:= \frac{\eta_{h,a}^{(2)}}{J(u)-J(\tilde{u})}.\]

In Subsection 3.4, we replaced \(u\) and \(z\) in 5 by solutions on the enriched space. In this section we will investigate replacing \(u\) and \(z\) by some arbitrary interpolations \(I_{h,u}: U_h \mapsto U_h^{(2)}\) and \(I_{h,z}: V_h \mapsto V_h^{(2)}\) of \(\tilde{u}\) and \(\tilde{z}\), respectively. One such interpolation is presented in the work [9], [87]. For instance, for \(Q_1\) or \(P_1\) elements, the nodes coincide with the nodes of \(Q_2\) and \(P_2\) finite elements on a coarser mesh, see Figure [fig:32Interpol32Q1] and Figure [fig:32Interpol32P1] for visualization, respectively.

Now, let \(I_{h,u}\tilde{u} \in U_h^{(2)}\) be an arbitrary interpolation, approximating \({u}_h^{(2)} \in U_h^{(2)}\) which solves 6 . Furthermore let, \(I_{h,z}\tilde{z} \in V_h^{(2)}\) be an interpolation, approximating \({z}_h^{(2)} \in V_h^{(2)}\), which solves 7 .

**Theorem 7**. *Let us assume the assumptions of Theorem 2 and let \(\tilde{u} \in U_h\)
and \(\tilde{z} \in V_h\) be arbitrary but fixed. Then for \(I_{h,u}\tilde{u} \in U_h^{(2)}\) and \(I_{h,z}\tilde{z} \in V_h^{(2)}\) holds \[\begin{align} \label{eq:32discrete32error32representation32interpolation32} \begin{aligned} J(I_{h,u}\tilde{u})-J(\tilde{u})&=
\frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right) -\rho (\tilde{u})(\tilde{z}) \\ &-\rho(I_{h,u}\tilde{u})(\frac{I_{h,z}\tilde{z}+
\tilde{z}}{2})+\frac{1}{2}\rho^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{aligned}
\end{align}\qquad{(12)}\] The term \(\tilde{\mathcal{R}}^{(3)(2)}\) is given by \[\begin{align}
\label{Error32Estimator:32Remainderterm1} \tilde{\mathcal{R}}^{(3)(2)}:=\frac{1}{2}\int_{0}^{1}[&J_h'''(\tilde{u}+se_{I_{h,u}})(e_{I_{h,u}},e_{I_{h,u}},e_{I_{h,u}}) \\
-&(\mathcal{A}^{(2)})'''(\tilde{u}+se_{I_{h,u}})(e_{I_{h,u}},e_{I_{h,u}},e_{I_{h,u}},\tilde{z}+se_{I_{h,z}}^*)\\ -&3(\mathcal{A}^{(2)})''(\tilde{u}+se_{I_{h,u}})(e_{I_{h,u}},e_{I_{h,u}},e_{I_{h,u}})]s(s-1)\,ds,
\end{align}\qquad{(13)}\] with \[e_{I_{h,u}}=I_{h,u}\tilde{u}-\tilde{u} \qquad \text{and} \qquad e_{I_{h,z}}^* =I_{h,z}\tilde{z}-\tilde{z}.\]*

**Proof.* The proof is similar to the proof of [63], [87] and Theorem 2 and Theorem 1. First, we define \(x_I := (I_{h,u}\tilde{u},I_{h,z}\tilde{z}) \in X_h^{(2)}:=U_h^{(2)} \times V_h^{(2)}\) and \(\tilde{x}:=(\tilde{u},\tilde{z}) \in X_h^{(2)}\). Since \(\mathcal{A}^{(2)} \in \mathcal{C}^3(U_h^{(2)},V_h^{(2)})\) and the \(J \in \mathcal{C}^3(U_h^{(2)},\mathbb{R})\), we can define a discrete Lagrange functional \[\mathcal{L}_h(\hat{x}):= J_h(\hat{u})-\mathcal{A}^{(2)}(\hat{u})(\hat{z}) \quad \forall (\hat{u},\hat{z})=:\hat{x} \in X_h^{(2)},\] which belongs to \(\mathcal{C}^3(X_h^{(2)},\mathbb{R})\). Following the same steps as in Theorem 1, we get \[\mathcal{L}_h(x_I)-\mathcal{L}(\tilde{x})=\int_{0}^{1} \mathcal{L}'(\tilde{x}+s(x_I-\tilde{x}))(x_I-\tilde{x})\,ds.\] Using the trapezoidal rule \[\int_{0}^{1}f(s)\,ds =\frac{1}{2}(f(0)+f(1))+ \frac{1}{2} \int_{0}^{1}f''(s)s(s-1)\,ds,\] with \(f(s):= \mathcal{L}_h'(\tilde{x}+s(x_I-\tilde{x}))(x_I-\tilde{x})\), cf. [87], we obtain \[\begin{align} \mathcal{L}_h(x_I)-\mathcal{L}(\tilde{x_I}) =& \frac{1}{2}(\mathcal{L}_h'(x_I)(x_I-\tilde{x}) +\mathcal{L}'(\tilde{x})(x_I-\tilde{x})) + \mathcal{R}^{(3)}. \end{align}\] Furthermore, it follows that \[\begin{align} J_h(I_{h,u}\tilde{u})-J_h(\tilde{u})=&\mathcal{L}_h(x_I)-\mathcal{L}(\tilde{x}) +{\mathcal{A}^{(2)}(I_{h,u}\tilde{u})(I_{h,z}\tilde{z}) }- \mathcal{A}^{(2)}(\tilde{u})(\tilde{z}) \\ =& \frac{1}{2}\left(\mathcal{L}_h'(x_I)(x_I-\tilde{x}) +\mathcal{L}_h'(\tilde{x})(x_I-\tilde{x})\right) \\ +&\underbrace{\mathcal{A}^{(2)}(I_{h,u}\tilde{u})(I_{h,z}\tilde{z})}_{=-\rho_h(I_{h,u}\tilde{u})(I_{h,z}\tilde{z})} \underbrace{-\mathcal{A}^{(2)}(\tilde{u})(\tilde{z})}_{=\rho_h(\tilde{u})(\tilde{z})}+ \tilde{\mathcal{R}}^{(3)(2)}. \end{align}\] Investigating the part \(\left(\mathcal{L}_h'(x_I)(x_I-\tilde{x}) +\mathcal{L}_h'(\tilde{x})(x_I-\tilde{x})\right)\), we observe \[\begin{align} \mathcal{L}_h'(x_I)(x_I-\tilde{x}) +\mathcal{L}_h'(\tilde{x})(x_I-\tilde{x}) = & \underbrace{J_h'(I_{h,u}\tilde{u})(I_{h,u}\tilde{u}-\tilde{u})-(\mathcal{A}^{(2)})'(I_{h,u}\tilde{u})(I_{h,u}\tilde{u}-\tilde{u},I_{h,z}\tilde{z})}_{=\rho_h^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})} \\ &\underbrace{-{\mathcal{A}^{(2)}(I_{h,u}\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})}}_{=\rho_h(I_{h,u}\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})}\\ &+\underbrace{J_h'(\tilde{u})(I_{h,u}\tilde{u}-\tilde{u})-(\mathcal{A}^{(2)})'(\tilde{u})(I_{h,u}\tilde{u}-\tilde{u},\tilde{z})}_{=\rho_h^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})} \\ &\underbrace{-\mathcal{A}^{(2)}(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})}_{=\rho_h(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})}. \end{align}\] To sum up, we get \[\begin{align} J_h(I_{h,u}\tilde{u})-J_h(\tilde{u}) &=\frac{1}{2}\left(\rho_h(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho_h^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right)\\ &+\frac{1}{2}\left(\rho_h(I_{h,u}\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho_h^*((I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right)\\ &-\rho_h(I_{h,u}\tilde{u})(I_{h,z}\tilde{z})+\rho_h(\tilde{u})(\tilde{z})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{align}\] After simplifications, we get \[\begin{align} J_h(I_{h,u}\tilde{u})-J_h(\tilde{u}) &=\frac{1}{2}\left(\rho_h(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho_h^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right)\\ &-\frac{1}{2}\rho_h(I_{h,u}\tilde{u})(I_{h,z}\tilde{z}+\tilde{z})+\frac{1}{2}\rho_h^*((I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\\ &+\rho_h(\tilde{u})(\tilde{z})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{align}\] On the discrete spaces \(U_h^{(2)}\) and \(V_h^{(2)}\), we have \(\rho_h=\rho\) and \(\rho_h^*=\rho^*\), hence we get \[\begin{align} \begin{aligned} J(I_{h,u}\tilde{u})-J(\tilde{u})&= \frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right) -\rho (\tilde{u})(\tilde{z}) \\ &-\rho(I_{h,u}\tilde{u})(\frac{I_{h,z}\tilde{z}+ \tilde{z}}{2})+\frac{1}{2}\rho^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{aligned} \end{align}\] This concludes the proof. ◻*

With this identity, we can define the error estimator for interpolation \(\eta_{I}\) as \[\label{eq:32def32etaI} \begin{align} \eta_{I}&:= \frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right) -\rho (\tilde{u})(\tilde{z}) \\ &-\rho(I_{h,u}\tilde{u})(\frac{I_{h,z}\tilde{z}+ \tilde{z}}{2})+\frac{1}{2}\rho^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{align}\tag{15}\] For further results, we require again a saturation assumption for the interpolation.

**Assumption 2** (Saturation assumption for interpolations). *Let \(I_{h,u}\) be an interpolation. Then the saturation assumption is fulfilled if there exists a \(b_0 \in
(0,1)\) and a \(b_h^I\in (0,b_0)\) with \[\label{eq:32saturation32interpol46} |J(u)-J(I_{h,u}\tilde{u})| \leq b_h^I
|J(\tilde{u})-J(u)|.\qquad{(14)}\] *

**Theorem 8**. *Let Assumption 2 be satisfied. Additionally let all assumptions of Theorem 7 be fulfilled. Then the error estimator \(\eta_{I}\) defined
in 15 , i.e. \[\begin{align} \eta_{I}&:= \frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right) -\rho
(\tilde{u})(\tilde{z}) \\ &-\rho(I_{h,u}\tilde{u})(\frac{I_{h,z}\tilde{z}+ \tilde{z}}{2})+\frac{1}{2}\rho^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})+ \tilde{\mathcal{R}}^{(3)(2)}. \end{align}\] is efficient and reliable
with the constants \(\underline{c}=\frac{1}{1+b_h^I}\) and \(\overline{c}=\frac{1}{1-b_h^I}\).*

**Proof.* The proof follows the same steps as in Theorem 3. ◻*

If Assumption 2 is fulfilled, Theorem 8 shows us that the error estimator \(\eta_I\) defined in 8 is efficient and reliable as in the case of enriched spaces. We split \(\eta_I\) into the following parts \[\label{def:32etaI32Parts} \begin{align} \eta_{I}&:= \underbrace{\frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right) }_{:=\eta_{I,h}} \underbrace{-\rho (\tilde{u})(\tilde{z})}_{:=\eta_k} \\ &\underbrace{-\rho(I_{h,u}\tilde{u})(\frac{I_{h,z}\tilde{z}+ \tilde{z}}{2})}_{:=\eta_{I_u}}+\underbrace{\frac{1}{2}\rho^*(I_{h,u}\tilde{u},I_{h,z}\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})}_{:=\eta_{I_z}}+ \underbrace{\tilde{\mathcal{R}}^{(3)(2)}}_{:=\eta_{\mathcal{R}_I}}. \end{align}\tag{16}\] Additionally, we mention that if \(I_{h,u}\tilde{u}=u_h^{(2)}\) solving 6 and if \(I_{h,z}\tilde{z}=z_h^{(2)}\) solving 7 , then \(\eta_I=\eta^{(2)}\).

As for the enriched chase, \(\eta_{\mathcal{R}_I}\) is of higher order, and usually is neglected in literature. However, we would like to mention, that here, higher-order means with respect to \(I_{h,z}\tilde{z}-\tilde{z}\) and \(I_{h,u}\tilde{u}-\tilde{u}\), respectively.

The iteration error estimator for interpolation and enriched solutions is identical. Therefore we again use the same symbol. As already described in Section 3.5.2, this error estimator can be used to stop the linear or non-linear solver.

This part resembles the error which is introduced by the interpolation of the primal variable \(u\) in the enriched space. If \(I_{h,u}\tilde{u}=u_h^{(2)}\), where \(u_h^{(2)}\) solves 6 , then \(\eta_{I_u}=0\). This part can be used to decide whether \(I_{h,u}\tilde{u}\) or \(u_h^{(2)}\) should be used to construct the error estimator.

The adjoint interpolation error estimator \(\eta_{I_z}\) represents the error variable \(z\) in the enriched. Again, if \(I_{h,z}\tilde{z}=z_h^{(2)}\), where \(z_h^{(2)}\) solves 7 , then \(\eta_{I_z}=0\). As for the primal interpolation error estimator, we can decide whether \(I_{h,z}\tilde{z}\) or \(z_h^{(2)}\) is used during the computation.

The part \[\label{eq:32def32primal32and32adjoint32error32estimator32Interpolation} \eta_{I,h}:=\frac{1}{2}\left(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})+\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\right),\tag{17}\] represents the discretization error of the error estimator as discussed in [74], [87]. Here, the part \(\rho(\tilde{u})(I_{h,z}\tilde{z}-\tilde{z})\) and \(\rho^*(\tilde{u},\tilde{z})(I_{h,u}\tilde{u}-\tilde{u})\) do not necessarily coincide, even if the given problem is affine linear. As in the enriched case, this part of the error estimator can be localized, which is shown in Section 3.10. Furthermore, this error estimator can be used to adapt the finite dimensional subspace \(U_h\) or, in case of finite elements, the mesh \(\mathcal{T}_h\).

**Theorem 9** (Efficiency and Reliability Result for \(\eta_{I,h}\)). *For the discretization error estimator it holds \[\frac{1}{1+b_h^I}\left(|\eta_{I,h}|-|\eta_k
+\eta_{\mathcal{R}_I}+\eta_{I_u}+\eta_{I_z}|\right)\leq |J(u)-J(\tilde{u})| \leq\frac{1}{1+b_h^I}\left(|\eta_{I,h}|+|\eta_k +\eta_{\mathcal{R}_I}+\eta_{I_u}+\eta_{I_z}|\right).\] Furthermore, if there exists a \(\alpha_{\eta_{I,h}} \in (0,1)\) such that \[|\eta_k +\eta_{\mathcal{R}_I}+\eta_{I_u}+\eta_{I_z}| \leq \alpha_{\eta_{I,h}} |\eta_{I,h}|,\] then the discretization error estimator is efficient
and reliable, i.e. \[\underline{c}_h|\eta_{I,h}|\leq |J(u)-J(\tilde{u})| \leq\overline{c}_h|\eta_{I,h}|,\] with \(\underline{c}_h:=(1-\alpha_{\eta_{I,h}})\frac{1}{1+b_h^I}\) and \(\overline{c}_h:=(1+\alpha_{\eta_{I,h}})\frac{1}{1+b_h^I}\).*

**Proof.* The proof follows from the same arguments as the proof of Theorem 5. ◻*

In this subsection, we define and investigate effectivity indices for interpolations. For background and motivation of effectivity indices, we refer the reader to Section 3.6.

**Theorem 10** (Bounds on the effectivity index). *We assume that \(|J(u)-J(\tilde{u})|\not=0\). Additionally, let the assumptions of Theorem 3 be fulfilled. Then, for the effectivity index \(I_{\mathrm{eff},+}^I\) defined by \[\label{def:32Ieff4332Int} I_{\mathrm{eff},+}^I:= \frac{\eta_I}{J(u)-J(\tilde{u})},\qquad{(15)}\] we have the bounds \[\label{eq:32Ieff4332bounds32Int}
1-b_h^I\leq |I_{\mathrm{eff},+}^I| \leq 1+b_h^I.\qquad{(16)}\] If the assumptions of Theorem 5 are
fulfilled then for the effectivity index \(I^I_{\mathrm{eff}}\) defined by \[\label{def:32Ieff32Int} I^I_{\mathrm{eff}}:=
\frac{\eta_{I,h}}{J(u)-J(\tilde{u})},\qquad{(17)}\] we have the bounds \[\label{eq:32Ieff32bounds32Int} \frac{1-b_h^I}{1+\alpha_{\eta_{I,h}}}\leq |I^I_{\mathrm{eff}}| \leq
\frac{1+b_h^I}{1-\alpha_{\eta_{I,h}}}.\qquad{(18)}\] *

**Proof.* The proof follows the same steps as the proof of Theorem 6. ◻*

In this subsection, we address error localization such that globally defined error estimators \(\eta\) can be applied locally to mesh elements in order to steer adaptive algorithms. To this end, the previous a posteriori error estimators \(\eta\) need to be split into element-wise or DoF-wise contributions \(\eta_i, i=1,\ldots,N\), where \(N\) is the number of degrees of freedom. Three known approaches are the classical integration by parts [46], a variational filtering operator over patches of elements [41] and a variational partition-of-unity localization [62]. For stationary problems, the effectivity of these localizations was established and numerically substantiated in [62]. The quality measure of the localization process is the so-called indicator index that was as well introduced in [62]:

**Definition 2**. *Let \(u\in U\) be the solution of 1 and \(\tilde{u}\in U_h\). For the definition of the discretization
error \(\eta_h^{(2)}\) in 13 , and localizing to all degrees of freedom \(i=1,\ldots,N\) of the governing triangulation,
the indicator index is defined as \[I_{ind} := \frac{\sum_{i=1}^N |\eta_h^{(2)}(i)|}{|J(u)-J(\tilde{u})|}.\] Here, \(\eta_h^{(2)}(i)\) is the localization of \(\eta_h^{(2)}\) to the degree of freedom \(i\) for \(i=1,\ldots,N\).*

For good effectivity and indicator indices, an important point is the influence of neighboring mesh elements; see e.g., [3]. Consequently, in the traditional method, integration by parts ensured gathering information from neighbor faces and edges [9]. For coupled problems and multiphysics applications, this procedure is error-prone and might be computationally expensive as higher-order operators need to be evaluated. Consequently, the objective is to stay in the weak formulation as done with the filtering approach in [41]. A further simplification was then the so-called PU-DWR approach [62], which we shall explain in the following. Using the PU, we touch different mesh elements per PU-node, and consequently, we gather the required information from neighboring elements. In our adaptive procedures in the remainder of this work, we always employ the PU-DWR method, namely for stationary problems with single goals, multiple goals, and space-time situations with both single and multiple goals.

In this section, we derive abstract results and show that the PU-DWR localization enters simply as a modified test function into the semi-linear/bilinear forms of the operator equation and the right hand side. To start, we first introduce the PU space and its fundamental property. Let us assume \(\{\chi_1,\ldots,\chi_M\}\) is a basis of the PU space \(V_{PU}\) such that \[\begin{align} \label{eq95PU95one} \sum_{i=1}^{M} \chi_{i} \equiv1. \end{align}\tag{18}\] holds. Common choices for the PU spaces are low-order finite element spaces such as \(V_{PU} = Q_1\) [62], \(V_{PU} = P_1\), \(V_{PU} = \widetilde{X}_{k,h}^{0,1}\), i.e., a \(cG(1)dG(0)\) space-time discretization [64] or \(V_{PU} = X_{kh}^{1,1}\), i.e. a \(cG(1)cG(1)\) discretization, for example in \(d+1\)-dimensional space-time discretizations [63]. In general, this ensures a coupling between neighboring temporal elements to address the problem shown in [3]. However, for discontinuous Galerkin discretizations, the dominating edge residuals, i.e. jump terms, are explicitly included in the estimator.

**Proposition 1** (Localized error estimator). *Let the previous PU be given. The localized form of the error estimator \(\eta\) 5 is \[J(u)-J(\tilde{u}) = \eta = \sum_{i=1}^M
\frac{1}{2}\rho(\tilde{u})((z-\tilde{z})\chi_i)
+\frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u-\tilde{u})\chi_i)
+\rho(\tilde{u})(\tilde{z}) + \mathcal{R}^{(3)}(\chi_i).\] As before, the part \(\rho(\tilde{u})(\tilde{z})\) determines the deviation of the approximate solution \(\tilde{u}\) in
comparison to the ‘exact’ discrete solution \(u_h\). This can be iteration errors due iterative linear or non-linear solutions. Since they act globally on the entire solution, no localization to \(i\) with the PU function is required here.*

*Proof.* Inserting the PU function \(\chi_i\) into 5 and using property 18 , we immediately establish the result. ◻

**Corollary 1** (Localized error estimator). *Neglecting the remainder term yields a computable form and it holds \[J(u)-J(\tilde{u}) \approx \eta = \sum_{i=1}^M
\frac{1}{2}\rho(\tilde{u})((z-\tilde{z})\chi_i)
+\frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u-\tilde{u})\chi_i)
+\rho(\tilde{u})(\tilde{z}).\]*

**Definition 3**. *The error estimator in Corollary 1 is composed by the following parts: \[\eta = \eta_h +
\eta_k := \sum_{i=1}^M (\eta_p + \eta_a) + \eta_k,\] where \(\eta_h\) denotes the discretization error and \(\eta_k\) the non-linear iteration error. Specifically, we have \[\begin{align}
\eta_p &:= \eta_p(i) := \frac{1}{2}\rho(\tilde{u})((z-\tilde{z})\chi_i), \\
\eta_a&:= \eta_a(i) := \frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u-\tilde{u})\chi_i),\\
\eta_k &:= \rho(\tilde{u})(\tilde{z}).
\end{align}\]*

In the computational realization it is immediately clear from Section 3 that \(z\) and \(\tilde{z}\) as well as \(u\) and \(\tilde{u}\) are approximated through discrete unknowns from spaces such as \(U_h\). Here, it is important that \(z_h \approx z\) and \(\tilde{z}_h \approx \tilde{z}\) come from different discrete spaces since otherwise \(z_h-\tilde{z}_h \equiv 0\) and \(u_h-\tilde{u}_h \equiv 0\). A practical version of the previous localized form reads:

**Proposition 2** (Practical error estimator). *Let \(\tilde{u}\in U_h\) be a low-order approximation to 2 , \(u_h^{(2)}\in
U_h^{(2)}\) the higher-order solution to 6 , and \(\tilde{z} \in U_h\) be a low-order approximation to 4 \(z_h^{(2)}\in V_h^{(2)}\) the higher-order adjoint solutions 7 , respectively. The practical localized PU error estimator reads \[J(u)-J(\tilde{u})
\approx \eta = \sum_{i=1}^M
\frac{1}{2}\rho(\tilde{u})((z_h^{(2)} - \tilde{z})\chi_i)
+\frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u_h^{(2)} - \tilde{u})\chi_i)
+\rho(\tilde{u})(\tilde{z}),\] where we now re-define the previous notation and obtain as error parts \[\eta = \eta_h + \eta_k := \sum_{i=1}^M (\eta_p + \eta_a) + \eta_k\] with \[\begin{align}
\eta_p &:= \eta_p(i) := \frac{1}{2}\rho(\tilde{u})((z_h^{(2)} - \tilde{z})\chi_i),\\
\eta_a &:= \eta_a(i) := \frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u_h^{(2)} - \tilde{u})\chi_i),\\
\eta_k &:= \rho(\tilde{u})(\tilde{z}).
\end{align}\]*

In this part, we use \(\tilde{u}_h^{(2)}=u_h^{(2)}\) in the case of enriched spaces and \(\tilde{u}_h^{(2)}=I_{h,u}\tilde{u}\) in the case of interpolation. For unified notation we use \(\tilde{u}_h^{(2)}\) for the interpolation and enriched solutions. Furthermore, we use \(\tilde{z}_h^{(2)}=z_h^{(2)}\) in the case of enriched spaces and \(\tilde{z}_h^{(2)}=I_{h,z}\tilde{z}\) in the case of interpolation. As before, let \(\chi_i \in S_1\), where \(S_1 \in \{P_1,Q_1\}\). Then, we know that \(\sum_{i=1}^{M}\chi_i=1\), where \(M:=\operatorname{dim}(S_1)\). It holds \[\begin{align} &\frac{1}{2}\left(\eta_p(\tilde{u})(\tilde{z}_h^{(2)}-\tilde{z})+\eta_a(\tilde{u},\tilde{z})(\tilde{u}_h^{(2)}-\tilde{u})\right), \nonumber\\ =&\frac{1}{2}\left(\eta_p(\tilde{u})((\tilde{z}_h^{(2)}-\tilde{z})\sum_{i=1}^{M}\chi_i)+\eta_a(\tilde{u},\tilde{z})((\tilde{u}_h^{(2)}-\tilde{u})\sum_{i=1}^{M}\chi_i)\right) \nonumber\\ =&\sum_{i=1}^{M}\underbrace{\frac{1}{2}\left(\eta_p(\tilde{u})((\tilde{z}_h^{(2)}-\tilde{z})\chi_i)+\eta_a(\tilde{u},\tilde{z})((\tilde{u}_h^{(2)}-\tilde{u})\chi_i)\right)}_{:=\eta_i^{PU}} \label{eq:32eta95iPU}\\ =&\sum_{i=1}^{M}\eta_i^{PU}. \nonumber \end{align}\tag{19}\] These indicators \(\eta_i^{PU}\) represent an error distribution of the PU or the nodal error contribution. This nodal error estimator is equally distributed to all elements sharing that node. If \(S_1=Q_1\), then adaptive refinement will introduce hanging nodes. The nodal error estimator on these hanging nodes is distributed to the corresponding neighboring nodes as visualized in Figure [fig:32NiPU32hanging32nodes].

**Remark 3**. *Of course the partition of unity technique is not restricted to the finite element method and can be applied to other discretization techniques like isogeometric analysis[27] as well.*

In this section, we briefly describe the fundamental algorithms for enriched approximation and interpolation.

Using the result of the previous subsections of Section 3, we can construct the following algorithm if we use enriched spaces.

**Remark 4**. *In Algorithm 4,the Line [Alg:32Single:32Step32uh2] and [Alg:32Single:32Step32uh] can be swapped.*

From the previous subsections of Section 3, we derive the following algorithm featuring interpolations.

**Remark 5**. *In Algorithm 5, we cannot swap the solution steps with their corresponding interpolation steps, e.g.Line [Alg:32Single:32Step32Iuh2] cannot be done before Line [Alg:32Single:32Step32uh32I].*

The basic error identity sketched in Section 3.3 poses very little assumptions on the underlying problem and in particular on \(u,\tilde{u}\in U\) and the adjoints \(z,\tilde{z}\in V\). It can directly be applied to any \(U\)-conforming discretization, i.e. \(u_h\in U_h\subset U\) and \(z_h\in V_h\subset V\) yielding the identity \[\label{ns:id} J(u)-J(u_h) = \frac{1}{2}\big(\rho(u_h)(z-z_h) + \rho^*(u_h,z_h)(u-u_h)\big) - \rho(u_h)(z_h) + \mathcal{R}^{(3)}.\tag{20}\] If we further consider consistent discretizations, i.e. discrete solutions satisfying \[\label{ns:cons} \begin{align} \mathcal{A}(u_h)(\phi_h) &= 0\quad \forall \phi_h\in V_h,\\ \mathcal{A}'(u_h)(\psi_h,z_h) &= J'(u_h)(\psi_h)\quad \forall \psi_h\in U_h, \end{align}\tag{21}\] and assume for now that the discrete solution \(u_h\in U_h\) is indeed a solution and no iteration error remains, 20 simplifies to \[\label{ns:id1} J(u)-J(u_h) = \frac{1}{2}\big(\rho(u_h)(z-z_h) + \rho^*(u_h,z_h)(u-u_h)\big) + \mathcal{R}^{(3)}.\tag{22}\] This gives rise to the classical formulation of the adjoint error identity [9] that, using Galerkin orthogonality, allows to localize the error by replacing the approximation errors \(u-u_h\) and \(z-z_h\) by any interpolation \(I_h u\in U_h\) and \(I_h z\in V_h\) \[\label{ns:id2} J(u)-J(u_h) = \frac{1}{2}\big(\rho(u_h)(z-I_h z) + \rho^*(u_h,z_h)(u-I_h u)\big) + \mathcal{R}^{(3)}.\tag{23}\] Besides conformity and consistency of the discretization, the fundamental assumption is the variational principle defining the solution \(u\in U\) and its adjoint \(z\in V\).

In the following paragraphs, we examine various problems in which one or more of these assumptions are violated. One straightforward example is the realization of Dirichlet conditions. Assume that the proper solution space would be \[u\in U = u^D+H^1_0(\Omega),\] where \(u^D\in H^1(\Omega)\) is the extension of some boundary data \(g\in H^\frac{1}{2}(\partial\Omega)\). The discrete finite element realization would find \[u_h\in U_h = I_h u^D + U_h^0,\] where \(U_h^0\subset H^1_0(\Omega)\) has homogeneous Dirichlet data, but where \(I_h u^D\neq u^D\) such that \(U_h\not\subset U\). An even simpler example is the weak enforcement of boundary conditions by Nitsche’s method [190], where, in general, \(u_h\neq 0\) on the discrete boundary and hence \(u_h\not\in H^1_0(\Omega)\). Further examples leading to non-conforming discretizations are approximations on curved domains, where \(\Omega\neq \Omega_h\) [191], non-conforming finite elements such as the Crouzeix-Raviart element [192] or discontinuous Galerkin methods [193], [194] for elliptic problems.

Another potential source of problems lies in the non-consistency of the discrete formulation. While consistency is required to use Galerkin orthogonality for localization in the spirit of 23 , it is not essential for the path outlined in Section 3 as the error identity contains the term \(\rho(\tilde{u})(\tilde{z})\), which, on the one hand, stands for iteration errors, but which also includes all non-conformity errors. Classical sources of non-consistency are stabilized finite element methods, where the discrete variational formulation must be enriched. Examples are transport stabilizations such as streamline diffusion for a simple diffusion transport problem \[\label{ns:stab} \mathcal{A}_h(\cdot)(\cdot) = \mathcal{A}(\cdot)(\cdot)+\mathcal{S}(\cdot)(\cdot),\tag{24}\] where \[\label{ns:stab:sd} \mathcal{A}_h(u)(\phi) = \underbrace{(\nabla u,\nabla \phi) + (\beta\cdot\nabla u,\phi)}_{={\mathcal{A}}(u)(\phi)} + \underbrace{(\delta_h\beta\cdot\nabla u,\beta\cdot\nabla \phi)}_{={\mathcal{S}}(u)(\phi)}\tag{25}\] and where the solution \(u\in U\) to \(\mathcal{A}(u)(\phi)=0\) for all \(\phi\in V\) does not satisfy the discrete problem \(\mathcal{A}_h(u)(\phi)=0\). Other examples are stabilizations of saddle-point problems such as local projections [195]. Non-consistency is also relevant for time-dependent problems. While the DWR estimator can directly be applied to space-time Galerkin methods as discussed in Section 3, Galerkin time discretizations might not be favorable in terms of computational complexity. Instead, efficient time-stepping methods are used for simulation and their similarity to certain Galerkin space-time discretization is utilized for error estimation only [61], [196].

The idea of different methods in simulation and theoretical error analysis can be taken even further. For example, the error estimator can in principle be applied to discrete solutions that are not based on variational principles at all, but are obtained using finite volume methods [197] or neural networks [165].

As before, let \(\mathcal{A}(\cdot)(\cdot)\) be the variational semi-linear form describing the problem at hand. Then, we have: Find \[u\in U: \quad \mathcal{A}(u)(v) = 0\quad \forall v\in V.\] Then, \(u_h\in U_h\subset U\) shall be the discrete approximation that is given in the modified variational form \[u_h\in U_h\subset U\quad \mathcal{A}_h(u_h)(v_h) = 0\quad \forall v_h\in V_h,\] where we assume that consistency does not hold, i.e. in the general case we have \[\mathcal{A}_h(u)(v_h) \neq 0\quad\forall v_h\in V_h.\] For the following we assume that the discrete variational form can be written as \(\mathcal{A}_h(\cdot)(\cdot)=\mathcal{A}(\cdot)(\cdot)+\mathcal{S}_h(\cdot)(\cdot)\), and that the discrete adjoint problem is given by: \[\text{Find z_h\in V_h:} \quad \mathcal{A}_h'(u_h)(v_h,z_h) = J'(u_h)(v_h)\quad\forall v_h\in U_h,\] with \(\mathcal{A}_h'(\cdot)(\cdot,\cdot) = \mathcal{A}'(\cdot)(\cdot,\cdot) + \mathcal{S}^*_h(\cdot)(\cdot,\cdot)\), where \(\mathcal{S}^*_h\) is not necessarily the derivative of \(\mathcal{S}_h(\cdot,\cdot)\).

**Theorem 11**. *Let \(u\in U\) and \(z\in V\) be primal and adjoint solutions such that \[\mathcal{A}(u)(v) = 0\quad\forall v\in V,\quad
\mathcal{A}'(u)(v,z) = J'(u)(v)\quad\forall v\in U.\] Further, let \(u_h\in U_h\subset U\) and \(z_h\in V_h\subset V\) be primal and adjoint discrete approximations such that
\[\mathcal{A}_h(u_h)(v_h) \approx 0\quad\forall v_h\in V_h,\quad \mathcal{A}_h'(u_h)(v_h,z_h) \approx J'(u_h)(v_h)\quad\forall v_h\in U_h,\] where \[\label{ns:stabform} \mathcal{A}_h(u)(v) = \mathcal{A}(u)(v) + \mathcal{S}_h(u)(v),\quad \mathcal{A}'_h(u)(v,z) = \mathcal{A}'(u)(v,z) + \mathcal{S}_h^*(u)(v,z).\qquad{(19)}\] Then, the following error
representation holds \[\label{eq:ConsistencyError} J(u)-J(u_h)= \frac{1}{2}\left(\rho(u_h)(z-z_h)+\rho^*(u_h,z_h)(u-u_h)\right) - \rho_h(u_h)(z_h)-\mathcal{S}_h(u_h)(z_h) +
\mathcal{R}^{(3)}.\qquad{(20)}\] Here, \(\rho(\cdot)(\cdot),\rho^*(\cdot,\cdot)(\cdot)\) and \(\mathcal{R}^{(3)}\) are defined as in Theorem 1, whereas \[\rho_h(u)(v) = -A_h(u)(v),\quad \rho_h^*(u,z)(v) = J'(u)(v)-A'(u)(v,z).\] Given interpolations \(I_h:U\to
U_h\) and \(I_h^*:V\to V_h\), the localization of this error identity reads \[\begin{gather} \label{eq:ConsistencyErrorLocal} J(u)-J(u_h)=
\frac{1}{2}\left(\rho(u_h)(z-I_h^*z)+\rho^*(u_h,z_h)(u-I_hu) \right)\\ -\rho_h(u_h)(z_h) +\frac{1}{2}\left(\rho_h(u_h)(I_h^*z - z_h)+\rho_h^*(u_h,z_h)(I_h^*u - u_h)\right)\\ -\frac{1}{2}\left(\mathcal{S}_h(u_h)(I_h^*z - z_h)+\mathcal{S}_h^*(u_h)(I_h^*u -
u_h,z_h)\right) - \mathcal{S}_h(u_h)(z_h) + \mathcal{R}^{(3)}.
\end{gather}\qquad{(21)}\] *

**Proof.* The error identity ?? of Theorem 1 does not require consistency and is directly applicable \[J(u)-J(u_h)= \frac{1}{2}\left(\rho(u_h)(z-z_h)+\rho^*(u_h,z_h)(u-u_h)\right) - \rho(u_h)(z_h) + \mathcal{R}^{(3)}.\] Using ?? , ?? directly follows.*

*If we want to exploit Galerkin orthogonality to localize by means of replacing the approximation errors \(u-u_h\) and \(z-z_h\) by interpolation weights, we introduce \(\pm \rho(u_h)(z_h-I_h^*z)\) and \(\pm \rho^*(u_h)(u_h-I_h u,z_h)\), and we obtain \[\begin{gather}
\label{ns:7} J(u)-J(u_h)= \frac{1}{2}\left(\rho(u_h)(z-I_h^*z)+\rho^*(u_h,z_h)(u-I_hu) \right) -\rho_h(u_h)(z_h) - \mathcal{S}_h(u_h)(z_h)\\ +\frac{1}{2}\left(\rho(u_h)(I_h^*z - z_h)+\rho^*(u_h,z_h)(I_h^*u - u_h)\right) + \mathcal{R}^{(3)}.
\end{gather}\qquad{(22)}\] We further use ?? to split \(\rho(\cdot)(\cdot) = \rho_h(\cdot)(\cdot)-\mathcal{S}_h(\cdot)(\cdot)\), and likewise \(\rho^*(\cdot,\cdot)(\cdot)\) to
separate the error identity into weighted residuals, iteration errors and consistency errors to reach ?? . ◻*

The error identities ?? and ?? now consist of five parts: primal and dual weighted residuals, the discrete iteration error, the non-consistency error and, finally, the remainder. Apart from the remainder, all other terms can be evaluated numerically, where for the residuals one of the approximations discussed in Section 3.4 or Section 3.7 has to be used. The iteration errors vanish if the algebraic problems are solved to sufficient precision.

**Remark 6** (Consistency remainders). *The two additional non-consistency terms appearing in ?? can not be directly evaluated, as they depend on the unknown exact solutions \(u\in U\) and \(z\in V\). However, they are usually negligible, as they carry an additional order compared to the consistency term \(S_h(u_h)(z_h)\). As an example we consider the streamline diffusion
stabilization 25 , where \[\mathcal{S}_h(u_h,v_h) = (\delta_h \beta\cdot\nabla u_h,\beta\cdot\nabla v_h)\] and \(\mathcal{S}^*_h(u_h)(v_h,z_h) =
\mathcal{S}_h(v_h,z_h)\). While \(\mathcal{S}_h(u_h)(z_h)\) can simply be estimated as \[\label{ns:x} \big|\mathcal{S}_h(u_h)(z_h)\big| \le c |\delta_h|\,\|\nabla
u_h\|\cdot \|\nabla z_h\|,\qquad{(23)}\] primal and adjoint consistency remainders give rise to \[\big|\mathcal{S}_h(u_h)(I_h^*-z_h)\big| \le c |\delta_h|\,\|\nabla u_h\|\cdot \big(\|\nabla (z-z_h)\|+\|\nabla
(z-I_h^*z)\|\big),\] which is of higher order as compared to the main consistency term ?? .*

Besides stabilization techniques, another typical application of tracking consistency errors in error estimation is found in time-stepping methods. To illustrate this, let \(u'(t) = f(t,u(t))\) be a given scalar initial value problem on \(I=(0,T)\) with \(u(0)=u_0\). By \(0=t_0<t_1<\dots<t_N=T\) we introduce discrete points in time, where we uniformly assume \(\Delta t=t_n-t_{n-1}\) for simplicity of notation only. Simple time-stepping schemes like the backward Euler method \[\label{ns:ie} u^n + \Delta t f(t_n,u^n) = u^{n-1}\tag{26}\] find their counterpart in Galerkin time-stepping methods such as introduced by Eriksson, Estep, Hansbo, and Johnson [8] as well as Thomée [30]. Time-stepping methods and Galerkin methods must still be seen as two distinct methods. Introducing the discontinuous space \[U_h = \{ v\in L^2(I) | \; v(t_0)\in\mathbb{R},\; v\big|_{(t_{n-1},t_n]}\in\mathbb{R},\; n=1,\dots,N\},\] the backward Euler method 26 can be formulated as to find \(u_h\in U_h\) such that \[A_h(u_h)(v_h) = 0\quad\forall v_h\in V_h=U_h\] with the fully discrete variational formulation \[\label{var:be} A_h(u_h)(v_h) =\left(u_h(t_0)-u_0\right)\cdot v_h(t_0) + \sum_{n=1}^N \left(u_h(t_n)-u_h(t_{n-1}) + \Delta t f\left(t_n,u(t_n)\right)\right)\cdot v_h(t_n).\tag{27}\] In contrast, the corresponding discontinuous Galerkin approach also defines \(u_h\in U_h\) given by \[A(u_h)(v_h) = 0\quad\forall v_h\in U_h,\] where the continuous variational form is \[\begin{gather} \label{var:dg} A(u_h)(v_h) =\left(u_h(t_0)-u_0\right)\cdot v_h(t_0)\\ + \sum_{n=1}^N\left\{ \left(u_h(t_n)-u_h(t_{n-1}\right)\cdot v_h(t_n) +\int_{t_{n-1}}^{t_n} f(t,u_h(t))\cdot v_h(t)\right\}. \end{gather}\tag{28}\] Both forms, 27 and 28 , only differ by numerical quadrature. However, as 27 is based on the box rule which has the same order of convergence as the backward Euler method itself, both discrete approaches must be considered substantially different.

The consistency error can be introduced as \(\mathcal{S}_h(u)(z) = A_h(u)(z)-A(u)(z)\) and for the Euler method we obtain \[\label{var:dg95consistency} \begin{align} |S(u_h)(z_h)| &= \sum_{n=1}^N \int_{t_{n-1}}^{t_n} \Delta t\cdot |z_h(t_n)|\cdot \left| f(t_n,u_h(t_n)) - \int_{t_{n-1}}^{t_n} f(t,u_h(t))\,\text{d}t\right|\\ &= \sum_{n=1}^N \int_{t_{n-1}}^{t_n} \Delta t^2\cdot |z_h(t_n)|\cdot \left| f'(\xi_n) + O(\Delta t)\right|, \end{align}\tag{29}\] as \(u_h(t)=u_h(t_n)\) for \(t\in (t_{n-1},t_n]\) and assuming that \(f(\cdot)\) is differentiable. This yields the same order as the truncation error of the backward Euler method itself. The consistency term cannot be neglected, but its evaluation must be included as part of the error identity.

If efficient time-stepping methods are used in simulations and a posteriori error estimation based on the DWR method is applied, the consistency error hence has to be tracked. For details and applications to flow problems discretized by the Crank-Nicolson scheme and Fractional-Step-\(\theta\) method, we refer to [61], [130], [196].

The application to non-conforming discretizations is more cumbersome and there is no one standard approach. This is already due to the fact that even the definition of the error \(u-u_h\) can be difficult. This is the case, for example, if \(u\) is given on a domain \(\Omega\) and \(u_h\) on a discretized domain \(\Omega_h\). In this case, one remedy could be to limit all considerations to the intersection \(\Omega\cap\Omega_h\), see for instance [191].

Since the sources of non-conformity are manifold, we discuss two examples in the following, first classical non-conforming finite element approaches and then discretizations that are not based on variational principles at all.

Let \(u\in U\) be the solution to \[\mathcal{A}(u)(v) = 0\quad\forall v\in V.\] By \(U_h\not\subset U\) and \(V_h\not\subset V\) we denote a non-conforming finite element discretization. As specific problem, consider a \(H^1\)-elliptic problem in mind with \(U=V=H^1_0(\Omega)\) and where \(U_h=V_h\) is the non-conforming Crouzeix-Raviart element. We must assume that the variational form \(\mathcal{A}(\cdot)(\cdot)\) is also defined in \(U_h\times U_h\) and elliptic. Hence, we introduce \(X=U\cup U_h\) and consider \(\mathcal{A}:X\times X\to\mathbb{R}\) and, likewise, we assume that the functional is defined on \(J:X\to\mathbb{R}\), which is a limitation as, for instance, line integrals along mesh edges are not well defined on \(U_h\).

The direct application of Theorem 1 is not possible, since the fact that the residual of the continuous solution disappears in a conformal discretization was exploited here, namely \(\mathcal{L}'(u,z)(u-u_h,z-z_h)=0\). This is not necessarily the case for a non-conforming discretization. Instead, what remains is the term \[\mathcal{L}'(u,z)(u-u_h,z-z_h)= \frac{1}{2}\left( -\mathcal{A}(u)(z-z_h) + J'(u)(u-u_h) - \mathcal{A}'(u)(u-u_h,z)\right).\] Its meaning must be considered on a case-by-case basis. For Poisson’s problem \(\mathcal{A}(u)(v) = (\nabla u,\nabla v)-(f,v)\), the error identity from Theorem 1 reads as \[\begin{gather} J(u-u_h) = \frac{1}{2}\left(\rho(u_h)(u-u_h)+\rho^*(u_h,z_h)(z-z_h)\right) - \rho(u_h)(z_h)\\ \frac{1}{2}\left( (f,z-z_h) - \left(\nabla u,\nabla (z-z_h)\right)+ J'(u)(u-u_h) - \left(\nabla (u-u_h),\nabla z\right)\right) +\mathcal{R}^{(3)}. \end{gather}\] Using integration by parts, these terms get \[\begin{align} (f,z-z_h) - \left(\nabla u,\nabla (z-z_h)\right) &= \sum_{{K\in \mathcal{T}_h}} \langle \partial_n u,[z_h]\rangle_{{\partial K}}\\ J'(u)(u-u_h) - \left(\nabla (u-u_h),\nabla z\right)&= \sum_{{K\in \mathcal{T}_h}} \langle \partial_n z,[u_h]\rangle_{{\partial K}}, \end{align}\] where \([\cdot]\) denotes the jump over the element’s edge \(\partial K\) and \(\mathcal{T}_h\) is the decomposition of \(\Omega_h\). This term exactly measures the non-conformity of the discrete approximations \(u_h\) and \(z_h\). Instead of a rigorous bound, it can be estimated based on replacing \(\partial_n u\) and \(\partial_n z\) by discrete reconstructions, see [198] for a similar procedure. We refer to [114] for a detailed presentation of a goal-oriented error estimator for different conforming and non-conforming discretizations.

Finally, we explore to which extend the error estimator can be applied to discretizations that cannot be written in the form of a variational formulation at all. These could be, for example, finite difference methods, or the approximation of
differential equations with neural networks like *Physics Informed Neural Networks (PINNs)* such as the *Deep Ritz* method [199]. In both
cases, the discrete approximation \(u_h\) cannot be represented by a variational formulation. Once again, we cannot rely on a uniform theoretical principle, but must argue on a case-by-case basis.

What we require in any case is the embedding of the discrete solution \(u_h\) into a function space \(\mathcal{E}:u_h\mapsto U_h\) that is in some sense compatible, meaning that it is either conforming \(U_h\subset U\) or that the variational form can be extended onto \(X:=U\cup U_h\) such as in the case of the non-conforming Crouzeix-Raviart element. Considering finite difference approximations \(\mathcal{E}(u_h) = I_h u_h\), the interpolation into a usual finite element space is one conforming option. In this setting, the original error identity from Theorem 1 can directly be applied:

**Theorem 12** (Error identity for conforming embeddings). *Let \(u\in U\) and \(z\in V\) be the primal and dual solutions. Let \(u_h\)
and \(z_h\) be any discrete approximations. Further, let \[\mathcal{E}:u_h\mapsto U_h\subset U,\quad
\mathcal{E}^*:z_h\mapsto V_h\subset V,\] be embeddings of the discrete approximations into subspaces. Then, it holds \[\begin{gather}
J(u)-J(\mathcal{E}(u_h)) = \frac{1}{2}
\left(
\rho\left(\mathcal{E}(u_h)\right)\left(z-\mathcal{E}^*(z_h)\right)+
\rho^*\left(\mathcal{E}(u_h),\mathcal{E}^*(z_h)\right)\left(u-\mathcal{E}(u_h)\right)\right)\\
-\rho\left(\mathcal{E}(u_h)\right)\left(\mathcal{E}^*(z_h)\right) + \mathcal{R}^{(3)}.
\end{gather}\]*

Like all previous error identities, this one suffers from the need to approximate numerically the primal and adjoint solutions \(u\in U\) and \(z\in V\). Furthermore, the embeddings \(\mathcal{E}(u_h)\) and \(\mathcal{E}^*(z_h)\) must be available for numerical evaluation. But most important, the error estimator is still based on variational principles although the discrete approximations \(u_h\) and \(z_h\) can be obtained in completely different ways. To evaluate the error, integrals over the computational domains containing the embeddings \(\mathcal{E}(u_h)\) and \(\mathcal{E}^*(z_h)\) must be computed. This will either require a mesh, which is all but trivial, if the choice of discretization method - e.g. smoothed particle hydrodynamics - was motivated as being mesh-free, or, if, for instance, neural networks are chosen for high-dimensional problems. Resorting to simple Monte Carlo quadrature avoids this problem, but it brings along a substantial quadrature error. We refer to [165] for a first discussion on this mostly open topic.

In the past two decades, numerous efforts regarding the efficient and robust numerical solution of time-dependent, non-linear, coupled partial differential equations (PDE) have been made. Often, such PDE systems arise in continuum mechanics from conservation laws such as conservation of mass, momentum, angular momentum, and energy. Examples are porous media applications [200]–[202], multi-phase flow and heat transfer processes [203], and fluid-structure interaction [130], [204], [205]. In addition, such multiphysics problems may be subject to inequality constraints, which result into coupled variational inequality systems (CVIS), for which multiphysics phase-field fracture is an example [134]. In such problems, often more than one single goal functional shall be controlled and measured up to a certain accuracy. This is the motivation for multigoal-oriented error control and adaptivity that we will explain in this section.

Let us assume that we are interested in the evaluation of \(N\) functionals, which we denote by \(J_1, J_2, \ldots,J_N\). This can also be seen as a vector valued quantity of interest \[\label{eq:32def32vecJ} \vec{J}(v):=(J_1(v),J_2(v),\cdots, J_{N}(v) ).\tag{30}\] There are a lot of works on multi-goal oriented adaptivity [48]–[52], [56], [56], [101], [103], [132], [206]–[208]. The aim of this section is to combine the functionals \(J_1, J_2, \ldots,J_N\) into a single functional and apply the results from Section 3. A simple idea to combine this functionals is to add them up. We define \(J_A\) as \[\label{eq:32def32Addfunctionals} J_A := \sum_{i=1}^N J_i.\tag{31}\] This functional allows us to combine the functionals. Let us assume that \(u\) solves the primal problem 1 , and let \(\tilde{u} \in U_h^{(2)}\) be approximation to \(u_h\) solving 2 . For instance, let us assume \(N=3\) and the error \(J_1(u)-J_1(\tilde{u})=1\), \(J_2(u)-J_2(\tilde{u})=1\) and \(J_3(u)-J_3(\tilde{u})=-2\). Then \(J_A(u)-J_A(\tilde{u})=0\). This shows that \(J_A\) is not suitable, since we get that the error vanishes even though the error does not vanish in any of the original functionals.

In this part, we introduce a weighting function to balance the sum of the different goal functionals. This function includes a sign evaluation in order to avoid cancellation of goal functionals with similar values, but different signs.

**Definition 4** (Error-Weighting function). *Let \(M \subseteq \mathbb{R}^N\). The function \(\mathfrak{E}: (\mathbb{R}^+_0)^N \times M \mapsto \mathbb{R}^+_0\) is an
*error-weighting function* if \(\mathfrak{E}(\cdot,m) \in \mathcal{C}^3((\mathbb{R}^+_0)^N,\mathbb{R}^+_0)\) is strictly monotonically increasing in each component and \(\mathfrak{E}
(0,m)=0\) for all \(m \in M\).*

Examples for such error weighting functions are \[\label{eq:error32weighting32function:example321} \mathfrak{E}(x,m):=\sum_{i=1}^N \frac{x_i}{|m_i|},\tag{32}\] \[\label{eq:error32weighting32function:example322} \mathfrak{E}(x,m):=\sum_{i=1}^N x_i,\tag{33}\] \[\label{eq:error32weighting32function:example323} \mathfrak{E}(x,m):=\sum_{i=1}^N \frac{x_i^p}{|m_i|^p} \qquad p \in (1,\infty),\tag{34}\] and \[\label{eq:error32weighting32function:example324} \mathfrak{E}(x,m):=\sum_{i=1}^N \sqrt{x_i},\tag{35}\] for \(x,m \in (\mathbb{R}^+_0)^N \times M\). These functions should mimic some kind of norm or metric. Here, the elements \(m\in M\) are used to weight the contributions. These are user chosen weights, or weights balancing the relative errors instead of the absolute ones. Finally, we define the error functional as follows \[\begin{align} \label{eq:32TrueErrorrepresentationFunctional} \tilde{J}_{\mathfrak{E}}(v):=\mathfrak{E}(|\vec{J}(u)-\vec{J}(v)|_N, m) \qquad \forall v \in U, \end{align}\tag{36}\] where, \(|\cdot|_N\) describes the component wise absolute value. It follows from the definition of \(\mathfrak{E}\) that \(J_\mathfrak{E}(v) \in \mathbb{R}^+_0\) for all \(v \in V\).

**Remark 7**. *The error functional \(\tilde{J}_{\mathfrak{E}}(v)\) mimics a semi-metric (as in [209], [210]) for the errors in the functionals \(J_i\) Hence, \(\tilde{J}_{\mathfrak{E}}(v)\)
represents a semi-metric, which ensures that \(\tilde{J}_{\mathfrak{E}}\) is monotonically increasing if \(|J_i(u)-J_i(\tilde{u})|\) is monotonically increasing.*

We notice that \(\tilde{J}_{\mathfrak{E}}(v)\) as defined in 36 is not computable, since we do not know the exact solution \(u\). And if we were to know \(u\), we also would know all \(J_i(u)\), for which we would not need any finite element simulations. In Section 3, we had a similar problem. There, we introduced the enriched solutions \(u_h^{(2)}\) solving 6 , which was used to replace all \(u\) in the error identity of Theorem 2. To this end, we apply the same approach and we define the computable error functional \[\begin{align} \label{eq:32ErrorrepresentationFunctional} {J}_{\mathfrak{E}}(v):=\mathfrak{E}(|\vec{J}(u_h^{(2)})-\vec{J}(v)|_N,m) \qquad \forall v \in U, \end{align}\tag{37}\] where \(u_h^{(2)}\) solves the enriched primal problem 6 . This functional allows us to use the ideas from Section 3 with \(J=J_\mathfrak{E}\). From 37 we conclude that \[\label{eq:32Errorweighting-function32def} {J}_{\mathfrak{E}}'(v)(\delta v)=\sum_{i=1}^{N}\operatorname{sign}\left(J_i(u_h^{(2)})-J_i(v)\right)\frac{\partial \mathfrak{E}}{\partial x_i}(|\vec{J}(u_h^{(2)})-\vec{J}(v)|_N,m) \qquad \forall v,\partial v \in U,\tag{38}\] The value \(\operatorname{sign}\left(J_i(u)-J_i(v)\right)\) which is approximated by the \(\operatorname{sign}\left(J_i(u_h^{(2)})-J_i(v)\right)\) in 38 is important to avoid error cancellations. In [48], [49], the sign was approximated using an adjoint to adjoint (dual to dual) problem, whereas in a prior work [52] the sign was approximated using the difference \(J_i(u_h^{(2)})-J_i(\tilde{u})\). For a linear partial differential equation and linear goal functionals those two methods coincide.

Common choices are \(m=|\vec{J}(\tilde{u})|_N\) or \(m=|\vec{J}(u_h^{(2)})|_N\). Together with the error function 32 , this choice aims for a similar relative error in all the functionals. In contrast the error function 33 aims for a similar absolute error. Furthermore, the error function 34 penalizes larger relative errors, and the error function 35 aims for a similar decrease of the error in all the functionals. The choice \(m_i=\omega_i |J_i(\tilde{u})|\), where \(\omega_i \in\mathbb{R}^+_0\) are user chosen weights, in combination with the error functional 32 , leads to an error functional with almost the same properties as \(J_c\) in [48], [49], [51], [52], [103], [145], [206], [211] defined as \[\label{eq:32def32J95c} J_c(v):=\sum_{i=1}^{N}w_i J_i(v),\tag{39}\] where \(w_i:=\omega_i\operatorname{sign}\left(J_i(u_h^{(2)})-J_i(\tilde{u})\right)/|J_i(\tilde{u})|\).

Since we work with the combined goal functional, the localization procedure can be carried out in the same fashion as previously described in Section 3.10. An immediate consequence of Proposition 2 is

**Proposition 3** (Practical error estimator for the functional \(J_\mathfrak{E}\)). *Let \(\tilde{u}\in U_h\) be a low-order approximation to 2 , \(u_h^{(2)}\in U_h^{(2)}\) the higher-order solution to 6 , and \(\tilde{z} \in U_h\) be a low-order
approximation to 4 \(z_h^{(2)}\in V_h^{(2)}\) the higher-order adjoint solutions 7 , respectively. The practical localized PU error
estimator reads for the error functional reads \[J_\mathfrak{E}(u)-J_\mathfrak{E}(\tilde{u}) \approx \eta^{(2)} = \sum_{i=1}^M
\frac{1}{2}\rho(\tilde{u})((z_h^{(2)} - \tilde{z})\chi_i)
+\frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u_h^{(2)} - \tilde{u})\chi_i)
+\rho(\tilde{u})(\tilde{z})\] where we now re-define the previous notation and obtain as error parts \[\eta^{(2)} = \eta_h^{(2)} + \eta_k := \sum_{i=1}^M (\eta_p + \eta_a) + \eta_k\] with \[\begin{align}
\eta_p &:= \eta_p(i):= \frac{1}{2}\rho(\tilde{u})((z_h^{(2)} - \tilde{z})\chi_i),\\
\eta_a &:= \eta_a(i):= \frac{1}{2}\rho^*(\tilde{u},\tilde{z})((u_h^{(2)} - \tilde{u})\chi_i),\\
\eta_k &:= \rho(\tilde{u})(\tilde{z}).
\end{align}\]*

Summarizing the previous ingredients allows us to formulate multigoal algorithms for adaptivity in which non-linear iteration errors are balanced with discretization errors.

**Remark 8**. *Since we have to approximate \(\operatorname{sign}\left(J_i(u_h^{(2)})-J_i(v)\right)\), we only provide an algorithm using enriched spaces.*

**Remark 9**. *If the non-linear solver is stopped using the iteration error estimator \(\eta_k\), we require \(u_h^{(2)}\) before \(\tilde{u}\). Therefore, in contrast to Algorithm 4, in Algorithm 6 Line [Alg:32Multi:32Step32uh2] and [Alg:32Multi:32Step32uh] must not be swapped.*

In this section, we substantiate our theoretical results and numerical algorithms from the previous sections with the help of four numerical examples. These include linear and non-linear stationary settings as well as a non-linear space-time test. All
examples are evaluated with the help of state-of-the-art measures by observing error reductions, reductions in the estimators, effectivity indices and indicator indices. These demonstrate the performance of goal-oriented adaptivity for different types of
discretizations. The numerical tests are computed with the open-source finite element libraries deal.II [75]–[77] and MFEM [212], [213]. Parts of our code developments are open-source on GitHub^{9}.

The purpose of this first numerical example is to illustrate some of the previous algorithmic and theoretical developments with the help of a numerical experiment including open-source code developments. We span from the problem statement, over discretization, numerical solution to goal functional evaluations on adaptively refined meshes using algorithms from Section 3.11. To this end, we employ a simple example, namely Poisson’s problem in two dimensions. In the numerical simulations, we also demonstrate trivial effects such as that the dual space must be richer than the primal function space for the error interpolations, as otherwise the error estimator is identically zero due to Galerkin orthogonality; cf. 23 .

The code is based on deal.II [75], [77] and the
current release version 9.5.1 [214]. Furthermore, this code is published open-source on GitHub^{10}. Conceptually, this code builds directly upon [62].

Let \(\Omega:=(0,1)^2\) be the domain with the Dirichlet boundary \(\partial\Omega\). Let \(f:\Omega\to\mathbb{R}\) be some volume force. Find \(u:\bar{\Omega}\to\mathbb{R}\) such that \[\begin{align} - \Delta u &= f \quad\text{in } \Omega, \\ u &= 0 \quad\text{on } \partial \Omega, \end{align}\] where \(f=-1\).

We define the function spaces, here \(U = V:= H^1_0(\Omega)\). Then, the weak form reads: Find \(u\in U\) such that \[\mathcal{A}(u)(\psi) = l(\psi) \quad\forall \psi\in V,\] with \[\mathcal{A}(u)(\psi):= \int_{\Omega} \nabla u\cdot \nabla \psi\, dx, \quad and \quad l(\psi) := \int_{\Omega} f\psi\, dx.\] The discrete problem is formulated on quadrilateral elements (Section 2.2.2) as decomposition of the domain \(\Omega\) and reads: Find \(u_{h}\in U_h\) such that \[\mathcal{A}(u_{h})(\psi_{h}) = l(\psi_{h}) \quad\forall \psi_{h}\in V_h,\] with \[\mathcal{A}(u_{h})(\psi_{h}):= \int_{\Omega} \nabla u_{h}\cdot \nabla \psi_{h}\, dx, \quad and \quad l(\psi_{h}) := \int_{\Omega} f\psi_{h}\, dx.\] For implementation reasons, the numerical solution is obtained within a Newton scheme in which the linear equations are solved with a direct solver (UMFPACK [215]). Clearly, the problem is linear and for this reason Newton’s method converges within one single iteration. The main reason using a Newton method is that the code can easily be extended to non-linear problems.

Furthermore, the goal functional is given as a point evaluation in the middle point: \[J(u) = u(0.5,0.5).\] For the later comparison, the reference value is determined as \[\label{eq95ref95value95736} u_{ref}(0.5,0.5) := -7.3671353258859554e-02,\tag{40}\] and was obtained on a globally refined super-mesh.

The adjoint problem is derived as in Equation 3 \[\text{Find } z\in V: \quad \mathcal{A}(\varphi,z) = J(\varphi) \quad\forall \varphi\in U,\] where we notice that both the left hand side and right sides are linear. Specifically, these are given as: \[\mathcal{A}(\varphi,z) = \int_{\Omega} \nabla \varphi\cdot \nabla z\, dx, \quad J(\varphi) = \varphi(0.5,0.5).\] The adjoint is solved similar to the primal problem, namely with Newton’s method (converging in one single iteration) and UMFPACK are employed again. Here, Newton’s method is not at all necessary since the adjoint is always linear as previously discussed. However, for implementation convenience, we employ the same Newton solver for the primal and adjoint problem, which constitutes solely our reason.

The objectives of our studies are:

Evaluating goal-oriented exact error \(J(u_{ref}) - J(u_{h})\);

Evaluating error estimator \(\eta\);

Computing \(I_\mathrm{eff}\) and \(I_{ind}\);

Employing PU-DWR for local mesh adaptivity;

Showing necessity of higher-order information of the adjoint \(z\) in \(z_{h}^{high} - z_{h}\) in the dual-weighted residual estimator;

Employing higher-order shape functions;

Employing higher-order PU function.

In this section, we conduct in total ten numerical experiments to investigate our previous objectives.

In this first numerical test (Table [tab951456951]), the adjoint has the same order as the primal solution, and due to Galerkin orthogonality (see e.g., 23 ), the estimator \(\eta\) is identical to zero and consequently, the adaptive algorithm stops after the first iteration, which is certainly not what we are interested in.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00
================================================================================
```

In this second experiment (Table [tab951456952]), we choose now a higher-order adjoint solution and obtain a typical adaptive loop. In the true (exact) error and the estimator \(\eta\), we observe quadratic convergence as to be expected. The \(I_\mathrm{eff}\) and \(I_{ind}\) are around the optimal value \(1\). These findings are in excellent agreement with similar results published in the literature.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 2.00e-02 2.00e-02 9.98e-01 9.98e-01
50 4.01e-03 4.03e-03 4.75e-03 1.00e+00 1.18e+00
162 9.27e-04 9.28e-04 1.17e-03 1.00e+00 1.26e+00
482 2.37e-04 2.44e-04 3.08e-04 1.03e+00 1.30e+00
1794 5.82e-05 5.91e-05 7.52e-05 1.02e+00 1.29e+00
4978 1.73e-05 1.92e-05 2.34e-05 1.11e+00 1.35e+00
14530 5.13e-06 6.02e-06 7.74e-06 1.18e+00 1.51e+00
================================================================================
```

In this third example (Table [tab951456953]), we increase the adjoint polynomial order and observe that we obtain roughly the same error tolerance of about \(10^{-6}\) with less degrees of freedom. The \(I_\mathrm{eff}\) is still optimal, while \(I_{ind}\) shows a slight overestimation with \(I_{ind} \approx 2\).

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 2.01e-02 3.43e-02 9.99e-01 1.71e+00
50 4.01e-03 4.01e-03 8.57e-03 1.00e+00 2.14e+00
162 9.27e-04 9.27e-04 1.99e-03 1.00e+00 2.14e+00
274 3.33e-04 3.98e-04 7.02e-04 1.20e+00 2.11e+00
994 8.45e-05 9.16e-05 1.67e-04 1.08e+00 1.98e+00
2578 2.39e-05 2.65e-05 4.64e-05 1.11e+00 1.94e+00
8482 6.00e-06 6.35e-06 1.21e-05 1.06e+00 2.01e+00
================================================================================
```

In this fourth numerical test (Table [tab951456954]), the results are close to the previous setting, which basically shows that such high-order adjoint solutions, do not increase necessarily anymore the error estimator and adaptivity.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 2.01e-02 4.64e-02 9.99e-01 2.31e+00
50 4.01e-03 4.01e-03 1.02e-02 1.00e+00 2.54e+00
162 9.27e-04 9.27e-04 2.45e-03 1.00e+00 2.65e+00
482 2.37e-04 2.44e-04 6.25e-04 1.03e+00 2.63e+00
802 9.72e-05 1.12e-04 2.10e-04 1.16e+00 2.16e+00
2578 2.39e-05 2.65e-05 5.21e-05 1.11e+00 2.18e+00
8290 6.12e-06 6.57e-06 1.37e-05 1.07e+00 2.24e+00
================================================================================
```

The fifth numerical experiment (Table [tab951456955]) has the inverted polynomial order in which Galerkin orthogonality is even more violated than in the \(Q_1/Q_1\) case. Our numerical results confirm the theory.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
50 4.66e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00
================================================================================
```

In the sixth test (Table [tab951456956]), we are in the same situation with equal-order polynomials as in the \(Q_1/Q_1\) case, and consequently, again due to Galerkin orthogonality the estimator \(\eta\) is zero.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
50 4.66e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00
================================================================================
```

We perform now a seventh experiment (Table [tab951456957]), which again works. In the last row, the \(I_\mathrm{eff}\) is off, likely to the reason that the reference value 40 is not accurate enough anymore. This shows that numerically obtained reference values must be computed with care.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
50 4.66e-05 2.31e-05 1.27e-03 4.96e-01 2.72e+01
162 1.98e-05 1.96e-05 1.09e-04 9.88e-01 5.47e+00
578 1.45e-06 1.44e-06 7.62e-06 9.96e-01 5.27e+00
1826 2.54e-08 3.31e-08 7.26e-07 1.31e+00 2.86e+01
6978 1.96e-09 7.72e-11 5.27e-08 3.94e-02 2.69e+01
================================================================================
```

In this eighth example (Table [tab951456958]), the results are similar to the previous test case, which basically confirms the fourth example, that higher-order adjoint do not contribute to better findings anymore for this specific configuration.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
50 4.66e-05 3.31e-05 1.33e-03 7.09e-01 2.85e+01
162 1.98e-05 1.98e-05 9.63e-05 1.00e+00 4.86e+00
578 1.45e-06 1.45e-06 6.88e-06 1.00e+00 4.75e+00
1826 2.54e-08 3.38e-08 6.88e-07 1.33e+00 2.71e+01
6978 1.96e-09 4.32e-11 5.00e-08 2.20e-02 2.55e+01
================================================================================
```

We finally conduct two tests with higher-order PU polynomial degrees. In Table [tab951456959], the \(I_\mathrm{eff}\) performs very well, while the indicator index shows over estimation of a factor \(4\).

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 2.00e-02 3.61e-02 9.98e-01 1.80e+00
50 4.01e-03 4.03e-03 1.30e-02 1.00e+00 3.24e+00
162 9.27e-04 9.28e-04 3.70e-03 1.00e+00 4.00e+00
482 2.37e-04 2.40e-04 1.05e-03 1.01e+00 4.40e+00
1106 7.73e-05 8.52e-05 3.27e-04 1.10e+00 4.23e+00
3810 2.01e-05 2.09e-05 8.62e-05 1.04e+00 4.29e+00
13250 5.39e-06 6.01e-06 2.39e-05 1.11e+00 4.43e+00
================================================================================
```

In this final test (Table [tab9514569510]), the indicator index shows an overestimation of a factor about \(8\). In terms of the true error and estimated error as well as the \(I_\mathrm{eff}\), the results are close to being optimal, while the indicator index shows an over estimation. With respect to a higher computational cost in computing the PU, these findings suggest that a low-order PU for this configuration is sufficient.

```
================================================================================
Dofs Exact err Est err Est ind Eff Ind
--------------------------------------------------------------------------------
18 2.01e-02 2.00e-02 6.22e-02 9.98e-01 3.10e+00
50 4.01e-03 4.03e-03 2.18e-02 1.00e+00 5.43e+00
162 9.27e-04 9.28e-04 6.14e-03 1.00e+00 6.62e+00
482 2.37e-04 2.32e-04 1.69e-03 9.78e-01 7.11e+00
1602 6.21e-05 6.95e-05 4.61e-04 1.12e+00 7.42e+00
5618 1.59e-05 1.74e-05 1.22e-04 1.10e+00 7.71e+00
19602 4.32e-06 5.16e-06 3.31e-05 1.19e+00 7.67e+00
================================================================================
```

Finally, from the `vtk`

raw data, using `visit`

[216], we show some graphical solutions in Figure 7 that
illustrate the performance of our algorithms.

In the second example, we consider the homogeneous Dirichlet boundary value problem for a non-linear elliptic partial differential equation in the plane domain \(\Omega=(0,5)\times (0,3) \setminus \left((1,2)^2 \cup (3,4)\times (1,2) \right) \subset \mathbb{R}^2\) visualized together with the initial mesh and the quantities of interest in Figure 8.

We look for some function \(u\) such that \[\label{eq:Ex2:pde} -\text{div}(\nu(|\nabla u|)\nabla u)=f\; \text{in}\;\Omega \quad \text{and} \quad u=0 \; \text{on}\; \partial\Omega,\tag{41}\] with \(\nu(s)=2+\text{arctan}(s^2)\) and given right-hand side \(f=10\). Such kind of non-linear partial differential equations arise, for instance, in magneto-static where \(\nu\) is the reluctivity that non-linearly depends on the gradient of \(u\); see, e.g., [217]. The weak form of 41 reads as follows: Find \(u \in U = V = H_0^1(\Omega) =\mathring{W}_2^1(\Omega)\) such that \[\label{eq:Ex2:WeakForm} \mathcal{A}(u)(v) = 0 \quad \forall \,v \in V\tag{42}\] that can equivalently be written as non-linear operator equation \[\label{eq:Ex2:OperatorEquation} \mathcal{A}(u) = 0 \quad \text{in} \; V^* = U^* = H^{-1}(\Omega),\tag{43}\] where the non-linear operator \(\mathcal{A}(u): H_0^1(\Omega) \rightarrow H^{-1}(\Omega)\) is well defined by the variational identity \[\label{eq:Ex2:OperatorDefinition} \mathcal{A}(u)(v) := \langle \nu(|\nabla u|)\nabla u,\nabla v \rangle-\langle f,v \rangle \quad \forall\, v, u \in H_0^1(\Omega).\tag{44}\]

Since the non-linear operator is strongly monotone and Lipschitz continuous, the operator equation 43 and the equivalent weak or variational formulation 42 have a unique solution \(u \in H_0^1(\Omega)\); see, e.g., [171], [217]. The differentiability properties of \(\nu\) yield the corresponding Fréchet differentiability properties of the non-linear operator \(\mathcal{A}\). For instance, the first Fréchet derivative of \(\mathcal{A}\) at some \(u \in U\) is nothing but the bounded, linear operator \(\mathcal{A}'(u): U \rightarrow V^*\) defined by the variational identity \[\label{eq:Ex2:FrechetDerivative} \mathcal{A}'(u)(v,w) := \langle (\nu(|\nabla u|) I + \frac{\nu'(|\nabla u|)}{|\nabla u|}\nabla u (\nabla u)^T) \nabla v,\nabla w \rangle \quad \forall\, v, w \in H_0^1(\Omega).\tag{45}\]

In the following, we are interested in the following quantities of interest:

the flux on \(\Gamma_{\text{Flux}}:=\{0\}\times (0,3)\): \(J_1(u):=\int_{\Gamma_{\text{Flux}}} \nabla u \cdot n \text{ d}s_x,\)

the point evaluation at \(x_0=(0.2,0.2)\): \(J_2(u):=u(x_0),\)

the point evaluation at \(x_1=(0.9,0.1)\): \(J_3(u):=u(x_1).\)

We can assume that the functionals are well defined at the solution \(u\). As an error weighting function, we choose the function presented in 32 , i.e. \[\mathfrak{E}(x,m)=\sum_{i=1}^N \frac{x_i}{|m_i|},\] with \(m=(J_1(u_h),J_2(u_h),J_3(u_h))\). Then, the combined functional follows from 36 .

The basic finite element discretization is performed by \(Q_1\) finite elements, whereas \(Q_2\) finite elements are used to construct the enriched spaces. For instance, the finite element scheme for solving the primal problem 42 reads as follows: Find \(u_h \in U_h \subset U\) such that \[\label{eq:Ex2:FinitElementScheme} \mathcal{A}(u_h)(v_h) = 0 \quad \forall \,v_h \in V_h=U_h,\tag{46}\] where \(U_h\) is spanned by the \(Q_1\) basis functions corresponding to the internal nodes.

Figure 8 shows the initial mesh with three quantities of interest, and the adaptively refined mesh after 25 refinement steps when the adaptivity is driven by the combined functional \(J_\mathfrak{E}\). As expected, we observe heavy mesh refinement where the three functionals are computed. The 8 interior corners with a re-entrant angle of \(\frac{3}{2}\pi\) give rise to singularities that pollute the solution elsewhere. However, the mesh refinement around these singularities depends on the distance to the places where we compute the functionals. The strongest refinement is observed in the lower left re-entrant corner, whereas the initial mesh is almost not refined in the upper right re-entrant corner that has the largest distance from the place where we evaluate the functionals. Figure 9 shows the finite element solution \(u_h\) and the localized error estimator after 25 adaptive refinements of the initial mesh. We observe that \(u_h\) is recovered more accurately in regions where the accuracy is needed for computing accurate values of the functionals. This goal-oriented adaptive refinement aiming at the joint accurate computation of local functionals is very different from the adaptive mesh refinement driven by global values like the \(H^1\)-norm of the discretization error in the residual a posteriori error estimator; see, e.g. [7]. Indeed, in Figure 13, we display the finite element solution \(u_h\) and the mesh after 25 adaptive refinements of the initial mesh when the adaptive refinement is driven by the \(L_2\)-norm, i.e., by the functional \(J(u) = \|u\|^2_{L_2(\Omega)}\) that is a global functional too. As expected, the refinement is located at the 8 singularities in the 8 re-entrant corners. In addition to this, Figures [fig:32Ex2:32IeffL2] and 14 illustrates the numerical behavior of the effectivity indices \(I_{\mathrm{eff}}\), \(I_{\mathrm{eff,a}}\), and \(I_{\mathrm{eff,p}}\) for \(J(u)\) and the decay of the absolute error for the functional \(J\) in the cases of adaptive and uniform refinements, respectively. This error decay behaves like \(O(h^{3/2}) = O(\text{DoFs}^{-3/4})\) for uniform mesh refinement due to the singularities, whereas \(O(\text{DoFs}^{-1})\), that is equivalent to \(O(h^{2})\) on a uniform mesh, is observed for the adaptive finite element procedure driven by the functional \(J\).

Furthermore, Figure [fig:32Ex2:32Ieff] shows the numerical behavior of the effectivity indices \(I_{\mathrm{eff}}\), \(I_{\mathrm{eff,a}}\), and \(I_{\mathrm{eff,p}}\) for the combined functional \(J_\mathfrak{E}\). We see that all three effectivity indices are close to \(1\), and \(I_{\mathrm{eff}}\) is practically almost \(1\) for more than \(10000\) DoFs. Figure 10 provides the error decay of the relative errors for \(J_1\), \(J_2\), \(J_3\), and the absolute error for \(J_\mathfrak{E}\). Figures [fig:32Ex2:32Error:32J1], 11, and [fig:32Ex2:32Error:32J3] compare the error decay of the uniform and adaptive refinements for the functionals \(J_1\), \(J_2\), and \(J_3\), respectively, whereas Figure 12 shows that the error estimator \(\eta_h\) is practically identical with the error for \(J_\mathfrak{E}\), and both decay like \(O((\text{DoFs})^{-1})\) as expected in the adaptive case.

In this third example, we consider the flow around a cylinder benchmark from [218]. We notice that this example is an extension of our prior study [73].

The domain \(\Omega:=(0,L)\times (0,H)\setminus \mathcal{B}\) where \(L=2.2\), \(H=0.41\) and \(\mathcal{B}:=\{x \in \mathbb{R}:|x-(0.2,0.2)|<0.05\}\). The domain as well as the boundary conditions are depicted in Figure [fig:32Ex3:32Omega43Boundary].

Find \(u:=(v,p)\) such that \[\begin{align} - \nabla\cdot\left (\nu \nabla v\right)+(v \cdot \nabla) v + \nabla p =& 0 \qquad\text{in }\Omega,\nonumber\\ \nabla \cdot v =& 0 \qquad\text{in }\Omega,\nonumber\\ \end{align}\] with \(\nu =10^{-3}\). Furthermore, we consider the following boundary conditions \[\begin{align} v=&0 \qquad\text{on } \Gamma_{\mathrm{no-slip}}:=\partial \mathcal{B} \cup \Gamma_{0}, \\ v=&\hat{v} \qquad \text{on } \Gamma_{\mathrm{inflow}}, \nonumber\\ \nu \frac{\partial v}{\partial \vec{n}} - p \cdot \vec{n }=& 0 \qquad\text{on } \Gamma_{\mathrm{outflow}}, \nonumber\\ \end{align}\] where \(\Gamma_{\mathrm{inflow}}:=\{0\}\times (0,H)\), \(\Gamma_{\mathrm{outflow}}:=\{L\}\times (0,H)\), and \(\Gamma_{0}:=(0,L)\times\{0,H\}\). The inflow profile \(\hat{v}\) is \(\hat{v}(x_1,x_2):=1.2x_2(H-x_2)/H^2\). Additionally, we mention that the pressure is uniquely determined due to the do-nothing condition on \(\Gamma_{\mathrm{outflow}}\); see [219].

The weak form of the problem is given by: Find \(u:=(v,p) \in (\hat{v},0) + U:=H^1_{\Gamma_\mathrm{no-slip}\cup\Gamma_{\mathrm{inflow}}}(\Omega) \times L^2(\Omega)\) such that \[\mathcal{A}(u)(\psi):=(\nu \nabla v, \nabla \psi_v)+(\nabla v \cdot v, \psi_v)+(p, \operatorname{div}(\psi_v))+(\operatorname{div}(u),\psi_p)\] for all test functions \((\psi_v,\psi_p)=:\psi \in V=U\), where \(H^1_{\Gamma_\mathrm{no-slip}\cup\Gamma_{\mathrm{inflow}}}(\Omega):=\{v \in H^1(\Omega): u =0 \text{ on }\Gamma_\mathrm{no-slip}\cup\Gamma_{\mathrm{inflow}} \}\).

The domain is decomposed into quadrilateral elements; cf. Figure 16 (left). As discretization, we use \(Q_2^d\) finite elements for the velocity \(v\) and \(Q_1\) finite elements for the pressure \(p\). For the enriched space we use \(Q_4^d\) finite elements for the velocity \(v\) and \(Q_2\) finite elements for the pressure \(p\). This discretizations are inf-sup stable as shown in [220].

We are interested in the following functionals of interest:

the pressure difference between \(x_1=(0.15,0.2)\) and \(x_2=(0.25,0.2)\): \[J_1(u):=p(x_1)-p(x_2),\]

the drag at the ball \(\mathcal{B}\): \[J_2({u}):= 500\int_{\partial \mathcal{B}} \left[\nu \frac{\partial u}{\partial \vec{n}} - p \vec{n }\right]\cdot \vec{e}_1\,\mathrm{ d}s_{(x)},\]

the lift at the ball \(\mathcal{B}\): \[J_3({u}):= 500\int_{\partial \mathcal{B}} \left[\nu \frac{\partial u}{\partial \vec{n}} - p \vec{n }\right]\cdot \vec{e}_2\,\mathrm{ d}s_{(x)},\]

the square of the magnitude of velocity at \(x_3=(2.2,0.205)\): \[J_4(u):=|v|^2(x_3),\]

where \(\vec{e}_1 = (1,0)\) and \(\vec{e}_2 = (0,1)\). For \(J_1\), \(J_2\), \(J_3\), we use the reference values from [221], and for \(J_4\) we used uniform \(p\)-refinement on an adaptively refined grid. The reference values are :

\(J_1(u)=0.11752016697\),

\(J_2(u)=5.57953523384\),

\(J_3(u)=0.010618948146\),

\(J_4(u)=0.088364291405\).

As in the previous example, we choose the error weighting function 32 , i.e. \[\mathfrak{E}(x,m)=\sum_{i=1}^N \frac{x_i}{|m_i|},\] with \(m=(J_1(u_h),J_2(u_h),J_3(u_h),J_4(u_h))\).

In Figure 18, we can see that the error functional \(J_\mathfrak{E}\) dominates the relative error in the single functionals. At the beginning of the refinement procedure, the error in \(J_4\) is the smallest of all the functionals. However, it requires different local refinement needs than the other functionals. Therefore, the error in \(J_4\) is not reduced in the first \(10\) refinement steps. In fact, the elements around the point \(x_3\) are not refined a single time during the first \(9\) refinement steps, cf. Figure 16. In Figure 17, there is already refinement in the area around \(x_3\) but not at \(x_3\) after \(11\) adaptive refinements, but a lot at \(x_3\) after \(17\) adaptive refinements. We also observe that the error of \(J_4\) is reduced in those refinement steps, cf. Figure 18 and Figure 20. The effectivity index \(I_\mathrm{eff}\) in Figure [fig:32Ex3:32Ieff] is between 0.5 and 2.5. The adjoint and primal effectivity indices have a similar value. In Figures [fig:32Ex3:32Error32J1] - 20, we observe that adaptive refinement has a smaller error than uniform refinement after \(30\, 000\) DoFs. Additionally, the convergence rate is better for adaptive refinement for all \(4\) functionals. At \(100\, 000\) DoFs, all errors are approximately one order of magnitude lower for adaptive refinement than for uniform refinement.

For our final application, we consider a non-linear evolutionary PDE, in particular, the regularized parabolic \(p\)-Laplace PDE. Let \(\Omega \subset \mathbb{R}^d\) again denote our spatial domain, and let \(T > 0\) be the final time horizon. To illustrate the performance of the space-time multi-goal oriented finite element method, we perform numerical experiments for \(d=2\) and \(d=3\). Here we consider problems with spatial singularities, but that are smooth in time. Depending on the goal functionals, we expect refinement not only in the domains of interest, but also towards the spatial singularities. In particular, we consider the classical L-shaped domain for the 2D+1 (\(d=2\)) case, and the Fichera corner domain for the 3D+1 (\(d=3\)) case.

The problem reads: Find \(u\) such that \[\begin{gather} \tag{47} \partial_t u - \operatorname{div}_x( (|\nabla_x u|^2 + \epsilon^2)^{(p-2)/2} \nabla_x u) = f\text{ in } Q = \Omega \times (0,T),\\ {u = 0 \text{ on } \Sigma_0=\Omega\times\{0\} \quad\text{and}\quad u = 0\text{ on } \Sigma=\partial \Omega\times(0,T).} \tag{48} \end{gather}\] Note that this form of the regularized \(p\)-Laplacian is motivated by the Carreau model. In the numerical experiment, we choose \(p=4\) and \(\epsilon = 10^{-2}\). We can formulate the above problem again in the abstract framework introduced in Section 2. Indeed, with the choices \(V = L_p(0,T; \mathring{W}^{1}_{p}(\Omega))\) and \(U = \{ v \in V: \partial_t v \in V^*, v = 0\;\text{on}\;\Sigma_0\}\), we arrive at the operator equation: Find \(u \in U\) such that \[\label{eq:ex4:OperatorEquation} \mathcal{A}(u) = 0\quad in\; V^*,\tag{49}\] where \[\langle \mathcal{A}(u), v \rangle = \langle\partial_t u, v\rangle + \langle (|\nabla_x u|^2 + \epsilon^2)^{(p-2)/2} \nabla_x u, \nabla_x v \rangle - \langle f, v\rangle \quad for all\, v\in V,\] \(L_p(0,T; X)\) denotes the Bochner space of \(L_p\)-functions that map from the interval \((0,T)\) to \(X\), and \(\mathring{W}^{1}_{p}(\Omega) = \{v \in L_p(\Omega): \nabla v \in (L_p(\Omega))^d, v = 0 \;\text{on}\;\partial \Omega \}\). The existence and uniqueness of the solution of the operator equation 49 follow from standard monotonicity arguments; see, e.g.,[171], [222].

Evolutionary problems like 47 –48 are usually solved by the method of lines, i.e.we perform semi-discretization separately for the spatial domain \(\Omega\) and the time interval \((0,T)\). An alternative ansatz is an all-at-once discretization, i.e.we directly discretize the space-time cylinder \(Q\subset \mathbb{R}^{d+1}\). In this work, we will focus on the latter approach. We will briefly describe the key points, for more details on space-time methods see, e.g., [34]. Starting with a space-time Galerkin ansatz with finite dimensional subspaces \(U_h \subset U\) and \(V_h = U_h \subset V\), we arrive at the discretized problem: Find \(u_h \in U_h\) such that \[\langle \mathcal{A}(u_h), v_h \rangle = 0 \quad \text{for all}\;v_h \in V_h.\] The above discrete problem has a unique solution; see [63] and the references therein. Finally, in order to solve the non-linear systems of equations, we again apply Newton’s method, where in each iteration, we have to solve the linear problem: Find \(w_h \in U_h\) such that \[\langle \mathcal{A}'(u_h)w_h, v_h \rangle = -\langle \mathcal{A}(u_h), v_h \rangle \quad \text{for all}\;v_h \in V_h.\] The above equation is uniquely solvable. Hence, Newton’s method is well defined; see [63] and the references therein.

We apply Algorithm 6, i.e.the enriched primal and adjoint problems are solved by Newton’s method and preconditioned GMRES, respectively.

We implemented the space-time finite element method using the finite element library MFEM [212]. We use a \(P_1\) conforming ansatz space for the 2D+1 (\(d=2\)) and 3D+1 (\(d=3\)) case, i.e.we use a simplicial decomposition \(\mathcal{T}_h\) of the \((d+1)\)-dimensional space-time cylinder \(Q \subset \mathbb{R}^{d+1}\).

We consider two goal functionals:

Let \(\Omega_I := (-15/16,-11/16)\times (11/16,15/16)\) and \(\Omega_I := (11/16,15/16)^3\) for \(d=2\) and \(d=3\), respectively, and \[J_1(u) := \int_{0.5}^{0.75}\!\int_{\Omega_I}\! |\nabla u(x,t)|^4 \;\mathrm{d}x\mathrm{d}t,\] i.e.we are interested in the energy in some subdomain far away from the corner singularity, and

the “average” of the solution at final time \(T\), i.e. \[J_2(u) = \int_\Omega\! u(x,T) \;\mathrm{d}x.\]

As an error weighting function, we choose the function presented in 32 . Thus, our combined functional has the form \[J_{\mathfrak{E}}(v) = \frac{|J_1(u_h^{(2)}) - J_1(v)|}{J_1(u_h)} + \frac{|J_2(u_h^{(2)}) - J_2(v)|}{J_2(u_h)}.\]

First we consider the L-shaped domain \(\Omega = (-1,1)^2 \setminus \{ [0,1)\times(-1,0]\}\) as depicted in the left of Figure 21. Here, we know the exact solution of the corresponding
elliptic problem^{11} and we manufacture a time-dependent solution for the parabolic \(p\)-Laplace initial-boundary value problem, which is
given by \[u(x_1,x_2,t) = \frac{3}{2} (1-x_1)^2(1-x_2)^2 r^{2/3} \sin(\frac{2}{3}\theta)\,\sin(t),\] where \((r,\theta)\) denote the polar representation of the spatial coordinates \((x_1,x_2)\). The right-hand side is computed accordingly.

We first consider the efficiency indices as presented in the left of Figure 22. Here, we observe that after some initial oscillations, the combined efficiency index \(I_{\mathrm{eff}}\) is very close to \(1\). These initial oscillations most likely occur as the integration domain \(Q_I = \Omega_I \times (0.5,0.75)\) is not exactly captured in the initial mesh; see Figure 24 (top right). Next, we consider the convergence of the error for the functionals. In the right of Figure 22, we observe that the error in the combined functional \(J_{\mathfrak{E}}\) and the error in the “average” \(J_2\) converge with almost \(\mathcal{O}(\mathrm{DoFs}^{-2/3})\), i.e.quadratic. While the error for the “energy” integral \(J_1\) is more oscillatory, we observe an overall almost quadratic convergence rate of around \(\mathcal{O}(\mathrm{DoFs}^{-2/3})\). In Figure 23, we present the error of the individual functionals, and compare the adaptive refinements with uniform refinements. We again observe the almost quadratic convergence of the adaptive refinements, while uniform refinements on the one hand need more DoFs to reach a specific error, and on the other hand have a worse convergence rate in the pre-asymptotic range.

Next, we consider the meshes obtained from the adaptive algorithm. In Figure 24, we present the space-time mesh in the first two rows, and three cuts with the \((x_1,x_2)\)-plane in the bottom row. In the upper two rows, we additionally present intersections of the space-time cylinder with the \((x_1,t)\)-plane at \(x_2 = 13/16\) and with the \((x_2,t)\)-plane at \(x_1 = -13/16\) in the right column. We indicate these cuts also in the left columns by thick red lines. In the upper row, we show the initial meshes, and in the middle row the meshes obtained after 18 adaptive refinements using Algorithm 6. As expected, we observe that the top of the space-time cylinder is refined, with special focus on the corner singularity. Moreover, we observe refinements toward the integration domain \(Q_I\) of \(J_2\), which is indicated with red shading. Since \(\Omega_I\) is quite far away from the singularity, the pollution effect is quite mild. This can also be observed in the bottom row of Figure 24, as the refinements towards the re-entrant corner are far less in the first two columns compared with the third column. From left to right, the cuts are located at \(t=0.5\), \(t = 0.625\), and \(t = 1\), respectively.

For the Fichera corner \(\Omega = (-1,1)^3 \setminus(-1,0]^3\) (see Figure 21 right), in addition to the corner singularity at the origin, there can also be edge singularities along the
re-entrant edges. In particular, we choose the right hand side ^{12} \[f(x,t) = \frac{\sin(t)}{\sqrt{x_1^2 + x_2^2 + x_3^3}}.\] In this case,
to the best of our knowledge, there is no explicit form of the solution \(u\). Nevertheless, we will again apply our adaptive multi-goal Algorithm 6, with the same
functionals as for the L-shape domain. For the first functional, we now use the three-dimensional \(\Omega_I = (11/16, 15/16)^3\), see the right of Figure 21 for a visualization. Instead
of presenting the convergence behavior of the error, we plot the convergence behavior of the error indicator \(\eta_h^{(2)}\). As we have observed in the previous examples, after some refinements \(\eta_h^{(2)}\) provides a good estimate for the actual error.

In Figure 25, we present the parts of the error estimator we use to drive the adaptive mesh refinement. We observe that the primal and adjoint parts of the error estimator both decay with a rate of \(\mathcal{O}(\mathrm{DoFs}^{-2/4})\), i.e. quadratic convergence. In Figure 26, we compare the convergence history of the functional values obtained by uniform refinements and adaptive refinements. As a reference value for the “exact” functional values, we plot the value of the functionals evaluated at the enriched solution \(u_h^{(2)}\) on the finest adaptive mesh.

Finally, we present intersections of the four dimensional space-time mesh with the \((x_1,x_2,x_3)\)-hyper-plane in Figure 27. From top to bottom, the cuts were performed at \(t=0\), \(t=0.5\) and \(t=1\), respectively. In the left column, we present the initial mesh, and in the right column the mesh after \(18\) adaptive refinements. We observe that the mesh refinements toward the edge and corner singularities are much more prominent at final time \(t=1\) as, e.g.,at time \(t=0.5\). This behavior is similar to our earlier observations for the 2D+1 case.

In this work, we reviewed single- and multigoal-oriented error control and adaptivity. Our notation was kept in general terms such that stationary and non-stationary (space-time) situations are covered. First, we explained in detail single goal-oriented error estimation with the help of the dual-weighted residual method. The error estimators cover both discretization and non-linear iteration errors. Therein, some new theoretical results were presented as well. Prior efficiency and reliability results only require one saturation assumption, rather a strengthened condition. Next, we considered non-standard discretizations such as stabilization terms and non-consistencies (e.g., classical finite difference based time-stepping schemes) or non-conformal methods (e.g., finite differences or neural network approximations) of numerical schemes in the frame of goal-oriented error estimation. Third, we concentrated on multigoal-oriented error estimation in which we have put a lot of efforts in the last eight years due to advances and increasing interest of solving multiphysics partial differential equations and coupled variational inequality systems. Besides theoretical results, we provided several adaptive algorithms for single and multiple goal functional evaluations. In order to substantiate our developments, four numerical examples were designed and computationally analyzed. In the first example, namely Poisson’s problem, the implementation is fully open-source and follows the classical structure of deal.II tutorial steps. Therefore, this example has a flavor of educational purpose. In the last example, recent work on the space-time PU-DWR method is further extended to multiple-goal functionals and singularities in the solution. Ongoing and future work is concerned with the extension to cover more terms in the error estimators that besides discretization and non-linear iteration errors, as well linear iteration errors, model errors, and model order reduction techniques are covered. In conjunction with high-performance parallel computing these yield important components to continue to solve efficiently and with desired accuracies multiphysics problems with physics-based discretizations and fast physics-based numerical solvers.

BE and TW have been supported by the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453). Furthermore, BE greatly acknowledges his current Humboldt Postdoctoral Fellowship funding. Furthermore, we would like to thank the Johann Radon
Institute for Computational and Applied Mathematics (RICAM) for providing computing resource support of the high-performance computing cluster Radon1 ^{13}.

[1]

I. Babuška and W. C. Rheinboldt, “Error estimates for adaptive finite element computations,” *SIAM J. Numer. Anal.*, vol. 15, no. 4, pp. 736–754, OPT, 1978 , OPT, doi:
0036-1429/78/1504-0007 , OPTisbn = {}, OPT.

[2]

W. Dörfler, “A convergent adaptive algorithm for poisson’s equation. , f,” *SIAM Journal on Numerical Analysis , journal = SIAM J. Numer. Anal.*, vol. 33,
no. 3, pp. 1106–1124, 1996.

[3]

C. Carstensen and R. Verfürth, “Edge residuals dominate a posteriori error estimates for low order finite element methods,” *SIAM J. Numer. Anal. , FJOURNAL = SIAM
Journal on Numerical Analysis*, vol. 36, no. 5, pp. 1571–1587, 1999, doi: 10.1137/S003614299732334X.

[4]

P. Morin, R. H. Nochetto, and K. G. Siebert, “Data oscillation and convergence of adaptive FEM,” *SIAM J. Numer. Anal. , FJOURNAL = SIAM Journal on
Numerical Analysis*, vol. 38, no. 2, pp. 466–488, 2000, doi: 10.1137/S0036142999360044.

[5]

P. Binev, W. Dahmen, and R. DeVore, “Adaptive finite element methods with convergence rates,” *Numer. Math. , FJOURNAL = Numerische Mathematik*, vol. 97,
no. 2, pp. 219–268, 2004, doi: 10.1007/s00211-003-0492-7.

[6]

R. Stevenson, “Optimality of a standard adaptive finite element method,” *Found. Comput. Math. , FJOURNAL = Foundations of Computational Mathematics. The Journal of
the Society for the Foundations of Computational Mathematics*, vol. 7, no. 2, pp. 245–269, 2007, doi: 10.1007/s10208-005-0183-0.

[7]

isbn=9780471967958. Verfürth R., *A review of a posteriori error estimation and adaptive mesh-refinement techniques*.
Wiley-Teubner, 1996.

[8]

K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, “Introduction to adaptive methods for differential equations , BOOKTITLE =
Acta numerica, 1995, SERIES = Acta Numer.” Cambridge Univ. Press, Cambridge, 1995 , MRCLASS = {65L70}, MR, pp. 105–158.

[9]

R. Becker and R. Rannacher, “An optimal control approach to a posteriori error estimation in finite element methods,” *Acta Numer. , FJOURNAL = Acta
Numerica*, vol. 10, no. 2009692 , MRREVIEWER = Wen Bin Liu, pp. 1–102, 2001, [Online]. Available: https://doi.org/10.1017/S0962492901000010.

[10]

L. Chamoin and F. Legoll, “An introductory review on a posteriori error estimation in finite element computations,” *SIAM Review*, vol. 65, no. 4, pp. 963–1028, 2023,
doi: 10.1137/21M1464841.

[11]

M. Ainsworth and J. T. Oden, *A posteriori error estimation in finite element analysis , SERIES = Pure and Applied
Mathematics*. John Wiley & Sons, New York, 2000, pp. xx+240, ISBN = 0-471-29411-X, MRCLASS = 65–02 (65N15), MR.

[12]

I. Babuška and O. T. Strouboulis ALTeditor =, *The finite element method and its reliability*, vol., OPT. Oxford University
Press, 2001.

[13]

[14]

O. W. Han ALTeditor =, *A posteriori error analysis via duality theory: With applications in modeling and numerical
approximations*, vol., OPT. Springer, 2004.

[15]

P. Neittaanmäki and O. S. Repin ALTeditor =, *Reliable methods for computer simulation: Error control and posteriori estimates*,
vol., OPT. Elsevier, 2004.

[16]

S. Repin, *A posteriori estimates for partial differential equations , series = RSCAM, OPTnote = Radon Series on Computational and Applied
Mathematics*, vol. 4. de Gruyter, 2008 , address = {Berlin}, OPT, p. xii+316.

[17]

K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, *Computational differential equations*. Cambridge University Press, 2009 ,
note={http://www.csc.kth.se/~jjan/private/cde.pdf}.

[18]

[19]

S. C. Brenner and L. R. Scott, *The mathematical theory of finite element methods , series = Texts in Applied
Mathematics*, vol. 15 , edition = Third. Springer, New York, 2008, pp. xviii+397, isbn = 978-0-387-75933-3, mrclass = 65-01 (65-02), mr.

[20]

[21]

P. G. Ciarlet, *Finite element method for elliptic problems*. Society for Industrial; Applied Mathematics , address = Philadelphia, PA, USA, 2002 , isbn =
{0898715148}.

[22]

T. J. R. Hughes, *The finite element method*. Dover Publications , address = , isbn = , series =, 2000, edition = {}.

[23]

[24]

A. Ern and J.-L. Guermond, *Finite elements i-III , OPTsubtitle = Approximation and Interpolation*. Springer, 2021.

[25]

T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,”
*Comput. Methods Appl. Mech. Eng.*, vol. 194 , OPT, pp. 4135–4195, OPT, 2005 , OPT, doi:, OPT.

[26]

J. A. Cottrell, T. J. R. Hughes, and citeulike =. 10673385. Bazilevs Y., *Isogeometric analysis: Toward integration of CAD and FEA*. John Wiley &
Sons, 2009.

[27]

“Mathematical analysis of variational isogeometric methods,” *Acta Numer. , FJOURNAL = Acta Numerica*, vol. 23, no. 3202239 , MRREVIEWER = Marius Ghergu,
pp. 157–287, 2014, doi: 10.1017/S096249291400004X.

[28]

C. Johnson, *Numerical solution of partial differential equations by the finite element method , NOTE = Reprint of the 1987 edition*. Dover
Publications, Inc., Mineola, NY, 2009, pp. ii+279, ISBN = 978-0-486-46900-3; 0-486-46900-X, MRCLASS = 65–01 (65M60 65N30), MR.

[29]

T. J. R. Hughes and G. M. Hulbert, “Space-time finite element methods for elastodynamics: Formulations and error estimates,” *Comp. Methods Appl. Mech. Engrg.*, vol.
66 , OPT, pp. 339–363, OPT, 1988 , OPT, doi:, OPTisbn = {}, OPT.

[30]

V. Thomée, *G alerkin finite element methods for parabolic problems*. Springer, 1997.

[31]

O. Steinbach, “Space-time finite element methods for parabolic problems,” *Comput. Methods Appl. Math.*, vol. 15, no. 4, pp. 551–566, OPT, 2015 , OPT, doi:, OPTisbn = {}, OPT.

[32]

R. H. Nochetto, S. A. Sauter, and Ch. Wieners, “Space-time methods for time-dependent partial differential equations,” *Oberwolfach reports*, vol. 14, no. 1, pp.
863–947, 2017, doi: 10.4171/OWR/2017/15.

[33]

O. Steinbach and H. Yang, “Space-time finite element methods for parabolic evolution equations: Discretization, a posteriori error
estimation, adaptivity and solution , booktitle = Space-Time Methods: Application to Partial Differential Equations, OPTcrossref = , OPT,” vol. 25 , OPT, de Gruyter, 2019 , OPTeditor = {}, pp. 207–248.

[34]

[35]

T. E. Tezduyar and K. Takizawa, “Space?time computations in practical engineering applications: A summary of the 25-,” *Comp. Mech.*, vol. 63, pp. 747–753, year=2019,
history.

[36]

S. Larsson, R. H. Nochetto, S. A. Sauter, and Ch. Wieners, “Space-time methods for time-dependent partial differential equations,” *Oberwolfach reports*, vol. 6, no.
1, pp. 1–80, 2022, doi: 10.4171/OWR/2022/6.

[37]

R. H. Nochetto, A. Schmidt, K. G. Siebert, and A. Veeser, “Pointwise a posteriori error estimates for monotone semi-linear equations,” *Numer. Math. , FJOURNAL =
Numerische Mathematik*, vol. 104, no. 4, pp. 515–538, 2006, doi: 10.1007/s00211-006-0027-0.

[38]

N. A. Pierce and M. B. Giles, “Adjoint and defect error bounding and correction for functional estimates,” *J. Comput. Phys. , issue_date = November
2004*, vol. 200, no. 2, pp. 769–794, Nov. 2004, doi: 10.1016/j.jcp.2004.05.001 , acmid = {1044303}.

[39]

S. Prudhomme, J. T. Oden, T. Westermann, J. Bass, and M. E. Botkin, “Practical methods for a posteriori error estimation in engineering applications,” *Int. J. Numer.
Methods Eng. , FJOURNAL = International Journal for Numerical Methods in Engineering*, vol. 56, no. 8, pp. 1193–1224, 2003, doi: 10.1002/nme.609.

[40]

M. B. Giles and E. Süli, “Adjoint methods for PDEs: A posteriori error analysis and postprocessing by duality,” *Acta Numer. , FJOURNAL = Acta
Numerica*, vol. 11, no. 2009374 , MRREVIEWER = Lucas Jódar, pp. 145–236, 2002, [Online]. Available: https://doi.org/10.1017/S096249290200003X.

[41]

M. Braack and A. Ern, “A posteriori control of modeling errors and discretization errors,” *Multiscale Model. Simul. , FJOURNAL = Multiscale Modeling &
Simulation. A SIAM Interdisciplinary Journal*, vol. 1, no. 2, pp. 221–238, 2003, [Online]. Available: https://doi.org/10.1137/S1540345902410482.

[42]

M. Giles and N. A. Pierce, “Analysis of adjoint error correction for superconvergent functional estimates.” 2001 , note = {Oxford University Computing Laboratory Report NA
01/14}.

[43]

N. A. Pierce and M. B. Giles, “Adjoint recovery of superconvergent functionals from PDE approximations,” *SIAM Rev. , FJOURNAL = SIAM
Review*, vol. 42, no. 2, pp. 247–264, 2000, doi: 10.1137/S0036144598349423.

[44]

S. Prudhomme and J. T. Oden, “On goal-oriented error estimation for elliptic problems: Application to the control of pointwise errors , ANNOTE = New advances
in computational methods (Cachan, 1997),” *Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 176, no. 1–4, pp. 313–331, 1999, [Online]. Available: https://doi.org/10.1016/S0045-7825(98)00343-0.

[45]

R. Rannacher and F.-T. Suttmeier, “A feed-back approach to error control in finite element methods: Application to linear elasticity,” *Comput. Mech. , FJOURNAL =
Computational Mechanics*, vol. 19, no. 5, pp. 434–446, 1997, doi: 10.1007/s004660050191.

[46]

R. Becker and R. Rannacher, “A feed-back approach to error control in finite element methods: Basic analysis and examples,” *East-West J. Numer. Math.*, vol. 4, pp.
237–264, 1996.

[47]

M. Schmich and B. Vexler, “Adaptivity with dynamic meshes for space-time finite element discretizations of parabolic equations,” *SIAM J. Sci. Comput.*, vol. 30, no.
1, pp. 369–393, 2008.

[48]

R. Hartmann and P. Houston, “Goal-oriented a posteriori error estimation for multiple target functionals , BOOKTITLE = Hyperbolic problems: theory, numerics,
applications,” Springer, Berlin, 2003 , MRCLASS = {35A15 (49M15 65M15 65M60)}, MR, pp. 579–588.

[49]

R. Hartmann, “Multitarget error estimation and adaptivity in aerodynamic flow simulations,” *SIAM J. Sci. Comput. , FJOURNAL = SIAM Journal on Scientific
Computing*, vol. 31, no. 1, pp. 708–731, 2008, doi: 10.1137/070710962.

[50]

F. E. M. (FEM). D. Pardo" keywords = "Multiple goals, “Multigoal-oriented adaptivity for hp-finite element methods,” *Procedia Comput. Sci.", optjournal = "Procedia
Computer Science*, vol. 1, no. 1, pp. 1953–1961, 2010", OPTnote = "ICCS 2010, doi: https://doi.org/10.1016/j.procs.2010.04.219.

[51]

B. Endtmayer and T. Wick, “A Partition-of-Unity Dual-Weighted Residual Approach for
Multi-Objective Goal Functional Error Estimation Applied to Elliptic Problems,” *Comput. Methods Appl. Math. , FJOURNAL =
Computational Methods in Applied Mathematics*, vol. 17, no. 4, pp. 575–599, 2017, [Online]. Available: https://doi.org/10.1515/cmam-2017-0001.

[52]

B. Endtmayer, U. Langer, and T. Wick, “Multigoal-oriented error estimates for non-linear problems,” *J. Numer. Math. , FJOURNAL = Journal of Numerical
Mathematics*, vol. 27, no. 4, pp. 215–236, 2019, [Online]. Available: https://doi.org/10.1515/jnma-2018-0038.

[53]

R. H. Nochetto, A. Veeser, and M. Verani, “A safeguarded dual weighted residual method,” *IMA J. Numer. Anal. , FJOURNAL = IMA Journal of Numerical
Analysis*, vol. 29, no. 1, pp. 126–140, 2009, doi: 10.1093/imanum/drm026.

[54]

M. Ainsworth and R. Rankin, “Guaranteed computable bounds on quantities of interest in finite element computations,” *Int. J. Numer. Methods Eng. , FJOURNAL =
International Journal for Numerical Methods in Engineering*, vol. 89, no. 13, pp. 1605–1634, 2012, doi: 10.1002/nme.3276.

[55]

M. Feischl, D. Praetorius, and K. G. van der Zee, “An,” *SIAM J. Numer. Anal. , FJOURNAL = SIAM Journal on Numerical Analysis*, vol. 54, no. 3, pp.
1423–1448, 2016, doi: 10.1137/15M1021982.

[56]

K. Kergrene, S. Prudhomme, L. Chamoin, and M. Laforest, “A new goal-oriented formulation of the finite element method,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL =
Computer Methods in Applied Mechanics and Engineering*, vol. 327, no. 3725770, pp. 256–276, 2017, [Online]. Available: https://doi.org/10.1016/j.cma.2017.09.018.

[57]

M. S. Mommer and R. Stevenson, “A goal-oriented adaptive finite element method with convergence rates,” *SIAM J. Numer. Anal. , FJOURNAL = SIAM Journal on Numerical
Analysis*, vol. 47, no. 2, pp. 861–886, 2009, doi: 10.1137/060675666.

[58]

R. Becker, E. Estecahandy, and D. Trujillo, “Weighted marking for goal-oriented adaptive finite element methods,” *SIAM J. Numer. Anal. , issue_date = November
2011*, vol. 49, no. 6, pp. 2451–2469, 2011, doi: 10.1137/100794298 , acmid = {2340490}.

[59]

P. Ladevèze, F. Pled, and L. Chamoin, “New bounding techniques for goal-oriented error estimation applied to linear problems,” *Int. J. Numer. Methods Eng. , FJOURNAL =
International Journal for Numerical Methods in Engineering*, vol. 93, no. 13, pp. 1345–1380, 2013, [Online]. Available: https://doi.org/10.1002/nme.4423.

[60]

I. Mozolevski and S. Prudhomme, “Goal-oriented error estimation based on equilibrated-flux reconstruction for finite element approximations of elliptic problems,”
*Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 288, no. 3327021 , MRREVIEWER = Lucas Jódar, pp. 127–145, 2015, [Online]. Available: https://doi.org/10.1016/j.cma.2014.09.025.

[61]

D. Meidner and T. Richter, “Goal-oriented error estimation for the fractional step theta scheme,” *Comput. Methods Appl. Math. , FJOURNAL = Computational Methods in
Applied Mathematics*, vol. 14, no. 2, pp. 203–230, 2014, doi: 10.1515/cmam-2014-0002.

[62]

T. Richter and T. Wick, “Variational localizations of the dual weighted residual estimator,” *J. Comput. Appl. Math. , FJOURNAL = Journal of Computational and
Applied Mathematics*, vol. 279, no. 3293320, pp. 192–208, 2015, doi: 10.1016/j.cam.2014.11.008.

[63]

B. Endtmayer, U. Langer, and A. Schafelner, “Goal-oriented adaptive space-time finite element methods for regularized parabolic \(p\)-Laplace problems,” *arXiv preprint arXiv:2306.07167*, 2023.

[64]

J. P. Thiele and T. Wick, “Numerical modeling and open-source implementation of variational partition-of-unity localizations of space-time dual-weighted residual estimators for
parabolic problems,” *J. Sci. Comput.*, vol. 99, no. 25 , OPT, pp. 237–264, 2024.

[65]

D. Braess, V. Pillwein, and J. Schöberl, “Equilibrated residual error estimates are \(p\)-robust,” *Comput. Methods Appl.
Mech. Eng. , FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 198, no. 13–14, pp. 1189–1197, 2009, [Online]. Available: https://doi.org/10.1016/j.cma.2008.12.010.

[66]

B. N. Granzow, D. T. Seidl, and S. D. Bond, “Linearization errors in discrete goal-oriented error estimation,” *Comput. Methods Appl. Mech. Eng. ,
optjournal=Computer Methods in Applied Mechanics and Engineering*, vol. 416, p. 116364, 2023.

[67]

C. Carstensen, M. Feischl, M. Page, and D. Praetorius, “Axioms of adaptivity,” *Comput. Math. Appl. , FJOURNAL = Computers & Mathematics with
Applications. An International Journal*, vol. 67, no. 6, pp. 1195–1253, 2014, doi: 10.1016/j.camwa.2013.12.003.

[68]

M. Holst, S. Pollock, and Y. Zhu, “Convergence of goal-oriented adaptive finite element methods for semilinear problems,” *Comput. Vis. Sci. , FJOURNAL = Computing
and Visualization in Science*, vol. 17, no. 1, pp. 43–63, 2015, doi: 10.1007/s00791-015-0243-1.

[69]

M. Holst and S. Pollock, “Convergence of goal-oriented adaptive finite element methods for nonsymmetric problems,” *Numer. Methods Partia. Diff. Equ. , FJOURNAL =
Numerical Methods for Partial Differential Equations. An International Journal*, vol. 32, no. 2, pp. 479–509, 2016, doi: 10.1002/num.22002.

[70]

R. Becker, M. Innerberger, and D. Praetorius, “Optimal convergence rates for goal-oriented FEM with quadratic goal functional,” *Comput. Methods Appl. Math.
, FJOURNAL = Computational Methods in Applied Mathematics*, vol. 21, no. 2, pp. 267–288, 2021, doi: 10.1515/cmam-2020-0044.

[71]

R. Becker, M. Brunner, M. Innerberger, J. M. Melenk, and D. Praetorius, “Rate-optimal goal-oriented adaptive FEM for semilinear elliptic PDEs,” *Comput. Math. Appl. ,
optjournal = Computers & Mathematics with Applications*, vol. 118, pp. 18–35, 2022.

[72]

R. Becker, G. Gantner, M. Innerberger, and opt Praetorius Dirk, “Goal-oriented adaptive finite element methods with optimal computational complexity,” *Numerische
Mathematik , journal=Numer. Math.*, vol. 153, no. 1, pp. 111–140, 2023.

[73]

B. Endtmayer, U. Langer, and T. Wick, “Two-Side a Posteriori Error Estimates for the Dual-Weighted
Residual Method,” *SIAM J. Sci. Comput. , FJOURNAL = SIAM Journal on Scientific Computing*, vol. 42, no. 1, pp. A371–A394, 2020, doi: 10.1137/18M1227275.

[74]

B. Endtmayer, U. Langer, and T. Wick, “Reliability and efficiency of DWR-type a posteriori error estimates with smart sensitivity weight recovering,” *Comput. Methods
Appl. Math. , FJOURNAL = Computational Methods in Applied Mathematics*, vol. 21, no. 2, pp. 351–371, 2021.

[75]

W. Bangerth, R. Hartmann, and G. Kanschat, “Deal.II – a general purpose object oriented finite element library,” *ACM Trans. Math. Softw.*, vol. 33, no. 4, pp.
24/1–24/27, 2007.

[76]

G. Alzetta *et al.*, “The *J. Numer. Math. , FJOURNAL = Journal of Numerical
Mathematics*, vol. 26, no. 4, pp. 173–183, 2018, doi: 10.1515/jnma-2018-0054.

`deal.`

`II`

library, version 9.0,” [77]

D. Arndt *et al.*, “The deal.II finite element library: Design, features, and insights,” *Comput. Math. Appl. , optjournal =
Computers & Mathematics with Applications*, vol. 81, pp. 407–422, 2021, doi: 10.1016/j.camwa.2020.02.022.

[78]

J. M. Melenk and I. Babuška, “The partition of unity finite element method: Basic theory and applications,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer
Methods in Applied Mechanics and Engineering*, vol. 139, no. 1–4, pp. 289–314, 1996, doi: 10.1016/S0045-7825(96)01087-0.

[79]

I. Babuška and J. M. Melenk, “The partition of unity method,” *Int. J. Numer. Methods Eng. , FJOURNAL = International Journal for Numerical Methods in
Engineering*, vol. 40, no. 4, pp. 727–758, 1997, doi: 10.1002/(SICI)1097-0207(19970228)40:4<727::AID-NME86>3.3.CO;2-E.

[80]

J. P. Thiele and T. Wick, “Jpthiele/pu-dwr-diffusion: v1.0.0.” Zenodo , version = v1.0.0, Feb. 2024, doi: 10.5281/zenodo.10641119.

[81]

J. P. Thiele, “Jpthiele/pu-dwr-combustion: v1.0.0.” Zenodo , version = v1.0.0, Feb. 2024, doi: 10.5281/zenodo.10641104.

[82]

J. Roth, J. P. Thiele, U. K?cher, and T. Wick, “Tensor-Product Space-Time Goal-Oriented Error
Control and Adaptivity With Partition-of-Unity Dual-Weighted Residuals for Nonstationary Flow Problems,”
*Comput. Methods Appl. Math.*, vol. 24, no. 1, pp. 185–214, 2024 , lastchecked = {2024-01-15}, doi: doi:10.1515/cmam-2022-0200.

[83]

R. Becker and R. Rannacher, “Weighted a posteriori error control in FE methods , booktitle = Lecture ENUMATH-95, Paris, Sept. 18-22, 1995, in:
Proc. ENUMATH-97, Heidelberg, Sept. 28 - Oct.3, 1997 (H.G. Bock, et al., eds),” 1998.

[84]

A. Rademacher, “Adaptive finite element methods for nonlinear hyperbolic problems of second order , school = Technische Universität
Dortmund,” PhD thesis, 2009.

[85]

W. Bangerth, M. Geiger, and R. Rannacher, “Adaptive Galerkin finite element methods for the wave equation,” *Comput. Methods Appl. Math.*, vol. 10, pp.
3–48, 2010.

[86]

M. Soszyńska and T. Richter, “Adaptive time-step control for a monolithic multirate scheme coupling the heat and wave equation,” *BIT Numerical Mathematics*, 2021,
doi: 10.1007/s10543-021-00854-3 , eprint = {2007.05372}, archiveprefix = {arxiv}.

[87]

R. Rannacher and J. Vihharev, “Adaptive finite element analysis of nonlinear problems: Balancing of discretization and iteration errors,” *J. Numer. Math. , FJOURNAL =
Journal of Numerical Mathematics*, vol. 21, no. 1, pp. 23–61, 2013, doi: 10.1515/jnum-2013-0002.

[88]

R. Rannacher and F.-T. Suttmeier, “A posteriori error control in finite element methods via duality techniques: Application to perfect plasticity,” *Comput. Mech. ,
FJOURNAL = Computational Mechanics*, vol. 21, no. 2, pp. 123–133, 1998, doi: 10.1007/s004660050288.

[89]

R. Rannacher and F.-T. Suttmeier, “A posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity,” *Comput. Methods Appl. Mech. Eng. ,
optjournal = "Computer Methods in Applied Mechanics and Engineering*, vol. 176, no. 1–4, pp. 333–361, 1999", note =, doi: http://dx.doi.org/10.1016/S0045-7825(98)00344-2.

[90]

R. Rannacher and F.-T. Suttmeier, “Error estimation and adaptive mesh design for FE models in elasto-plasticity , booktitle = Error-Controlled
Adaptive FEMs in Solid Mechanics,” John Wiley, 2000 , editor = {E. Stein}.

[91]

C. Mehlmann and T. Richter, “A goal oriented error estimator and mesh adaptivity for sea ice simulations,” *Ocean Model. , optjournal = Ocean Modelling*,
vol. 154, p. 101684, 2020, doi: https://doi.org/10.1016/j.ocemod.2020.101684.

[92]

R. Becker, V. Heuveline, and R. Rannacher, “An optimal control approach to adaptivity in computational fluid mechanics,” *Int. J. Numer. Methods Fluids , fjournal =
Int. J. Numer. Methods Fluids, OPTfjournal = International Journal for Numerical Methods in Fluids*, vol. 40, no. 1–2, pp. 105–120, 2002.

[93]

M. Braack and T. Richter, “Solutions of 3D navier–stokes benchmark problems with adaptive finite elements,” *Comput. Fluids , optjournal=Computers &
Fluids*, vol. 35, pp. 372–392, May 2006, doi: 10.1016/j.compfluid.2005.02.001.

[94]

M. Besier, “Adaptive finite element methods for computing nonstationary incompressible flows , school = University of Heidelberg,” PhD thesis,
2009.

[95]

M. Besier and R. Rannacher, “Goal-oriented space-time adaptivity in the finite element galerkin method for the computation of nonstationary incompressible flow,” *Int. J.
Num. Meth. Fluids*, vol. 70, pp. 1139–1166, 2012.

[96]

B. Endtmayer, U. Langer, J. P. Thiele, and booktitle=Numerical. M. and A. A. E. 2019. Wick Thomas, “Hierarchical DWR error estimates for the navier-stokes equations: \(h\) and \(p\) enrichment,” Springer, 2021, pp. 363–372.

[97]

D. Kuzmin and S. Korotov, “Goal-oriented a posteriori error estimates for transport problems,” *Math. Comput. Simulation , FJOURNAL = Mathematics and Computers in
Simulation*, vol. 80, no. 8, pp. 1674–1683, 2010, doi: 10.1016/j.matcom.2009.03.008.

[98]

M. Bause, M. P. Bruchhäuser, and U. Köcher, “Flexible goal-oriented adaptivity for higher-order space–time discretizations of transport problems with coupled flow,”
*Comput. Math. Appl. , optjournal = Computers & Mathematics with Applications*, vol. 91, pp. 17–35, 2021.

[99]

M. P. Bruchhäuser, “Goal-oriented space-time adaptivity for a multirate approach to coupled flow and transport,” PhD thesis, 2022 , school={Helmut Schmidt University
Hamburg}.

[100]

M. P. Bruchhäuser and M. Bause, “A cost-efficient space-time adaptive algorithm for coupled flow and transport,” *Comput. Methods Appl. Math. ,
optjournal=Computational Methods in Applied Mathematics*, vol. 23, no. 4, pp. 849–875, 2023.

[101]

S. Beuchler, B. Endtmayer, J. Lankeit, and T. Wick, “Multigoal-oriented a posteriori error control for heated material processing using a generalized Boussinesq
model,” *Comptes Rendus. Mécanique*, vol. Special, 2023, doi: 10.5802/crmeca.160 , language = {en}, note = {Online first}.

[102]

S. Beuchler, A. Demircan, B. Endtmayer, U. Morgner, and T. Wick, “Mathematical modeling and numerical multigoal-oriented a posteriori error control and adaptivity for a
stationary, nonlinear, coupled flow temperature model with temperature dependent density.” 2024 , note = {submitted}.

[103]

B. Endtmayer, A. Demircan, D. Perevoznik, U. Morgner, S. Beuchler, and T. Wick, “Adaptive finite element simulations of laser-heated material flow using a boussinesq
model,” *Proc. Appl. Math. Mech.*, vol. 23, no. 1, p. e202200219, 2023, doi: https://doi.org/10.1002/pamm.202200219.

[104]

R. Becker, M. Braack, and R. Rannacher, “On error control for reactive flow problems , booktitle = Scientific Computing in Chemical Engineering II,
Computational Fluid Dynamics, Reaction Engineering, and Molecular Properties,” vol. 1, Springer Berlin, 1999, pp. 320–327.

[105]

T. Richter, “Parallel multigrid method for adaptive finite elements with application to 3D flow problems , school = University of Heidelberg,”
PhD thesis, 2005.

[106]

M. Braack and book T. Richter, “Reactive flows, diffusion and
transport , title = Mesh and model adaptivity for flow problems,” Springer Berlin Heidelberg, 2006, pp. 47–75.

[107]

A. Bespalov, D. Praetorius, L. Rocchi, and M. Ruggeri, “Goal-oriented error estimation and adaptivity for elliptic PDEs with parametric or uncertain inputs,” *Comput.
Methods Appl. Mech. Eng. , optjournal=Computer Methods in Applied Mechanics and Engineering*, vol. 345, pp. 951–982, 2019.

[108]

P. Ingelstr?m and keywords =. M. equations,. A. G. error estimation,. H. basis functions,. W. simulation A. Bondeson, “Goal-oriented error estimation and h-adaptivity for
maxwell?s equations,” *Comput. Methods Appl. Mech. Eng. , optjournal = Computer Methods in Applied Mechanics and Engineering*, vol. 192, no. 22, pp. 2597–2616, 2003, doi: https://doi.org/10.1016/S0045-7825(03)00295-0.

[109]

V. Heuveline and R. Rannacher, “A posteriori error control for finite element approximations of elliptic eigenvalue problems,” *Adv. Comput. Math. , optjournal="Advances
in Computational Mathematics*, vol. 15, no. 1, pp. 107–138, 2001, doi: 10.1023/A:1014291224961.

[110]

B. Cockburn and S. Xia, “An adjoint-based super-convergent galerkin approximation of eigenvalues,” *J. Comput. Phys. , optjournal=Journal of Computational
Physics*, vol. 449, p. 110816, 2022.

[111]

J. T. Oden and S. Prudhomme, “Estimation of modeling error in computational mechanics,” *J. Comput. Phys. , optjournal = "Journal of Computational Physics*, vol.
182, no. 2, pp. 496–515, 2002.

[112]

A. Rademacher, “Mesh and model adaptivity for frictional contact problems , opt,” *Numerische Mathematik , journal =Numer. Math*, vol. 142 , OPT, no. 2,
pp. 465–523, 2019.

[113]

S. Weißer and T. Wick, “The Dual-Weighted Residual Estimator Realized on Polygonal
Meshes,” *Comput. Methods Appl. Math. , FJOURNAL = Computational Methods in Applied Mathematics*, vol. 18, no. 4, pp. 753–776, 2018, doi: 10.1515/cmam-2017-0046.

[114]

G. Mallik, M. Vohralik, and A. posteriori error estimate Soleiman Yousef" keywords = "Quantity of interest, “Goal-oriented a posteriori error estimation for conforming and
nonconforming approximations with inexact solvers,” *J. Comput. Appl. Math.*, vol. 366, p. 112367, 2020, doi: https://doi.org/10.1016/j.cam.2019.112367.

[115]

T. Richter, “A posteriori error estimation and anisotropy detection with the dual-weighted residual method,” *Int. J. Numer. Methods Fluids , FJOURNAL =
International Journal for Numerical Methods in Fluids*, vol. 62, no. 1, pp. 90–118, 2010, doi: 10.1002/fld.2016.

[116]

E. T. Chung, W. T. Leung, and S. Pollock, “Goal-oriented adaptivity for GMsFEM,” *J. Comput. Appl. Math. , optjournal=Journal of Computational and Applied
Mathematics*, vol. 296, pp. 625–637, 2016.

[117]

M. Maier and R. Rannacher, “Duality-based adaptivity in finite element discretization of heterogeneous multiscale problems.” *J. Numer. Math. , OPTFJOURNAL
=Journal of Numerical Mathematics*, vol. 24, no. 3, pp. 167–187, 2016.

[118]

Matthias. Maier and Rolf. Rannacher, “A duality-based optimization approach for model adaptivity in heterogeneous multiscale problems,” *Multiscale Model. Simul.*,
vol. 16, no. 1, pp. 412–428, 2018, doi: 10.1137/16M1105670.

[119]

L. Lautsch and T. Richter, “Error estimation and adaptivity for differential equations with multiple scales in time,” *Comput. Methods Appl. Math. , OPTFJOURNAL =
Computational Methods in Applied Mathemacics*, 2021, doi: 10.1515/cmam-2021-0030 , eprint = {2006.05420}.

[120]

D. Dominguez, L. Lautsch, and T. Richter, “A variational approach for temporal multiscale problems and its application to adaptivity and optimization , eid = e202300193,” *Proc. Appl. Math. Mech.*, 2023, doi: https://doi.org/10.1002/pamm.202300193.

[121]

T. Grätsch and K.-J. Bathe, “Goal-oriented error estimation in the analysis of fluid flows with structural interactions,” *Comp. Methods Appl. Mech. Engrg.*, vol.
195, pp. 5673–5684, 2006.

[122]

T. Richter, “Goal-oriented error estimation for fluid-structure interaction problems,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer Methods in Applied
Mechanics and Engineering*, vol. 223/224, no. 2917479 , MRREVIEWER = Xiaoguang Allan Zhong, pp. 28–42, 2012, doi: 10.1016/j.cma.2012.02.014.

[123]

“An Eulerian approach to fluid-structure interaction and goal-oriented mesh adaption,” *Int. J. Numer. Methods Fluids*, vol. 51, pp. 1017–1039,
2006.

[124]

T. Dunne, “Adaptive finite element approximation of fluid-structure interaction based on Eulerian and arbitrary Lagrangian-Eulerian
variational formulations , school = University of Heidelberg,” PhD thesis, 2007.

[125]

T. Dunne, T. Richter, and R. Rannacher, “Numerical simulation of fluid-structure interaction based on monolithic variational formulations,” Springer, 2010 , series =
{Comtemporary Challenges in Mathematical Fluid Mechanics}, address = {World Scientific, Singapore}, pp. 1–75.

[126]

T. Wick, “Adaptive Finite Element Simulation of
Fluid-Structure Interaction with Application to Heart-Valve Dynamics , school = University of Heidelberg,” PhD thesis, 2011 ,
OPT.

[127]

T. Wick, “Goal-oriented mesh adaptivity for fluid-structure interaction with application to heart-valve settings,” *Arch. Mech. Eng. , optjournal = Archive of
Mechanical Engineering *, vol. LIX, no. 1, pp. 73–99, 2012, doi: 10.2478/v10180-012-0005-2.

[128]

P. W. Fick, E. H. van Brummelen, and K. G. van der Zee, “On the adjoint-consistent formulation of interface conditions in goal-oriented error estimation and adaptivity for
fluid-structure interaction,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 199, no. 49–52, pp. 3369–3385, 2010, doi: 10.1016/j.cma.2010.07.009.

[129]

K. G. van der Zee, E. H. van Brummelen, I. Akkerman, and R. de Borst, “Goal-oriented error estimation and adaptivity for fluid-structure interaction using exact linearized
adjoints,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 200, no. 37–40, pp. 2738–2757, 2011, doi: 10.1016/j.cma.2010.12.010.

[130]

T. Richter, *Fluid-structure interactions. Models, analysis and finite elements*, vol. 118 , series = Lecture notes in computational science and engineering. Springer,
2017.

[131]

L. Failer and T. Wick, “Adaptive time-step control for nonlinear fluid-structure interaction,” *J. Comp. Phys.*, vol. 366, pp. 448–477, 2018.

[132]

K. Ahuja, B. Endtmayer, M. C. Steinbach, and f. Wick Thomas, “Multigoal-oriented error estimation and mesh adaptivity for fluid–structure interaction,” *Journal of
Computational and Applied Mathematics , journal=J. Comput. Appl. Math.*, vol. 412, p. 114315, 2022.

[133]

T. Wick, “Goal functional evaluations for phase-field fracture using PU-based DWR mesh adaptivity,” *Comput. Mech. , FJOURNAL =
Computational Mechanics*, vol. 57, no. 6, pp. 1017–1035, 2016, doi: 10.1007/s00466-016-1275-1.

[134]

[135]

T. Wick, “Dual-weighted residual a posteriori error estimates for a penalized phase-field slit discontinuity problem,” *Comput. Methods Appl. Math. , optjournal =
Computational Methods in Applied Mathematics*, vol. 21, no. 3, pp. 693–707, 2021 , lastchecked = {2024-03-25}, doi: doi:10.1515/cmam-2020-0038.

[136]

R. Becker, H. Kapp, and R. Rannacher, “Adaptive finite element methods for optimal control of partial differential equations: Basic concepts,” *SIAM J. Optim.
Control*, vol. 39, pp. 113–132, 2000.

[137]

R. Becker, M. Braack, D. Meidner, R. Rannacher, and B. Vexler, “Adaptive finite element methods for
PDE-constrained optimal control problems , BOOKTITLE = Reactive flows, diffusion and transport,” Springer, Berlin, 2007 , MRCLASS = {49M25 (35J60 49K20 65N30)}, MR, pp. 177–205.

[138]

D. Meidner and B. Vexler, “Adaptive space-time finite element methods for parabolic optimization problems,” *SIAM J. Control Optim. , optjournal = SIAM Journal on
Control and Optimization*, vol. 46, no. 1, pp. 116–142, 2007.

[139]

D. Meidner, “Adaptive space-time finite element methods for optimization problems governed by nonlinear parabolic systems , school = University of
Heidelberg,” PhD thesis, 2008.

[140]

M. Hintermüller and R. H. W. Hoppe, “Goal-oriented adaptivity in control constrained optimal control of partial differential equations , opt,” *SIAM Journal on Control
and Optimization , journal = SIAM J. Control Opt.*, vol. 47, no. 4, pp. 1721–1743, 2008, doi: 10.1137/070683891.

[141]

B. Vexler and W. Wollner, “Adaptive finite elements for elliptic optimization problems with control constraints , opt,” *SIAM Journal on Control and Optimization ,
journal = SIAM J. Control Opt.*, vol. 47, no. 1, pp. 509–534, 2008, doi: 10.1137/070683416.

[142]

M. Hintermüller and R. H. W. Hoppe, “Goal-oriented adaptivity in pointwise state constrained optimal control of partial differential equations , opt,” *SIAM Journal on
Control and Optimization , journal = SIAM J. Control Opt.*, vol. 48, no. 8, pp. 5468–5487, 2010, doi: 10.1137/090761823.

[143]

R. Rannacher and B. Vexler, “Adaptive finite element discretization in PDE-based optimization,” *GAMM-Mitteilungen*, vol. 33, no. 2, pp. 177–193, keywords =
PDE–based optimization, finite element method, goal–oriented adaptivity, 2010, doi: 10.1002/gamm.201010014.

[144]

M. Hintermüller, M. Hinze, C. Kahle, and T. Keil, “A goal-oriented dual-weighted adaptive finite element approach for the optimal control of a nonsmooth
Cahn-Hilliard-Navier-Stokes system,” *Optim. Eng. , FJOURNAL = Optimization and Engineering. International Multidisciplinary Journal to Promote Optimization Theory & Applications in
Engineering Sciences*, vol. 19, no. 3, pp. 629–662, 2018, doi: 10.1007/s11081-018-9393-6.

[145]

B. Endtmayer, U. Langer, I. Neitzel, T. Wick, and W. Wollner, “Multigoal-oriented optimal control problems with nonlinear PDE constraints,” *Comput. Math.
Appl.*, vol. 79, no. 10, pp. 3001–3026, 2020, doi: 10.1016/j.camwa.2020.01.005.

[146]

H. Blum and F.-T. Suttmeier, “An adaptive finite element discretisation for a simplified signorini problem,” *Calcolo*, vol. 37, no. 2, pp. 65–77, 1999.

[147]

H. Blum and F.-T. Suttmeier, “Weighted error estimates for finite element solutions of variational inequalities,” *Computing*, vol. 65, no. 2, pp. 119–134,
language=English, 2000, doi: 10.1007/s006070070015.

[148]

H. Blum, A. Schröder, and booktitle=Numerical. M. and A. A. Suttmeier Franz-Theo, “A posteriori estimates for FE-solutions of variational inequalities,” Springer, 2003, pp.
669–680.

[149]

A. Schröder and A. Rademacher, “Goal-oriented error control in adaptive mixed FEM for Signorini’s problem,” *Comput. Methods Appl. Mech. Eng. ,
FJOURNAL = Computer Methods in Applied Mechanics and Engineering*, vol. 200, no. 1–4, pp. 345–355, 2011, doi: 10.1016/j.cma.2010.08.015.

[150]

A. Rademacher and A. Schröder, “Dual weighted residual error control for frictional contact problems,” *Comput. Methods Appl. Math. , FJOURNAL = Computational
Methods in Applied Mathematics*, vol. 15, no. 3, pp. 391–413, 2015.

[151]

[152]

P. Stolfo, A. Rademacher, and A. Schröder, “Dual weighted residual error estimation for the finite cell method,” *J. Numer. Math.*, vol. 27, no. 2, 2019.

[153]

P. Di Stolfo and O. J. and W. P. Schröder Andreas", “Error control and adaptivity for the finite cell method",
bookTitle="non-standard discretisation methods in solid mechanics,” Springer International Publishing", address="Cham, 2022, pp. 377–403.

[154]

D. Meidner, R. Rannacher, and J. Vihharev, “Goal-oriented error control of the iterative solution of finite element equations,” *J. Numer. Math. , FJOURNAL =
Journal of Numerical Mathematics*, vol. 17, no. 2, pp. 143–172, 2009, doi: 10.1515/JNUM.2009.009.

[155]

R. Rannacher, A. Westenberger, and W. Wollner, “Adaptive finite element solution of eigenvalue problems: Balancing of discretization and iteration error,” *J. Numer.
Math. , FJOURNAL = Journal of Numerical Mathematics*, vol. 18, no. 4, pp. 303–327, 2010, [Online]. Available: https://doi.org/10.1515/JNUM.2010.015.

[156]

V. Dolejšı́, O. Bartoš, and f. Roskovec Filip, “Goal-oriented mesh adaptation method for nonlinear problems including algebraic errors,” *Comput. Math. Appl. , optjournal
= Computers & Mathematics with Applications, journal = Comput. Math. Appl.*, vol. 93, pp. 178–198, 2021.

[157]

V. Dolejšı́ and S. Congreve, “Goal-oriented error analysis of iterative galerkin discretizations for nonlinear problems including linearization and algebraic errors,” *J.
Comput. Appl. Math. , optjournal=Journal of Computational and Applied Mathematics*, vol. 427, p. 115134, 2023.

[158]

Ch. Goll, R. Rannacher, and W. Wollner, “The damped Crank-Nicolson time-marching scheme for the adaptive solution of the
Black-Scholes equation,” *J. Comput. Finance*, vol. 18, no. 4, pp. 1–37, 2015.

[159]

M. Meyer and H. G. Matthies, “Efficient model reduction in non-linear dynamics using the Karhunen-Loève expansion and dual-weighted-residual
methods,” *Comput. Mech.*, vol. 31, no. 1, pp. 179–191, 2003, doi: 10.1007/s00466-002-0404-1.

[160]

H. Fischer, J. Roth, L. Chamoin, A. Fau, M. F. Wheeler, and T. Wick, “Adaptive space-time model order reduction with dual-weighted residual (MORe DWR) error control for
poroelasticity,” *for advanced modeling and simulations in engineering sciences*. 2024 , note = {accepted for publication in}.

[161]

H. Fischer, J. Roth, T. Wick, L. Chamoin, and keywords =. T. space reduced modeling,. D. residual method,. G. error control,. I. proper orthogonal decomposition Amelie Fau, “MORe
DWR: Space-time goal-oriented error control for incremental POD-based ROM for time-averaged goal functionals,” *J. Comput. Phys. , optjournal = Journal of Computational Physics*, vol. 504, p. 112863, 2024, doi: https://doi.org/10.1016/j.jcp.2024.112863.

[162]

K. Kergrene, L. Chamoin, M. Laforest, and S. Prudhomme, “On a goal-oriented version of the proper generalized decomposition method,” *J. Sci. Comput. , FJOURNAL =
Journal of Scientific Computing*, vol. 81, no. 1, pp. 92–111, 2019, doi: 10.1007/s10915-019-00918-1.

[163]

L. Chamoin and F. Legoll, “Goal-oriented error estimation and adaptivity in MsFEM computations,” *Comput. Mech. , FJOURNAL =
Computational Mechanics*, vol. 67, no. 4, pp. 1201–1228, 2021, doi: 10.1007/s00466-021-01990-x.

[164]

I. Brevis, I. Muga, and keywords =. G. finite elements,. M. acceleration,. R. minimization,. P. method,. W. inner D. algorithms Kristoffer G. van der Zee, “A machine-learning
minimal-residual (ML-MRes) framework for goal-oriented finite element discretizations,” *Comput. Math. Appl. , optjournal = Computers & Mathematics with Applications*, vol. 95, pp. 186–199, 2021 , note = {Recent
Advances in Least-Squares and Discontinuous Petrov?Galerkin Finite Element Methods}, doi: https://doi.org/10.1016/j.camwa.2020.08.012.

[165]

P. Minakowski and T. Richter, “A priori and a posteriori error estimates for the deep Ritz method applied to the Laplace and Stokes
problem,” *J. Comput. Appl. Math. , optjournal = Journal of Computational and Applied Mathematics*, vol. 421, p. 114845, 2023, doi: 10.1016/j.cam.2022.114845 , eprint = {2107.11035}.

[166]

J. Roth, M. Schröder, and T. Wick, “Neural network guided adjoint computations in dual weighted residual error estimation,” *SN Applied Sciences*, vol. 4, no. 2, p.
62, 2022, doi: 10.1007/s42452-022-04938-9.

[167]

J. T. Oden, “Adaptive multiscale predictive modelling,” *Acta Numer.*, vol. 27, pp. 353–450, 2018, doi: 10.1017/S096249291800003X.

[168]

C. Goll, T. Wick, and W. Wollner, “DOpElib : Differential equations and optimization environment; A goal oriented software library for solving PDEs and
optimization problems with PDEs,” *Archive of Numerical Software*, vol. 5, no. 2, pp. 1–14, 2017, doi: 10.11588/ans.2017.2.11815.

[169]

U. Köcher, M. P. Bruchhäuser, and M. Bause, “Efficient and scalable data structures and algorithms for goal-oriented adaptivity of space–time FEM codes,”
*SoftwareX*, vol. 10, p. 100239, 2019.

[170]

J. P. Thiele, “Error-controlled space-time finite elements, algorithms, and implementations for nonstationary problems , school = Leibniz University Hannover,”
PhD thesis, 2024.

[171]

O. E. Zeidler ALTeditor =, *Nonlinear functional analysis and its applications II/b: Nonlinear monotone
operators*, vol., OPT. Springer, 1990.

[172]

P. Deuflhard, *Newton methods for nonlinear problems*, vol. 35 , series = Springer Series in Computational Mathematics. Springer Berlin Heidelberg, 2011.

[173]

M. P. Bruchhäuser, K. Schwegler, and M. Bause, “Numerical study of goal-oriented error control for stabilized finite element methods , BOOKTITLE = Advanced
Finite Element Methods with Applications: Selected Papers from the 30th Chemnitz Finite Element Symposium 2017,” Springer International Publishing, 2019, pp. 85–106.

[174]

E. Creusé, S. Nicaise, and Z. Tang, “Goal-oriented error estimation based on equilibrated flux and potential reconstruction for the approximation of elliptic and parabolic
problems,” *Comput. Math. Appl. , optjournal = Computers & Mathematics with Applications*, vol. 146, pp. 323–338, 2023.

[175]

K. Schwegler, M. P. Bruchhäuser, and M. Bause, “Goal-oriented a posteriori error control for nonstationary convection-dominated transport problems,” *ArXiv e-prints ,
archivePrefix = "arXiv", eprint = 1601.06544, primaryClass = "math.NA", keywords = Mathematics - Numerical Analysis, 65M12, 65M60*, 2016, [Online]. Available: http://adsabs.harvard.edu/abs/2016arXiv160106544S , adsnote = {Provided by the SAO/NASA Astrophysics Data System}.

[176]

C. Erath, G. Gantner, and D. Praetorius, “Optimal convergence behavior of adaptive FEM driven by simple (h- h/ 2)-type error estimators,” *Comput. Math. Appl.*, vol.
79, no. 3, pp. 623–642, 2020.

[177]

R. E. Bank and A. Weiser, “Some a posteriori error estimators for elliptic partial differential equations,” *Math. Comp. , FJOURNAL = Mathematics of
Computation*, vol. 44, no. 170, pp. 283–301, 1985, doi: 10.2307/2007953.

[178]

F. A. Bornemann, B. Erdmann, and R. Kornhuber, “A posteriori error estimates for elliptic problems in two and three space dimensions,” *SIAM J. Numer. Anal. , FJOURNAL =
SIAM Journal on Numerical Analysis*, vol. 33, no. 3, pp. 1188–1204, 1996, doi: 10.1137/0733059.

[179]

R. E. Bank and R. K. Smith, “A posteriori error estimates based on hierarchical bases,” *SIAM J. Numer. Anal. , FJOURNAL = SIAM Journal on Numerical
Analysis*, vol. 30, no. 4, pp. 921–935, 1993, doi: 10.1137/0730048.

[180]

A. De Rossi, “Saturation assumption and finite element method for a one-dimensional model,” *RGMIA Research Report Collection*, vol. 5, no. 2, pp. Article 13, 1–6,
2002.

[181]

A. Agouzal, “On the saturation assumption and hierarchical a posteriori error estimator,” *Comput. Methods Appl. Math. , FJOURNAL = Computational Methods in
Applied Mathematics*, vol. 2, no. 2, pp. 125–131, 2002.

[182]

W. Dörfler and R. H. Nochetto, “Small data oscillation implies the saturation assumption,” *Numer. Math. , FJOURNAL = Numerische Mathematik*, vol. 91,
no. 1, pp. 1–12, 2002, doi: 10.1007/s002110100321.

[183]

B. Achchab, S. Achchab, and A. Agouzal, “Some remarks about the hierarchical a posteriori error estimate,” *Numer. Methods Partia. Diff. Equ. , FJOURNAL =
Numerical Methods for Partial Differential Equations. An International Journal*, vol. 20, no. 6, pp. 919–932, 2004, doi: 10.1002/num.20016.

[184]

C. Carstensen, D. Gallistl, and J. Gedicke, “Justification of the saturation assumption,” *Numer. Math. , FJOURNAL = Numerische Mathematik*, vol. 134,
no. 1, pp. 1–25, 2016, doi: 10.1007/s00211-015-0769-7.

[185]

R. E. Bank, A. Parsania, and S. Sauter, “Saturation estimates for \(hp\)-finite element methods,” *Comput. Vis. Sci. ,
FJOURNAL = Computing and Visualization in Science*, vol. 16, no. 5, pp. 195–217, 2013, doi: 10.1007/s00791-015-0234-2.

[186]

S. Ferraz-Leite, C. Ortner, and D. Praetorius, “Convergence of simple adaptive Galerkin schemes based on \(h-h/2\) error
estimators,” *Numer. Math. , FJOURNAL = Numerische Mathematik*, vol. 116, no. 2, pp. 291–316, 2010, doi: 10.1007/s00211-010-0292-9.

[187]

H. Kim and S.-G. Kim, “Saturation assumptions for a 1d convection-diffusion model,” *Korean J. Math. , fjournal = The Korean Journal of Mathematics*,
vol. 22, no. 4, pp. 599–609, keywords = saturatoin assumption, hierarchical basis, convection–diffusion model, 2014, [Online]. Available: http://kkms.org/index.php/kjm/article/view/318.

[188]

B. Endtmayer, “Multi-goal oriented a posteriori error estimates for nonlinear partial differential equations,” PhD thesis, 2020 , school={Johannes Kelper University Linz},
OPTnote={Available at: \url{https://epub.jku.at/obvulihs/download/pdf/5767444}}.

[189]

I. Babuška and W. C. Rheinboldt, “A-posteriori error estimates for the finite element method , opt,” *International Journal for Numerical Methods in Engineering ,
journal = Int. J. Numer. Methods Eng.*, vol. 12, no. 10, pp. 1597–1615, 1978, doi: 10.1002/nme.1620121010.

[190]

J. Nitsche, “Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von
Teilräumen, die keinen Randbedingungen unterworfen sind , NOTE = Collection of articles dedicated to Lothar Collatz on his sixtieth birthday,” *Abh. Math. Sem. Univ. Hamburg , FJOURNAL =
Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg*, vol. 36, no. 341903 , MRREVIEWER = I. Aganovic, pp. 9–15, 1971, [Online]. Available: https://doi.org/10.1007/BF02995904.

[191]

P. Minakowski and T. Richter, “Finite element error estimates on geometrically perturbed domains,” *J. Sci. Comput. , optjournal = Journal of Scientific
Computing*, vol. 84, no. 30, 2020, doi: 10.1007/s10915-020-01285-y.

[192]

M. Crouzeix and Raviart, P.-A., “Conforming and nonconforming finite element methods for solving the stationary stokes equations i,” *R.A.I.R.O.*, vol. 7, no. R3,
pp. 33–75, 1973, doi: 10.1051/m2an/197307R300331.

[193]

[194]

D. A. Di Pietro and A. Ern, *Mathematical aspects of discontinuous Galerkin methods , SERIES = Mathématiques
& Applications (Berlin) {Mathematics & Applications}*, vol. 69. Springer, Heidelberg, 2012, pp. xviii+384, ISBN = 978-3-642-22979-4, MRCLASS = 65–02 (35A35 35F15 35J25 35Q35 65M60 65N30), MR.

[195]

R. Becker and M. Braack, “A finite element pressure gradient stabilization for the Stokes equations based on local projections,” *Calcolo*, vol. 38 ,
__markedentry = {richter:6}, Library = private, Timestamp = 2015.01.13, no. 4, pp. 173–199, 2001.

[196]

D. Meidner and T. Richter, “A posteriori error estimation for the fractional step theta discretization of the incompressible Navier-Stokes
equations,” *Comp. Meth. Appl. Mech. Engrg.*, vol. 288 , OPT, pp. 45–59, 2015 , OPT.

[197]

Q. Chen and M. Gunzburger, “Goal-oriented a posteriori error estimation for finite,” *J. Comput. Appl. Math. , optjournal = Journal of Computational and Applied
Mathematics, volume = 265*, vol. methods, pp. 69–82, 2014, doi: 10.1016/j.cam.2013.10.004.

[198]

G. Grajewski, J. Hron, and S. Turek, “Dual weighted a posteriori error estimation for a new nonconforming linear finite element on quadrilaterals,” *Appl. Numer. Math. ,
optjournal = Applied Numerical Mathematics*, vol. 54, pp. 504–518, 2005.

[199]

W. E and B. Yu, “The Deep Ritz method: A deep learning-based numerical algorithm for solving variational problems,” *Commun. Math. Stat. ,
optjournal = Communications in Mathematics and Statistics*, vol. 6, no. 1, pp. 1–12, Feb. 2018, doi: 10.1007/s40304-018-0127-z , file =
{:Literature/E2018.pdf:PDF}.

[200]

G. Chavent and J. Jaffré, *Mathematical models and finite elements for reservoir simulation: Single phase, multiphase and multicomponent flows through porous media*, vol. 17.
Elsevier, 1986.

[201]

R. W. Lewis and B. A. Schrefler, *The finite element method in the static and dynamic deformation and consolidation of porous media, 2nd edition*. Wiley, 1999.

[202]

O. Coussy, *Poromechanics*. Wiley , address = , isbn = , series =, 2004.

[203]

[204]

[205]

Y. Bazilevs, K. Takizawa, and T. E. Tezduyar, *Computational fluid-structure interaction: Methods and Applications*, vol., series = , address =. Wiley,
2013.

[206]

B. Endtmayer, U. Langer, and T. Wick, “Multiple goal-oriented error estimates applied to 3d non-linear problems,” *Proc. Appl. Math. Mech.*, vol. 18, no. 1, p.
e201800048, 2018, doi: 10.1002/pamm.201800048.

[207]

E. H. van Brummelen, S. Zhuk, and G. J. van Zwieten, “Worst-case multi-objective error estimation and adaptivity,” *Comput. Methods Appl. Mech. Eng. , FJOURNAL =
Computer Methods in Applied Mechanics and Engineering*, vol. 313, no. 3577958, pp. 723–743, 2017, [Online]. Available: https://doi.org/10.1016/j.cma.2016.10.007.

[208]

S. Beuchler, B. Endtmayer, and T. Wick, “Goal oriented error control for stationary incompressible flow coupled to a heat equation,” *Proc. Appl. Math. Mech.*, vol.
21, no. 1, p. e202100151, 2021, doi: https://doi.org/10.1002/pamm.202100151.

[209]

W. Sierpiński, *General topology , SERIES = Mathematical Expositions, No. 7, NOTE = Translated by C. Cecilia Krieger*. University of
Toronto Press, Toronto, 1952, pp. xii+290, MRCLASS = 56.0X, MR.

[210]

M. A. Khamsi and W. A. Kirk, *An introduction to metric spaces and fixed point theory , SERIES = Pure and Applied
Mathematics (New York)*. Wiley-Interscience, New York, 2001, pp. x+302, ISBN = 0-471-41825-0, MRCLASS = 46-01 (46Bxx 47-01 47H10 54-01 54E35 54H25), MR.

[211]

B. Endtmayer, U. Langer, I. Neitzel, T. Wick, and W. Wollner, “Mesh adaptivity and error estimates applied to a regularized \(p\)-Laplacian constrainted optimal control problem for multiple quantities of interest,” *Proc. Appl. Math. Mech.*, vol. 19, no. 1, p. e201900231, 2019, doi: 10.1002/pamm.201900231.

[212]

R. Anderson *et al.*, “MFEM : A modular finite element methods library,” *Comput. Math. Appl. , optjournal = Computers & Mathematics with
Applications, journal = Comput. Math. Appl.*, vol. 81, pp. 42–74, 2021, doi: 10.1016/j.camwa.2020.06.009.

[213]

“MFEM : Modular finite element methods {Software} , howpublished = mfem.org.” doi: 10.11578/dc.20171025.1248.

[214]

D. Arndt *et al.*, “The *J. Numer. Math. , OPTFJOURNAL =Journal of Numerical Mathematics*, vol. 31,
no. 3, pp. 231–246, 2023, doi: 10.1515/jnma-2023-0089.

`deal.II`

library, version 9.5,” [215]

T. A. Davis, “Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method,” *ACM Trans. Math. Softw. , issue_date = June 2004*, vol. 30, no.
2, pp. 196–199, Jun. 2004, doi: 10.1145/992200.992206 , acmid = {992206}.

[217]

B. Heise, “Analysis of a fully discrete finite element method for a nonlinear magnetic field problem,” *SIAM J. Numer. Anal.*, vol. 31, no. 3, pp. 745–759, OPT, 1994
, OPT, doi:, OPTisbn = {}, OPT.

[218]

M. Schäfer, S. Turek, F. Durst, E. Krause, and booktitle=Flow. simulation with high computers I. Rannacher R., “Benchmark computations of laminar flow around a cylinder,”
Springer, 1996, pp. 547–566.

[219]

J. G. Heywood, R. Rannacher, and S. Turek, “Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes
equations,” *Int. J. Numer. Meth. Fl.*, vol. 22, no. 5, pp. 325–352, 1996.

[220]

W. Zulehner, “A short note on inf-sup conditions for the taylor-hood family \(Q_k\)-\(Q_{k-1}\),” *arXiv preprint arXiv:2205.14223*, 2022.

[221]

G. Nabh, “On high order methods for the stationary incompressible navier-stokes equations,” PhD thesis, 1998 , school={Interdisziplin{\"a}res Zentrum f{\"u}r Wiss. Rechnen
der Univ. Heidelberg}.

[222]

O. T. Roubı́c̆ek ALTeditor =, *Nonlinear partial differential equations with
applications*, vol. 153 , OPT. Springer, 2013.

[223]

T. Apel and S. Nicaise, “The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges , f,” *Mathematical Methods in
the Applied Sciences , Journal = Math. Methods Appl. Sci.*, vol. 21, no. 6, pp. 519–549, 1998 , Language = {English}, doi: 10.1002/(SICI)1099-1476(199804)21:6<519::AID-MMA962>3.0.CO;2-R , Keywords =
{65N30,65N50,65N12,35J25}, zbMATH = {1210575}, Zbl = {0911.65107}.

See https://math.nist.gov/amr-benchmark/index.html, and select “L-shaped Domain Homogeneous Boundary Conditions”.↩︎

See https://math.nist.gov/amr-benchmark/index.html and select “Fichera Corner with Vertex and Edge Singularities”, or see [223].↩︎