Tensor-Product Space-Time Goal-Oriented Error Control and Adaptivity With Partition-of-Unity Dual-Weighted Residuals for Nonstationary Flow Problems

Julian Roth; Jan Philipp Thiele; Uwe Köcher; Thomas Wick

doi:10.1515/cmam-2022-0200

Article Publicly Available

Tensor-Product Space-Time Goal-Oriented Error Control and Adaptivity With Partition-of-Unity Dual-Weighted Residuals for Nonstationary Flow Problems

Julian Roth , Jan Philipp Thiele , Uwe Köcher and Thomas Wick

Published/Copyright: May 31, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 24 Issue 1

Abstract

In this work, the dual-weighted residual method is applied to a space-time formulation of nonstationary Stokes and Navier–Stokes flow. Tensor-product space-time finite elements are being used to discretize the variational formulation with discontinuous Galerkin finite elements in time and inf-sup stable Taylor–Hood finite element pairs in space. To estimate the error in a quantity of interest and drive adaptive refinement in time and space, we demonstrate how the dual-weighted residual method for incompressible flow can be extended to a partition-of-unity based error localization. We substantiate our methodology on 2D benchmark problems from computational fluid mechanics.

Keywords: Tensor-Product Space-Time Finite Elements; Dual-Weighted Residuals; Mesh Adaptivity; Dynamically Changing Meshes; Incompressible Navier–Stokes Equations

MSC 2010: 65N30; 65N50; 76D05

1 Introduction

This work is devoted to the efficient numerical solution of the incompressible Navier–Stokes equations employing adaptive Galerkin finite elements in space and time. Adaptivity is based on a posteriori goal-oriented error control using a partition-of-unity dual-weighted residual (PU-DWR) approach. The closest studies to ours are on the one hand [15] as well as the PhD thesis [54] in which similar concepts for full space-time adaptivity for the Navier–Stokes equations were developed. The main differences are two-fold in that we work with tensor-product space-time finite elements and that we utilize a partition-of-unity for localization (originally proposed for stationary problems in [51]) rather than the filtering approach [18]. On the other hand, another closely related study is [58] in which the current PU-DWR method was developed and analyzed for parabolic problems. Using discontinuous Galerkin (dG) test functions in time decouples the formulation and yields a natural localization for which a PU in time is not mandatory. However, we view the PU as a systematic extension to full space-time settings that holds for discontinuous and continuous Galerkin time discretizations. For instance, full continuous Galerkin space-time settings in 4D as in the work of Schafelner [52], error control and adaptivity are naturally achieved with a space-time PU localization. Furthermore, higher-order time derivatives and integration by parts in time in the localization process can be included (as this was one of the main motivations of Richter–Wick [51]).

Employing (full) space-time adaptivity using the DWR method dates back to [55], which developed and tested the method for parabolic problems. Shortly afterward, the space-time DWR method was applied to the incompressible Navier–Stokes equations in [54, 15] and for second-order hyperbolic problems (i.e. the elastic wave equation) in [48, 7]. In this regard, the first space-time finite element formulation (without adaptivity) for incompressible flow was proposed in [45]. Employing space-time (therein called G ⁢ 2 ) for turbulent flow and goal-oriented adaptivity for the efficient computation of the mean drag was done in [31].

Transport problems with coupled flow with an emphasis on adaptivity of the transport part was investigated in [12], and the extension to a multirate framework was suggested in [21]. Goal-oriented error control with space-time formulations, but only applied to temporal error control for parabolic problems, the Navier–Stokes equations and fluid-structure interaction was undertaken in [42, 43, 26], respectively. Space-time methods with residual-based error control and adaptivity can be found for instance in [38, 52]. Finally, [44] proposed goal-oriented space-time adaptivity with optimal control and [36, 37] developed residual-based adaptivity in optimal control.

The main objectives of this work are to develop space-time PU-DWR for the incompressible Navier–Stokes equations. Therein, we emphasize the following novelties in comparison to existing work. Our approach allows for natural higher-order time discretizations by setting the polynomial degrees of the trial and test basis functions since a general variational discontinuous Galerkin (dG) formulation in time is used and there is no need for the explicit derivation of time-stepping schemes (as for instance done in [55, 15]). Therefore, our framework is more generally extendable to other PDEs since time stepping formulations, e.g. a piecewise linear and discontinuous in time dG ( 1 ) , neither have to be explicitly derived nor implemented for coupled multiphysics problems. To this end, we consider piecewise quadratic time discretizations ( dG ( 2 ) ) for the dual problems without major methodological efforts. Precisely, the defining nodes for the piecewise polynomial basis functions, which are called support points, can be flexibly, arbitrarily and independently chosen within the temporal elements. Such support points might be the quadrature points of the used quadrature for the temporal integration or some other distribution in the temporal element such that the temporal polynomial degree holds. For instance, in this work, we apply (and compare) schemes generated by the Gauss–Lobatto and the Gauss–Legendre formulae. Finally, our implementation is done in the modern finite element library deal.II [6, 5]. These developments are accompanied with a detailed computational analysis in form of error reductions and studies of effectivity indices. To better understand the behavior, the (linear) Stokes equations are considered as well. Therein, Galerkin discretizations using tensor products in time and space are employed. In time, discontinuous Galerkin dG ( r ) elements with polynomial degree r ≥ 0 are used for the primal forward problem. The adjoint problem is approximated with globally higher-order finite elements. In space, the classical Taylor–Hood element (i.e. continuous Galerkin cG ( s ) , with s ≥ 1 with s = 2 for the velocities and s = 1 for the pressure) is employed for the primal problem and a higher-order approximation for the adjoint. A delicate issue are pressure-robust discretizations on dynamically changing adaptive meshes for which we employ a divergence-free projection from [16], but we make several remarks to more modern methods. These developments yield expressions for error indicators in time and space, which are then utilized to formulate a final adaptive algorithm.

The outline of this paper is as follows. In Section 2, the space-time formulation and discretization are introduced. Next, in Section 3, the space-time PU-DWR is developed and error estimators and error indicators are derived. Afterward, in Section 4, several benchmark tests are adopted to substantiate our algorithmic developments. Our work is summarized in Section 5.

2 Space-Time Formulation and Discretization

We model viscous fluid flow with constant density and temperature by the Navier–Stokes equations [27, 57, 49, 28, 59, 29].

Formulation 1

Formulation 1 (Incompressible Isothermal Navier–Stokes Equations)

Find the vector-valued velocity

v : ( 0 , T ) × Ω → R d with ⁢ d ∈ { 2 , 3 }

and the scalar-valued pressure p : ( 0 , T ) × Ω → R such that

∂ t v − ∇ x ⋅ σ + ( v ⋅ ∇ x ) ⁢ v = 0 in ⁢ ( 0 , T ) × Ω , ∇ x ⋅ v = 0 in ⁢ ( 0 , T ) × Ω , v = v D on ⁢ ( 0 , T ) × ∂ Ω , v ⁢ ( 0 ) = v 0 in ⁢ Ω ,

where we use the non-symmetric Cauchy stress tensor for the stress 𝝈 defined as

σ := σ ⁢ ( v p ) = − p ⁢ I + ν ⁢ ∇ x v .

Here, 𝜈 is the kinematic viscosity, v 0 is the initial velocity and v D is a possibly time-dependent Dirichlet boundary value.

Remark 1

By omitting the nonlinear convection term ( v ⋅ ∇ x ) ⁢ v , we obtain the Stokes problem.

We employ the following notation for the function spaces. Spatial function spaces are denoted by V := V 1 × V 2 , where V 1 is a vector-valued function space and V 2 is a scalar-valued function space. Without indices, the function space of the continuous problem statement is V ⁢ ( Ω ) = H 0 1 ⁢ ( Ω ) d × L 2 ⁢ ( Ω ) . The space of lowest-order Taylor–Hood elements (quadratic FE for velocity, linear FE for pressure) is V h 2 / 1 ⁢ ( T h ) , and the space of fourth-order FE for velocity and second-order FE for pressure on a given spatial mesh T h is V h 4 / 2 ⁢ ( T h ) . Temporal function spaces are denoted by 𝑋. Without indices, 𝑋 and X ~ are the continuous level function spaces without and with discontinuities respectively. The space of discontinuous finite elements in time of 𝑟-th order is X k dG ( r ) , and the space of discontinuous finite elements in time of ( r + 1 ) -th order is X k dG ( r + 1 ) . We denote spatio-temporal function spaces by combining a spatial and a temporal function space,

X ⁢ ( I , V ) := { ( v p ) | v ∈ L 2 ⁢ ( I , V 1 ) , ∂ t v ∈ L 2 ⁢ ( I , V 1 ∗ ) , p ∈ L 2 ⁢ ( I , V 2 ) , ∫ Ω p ⁢ ( t ) ⁢ d x = 0 ⁢ for all ⁢ t ∈ I }

that will be the space-time Hilbert space on which we define the continuous weak formulation. Here, the dual space of V 1 is denoted by V 1 * = L ⁢ ( V 1 , R ) . This is a good choice for the continuous level function space since we have a continuous embedding into a function space where the fluid velocity is continuous on I ¯ (see [22]). This ensures that the initial condition is well-defined for the velocity.

Remark 2

When we have homogeneous Neumann boundary conditions σ ⋅ n = 0 on Γ out ⊂ ∂ Ω , the so-called do-nothing (outflow) condition [30], the pressure does not require normalization for its uniqueness. As the Dirichlet constraints are only set on part of the boundary, i.e. Γ D = ∂ Ω ∖ Γ out , the continuous spatial function space is V = H 0 , Γ D 1 ⁢ ( Ω ) d × L 2 ⁢ ( Ω ) . Then we have the spatio-temporal function space

X ⁢ ( I , V ) := { ( v p ) | v ∈ L 2 ⁢ ( I , V 1 ) , ∂ t v ∈ L 2 ⁢ ( I , V 1 ∗ ) , p ∈ L 2 ⁢ ( I , V 2 ) } .

Finally, to simplify notation, we introduce

( f , g ) := ( f , g ) L 2 ⁢ ( Ω ) := ∫ Ω f ⋅ g ⁢ d x , ( ( f , g ) ) := ( f , g ) L 2 ⁢ ( I , L 2 ⁢ ( Ω ) ) := ∫ I ( f , g ) ⁢ d t .

In this notation, f ⋅ g represents the Euclidean inner product if 𝑓 and 𝑔 are scalar- or vector-valued and the Frobenius inner product if 𝑓 and 𝑔 are matrices.

For a given initial value v 0 ∈ L 2 ⁢ ( Ω ) d , the weak formulation reads as follows: find

U := ( v p ) ∈ ( v D 0 ) + X ⁢ ( I , V ⁢ ( Ω ) )

such that

A ⁢ ( U ) ⁢ ( Φ ) := ( ( ∂ t v , Φ v ) ) + a ⁢ ( U ) ⁢ ( Φ ) + ( v ⁢ ( 0 ) − v 0 , Φ v ⁢ ( 0 ) ) = 0 for all ⁢ Φ := ( Φ v Φ p ) ∈ X ⁢ ( I , V ⁢ ( Ω ) )

where, for the Navier–Stokes equations, we have

a ⁢ ( U ) ⁢ ( Φ ) := − ( ( p , ∇ x ⋅ Φ v ) ) + ν ⁢ ( ( ∇ x v , ∇ x Φ v ) ) + ( ( ( v ⋅ ∇ x ) ⁢ v , Φ v ) ) + ( ( ∇ x ⋅ v , Φ p ) ) ,

and for the Stokes equations, we have

a ⁢ ( U ) ⁢ ( Φ ) := − ( ( p , ∇ x ⋅ Φ v ) ) + ν ⁢ ( ( ∇ x v , ∇ x Φ v ) ) + ( ( ∇ x ⋅ v , Φ p ) ) .

2.1 Discretization in Time

Let T k := { I m , m = 1 , … , M } , with I m := ( t m − 1 , t m ) , be a partitioning of time, i.e. I ¯ = [ 0 , T ] = ⋃ m = 1 M I ¯ m . Then, in an intermediate step similar to spatial dG methods, we can define a broken function space that allows for discontinuities at the temporal grid points,

(2.1) X ~ ⁢ ( T k , V ) := { ( v p ) ∈ L 2 ⁢ ( I , L 2 ⁢ ( Ω ) d + 1 ) , ( v p ) | I m ∈ X ⁢ ( I m , V ) ⁢ for all ⁢ I m ∈ T k } .

As local continuous embeddings into C ⁢ ( I ¯ m , V 1 ) exist for 𝒗, we can define the limits of 𝒇 at time t m from above and from below for a function f ∈ X ~ ⁢ ( T k , V ) as

f m ± := lim ϵ ↘ 0 f ⁢ ( t m ± ϵ ) ,

and the jump of the function value of 𝒇 at time t m as [ f ] m := f m + − f m − . With this, we can introduce a weak formulation with discontinuities on the continuous level by interval-wise evaluation of the temporal integrals and inserting jumps between the temporal elements.

Find

U := ( v p ) ∈ ( v D 0 ) + X ~ ⁢ ( T k , V ⁢ ( Ω ) )

such that

(2.2) A ~ ⁢ ( U ) ⁢ ( Φ k ) := ∑ m = 1 M ∫ I m ( ∂ t v , Φ v ) − ( p , ∇ x ⋅ Φ v ) + ν ⁢ ( ∇ x v , ∇ x Φ v ) ⁢ d ⁢ t + ∑ m = 1 M ∫ I m ( ( v ⋅ ∇ x ) ⁢ v , Φ v ) + ( ∇ x ⋅ v , Φ p ) ⁢ d ⁢ t + ∑ m = 1 M − 1 ( [ v ] m , Φ m v , + ) + ( v , 0 + ⁢ ( 0 ) − v 0 , Φ 0 v , + ⁢ ( 0 ) ) = 0 for all ⁢ Φ := ( Φ v Φ k p ) ∈ X ~ ⁢ ( T k , V ⁢ ( Ω ) ) .

As functions f ∈ X ⁢ ( I , V ⁢ ( Ω ) ) are continuous on 𝐼, [ f ] m = 0 holds for all m = 1 , … , M . Therefore, it follows that X ⁢ ( I , V ⁢ ( Ω ) ) ⊂ X ~ ⁢ ( T k , V ⁢ ( Ω ) ) holds as well and the above weak formulation is consistent. By discretizing the space-time function space (2.1) in time with the discontinuous Galerkin method of order r ∈ N 0 ( dG ( r ) ), we obtain the semi-discrete space as

X k dG ( r ) ( T k , V ) := { ( v k p k ) | v k | I m ∈ P r ⁢ ( I m , V 1 ) , p k | I m ∈ P r ⁢ ( I m , V 2 ) ⁢ for all ⁢ 1 ≤ m ≤ M , ∫ Ω p k ( t ) d x = 0 for all t ∈ I } ⊂ X ~ ( T k , V ) .

Here, P r ⁢ ( I m , Y ) is the space of polynomials of order 𝑟, which map from the time interval I m into the function space 𝑌. The dG ( r ) time discretization for the cases r = 0 and r = 1 is illustrated in Figure 1.

Figure 1

dG ( r ) time discretization.

$(a) dG ( 0 ) \operatorname*{dG}(0)$

(a)

dG ( 0 )

$(b) dG ( 1 ) \operatorname*{dG}(1)$

(b)

dG ( 1 )

It has been derived in [54] for the Navier–Stokes equations that we obtain the backward Euler scheme by using suitable numerical quadrature for the temporal integrals of the dG ( 0 ) discretization. Therefore, it can be seen as a variant of that scheme. Additionally, the corresponding time-stepping scheme of dG ( 1 ) has been formulated therein. The dG ( r ) time discretization of the Navier–Stokes equations reads as follows: find

U k := ( v k p k ) ∈ ( v D 0 ) + X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) )

such that

A k ⁢ ( U k ) ⁢ ( Φ k ) := ∑ m = 1 M ∫ I m ( ∂ t v k , Φ k v ) − ( p k , ∇ x ⋅ Φ k v ) + ν ⁢ ( ∇ x v k , ∇ x Φ k v ) ⁢ d ⁢ t + ∑ m = 1 M ∫ I m ( ( v k ⋅ ∇ x ) ⁢ v k , Φ k v ) + ( ∇ x ⋅ v k , Φ k p ) ⁢ d ⁢ t + ∑ m = 1 M − 1 ( [ v k ] m , Φ k , m v , + ) + ( v k , 0 + ⁢ ( 0 ) − v 0 , Φ k , 0 v , + ⁢ ( 0 ) ) = 0 for all ⁢ Φ k := ( Φ k v Φ k p ) ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) .

2.2 Discretization in Space

Next, we describe the spatial discretization of the dG ( r ) formulation for the (Navier–)Stokes equations with continuous finite elements. Let T h be a shape-regular mesh of the domain Ω ⊂ R d . The elements K ∈ T h are a non-overlapping covering of Ω ¯ , where the elements 𝐾 are quadrilaterals for d = 2 and hexahedrals for d = 3 . The discretization parameter ℎ is defined element-wise as h | K = h K = diam ⁡ ( K ) . After adaptive refinement in space, h K can vary across elements and we need to account for hanging nodes in the assembly of the matrices and vectors. These are degrees of freedom that are located in the middle of the neighboring elements’ edge ( d = 2 ) or face ( d = 3 ). Furthermore, on different time intervals, we allow for changing spatial meshes, i.e. possibly T h ⁢ ( I m ) ≠ T h ⁢ ( I n ) , but for simplicity, we drop the dependence on the time interval and only mention it when it is not obvious from the context.

During the spatial discretization of the mixed system of the (Navier–)Stokes equations, we need to account for the inf-sup or Ladyzhenskaya–Babuska–Brezzi (LBB) condition

inf q ∈ L 2 ⁢ ( Ω ) sup φ ∈ H 0 1 ⁢ ( Ω ) d ( q , ∇ ⋅ φ ) ∥ q ∥ L 2 ⁢ ( Ω ) ⁢ ∥ φ ∥ H 1 ⁢ ( Ω ) d ≥ γ > 0 ,

which guarantees well-posedness (existence, uniqueness and stability) of the pressure solution under suitable boundary conditions and assumptions on the domain regularity; see e.g. [28, 50]. In the spatial discretization, we now replace L 2 ⁢ ( Ω ) and H 0 1 ⁢ ( Ω ) d with conforming finite element spaces that still fulfill a discrete version of this stability estimate. For this, we define the mixed finite element space

V h s v / s p ⁢ ( T h ) := { ( v p ) ∈ C ⁢ ( Ω ¯ ) d + 1 | v | K ∈ Q s v ⁢ ( K ) d , p | K ∈ Q s p ⁢ ( K ) ⁢ for all ⁢ K ∈ T h } ,

where Q s ⁢ ( K ) is being constructed by mapping tensor-product polynomials of degree 𝑠 from the master element ( 0 , 1 ) d to the element 𝐾. In our computations, we use the function space V h 2 / 1 , the so-called Taylor–Hood elements [33] with quadratic finite elements for the velocity and linear finite elements for the pressure. Additionally, we use as V h 4 / 2 an enriched space with fourth-order finite elements for the velocity and second-order finite elements for the pressure. In R 2 , we could have also used V h 3 / 2 as an enriched space since Q 3 / Q 2 was shown to be an inf-sup stable Taylor–Hood pair [56], but this is not necessarily the case in 3D anymore. Hence, we approximate the enriched velocity by Q 4 instead of Q 3 finite elements in space. The fully discretized problem then reads as follows: find

U k ⁢ h := ( v k ⁢ h p k ⁢ h ) ∈ ( v D 0 ) + X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) )

such that

(2.3) A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Φ k ⁢ h ) := ∑ m = 1 M ∫ I m ( ∂ t v k ⁢ h , Φ k ⁢ h v ) − ( p k ⁢ h , ∇ x ⋅ Φ k ⁢ h v ) + ν ⁢ ( ∇ x v k ⁢ h , ∇ x Φ k ⁢ h v ) ⁢ d ⁢ t + ∑ m = 1 M ∫ I m ( ( v k ⁢ h ⋅ ∇ x ) ⁢ v k ⁢ h , Φ k ⁢ h v ) + ( ∇ x ⋅ v k ⁢ h , Φ k ⁢ h p ) ⁢ d ⁢ t + ∑ m = 1 M − 1 ( [ v k ⁢ h ] m , Φ k ⁢ h , m v , + ) + ( v k ⁢ h , 0 + ⁢ ( 0 ) − v 0 , Φ k ⁢ h , 0 v , + ⁢ ( 0 ) ) = 0 for all ⁢ Φ k ⁢ h := ( Φ k ⁢ h v Φ k ⁢ h p ) ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

2.3 Tensor-Product Space-Time Discretization

The above space-time analysis is required for the theory of error estimation for the time-dependent partial differential equations. However, this formulation is usually transformed into a time-stepping scheme in the literature by applying a suitable numerical quadrature to the temporal integrals (see also Section 2.1). Instead, we will follow the approach of Bruchhäuser, Köcher and Bause [21, Section 4], which uses tensor-product spaces to directly work with the fully discrete weak formulation (2.3). In this approach, the space-time cylinder ( 0 , T ) × Ω is divided and discretized into space-time slabs

S l n := ( ⋃ m = l n I m × T h m )

on which the weak formulation is solved sequentially. The spatial mesh is slabwise constant, i.e. T h l = ⋯ = T h n .

For Hilbert function spaces H 1 ⁢ ( I ) and H 2 ⁢ ( Ω ) , the closure of the Hilbert space tensor product H 1 ⁢ ( I ) ⊗ ^ H 2 ⁢ ( Ω ) is isometric to H 1 ⁢ ( I , H 2 ⁢ ( Ω ) ) , and the isometry also holds for tensor products of finite-dimensional subspaces [47, Chapter 1.2]. From this, it follows that f ⁢ ( t ) ⁢ g ⁢ ( x ) with f ∈ H 1 ⁢ ( I ) , g ∈ H 2 ⁢ ( Ω ) can be identified with a function h ⁢ ( t , x ) ∈ H 1 ⁢ ( I , H 2 ⁢ ( Ω ) ) . Hence, we can define the space-time tensor-product finite element basis on S l n by taking the tensor product of the temporal discontinuous Galerkin basis functions and the spatial continuous Galerkin basis functions. This allows for the extension of a simple finite element code for a stationary Stokes problem by creating a temporal mesh and adding temporal loops over elements, quadrature points and degrees of freedom in the assembly of the FEM matrices. Furthermore, unlike the time-stepping formulations of (2.3) where each polynomial degree in time needs to be derived and implemented separately, changing the temporal degree of the space-time discretization can be performed simply by changing the polynomial degree of the temporal finite elements and our code theoretically allows for hp-adaptivity in time since the temporal degree could be chosen on each slab individually. Another benefit of using space-time slabs for adaptive refinement is that the spatial mesh is fixed for the length of the slab which is sometimes desirable for numerical stability. On the other hand, using slabs with multiple time elements comes at the cost of solving larger linear equation systems.

Remark 3

Remark 3 (Support Points for the Temporal Elements)

As we use discontinuous elements in time, we are not limited to placing the outermost temporal degrees of freedom on the element edges. Instead, we can use arbitrary points without the need for recalculation of a time-stepping scheme. For spatial dG methods, Gauss–Legendre support points have been shown to provide a better discretization [23].

2.4 Linearization of the Semilinear Form

For the Navier–Stokes equations, we need to include the nonlinear convection term ( ( ( v ⋅ ∇ x ) ⁢ v , Φ v ) ) in our variational formulation A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Φ k ⁢ h ) = 0 . Therefore, we need to apply linearization to obtain a solvable linear equation system. In this paper, we decided to solve this variational root finding problem by using the residual-based Newton method. The resulting nonlinear iteration step j + 1 ∈ N reads as follows: find

δ ⁢ U k ⁢ h := ( δ ⁢ v k ⁢ h δ ⁢ p k ⁢ h ) ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) )

such that

A k ⁢ h , U ′ ⁢ ( U k ⁢ h j ) ⁢ ( δ ⁢ U k ⁢ h , Φ k ⁢ h ) = − A k ⁢ h ⁢ ( U k ⁢ h j ) ⁢ ( Φ k ⁢ h ) for all ⁢ Φ k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

Obtain the new iterate by

U k ⁢ h j + 1 = U k ⁢ h j + δ ⁢ U k ⁢ h .

Remark 4

Note that we will also need to assemble the directional derivative A k ⁢ h , U ′ for the adjoint problem (3.2) with the minor difference that there the trial and test functions are interchanged. This corresponds to transposing the system matrix from a Newton step.

The derivative of the complete semilinear form is given by

A k ⁢ h , U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h , Φ k ⁢ h ) = ∑ m = 1 M ∫ I m ( ∂ t δ ⁢ v k ⁢ h , Φ k ⁢ h v ) − ( δ ⁢ p k ⁢ h , ∇ x ⋅ Φ k ⁢ h v ) + ν ⁢ ( ∇ x δ ⁢ v k ⁢ h , ∇ x Φ k ⁢ h v ) ⁢ d ⁢ t + ∑ m = 1 M ∫ I m ( ( δ ⁢ v k ⁢ h ⋅ ∇ x ) ⁢ v k ⁢ h + ( v k ⁢ h ⋅ ∇ x ) ⁢ δ ⁢ v k ⁢ h , Φ k ⁢ h v ) + ( ∇ x ⋅ δ ⁢ v k ⁢ h , Φ k ⁢ h p ) ⁢ d ⁢ t + ∑ m = 1 M − 1 ( [ δ ⁢ v k ⁢ h ] m , Φ k ⁢ h , m v , + ) + ( δ ⁢ v k ⁢ h , 0 + ⁢ ( 0 ) , Φ k ⁢ h , 0 v , + ⁢ ( 0 ) ) .

The Newton algorithm was taken from [60] and the corresponding GitHub code^[1], which could naturally be reused in our space-time setting. For the numerical tests, we used max line search = 10 as the maximum number of line search steps and λ = 0.6 as the line search damping parameter. We set the maximum number of Newton steps to max Newton = 10 and the tolerance TOL to 10 − 10 .

Since we are dealing with a time-dependent problem, the initial guess can be chosen as the solution of the previous slab evaluated at its end time, i.e. for the slabs S m − 1 m − 1 := ( t m − 2 , t m − 1 ) × T h and S m m := ( t m − 1 , t m ) × T h , we define

U k ⁢ h m , 0 ≡ U k ⁢ h m − 1 ⁢ ( t m − 1 ) .

Remark 5

Alternatively one could create the initial guess U k ⁢ h m , 0 as the temporal extrapolant of the entire solution from the previous slab, e.g.

U k ⁢ h m , 0 ⁢ ( t ) := U k ⁢ h m − 1 ⁢ ( t m − 2 ) + ( t − t m − 2 ) ⋅ U k ⁢ h m − 1 ⁢ ( t m − 1 ) − U k ⁢ h m − 1 ⁢ ( t m − 2 ) t m − 1 − t m − 2

for linear extrapolation in time.

Additionally, we need to prescribe the non-homogeneous Dirichlet conditions on the initial guess. However, the boundary conditions for the update δ ⁢ U k ⁢ h m become homogeneous, i.e. we enforce δ ⁢ U k ⁢ h m = 0 on the Dirichlet boundary.

In order to save some computational cost, to compute the directional derivative A k ⁢ h , U ′ ⁢ ( U k ⁢ h m , j ) ⁢ ( δ ⁢ U k ⁢ h m , Φ k ⁢ h ) at each step of the Newton method, in our program, we use the approximation

A k ⁢ h , U ′ ⁢ ( U k ⁢ h m , j + 1 ) ⁢ ( ⋅ , Φ k ⁢ h ) ≈ A k ⁢ h , U ′ ⁢ ( U k ⁢ h m , j ) ⁢ ( ⋅ , Φ k ⁢ h )

when the reduction rate

θ j := ∥ R ⁢ ( U k ⁢ h m , j + 1 ) ∥ ∥ R ⁢ ( U k ⁢ h m , j ) ∥

is sufficiently small. We use these so-called simplified-Newton steps when θ j ≤ θ max := 0.1 .

2.5 Divergence-Free Projection

In fluid dynamics, working with dynamically changing meshes leads to non-physical oscillations in the values of goal functionals, which has been proven in [16] for a Stokes model problem. Further analysis of the numerical artifacts caused by the violation of the inf-sup condition, due to the usage of time-varying spatial meshes, can be found in [19, 9, 2, 39]. To avoid the resulting numerical artifacts, Besier and Wollner [16] proposed to either repeat the last time step on the finite element space of the new slab or to use a divergence-free projection for the initial condition when switching to a different spatial mesh. Here, we decided to focus on the latter and computationally analyze the impact of a divergence-free L 2 and a divergence-free H 0 1 projection [16].

Remark 6

We notice that this approach uses a global correction, which is expensive and requires the solution of an additional problem. Using the Helmholtz decomposition, Linke [40] showed that the continuous velocity solution is invariant to irrotational forcings. However, conforming mixed finite elements lose these properties on the discrete level leading to non-physical solutions. Various fixes were proposed and studied using local projections to H ⁢ ( div ) -conforming, i.e. discretely divergence-free, finite elements inside the weak formulation of the (Navier–)Stokes equations [40, 41, 35, 39]. However, these studies for fluid flow have all been done on triangular elements, and applying them with corresponding quadrilateral elements did not yield the desired results for our fluid flow formulations. In incompressible solid mechanics, quadrilateral elements with such pressure-robust schemes seem to work though [10, 11].

This problem can be formalized as follows. Let

S m m := ( t m − 1 , t m ) × T h m and S m + 1 m + 1 := ( t m , t m + 1 ) × T h m + 1

be two space-time slabs with different spatial triangulations T h m ≠ T h m + 1 . Then we need to find a projection Π of the initial condition

U k ⁢ h ⁢ ( t m ) := ( v k ⁢ h ⁢ ( t m ) p k ⁢ h ⁢ ( t m ) ) ∈ ( v D ⁢ ( t m ) 0 ) + V h s v / s p ⁢ ( T h m )

into the function space of the next slab

( v D ⁢ ( t m ) 0 ) + V h s v / s p ⁢ ( T h m + 1 )

such that the projected initial condition is divergence free, i.e.

( ∇ x ⋅ ( Π ⁢ v k ⁢ h ⁢ ( t m ) ) , Φ h p ) = 0 for all ⁢ Φ h p ∈ V h s p ⁢ ( T h m + 1 ) .

To achieve this goal, Besier and Wollner [16] presented and analyzed a divergence-free L 2 projection and a divergence-free H 0 1 projection.

Divergence-free L 2 projection: find
Π ⁢ U k ⁢ h ⁢ ( t m ) = ( Π ⁢ v k ⁢ h ⁢ ( t m ) Π ⁢ p k ⁢ h ⁢ ( t m ) ) ∈ ( v D ⁢ ( t m ) 0 ) + V h s v / s p ⁢ ( T h m + 1 )
such that
( Π ⁢ v k ⁢ h ⁢ ( t m ) , Φ h v ) − ( Π ⁢ p k ⁢ h ⁢ ( t m ) , ∇ x ⋅ Φ h v ) + ( ∇ x ⋅ ( Π ⁢ v k ⁢ h ⁢ ( t m ) ) , Φ h p ) = ( v k ⁢ h ⁢ ( t m ) , Φ h v ) for all ⁢ Φ h := ( Φ h v Φ h p ) ∈ V h s v / s p ⁢ ( T h m + 1 ) .
Divergence-free H 0 1 projection: find
Π ⁢ U k ⁢ h ⁢ ( t m ) = ( Π ⁢ v k ⁢ h ⁢ ( t m ) Π ⁢ p k ⁢ h ⁢ ( t m ) ) ∈ ( v D ⁢ ( t m ) 0 ) + V h s v / s p ⁢ ( T h m + 1 )
such that
( ∇ x ( Π ⁢ v k ⁢ h ⁢ ( t m ) ) , ∇ x Φ h v ) − ( Π ⁢ p k ⁢ h ⁢ ( t m ) , ∇ x ⋅ Φ h v ) + ( ∇ x ⋅ ( Π ⁢ v k ⁢ h ⁢ ( t m ) ) , Φ h p ) = ( ∇ x v k ⁢ h ⁢ ( t m ) , ∇ x Φ h v ) for all ⁢ Φ h := ( Φ h v Φ h p ) ∈ V h s v / s p ⁢ ( T h m + 1 ) .

3 Space-Time Dual-Weighted Residual Error Estimation

In many practical applications, we are in general not interested in the entire solution of the Navier–Stokes equations but only in some quantity of interest. In this chapter, we will describe how we can use the Dual-Weighted Residual (DWR) method [14, 13] (see also other references [8, 55, 15, 1, 46]) to estimate the overall error of our finite element solution in this quantity of interest and how we can localize the error with a space-time partition of unity (PU) [51, 58] to refine the spatial and temporal meshes.

3.1 DWR for Parabolic PDEs

In the following analysis, we will work with homogeneous boundary conditions to simplify the presentation. Let us consider a goal functional J : X ~ ⁢ ( I , V ⁢ ( Ω ) ) → R of the form

J ⁢ ( U ) = ∫ 0 T J 1 ⁢ ( U ⁢ ( t ) ) ⁢ d t + J 2 ⁢ ( U ⁢ ( T ) ) ,

which represents our physical quantity of interest. We now want to minimize the difference between the quantity of interest of the continuous solution 𝑼 and its FEM approximation U k ⁢ h , i.e.

J ⁢ ( U ) − J ⁢ ( U k ⁢ h )

subject to the constraint that the variational formulation is being satisfied. As in constrained optimization, we define the Lagrange functionals

L : X ~ ⁢ ( I , V ⁢ ( Ω ) ) × X ~ ⁢ ( I , V ⁢ ( Ω ) ) → R , ( U , Z ) ↦ J ⁢ ( U ) − A ⁢ ( U ) ⁢ ( Z ) , L k : X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) × X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) → R , ( U k , Z k ) ↦ J ⁢ ( U k ) − A k ⁢ ( U k ) ⁢ ( Z k ) , L k ⁢ h : X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) × X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) → R , ( U k ⁢ h , Z k ⁢ h ) ↦ J ⁢ ( U k ⁢ h ) − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k ⁢ h ) .

The stationary points ( U , Z ) , ( U k , Z k ) and ( U k ⁢ h , Z k ⁢ h ) of the Lagrange functionals ℒ, L k and L k ⁢ h need to satisfy the Karush–Kuhn–Tucker first-order optimality conditions. Firstly, the stationary points are solutions to the equations

L Z ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ Z ) = 0 for all ⁢ δ ⁢ Z ∈ X ~ ⁢ ( I , V ⁢ ( Ω ) ) , L k , Z ′ ⁢ ( U k , Z k ) ⁢ ( δ ⁢ Z k ) = 0 for all ⁢ δ ⁢ Z k ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) , L k ⁢ h , Z ′ ⁢ ( U k ⁢ h , Z k ⁢ h ) ⁢ ( δ ⁢ Z k ⁢ h ) = 0 for all ⁢ δ ⁢ Z k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

We call these equations the primal problems and their solutions 𝑼, U k and U k ⁢ h the primal solutions. Secondly, the stationary points must also satisfy the equations

L U ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ U ) = 0 for all ⁢ δ ⁢ U ∈ X ~ ⁢ ( I , V ⁢ ( Ω ) ) , L k , U ′ ⁢ ( U k , Z k ) ⁢ ( δ ⁢ U k ) = 0 for all ⁢ δ ⁢ U k ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) , L k ⁢ h , U ′ ⁢ ( U k ⁢ h , Z k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h ) = 0 for all ⁢ δ ⁢ U k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

These equations are called the adjoint or dual problems and their solutions 𝒁, Z k and Z k ⁢ h are the adjoint solutions. To be able to decide whether the spatial or temporal mesh should be refined, we split the total discretization error into

J ⁢ ( U ) − J ⁢ ( U k ⁢ h ) = ( J ⁢ ( U ) − J ⁢ ( U k ) ) + ( J ⁢ ( U k ) − J ⁢ ( U k ⁢ h ) ) ≈ η k + η h ,

where η k stands for the temporal error estimate and η h stands for the spatial error estimate. Following [15, Theorem 5.2] and neglecting the remainder terms, which are of higher order, we get the error estimates

(3.1) J ⁢ ( U ) − J ⁢ ( U k ) ≈ 1 2 ( ρ ( U k ) ( Z − Z k ) + ρ * ( U k , Z k ) ( U − U k ) ) = : η k , J ⁢ ( U k ) − J ⁢ ( U k ⁢ h ) ≈ 1 2 ( ρ ( U k ⁢ h ) ( Z k − Z k ⁢ h ) + ρ * ( U k ⁢ h , Z k ⁢ h ) ( U k − U k ⁢ h ) ) = : η h ,

where we denote the primal and adjoint residuals as

ρ ⁢ ( U ) ⁢ ( δ ⁢ Z ) := L k ⁢ h , Z ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ Z ) , ρ * ⁢ ( U , Z ) ⁢ ( δ ⁢ U ) := L k ⁢ h , U ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ U ) .

We now have two problems to solve: the primal and the adjoint problem. All the remaining variables in the error estimator can be computed via higher-order interpolation or via interpolation to a lower-order space.

3.2 Primal Problem

To derive the primal problem, we need to compute the Gâteaux derivatives of the Lagrange functionals ℒ, L k and L k ⁢ h with respect to the adjoint solution 𝒁. Using the definition of the continuous (𝐴), the time-discrete ( A k ) and the fully discrete ( A k ⁢ h ) variational formulations, we observe that they are linear in the test functions, and hence we have

L Z ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ Z ) = − A ⁢ ( U ) ⁢ ( δ ⁢ Z ) = 0 for all ⁢ δ ⁢ Z ∈ X ~ ⁢ ( I , V ⁢ ( Ω ) ) , L k , Z ′ ⁢ ( U k , Z k ) ⁢ ( δ ⁢ Z k ) = − A k ⁢ ( U k ) ⁢ ( δ ⁢ Z k ) = 0 for all ⁢ δ ⁢ Z k ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) , L k ⁢ h , Z ′ ⁢ ( U k ⁢ h , Z k ⁢ h ) ⁢ ( δ ⁢ Z k ⁢ h ) = − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ Z k ⁢ h ) = 0 for all ⁢ δ ⁢ Z k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

To get the primal solution, we simply need to solve the weak formulation of the Navier–Stokes problem.

3.3 Adjoint Problem

To derive the adjoint problem, we need to compute the Gâteaux derivatives of the Lagrange functionals ℒ, L k and L k ⁢ h with respect to the primal solution 𝑼. We then get

L U ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ U ) = J U ′ ⁢ ( U ) ⁢ ( δ ⁢ U ) − A U ′ ⁢ ( U ) ⁢ ( δ ⁢ U , Z ) = 0 for all ⁢ δ ⁢ U ∈ X ~ ⁢ ( I , V ⁢ ( Ω ) ) , L k , U ′ ⁢ ( U k , Z k ) ⁢ ( δ ⁢ U k ) = J U ′ ⁢ ( U k ) ⁢ ( δ ⁢ U k ) − A k , U ′ ⁢ ( U k ) ⁢ ( δ ⁢ U k , Z k ) = 0 for all ⁢ δ ⁢ U k ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) , L k ⁢ h , U ′ ⁢ ( U k ⁢ h , Z k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h ) = J U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h ) − A k ⁢ h , U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h , Z k ⁢ h ) = 0 for all ⁢ δ ⁢ U k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

Therefore, to obtain the adjoint solution, we need to solve an additional equation, the adjoint problem

(3.2) A k ⁢ h , U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h , Z k ⁢ h ) = J U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h ) .

Remark 7

For linear PDEs, like the Stokes equations, and linear goal functionals, like the mean drag functional, the adjoint problems simplify to

L U ′ ⁢ ( U , Z ) ⁢ ( δ ⁢ U ) = J ⁢ ( δ ⁢ U ) − A ⁢ ( δ ⁢ U ) ⁢ ( Z ) = 0 for all ⁢ δ ⁢ U ∈ X ~ ⁢ ( I , V ⁢ ( Ω ) ) , L k , U ′ ⁢ ( U k , Z k ) ⁢ ( δ ⁢ U k ) = J ⁢ ( δ ⁢ U k ) − A k ⁢ ( δ ⁢ U k ) ⁢ ( Z k ) = 0 for all ⁢ δ ⁢ U k ∈ X k dG ( r ) ⁢ ( T k , V ⁢ ( Ω ) ) , L k ⁢ h , U ′ ⁢ ( U k ⁢ h , Z k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h ) = J ⁢ ( δ ⁢ U k ⁢ h ) − A k ⁢ h ⁢ ( δ ⁢ U k ⁢ h ) ⁢ ( Z k ⁢ h ) = 0 for all ⁢ δ ⁢ U k ⁢ h ∈ X k dG ( r ) ⁢ ( T k , V h s v / s p ⁢ ( T h ) ) .

In this scenario, we have the error identities

(3.3) J ⁢ ( U ) − J ⁢ ( U k ) = − A k ⁢ ( U k ) ⁢ ( Z − Z k ) , J ⁢ ( U k ) − J ⁢ ( U k ⁢ h ) = − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k − Z k ⁢ h ) .

We show this error identity for the spatial error contribution. Using the linearity of the goal functional and the definition of the adjoint problem, we have

J ⁢ ( U k ) − J ⁢ ( U k ⁢ h ) = J ⁢ ( U k − U k ⁢ h ) = A k ⁢ ( U k − U k ⁢ h ) ⁢ ( Z k ) = A k ⁢ h ⁢ ( U k − U k ⁢ h ) ⁢ ( Z k ) = − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k ) = − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k − Z k ⁢ h ) .

For the temporal error identity, we can follow the same steps with the continuous level weak formulation including discontinuities A ~ (2.2).

By Remark 7, we can thus write the adjoint problem (3.2) for the Stokes problem as

A k ⁢ h ⁢ ( δ ⁢ U k ⁢ h ) ⁢ ( Z k ⁢ h ) = J ⁢ ( δ ⁢ U k ⁢ h ) ⇔ ∑ m = 1 M ∫ I m ( ∂ t δ ⁢ v k ⁢ h , Z k ⁢ h v ) − ( δ ⁢ p k ⁢ h , ∇ x ⋅ Z k ⁢ h v ) + ν ⁢ ( ∇ x δ ⁢ v k ⁢ h , ∇ x Z k ⁢ h v ) + ( ∇ x ⋅ δ ⁢ v k ⁢ h , Z k ⁢ h p ) ⁢ d ⁢ t + ∑ m = 1 M − 1 ( [ δ ⁢ v k ⁢ h ] m , Z k ⁢ h , m v , + ) + ( δ ⁢ v k ⁢ h , 0 + ⁢ ( 0 ) , Z k ⁢ h , 0 v , + ⁢ ( 0 ) ) = J ⁢ ( δ ⁢ U k ⁢ h ) .

To move the time derivative from the test function to the adjoint solution, we apply integration by parts in time and get

∑ m = 1 M ∫ I m ( δ ⁢ v k ⁢ h , − ∂ t Z k ⁢ h v ) − ( δ ⁢ p k ⁢ h , ∇ x ⋅ Z k ⁢ h v ) + ν ⁢ ( ∇ x δ ⁢ v k ⁢ h , ∇ x Z k ⁢ h v ) + ( ∇ x ⋅ δ ⁢ v k ⁢ h , Z k ⁢ h p ) ⁢ d ⁢ t − ∑ m = 1 M ( δ ⁢ v k ⁢ h , m − , [ Z k ⁢ h v ] m ) + ( δ ⁢ v k ⁢ h , M − ⁢ ( T ) , Z k ⁢ h , M v , − ⁢ ( T ) ) = J ⁢ ( δ ⁢ U k ⁢ h ) .

Similarly, to solve the adjoint problem of the Navier–Stokes equations, we now need to solve the general adjoint problem (3.2), which due to the linearity of our goal functional can be simplified to

A k ⁢ h , U ′ ⁢ ( U k ⁢ h ) ⁢ ( δ ⁢ U k ⁢ h , Z k ⁢ h ) = J ⁢ ( δ ⁢ U k ⁢ h ) .

Note that we still need to access the primal solution U k ⁢ h to be able to solve the adjoint equation. The adjoint problem for the Navier–Stokes equations reads

∑ m = 1 M ∫ I m ( δ ⁢ v k ⁢ h , − ∂ t Z k ⁢ h v ) − ( δ ⁢ p k ⁢ h , ∇ x ⋅ Z k ⁢ h v ) + ν ⁢ ( ∇ x δ ⁢ v k ⁢ h , ∇ x Z k ⁢ h v ) ⁢ d ⁢ t + ∑ m = 1 M ∫ I m ( ( δ ⁢ v k ⁢ h ⋅ ∇ x ) ⁢ v k ⁢ h + ( v k ⁢ h ⋅ ∇ x ) ⁢ δ ⁢ v k ⁢ h , Z k ⁢ h v ) + ( ∇ x ⋅ δ ⁢ v k ⁢ h , Z k ⁢ h p ) ⁢ d ⁢ t − ∑ m = 1 M ( δ ⁢ v k ⁢ h , m − , [ Z k ⁢ h v ] m ) + ( δ ⁢ v k ⁢ h , M − ⁢ ( T ) , Z k ⁢ h , M v , − ⁢ ( T ) ) = J ⁢ ( δ ⁢ U k ⁢ h ) .

3.4 PU-DWR Navier–Stokes Error Estimator

For linear PDEs, like the Stokes problem, and linear goal functionals, we have the error identity

J ( U ) − J ( U k ⁢ h ) = [ J ( U ) − J ( U k ) ] + [ J ( U k ) − J ( U k ⁢ h ) ] = − A k ⁢ ( U k ) ⁢ ( Z − Z k ) ⏟ = ⁣ : η k − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k − Z k ⁢ h ) ⏟ = ⁣ : η h = : η ,

which was derived in (3.3). For the nonlinear Navier–Stokes equations, we do not have the error identity (3.3), but need to use approximation (3.1) which contains the primal residual 𝜌 and the adjoint residual ρ ∗ . Including the adjoint residual yields an estimator which is accurate up to a remainder term of third order, compared to a purely primal residual based error estimator (see [14, Section 2] or [8, Chapter 6]) with a second-order remainder term only. For mildly nonlinear problems, such as incompressible Navier–Stokes (being semilinear in a laminar flow regime, i.e. nonlinearity in a low-order term) and linear goal functionals, the estimator asymptotically yields the same primal and adjoint residuals. As verification in the stationary setting, we adopted [25, Section 7.2] and evaluated 𝜌 and ρ ∗ separately, getting I eff primal ≈ I eff adjoint ≈ 1 . Clearly, for highly nonlinear problems (e.g. quasi-linear or fully nonlinear) and nonlinear goal functionals, both estimator parts are necessary as demonstrated in [25, Figure 4] and [24, Section 6.5]. Since being sufficiently accurate for such settings and being cheaper, the primal estimator only has been used in fluid applications; see e.g. [14, Section 8], [20] or [32]. Hence, for the estimation of the temporal and spatial error contributions, we employ

J ⁢ ( U ) − J ⁢ ( U k ) ≈ − A k ⁢ ( U k ) ⁢ ( Z − Z k ) , J ⁢ ( U k ) − J ⁢ ( U k ⁢ h ) ≈ − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( Z k − Z k ⁢ h ) .

To localize the spatial and temporal error, we now include the partition of unity { χ i m } i , m with ∑ m = 1 M ∑ i ∈ T h m χ i m ≡ 1 in the adjoint weights to localize the error. For the partition of unity, we choose linear finite elements in space ( cG ( 1 ) ) and piecewise discontinuous finite elements in time ( dG ( 0 ) ). This allows for the localization of the error to the individual temporal elements and the vertices of the spatial elements. We have

(3.4) η k = ∑ m = 1 M ∑ i ∈ T h m − A k ⁢ ( U k ) ⁢ ( ( Z − Z k ) ⁢ χ i m ) ⏟ = ⁣ : η k i , m , η h = ∑ m = 1 M ∑ i ∈ T h m − A k ⁢ h ⁢ ( U k ⁢ h ) ⁢ ( ( Z k − Z k ⁢ h ) ⁢ χ i m ) ⏟ = ⁣ : η h i , m .

Exemplarily, we will show the formula for the computation of the temporal error indicator η k i , m by applying the product rule

η k i , m = − A k ⁢ ( U k ) ⁢ ( ( Z − Z k ) ⁢ χ i m ) = − ∫ I m ( ∂ t v k , ( Z v − Z k v ) ⁢ χ i m ) − ( p k , ∇ x ⋅ { ( Z v − Z k v ) ⁢ χ i m } ) + ν ⁢ ( ∇ x v k , ∇ x { ( Z v − Z k v ) ⁢ χ i m } ) ⁢ d ⁢ t − ∫ I m ( ( v k ⋅ ∇ x ) ⁢ v k , ( Z v − Z k v ) ⁢ χ i m ) + ( ∇ x ⋅ v k , ( Z p − Z k p ) ⁢ χ i m ) ⁢ d ⁢ t − ( [ v k ] m − 1 , ( Z m − 1 v , + − Z k , m − 1 v , + ) ⁢ χ i m ) = − ∫ I m ( ∂ t v k , ( Z v − Z k v ) ⁢ χ i m ) − ( p k , { ∇ x ⋅ ( Z v − Z k v ) } ⁢ χ i m + ( Z v − Z k v ) ⋅ ( ∇ x χ i m ) ) ⁢ d ⁢ t − ∫ I m ν ⁢ ( ∇ x v k , { ∇ x ( Z v − Z k v ) } ⁢ χ i m + ( Z v − Z k v ) ⊗ ( ∇ x χ i m ) ) ⁢ d t − ∫ I m ( ( v k ⋅ ∇ x ) ⁢ v k , ( Z v − Z k v ) ⁢ χ i m ) + ( ∇ x ⋅ v k , ( Z p − Z k p ) ⁢ χ i m ) ⁢ d ⁢ t − ( [ v k ] m − 1 , ( Z m − 1 v , + − Z k , m − 1 v , + ) ⁢ χ i m ) .

Here, v ⊗ w stands for the tensor product of two vectors 𝒗 and 𝒘, i.e. ( v ⊗ w ) i , j := v i ⁢ w j and the initial condition is included as [ v k ] 0 := v + ⁢ ( 0 ) − v 0 . With this partition of unity approach, the spatial and temporal error can now be localized to each space-time degree of freedom of the cG ( 1 ) ⁢ dG ( 0 ) space-time discretization and will later indicate which spatial or temporal elements have the largest error and should be refined. However, to test the quality of our error estimator, we will initially work with globally refined space-time meshes and monitor how well the estimated error coincides with the true error by computing the effectivity index

I eff := | η J ⁢ ( U ) − J ⁢ ( U k ⁢ h ) | = | η k + η h J ⁢ ( U ) − J ⁢ ( U k ⁢ h ) | .

We asymptotically expect I eff → k , h → 0 1 for problems that satisfy a saturation assumption with a constant in the upper bound that converges to 0 for k , h → 0 (see [58]) (based on the ideas from the spatial case [25]).

3.5 Adaptive Algorithm

We will use the adaptive algorithm proposed in [58], which is a time-slabbing adaptation of the algorithm proposed in [55] to balance the temporal and the spatial discretization error.

Algorithm 1

Algorithm 1 (MARK and REFINE for the time-slabbing approach)

In our first numerical experiments, we chose a sufficiently big equilibration factor 𝑐 ( c = 10 7 ) to always refine both in space and time. In the future, we want to use a smaller equilibration constant 𝑐 to only refine either the temporal or the spatial discretization which dominates the total error. A computational analysis of the choice of the equilibration factor 𝑐 can be found in Appendix A.

After having computed the temporal error indicators for each temporal element and the spatial error indicators for each spatial element for each slab, we mark spatial and temporal elements for refinement according to Algorithm 1.

3.6 Evaluation of the Error Indicators

Until now, we have not discussed the practical realization of the PU-DWR error estimator (3.4). To make the error indicators computable, we need to make some approximations and get the unknown solutions by post-processing of the primal and adjoint solutions. In the works of Besier (and co-workers) [55, 54, 15], the linearization point in the temporal error has been replaced by the primal solution, i.e.

η k = − A k ⁢ ( U k ) ⁢ ( Z − Z k ) ≈ − A k ⁢ ( U k ⁢ h ) ⁢ ( Z − Z k ) ,

and he has discussed that this additional approximation error is small in comparison to the total discretization error [55, Remark 3.2]. Therefore, we will use this simplification in the same way. We still need to think about the approximation of the weights Z − Z k and Z k − Z k ⁢ h . This approximation depends on the choice of our finite element space for the adjoint problem, and in this work, we consider only a mixed-order error estimator.

For mixed order, the adjoint solution is of higher order in space and time than the primal solution. More concretely, we choose the primal solution to be Q 2 / Q 1 in space and dG ( 1 ) in time, whereas the adjoint solution is Q 4 / Q 2 in space and dG ( 2 ) in time. This setting has already been discussed in [12, Section 4], and in this work, we approximate the adjoint weights Z − Z k and Z k − Z k ⁢ h using the adjoint solution Z k ⁢ h ( 2 ) ∈ X k dG ( 2 ) ⁢ ( T k , V h 4 / 2 ⁢ ( T h ) ) by

Z − Z k ≈ ( id − I k dG ( 1 ) ) ⁢ Z k ⁢ h ( 2 ) , Z k − Z k ⁢ h ≈ ( id − I h 2 / 1 ) ⁢ I k dG ( 1 ) ⁢ Z k ⁢ h ( 2 ) .

Here, I h 2 / 1 : X ⁢ ( I , V h 4 / 2 ⁢ ( T h ) ) → X ⁢ ( I , V h 2 / 1 ⁢ ( T h ) ) denotes the interpolation in space from the enriched spatial function space to the primal spatial function space and I k dG ( 1 ) : X k dG ( 2 ) ⁢ ( T k , V ⁢ ( Ω ) ) → X k dG ( 1 ) ⁢ ( T k , V ⁢ ( Ω ) ) denotes the interpolation in time from the enriched temporal function space to the primal temporal function space.

From an implementational point of view, this approach is very convenient since it only requires interpolation of the adjoint solution to a lower-order function in space. This can be done by evaluating the adjoint solution at the space-time degrees of freedom of the lower-order space.

4 Numerical Tests

In this section, we perform some numerical experiments to verify the validity of our space-time error estimation framework. In the first two experiments, we consider two-dimensional time-dependent laminar flow around a cylinder, which is being modeled by the Stokes or Navier–Stokes equations. In a further experiment, we model the two-dimensional time-dependent laminar backward-facing step problem and do adaptive refinement using a nonlinear functional. The code is written in the C++ based FEM library deal.II [6, 5].

4.1 Stokes Flow around a Cylinder

In the first numerical experiment, we consider the test case 2D-3^[2] from the Schäfer–Turek Navier–Stokes benchmarks [53]. The 2D-3 benchmark problem describes the two-dimensional laminar flow around a circular cylinder with a time-dependent inflow condition. In Figure 2, the domain Ω and the different boundary components Γ in , Γ out , Γ wall and Γ circle are depicted.

Figure 2

Domain of the 2D Navier–Stokes benchmark problem.

Here, the kinematic viscosity is ν = 10 − 3 . The domain is defined as Ω := ( 0 , 2.2 ) × ( 0 , 0.41 ) ∖ B r ⁢ ( 0.2 , 0.2 ) with r = 0.05 , where B r ⁢ ( x , y ) denotes the ball around ( x , y ) with radius 𝑟 and we denote the time interval by I := ( 0 , 8 ) . On the boundary ∂ Ω , we prescribe Dirichlet and Neumann boundary conditions. We apply the no-slip condition, i.e. homogeneous Dirichlet boundary conditions v = 0 on Γ wall ∪ Γ circle , where

Γ wall = ( 0 , 2.2 ) × { 0 } ∪ ( 0 , 2.2 ) × { 0.41 } and Γ circle = ∂ B r ⁢ ( 0.2 , 0.2 ) .

On the boundary Γ in = { 0 } × ( 0 , 0.41 ) , we enforce a time-dependent parabolic inflow profile, which is a type of inhomogeneous Dirichlet boundary condition, i.e. v = v D on Γ in . For the 2D-3 benchmark, the inflow parabola is given by

v D ⁢ ( t , 0 , y ) = ( sin ⁡ ( π ⁢ t 8 ) ⁢ 6 ⁢ y ⁢ ( 0.41 − y ) 0.41 2 0 ) .

For brevity, we denote the Dirichlet boundary as

Γ D := Γ in ∪ Γ wall ∪ Γ circle .

The outflow boundary Γ out = { 2.2 } × ( 0 , 0.41 ) has Neumann boundary conditions. To get a correct outflow profile of the pressure, it has been discussed in [30] that the (modified) do-nothing condition ν ⁢ ( ∇ x v ) ⋅ n − p ⁢ n = 0 needs to be enforced on Γ out . For the non-symmetric stress tensor, this corresponds to homogeneous Neumann boundary conditions, i.e. σ ⋅ n = 0 on Γ out . It can be shown that, for this problem, the Reynolds number satisfies 0 ≤ Re ⁡ ( t ) ≤ 100 .

At first, we consider a simplification of this benchmark problem, where the nonlinearity is omitted such that we solve the Stokes equations. As mentioned in Remark 7, this results in an error identity for our primal error estimator and we expect I eff → 1 under uniform refinement of the space-time meshes. For our calculations, we use the mean drag goal functional

J drag ⁢ ( U ) := 1 8 ⁢ ∫ 0 8 c D ⁢ ( t ) ⁢ d t = 20 8 ⁢ ∫ 0 8 ∫ Γ circle σ ⁢ ( U ) ⋅ n ⋅ e 1 ⁢ d s ⁢ d t ≈ 0.40284197 ,

where e 1 = ( 1 , 0 ) T , c D denotes the drag coefficient and the reference value has been calculated on a fine spatio-temporal grid with 10,240 temporal and 1,480,000 spatial degrees of freedom.

4.1.1 Uniform Refinement

Using our mixed-order error estimator, we perform a grid search over the number of primal temporal degrees of freedom N k ∈ { 20 , 40 , 80 , 160 , 320 , 640 } and over the number of primal spatial degrees of freedom

N h ∈ { 445 , 1610 , 6100 , 23720 , 93520 } .

Furthermore, we compare Gauss–Legendre and Gauss–Lobatto quadrature for the temporal support points of the dG ( r ) method.

Gauss–Legendre Support Points in Time

We observe in Table 1 that the effectivity index for the Stokes test case is almost constant under temporal refinement. This can be explained by Table 2 and Table 3, where we find that the discretization error consists mainly of the spatial discretization error, while the temporal discretization error is a few orders of magnitude smaller for this test case.

Looking at the columns of Table 2, i.e. keeping the number of temporal degrees of freedom fixed and only refining uniformly in space, we see that the temporal error estimator remains almost constant. This verifies that our temporal error estimator does not depend on spatial refinement.

Looking at the rows of Table 3, i.e. keeping the number of spatial degrees of freedom fixed and only refining uniformly in time, we see that the spatial error remains almost constant. This verifies that our spatial error estimator does not depend on temporal refinement. Additionally, we observe quadratic convergence of the spatial error with respect to the spatial mesh size ℎ.

Table 1

Effectivity indices of mixed order for Stokes 2D-3 with Gauss–Legendre support points in time.

N h \ N k	20	40	80	160	320
445	0.61	0.61	0.61	0.61	0.61
1,610	0.68	0.68	0.68	0.68	0.68
6,100	0.73	0.72	0.72	0.72	0.72
23,720	0.77	0.75	0.75	0.75	0.75
93,520	0.90	0.81	0.80	0.80	0.80

Table 2

Temporal error estimator of mixed order for Stokes 2D-3 with Gauss–Legendre support points in time.

N h \ N k	20	40	80	160	320
445	− 2.7901 ⋅ 10 − 6	− 5.2875 ⋅ 10 − 7	− 1.0091 ⋅ 10 − 7	− 1.7218 ⋅ 10 − 8	− 2.5813 ⋅ 10 − 9
1,610	− 1.0153 ⋅ 10 − 5	− 1.5887 ⋅ 10 − 6	− 2.2938 ⋅ 10 − 7	− 3.1912 ⋅ 10 − 8	− 4.3840 ⋅ 10 − 9
6,100	− 1.1674 ⋅ 10 − 5	− 1.9094 ⋅ 10 − 6	− 2.8545 ⋅ 10 − 7	− 4.0021 ⋅ 10 − 8	− 5.3819 ⋅ 10 − 9
23,720	− 1.1797 ⋅ 10 − 5	− 1.9447 ⋅ 10 − 6	− 2.9447 ⋅ 10 − 7	− 4.1974 ⋅ 10 − 8	− 5.7309 ⋅ 10 − 9
93,520	− 1.1806 ⋅ 10 − 5	− 1.9471 ⋅ 10 − 6	− 2.9515 ⋅ 10 − 7	− 4.2178 ⋅ 10 − 8	− 5.7855 ⋅ 10 − 9

Table 3

Spatial error estimator of mixed order for Stokes 2D-3 with Gauss–Legendre support points in time.

N h \ N k	20	40	80	160	320
445	3.7394 ⋅ 10 − 2	3.7397 ⋅ 10 − 2	3.7398 ⋅ 10 − 2	3.7398 ⋅ 10 − 2	3.7398 ⋅ 10 − 2
1,610	1.3111 ⋅ 10 − 2	1.3110 ⋅ 10 − 2	1.3110 ⋅ 10 − 2	1.3110 ⋅ 10 − 2	1.3110 ⋅ 10 − 2
6,100	4.0683 ⋅ 10 − 3	4.0674 ⋅ 10 − 3	4.0673 ⋅ 10 − 3	4.0673 ⋅ 10 − 3	4.0673 ⋅ 10 − 3
23,720	1.1460 ⋅ 10 − 3	1.1457 ⋅ 10 − 3	1.1457 ⋅ 10 − 3	1.1457 ⋅ 10 − 3	1.1457 ⋅ 10 − 3
93,520	3.0495 ⋅ 10 − 4	3.0486 ⋅ 10 − 4	3.0485 ⋅ 10 − 4	3.0485 ⋅ 10 − 4	3.0484 ⋅ 10 − 4

Table 4

Effectivity indices of mixed order for Stokes 2D-3 with Gauss–Lobatto support points in time.

N h \ N k	20	40	80	160	320
445	0.58	0.61	0.61	0.61	0.61
1610	0.58	0.65	0.67	0.68	0.68
6100	0.45	0.63	0.70	0.72	0.72
23720	0.23	0.49	0.66	0.73	0.75
93520	0.08	0.25	0.52	0.70	0.77

Table 5

Temporal error estimator of mixed order for Stokes 2D-3 with Gauss–Lobatto support points in time.

N h \ N k	20	40	80	160	320
445	− 1.0450 ⋅ 10 − 3	− 3.3808 ⋅ 10 − 4	8.4236 ⋅ 10 − 5	2.6087 ⋅ 10 − 4	3.2106 ⋅ 10 − 4
1,610	− 2.4703 ⋅ 10 − 4	− 1.6773 ⋅ 10 − 4	− 8.1390 ⋅ 10 − 5	− 1.2743 ⋅ 10 − 5	2.5248 ⋅ 10 − 5
6,100	− 4.2596 ⋅ 10 − 5	− 2.9024 ⋅ 10 − 5	− 2.2523 ⋅ 10 − 5	− 1.5217 ⋅ 10 − 5	− 7.6783 ⋅ 10 − 6
23,720	− 1.7540 ⋅ 10 − 5	− 4.5460 ⋅ 10 − 6	− 2.5516 ⋅ 10 − 6	− 2.0765 ⋅ 10 − 6	− 1.6699 ⋅ 10 − 6
93,520	− 1.5636 ⋅ 10 − 5	− 2.4541 ⋅ 10 − 6	− 4.8653 ⋅ 10 − 7	− 2.0288 ⋅ 10 − 7	− 1.5701 ⋅ 10 − 7

Table 6

Spatial error estimator of mixed order for Stokes 2D-3 with Gauss–Lobatto support points in time.

N h \ N k	20	40	80	160	320
445	3.8122 ⋅ 10 − 2	3.7657 ⋅ 10 − 2	3.7294 ⋅ 10 − 2	3.7132 ⋅ 10 − 2	3.7076 ⋅ 10 − 2
1,610	1.3234 ⋅ 10 − 2	1.3249 ⋅ 10 − 2	1.3185 ⋅ 10 − 2	1.3121 ⋅ 10 − 2	1.3085 ⋅ 10 − 2
6,100	4.0609 ⋅ 10 − 3	4.0857 ⋅ 10 − 3	4.0874 ⋅ 10 − 3	4.0820 ⋅ 10 − 3	4.0749 ⋅ 10 − 3
23,720	1.1383 ⋅ 10 − 3	1.1456 ⋅ 10 − 3	1.1473 ⋅ 10 − 3	1.1476 ⋅ 10 − 3	1.1473 ⋅ 10 − 3
93,520	3.0247 ⋅ 10 − 4	3.0437 ⋅ 10 − 4	3.0485 ⋅ 10 − 4	3.0496 ⋅ 10 − 4	3.0499 ⋅ 10 − 4

Gauss–Lobatto Support Points in Time

In contrast to Table 1, we observe in Table 4 that the effectivity indices with Gauss–Lobatto support points in time are smaller for coarse temporal meshes. However, for fine temporal scales, the effectivity indices from Table 4 converge to the effectivity indices of Table 1.

In Table 5, we observe that the temporal error estimator with Gauss–Lobatto support points is larger than in Table 2 for Gauss–Legendre. Therefore, from now on, we only use Gauss–Legendre support points in time.

Finally, in Table 6, we get similar spatial error estimates as in Table 3. The only difference is that the spatial error estimator changes slightly less under temporal refinement for Gauss–Legendre support points in time, which further motivates our preferred usage of this quadrature formula.

4.1.2 Adaptive Refinement

We now apply the space-time adaptivity, Algorithm 1, starting from the coarsest spatial mesh with N h = 445 primal DoFs and M = 20 slabs with one temporal element each, i.e. we have N k = 40 initial temporal degrees of freedom due to the dG ( 1 ) discretization. This results in 17,800 space-time degrees of freedom.

As a marking strategy in space, we use the averaging strategy, which marks all elements K ∈ T h such that

η h | K > α | T h | ⁢ ∑ K ′ ∈ T h η h | K ′ .

In our computations, we used α = 1.1 . As a marking strategy in time, we use a fixed rate strategy (fixed number), where the ε % of temporal cells with the largest error are being refined. In our computations, we used ε = 75 .

Table 7

Adaptive refinement of mixed order on dynamic meshes for Stokes 2D-3.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 5.2875 ⋅ 10 − 7	3.7397 ⋅ 10 − 2	3.7397 ⋅ 10 − 2	6.0872 ⋅ 10 − 2	0.61
55,328	303,288	36	5.6502 ⋅ 10 − 8	1.5619 ⋅ 10 − 2	1.5619 ⋅ 10 − 2	2.1453 ⋅ 10 − 2	0.73
259,712	1,468,032	64	− 3.8869 ⋅ 10 − 8	4.1840 ⋅ 10 − 3	4.1840 ⋅ 10 − 3	5.7545 ⋅ 10 − 3	0.73
902,838	5,137,110	113	− 8.8754 ⋅ 10 − 9	1.1668 ⋅ 10 − 3	1.1668 ⋅ 10 − 3	1.5598 ⋅ 10 − 3	0.75
3,438,496	19,668,018	199	− 7.5194 ⋅ 10 − 9	3.1632 ⋅ 10 − 4	3.1632 ⋅ 10 − 4	4.2411 ⋅ 10 − 4	0.75

Table 8

Adaptive refinement of mixed order on dynamic meshes for Stokes 2D-3 with divergence-free L 2 projection.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 1.7268 ⋅ 10 − 6	3.7635 ⋅ 10 − 2	3.7634 ⋅ 10 − 2	5.9702 ⋅ 10 − 2	0.63
55,092	301,932	36	− 4.7644 ⋅ 10 − 8	1.5736 ⋅ 10 − 2	1.5736 ⋅ 10 − 2	2.1001 ⋅ 10 − 2	0.75
259,712	1,468,032	64	− 4.2556 ⋅ 10 − 8	4.2011 ⋅ 10 − 3	4.2011 ⋅ 10 − 3	5.3248 ⋅ 10 − 3	0.79
902,838	5,137,110	113	− 1.0787 ⋅ 10 − 8	1.1719 ⋅ 10 − 3	1.1718 ⋅ 10 − 3	1.3084 ⋅ 10 − 3	0.90
3,437,792	19,663,926	199	− 7.5608 ⋅ 10 − 9	3.1791 ⋅ 10 − 4	3.1790 ⋅ 10 − 4	2.7098 ⋅ 10 − 4	1.17

Table 9

Adaptive refinement of mixed order on dynamic meshes for Stokes 2D-3 with divergence-free H 0 1 projection.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 1.9364 ⋅ 10 − 6	3.7704 ⋅ 10 − 2	3.7702 ⋅ 10 − 2	5.9309 ⋅ 10 − 2	0.64
55,092	301,932	36	− 9.5524 ⋅ 10 − 8	1.5745 ⋅ 10 − 2	1.5745 ⋅ 10 − 2	2.0791 ⋅ 10 − 2	0.76
259,712	1,468,032	64	− 4.7248 ⋅ 10 − 8	4.1966 ⋅ 10 − 3	4.1966 ⋅ 10 − 3	5.0362 ⋅ 10 − 3	0.83
902,838	5,137,110	113	− 1.1303 ⋅ 10 − 8	1.1696 ⋅ 10 − 3	1.1696 ⋅ 10 − 3	1.1248 ⋅ 10 − 3	1.04
3,439,256	19,672,314	199	− 7.3748 ⋅ 10 − 9	3.1672 ⋅ 10 − 4	3.1672 ⋅ 10 − 4	1.5370 ⋅ 10 − 4	2.06

In Table 7, we show the results from five adaptive refinement loops with a mixed-order error estimator for the Stokes 2D-3 problem. We do not use a divergence-free projection here, which leads to oscillations in the goal functionals (cf. Figure 4). However, the error estimates for the Stokes problem are still acceptable and the effectivity index does not change much.

In Table 8, we have the same setup as in Table 7 and additionally employ a divergence-free L 2 projection from Section 2.5. We observe that ensuring the incompressibility of the snapshot on the new mesh yields a smaller true error in the goal functional between the reference solution and the finite element solution than in Table 7. In contrast to Table 8, we use a divergence-free H 0 1 projection in Table 9, but we make similar observations when it comes to the true error.

Figure 3

Comparison of the error in the mean drag for the Stokes problem for uniform and adaptive spatio-temporal refinement.

In Figure 3, we start with a common coarse spatio-temporal mesh and iteratively refine uniformly in space and time or adaptively refine in space and time based on Algorithm 1. As expected the dual-weighted residual method driven mesh refinement leads to faster convergence in the error in the mean drag functional. In particular, for adaptive refinement, a lot fewer degrees of freedom are needed to reach a prescribed accuracy in the goal functional, e.g. to achieve a tolerance of 2 ⋅ 10 − 3 , using uniform refinement requires 7,590,400 space-time degrees of freedom, whereas adaptive refinement only needs 902,838 space-time degrees of freedom for the primal problem.

Figure 4

Trajectory of the drag/lift coefficient for the Stokes problem for adaptive refinement.

(a)

Drag coefficient

(b)

Lift coefficient

In Figure 4, we compare the trajectory for the drag and lift coefficients over time for the Stokes 2D-3 problem and the different adaptive refinement methods. The drag trajectory is mostly smooth with only minor sawtooth-like oscillations near the maximal drag. However, the lift coefficient displays large spikes when using adaptive refinement without ensuring the incompressibility of the solution on the next spatial mesh. This underlines the necessity of methods like the divergence-free L 2 and H 0 1 projections from Besier and Wollner [16], which lead to smoother lift coefficient trajectories that more closely match the uniformly refined solution.

4.2 Navier–Stokes Flow around a Cylinder

We now include the nonlinear term in the variational formulation and solve the Navier–Stokes equations. In contrast to the Stokes problem from Section 4.1, there are plenty of reference results for the Navier–Stokes 2D-3 benchmark problem in the literature [53, 34], and we decided to use the publicly available results of the FeatFlow software^[3] of the group at the TU Dortmund around Stefan Turek. The FeatFlow authors discretize in space with Q 2 / P 1 disc elements, and in time, they discretize with the Crank–Nicolson scheme. For their finest calculations, they use a fixed time step size of k = 0.000625 and 667,264 degrees of freedom in space. We use their mean drag

J drag ⁢ ( U ) ≈ 1.6031368

as our reference solution. As marking strategies, we again used the averaging strategy in space and the fixer rate strategy in time with the same parameters as for the linear case ( α = 1.1 and ε = 75 ).

Table 10

Uniform refinement of mixed order on dynamic meshes for Navier–Stokes 2D-3.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
128,800	732,000	40	1.3424 ⋅ 10 − 3	2.8297 ⋅ 10 − 1	2.8431 ⋅ 10 − 1	2.9690 ⋅ 10 − 1	0.96
976,000	5,692,800	80	6.8970 ⋅ 10 − 2	1.3391 ⋅ 10 − 1	2.0288 ⋅ 10 − 1	1.3978 ⋅ 10 − 1	1.45
7,590,400	44,889,600	160	5.5834 ⋅ 10 − 3	2.4517 ⋅ 10 − 2	3.0100 ⋅ 10 − 2	2.5428 ⋅ 10 − 2	1.18

Table 11

Adaptive refinement of mixed order on dynamic meshes for Navier–Stokes 2D-3.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
63,336	350,244	36	2.3497 ⋅ 10 − 7	2.6516 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6719 ⋅ 10 − 1	0.99
228,568	1,285,878	64	1.7174 ⋅ 10 − 3	− 1.0803 ⋅ 10 − 2	9.0853 ⋅ 10 − 3	1.2438 ⋅ 10 − 1	0.07
835,942	4,749,933	113	− 1.3694 ⋅ 10 − 1	1.7656 ⋅ 10 − 1	3.9619 ⋅ 10 − 2	2.7110 ⋅ 10 − 2	1.46
3,008,930	17,274,966	199	3.9451 ⋅ 10 − 3	1.1779 ⋅ 10 − 2	1.5724 ⋅ 10 − 2	1.1074 ⋅ 10 − 2	1.42

Table 12

Adaptive refinement of mixed order on dynamic meshes for Navier–Stokes 2D-3 with divergence-free L 2 projection.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 9.3298 ⋅ 10 − 6	3.8542 ⋅ 10 − 1	3.8541 ⋅ 10 − 1	5.5752 ⋅ 10 − 1	0.69
63,454	350,922	36	− 1.4015 ⋅ 10 − 7	2.6505 ⋅ 10 − 1	2.6505 ⋅ 10 − 1	2.6642 ⋅ 10 − 1	0.99
230,032	1,294,482	64	8.9182 ⋅ 10 − 4	− 1.2571 ⋅ 10 − 2	1.1679 ⋅ 10 − 2	1.2586 ⋅ 10 − 1	0.09
828,744	4,706,883	113	− 1.1615 ⋅ 10 − 1	7.6888 ⋅ 10 − 2	3.9265 ⋅ 10 − 2	2.5449 ⋅ 10 − 2	1.54
3,004,686	17,251,722	199	4.3194 ⋅ 10 − 3	1.9094 ⋅ 10 − 2	2.3414 ⋅ 10 − 2	1.9674 ⋅ 10 − 2	1.19

Table 13

Adaptive refinement of mixed order on dynamic meshes for Navier–Stokes 2D-3 with divergence-free H 0 1 projection.

#DoF(primal)	#DoF(adjoint)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 9.1676 ⋅ 10 − 6	3.8548 ⋅ 10 − 1	3.8547 ⋅ 10 − 1	5.5768 ⋅ 10 − 1	0.69
63,690	352,278	36	1.4710 ⋅ 10 − 6	2.6493 ⋅ 10 − 1	2.6493 ⋅ 10 − 1	2.6667 ⋅ 10 − 1	0.99
234,878	1,322,502	64	7.5795 ⋅ 10 − 4	− 4.4232 ⋅ 10 − 3	3.6653 ⋅ 10 − 3	1.2633 ⋅ 10 − 1	0.03
834,710	4,741,881	113	− 2.3546 ⋅ 10 − 3	− 7.6972 ⋅ 10 − 2	7.9327 ⋅ 10 − 2	1.9651 ⋅ 10 − 2	4.04
3,044,708	17,485,449	199	3.0977 ⋅ 10 − 3	9.0900 ⋅ 10 − 3	1.2188 ⋅ 10 − 2	6.5227 ⋅ 10 − 3	1.87

In Table 10, we observe that the effectivity indices are close to 1 under uniform spatio-temporal refinement for the Navier–Stokes 2D-3 problem. We want to remark again that, for this nonlinear problem, our primal error estimator no longer satisfies the error identity, but we can computationally verify its accuracy.

In Table 11, we display the results from five DWR refinement loops with a mixed-order error estimator for the Navier–Stokes 2D-3 problem. Although we do not employ a divergence-free projection, the effectivity indices hint at a good match between the true error and the error estimator. However, in the third loop, a significant underestimation of the error in the mean drag can be observed which nevertheless does not persist after further refinement.

In Table 12, we observe that the accuracy of the error estimator is improved compared to Table 11 by employing a divergence-free L 2 projection. Additionally, this reduces the true error and thus yields a faster convergence toward the true solution.

In Table 13, we have the same setup as in Table 12 but now use a divergence-free H 0 1 projection. In comparison to Table 12, the effectivity indices are now slightly worse, but the true error is lower.

Figure 5

Comparison of the error in the mean drag for the Navier–Stokes problem for uniform and adaptive spatio-temporal refinement.

In Figure 5, we visualize the true error from Tables 10, 11, 12 and 13. Analogous to Figure 3, we detect that the error of our adaptively refined solution converges faster than a uniformly refined solution.

Figure 6

Trajectory of the drag/lift coefficient for the Navier–Stokes problem for adaptive refinement.

(a)

Drag coefficient

(b)

Lift coefficient

In Figure 6, we compare the trajectory of the drag and lift coefficient for the Navier–Stokes 2D-3 benchmark problem under different kinds of adaptive refinement. Like we observed in Figure 4, the divergence-free projections create less oscillatory drag and lift trajectories. This can be seen especially for the lift coefficient, where the projections help to remove the large spikes caused by using dynamic meshes. At the same time, this causes a slight damping between t = 4 ⁢ s and t = 6 ⁢ s .

$Figure 7 Snapshots of the velocity magnitude of the primal Navier–Stokes problem after four adaptive refinements with divergence-free L 2 L^{2} projection.$

Figure 7

Snapshots of the velocity magnitude of the primal Navier–Stokes problem after four adaptive refinements with divergence-free L 2 projection.

In Figure 7, we show some dynamic adaptive spatial meshes together with the snapshots of the velocity magnitude of the primal Navier–Stokes problem after four adaptive refinements with divergence-free L 2 projection. We observe that mainly the region around the cylinder is being refined, which is to be expected given that the drag coefficient is a boundary integral. Further away from the cylinder, e.g. close to the outflow boundary, the solution does not need to be resolved with the same accuracy. Moreover, since we are using dynamic meshes, the spatial mesh refinement can change over time as shown in Figure 7.

$Figure 8 Temporal mesh in the final DWR loop for Navier–Stokes 2D-3 with divergence-free L 2 L^{2} projection.$

Figure 8

Temporal mesh in the final DWR loop for Navier–Stokes 2D-3 with divergence-free L 2 projection.

In Figure 8, we show the temporal mesh in the fifth DWR loop for the Navier–Stokes 2D-3 benchmark when using a divergence-free L 2 projection. The temporal mesh is also being refined adaptively, and we can see that, in some regions, the diameter of the temporal elements is smaller, especially in the second half of the temporal domain.

4.3 Navier–Stokes Flow over Backward Facing Step

In the third numerical experiment, we consider a version of the two-dimensional laminar backward-facing step benchmark [4, 17] with a time-dependent inflow condition. The domain is shown in Figure 9.

Figure 9

Domain of the 2D backward-facing step benchmark.

The domain is defined as Ω := ( − 1 2 , 0 ) × ( 1 2 , 1 ) ∪ ( 0 , 2 ) × ( 0 , 1 ) . The kinematic viscosity is ν = 1 , and we consider the time interval I = ( 0 , 8 ) . We prescribe homogeneous Neumann boundary conditions on the outflow boundary Γ out := { 2 } × ( 0 , 1 ) . On the inflow boundary Γ in := { − 1 2 } × ( 1 2 , 1 ) , we enforce the time-dependent parabolic inflow profile

v D ⁢ ( t , − 1 2 , y ) = ( sin ⁡ ( π ⁢ t 8 ) ⁢ 80 ⁢ ( − 1 + ( 3 − 2 ⁢ y ) ⁢ y ) 0 ) .

On all remaining boundaries, homogeneous Dirichlet boundary conditions are being applied. In contrast to the previous test cases, we now consider the mean squared vorticity

J vorticity ⁢ ( U ) = 1 8 ⁢ ∫ 0 8 ∫ Ω ( ∇ × v ) 2 ⁢ d x ⁢ d t = 1 8 ⁢ ∫ 0 8 ∫ Ω ( ∂ 1 v 2 − ∂ 2 v 1 ) 2 ⁢ d x ⁢ d t

as our goal functional. Note that this is a nonlinear goal functional and the right-hand side of the dual problem is given by

J vorticity , U ′ ⁢ ( U ) ⁢ ( δ ⁢ U ) = 2 8 ⁢ ∫ 0 8 ∫ Ω ( ∇ × v ) ⋅ ( ∇ × δ ⁢ v ) ⁢ d x ⁢ d t .

Using 80 initial temporal degrees of freedom and a spatial grid with 1,439 primal degrees of freedom, we obtained the series of functional values from Table 14. We can see that this series converges from below and that the accuracy is increasing with each refinement step.

Table 14

Mean vorticity value under uniform spatio-temporal refinement for the backward-facing step test case.

#temporal DoF(primal)	#spatial DoF(primal)	J ⁢ ( U k ⁢ h )
80	1,439	506.20900857749609
160	5,467	510.27887058685025
320	21,299	512.24259375586087
640	84,067	513.18333123881553
1,280	334,019	513.63122326525161
2,560	1,331,587	513.84343972465513

For the sake of brevity, we will only discuss adaptive results in detail. For this test case, we used the fixed rate marking strategy both in space and time with ε = 30 .

Table 15

Signed effectivity indices for the different projections with a once coarsened spatio-temporal mesh (coarse) compared to the finest reference mesh (fine). Positive I eff means the computed functional value is below the reference value.

DWR loop	Coarse I eff	Fine I eff	Coarse I eff L 2	Fine I eff L 2	Coarse I eff H 0 1	Fine I eff H 0 1
1	1.06	1.03	1.05	1.03	1.06	1.03
2	1.15	1.08	1.15	1.08	1.10	1.04
3	1.34	1.17	1.34	1.17	1.83	1.53
4	2.02	1.37	2.21	1.47	−2.51	−3.97
5	−1.11	2.06	−1.53	2.47	−6.46	−0.70

Table 15 shows the resulting effectivity indices for the different projection strategies (no projection, L 2 projection, H 0 1 projection). To investigate the effect of the reference value on the approximated real error, we compare the results with the last (fine) and second-to-last (coarse) mean vorticity as reference values J vorticity ⁢ ( U ) . As long as the vorticity on the adaptive mesh is below the reference value, we can see that we get a better effectivity with the finer reference value, suggesting that our estimator is a good predictor of the true error. Additionally, we see a bad effectivity index for the loop in which we have the same h min as in the reference grid. However, this is not surprising as it is generally not advisable to analyze adaptive results on the reference mesh level.

To overcome the issue of the unreliable reference solution and since the functional values in Table 14 converge nicely, we extrapolated a new value from the two finest values using the formula in [3, Section 17.7.5] with h 2 = h 1 2 , obtaining

J extr. := ( J 2 − 1 4 ⁢ J 1 ) ⁢ 4 3 = 513.9141785444563 .

Additionally, we analyzed the differences between consecutive functional values, i.e. d ⁢ J n := J n − J n − 1 , where J n is the functional value at grid level 𝑛 starting with J 0 = 506.209 . As these differences are roughly halved in each uniform refinement step, we applied linear regression to the five log ⁡ ( J n ) values and obtained

log ⁡ ( d ⁢ J ⁢ ( n ) ) ≈ 2.148469873234862 − 0.738556050505728 ⋅ n .

With this, we calculated new functional values by

J reg. ⁢ ( n ) = J 0 + ∑ k = 1 n exp ⁡ ( log ⁡ ( d ⁢ J ⁢ ( k ) ) ) .

Using this formula, we obtained two new functional values

J reg. ⁢ ( 6 ) = 513.9587146059994 and J reg. ⁢ ( 25 ) = 514.0520354056016 .

Table 16

Signed effectivity indices for the different projections with an extrapolated reference value (extr.) vs. reference values J reg. ⁢ ( 6 ) and J reg. ⁢ ( 25 ) obtained by regression of the log-differenced function values. Positive I eff means the computed functional value is below the reference value.

Loop	Extr. I eff	Reg. 6 I eff	Reg. 25 I eff	Extr. I eff L 2	Reg. 6 I eff L 2	Reg. 25 I eff L 2	Extr. I eff H 0 1	Reg. 6 I eff H 0 1	Reg. 25 I eff H 0 1
1	1.02	1.02	1.00	1.01	1.01	1.00	1.02	1.02	1.00
2	1.06	1.05	1.02	1.06	1.04	1.02	1.02	1.01	0.98
3	1.12	1.09	1.03	1.12	1.09	1.03	1.45	1.40	1.31
4	1.23	1.16	1.04	1.32	1.24	1.11	−4.94	−5.81	−9.30
5	1.54	1.33	1.03	1.78	1.51	1.15	−0.72	−0.73	−0.76

Table 16 shows that extrapolation of the reference value greatly improves the effectivity indices. In particular, recursively applying the extrapolation formula a few times (cf. reg. 25) yields effectivity indices much closer to 1 when using no projection or an L 2 projection. However, the H 0 1 projection results still do not look good on the finer meshes.

$Figure 10 Comparison of the true error in the mean squared vorticity for the backward-facing step for uniform and adaptive spatio-temporal refinement using J reg. ⁢ ( 25 ) J_{\textup{reg.}}(25) as the reference value.$

Figure 10

Comparison of the true error in the mean squared vorticity for the backward-facing step for uniform and adaptive spatio-temporal refinement using J reg. ⁢ ( 25 ) as the reference value.

In Figure 10, we visualize the true error in the mean squared vorticity for the backward-facing step for uniform and adaptive spatio-temporal refinement. Analogous to Figure 3 and Figure 5 for the (Navier–)Stokes 2D-3 benchmark, we detect that the error of our adaptively refined solution converges with a higher rate than a uniformly refined solution. In particular, when comparing the fourth refinement cycle, we observe that adaptive refinement uses only 1 6 of the number of space-time degrees of freedom as uniform refinement to reach a tolerance of 0.7 in the goal functional.

$Figure 11 Snapshot of the velocity magnitude and streamlines of the primal backward-facing step problem at t ≈ 4 t\approx 4 .$

Figure 11

Snapshot of the velocity magnitude and streamlines of the primal backward-facing step problem at t ≈ 4 .

$Figure 12 Snapshot of the vorticity magnitude of the primal backward-facing step problem at t ≈ 4 t\approx 4 .$

Figure 12

Snapshot of the vorticity magnitude of the primal backward-facing step problem at t ≈ 4 .

$Figure 13 Spatial grid at t ≈ 4 t\approx 4 for the backward-facing step problem after 4 adaptive refinements with divergence-free L 2 L^{2} projection.$

Figure 13

Spatial grid at t ≈ 4 for the backward-facing step problem after 4 adaptive refinements with divergence-free L 2 projection.

In Figure 11, we show the snapshot of the velocity magnitude of the primal backward-facing step problem at time t ≈ 4 and its streamlines. Looking at the streamlines, we observe a corner vortex in the lower left corner after the step. In Figure 12, we plot the snapshot of the vorticity magnitude of the primal backward-facing step problem at time t ≈ 4 . It can be seen that the vorticity magnitude is largest at the reentrant corner x = ( 0 , 1 2 ) .

In Figure 13, we display a four times adaptively refined grid with divergence-free L 2 projection. On the one hand, we observe that the grid is mostly being refined close to the reentrant corner where the vorticity magnitude is at its maximum. On the other hand, closer to the inflow or outflow boundaries, the mesh is much coarser.

5 Conclusions and Outlook

In this work, we formulated a space-time PU-DWR approach for goal-oriented error control and local mesh adaptivity for the incompressible Navier–Stokes equations. Therein, our focus was temporal and spatial error control, which has been computationally analyzed for different benchmark problems such as the Schäfer–Turek 1996 benchmark. In terms of convergence as well as evaluation of effectivity indices, we observed excellent performances for different goal functionals. In ongoing work, an interesting aspect would be to implement Linke’s pressure-robust discretizations rather than the global corrections proposed by Besier and Wollner.

Funding source: Deutsche Forschungsgemeinschaft

Award Identifier / Grant number: 433082294

Funding statement: The first and fourth authors acknowledge the funding of the German Research Foundation (DFG; http://dx.doi.org/10.13039/501100001659) within the framework of the International Research Training Group on Computational Mechanics Techniques in High Dimensions GRK 2657 under Grant Number 433082294.

A Computational Tests of the Equilibration Factor in Algorithm 1

A.1 Changing Equilibration Factor with a Fixed Number of Initial Temporal Elements

To test the impact of the equilibration constant 𝑐 on the adaptive results, we repeated the simulation from Table 11 – the adaptive refinement of Navier–Stokes 2D-3 without divergence-free projection – with

c ∈ { 1 , 5 , 10 3 , 10 5 , 10 7 } .

As shown in these tables, we need to increase the equilibration constant 𝑐 to refine both the spatial and temporal meshes. In particular, for c = 1 (cf. Table 17) and c = 5 (cf. Table 18), the temporal mesh is not being refined at all since the spatial error dominates the temporal error. If we then increase the equilibration constant to c = 10 3 (cf. Table 19) or c = 10 5 (cf. Table 20), then we get some temporal refinement, but not in every refinement cycle. This motivates larger values for the equilibration constant, e.g. c = 10 7 (cf. Table 21 and Table 11), to refine both the spatial and temporal meshes.

Table 17

Adaptive refinement of Navier–Stokes 2D-3 with c = 1 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
35,428	195,978	20	− 3.7918 ⋅ 10 − 5	2.6520 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6716 ⋅ 10 − 1	0.99
72,440	407,643	20	− 3.9093 ⋅ 10 − 4	1.4886 ⋅ 10 − 2	1.4495 ⋅ 10 − 2	1.2598 ⋅ 10 − 1	0.12
143,692	815,421	20	− 1.5155 ⋅ 10 − 4	2.2525 ⋅ 10 − 2	2.2374 ⋅ 10 − 2	5.1573 ⋅ 10 − 2	0.43
285,260	1,633,776	20	− 2.1708 ⋅ 10 − 5	3.1083 ⋅ 10 − 3	3.0866 ⋅ 10 − 3	3.1490 ⋅ 10 − 2	0.10

Table 18

Adaptive refinement of Navier–Stokes 2D-3 with c = 5 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
35,428	195,978	20	− 3.7918 ⋅ 10 − 5	2.6520 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6716 ⋅ 10 − 1	0.99
72,440	407,643	20	− 3.9093 ⋅ 10 − 4	1.4886 ⋅ 10 − 2	1.4495 ⋅ 10 − 2	1.2598 ⋅ 10 − 1	0.12
143,692	815,421	20	− 1.5155 ⋅ 10 − 4	2.2525 ⋅ 10 − 2	2.2374 ⋅ 10 − 2	5.1573 ⋅ 10 − 2	0.43
285,260	1,633,776	20	− 2.1708 ⋅ 10 − 5	3.1083 ⋅ 10 − 3	3.0866 ⋅ 10 − 3	3.1490 ⋅ 10 − 2	0.10

Table 19

Adaptive refinement of Navier–Stokes 2D-3 with c = 10 3 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
35,428	195,978	20	− 3.7918 ⋅ 10 − 5	2.6520 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6716 ⋅ 10 − 1	0.99
72,440	407,643	20	− 3.9093 ⋅ 10 − 4	1.4886 ⋅ 10 − 2	1.4495 ⋅ 10 − 2	1.2598 ⋅ 10 − 1	0.12
263,200	1,494,210	36	− 2.7724 ⋅ 10 − 1	2.4822 ⋅ 10 − 1	2.9025 ⋅ 10 − 2	5.1403 ⋅ 10 − 2	0.56
995,472	5,718,699	64	1.4090 ⋅ 10 − 1	− 2.5893 ⋅ 10 − 1	1.1804 ⋅ 10 − 1	3.8993 ⋅ 10 − 2	3.03

Table 20

Adaptive refinement of Navier–Stokes 2D-3 with c = 10 5 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
63,336	350,244	36	2.3497 ⋅ 10 − 7	2.6516 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6719 ⋅ 10 − 1	0.99
129,164	726,756	36	4.0473 ⋅ 10 − 4	− 4.1378 ⋅ 10 − 2	4.0973 ⋅ 10 − 2	1.2605 ⋅ 10 − 1	0.33
468,778	2662221	64	− 1.7039 ⋅ 10 + 0	2.8312 ⋅ 10 + 0	1.1274 ⋅ 10 + 0	5.2671 ⋅ 10 − 2	21.4
1,719,320	9,873,924	113	2.4233 ⋅ 10 − 2	− 7.5433 ⋅ 10 − 4	2.3479 ⋅ 10 − 2	2.2151 ⋅ 10 − 2	1.06

Table 21

Adaptive refinement of Navier–Stokes 2D-3 with c = 10 7 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
17,800	96,600	20	− 8.4213 ⋅ 10 − 6	3.8557 ⋅ 10 − 1	3.8557 ⋅ 10 − 1	5.5869 ⋅ 10 − 1	0.69
63,336	350,244	36	2.3497 ⋅ 10 − 7	2.6516 ⋅ 10 − 1	2.6516 ⋅ 10 − 1	2.6719 ⋅ 10 − 1	0.99
228,568	1,285,878	64	1.7174 ⋅ 10 − 3	− 1.0803 ⋅ 10 − 2	9.0853 ⋅ 10 − 3	1.2438 ⋅ 10 − 1	0.07
835,942	4,749,933	113	− 1.3694 ⋅ 10 − 1	1.7656 ⋅ 10 − 1	3.9619 ⋅ 10 − 2	2.7110 ⋅ 10 − 2	1.46
3,008,930	17,274,966	199	3.9451 ⋅ 10 − 3	1.1779 ⋅ 10 − 2	1.5724 ⋅ 10 − 2	1.1074 ⋅ 10 − 2	1.42

A.2 Changing Number of Initial Temporal Elements with Fixed Equilibration Factor

We deliberately started with a very coarse temporal mesh, which might lead to an underestimation of the temporal error. Therefore, we performed additional tests with finer temporal meshes and c = 1 as well as c = 5 to see whether the equilibration plays a more important role when the error is better resolved in time. Overall, the behavior is similar to M = 20 (cf. Table 17 and Table 18) with most steps only refining the spatial grids. For M = 40 (cf. Table 22 and Table 23) and M = 80 (cf. Table 24 and Table 25), the last step is different with only temporal refinement for c = 1 and temporal as well as spatial refinement for c = 5 . Comparing actual errors as well as the number of DoFs, we can see that sufficient temporal refinement (roughly more than 120 intervals; see Tables 24, 25, 26, 27) is needed to get an error of around 0.01 in the fourth refinement step. So both enforced global refinement as well as equilibration with a finer initial grid led to good error reduction. However, while the number of degrees of freedom in the final step is smaller with equilibration, the total cost is higher in comparison since solving and estimating have to be performed for more time steps. Additionally, the effectivity indices for enforced global refinement seem to stabilize in the last two steps.

Table 22

Adaptive refinement of Navier–Stokes 2D-3 with c = 1 and M = 40 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
35,600	193,200	40	3.3813 ⋅ 10 − 8	3.8553 ⋅ 10 − 1	3.8553 ⋅ 10 − 1	5.5866 ⋅ 10 − 1	0.69
70,880	392,145	40	2.5936 ⋅ 10 − 6	2.6511 ⋅ 10 − 1	2.6511 ⋅ 10 − 1	2.6762 ⋅ 10 − 1	0.99
147,986	833,289	40	2.3636 ⋅ 10 − 4	− 3.8774 ⋅ 10 − 2	3.8538 ⋅ 10 − 2	1.2625 ⋅ 10 − 1	0.31
292,392	1661,022	40	3.1185 ⋅ 10 − 1	− 1.0598 ⋅ 10 − 1	2.0587 ⋅ 10 − 1	5.5182 ⋅ 10 − 2	3.73
513,362	2915,817	71	− 2.8892 ⋅ 10 − 1	2.3645 ⋅ 10 − 1	5.2464 ⋅ 10 − 2	3.4048 ⋅ 10 − 2	1.54

Table 23

Adaptive refinement of Navier–Stokes 2D-3 with c = 5 and M = 40 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
35,600	193,200	40	3.3813 ⋅ 10 − 8	3.8553 ⋅ 10 − 1	3.8553 ⋅ 10 − 1	5.5866 ⋅ 10 − 1	0.69
70,880	392,145	40	2.5936 ⋅ 10 − 6	2.6511 ⋅ 10 − 1	2.6511 ⋅ 10 − 1	2.6762 ⋅ 10 − 1	0.99
147,986	833,289	40	2.3636 ⋅ 10 − 4	− 3.8774 ⋅ 10 − 2	3.8538 ⋅ 10 − 2	1.2625 ⋅ 10 − 1	0.31
292,392	1,661,022	40	3.1185 ⋅ 10 − 1	− 1.0598 ⋅ 10 − 1	2.0587 ⋅ 10 − 1	5.5182 ⋅ 10 − 2	3.73
1,092,116	6,271,314	71	8.2739 ⋅ 10 − 2	− 6.5925 ⋅ 10 − 2	1.6814 ⋅ 10 − 2	3.4252 ⋅ 10 − 2	0.49

Table 24

Adaptive refinement of Navier–Stokes 2D-3 with c = 1 and M = 80 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
71,200	386,400	80	4.9319 ⋅ 10 − 7	3.8556 ⋅ 10 − 1	3.8556 ⋅ 10 − 1	5.5865 ⋅ 10 − 1	0.69
141,360	782,082	80	2.4853 ⋅ 10 − 6	2.6658 ⋅ 10 − 1	2.6658 ⋅ 10 − 1	2.6925 ⋅ 10 − 1	0.99
295,970	1,666,374	80	6.3941 ⋅ 10 − 3	8.3012 ⋅ 10 − 2	8.9406 ⋅ 10 − 2	1.2361 ⋅ 10 − 1	0.72
602,996	3,427,926	80	− 2.4988 ⋅ 10 − 1	1.3495 ⋅ 10 − 1	1.1493 ⋅ 10 − 1	3.1269 ⋅ 10 − 2	3.68
1,053,362	5,987,769	141	2.9715 ⋅ 10 − 3	− 2.0812 ⋅ 10 − 2	1.7841 ⋅ 10 − 2	1.2436 ⋅ 10 − 2	1.43

Table 25

Adaptive refinement of Navier–Stokes 2D-3 with c = 5 and M = 80 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
71,200	386,400	80	4.9319 ⋅ 10 − 7	3.8556 ⋅ 10 − 1	3.8556 ⋅ 10 − 1	5.5865 ⋅ 10 − 1	0.69
141,360	782,082	80	2.4853 ⋅ 10 − 6	2.6658 ⋅ 10 − 1	2.6658 ⋅ 10 − 1	2.6925 ⋅ 10 − 1	0.99
295,970	1,666,374	80	6.3941 ⋅ 10 − 3	8.3012 ⋅ 10 − 2	8.9406 ⋅ 10 − 2	1.2361 ⋅ 10 − 1	0.72
602,996	3,427,926	80	− 2.4988 ⋅ 10 − 1	1.3495 ⋅ 10 − 1	1.1493 ⋅ 10 − 1	3.1269 ⋅ 10 − 2	3.68
2,160,878	12,413,085	141	6.3265 ⋅ 10 − 3	6.4954 ⋅ 10 − 3	1.2822 ⋅ 10 − 2	9.8795 ⋅ 10 − 3	1.30

Table 26

Adaptive refinement of Navier–Stokes 2D-3 with c = 1 and M = 160 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
142,400	772,800	160	1.2458 ⋅ 10 − 7	3.8558 ⋅ 10 − 1	3.8558 ⋅ 10 − 1	5.5865 ⋅ 10 − 1	0.69
282,110	1,560,498	160	1.2875 ⋅ 10 − 6	2.6660 ⋅ 10 − 1	2.6660 ⋅ 10 − 1	2.6935 ⋅ 10 − 1	0.99
589,760	3,320,313	160	− 8.6094 ⋅ 10 − 5	1.4437 ⋅ 10 − 1	1.4428 ⋅ 10 − 1	1.1088 ⋅ 10 − 1	1.30
1,221,266	6,944,181	160	3.6799 ⋅ 10 − 3	− 2.8981 ⋅ 10 − 2	2.5301 ⋅ 10 − 2	4.8499 ⋅ 10 − 3	5.22
2,439,656	14,009,478	160	4.1653 ⋅ 10 − 3	− 9.1947 ⋅ 10 − 3	5.0294 ⋅ 10 − 3	1.2499 ⋅ 10 − 2	0.40

Table 27

Adaptive refinement of Navier–Stokes 2D-3 with c = 5 and M = 160 .

#DoF(primal)	#DoF(dual)	𝑀	η k	η h	𝜂	J ⁢ ( U ) − J ⁢ ( U k ⁢ h )	I eff
142,400	772,800	160	1.2458 ⋅ 10 − 7	3.8558 ⋅ 10 − 1	3.8558 ⋅ 10 − 1	5.5865 ⋅ 10 − 1	0.69
282,110	1,560,498	160	1.2875 ⋅ 10 − 6	2.6660 ⋅ 10 − 1	2.6660 ⋅ 10 − 1	2.6935 ⋅ 10 − 1	0.99
589,760	3,320,313	160	− 8.6094 ⋅ 10 − 5	1.4437 ⋅ 10 − 1	1.4428 ⋅ 10 − 1	1.1088 ⋅ 10 − 1	1.30
1,221,266	6,944,181	160	3.6799 ⋅ 10 − 3	− 2.8981 ⋅ 10 − 2	2.5301 ⋅ 10 − 2	4.8499 ⋅ 10 − 3	5.22
2,439,656	14,009,478	160	4.1653 ⋅ 10 − 3	− 9.1947 ⋅ 10 − 3	5.0294 ⋅ 10 − 3	1.2499 ⋅ 10 − 2	0.40

Acknowledgements

We thank Hendrik Fischer (Hannover/Paris), Marius Paul Bruchhäuser (Hamburg), Bernhard Endtmayer (Hannover), Johannes Lankeit (Hannover) and Ludovic Chamoin (Paris) for fruitful discussions and comments.

References

[1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Pure Appl. Math. (New York), Wiley-Interscience, New York, 2000. 10.1002/9781118032824Search in Google Scholar

[2] F. Alauzet and M. Mehrenberger, P 1 -conservative solution interpolation on unstructured triangular meshes, Internat. J. Numer. Methods Engrg. 84 (2010), no. 13, 1552–1588. 10.1002/nme.2951Search in Google Scholar

[3] S. Amstutz and T. Wick, Refresher course in maths and a project on numerical modeling done in twos, Institutionelles Repositorium der Leibniz Universität Hannover, Universität Hannover, Hannover (2022), https://doi.org/10.15488/11629. Search in Google Scholar

[4] B. Armaly, F. Durst, Jose Pereira and B. Schönung, Experimental and theoretical investigation of backward-facing step flow, J. Fluid Mech. 127 (1983), 473–496. 10.1017/S0022112083002839Search in Google Scholar

[5] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin and D. Wells, The deal.II finite element library: Design, features, and insights, Comput. Math. Appl. 81 (2021), 407–422. 10.1016/j.camwa.2020.02.022Search in Google Scholar

[6] D. Arndt, W. Bangerth, M. Feder, M. Fehling, R. Gassmöller, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, S. Sticko, B. Turcksin and D. Wells, The deal.II library, version 9.4., J. Numer. Math. 30 (2022), no. 3, 231–246. 10.1515/jnma-2022-0054Search in Google Scholar

[7] W. Bangerth, M. Geiger and R. Rannacher, Adaptive Galerkin finite element methods for the wave equation, Comput. Methods Appl. Math. 10 (2010), no. 1, 3–48. 10.2478/cmam-2010-0001Search in Google Scholar

[8] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Lectures in Math. ETH Zürich, Birkhäuser, Basel, 2003. 10.1007/978-3-0348-7605-6Search in Google Scholar

[9] E. Bänsch, F. Karakatsani and C. G. Makridakis, A posteriori error estimates for fully discrete schemes for the time dependent Stokes problem, Calcolo 55 (2018), no. 2, Paper No. 19. 10.1007/s10092-018-0259-2Search in Google Scholar

[10] S. Basava, K. Mang, M. Walloth, T. Wick and W. Wollner, Adaptive and pressure-robust discretization of incompressible pressure-driven phase-field fracture, Non-Standard Discretisation Methods in Solid Mechanics, Lect. Notes Appl. Comput. Mech. 98, Springer, Cham, (2022), 191–215. 10.1007/978-3-030-92672-4_8Search in Google Scholar

[11] S. R. Basava and W. Wollner, Gradient robust mixed methods for nearly incompressible elasticity, J. Sci. Comput. 95 (2023), Article No. 93. 10.1007/s10915-023-02227-0Search in Google Scholar

[12] 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. 91 (2021), 17–35. 10.1016/j.camwa.2020.08.028Search in Google Scholar

[13] 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. 4 (1996), no. 4, 237–264. Search in Google Scholar

[14] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer. 10 (2001), 1–102. 10.1017/S0962492901000010Search in Google Scholar

[15] M. Besier and R. Rannacher, Goal-oriented space-time adaptivity in the finite element Galerkin method for the computation of nonstationary incompressible flow, Internat. J. Numer. Methods Fluids 70 (2012), no. 9, 1139–1166. 10.1002/fld.2735Search in Google Scholar

[16] M. Besier and W. Wollner, On the pressure approximation in nonstationary incompressible flow simulations on dynamically varying spatial meshes, Internat. J. Numer. Methods Fluids 69 (2012), no. 6, 1045–1064. 10.1002/fld.2625Search in Google Scholar

[17] G. Biswas, M. Breuer and F. Durst, Backward-facing step flows for various expansion ratios at low and moderate Reynolds numbers, J. Fluids Eng.-Trans. ASME 126 (2004), 362–374. 10.1115/1.1760532Search in Google Scholar

[18] M. Braack and A. Ern, A posteriori control of modeling errors and discretization errors, Multiscale Model. Simul. 1 (2003), no. 2, 221–238. 10.1137/S1540345902410482Search in Google Scholar

[19] M. Braack, J. Lang and N. Taschenberger, Stabilized finite elements for transient flow problems on varying spatial meshes, Comput. Methods Appl. Mech. Engrg. 253 (2013), 106–116. 10.1016/j.cma.2012.08.006Search in Google Scholar

[20] M. Braack and T. Richter, Solutions of 3D Navier–Stokes benchmark problems with adaptive finite elements, Comput. Fluids 35 (2006), no. 4, 372–392. 10.1016/j.compfluid.2005.02.001Search in Google Scholar

[21] M. P. Bruchhäuser, U. Köcher and M. Bause, On the implementation of an adaptive multirate framework for coupled transport and flow, J. Sci. Comput. 93 (2022), no. 3, Paper No. 59. 10.1007/s10915-022-02026-zSearch in Google Scholar

[22] R. Dautray and J.-L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology. Vol. 5, Springer, Berlin, 2000. Search in Google Scholar

[23] J. D. De Basabe, M. K. Sen and M. F. Wheeler, The interior penalty discontinuous Galerkin method for elastic wave propagation: Grid dispersion, Geophys. J. Internat. 175 (2008), no. 1, 83–93. 10.1111/j.1365-246X.2008.03915.xSearch in Google Scholar

[24] B. Endtmayer, U. Langer and T. Wick, Multigoal-oriented error estimates for non-linear problems, J. Numer. Math. 27 (2019), no. 4, 215–236. 10.1515/jnma-2018-0038Search in Google Scholar

[25] B. Endtmayer, U. Langer and T. Wick, Two-side a posteriori error estimates for the dual-weighted residual method, SIAM J. Sci. Comput. 42 (2020), no. 1, A371–A394. 10.1137/18M1227275Search in Google Scholar

[26] L. Failer and T. Wick, Adaptive time-step control for nonlinear fluid-structure interaction, J. Comput. Phys. 366 (2018), 448–477. 10.1016/j.jcp.2018.04.021Search in Google Scholar

[27] G. P. Galdi, An Introduction to the Mathematical Theory of the Navier–Stokes Equations: Steady-State Problems, Springer Monogr. Math., Springer, New York, 2011. 10.1007/978-0-387-09620-9Search in Google Scholar

[28] V. Girault and P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Springer Ser. Comput. Math. 5, Springer, Berlin, 1986. 10.1007/978-3-642-61623-5Search in Google Scholar

[29] R. Glowinski, Finite element methods for incompressible viscous flow, Numerical Methods for Fluids (Part 3), Handb. Numer. Anal. 9, North-Holland, Amsterdam (2003), 3–1176. 10.1016/S1570-8659(03)09003-3Search in Google Scholar

[30] J. G. Heywood, R. Rannacher and S. Turek, Artificial boundaries and flux and pressure conditions for the incompressible Navier–Stokes equations, Internat. J. Numer. Methods Fluids 22 (1996), no. 5, 325–352. 10.1002/(SICI)1097-0363(19960315)22:5<325::AID-FLD307>3.0.CO;2-YSearch in Google Scholar

[31] J. Hoffman, Efficient computation of mean drag for the subcritical flow past a circular cylinder using general Galerkin G2, Internat. J. Numer. Methods Fluids 59 (2009), no. 11, 1241–1258. 10.1002/fld.1865Search in Google Scholar

[32] J. Hoffman and C. Johnson, Adaptive finite element methods for incompressible fluid flow, Error Estimation and Adaptive Discretization Methods in Computational Fluid Dynamics, Lect. Notes Comput. Sci. Eng. 25, Springer, Berlin (2003), 97–157. 10.1007/978-3-662-05189-4_3Search in Google Scholar

[33] P. Hood and C. Taylor, Navier–Stokes equations using mixed interpolation, Finite Element Methods in Flow Problems, John Wiley, New York (1974), 121–132. Search in Google Scholar

[34] V. John, Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder, Internat. J. Numer. Methods Fluids 44 (2004), no. 7, 777–788. 10.1002/fld.679Search in Google Scholar

[35] V. John, A. Linke, C. Merdon, M. Neilan and L. G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM Rev. 59 (2017), no. 3, 492–544. 10.1137/15M1047696Search in Google Scholar

[36] U. Langer and A. Schafelner, Adaptive space-time finite element methods for non-autonomous parabolic problems with distributional sources, Comput. Methods Appl. Math. 20 (2020), no. 4, 677–693. 10.1515/cmam-2020-0042Search in Google Scholar

[37] U. Langer and A. Schafelner, Adaptive space-time finite element methods for parabolic optimal control problems, J. Numer. Math. 30 (2022), no. 4, 247–266. 10.1515/jnma-2021-0059Search in Google Scholar

[38] U. Langer and O. Steinbach, Space-Time Methods—Applications to Partial Differential Equations, Radon Ser. Comput. Appl. Math. 25, De Gruyter, Berlin, 2019. 10.1515/9783110548488Search in Google Scholar

[39] P. L. Lederer, A. Linke, C. Merdon and J. Schöberl, Divergence-free reconstruction operators for pressure-robust Stokes discretizations with continuous pressure finite elements, SIAM J. Numer. Anal. 55 (2017), no. 3, 1291–1314. 10.1137/16M1089964Search in Google Scholar

[40] A. Linke, On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime, Comput. Methods Appl. Mech. Engrg. 268 (2014), 782–800. 10.1016/j.cma.2013.10.011Search in Google Scholar

[41] A. Linke, G. Matthies and L. Tobiska, Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors, ESAIM Math. Model. Numer. Anal. 50 (2016), no. 1, 289–309. 10.1051/m2an/2015044Search in Google Scholar

[42] D. Meidner and T. Richter, Goal-oriented error estimation for the fractional step theta scheme, Comput. Methods Appl. Math. 14 (2014), no. 2, 203–230. 10.1515/cmam-2014-0002Search in Google Scholar

[43] D. Meidner and T. Richter, A posteriori error estimation for the fractional step theta discretization of the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 288 (2015), 45–59. 10.1016/j.cma.2014.11.031Search in Google Scholar

[44] D. Meidner and B. Vexler, Adaptive space-time finite element methods for parabolic optimization problems, SIAM J. Control Optim. 46 (2007), no. 1, 116–142. 10.1137/060648994Search in Google Scholar

[45] S. Mittal, A. Ratner, D. Hastreiter and T. E. Tezduyar, Space-time finite element computation of incompressible flows with emphasis on flow involving oscillating cylinders, Internat. Video J. Engrg. Res. 1 (1991), 83–86. Search in Google Scholar

[46] J. T. Oden, Adaptive multiscale predictive modelling, Acta Numer. 27 (2018), 353–450. 10.1017/S096249291800003XSearch in Google Scholar

[47] R. Picard and D. McGhee, Partial Differential Equations, De Gruyter Exp. Math. 55, Walter de Gruyter, Berlin, 2011. 10.1515/9783110250275Search in Google Scholar

[48] A. Rademacher, Adaptive finite element methods for nonlinear hyperbolic problems of second order, PhD thesis, Technische Universität Dortmund, 2009. Search in Google Scholar

[49] R. Rannacher, Finite element methods for the incompressible Navier–Stokes equations, Fundamental Directions in Mathematical Fluid Mechanics, Adv. Math. Fluid Mech., Birkhäuser, Basel (2000), 191–293. 10.1007/978-3-0348-8424-2_6Search in Google Scholar

[50] R. Rannacher, Probleme der Kontinuumsmechanik und ihre numerische Behandlung, Heidelberg University, Heidelberg, 2017. Search in Google Scholar

[51] T. Richter and T. Wick, Variational localizations of the dual weighted residual estimator, J. Comput. Appl. Math. 279 (2015), 192–208. 10.1016/j.cam.2014.11.008Search in Google Scholar

[52] A. Schafelner, Space-time finite element methods, PhD thesis, Johannes Kepler University, Linz, 2021. Search in Google Scholar

[53] M. Schäfer, S. Turek, F. Durst, E. Krause and R. Rannacher, Benchmark computations of laminar flow around a cylinder, Flow Simulation with High-Performance Computers II, Vieweg & Teubner, Wiesbaden (1996), 547–566. 10.1007/978-3-322-89849-4_39Search in Google Scholar

[54] M. Schmich, Adaptive finite element methods for computing nonstationary incompressible flows, Doctoral Thesis, Heidelberg University, 2009. Search in Google Scholar

[55] M. Schmich and B. Vexler, Adaptivity with dynamic meshes for space-time finite element discretizations of parabolic equations, SIAM J. Sci. Comput. 30 (2007/08), no. 1, 369–393. 10.1137/060670468Search in Google Scholar

[56] R. Stenberg, Error analysis of some finite element methods for the Stokes problem, Math. Comp. 54 (1990), no. 190, 495–508. 10.1090/S0025-5718-1990-1010601-XSearch in Google Scholar

[57] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, AMS Chelsea, Providence, 2001. 10.1090/chel/343Search in Google Scholar

[58] J. P. Thiele and T. Wick, Variational partition-of-unity localizations of space-time dual weighted residual estimators for parabolic problems, preprint (2022), https://arxiv.org/abs/2207.04764. Search in Google Scholar

[59] S. Turek, Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach, Lect. Notes Comput. Sci. Eng. 6, Springer, Berlin, 1999. 10.1007/978-3-642-58393-3Search in Google Scholar

[60] T. Wick, Solving monolithic fluid-structure interaction problems in arbitrary Lagrangian Eulerian coordinates with the deal.II library, Arch. Numer. Softw. 1 (2013), 10.11588/ans.2013.1.10305. 10.11588/ans.2013.1.10305Search in Google Scholar

Received: 2022-10-06

Revised: 2023-02-20

Accepted: 2023-04-26

Published Online: 2023-05-31

Published in Print: 2024-01-01

Articles in the same Issue

https://doi.org/10.1515/cmam-2022-0200

Keywords for this article

Tensor-Product Space-Time Finite Elements; Dual-Weighted Residuals; Mesh Adaptivity; Dynamically Changing Meshes; Incompressible Navier–Stokes Equations