Keywords: Navier–Stokes Equations \(\cdot\) Goal-Oriented Error Control \(\cdot\) Dual Weighted Residual Method \(\cdot\) Space-Time Adaptivity \(\cdot\) Geometric Multigrid Method \(\cdot\) Local Vanka Smoother \(\cdot\) Nitsche Method

1 Introduction↩︎

The accurate numerical simulation of incompressible viscous fluid flow governed by the Navier–Stokes equations remains challenging, in particular when strong dynamics and fine-scale structures in space and time have to be resolved. The robust and efficient solution of the arising algebraic systems adds a further layer of complexity. Moreover, algorithms that address these challenges have to be implemented in a software framework that is suitable for computations on modern hardware architectures with tens of thousands of compute cores, including hybrid CPU and GPU systems. This places high demands on the underlying data structures. In this work, we design an approach that addresses these requirements. The implementation is based on the deal.II library [1].

Space-time mesh adaptation based on rigorous error representations, or on indicators derived from such representations, has been shown to yield accurate results on computationally feasible grids at moderate numerical cost, even for problems with complex solution structures [2], [3]. Among these approaches, the Dual Weighted Residual method (DWR) [4] offers the possibility of controlling the discretization error in user-defined goal quantities of practical interest. The key idea of this approach is to compute dual weights for the local residuals entering the error representation. This requires the solution of the generally nonlinear primal problem and its linearized dual problem by space-time finite element methods. The DWR method has been successfully applied to numerous problem classes, including fluid flow [5], [6], wave propagation [7], fluid-structure interaction [8], [9], among others. For a comprehensive review with detailed references on various model problems, we refer, for instance, to [4], [8], [10]. Error control in norms associated with the weak formulation constitutes a conceptually different class of techniques for a posteriori error estimation and automatic mesh adaptation. For this, an approach using equilibrated flux reconstructions is proposed in [11].

The DWR concept relies on solving the primal Navier–Stokes system and the associated linearized dual problem by space-time finite element methods. Both, lower- and higher-order approximations, in space and time are of interest. They should maintain stability and preserve as much of the structure of the continuous problem as possible. For the time discretization we use the discontinuous Galerkin method [12], [13]. In space, we employ inf-sup stable finite element pairs with discontinuous pressure. Discontinuous pressure approximations have demonstrated advantages over continuous ones [14]. Moreover, they are advantageous for the construction and application of the slab-wise local Vanka smoother within the geometric multigrid preconditioner. Nonhomogeneous Dirichlet boundary conditions are imposed by Nitsche’s method [15], which leads to a strictly variational formulation rather than enforcing the boundary data through the definition of the function spaces and algebraic systems.

We localize the space-time finite element approximation of the primal and dual DWR problems to space-time slabs \(\Omega\times (t_{n-1},t_n]\), where \(\Omega \subset \mathbb{R}^d\) with \(d \in \{2,3\}\). This yields a time-marching approach. For higher polynomial degrees in time and high spatial resolution, the arising algebraic systems already become large, and the Newton-linearized Jacobian has a complex block structure that complicates the construction of efficient preconditioners. Holistic space-time finite element frameworks, in which all time steps are solved fully coupled and simultaneously, have also been studied yet; cf., e.g., [16]–[19]. Their feasibility continues to remain an open problem if three space variables and a huge number of time steps are involved. The lack of efficient preconditioning techniques for such systems still limits the application of holistic approaches. For this reason, we use a slab-wise approach with a tensor-product mesh structure in this work.

For the slab-wise solution of the Newton-linearized primal problems and of the linear dual problems, we apply a Flexible Generalized Minimum RESidual (FGMRES) method preconditioned by a single \(V\)-cycle of a geometric multigrid (GMG) method. Nowadays, multigrid techniques are typically used as preconditioners for Krylov subspace methods rather than as stand-alone solvers. GMG methods have been studied extensively for the iterative solution of discrete systems arising from stationary partial differential equations. For this we refer to the pioneering work [20]. Owing to their optimal linear complexity, GMG methods are highly efficient and robust for this class of problems [21]. They have also been generalized in several directions, including \(hp\) multigrid methods based on hierarchies of polynomial degrees and meshes. For space-time finite element approximations of the Navier–Stokes equations, block systems with generalized saddle-point structure arise, where the precise block structure depends on the ordering of the slab unknowns [22], [23]. For such systems, the smoothing procedure is decisive for the overall efficiency of the multigrid preconditioner. In the present work, we use a local Vanka smoother, which can be interpreted as a non-overlapping Schwarz method. It is among the most effective smoothers for saddle-point problems, although its cost increases with the polynomial degree and with the number of space-time unknowns [23]. Algebraic MultiGrid (AMG) methods, a further class of multigrid techniques, perform damping of high frequency errors in the solution to linear systems and coarse grid corrections on the algebraic level and do not need any geometric information for building a hierarchy of problems [24], [25]. However, memory usage limits the application of AMG if large scale systems are encountered.

In this work, we study the combination of DWR error estimation, automatic space-time mesh adaptation, and slab-wise geometric multigrid preconditioning for the Navier–Stokes equations. The focus is on goal-oriented adaptivity and its interaction with a GMG-preconditioned Newton–Krylov method. For this purpose, we combine a DWR a posteriori error estimator with adaptive space-time refinement and GMRES iterations preconditioned by slab-wise geometric multigrid. The implementation uses MPI parallelization in deal.II, discontinuous Galerkin time discretization, and inf-sup stable finite element pairs with discontinuous pressure on tensor-product meshes. On benchmark computations, we assess estimator quality, reduction of the goal error relative to total computational work, and solver performance for both primal and adjoint problems. A central question we address is whether the gain from goal-oriented adaptivity justifies the additional cost of the adjoint solve and how this trade-off interacts with the robustness of the solvers.

This work is organized as follows. In Sec. 2, we present the weak form as well as the space-time discretization schemes. In Sec. 3, we derive an a posteriori error representation for the model problem. In Sec. 4, we discuss practical aspects of the implementation as well as the underlying space-time adaptive algorithm. In Sec. 5, we present the results of our numerical examples. Finally, Sec. 6 summarizes the results and outlines directions for future research.

2 Model Problem and Discretization in Space and Time↩︎

In this work, we investigate the accurate and efficient numerical simulation of incompressible viscous flow modeled by the nonstationary Navier–Stokes equations \[\label{eq:1:navier} \begin{array}{rcl @{\,\,}l @{\,\,}l @{\,}l} \partial_t \boldsymbol{v} - \nu \Delta \boldsymbol{v} +\boldsymbol{v}\cdot\nabla\boldsymbol{v} + \nabla p &=& \boldsymbol{f} & \textrm{in} & Q & = \Omega \times I\,,\\[.5ex] \nabla \cdot \boldsymbol{v}&=& 0 & \textrm{in} & Q & = \Omega \times I\,,\\[.5ex] \boldsymbol{v}&=& \boldsymbol{v}_D & \textrm{on} & \Sigma_{D} & = \Gamma_{D} \times I\,,\\[.5ex] \nu\partial_{\boldsymbol{n}}\boldsymbol{v}- p\,\boldsymbol{n} &=& \boldsymbol{0} & \textrm{on} & \Sigma_{\textrm{outflow}} & = \Gamma_{\textrm{outflow}} \times I\,,\\[.5ex] \boldsymbol{v}(0) &=& \boldsymbol{v}_0 & \textrm{on} & \Sigma_0 & = \Omega \times \{ 0 \}\,,\\[.5ex] \end{array}\tag{1}\] In 1 , we denote by \(Q=\Omega \times I\) the space-time domain, where \(\Omega\subset \mathbb{R}^d\), with \(d=2\) or \(d=3\), is a polygonal or polyhedral bounded domain with Lipschitz boundary \(\partial\Omega\) and \(I=(0,T)\), with \(0 < T < \infty\), is a finite time interval. The flow behavior at the boundary \(\partial\Omega = \Gamma_{D} \cup \Gamma_{\textrm{outflow}}\) is modeled by a function \(\boldsymbol{v}_D\) on the segments \(\Gamma_D\coloneq\Gamma_{\textrm{inflow}} \cup \Gamma_{\textrm{wall}}\), that prescribes an inflow profile \(\boldsymbol{v}_{\textrm{in}}\) on \(\Gamma_{\textrm{inflow}}\) and a no-slip condition \(\boldsymbol{v}=\boldsymbol{0}\) on the fixed wall segment \(\Gamma_{\textrm{wall}}\), as well as a so-called “do-nothing” condition on the artificial boundary segment \(\Gamma_{\textrm{outflow}}\). For this refer to [26], [27]. The velocity field \(\boldsymbol{v}\) and the pressure variable \(p\) are the unknown variables. The parameter \(\nu>0\) denotes the fluid’s viscosity, \(\boldsymbol{f}\) is a given volume force density, and \(\boldsymbol{v}_{0}\) is the initial value.

For the sake of brevity, we present the continuous and semidiscrete in time variational problem for homogeneous Dirichlet boundary conditions only. The mixed boundary conditions in 1 are incorporated in the fully discrete system and derivation of the duality based a posteriori error estimation; cf. Remark 2. Further, they are used in one of our numerical experiments. Throughout, we use standard notation for Sobolev and Bochner spaces. In particular, by \(L_\sigma^2(\Omega)^d\) we denote the \(L^2\)-closure of the set of divergence-free test functions on \(\Omega\). For the existence, uniqueness and regularity of weak solutions to the Navier–Stokes equations such that all of the arguments and terms used below are well-defined, we refer to the literature, cf., e.g., [26], [28], [29].

2.1 Weak Form and Discretization in Space and Time↩︎

For the discretization in space and time we follow the lines of [5], [22]. The variational formulation of 1 reads as:
For \(\boldsymbol{f} \in L^2(I;H^{-1}(\Omega)^{d})\) and \(\boldsymbol{v}_{0}\in L_\sigma^2(\Omega)^d,\) find \(\boldsymbol{u}\coloneq\{\boldsymbol{v},p\}\in X\) such that \[\label{eq:2:weakform} A(\boldsymbol{u})(\boldsymbol{\varphi}) = F(\boldsymbol{\psi}) \quad \forall \boldsymbol{\varphi}\coloneq\{\boldsymbol{\psi},\xi\} \in X\,,\tag{2}\] where the semi-linear form \(A(\cdot)(\cdot)\) and the linear form \(F(\cdot)\) are defined by \[\begin{align} \tag{3} A(\boldsymbol{u})(\boldsymbol{\varphi}) & \coloneq & \int_{I}\big\{(\partial_t \boldsymbol{v},\boldsymbol{\psi}) + a(\boldsymbol{u})(\boldsymbol{\varphi}) \big\} \mathrm{d} t\; + (\boldsymbol{v}(0),\boldsymbol{\psi}(0) )\,, \\[1.5ex] \tag{4} F(\boldsymbol{\psi}) & \coloneq & \int_I(\boldsymbol{f},\boldsymbol{\psi})\;\mathrm{d}t + (\boldsymbol{v}_{0},\boldsymbol{\psi}(0) ) \,, \end{align}\] with the inner semi-linear form \(a(\cdot)(\cdot)\) given by \[\label{eq:5:inSLFa} a(\boldsymbol{u})(\boldsymbol{\varphi}) \coloneq \nu(\nabla \boldsymbol{v}, \nabla \boldsymbol{\psi}) +(\boldsymbol{v}\cdot \nabla \boldsymbol{v}, \boldsymbol{\psi}) -(p,\nabla\cdot\boldsymbol{\psi}) + (\nabla\cdot\boldsymbol{v},\xi)\,.\tag{5}\] Here, \((\cdot,\cdot)\) denotes the inner product on \(L^2(\Omega)\) (or \(L^2(\Omega)^d\)) and the underlying function space \(X\) is defined as \[\label{eq:6:X} X\coloneq\{\boldsymbol{u}=\{\boldsymbol{v},p\}, \boldsymbol{v}\in L^2\big(I;H_0^1(\Omega)^d\big), \partial_t\boldsymbol{v}\in L^2\big(I;H^{-1}(\Omega)^d\big), p\in L^2\big(I;L_0^2(\Omega)\big)\}\,.\tag{6}\] For the semi-discretization in time we use a discontinuous Galerkin method dG(\(r\)) with an arbitrary polynomial degree \(r\ge 0\). Let \(\mathcal{T}_{\tau}\,=\{0\}\cup I_1 \cup \cdots \cup I_N\) be a partition of the closure of the time domain \(\bar{I}=[0,T]\) with left-open subintervals \(I_n\coloneq(t_{n-1},t_n]\), \(n=1,\dots,N\), with \(0=:t_0<t_1<\dots<t_N\coloneq T\) and time step sizes \(\tau_n=t_n-t_{n-1}\) and the global time discretization parameter \(\tau=\max_{n}\,\tau_{n}\). On the subintervals \(I_n\), we define the time-discrete function space \(X_{\tau}^{r}\) as \[\label{eq:7:Def95X95tau95r} \begin{align} X_{\tau}^{r} \coloneq \Big\{ & \boldsymbol{u}_{\tau}=\{\boldsymbol{v}_{\tau},p_{\tau}\}\in L^{2}\big(I; H_0^1(\Omega)^d\times L_0^2(\Omega)\big)\big| \boldsymbol{v}_{\tau}|_{I_{n}}\in \mathcal{P}_{r}\big(I_{n}; H_0^1(\Omega)^d\big),\\ & \boldsymbol{v}_{\tau}(0)\in L^2(\Omega)^d, p_{\tau} \in \mathcal{P}_{r}\big(I_{n}; L_0^2(\Omega)\big)\,, n=1,\dots,N \Big\}\,, \end{align}\tag{7}\] where \(\mathcal{P}_{r}(I_{n}; Y)\) denotes the space of all polynomials in time up to degree \(r\geq0\) on \(I_n\) with values in a Banach space \(Y\). For some discontinuous in time function \(\boldsymbol{v}_{\tau}\in X_{\tau}^{r}\) we define the limits \(\boldsymbol{v}_{\tau,n}^{\pm}\) from above and below of \(\boldsymbol{v}_{\tau}\) at \(t_n\) as well as their jump at \(t_n\) by \[\begin{array}{lcrclcr} \boldsymbol{v}_{\tau,n}^{\pm} & \coloneq & \displaystyle\lim_{t\mapsto t_n\pm0} \boldsymbol{v}_\tau(t) \,, & [\boldsymbol{v}_{\tau}]_{n} & \coloneq & \boldsymbol{v}_{\tau,n}^{+} -\boldsymbol{v}_{\tau,n}^{-} \,. \end{array}\] Then, the semi-discretization in time scheme of the incompressible Navier–Stokes equations 2 reads as follows: Find \(\boldsymbol{u}_{\tau}=\{\boldsymbol{v}_{\tau},p_{\tau}\} \in X_{\tau}^{r}\) such that \[\label{eq:8:time-discrete} A_{\tau}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau}) = F_\tau(\boldsymbol{\psi}_{\tau}) \quad \forall \boldsymbol{\varphi}_{\tau}=\{\boldsymbol{\psi}_{\tau},\xi_{\tau}\} \in X_{\tau}^{r}\,,\tag{8}\] where the semi-discrete semi-linear form \(A_{\tau}(\cdot)(\cdot)\) and linear form \(F_{\tau}(\cdot)\) are defined by \[\tag{9} \begin{align} \tag{10} A_{\tau}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau}) & \coloneq \sum_{n=1}^{N}\int_{I_n}\big\{(\partial_{t} \boldsymbol{v}_{\tau},\boldsymbol{\psi}_{\tau}) + a(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau}) \big\} \mathrm{d} t\; \\ \nonumber & + \sum_{n=2}^{N}([\boldsymbol{v}_\tau]_{n-1},\boldsymbol{\psi}_{\tau,n-1}^+ ) + (\boldsymbol{v}_{\tau,0}^+,\boldsymbol{\psi}_{\tau,0}^+)\,, \\ \tag{11} F_\tau (\boldsymbol{\psi}_{\tau}) & \coloneq \int_I(\boldsymbol{f},\boldsymbol{\psi}_{\tau})\;\mathrm{d}t + (\boldsymbol{v}_{0},\boldsymbol{\psi}_{\tau,0}^+)\,. \end{align}\]

For the discretization in space, let \(\mathcal{T}_h^{n}=\{K\}\) be a triangulation of mesh size \(h_n\) of the domain \(\Omega\) on subinterval \(I_n\) such that \(\overline{\Omega}=\bigcup_{K\in\mathcal{T}_h^{n}}\overline{K}\). For \(t=0\), we let \(\mathcal{T}_h^0\) denote a decomposition of \(\Omega\) with mesh size \(h_0\). The maximal mesh size is denoted by \(h\coloneq \max_{n} h_n\). We define the space-time mesh slabs \(Q_n \coloneq \mathcal{T}_h^n \times I_n\) for \(n=1,\dots, N\). Our mesh adaptation process yields locally refined elements, which is enabled by using hanging nodes. We point out that the global conformity of the finite element approach is preserved since the unknowns at such hanging nodes are eliminated by interpolation between the neighboring ‘regular’ nodes; cf. [10] and [33].

In our GMG approach we apply a local (element-based) Vanka smoother [22]. Numerical studies have shown, that this smoother performs efficiently for local velocity-pressure couplings by using discontinuous pressure approximations. Therefore, we apply the inf-sup stable family of finite element pairs \(Q_k/P_{k-1}^{\textrm{disc}}\), for \(k\geq 2\) [26], [34]). Thus, we built our spatial discretization on the following pair of spaces: \[\label{eq:10:VhQh} \begin{align} \boldsymbol{V}_h^{k}(\mathcal{T}_h^n) &\coloneq \{\boldsymbol{v}_h\in [H^1(\Omega)]^{d}: (\boldsymbol{v}_h)_{|K}\in V^{k}(K)\; \forall K\in\mathcal{T}_h^n\}\,,\\ Q_h^{k-1}(\mathcal{T}_h^n) &\coloneq \{ q_h\in L_0^2(\Omega): (q_h)_{|K}\in Q^{k-1}(K)\; \forall K\in\mathcal{T}_h^n\}\,, \end{align}\tag{12}\] where the finite elements \(\boldsymbol{V}^{k}(K)\) and \(Q^{k-1}(K)\) are defined by \[\label{eq:11:VH} \boldsymbol{V}^{k}(K)\coloneq [Q_k]^d\circ F_K^{-1}\,,\qquad Q^{k-1}(K)\coloneq P_{k-1}^{\textrm{disc}}\circ F_K^{-1}\,,\tag{13}\] and \(F_K\) is the standard multilinear map from the reference element to \(K\). We employ the mapped variant of \(P_{k-1}^{\mathrm{disc}}\) for geometric consistency on curved and non-affine meshes and improved conditioning of the algebraic systems. Here, \({Q}_k, k\geq 2\), denotes the space of polynomials defined on the reference element with maximum degree \(k\) in each variable and \({P}_{k-1}^{\textrm{disc}}\) is the space of polynomials defined on the reference element with maximum degree \(k-1\). Now, by replacing \(H_0^1(\Omega)\) and \(L_0^2(\Omega)\) in the definition of the semi-discrete function space \(X_{\tau}^{r}\) in (7 ) by the spatial finite element spaces defined above, we obtain the fully discrete function space \[\begin{align} \label{eq:12:Def95X95tau95h95r95k} X_{\tau h}^{r,k} \coloneq \Big\{ & \boldsymbol{u}_{\tau h}=\{\boldsymbol{v}_{\tau h},p_{\tau h}\}\, \big|\, \boldsymbol{v}_{\tau h}|_{I_n} \in \mathcal{P}_r\big(I_n; \boldsymbol{V}_{h}^{k}(\mathcal{T}_h^n) \big) \,, \boldsymbol{v}_{\tau h}(0) \in \boldsymbol{V}_{h}^{k}(\mathcal{T}_h^0) \,,\\ & p_{\tau h}|_{I_n} \in \mathcal{P}_r\big(I_n; Q_{h}^{k-1}(\mathcal{T}_h^n) \big)\,, n=1,\dots,N \Big\}\,. \end{align}\tag{14}\] For the treatment of Dirichlet boundary conditions, we use Nitsche’s method (cf. [15]), where the Dirichlet boundary conditions are imposed in weak form. For this, we let \(N_{\Gamma_D}: L^2\big(I;H^{1/2}(\Gamma_D)^d\big)\times X_{\tau h}^{r,k} \to \mathbb{R}\) be defined by \[N_{\Gamma_D}(\boldsymbol{v}, \boldsymbol{\varphi}) x\coloneq \sum_{n=1}^{N}\int_{I_n}\big\{ -(\boldsymbol{v}, \nu\partial_{\boldsymbol{n}}\boldsymbol{\psi}+\xi\,\boldsymbol{n})_{\Gamma_D} +\nu\dfrac{\gamma_1}{h}(\boldsymbol{v}, \boldsymbol{\psi})_{\Gamma_D} +\dfrac{\gamma_2}{h}(\boldsymbol{v}\cdot\boldsymbol{n}, \boldsymbol{\psi}\cdot\boldsymbol{n})_{\Gamma_D} \big\}\;\mathrm{d} t\] for \(\boldsymbol{v}\in L^2\big(I;H^{1/2}(\Gamma_D)^d\) and \(\boldsymbol{\varphi}\in X_{\tau h}^{r,k}\) and some numerical (tuning) parameters \(\gamma_1>0\) and \(\gamma_2>0\) ensuring stability of the discrete system . Further, we define \({B}_{\Gamma_D}: X_{\tau h}^{r,k} \times X_{\tau h}^{r,k}\to \mathbb{R}\) by \[\label{eq:13:nitsche} {B}_{\Gamma_D}({\boldsymbol{u}},\boldsymbol{\varphi})\coloneq \sum_{n=1}^{N}\int_{I_n} -(\nu\partial_{\boldsymbol{n}}{\boldsymbol{v}}-{p}\,\boldsymbol{n},\boldsymbol{\psi})_{\Gamma_D}\,\mathrm{d} t +N_{\Gamma_D}({\boldsymbol{v}}, \boldsymbol{\varphi})\tag{15}\] for \((\boldsymbol{u},\boldsymbol{\varphi})\in X_{\tau h}^{r,k} \times X_{\tau h}^{r,k}\). Finally, the fully discrete problem to Equation 2 reads as follows: Find \(\boldsymbol{u}_{\tau h}=\{\boldsymbol{v}_{\tau h},p_{\tau h}\} \in X_{\tau h}^{r,k}\) such that \[\label{eq:14:fully-discrete} A_{\tau h}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h}) = F_{\tau h}(\boldsymbol{\varphi}_{\tau h}) \quad \forall \boldsymbol{\varphi}_{\tau h}=\{\boldsymbol{\psi}_{\tau h},\xi_{\tau h}\} \in X_{\tau h}^{r,k}\,,\tag{16}\] where the semi-discrete semi-linear form \(A_{\tau h}(\cdot)(\cdot)\) and the linear form \(F_{\tau h}(\cdot)\) are defined by \[\tag{17} \begin{align} \tag{18} A_{\tau h}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h}) \coloneq & A_{\tau}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h}) +{B}_{\Gamma_D}(\boldsymbol{u}_{\tau h},\boldsymbol{\varphi}_{\tau h}) \\ \tag{19} F_{\tau h}(\boldsymbol{\varphi}_{\tau h}) \coloneq & F_{\tau} (\boldsymbol{\psi}_{\tau h}) +N_{\Gamma_D}(\boldsymbol{v}_D, \boldsymbol{\varphi}_{\tau h})\,. \end{align}\] We note that in the case of homogeneous Dirichlet boundary conditions the last term in 19 vanishes since \(\boldsymbol{v}_D=\boldsymbol{0}\).

3 DWR-Based A Posteriori Error Estimation↩︎

In this section, we derive an a posteriori error representation based on the DWR method. This representation is given in terms of a user-chosen goal functional \(J:X\mapsto\mathbb{R}\), in general given as \[\label{eq:3:1:Jgeneral} J(\boldsymbol{u})=\int_0^T J_1(\boldsymbol{u}(t))\;\mathrm{d}t + J_2(\boldsymbol{u}(T))\,,\tag{20}\] where \(J_1\) and \(J_2\) are three times differentiable functionals and either of them may vanish; cf., e.g., [8]. Since we aim at controlling the error due to the discretization in time and in space, we split the a posteriori error representation with respect to \(J\) into the contributions \[\label{eq:3:2:Jsplitted} J(\boldsymbol{u})-J(\boldsymbol{u}_{\tau h}) = J(\boldsymbol{u})-J(\boldsymbol{u}_{\tau}) + J(\boldsymbol{u}_{\tau})-J(\boldsymbol{u}_{\tau h})\,.\tag{21}\] To derive these representations, we introduce the Lagrangian functionals \(\mathcal{L}: X\times X \rightarrow \mathbb{R}\), \(\mathcal{L}_\tau: X_{\tau}^{r} \times X_{\tau}^{r} \rightarrow \mathbb{R}\), and \(\mathcal{L}_{\tau h}: X_{\tau h}^{r,k} \times X_{\tau h}^{r,k} \rightarrow \mathbb{R}\) based on the Euler Lagrange method of constraint optimization: \[\tag{22} \begin{align} \tag{23} \mathcal{L}(\boldsymbol{u},\boldsymbol{z}) & \coloneq J(\boldsymbol{u}) + F(\boldsymbol{w}) - A(\boldsymbol{u})(\boldsymbol{z})\,, \\ \tag{24} \mathcal{L}_{\tau}(\boldsymbol{u}_{\tau},\boldsymbol{z}_{\tau}) & \coloneq J(\boldsymbol{u}_{\tau}) + F_\tau(\boldsymbol{w}_{\tau}) - A_{\tau}(\boldsymbol{u}_{\tau})(\boldsymbol{z}_{\tau})\,, \\ \tag{25} \mathcal{L}_{\tau h}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h}) & \coloneq J(\boldsymbol{u}_{\tau h}) + F_{\tau h} (\boldsymbol{z}_{\tau h}) - A_{\tau h}(\boldsymbol{u}_{\tau h})(\boldsymbol{z}_{\tau h}) \\ \nonumber &=\mathcal{L}_{\tau}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h}) -{B}_{\Gamma_D}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h}) +N_{\Gamma_D}(\boldsymbol{v}_D,\boldsymbol{z}_{\tau h}) \,. \end{align}\] Here, the Lagrange multipliers \(\boldsymbol{z}=\{\boldsymbol{w},q\}\), \(\boldsymbol{z}_{\tau}=\{\boldsymbol{w}_{\tau},q_{\tau}\},\) and \(\boldsymbol{z}_{\tau h}=\{\boldsymbol{w}_{\tau h},q_{\tau h}\}\) are called dual variables in contrast to the primal variables \(\boldsymbol{u}=\{\boldsymbol{v},p\}\), \(\boldsymbol{u}_{\tau}=\{\boldsymbol{v}_{\tau},p_{\tau}\},\) and \(\boldsymbol{u}_{\tau h}=\{\boldsymbol{v}_{\tau h},p_{\tau h}\}\); cf. [4], [5].

Considering the directional derivatives of the Lagrangian functionals, also known as Gâteaux derivatives, with respect to their first argument, i.e. \[\label{eq:3:4:Def95Gateaux95derivative} \mathcal{L}^{\prime}_{\boldsymbol{u}}(\boldsymbol{u},\boldsymbol{z})(\boldsymbol{\varphi}) \coloneq \lim_{t\neq0,t\rightarrow 0} t^{-1}\big\{\mathcal{L}(\boldsymbol{u}+ t\boldsymbol{\varphi},\boldsymbol{z}) -\mathcal{L}(\boldsymbol{u},\boldsymbol{z})\big\}, \quad \boldsymbol{\varphi}\in X\,,\tag{26}\] leads to the so-called dual problems, cf., e.g., [5]. The continuous, semi-discrete, and fully discrete dual solutions \(\boldsymbol{z}~\in~X\), \(\boldsymbol{z}_\tau~\in~X_{\tau}^{r},\) and \(\boldsymbol{z}_{\tau h} \in X_{\tau h}^{r,k}\) are determined by the optimality conditions \[\tag{27} \begin{align} \tag{28} \mathcal{L}^{\prime}_{\boldsymbol{u}}(\boldsymbol{u},\boldsymbol{z})(\boldsymbol{\varphi}) & = 0 \quad \forall \boldsymbol{\varphi}\in X\,, \\ \tag{29} \mathcal{L}^{\prime}_{\tau, \boldsymbol{u}}(\boldsymbol{u}_\tau,\boldsymbol{z}_\tau)(\boldsymbol{\varphi}_{\tau}) & = 0 \quad \forall \boldsymbol{\varphi}_{\tau} \in X_{\tau}^{r}\,, \\ \tag{30} \mathcal{L}^{\prime}_{\tau h, \boldsymbol{u}}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h})(\boldsymbol{\varphi}_{\tau h}) & = 0 \quad \forall \boldsymbol{\varphi}_{\tau h} \in X_{\tau h}^{r,k}\,. \end{align}\] More precisely, the continuous dual solution \(\boldsymbol{z}\in X\) is the solution of \[\label{eq:3:6:dual95continuous} A^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi},\boldsymbol{z}) = J^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi}) \quad \forall \boldsymbol{\varphi}\in X\,,\tag{31}\] where the adjoint semi-linear form \(A^{\prime}(\cdot)(\cdot, \cdot)\) is given by \[\label{eq:3:7:Def95A95prime95u95phi95z} A^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi},\boldsymbol{z}) \coloneq \int_I\big\{(\boldsymbol{\psi},-\partial_{t} \boldsymbol{w}) +a^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi},\boldsymbol{z})\big\}\mathrm{d}t +(\boldsymbol{\psi}(T),\boldsymbol{w}(T))\,.\tag{32}\] We note that for the representation 32 of \(A^{\prime}(\cdot)(\cdot, \cdot)\) integration by parts in time is applied, which is allowed for weak solutions \(\boldsymbol{w}\in X\); cf., e.g., [8]. The derivative \(a^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi},\boldsymbol{z})\) of the inner semi-linear form \(a(\boldsymbol{u})(\boldsymbol{z})\) in \(A^{\prime}\) admits the explicit form \[a^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi},\boldsymbol{z})= \nu(\nabla \boldsymbol{\psi},\nabla \boldsymbol{w}) + (\boldsymbol{v}\cdot \nabla \boldsymbol{\psi},\boldsymbol{w}) + (\boldsymbol{\psi}\cdot\nabla\boldsymbol{v},\boldsymbol{w}) -(\xi,\nabla\cdot\boldsymbol{w}) +(\nabla\cdot\boldsymbol{\psi},q)\,.\] The right-hand side of Eq. 31 is given by \[\label{eq:3:8:Def95J95prime95u95phi} J^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi}) \coloneq \int_I J_1^{\prime}(\boldsymbol{u})(\boldsymbol{\varphi})\mathrm{d}t+J_2^{\prime}(\boldsymbol{u}(T))(\boldsymbol{\varphi}(T))\,.\tag{33}\] Further, the semi-discrete dual solution \(\boldsymbol{z}_\tau\in X_{\tau}^{r}\) and the fully discrete dual solution \(\boldsymbol{z}_{\tau h} \in X_{\tau h}^{r,k}\) satisfy the equations \[\tag{34} \begin{align} \tag{35} A_{\tau}^{\prime}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau},\boldsymbol{z}_{\tau}) & = J^{\prime}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau}) \quad\forall \boldsymbol{\varphi}_{\tau}\in X_{\tau}^{r}\,, \\ \tag{36} A_{\tau h}^{\prime}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h},\boldsymbol{z}_{\tau h}) & = J^{\prime}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h}) \quad \forall \boldsymbol{\varphi}_{\tau h}\in X_{\tau h}^{r,k}\,, \end{align}\] where \(A_{\tau}^{\prime}(\cdot)(\cdot, \cdot)\) and \(A_{\tau h}^{\prime}(\cdot)(\cdot, \cdot)\) are given by \[\begin{align} A_{\tau}^{\prime}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau},\boldsymbol{z}_{\tau}) \coloneq & \sum_{n=1}^N\int_{I_n}\big\{ (\boldsymbol{\psi}_{\tau},-\partial_{t} \boldsymbol{w}_{\tau}) + a^{\prime}(\boldsymbol{u}_{\tau})(\boldsymbol{\varphi}_{\tau},\boldsymbol{z}_{\tau}) \big\}\mathrm{d}t \\ & - \sum_{n=1}^{N-1}(\boldsymbol{\psi}_{\tau,n}^-,[\boldsymbol{w}_{\tau}]_{n}) +(\boldsymbol{\psi}_{\tau,N}^-,\boldsymbol{w}_{\tau,N}^-)\,, \end{align}\] and \[\begin{align} A_{\tau h}^{\prime}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h},\boldsymbol{z}_{\tau h}) \coloneq & A_{\tau}^{\prime}(\boldsymbol{u}_{\tau h})(\boldsymbol{\varphi}_{\tau h},\boldsymbol{z}_{\tau h}) +{B}_{\Gamma_D}(\boldsymbol{\varphi}_{\tau h},\boldsymbol{z}_{\tau h})\,. \end{align}\]

Now, using an abstract error representation result in terms of Lagrangian functionals defined on different discretization levels given in [5], we get the following error representation formulas for the discretization errors in space and time.

3.0.0.1 Proof.

The technique to prove the temporal and spatial error representation formulas goes back to an abstract error result in terms of Lagrangian functionals defined on different discretization levels given by Besier and Rannacher applied to the incompressible Navier–Stokes equations using local projection stabilization (LPS) in[5]. More precisely, we are using the abstract result given in [5] with the following settings: \[\begin{array}{l@{\,}c@{\,}l@{\,}l@{\,}c@{\,}l@{\,}l@{\,}c@{\,}l@{\,}l@{\,}c@{\,}l@{\,}l@{\,}c@{\,}l@{\,}} L & = & \mathcal{L}\,,\;\; & \tilde{L}&=&\mathcal{L}_{\tau}\,,\;\;& Y_1 &=& X \times X\,,\;\;& Y_2 &=& X_{\tau}^{r} \times X_{\tau}^{r}\,,\;\;& Y&\coloneq&Y_1+Y_2\,, \\ L & = & \mathcal{L}_{\tau}\,,\;\; & \tilde{L}&=&\mathcal{L}_{\tau h}\,,\;\;& Y_1 &=& X_{\tau}^{r} \times X_{\tau}^{r}\,,\;\;& Y_2 &=& X_{\tau h}^{r,k} \times X_{\tau h}^{r,k}\,,\;\;& Y&\coloneq&Y_1\,, \end{array}\] where \(\mathcal{L},\mathcal{L}_{\tau},\) and \(\mathcal{L}_{\tau h}\) are the Lagrangian functionals defined in (22 ) and \(Y,Y_1\) and \(Y_2\) are function spaces defined in [5]. The resulting Nitsche terms are treated in the same way as the stabilization terms in the proof of [5]. \(\blacksquare\)

4 Practical Aspects↩︎

In this section, we introduce local error indicators in space and time, present the space-time adaptive algorithm, and describe the geometric multigrid preconditioned Newton–Krylov solver together with selected aspects of the parallel implementation.

4.1 Local Error Indicators in Space and Time↩︎

For the marking step in the adaptive mesh refinement process, we require local versions of the error representations ?? from Theorem 1. Neglecting the higher-order remainder terms \(\mathcal{R}_\tau\) and \(\mathcal{R}_h\) in ?? yields the global a posteriori estimators \[\begin{align} J(\boldsymbol{u})-J(\boldsymbol{u}_{\tau})&\approx \eta_{\tau}\,, \\ J(\boldsymbol{u}_{\tau})-J(\boldsymbol{u}_{\tau h})&\approx \eta_{h}\,, \end{align}\] with \[\begin{align} \eta_{\tau}\coloneq& \frac{1}{2}\rho_{\tau}(\boldsymbol{u}_{\tau})(\boldsymbol{z}-\tilde{\boldsymbol{z}}_{\tau}) + \frac{1}{2}\rho_{\tau}^{\ast}(\boldsymbol{u}_{\tau},\boldsymbol{z}_{\tau}) (\boldsymbol{u}-\tilde{\boldsymbol{u}}_{\tau})\,, \\[1ex] \eta_{h}\coloneq& \frac{1}{2}\rho_{\tau}(\boldsymbol{u}_{\tau h})(\boldsymbol{z}_{\tau}-\tilde{\boldsymbol{z}}_{\tau h}) + \frac{1}{2} \rho_{\tau}^{\ast}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h}) (\boldsymbol{u}_{\tau}-\tilde{\boldsymbol{u}}_{\tau h}) \\ & -\frac{1}{2}N_{\Gamma_D}(\boldsymbol{v}_{D,h}, \tilde{\boldsymbol{z}}_{\tau h}+\boldsymbol{z}_{\tau h}) + \frac{1}{2}{B}_{\Gamma_D}^{\prime}(\boldsymbol{u}_{\tau h}, \boldsymbol{z}_{\tau h})(\tilde{\boldsymbol{u}}_{\tau h}-\boldsymbol{u}_{\tau h},\tilde{\boldsymbol{z}}_{\tau h}-\boldsymbol{z}_{\tau h}) \\ &+{B}_{\Gamma_D}(\boldsymbol{u}_{\tau h}, \boldsymbol{z}_{\tau h}) \,. \end{align}\]

To obtain implementable estimators, we replace the unknown exact and semi-discrete solutions by the fully discrete approximations \(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h}\in X_{\tau h}^{r,k}\). The temporal weights \(\boldsymbol{u}-\tilde{\boldsymbol{u}}_{\tau}\) and \(\boldsymbol{z}-\tilde{\boldsymbol{z}}_{\tau}\), and the spatial weights \(\boldsymbol{u}_{\tau}-\tilde{\boldsymbol{u}}_{\tau h}\) and \(\boldsymbol{z}_{\tau}-\tilde{\boldsymbol{z}}_{\tau h}\), are then approximated as follows.

4.1.0.1 Temporal weights.

We approximate the temporal weights by higher-order reconstructions in time: \[\boldsymbol{u}-\tilde{\boldsymbol{u}}_{\tau} \approx \mathrm{E}_{\tau}^{(r+1)}\boldsymbol{u}_{\tau h}-\boldsymbol{u}_{\tau h}\,, \qquad \boldsymbol{z}-\tilde{\boldsymbol{z}}_{\tau} \approx \mathrm{E}_{\tau}^{(r+1)}\boldsymbol{z}_{\tau h}-\boldsymbol{z}_{\tau h}\,.\] Here \(\mathrm{E}_{\tau}^{(r+1)}\) denotes a reconstruction operator in time that lifts the discrete solution to a piecewise polynomial of degree \(r+1\), see [35], [36].

4.1.0.2 Spatial weights.

We approximate the spatial weights by patch-wise higher-order interpolation: \[\boldsymbol{u}_{\tau}-\tilde{\boldsymbol{u}}_{\tau h} \approx \mathrm{I}_{2h}^{(2k)}\boldsymbol{u}_{\tau h}-\boldsymbol{u}_{\tau h}\,, \qquad \boldsymbol{z}_{\tau}-\tilde{\boldsymbol{z}}_{\tau h} \approx \mathrm{I}_{2h}^{(2k)}\boldsymbol{z}_{\tau h}-\boldsymbol{z}_{\tau h}\,.\] Here \(\mathrm{I}_{2h}^{(2k)}\) denotes a patch-wise interpolation operator that acts on a patch of size \(2h\) and lifts the discrete solution to a piecewise polynomial of degree \(2k\), see [37].

With these approximations, the implemented temporal and spatial estimators are given by \[\tag{37} \begin{align}\tag{38} \tilde{\eta}_{\tau}\approx&\frac{1}{2}\rho_{\tau}(\boldsymbol{u}_{\tau h})(\mathrm{E}_{\tau}^{(r+1)}\boldsymbol{z}_{\tau h}-\boldsymbol{z}_{\tau h}) + \frac{1}{2}\rho_{\tau}^{\ast}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h})(\mathrm{E}_{\tau}^{(r+1)}\boldsymbol{u}_{\tau h}-\boldsymbol{u}_{\tau h})\,, \\ \tag{39} \tilde{\eta}_{h}\approx&\frac{1}{2}\rho_{\tau}(\boldsymbol{u}_{\tau h})(\mathrm{I}_{2h}^{(2k)}\boldsymbol{z}_{\tau h}-\boldsymbol{z}_{\tau h}) + \frac{1}{2}\rho_{\tau}^{\ast}(\boldsymbol{u}_{\tau h},\boldsymbol{z}_{\tau h})(\mathrm{I}_{2h}^{(2k)}\boldsymbol{u}_{\tau h}-\boldsymbol{u}_{\tau h}) \\\nonumber & -\frac{1}{2}N_{\Gamma_D}(\boldsymbol{v}_{D,h}, \mathrm{I}_{2h}^{(2k)}\boldsymbol{z}_{\tau h}+\boldsymbol{z}_{\tau h}) + \frac{1}{2}{B}_{\Gamma_D}^{\prime}(\boldsymbol{u}_{\tau h}, \boldsymbol{z}_{\tau h})(\mathrm{I}_{2h}^{(2k)}\boldsymbol{u}_{\tau h}-\boldsymbol{u}_{\tau h},\mathrm{I}_{2h}^{(2k)}\boldsymbol{z}_{\tau h}-\boldsymbol{z}_{\tau h}) \\ \nonumber & +{B}_{\Gamma_D}(\boldsymbol{u}_{\tau h}, \boldsymbol{z}_{\tau h}) \,. \end{align}\]

For mesh marking, these implemented estimators are decomposed into local contributions applying integration by parts on each element \(K\in\mathcal{T}_h^n\) to the viscous term. Following [4], [5], [10], we write \[\label{eq:localized-errors-iso} \tilde{\eta}_{\tau}=\sum_{n=1}^{N}\sum_{K\in\mathcal{T}_h^n}\tilde{\eta}_{\tau}^{K,n}, \qquad \tilde{\eta}_h=\sum_{n=1}^{N}\sum_{K\in\mathcal{T}_h^n}\tilde{\eta}_{h}^{K,n}\,.\tag{40}\] The local indicators \(\tilde{\eta}_{\tau}^{K,n}\) and \(\tilde{\eta}_{h}^{K,n}\) are then used for slab and element marking.

4.2 Space-Time Adaptive Algorithm↩︎

We present a condensed version of the adaptive solution algorithm based on tensor-product space-time slabs \(Q_n\). For further details on the slab data structures, we refer to [38] and [35].

In Remark 4, we introduced slab-wise refinement fractions. In this work, we either keep these fractions constant, as in previous work, or adapt them according to the following rule.

4.3 Geometric Multigrid Preconditioned Newton-Krylov Method↩︎

The space-time finite element discretization of the incompressible Navier–Stokes equations leads to the nonlinear algebraic system 16 .

To solve this problem, we employ a Newton–Krylov method preconditioned by geometric multigrid. Newton’s method is applied in the weak form \[\label{newton} A_{\tau h}'(\boldsymbol{u}^{(l-1)})(\boldsymbol{w}^{(l)},\boldsymbol{\varphi}_{\tau h}) = F_{\tau h}(\boldsymbol{\varphi}_{\tau h}) - A_{\tau h}(\boldsymbol{u}^{(l-1)})(\boldsymbol{\varphi}_{\tau h}), \quad \boldsymbol{u}^{(l)}=\boldsymbol{u}^{(l-1)}+\boldsymbol{w}^{(l)}.\tag{41}\] Here \(A_{\tau h}'(\boldsymbol{u}^{(l-1)})\) denotes the Fréchet derivative of the semilinear form \(A_{\tau h}\) at \(\boldsymbol{u}^{(l-1)}\). Algebraically, this corresponds to the Jacobian matrix evaluated at \(\boldsymbol{u}^{(l-1)}\) and applied to the Newton correction \(\boldsymbol{w}^{(l)}\). The Jacobian can be computed analytically, see [8]. To reduce numerical costs, the Jacobian is reused over several Newton iterations and consecutive slabs, and is reassembled only when the observed residual reduction factor \(\|\boldsymbol{r}^{(l)}\|/\|\boldsymbol{r}^{(l-1)}\|\) exceeds a prescribed threshold \(\rho_{\max}\). In our computations, \(\rho_{\max}\) is chosen empirically and typically lies between \(0.05\) and \(0.1\).

Each Newton step requires the solution of a large, sparse, nonsymmetric, and indefinite linear system arising from the saddle-point structure of the Navier–Stokes equations. We use the Generalized Minimal Residual method (GMRES) [40], [41] to solve these linear systems. To accelerate convergence, GMRES is preconditioned by geometric multigrid. On each slab \(I_n\), the multigrid method is built on a hierarchy of nested discrete saddle-point spaces \[\mathcal{X}_0^n \subset \mathcal{X}_1^n \subset \cdots \subset \mathcal{X}_L^n,\] associated with nested spatial meshes \(\mathcal{T}_0^n,\dots,\mathcal{T}_L^n=\mathcal{T}_h^n\). On each level, high-frequency errors are damped by a smoothing operator, while low-frequency components are transferred to coarser levels by \(L^2\)-projection-based restriction and interpolative prolongation. On the coarsest level, a direct solver is employed. Under the standard multigrid assumptions, one V-cycle has linear complexity in the number of degrees of freedom.

To treat the saddle-point structure efficiently and to obtain a parallelizable smoother, we employ a Vanka-type smoother [42] within the multigrid cycle. It is defined by \[\label{vanka0} \mathcal{V}_\ell(A_\ell) = \sum_{K\in\mathcal{T}_\ell} R_K^\top \bigl[R_KA_\ell R_K^\top\bigr]^{-1}R_K,\tag{42}\] where \(A_\ell\) denotes the system matrix on level \(\ell\), and \(R_K\) is the restriction operator to the degrees of freedom associated with the element \(K\in\mathcal{T}_\ell\). One smoothing step is then given by \[\label{vanka} \mathcal{S}(\boldsymbol{A}_\ell,\boldsymbol{b}_\ell,\boldsymbol{x}_k) = \boldsymbol{x}_k + \omega\,\mathcal{V}_\ell(A_\ell)\, (\boldsymbol{b}_\ell-\boldsymbol{A}_\ell\boldsymbol{x}_k),\tag{43}\] with damping parameter \(\omega\in(0,1]\). In our experience, two to five pre- and post-smoothing steps are sufficient for robust convergence. The discontinuous pressure approximation is beneficial for the Vanka smoother, since it localizes the pressure degrees of freedom and thereby improves the efficiency of the element-wise local solves.

In summary, the nonlinear space-time finite element system is linearized by Newton’s method and the resulting linear systems are solved by GMRES preconditioned by geometric multigrid with a Vanka-type smoother. For a more detailed description of the solver architecture and its parallelization, we refer to [22].

4.4 Parallel implementation↩︎

The implementation of the adaptive algorithm is distributed by MPI and uses the parallel infrastructure of deal.II for error estimation, marking, refinement, and multigrid setup across all processes. No stage is executed on a dedicated master rank. Collective communication is still required for the assembly of global indicators and for the coarse-grid solve.

At present, load balancing and repartitioning are not used. To maintain scalability after several adaptive cycles, we plan to incorporate dynamic repartitioning in future work.

5 Numerical Examples↩︎

In this section, we investigate the performance properties of the proposed approach, which integrates goal-oriented error control with iterative solver techniques. In particular, we study the convergence behavior, computational efficiency, as well as stability properties of the underlying space-time adaptive algorithm (Algorithm 1) introduced in Sect. 4.2. These properties are examined by means of the following three numerical examples.

1. Exact Analytical Solution: The first example is an academic test problem including a given analytical solution that goes back to [5]. It is used to validate accuracy and efficiency of the derived approach for different discretizations up to polynomial order 3 in time and 4 in space.

2. DFG Benchmark 2D-3: The second example is the famous DFG flow around cylinder benchmark 2D-3 that simulates the time-periodic behavior of a fluid in a pipe with a circular obstacle that goes back to [43]. The goal functional is chosen to control the mean drag.

We evaluate the accuracy and efficiency of our proposed approach in terms of several quantities. Firstly, to assess the accuracy of the duality based space-time error estimator, we compute the effectivity index defined as the ratio of the estimated to the exact error, given by \[\label{eq:Ieff} I_{\mathrm{eff}}\coloneq \left|\frac{ \tilde{\eta}_{\tau} + \tilde{\eta}_{h}}{J(\boldsymbol{u})-J(\boldsymbol{u}_{\tau h})}\right|\,.\tag{44}\] Secondly, to measure the gain in accuracy over the increase of work, we determine the empirical convergence rates, averaged over all DWR-loops, \[\label{eq:mean-rates} \overline{\operatorname{rate}}_{N_{\mathrm{ST}}}^{\|e_T\|} \quad \text{and} \quad \overline{\operatorname{rate}}_{W_{\mathrm{tot}}}^{\|e_T\|}\tag{45}\] of the reduction of the error in the goal quantity with respect to the increase of the number of space-time degrees of freedom \(N_{\mathrm{ST}}\) and the total work \(W_{\mathrm{tot}}\). The latter is measured by the cumulative linear solver work \[\label{eq:work} W_{\mathrm{tot}} = W_p + W_d,\quad\text{with primal and dual contributions}\tag{46}\] \[W_p \coloneq \sum_{n=1}^{N_{\mathrm{time}}} n_{\mathrm{L},n}^{p} N_{\mathrm{DoF},n}^{p}, \; W_d \coloneq \sum_{n=1}^{N_{\mathrm{time}}} n_{\mathrm{L},n}^{d} N_{\mathrm{DoF},n}^{d}\,,\] where \(n_{\mathrm{L},n}^{p}\) and \(n_{\mathrm{L},n}^{d}\) denote the numbers of primal and dual linear iterations on slab \(Q_n = \mathcal{T}_h^n \times I_n\), and \(N_{\mathrm{DoF},n}^{p}\), \(N_{\mathrm{DoF},n}^{d}\) the corresponding numbers of spatial degrees of freedom. The nonlinear primal problem is solved by Newton’s method. We therefore also report the number of nonlinear iterations on slab \(Q_n\), denoted by \(n_{\mathrm{NL},n}\). Accordingly, we let \(n_{\mathrm{L},n}^{p}\) denote the total number of linear iterations accumulated over the \(n_{\mathrm{NL},n}\) Newton steps on slab \(Q_n\). By \(n_{\mathrm{NL}}\), \(n_{\mathrm{L}}^{p}\), and \(n_{\mathrm{L}}^{d}\) we denote the total numbers of nonlinear, primal linear, and dual linear iterations over all slabs within one DWR loop \(\ell\). In addition, we record \[\label{eq:iteration-indicators} \min_{\ell=1,\dots,\ell_{\max}} \Bigl(\frac{n_{\mathrm{NL}}}{N_{\mathrm{time}}}\Bigr),\;\; \max_{\ell=1,\dots,\ell_{\max}} \Bigl(\frac{n_{\mathrm{NL}}}{N_{\mathrm{time}}}\Bigr),\;\; \min_{\ell=1,\dots,\ell_{\max}} \Bigl(\frac{n_{\mathrm{L}}^{p}}{n_{\mathrm{NL}}}\Bigr),\;\; \max_{\ell=1,\dots,\ell_{\max}} \Bigl(\frac{n_{\mathrm{L}}^{p}}{n_{\mathrm{NL}}}\Bigr),\tag{47}\] together with the relative dual contribution \[\label{eq:dual-work-fraction} \frac{W_d}{W_{\mathrm{tot}}}.\tag{48}\] Throughout, by \(\ell\) we denote the refinement level or DWR loop number, \(N_{\text{time}}\) the number of subintervals in time, \(N_{\text{space}}^{\text{max}}\) the maximum number of degrees of freedom associated with a spatial triangulation within one loop and \(N_{\text{tot}}\) the number of total degrees of freedom in space and time. We consider different finite element approximations dG(\(r\))-\(Q_k/P_{k-1}^{\mathrm{disc}}\) denoting a discontinuous in time approximation of polynomial degree \(r\geq 0\) as well as a continuous velocity and discontinuous pressure approximations in space of polynomial degree \(k\geq2\).

The implementation is based on the deal.II finite element library [1]. The tests are run on a single node with 2 Intel Xeon Platinum 8360Y CPUs and 1024   RAM of the HPC cluster HSUper at Helmut Schmidt University.

5.1 Example 1: Exact Analytical Solution↩︎

As a first test case, we consider a model problem on the space-time domain \(\Omega\times I = [0,1]^2\times [0, 1]\) with prescribed solution given for the velocity \(\boldsymbol{v} \colon \Omega\times I \to \mathbb{R}^2\) and pressure \(p \colon \Omega\times I \to \mathbb{R}\) by \[\label{eq:conv-test} \begin{align} \boldsymbol{v}(\boldsymbol{x},\,t) &= \sin(t) \begin{pmatrix} \sin^2(\pi x) \sin(\pi y) \cos(\pi y) \\ \sin(\pi x) \cos(\pi x) \sin^2(\pi y) \end{pmatrix}, \\ p(\boldsymbol{x},\,t) &= \sin(t) \sin(\pi x) \cos(\pi x) \sin(\pi y) \cos(\pi y)\,. \end{align}\tag{49}\] We choose the kinematic viscosity as \(\nu \in \{1, 10^{-3}\}\) and choose the external force \(\boldsymbol{f}\) such that the solution 49 satisfies 1 . The initial velocity is prescribed as zero and homogeneous Dirichlet boundary conditions are imposed on \(\partial\Omega\) for all times \[\boldsymbol{v} = \boldsymbol{0}\text{ on }\Omega\times \{0\},\quad \boldsymbol{v} = \boldsymbol{0}, \text{ on } \partial \Omega\times (0, T]\,.\] The initial space-time mesh \(\mathcal{T}_{h}\otimes\mathcal{M}_{\tau}\) is a uniform tensor-product partition of the space-time domain \(\Omega\times I\). We use discretizations with varying polynomial degrees \(r \in \{0,\dots,\,3\}\) in time, where we use \(\mathbb{Q}_{r+1}/\mathbb{P}_{r}^{\mathrm{disc}}\) elements in space, except for the \(dG(0)\) case where we use \(\mathbb{Q}_{2}/\mathbb{P}_{1}^{\mathrm{disc}}\) elements. The goal quantity is chosen to control the \(L^2\)-error \(e_N^-, e_N^-=u(\boldsymbol{x},T)-u_{\tau h}(T^{-})\), at the final time point \(T=1\), given by \[\label{eq:JL2final} J_T(u)= \frac{\big(u(\boldsymbol{x},T),e_N^-\big)}{\|e_N^-\|}\,.\tag{50}\] In Tables ¿tbl:tab:dG0?–¿tbl:tab:dG3? we present the development of the total discretization error \(J(e)=\|e_T\|\), the approximated spatial and temporal error estimators \(\eta_h\) and \(\eta_\tau\) as well as the effectivity index \(I_{\mathrm{eff}}\) during an adaptive refinement process for various discretizations. For all discretizations considered, we observe a consistent and satisfactory reduction of the discretization error \(\|e_T\|\). Concerning the accuracy of the employed error estimators, all cases exhibit reliable quantitative agreement with the discretization error, as indicated by effectivity indices \(I_{\mathrm{eff}}\) close to one (cf. the last column of Tables ¿tbl:tab:dG0?–¿tbl:tab:dG3?). From the perspective of computational efficiency in space–time adaptive algorithms, it is crucial to ensure an equilibrated reduction of the temporal as well as spatial discretization error. Referring to this, we point out for all test cases a good equilibration of the spatial and temporal error indicators \(\eta_h\) and \(\eta_\tau\) in the course of the refinement process (cf. column 6 and 7 of Tables ¿tbl:tab:dG0?–¿tbl:tab:dG3?).

5.5pt

4.5pt

4.5pt

4.5pt

Across all four fully adaptive runs, the DWR estimator in combination with adaptive refinement yields a consistent reduction of the goal error with increasing space-time resolution, while the corresponding effectivity indices remain in a reasonable range. The reported mean rates show that the decay of the goal error with respect to total work is close to its decay with respect to the space-time number of degrees of freedom, indicating that the adaptive refinement translates effectively into computational efficiency. A compact comparison of these approximation and cost quantities is given in Table 1.

The algebraic behavior is summarized in Figure 4, which shows stable nonlinear iteration counts per slab and only moderate growth of the primal linear iterations per nonlinear step over the adaptive loops. This indicates that the slab-wise geometric multigrid preconditioned Newton–Krylov method remains robust under successive adaptive refinement. At the same time, the relative dual work fraction stays at a comparable level across all discretizations, showing that the adjoint solve contributes a noticeable but controlled part of the overall computational cost.

In Fig. 5 we plot the share of the primal problem, the dual problem as well as the error estimation over the entire wallclock time and its percentage share for one adaptive run including all DWR loops for the different discretizations corresponding to Tables ¿tbl:tab:dG0?–¿tbl:tab:dG3?. While the total runtime from dG(\(1\))-\(Q_2/P_1^{\textrm{disc}}\) to dG(\(3\))-\(Q_4/P_3^{\textrm{disc}}\) increases progressively as expected, the percentage share of the three units remains stable and nearly unchanged across all discretizations. In comparison to the primal problem, the dual problem requires only one third of its computational effort, while the error estimator itself is negligible at less than one percent. The fact that the total runtime of the dG(\(0\))-\(Q_2/P_1^{\textrm{disc}}\) discretization is slightly higher to that of the dG(\(1\))-\(Q_2/P_1^{\textrm{disc}}\) one is due to the large number of slabs required in the first case.

Finally, in Fig. 6 and Fig. 7 we visualize the distribution of the adaptively determined time step sizes \(\tau_n\) and minimal spatial element diameter \(\min_{K\in\mathcal{T}_h^{n,\ell_{\max}}} h_K\) within the slab \(Q_n\) over the entire time interval \(I\), respectively, corresponding to the final DWR loop of the different discretizations listed in Tables ¿tbl:tab:dG0?–¿tbl:tab:dG3?. For all discretizations, the step and element sizes are relatively large at the beginning, but they decrease continuously over time, reaching their minimal values at the final time \(T=1\). We note that the initial steps correspond to the initial distribution of the temporal and spatial meshes. This behavior is consistent with the choice of the underlying goal functional 50 , which specifically targets the control of the \(L^2\)-error at the final time \(T\). Furthermore, we observe that the total length of the subintervals is significantly reduced when using higher-order discretizations compared to the lowest-order case.

5.2 Example 2: DFG Benchmark 2D-3↩︎

We consider the well-established benchmark problem of laminar flow around a cylinder introduced by Schäfer and Turek [43], precisely the unsteady 2D-3 test setting. In this case, the flow around a circular obstacle with Reynolds number \(\mathrm{Re}=100\) (the kinematic viscosity is set to \(\nu=0.001\)) is driven by a parabolic inflow profile whose amplitude is modulated in time by a \(\sin\) function. The geometry of the flow domain is shown in Fig. 8. We note that the coarse mesh used here is the same as in previous studies [5], [44]. The final time point is set to \(T=8\). The goal quantity is chosen to control the error of the mean drag coefficient. This coefficient is defined by \[\label{eq:Jmeandrag} J_{\textrm{drag}}(u)= \displaystyle-\int_I\{(\partial_t\boldsymbol{v},\hat{\boldsymbol{\psi}}_d)+a(\boldsymbol{u})(\hat{\boldsymbol{\varphi}}_d)\}\,\mathrm{d}t\,,\tag{51}\] where \(\hat{\boldsymbol{\varphi}}_d=(\hat{\boldsymbol{\psi}}_d,0)^\top\) fulfills \(\hat{\boldsymbol{\psi}}_d|_S=(|I|^{-1}C,0)^\top\,, \;\; \hat{\boldsymbol{\psi}}_d|_{\partial\Omega\textbackslash S}=\boldsymbol{0}\) with the constant \(C=2(U_{\mathrm{mean}}^2d)^{-1}\). Here, the mean inflow velocity is given by \[U_{\mathrm{mean}}=2/3\;\boldsymbol{v}_D(0,0.205,0)=\sin (\pi t/8)\] given by Eq. 52 and \(d\) denotes the diameter of the obstacle with boundary \(S=S=\partial B_r(0.2,0.2)\). We take the reference value for the goal functional from [5] and set it to 1.6072872.

The flow is driven by the Dirichlet inflow profile \(\boldsymbol{v}=\boldsymbol{v}_D\) on the left boundary \(\Gamma_{\mathrm{in}} \coloneq \{0\}\times[0,H],\) given by \[\label{eq:inflow} \boldsymbol{v}_D(\boldsymbol{x},t) = \begin{pmatrix} \sin\!\left(\frac{\pi}{8}t\right)\dfrac{6y(H-y)}{H^2}\\[0.3em] 0 \end{pmatrix} \qquad \text{on }\Gamma_{\mathrm{in}},\tag{52}\] where \(H=0.41\) denotes the channel height. On the wall boundary \(\Gamma_{\mathrm{wall}}\) and on the obstacle boundary we prescribe no-slip boundary conditions \(\boldsymbol{v}=\boldsymbol{0}\). On the outflow boundary \(\Gamma_{\mathrm{out}}\) we impose a do-nothing condition [27].

In Table ¿tbl:tab:dfg-drag-all? we present the results of an adaptive refinement run together with a comparison of the algebraic performance of a GMG-preconditioned Newton-Krylov method and an ILU method using a dG(\(1\))-\(Q_2/P_1^{\mathrm{disc}}\) discretization. The adaptive refinement reduces the mean drag error 51 monotonically. The drag value converges to the reference value and the effectivity indices stay close to one on all loops. This shows that the DWR estimator provides a reliable approximation of the goal error. The observed decay of the goal error with respect to the space-time number of degrees of freedom is at least linear and becomes higher on the finer loops. Measured against total algebraic work, the same trend is obtained more efficiently with GMG than with ILU.

This is consistent with the solver data. The GMG-preconditioned Newton-Krylov method remains robust over the full adaptive sequence, whereas ILU requires substantially more nonlinear iterations and much larger primal linear iteration counts. The larger number of nonlinear iterations for ILU can be explained by the lower accuracy of the linear solves in the Newton steps. Since the resulting Newton corrections are less accurate than with GMG, the convergence of the outer Newton iteration is correspondingly weaker. Overall, the results indicate that the adaptive strategy is effective for the drag goal quantity 51 and that GMG is the more robust and efficient preconditioner for this problem.

The temporal evolution of the drag coefficient is compared for the different DWR refinement loops in Fig. 9. We note that the curves converge by increasing the adaptive refinement (cf. loop 3 and 4) and are in good comparison to the results given in [5] and [44], respectively.

Finally, in Fig. 10 we present the discrete primal solution on the final DWR loop as well as the adaptively refined spatial meshes for all four DWR loops at time point \(t=5\). The spatial refinement remains concentrated close to the cylinder and in the wake, as confirmed by the selected meshes and solution snapshots. The drag and lift histories are smooth and consistent with the benchmark results, which indicates that the adaptive discretization resolves the relevant time dependent features and the Stokes projections avoids spurious oscillations.

6 Conclusion↩︎

We have developed a goal-oriented space-time adaptive finite element method for the incompressible Navier–Stokes equations based on the Dual Weighted Residual method. The approach combines discontinuous Galerkin time discretization, inf-sup stable mixed finite elements in space, adaptive refinement in space and time, and a slab-wise Newton–Krylov solver with geometric multigrid preconditioning. The numerical results show that this combination yields effective control of the goal quantities on dynamically adapted meshes.

The main computational benefit of the DWR approach is its problem-driven refinement. The estimator identifies those regions in space and time that are most relevant for the chosen output and thereby focuses resolution where it matters for the quantity of interest. At the same time, the results show that this accuracy comes at a substantial cost. In addition to the primal problem, the adaptive loop requires the repeated solution of a dual problem and the evaluation of weighted residual contributions. The efficiency of the overall approach therefore depends not only on the quality of the estimator, but equally on the robustness of the underlying iterative solvers.

The Newton iteration counts are essentially unaffected by the choice of linear preconditioner, whereas the GMRES iteration counts depend strongly on it. Thus, the additional cost of the DWR framework is governed primarily by the linear solver performance for the primal and dual subproblems. In this setting, geometric multigrid provides the required robustness, while ILU does not yield a competitive performance. This is especially important for higher-order discretizations, where the linear systems become more demanding and the dual solve adds a further substantial burden. For goal-oriented space-time adaptivity, robust preconditioning is therefore essential.

Overall, the results show that DWR-based space-time adaptivity is effective for goal-oriented error control, provided that it is paired with scalable linear algebra. The practical cost of the method is determined by the interplay between estimator quality and solver performance. Future work should therefore address both aspects together, in particular through cheaper dual solves, more selective dual updates, and improved multilevel solvers for higher-order and larger-scale space-time discretizations.

Acknowledgement↩︎

Computational resources (HPC cluster HSUper) have been provided by the project hpc.bw, funded by dtec.bw - Digitalization and Technology Research Center of the Bundeswehr. dtec.bw is funded by the European Union - NextGenerationEU.

References↩︎

---

1. \(^{,\ddagger}\) Helmut Schmidt University Hamburg, Faculty of Mechanical and Civil Engineering, Chair of Numerical Mathematics, Holstenhofweg 85, 22043 Hamburg, Germany,
   bruchhaeuser@hsu-hamburg.de (\(^\ast\)corresponding author), bause@hsu-hamburg.de↩︎

2. University of Magdeburg, Institute for Analysis and Numerics, Universitätsplatz 2, 39104 Magdeburg, Germany, nils.margenberg@ovgu.de.de↩︎