## 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 [14] as well as the PhD thesis [47] in which similar concepts for full space-time adaptivity for the Navier-Stokes equations were developed. The main differences are two-fold 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 [44]) rather than the filtering approach [17]. On the other hand, another closely related study is [53] in which the current PU-DWR method was developed and analyzed for parabolic problems.

Employing (full) space-time adaptivity using the DWR method dates back to [48], who developed and tested the method for parabolic problems. Shortly afterwards the space-time DWR method was applied to the incompressible Navier-Stokes equations in [47, 14] and for second-order hyperbolic problems (i.e., the elastic wave equation) in [41] and [6]. In this regard, the first space-time finite element formulation (without adaptivity) for incompressible flow was proposed in [52]. Employing space-time (therein called ) for turbulent flow and goal-oriented adaptivity for the efficient computation of the mean drag was done in [28].

Transport problems with coupled flow with an emphasis on adaptivity of the transport part was investigated in [11] and the extension to a multirate framework was suggested in [20]. 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 [37, 38, 23], respectively. Space-time methods with residual-based error control and adaptivity can be found for instance in [33, 45]. Finally, [39] proposed goal-oriented space-time adaptivity with optimal control and [31, 32] 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. These developments include 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 elements with polynomial degree 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 , with with for the velocities and 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 [15], 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 [24, 51, 42, 25, 54, 26]:

###### Formulation 2.1 (incompressible isothermal Navier-Stokes equations).

Find the vector-valued velocity

with and the scalar-valued pressure such that(1a) | |||||

(1b) | |||||

(1c) | |||||

(1d) |

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

(2) |

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

###### Remark 2.2.

By omitting the nonlinear convection term we obtain the Stokes problem.

We employ the following notation for the function spaces. Spatial function spaces are denoted by . Without indices is the analytical function space. is the space of lowest order Taylor-Hood elements (quadratic FE for velocity, linear FE for pressure) and is the space of fourth order FE for velocity and second order FE for pressure. Temporal function spaces are denoted by . Without indices is the analytical function space. is the space of discontinuous finite elements in time of .th order and is the space of discontinuous finite elements in time of .th order. We denote spatio-temporal function spaces by combining a spatial and a temporal function space:

(3) |

will be the space-time Hilbert space on which we define analytical weak formulation. Here denotes the dual space of , which contains bounded linear operators from to . This is a good choice for the analytical function space, since we have a continuous embedding into a function space where the fluid velocity is continuous on [14, 21]. This ensures that the initial condition is well-defined for the velocity.

###### Remark 2.3.

When we have homogeneous Neumann boundary conditions on , the so-called do-nothing (outflow) condition [27], the pressure does not require normalization for its uniqueness and we have the spatio-temporal function space

(4) |

where .

Finally, to simplify notation we introduce

(5) |

In this notation, 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 the weak formulation reads: Find such that

(6) |

where for the Navier-Stokes equations, we have

(7) |

and for the Stokes equations, we have

(8) |

### 2.1 Discretization in time

Let , with be a partitioning of time. Then, we define the semi-discrete space as

(9) |

where the space-time function space (3) has been discretized in time with the discontinuous Galerkin method of order (). is the space of polynomials of order , which map from the time interval into the space . Since functions in can have discontinuities between the time intervals, we define the limits of at time from above and from below for a function

(10) |

and the jump of the function value of at time as

(11) |

The time discretization for the cases and is illustrated in Figure 3.

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

Find such that

(12) | ||||

### 2.2 Discretization in space

Next, we describe the spatial discretization of the formulation for the (Navier-)Stokes equations with continuous finite elements. Let be a shape-regular mesh of the domain . The elements are a non-overlapping covering of , where the elements are quadrilaterals for and hexahedrals for . The discretization parameter is defined element-wise as . After adaptive refinement in space,

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 (

) or face (). Furthermore, on different time intervals we allow for changing spatial meshes, i.e. possibly , 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

(13) |

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., [25, 43]. In the spatial discretization we now replace and with conforming finite element spaces that still fulfill a discrete version of this stability estimate. For this we define the mixed finite element space

(14) |

where is being constructed by mapping tensor-product polynomials of degree from the master element to the element . In our computations, we use the function space , the so-called Taylor-Hood elements [50] with quadratic finite elements for the velocity and linear finite elements for the pressure. Additionally, we use as an enriched space with fourth order finite elements for the velocity and second order finite elements for the pressure. In , we could have also used as an enriched space, since was shown to be an inf-sup stable Taylor-Hood pair [49], but this must not be the case in 3D anymore. Hence, we approximate the enriched velocity by instead of finite elements in space. The fully discretized problem then reads: Find such that

(15) | ||||

### 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 et al. [20][Sec. 4], which uses tensor-product spaces to directly work with the fully discrete weak formulation (15). In this approach the space-time cylinder is divided into space-time slabs(16) |

on which the weak formulation is solved sequentially. The spatial mesh is slabwise constant, i.e. . Hence, we can define the space-time tensor-product finite element basis on 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 a loop over the temporal degrees of freedom in the assembly of the FEM matrices. Furthermore, unlike for the time-stepping formulations of (15) 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.

### 2.4 Linearization of the semilinear form

For the Navier-Stokes equations, we need to include the nonlinear convection term in our variational formulation .
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’s method.
The resulting nonlinear iteration step reads as:

Solve the problem:

Find such that

(17) |

Obtain the new iterate by:

(18) |

###### Remark 2.4.

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

Hence, we now need to compute the Fréchet derivative of the semi-linear form with respect to in the direction of . For a linear functional

(19) |

holds. Hence, it remains to compute the Fréchet derivative of the convection term

(20) |

which directly follows from the linearity of differentiation as

(21) |

The derivative of the complete semi-linear form is thus given by

(22) | ||||

Finally, we can formulate the full algorithm of Newton’s method, which we use for the numerical experiments.

Parameters:
, , ,

Input: Initial guess on slab

Output: FEM solution on slab .

(23) |

(24) |

(25) |

(26) |

For the numerical tests, we used as the maximum number of line search steps and as the line search damping parameter. We set the maximum number of Newton steps to and the tolerance to .

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 and , we define

(27) |

###### Remark 2.5.

Alternatively one could create the initial guess as the temporal extrapolant of the entire solution from the previous slab, e.g. 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 become homogeneous, i.e. we enforce on the Dirichlet boundary.

We observe that in Algorithm 1 we need to compute the directional derivative in each step of the Newton method. This is computationally expensive, since we need to assemble the Jacobian matrix in every step. Instead, in our program we use the approximation

(28) |

when the reduction rate

(29) |

is sufficiently small. We use these so called simplified-Newton steps, when .

### 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 [15] 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 [18, 8, 2, 34]. To avoid the resulting numerical artifacts, Besier and Wollner [15] 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 and a divergence free projection [15].

###### Remark 2.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 [35] 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 -conforming, i.e. discretely divergence free, finite elements inside the weak formulation of the (Navier-)Stokes equations [35],[36],[30],[34]. 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 [9, 10].

This problem can be formalized as follows:

Let and be two space-time slabs with different spatial triangulations . Then, we need to find a projection of the initial condition
into the function space of the next slab
such that the projected initial condition is divergence free, i.e.

(30) |

To achieve this goal Besier and Wollner [15] presented and analyzed a divergence free projection and a divergence free projection.

Divergence free projection

Find such that

(31) | ||||

Divergence free projection

Find such that

(32) | ||||

## 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 [13, 12] (see also other references [7, 48, 14, 1, 40]) 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) [44, 53] 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 of the form

(33) |

which represents our physical quantity of interest. We now want to minimize the difference between the quantity of interest of the analytical solution and its FEM approximation , i.e.

(34) |

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

(35) | ||||

(36) | ||||

(37) |

The stationary points , and of the Lagrange functionals , and need to satisfy the Karush-Kuhn-Tucker first order optimality conditions. Firstly, the stationary points are solutions to the equations

(38) | ||||

(39) | ||||

(40) |

We call these equations the primal problems and their solutions , and the primal solutions. Secondly, the stationary points must also satisfy the equations

(41) | ||||

(42) | ||||

(43) |

These equations are being called the adjoint or dual problems and their solutions , and 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

(44) |

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

(45) | ||||

where we denote the primal and adjoint residuals as

(46) |

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 , and with respect to the adjoint solution . Using the definition of the continuous (), the time-discrete () and the fully discrete () variational formulations, we observe that they are linear in the test functions and hence we have

(47) | ||||

(48) | ||||

(49) |

Therefore, 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 , and with respect to the primal solution . We then get