# An energy-efficient GMRES–Multigrid solver for space-time finite element
computation of dynamic poro- and thermoelasticity

Mathias Anselmann∗<EMAIL_ADDRESS>(corresponding author) Markus Bause∗
Nils Margenberg∗ Pavel Shamko∗
∗ Helmut Schmidt University, Faculty of Mechanical and Civil Engineering,
Holstenhofweg 85,
22043 Hamburg, Germany

###### Abstract

Abstract

We present families of space-time finite element methods (STFEMs) for a
coupled hyperbolic-parabolic system of poro- or thermoelasticity. Well-
posedness of the discrete problems is proved. Higher order approximations
inheriting most of the rich structure of solutions to the continuous problem
on computationally feasible grids are naturally embedded. However, the block
structure and solution of the algebraic systems become increasingly complex
for these members of the families. We present and analyze a robust geometric
multigrid (GMG) preconditioner for GMRES iterations. The GMG method uses a
local Vanka-type smoother. Its action is defined in an exact mathematical way.
Due to nonlocal coupling mechanisms of unknowns, the smoother is applied on
patches of elements. This ensures the damping of error frequencies. In a
sequence of numerical experiments, including a challenging three-dimensional
benchmark of practical interest, the efficiency of the solver for STFEMs is
illustrated and confirmed. Its parallel scalability is analyzed. Beyond this
study of classical performance engineering, the solver’s energy efficiency is
investigated as an additional and emerging dimension in the design and tuning
of algorithms and their implementation on the hardware.

Keywords. Poro-/thermoelasticity, dynamic Biot model, space-time finite
element approximation, geometric multigrid method, GMRES iterations,
computational efficiency, energy consumption.

##  1 Introduction

###  1.1 Mathematical model

We present and investigate numerically a geometric multigrid preconditioning
technique, based on a local Vanka-type smoother, for solving by GMRES
iterations the linear systems that arise from space-time finite element
discretizations of the coupled hyperbolic-parabolic system of dynamic
poroelasticity

$\displaystyle\rho\partial_{t}^{2}\boldsymbol{u}-\nabla\cdot(\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{u}))+\alpha\boldsymbol{\nabla}p$
$\displaystyle=\rho\boldsymbol{f}\,,$ $\displaystyle\quad\text{in
}\;\Omega\times(0,T]\,,$ (1.1a) $\displaystyle
c_{0}\partial_{t}p+\alpha\nabla\cdot\partial_{t}\boldsymbol{u}-\nabla\cdot(\boldsymbol{K}\nabla
p)$ $\displaystyle=g\,,$ $\displaystyle\quad\text{in }\;\Omega\times(0,T]\,,$
(1.1b)
$\displaystyle\boldsymbol{u}(0)=\boldsymbol{u}_{0}\,,\quad\partial_{t}\boldsymbol{u}(0)=\boldsymbol{u}_{1}\,,\quad
p(0)$ $\displaystyle=p_{0}\,,$ $\displaystyle\quad\text{in
}\;\Omega\times\\{0\\}\,,$ (1.1c) $\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}_{D}\,,$ $\displaystyle\quad\text{on
}\;\Gamma_{\boldsymbol{u}}^{D}\times(0,T]\,,$ (1.1d)
$\displaystyle-(\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{u})-\alpha
p\boldsymbol{E})\boldsymbol{n}$ $\displaystyle=\boldsymbol{t}_{N}\,,$
$\displaystyle\quad\text{on }\;\Gamma_{\boldsymbol{u}}^{{N}}\times(0,T]\,,$
(1.1e) $\displaystyle p$ $\displaystyle=p_{D}\,,$ $\displaystyle\quad\text{on
}\;\Gamma_{p}^{{D}}\times(0,T]\,,$ (1.1f) $\displaystyle-\boldsymbol{K}\nabla
p\cdot\boldsymbol{n}$ $\displaystyle=p_{N}\,,$ $\displaystyle\quad\text{on
}\;\Gamma_{p}^{{N}}\times(0,T]\,.$ (1.1g)

In (1.1), $\Omega\subset\mathbb{R}^{d}$, with $d\in\\{2,3\\}$, is an open
bounded Lipschitz domain with outer unit normal vector $\boldsymbol{n}$ to the
boundary $\partial\Omega$ and $T>0$ is the final time point. We let
$\partial\Omega={\Gamma_{\boldsymbol{u}}^{D}}\cup{\Gamma_{\boldsymbol{u}}^{N}}$
and $\partial\Omega={\Gamma_{p}^{D}}\cup{\Gamma_{p}^{N}}$ with closed portions
$\Gamma_{\boldsymbol{u}}^{D}$ and $\Gamma_{p}^{D}$ of non-zero measure.
Important applications of the model (1.1), that is studied as a prototype
system, arise in poroelasticity; cf. [57] and [13, 14, 15]. In poroelasticity,
Eqs. (1.1) are referred to as the dynamic Biot model. The system (1.1) is used
to describe flow of a slightly compressible viscous fluid through a deformable
porous matrix. The small deformations of the matrix are described by the
Navier equations of linear elasticity, and the diffusive fluid flow is
described by Duhamel’s equation. The unknowns are the effective solid phase
displacement $\boldsymbol{u}$ and the effective fluid pressure $p$. The
quantity
$\boldsymbol{\varepsilon}(\boldsymbol{u}):=(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{\top})/2$
denotes the symmetrized gradient or strain tensor and
$\boldsymbol{E}\in\mathbb{R}^{d,d}$ is the identity matrix. Further, $\rho$ is
the effective mass density, $\boldsymbol{C}$ is Gassmann’s fourth order
effective elasticity tensor, $\alpha$ is Biot’s pressure-storage coupling
tensor, $c_{0}$ is the specific storage coefficient and $\boldsymbol{K}$ is
the permeability field. For brevity, the positive quantities $\rho>0$,
$\alpha>0$ and $c_{0}>0$ as well as the tensors $\boldsymbol{C}$ and
$\boldsymbol{K}$ are assumed to be constant in space and time. The tensors
$\boldsymbol{C}$ and $\boldsymbol{K}$ are assumed to be symmetric and positive
definite,

$\displaystyle\exists
k_{0}>0\;\forall\boldsymbol{\xi}=\boldsymbol{\xi}^{\top}\in\mathbb{R}^{d,d}:$
$\displaystyle\quad\sum_{i,j,k,l=1}^{d}\xi_{ij}C_{ijkl}\xi_{kl}\geq
k_{0}\sum_{j,k=1}^{d}|\xi_{jk}|^{2}\,,$ (1.2a) $\displaystyle\exists
k_{1}>0\;\forall\boldsymbol{\xi}\in\mathbb{R}^{d}:$
$\displaystyle\quad\sum_{i,j,=1}^{d}\xi_{i}K_{ij}\xi_{j}\geq
k_{1}\sum_{i=1}^{d}|\xi_{i}|^{2}\,.$ (1.2b)

In order to enhance physical realism in poroelasticity, generalizations of the
model (1.1) have been developed and investigated in, e.g., [16, 49]. We note
that the system (1.1) is also formally equivalent to the classical coupled
thermoelasticity system which describes the flow of heat through an elastic
structure; cf. [35] and [20, 43]). In that context, $p$ denotes the
temperature, $c_{0}$ is the specific heat of the medium, and $\boldsymbol{K}$
is the conductivity. Then, the quantity $\alpha\nabla p$ arises from the
thermal stress in the structure, and the term
$\alpha\nabla\cdot\partial_{t}\boldsymbol{u}$ models the internal heating due
to the dilation rate. Well-posedness of (1.1) is ensured; cf., e.g., [35, 58,
55]. This can be shown by different mathematical techniques, by semigroup
methods [35, Thm. 2.2], Rothe’s method [58, Thm. 4.18 and Cor. 4.33] and
Picard’s theorem [55, Thm. 6.2.1]. In these works, boundary conditions
different to the ones in (1.1d) to (1.1d) are partly used.

###  1.2 Space-time finite element and multigrid techniques

The coupled hyperbolic-parabolic structure of the system (1.1) of partial
differential equations adds an additional facet of complexity onto its
numerical simulation. A natural and promising approach for the numerical
approximation of coupled systems is given by the application of space-time
finite element methods that are based on an uniform treatment of the space and
time variables by variational formulations. Therein, the discrete unknown
functions are defined on the entire space-time domain $\Omega\times[0,T]$ and
are represented in terms of basis functions. This facilitates the
discretization of even complex coupling terms that may involve combinations of
temporal and spatial derivatives as in (1.1b). Moreover, space-time finite
element methods enable the natural construction of higher order members in the
families of schemes. Higher order discretizations offer the potential to
achieve accurate results on computationally feasible grids with a minimum of
numerical costs. In particular for second-order hyperbolic problems, higher
order time discretization schemes have shown their superiority, measured in
terms of accuracy over temporal mesh resolution and compute time; cf., e.g.,
[5]. Finally, space-time adaptivity based on a-posteriori error control by
duality concepts and multi-rate in time approaches become feasible. For this
we refer to [9] for the presentation of the abstract framework of duality
based error control and, e.g., to [8, 19] for its application to physical
problems. For the variational time discretization of second-order hyperbolic
and parabolic equations, various families have been proposed and investigated
recently; cf., e.g., [5, 10, 41, 25, 60]. They differ by the approximation
order and regularity in time of the discrete solution as well as stability and
energy conservation properties. In this work, we use a discontinuous Galerkin
time discretization of arbitrary polynomial order and focus on the efficient
solution of the resulting complex algebraic block systems. Higher order
(variational) time discretizations for three-dimensional applications put a
facet of complexity on the algebraic solver and have rarely been studied yet
in the literature. For the discretization of the space variables, numerous
finite element approaches have been proposed and studied, in particular, for
the quasi-static Biot system of flow in deformable porous media. The latter
model results from (1.1) by neglecting the acceleration term
$\rho\partial_{t}^{2}\boldsymbol{u}$ in (1.1a). For the quasi-static Biot
model, three-field formulations with the fluid flux
$\boldsymbol{q}=-\boldsymbol{K}\nabla p$ as an additional explicit variable
and the application of mixed finite element methods for the spatial
discretization of $p$ and $\boldsymbol{q}$ attracted strong interest of
researchers; cf., e.g., [11, 30, 40, 51, 52] and the references therein. In
this work, we use an approximation of the variables $\boldsymbol{u}$ and $p$
by higher order inf-sup stable pairs of finite element spaces with continuous
or discontinuous discrete variable $p$. Precisely, either the well-known
Taylor–Hood family of finite element pairs for $\boldsymbol{u}$ and $p$ or a
continuous approximation of $\boldsymbol{u}$ combined with a discontinuous
approximation of $p$ is used; cf. (2.6). In the literature, the latter option
has shown its superiority for the approximation of incompressible viscous flow
problems [38]. Structure-preserving approximations with the fluid flux and
stress tensor as additional variables, that are based on a reformulation of
the coupled problem (1.1) as a first-order in space and time system (cf. [25,
Subsec. 7.1]), remain as a work for the future.

Figure 1.1: Space-time mesh for a piecewise linear ($k=1$) discontinuous
Galerkin time discretization and a Lagrange basis w.r.t. the $k+1$ Gauss-Radau
quadrature points of $I_{n}$.

In our space-time finite element approach, we choose a discontinuous temporal
test basis such a time-marching process is obtained. On each subinterval
$I_{n}=(t_{n-1},t_{n}]$ of the time mesh
$\mathcal{M}_{\tau}:=\\{I_{1},\ldots,I_{N}\\}$, we use a Lagrangian basis with
respect to the $(k+1)$ Gauss–Radau quadrature points of $I_{n}$ for the time
interpolation. Two discrete systems differing by the treatment of the term
$\nabla\cdot\partial_{t}\boldsymbol{u}$ in (1.1b) are then proposed. Dirichlet
boundary conditions are implemented in a weak form by means of Nitsche’s
method and/or a discontinuous (spatial) approximation of the scalar variable
$p$. Well-posedness of the discrete problems is carefully proved for arbitrary
polynomial degrees of the discretization. The degrees of freedom on $I_{n}$
are associated with the $(k+1)$ Gauss–Radau quadrature nodes of this
subinterval; cf. Fig. 1.1. On each subinterval $I_{n}$, the resulting
algebraic problem is a linear system of equations with a block matrix of
$(k+1)\times(k+1)$ subsystems where each of the subsystems exhibits a complex
structure itself. Here, $k$ denotes the piecewise polynomial order of the time
discretization; cf. [3, 2, 32]. The solution of these algebraic problems with
complex block structure demands for efficient and robust iteration schemes.
Geometric multigrid (GMG) methods are known as the most efficient iterative
techniques for the solution of large linear systems arising from the
discretization of partial differential equations, particularly of elliptic
type; cf. [17, 29, 62] for a general overview. In contrast to other iterative
solvers (like for instance BiCGStab or GMRES), multigrid solvers converge with
a rate which is independent of the (spatial) mesh size and require
computational costs that are only linearly dependent on the number of
unknowns. Massively parallel implementations of GMG methods on modern
architectures show excellent scalability properties and their high efficiency
has been recognized, e.g., in [45, 26, 26]. Today they are widely used and are
employed in many variants; cf. [23, 39, 63] for flow and saddle point
problems. Nowadays, GMG methods are employed as preconditioner in Krylov
subspace iterations, for instance GMRES iterations, to enhance their
robustness, which is also done here. Analyses of GMG methods (cf., e.g., [23,
31, 46]) have been done for linear systems in saddle point form, with matrix
$\boldsymbol{\mathcal{A}}=\begin{pmatrix}\boldsymbol{A}&\boldsymbol{B}^{\top}\\\
\boldsymbol{B}&-\boldsymbol{C}\end{pmatrix}$ and symmetric and positive
definite submatrices $\boldsymbol{A}$ and $\boldsymbol{C}$, arising for
instance from mixed discretizations of the Stokes problem. If higher order
discontinuous Galerkin time stepping is used, the resulting linear system
matrix is a $(k+1)\times(k+1)$ block matrix with each of its blocks being of
the form $\boldsymbol{\mathcal{A}}$. This imposes an additional facet of
complexity on the geometric multigrid preconditioner. To the best of our
knowledge, analyses of GMG methods arising form higher order (variational)
time discretizations of the coupled system (1.1) or related ones are still
missing in the literature.

GMG methods exploit different mesh levels of the underlying problem in order
to reduce different frequencies of the error by employing a relatively cheap
smoother on each grid level. Different iterative methods have been proposed in
the literature as smoothing procedure; cf. [23] and the references therein.
They range from low-cost methods like Richardson, Jacobi, and SOR applied to
the normal equation of the linear system to collective smoothers, that are
based on the solution of small local problems. Here, we use a Vanka-type
smoother [66, 46, 39] of the family of collective methods. Numerical
computations show that an elementwise application of the Vanka smoother fails
to reduce the high frequencies of the error on the multigrid levels. The
reason for this comes through spatial interelement couplings of degrees of
freedom of the scalar variable $p$ in (1.1). This holds for the discontinuous
and continuous spatial approximation of $p$. For the spatially continuous
approximation of $p$, the coupling of degrees of freedom is obvious. For the
spatially discontinuous approximation of $p$, the interelement coupling of
degrees of freedom of $p$ is due to the presence of the diffusion term in
(1.1b) and the surface integrals involving averages and jumps of the discrete
solution over interelement faces resulting from its discretization. This is in
contrast to the Navier–Stokes equations for which elementwise Vanka-type
smoother are known to perform excellently if discontinuous pressure
approximations are applied; cf. [3, 2, 39]. The difference comes through the
fact that in case of the Navier–Stokes system no interelement coupling of
pressure degrees of freedom arises. As a remedy, we propose the application of
the Vanka-type smoother on cell patches that are linked to the grid nodes and
built from four neighbored cells in two dimensions and eight neighbored
elements in three dimensions, with appropriate adaptations for grid nodes
close to or on the domain’s boundary. Further we employ some relaxation
strategy in the smoothing steps. Then an efficient damping of frequencies in
the error on the multigrid hierarchy is obtained. Here, this Vanka-type
smoother is presented in an mathematically exact way and its performance
properties are investigated by numerical experiments. Our numerical
experiments confirm that GMRES iterations that are preconditioned by the
proposed multigrid method converge at a rate which is (nearly) independent of
the mesh sizes in space and the time step size; cf. also [3, 2]. In the
smoothing steps, the application of the Vanka-type operator requires the
solution of patchwise linear systems of equations. For (1.1) and three space
dimensions the size of the local systems no longer continues to remain
negligible, in particular if a higher polynomial order in space and/or time is
applied. Their solution deserves attention. In our approach we compute of
inverse of the patchwise system matrix by LAPACK routines to solve to local
system of equations. Details of an data-structure efficient MPI-based parallel
implementation of the GMG algorithm in our in-house high-performance frontend
solver for the deal.II finite element library [6] can be found in [2] where
this is presented for the Navier–Stokes system. Multi-field formulations of
the system (1.1), that are not studied in this work, further increase the
complexity of the arising linear systems; cf. [25, Subsec. 7.1]. However, the
methodology of the linear solver, that is developed here, offers the potential
of its generalization to more advanced problem formulations and their
discretizations.

###  1.3 Energy efficiency

In the past, performance engineering and hardware engineering for large scale
simulations of physical phenomena have been eclipsed by the longing for ever
more performance where faster seemed to be the only paradigm. ”Classical”
performance engineering has been applied to enhance, firstly, the efficiency
of the current method on the target hardware or to find numerical alternatives
that might better fit to the hardware in use and/or, secondly, to develop
other numerical methods can be found to improve the numerical efficiency.
Tuning both simultaneously is called hardware-oriented numerics in the
literature; cf. [64, 65]. Since recently, a growing awareness of energy
consumption in computational science, particularly, in extreme scale computing
with a view to exascale computing has raised; cf., e.g., [54]. It has been
observed that as a consequence of decades of performance-centric hardware
development there is a huge gap between pure performance and energy
efficiency. An analysis of our algorithm’s parallel scaling and energy
consumption properties by performance models exceeds the scope of this work
and would overburden it. However, since energy consumption of application
codes on the available hardware is of growing awareness and a key for future
improvements, we study the energy consumption and parallel scaling properties
of our algorithm and its implementation by three-dimensional numerical
experiments. The development of a proper model that quantifies performance and
energy efficiency in some appropriate metric and can be used for a code
optimization still deserves research and is left as a work for the future.

###  1.4 Outline of the work

This work is organized as follows. In Sec. 2 we introduce our notation. In
Sec. 3 the space-time finite element approximation of arbitrary order of (1.1)
is derived and well-posedness of the fully discrete problem is proved. Our
GMRES–GMG solver is introduced in Sec. 4. In Sec. 5 our performed numerical
computations for analyzing the performance properties of the overall approach
are presented. In Sec. 6 we end with a summary and conclusions. In the
appendix, supplementary results are summarized.

##  2 Basic notation

In this work, standard notation is used. We denote by $H^{1}(\Omega)$ the
Sobolev space of $L^{2}(\Omega)$ functions with first-order weak derivatives
in $L^{2}(\Omega)$. Further, $H^{-1}(\Omega)$ is the dual space of
$H^{1}_{0}(\Omega)$, with the standard modification if the Dirichlet condition
is prescribed on a part $\Gamma^{D}\subset\Gamma$ of the boundary
$\partial\Omega$ only; cf. (1.1). The latter is not explicitly borne out by
the notation of $H^{-1}(\Omega)$. It is always clear from the context. For
vector-valued functions we write those spaces bold. By
$\langle\cdot,\cdot\rangle_{S}$ we denote the $L^{2}(S)$ inner product for a
domain $S$. For $S=\Omega$, we simply write $\langle\cdot,\cdot\rangle$. For
the norms of the Sobolev spaces the notation is

$\displaystyle\|\cdot\|:=\|\cdot\|_{L^{2}}\,,\qquad\|\cdot\|_{1}:=\|\cdot\|_{H^{1}}\,.$

For short, we put

$Q:=L^{2}(\Omega)\qquad\text{and}\qquad\boldsymbol{V}:=\left(H^{1}(\Omega)\right)^{d}\,.$

For a Banach space $B$, we let $L^{2}(0,T;B)$ be the Bochner space of
$B$-valued functions, equipped with its natural norm. For a subinterval
$J\subseteq[0,T]$, we will use the notation $L^{2}(J;B)$ for the corresponding
Bochner space. In what follows, the constant $c$ is generic and indepedent of
the size of the space and time meshes.

For the time discretization, we decompose the time interval $I:=(0,T]$ into
$N$ subintervals $I_{n}=(t_{n-1},t_{n}]$, $n=1,\ldots,N$, where
$0=t_{0}<t_{1}<\cdots<t_{N-1}<t_{N}=T$ such that $I=\bigcup_{n=1}^{N}I_{n}$.
We put $\tau:=\max_{n=1,\ldots,N}\tau_{n}$ with $\tau_{n}=t_{n}-t_{n-1}$.
Further, the set $\mathcal{M}_{\tau}:=\\{I_{1},\ldots,I_{N}\\}$ of time
intervals is called the time mesh. For a Banach space $B$ and any
$k\in\mathbb{N}_{0}$, we let

$\mathbb{P}_{k}(I_{n};B):=\bigg{\\{}w_{\tau}\,:\,I_{n}\to
B\,,\;w_{\tau}(t)=\sum_{j=0}^{k}W^{j}t^{j}\;\forall t\in I_{n}\,,\;W^{j}\in
B\;\forall j\bigg{\\}}\,.$ (2.1)

For $k\in\mathbb{N}_{0}$ we define the space of piecewise polynomial functions
in time with values in $B$ by

$Y_{\tau}^{k}(B):=\left\\{w_{\tau}:\overline{I}\rightarrow B\mid
w_{\tau}{}_{|I_{n}}\in\mathbb{P}_{k}(I_{n};B)\;\forall
I_{n}\in\mathcal{M}_{\tau},\,w_{\tau}(0)\in B\right\\}\subset L^{2}(I;B)\,.$
(2.2)

For any function $w:\overline{I}\to B$ that is piecewise sufficiently smooth
with respect to the time mesh $\mathcal{M}_{\tau}$, for instance for $w\in
Y^{k}_{\tau}(B)$, we define the right-hand sided and left-hand sided limit at
a mesh point $t_{n}$ by

$w^{+}(t_{n}):=\lim_{t\to
t_{n}+0}w(t),\quad\text{for}\;n<N\,,\qquad\text{and}\qquad
w^{-}(t_{n}):=\lim_{t\to t_{n}-0}w(t),\quad\text{for}\;n>0\,.$ (2.3)

For the integration in time of a discontinuous Galerkin approach it is natural
to use the right-sided $(k+1)$-point Gauss–Radau quadrature formula. On the
subinterval $I_{n}$, it reads as

$Q_{n}(w):=\frac{\tau_{n}}{2}\sum_{\mu=1}^{k+1}\hat{\omega}_{\mu}^{\operatorname{GR}}w(t_{n,\mu}^{\operatorname{GR}})\approx\int_{I_{n}}w(t)\,\mathrm{d}t\,,$
(2.4)

where
$t_{n,\mu}^{\operatorname{GR}}=T_{n}(\hat{t}_{\mu}^{\operatorname{GR}})$, for
$\mu=1,\ldots,k+1$, are the Gauss–Radau quadrature points on $I_{n}$ and
$\hat{\omega}_{\mu}^{\operatorname{GR}}$ the corresponding weights. Here,
$T_{n}(\hat{t}):=(t_{n-1}+t_{n})/2+(\tau_{n}/2)\hat{t}$ is the affine
transformation from $\hat{I}=[-1,1]$ to $I_{n}$ and
$\hat{t}_{\mu}^{\operatorname{GR}}$ are the Gauss–Radau quadrature points on
$\hat{I}$. Formula (2.4) is exact for all polynomials
$w\in\mathbb{P}_{2k}(I_{n};\mathbb{R})$. In particular, there holds that
$t_{n,k+1}^{\operatorname{GR}}=t_{n}$.

For the space discretization, let $\\{\mathcal{T}_{l}\\}_{l=0}^{L}$ be the
decomposition on every multigrid level of $\Omega$ into (open) quadrilaterals
or hexahedrals, with $\mathcal{T}_{l}=\\{K_{i}\mid
i=1,\ldots,N^{\text{el}}_{l}\\}$, for $l=0,\ldots,L$. These element types are
chosen for our implementation (cf. Sec. 5) that is based on the deal.II
library [6]. The finest partition is $\mathcal{T}_{h}=\mathcal{T}_{L}$. We
assume that all the partitions $\\{\mathcal{T}_{l}\\}_{l=0}^{L}$ are quasi-
uniform with characteristic mesh size $h_{l}$ and $h_{l}=\gamma h_{l-1}$,
$\gamma\in(0,1)$ and $h_{0}=\mathcal{O}(1)$. On the actual mesh level, the
finite element spaces used for approximating the unknowns $\boldsymbol{u}$ and
$p$ of (1.1) are of the form ($l\in\\{0,\ldots,L\\}$)

$\displaystyle\boldsymbol{V}_{h_{l}}^{l}$
$\displaystyle:=\\{\boldsymbol{v}_{h}\in\boldsymbol{V}\cap
C(\overline{\Omega})^{d}\;:\;\boldsymbol{v}_{h}{}_{|K}\in{\boldsymbol{V}(K)}\;\;\text{for
all}\;K\in\mathcal{T}_{l}\\}\,,$ (2.5a) $\displaystyle
Q_{h_{l}}^{l,\text{cont}}$ $\displaystyle:=\\{\boldsymbol{q}_{h}\in Q\cap
C(\overline{\Omega})\;:\;\boldsymbol{q}_{h}{}_{|K}\in{Q(K)}\;\;\text{for
all}\;K\in\mathcal{T}_{l}\\}\,,$ (2.5b) $\displaystyle
Q_{h_{l}}^{l,\text{disc}}$ $\displaystyle:=\\{\boldsymbol{q}_{h}\in
Q\;:\;\boldsymbol{q}_{h}{}_{|K}\in{Q(K)}\;\;\text{for
all}\;K\in\mathcal{T}_{l}\\}\,.$ (2.5c)

By an abuse of notation, we skip the index $l$ of the mesh level when it is
clear from the context and put

$\boldsymbol{V}_{h}:=\boldsymbol{V}_{h}^{l}\quad\text{and}\quad
Q_{h}:=Q_{h}^{l}\;\;\text{with}\;\;Q_{h}^{l}\in\\{Q_{h}^{l,\text{cont}},Q_{h}^{l,\text{disc}}\\}\,.$
(2.6)

For the local spaces $\boldsymbol{V}(K)$ and $Q(K)$ we employ mapped versions
of the pairs $\mathbb{Q}_{r}^{d}/\mathbb{Q}_{r-1}$ and
$\mathbb{Q}_{r}^{d}/\mathbb{P}_{r-1}^{\operatorname{disc}}$, for $r\geq 2$.
The pair $\mathbb{Q}_{r}^{d}/\mathbb{Q}_{r-1}$ with a (globally) continuous
approximation of the scalar variable $p$ in $Q_{h}^{l,\text{cont}}$ is the
well-known Taylor–Hood family of finite element spaces. The pair
$\mathbb{Q}_{r}^{d}/\mathbb{P}_{r-1}^{\operatorname{disc}}$ comprises a
discontinuous approximation of $p$ in the broken polynomial space
$Q_{h}^{l,\text{disc}}$. For the Navier–Stokes equations, the multigrid method
has shown to work best for higher-order finite element spaces with
discontinuous discrete pressure; cf. [38] and [2]. For a further discussion of
mapped and unmapped versions of the pair
$\mathbb{Q}_{r}^{d}/\mathbb{P}_{r-1}^{\operatorname{disc}}$ we refer to [36,
Subsec. 3.6.4]. For an analysis of stability properties of (spatial)
discretizations for the quasi-static Biot system we refer to, e.g., [53]. In
particular, both choices of the local finite element spaces,
$\mathbb{Q}_{r}^{d}/\mathbb{Q}_{r-1}$ and
$\mathbb{Q}_{r}^{d}/\mathbb{P}_{r-1}^{\operatorname{disc}}$, satisfy the
inf–sup stability condition,

$\inf_{q_{h}\in
Q_{h}\backslash\\{0\\}}\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}\backslash\\{\boldsymbol{0}\\}}\dfrac{b(\boldsymbol{v}_{h},q_{h})}{\|\boldsymbol{v}_{h}\|_{1}\,\|q_{h}\|}\geq\beta>0\,,$
(2.7)

for some constant $\beta$ independent of $h$; cf. [36, 48]. In [7, 47] optimal
interpolation error estimates for mapped finite elements on quadrilaterals and
hexahedra are studied. It turned out that the optimality is given for special
families of triangulations. In two and three dimensions, families of meshes,
which are obtained by a regular uniform refinement of an initial coarse grid,
are among these special families. Such a regular refinement that is natural
for the construction of the multigrid hierarchy is used in our computations.
Thus, for $\boldsymbol{v}\in\boldsymbol{H}^{r+1}(\Omega)$ and $q\in
H^{r}(\Omega)$ there exist approximations
$i_{h}\boldsymbol{v}\in\boldsymbol{V}_{h}$ and $j_{h}q\in Q_{h}$ such that,
with some generic constant $c>0$ independent of $h$,

$\|\boldsymbol{v}-i_{h}\boldsymbol{v}\|+h\|\nabla(\boldsymbol{v}-i_{h}\boldsymbol{v})\|\leq
ch^{r+1}\qquad\text{and}\qquad\|q-j_{h}q\|\leq ch^{r}\,.$ (2.8)

##  3 Space-time finite element approximation

For the discretization we rewrite (1.1) as a first-order in time system by
introducing the new variable $\boldsymbol{v}:=\partial_{t}\boldsymbol{u}$.
Then, we recover (1.1a) and (1.1b) as

$\displaystyle\partial_{t}\boldsymbol{u}-\boldsymbol{v}$
$\displaystyle=\boldsymbol{0}$ (3.1a)
$\displaystyle\rho\partial_{t}\boldsymbol{v}-\nabla\cdot(\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{u}))+\alpha\boldsymbol{\nabla}p$
$\displaystyle=\rho\boldsymbol{f}\,,$ (3.1b) $\displaystyle
c_{0}\partial_{t}p+\alpha\nabla\cdot\boldsymbol{v}-\nabla\cdot(\boldsymbol{K}\boldsymbol{\nabla}p)$
$\displaystyle=g$ (3.1c)

along with the initial and boundary conditions 1.1c to 1.1g. We employ
discontinuous Galerkin methods (cf. [61]) for the discretization of the time
variable and inf-sup stable pairs of finite elements (cf. Sec. 2) for the
approximation of the space variables in 3.1. The derivation of the discrete
scheme, presented below in Problem B.1, is standard and not explicitly
presented here. It follows the lines of [10], where continuous in time
Galerkin methods are applied to (1.1), and [3, 2, 34, 33], where discontinuous
in time Galerkin methods are used to discretize the Navier–Stokes system. In
contrast to [10], Dirichlet boundary conditions are implemented here by
Nitsche’s method [12, 24, 50]. This yields a strong link between two different
families of inf-stable finite element pairs for the space discretization. The
main reason for using Nitsche’s method here comes through our more general
software framework. It captures problems on evolving domains solved on fixed
computational background grids (cf. [4]), where Nitsche’s method is applied.
We note that Nitsche’s method does not perturb the convergence behavior of the
space-time discretization; cf. Sec. 5.

For the discrete scheme we need further notation. On the multigrid level $l$
with decomposition $\mathcal{T}_{l}$, for
$\boldsymbol{w}_{h},\boldsymbol{\chi}_{h}\in\boldsymbol{V}_{h}$ and
$q_{h},\psi_{h}\in Q_{h}$ we define

$\displaystyle A_{\gamma}(\boldsymbol{w}_{h},\boldsymbol{\chi}_{h})$
$\displaystyle:=\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{w}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{\chi}_{h})\rangle-\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{w}_{h})\boldsymbol{n},\boldsymbol{\chi}_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}+a_{\gamma}(\boldsymbol{w}_{h},\boldsymbol{\chi}_{h})\,,$
(3.2a) $\displaystyle C(\boldsymbol{\chi}_{h},q_{h})$
$\displaystyle:=-\alpha\langle\nabla\cdot\boldsymbol{\chi}_{h},q_{h}\rangle+\alpha\langle\boldsymbol{\chi}_{h}\cdot\boldsymbol{n},q_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\,,$
(3.2b) $\displaystyle B_{\gamma}(q_{h},\psi_{h})$
$\displaystyle:=\left\\{\begin{array}[]{@{}ll}\langle\boldsymbol{K}\nabla
q_{h},\nabla\psi_{h}\rangle-\langle\boldsymbol{K}\nabla
q_{h}\cdot\boldsymbol{n},\psi_{h}\rangle_{\Gamma_{p}^{D}}+b_{\gamma}(q_{h},\psi_{h})\,,&\text{for}\;\;Q_{h}=Q_{h}^{l,\text{cont}}\,,\\\\[12.91663pt]
\begin{array}[]{@{}l}\displaystyle\sum_{K\in\mathcal{T}_{l}}\langle\boldsymbol{K}\nabla
q_{h},\nabla\psi_{h}\rangle_{K}-\sum_{F\in\mathcal{F}_{h}}\big{(}\langle\\{\mskip-5.0mu\\{\boldsymbol{K}\nabla
q_{h}\\}\mskip-5.0mu\\}\cdot\boldsymbol{n},\lsem\psi_{h}\rsem\rangle_{F}\\\\[12.91663pt]
\displaystyle\quad+\langle\lsem
q_{h}\rsem,\\{\mskip-5.0mu\\{\boldsymbol{K}\nabla\psi_{h}\\}\mskip-5.0mu\\}\cdot\boldsymbol{n}\rangle_{F}\big{)}+\sum_{F\in\mathcal{F}_{h}}\frac{\gamma}{h_{F}}\langle\lsem
q_{h}\rsem,\lsem\psi_{h}\rsem\rangle_{F}\,,\end{array}&\text{for}\;\;Q_{h}=Q_{h}^{l,\text{disc}}\,,\\\\[4.30554pt]
\end{array}\right.$ (3.2g)

where, for
$\boldsymbol{w}\in\boldsymbol{H}^{1/2}(\Gamma^{D}_{\boldsymbol{u}})$ and $q\in
H^{1/2}(\Gamma^{D}_{p})$,

$\displaystyle a_{\gamma}(\boldsymbol{w},\boldsymbol{\chi}_{h})$
$\displaystyle:=-\langle\boldsymbol{w},\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{\chi}_{h})\boldsymbol{n}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}+\frac{\gamma_{a}}{h_{F}}\langle\boldsymbol{w},\boldsymbol{\chi}_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\,,$
(3.3a) $\displaystyle b_{\gamma}(q,\psi_{h})$ $\displaystyle:=-\langle
q,\boldsymbol{K}\nabla\psi_{h}\cdot\boldsymbol{n}\rangle_{\Gamma_{p}^{D}}+\frac{\gamma_{b}}{h_{F}}\langle
q,\psi_{h}\rangle_{\Gamma^{D}_{p}}\,.$ (3.3b)

The second of the options in (3.2g) amounts to a symmetric interior penalty
discontinuous Galerkin discretization of the scalar variable $p$; cf., e.g.,
[22, Sec. 4.2]. As usual, the average
$\\{\mskip-5.0mu\\{\cdot\\}\mskip-5.0mu\\}$ and jump $\lsem\cdot\rsem$ of a
function $w\in L^{2}(\Omega)$ on an interior face $F$ between two elements
$K^{+}$ and $K^{-}$, such that $F=\partial K^{+}\cap\partial K^{-}$, are

$\\{\mskip-5.0mu\\{w\\}\mskip-5.0mu\\}:=\frac{1}{2}(w^{+}+w^{-})\quad\text{and}\quad\lsem
w\rsem:=w^{+}-w^{-}\,.$

For boundary faces $F\subset\partial K\cap\partial\Omega$, we set
$\\{\mskip-5.0mu\\{w\\}\mskip-5.0mu\\}:=w_{|K}$ and $\lsem w\rsem:=w_{|K}$.
The set of all faces (interior and boundary faces) on the multigrid level
$\mathcal{T}_{l}$ is denoted by $\mathcal{F}_{h}$. In the second of the
options in (3.2g), the parameter $\gamma$ of the last term has to be chosen
sufficiently large, such that discrete coercivity on $Q_{h}$ of $B_{\gamma}$
is preserved. The local length $h_{F}$ is chosen as
$h_{F}=\\{\mskip-5.0mu\\{h_{F}\\}\mskip-5.0mu\\}:=\frac{1}{2}(|K^{+}|_{d}+|K^{-}|_{d})$
with Hausdorff measure $|\cdot|_{d}$; cf. [22, p. 125]. For boundary faces we
set $h_{F}:=|K|_{d}$. In (3.3), the quantities $\gamma_{a}$ and $\gamma_{b}$
are the algorithmic parameters of the stabilization terms in the Nitsche
formulation. To ensure well-posedness of the discrete systems the parameters
$\gamma_{a}$ and $\gamma_{b}$ have to be chosen sufficiently large; cf.
Appendix A. Based on our numerical experiments we choose the algorithmic
parameter $\gamma$, $\gamma_{a}$ and $\gamma_{b}$ in (3.2g) and (3.3) as

$\displaystyle\gamma_{a}=5\cdot 10^{4}\cdot
r\cdot(r+1)\quad\text{and}\quad\gamma=\gamma_{b}=\frac{1}{2}\cdot
r\cdot(r-1)\,,$

where $r$ is the polynomial degree of the finite element space (2.5a) for the
displacement variable.

Finally, for given $\boldsymbol{f}\in\boldsymbol{H}^{-1}(\Omega)$,
$\boldsymbol{u}_{D}\in\boldsymbol{H}^{1/2}(\Gamma^{D}_{\boldsymbol{u}})$,
$\boldsymbol{t}_{N}\in\boldsymbol{H}^{-1/2}(\Gamma^{N}_{\boldsymbol{u}})$ and
$g\in H^{-1}(\Omega)$, $p_{D}\in H^{1/2}(\Gamma^{D}_{p})$, $p_{N}\in
H^{-1/2}(\Gamma^{N}_{p})$ for $Q_{h}=Q_{h}^{l,\text{cont}}$, and properly
fitted assumptions about the data for $Q_{h}=Q_{h}^{l,\text{disc}}$, we put

$\displaystyle F_{\gamma}(\boldsymbol{\chi}_{h})$
$\displaystyle:=\langle\boldsymbol{f},\boldsymbol{\chi}_{h}\rangle-\langle\boldsymbol{t}_{N},\boldsymbol{\chi}_{h}\rangle_{\Gamma^{N}_{\boldsymbol{u}}}+a_{\gamma}(\boldsymbol{u}_{D},\boldsymbol{\chi}_{h})\,,$
(3.4a) $\displaystyle G_{\gamma}(\psi_{h})$
$\displaystyle:=\left\\{\begin{array}[]{@{}ll}\displaystyle\langle\langle
g,\psi_{h}\rangle-\alpha\langle\boldsymbol{v}_{D}\cdot\boldsymbol{n},\psi_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}-\langle
p_{N},\psi_{h}\rangle_{\Gamma^{N}_{p}}+b_{\gamma}(p_{D},\psi_{h})\,,&\text{for}\;\;Q_{h}=Q_{h}^{l,\text{cont}}\,,\\\\[12.91663pt]
\begin{array}[]{@{}l}\displaystyle\langle
g,\psi_{h}\rangle-\sum_{F\in\mathcal{F}_{h}^{D,\boldsymbol{u}}}\alpha\langle\boldsymbol{v}_{D}\cdot\boldsymbol{n},\psi_{h}\rangle_{F}-\sum_{F\in\mathcal{F}_{h}^{D,p}}\langle
p_{D},\\{\mskip-5.0mu\\{\boldsymbol{K}\nabla\psi_{h}\\}\mskip-5.0mu\\}\cdot\boldsymbol{n}\rangle_{F}\\\\[17.22217pt]
\displaystyle\;+\sum_{F\in\mathcal{F}_{h}^{D,p}}\frac{\gamma}{h}\langle
p_{D},\lsem\psi_{h}\rsem\rangle_{F}-\sum_{F\subset\mathcal{F}_{h}^{N,p}}\langle
p_{N},\lsem\psi_{h}\rsem\rangle_{F}\,,\end{array}&\text{for}\;\;Q_{h}=Q_{h}^{l,\text{disc}}\,.\end{array}\right.$
(3.4f)

In the second of the options in (3.4f), we denote by
$\mathcal{F}_{h}^{D,p}\subset\mathcal{F}_{h}$ and
$\mathcal{F}_{h}^{N,p}\subset\mathcal{F}_{h}$ the set of all element faces on
the boundary parts $\Gamma_{p}^{D}$ and $\Gamma_{p}^{N}$, respectively; cf.
(1.1). The second of the terms on the right-hand side of (3.4f), with
$\boldsymbol{v}_{D}=\partial_{t}\boldsymbol{u}_{D}$, is added to ensure
consistency of the form (3.2b) in the fully discrete formulation (3.5c) of
(1.1b), i.e., that the discrete equation (3.5c) is satisfied by the continuous
solution to (1.1).

We use a temporal test basis that is supported on the subintervals $I_{n}$;
cf., e.g., [2, 33]. Then, a time marching process is obtained. In that, we
assume that the trajectories $\boldsymbol{u}_{\tau,h}$,
$\boldsymbol{v}_{\tau,h}$ and $p_{\tau,h}$ have been computed before for all
$t\in[0,t_{n-1}]$, starting with approximations
$\boldsymbol{u}_{\tau,h}(t_{0}):=\boldsymbol{u}_{0,h}$,
$\boldsymbol{v}_{\tau,h}(t_{0}):=\boldsymbol{u}_{1,h}$ and
$p_{\tau,h}(t_{0}):=p_{0,h}$ of the initial values $\boldsymbol{u}_{0}$,
$\boldsymbol{v}_{0}$ and $p_{0,h}$. Then, we consider solving the following
local problem on $I_{n}$.

######  Problem 3.1 (Numerically integrated $I_{n}$-problem)

For given
$\boldsymbol{u}_{h}^{n-1}:=\boldsymbol{u}_{\tau,h}(t_{n-1})\in\boldsymbol{V}_{h}$,
$\boldsymbol{v}_{h}^{n-1}:=$
$\boldsymbol{v}_{\tau,h}(t_{n-1})\in\boldsymbol{V}_{h}$, and
$p_{h}^{n-1}:=p_{\tau,h}(t_{n-1})\in Q_{h}$ with
$\boldsymbol{u}_{\tau,h}(t_{0}):=\boldsymbol{u}_{0,h}$,
$\boldsymbol{v}_{\tau,h}(t_{0}):=\boldsymbol{u}_{1,h}$ and
$p_{\tau,h}(t_{0}):=p_{0,h}$, find
$(\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau,h},p_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$
such that

$\displaystyle\begin{aligned}
&Q_{n}\big{(}\langle\partial_{t}\boldsymbol{u}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle-\langle\boldsymbol{v}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle\big{)}+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle=\langle\boldsymbol{u}_{h}^{n-1},\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle\,,\\\\[4.30554pt]
\end{aligned}$ (3.5a) $\displaystyle\begin{aligned}
&Q_{n}\Big{(}\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{\chi}_{\tau,h}\rangle+A_{\gamma}(\boldsymbol{u}_{\tau,h},\boldsymbol{\chi}_{\tau,h})+C(\boldsymbol{\chi}_{\tau,h},p_{\tau,h})\Big{)}+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\chi}_{\tau,h}^{+}(t_{n-1})\rangle\\\\[4.30554pt]
&\quad=Q_{n}\Big{(}F_{\gamma}(\boldsymbol{\chi}_{\tau,h})\Big{)}+\langle\rho\boldsymbol{v}_{h}^{n-1},\chi_{\tau,h}^{+}(t_{n-1})\rangle\,,\\\\[4.30554pt]
\end{aligned}$ (3.5b) $\displaystyle\begin{aligned} &Q_{n}\Big{(}\langle
c_{0}\partial_{t}p_{\tau,h},\psi_{\tau,h}\rangle-C(\boldsymbol{v}_{\tau,h},\psi_{\tau,h})+B_{\gamma}(p_{\tau,h},\psi_{\tau,h})\Big{)}+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),\psi_{\tau,h}^{+}(t_{n-1})\rangle\\\\[4.30554pt]
&\quad=Q_{n}\Big{(}G_{\gamma}(\psi_{\tau,h})\Big{)}+\langle
c_{0}p_{h}^{n-1},\psi_{\tau,h}^{+}(t_{n-1})\rangle\end{aligned}$ (3.5c)

for all
$(\boldsymbol{\phi}_{\tau,h},\boldsymbol{\chi}_{\tau,h},\psi_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$.

The trajectories defined by Problem 3.1, for $n=1,\ldots,N$, satisfy that
$\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau,h}\in
Y_{\tau}^{k}(\boldsymbol{V}_{h})$ and $p_{\tau,h}\in Y_{\tau}^{k}(Q_{h})$. The
quadrature formulas on the left hand-side of (3.5) can be rewritten by time
integrals since the Gauss–Radau formula (2.4) is exact for all polynomials
$w\in P_{2k}(I_{n};\mathbb{R})$. Well-posedness of Problem 3.1 is ensured.

######  Lemma 3.2 (Existence and uniqueness of solutions to Problem 3.1)

Problem 3.1 admits a unique solution.

#### Proof

We prove Lem. 3.2 for $Q_{h}=Q_{h}^{l,\text{cont}}$ only, thus assuming the
first of the options in (3.2g) and (3.4f). For $Q_{h}=Q_{h}^{l,\text{disc}}$,
the proof can be done similarly by using, in addition, standard techniques of
error analysis for discontinuous Galerkin methods; cf., e.g., [22, Sec. 4].
Since Problem 3.5 is finite-dimensional, it suffices to prove uniqueness of
the solution. Existence of the solution then follows directly from its
uniqueness. Let now
$(\boldsymbol{u}^{(1)}_{\tau,h},\boldsymbol{v}^{(1)}_{\tau,h},p^{(1)}_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$
and
$(\boldsymbol{u}^{(2)}_{\tau,h},\boldsymbol{v}^{(2)}_{\tau,h},p^{(2)}_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$
denote two triples of solutions to (3.5a). Their differences
$(\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau,h},p_{\tau,h})=(\boldsymbol{u}^{(1)}_{\tau,h},\boldsymbol{v}^{(1)}_{\tau,h},p^{(1)}_{\tau,h})-(\boldsymbol{u}^{(2)}_{\tau,h},\boldsymbol{v}^{(2)}_{\tau,h},p^{(2)}_{\tau,h})$
then satisfy the equations

$\displaystyle
Q_{n}\big{(}\langle\partial_{t}\boldsymbol{u}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle-\langle\boldsymbol{v}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle\big{)}+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle=0\,,$
(3.6a) $\displaystyle
Q_{n}\Big{(}\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{\chi}_{\tau,h}\rangle+A_{\gamma}(\boldsymbol{u}_{\tau,h},\boldsymbol{\chi}_{\tau,h})+C(\boldsymbol{\chi}_{\tau,h},p_{\tau,h})\Big{)}+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\chi}_{\tau,h}^{+}(t_{n-1})\rangle=0\,,$
(3.6b) $\displaystyle Q_{n}\Big{(}\langle
c_{0}\partial_{t}p_{\tau,h},\psi_{\tau,h}\rangle-C(\boldsymbol{v}_{\tau,h},\psi_{\tau,h})+B_{\gamma}(p_{\tau,h},\psi_{\tau,h})\Big{)}+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),\psi_{\tau,h}^{+}(t_{n-1})\rangle=0$ (3.6c)

for all
$(\boldsymbol{\phi}_{\tau,h},\boldsymbol{\chi}_{\tau,h},\psi_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$.
We let $\boldsymbol{A}_{\gamma}:\boldsymbol{V}_{h}\mapsto\boldsymbol{V}_{h}$
be the discrete operator that is defined, for
$\boldsymbol{w}_{h}\in\boldsymbol{V}_{h}$ and all
$\boldsymbol{\phi}_{h}\in\boldsymbol{V}_{h}$, by

$\langle\boldsymbol{A}_{\gamma}\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}\rangle=A_{\gamma}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\,.$
(3.7)

In (3.6) we choose
$\boldsymbol{\phi}_{\tau,h}=\boldsymbol{A}_{\gamma}\boldsymbol{u}_{\tau,h}$,
$\boldsymbol{\chi}_{\tau,h}=\boldsymbol{v}_{\tau,h}$ and
$\psi_{\tau,h}=p_{\tau,h}$. Adding the resulting equations yields that

$\displaystyle
Q_{n}\big{(}\langle\partial_{t}\boldsymbol{u}_{\tau,h},\boldsymbol{A}_{\gamma}\boldsymbol{u}_{\tau,h}\rangle+\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{v}_{\tau,h}\rangle+\langle
c_{0}\partial_{t}p_{\tau,h},p_{\tau,h}\rangle+B_{\gamma}(p_{\tau,h},p_{\tau,h})\big{)}$
(3.8)
$\displaystyle\quad+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{A}_{\gamma}u_{\tau,h}^{+}(t_{n-1})\rangle+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{v}_{\tau,h}^{+}(t_{n-1})\rangle+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle=0\,.$

Recalling the exactness of the Gauss–Radau formula (2.4) for
$w\in\mathbb{P}_{2k}(I_{n};\mathbb{R})$, Eq. (3.8) yields that

$\displaystyle\frac{1}{2}\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}\big{(}\langle\boldsymbol{A}_{\gamma}\boldsymbol{u}_{\tau,h},\boldsymbol{u}_{\tau,h}\rangle+\langle\rho\boldsymbol{v}_{\tau,h},\boldsymbol{v}_{\tau,h}\rangle+\langle
c_{0}p_{\tau,h},p_{\tau,h}\rangle\big{)}\,\mathrm{d}t+Q_{n}\big{(}B_{\gamma}(p_{\tau,h},p_{\tau,h})\big{)}$
$\displaystyle\quad+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{A}_{\gamma}u_{\tau,h}^{+}(t_{n-1})\rangle+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{v}_{\tau,h}^{+}(t_{n-1})\rangle+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle=0\,.$

Using (3.7), this shows that

$\displaystyle
A_{\gamma}(\boldsymbol{u}_{\tau,h}(t_{n}),\boldsymbol{u}_{\tau,h}(t_{n}))+\langle\rho\boldsymbol{v}_{\tau,h}(t_{n}),\boldsymbol{v}_{\tau,h}(t_{n})\rangle+\langle
c_{0}p_{\tau,h}(t_{n}),p_{\tau,h}(t_{n})\rangle+2Q_{n}\big{(}B_{\gamma}(p_{\tau,h},p_{\tau,h})\big{)}$
(3.9)
$\displaystyle\quad+A_{\gamma}(\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{u}_{\tau,h}^{+}(t_{n-1}))+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{v}_{\tau,h}^{+}(t_{n-1})\rangle+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle=0\,.$

From (3.9) along with the discrete coercivity properties (A.3) of $A_{\gamma}$
and (A.5) of $B_{\gamma}$ we directly deduce that

$\boldsymbol{u}_{\tau,h}(t_{n})=\boldsymbol{u}_{\tau,h}^{+}(t_{n-1})=\boldsymbol{0}\,,\quad\boldsymbol{v}_{\tau,h}(t_{n})=\boldsymbol{v}_{\tau,h}^{+}(t_{n-1})=\boldsymbol{0}\quad\text{and}\quad
p_{\tau,h}(t_{n})=p_{\tau,h}^{+}(t_{n-1})=0\,,$ (3.10)

as well as

$p_{\tau,h}\big{(}t_{n,\mu}^{\text{GR}}\big{)}=0\,,\quad\text{for}\;\;\mu=1,\ldots,k+1\,.$
(3.11)

We note that $t_{n,k+1}^{\text{GR}}=t_{n}$. Relation (3.11) now implies that
$p_{\tau,h}\equiv 0$ on $I_{n}$.

To prove that $\boldsymbol{u}_{\tau,h}\equiv\boldsymbol{0}$ and
$\boldsymbol{v}_{\tau,h}\equiv\boldsymbol{0}$, by (3.10) it is sufficient to
show that $\boldsymbol{u}_{\tau,h}(t_{n,\mu}^{\text{G}})=\boldsymbol{0}$ and
$\boldsymbol{v}_{\tau,h}(t_{n,\mu}^{\text{G}})=\boldsymbol{0}$, for
$\mu=1,\ldots,k$, where $t_{n,\mu}^{\text{G}}$, for $\mu=1,\ldots,k$, are the
nodes of the $k$-point Gauss quadrature formula on $I_{n}$ that is exact for
all polynomials in $\mathbb{P}_{2k-1}(I_{n};\mathbb{R})$. This is done now.
Recalling (3.10), we conclude from (3.6a) by a suitable choice of test
functions that

$\partial_{t}\boldsymbol{u}_{\tau,h}(t_{n,\mu}^{\text{GR}})=\boldsymbol{v}_{\tau,h}(t_{n,\mu}^{\text{GR}})\,,\quad\text{for}\;\;\mu=1,\ldots,k\,.$
(3.12)

Next, choosing $\chi_{\tau,h}=\boldsymbol{v}_{\tau,h}$ in (3.6b) and recalling
(3.12) imply that

$Q_{n}\Big{(}\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{v}_{\tau,h}\rangle+A_{\gamma}(\boldsymbol{u}_{\tau,h},\partial_{t}\boldsymbol{u}_{\tau,h})\Big{)}=0$
(3.13)

By the exactness of the Gauss–Radau formula (2.4) for all
$w\in\mathbb{P}_{2k}(I_{n};\mathbb{R})$ we have from (3.13) that

$\int_{t_{n-1}}^{t_{n}}\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{v}_{\tau,h}\rangle\,\mathrm{d}t+\frac{1}{2}\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}A_{\gamma}(\boldsymbol{u}_{\tau,h},\boldsymbol{u}_{\tau,h})\,\mathrm{d}t=0\,.$
(3.14)

The second of the terms in (3.14) vanishes by (3.10). The stability result of
[42, Lem. 2.1] then implies that

$\boldsymbol{v}_{\tau,h}(t_{n,\mu}^{\text{G}})=\boldsymbol{0}\,,\quad\text{for}\;\;\mu=1,\ldots,k\,.$
(3.15)

From (3.15) along with (3.10) we then deduce that
$\boldsymbol{v}_{\tau,h}\equiv\boldsymbol{0}$ on $I_{n}$. Choosing the test
function $\boldsymbol{\phi}_{\tau,h}=\boldsymbol{u}_{\tau_{h}}$ in (3.6),
using $\boldsymbol{v}_{\tau,h}\equiv\boldsymbol{0}$ and applying the stability
result [42, Lem. 2.1], it follows that

$\boldsymbol{u}_{\tau,h}(t_{n,\mu}^{\text{G}})=\boldsymbol{0}\,,\quad\text{for}\;\;\mu=1,\ldots,k\,.$
(3.16)

From (3.16) along with (3.10) we then have that
$\boldsymbol{u}_{\tau,h}\equiv\boldsymbol{0}$ on $I_{n}$. Thus uniqueness of
solutions to (3.5) and, thereby, well-posedness of Problem 3.1 is thus
ensured.

$\blacksquare$

In Appendix B an alternative formulation for the system (3.6) is still
presented. It is based on using the time derivative
$\partial_{t}\boldsymbol{u}_{\tau,h}$ of the primal variable
$\boldsymbol{u}_{\tau,h}$ instead of the auxiliary variable
$\boldsymbol{v}_{\tau,h}$ in (3.5c). In this case, an additional surface
integral has to be included; cf. Eq. (B.1c) .

##  4 Algebraic solver by geometric multigrid preconditioned GMRES iterations

On the algebraic level, the variational problem (B.1) leads to linear systems
of equations with complex block structure, in particular if higher order
(piecewise) polynomial degrees $k$ for the approximation of the temporal
variable are involved. This demands for a robust and efficient linear solver,
in particular in the three-dimensional case $\Omega\subset\mathbb{R}^{3}$. For
solving (B.1) we consider using flexible GMRES iterations [56] that are
preconditioned by a $V$-cycle geometric multigrid method (GMG) based on a
local Vanka smoother. In [3], the combined GMG–GMRES approach is further
embedded in a Newton iteration for solving space-time finite element
discretizations of the Navier–Stokes system. Thus, nonlinear extensions of the
prototype model (1.1) are expected to become feasible by our approach as well.

To derive the algebraic form of (3.5), the discrete functions
$\boldsymbol{u}_{\tau,h}$, $\boldsymbol{v}_{\tau,h}$ and $p_{\tau,h}$ are
represented in a Lagrangian basis
$\\{\chi_{n,m}\\}_{m=1}^{k+1}\subset\mathbb{P}_{k}(I_{n};\mathbb{R})$ with
respect to the $(k+1)$ Gauss–Radau quadrature points of $I_{n}$, such that

$\displaystyle\boldsymbol{u}_{\tau,h}{}_{|I_{n}}(\boldsymbol{x},t)$
$\displaystyle=\sum_{m=1}^{k+1}\boldsymbol{u}_{n,m}(\boldsymbol{x})\chi_{n,m}(t)\,,\quad\boldsymbol{v}_{\tau,h}{}_{|I_{n}}(\boldsymbol{x},t)=\sum_{m=1}^{k+1}\boldsymbol{v}_{n,m}(\boldsymbol{x})\chi_{n,m}(t)\,,$
(4.1a) $\displaystyle p_{\tau,h}{}_{|I_{n}}(\boldsymbol{x},t)$
$\displaystyle=\sum_{m=1}^{k+1}p_{n,m}(\boldsymbol{x})\chi_{n,m}(t)\,.$ (4.1b)

The resulting coefficient functions
$(\boldsymbol{u}_{n,m},\boldsymbol{v}_{n,m},p_{n,m})\in\boldsymbol{V}_{h}\times\boldsymbol{V}_{h}\times
Q_{h}$, for $m=1,\ldots,k+1$, are developed in terms of the finite element
basis of $\boldsymbol{V}_{h}$ and $Q_{h}$, respectively. Letting
$\boldsymbol{V}_{h}=\operatorname{span}\\{\boldsymbol{\psi}_{1},\ldots,\boldsymbol{\psi}_{R}\\}$
and $Q_{h}=\operatorname{span}\\{\xi_{1},\ldots,\xi_{S}\\}$, we get that

$\boldsymbol{u}_{n,m}(\boldsymbol{x})=\sum_{r=1}^{R}u^{(r)}_{n,m}\boldsymbol{\psi}_{r}(\boldsymbol{x})\,,\;\;\boldsymbol{v}_{n,m}(\boldsymbol{x})=\sum_{r=1}^{R}v^{(r)}_{n,m}\boldsymbol{\psi}_{r}(\boldsymbol{x})\quad\text{and}\quad
p_{n,m}(\boldsymbol{x})=\sum_{s=1}^{S}p^{(s)}_{n,m}\,\xi_{s}(\boldsymbol{x})\,.$
(4.2)

For the coefficients of the expansions in (4.2) we define the subvectors

$\displaystyle\boldsymbol{U}_{n,m}$
$\displaystyle=\big{(}u^{(1)}_{n,m},\ldots,u^{(R)}_{n,m}\big{)}^{\top}\,,\;\;\boldsymbol{V}_{n,m}=\big{(}v^{(1)}_{n,m},\ldots,v^{(R)}_{n,m}\big{)}^{\top}\quad\text{and}$
(4.3a) $\displaystyle\boldsymbol{P}_{n,m}$
$\displaystyle=\big{(}p^{(1)}_{n,m},\ldots,p^{(S)}_{n,m}\big{)}^{\top}\,,\quad\text{for}\;\;m=1,\ldots,k+1\,,$
(4.3b)

of the degrees of freedom for all Gauss–Radau quadrature points and the global
solution vector on $I_{n}$ by

$\boldsymbol{X}_{n}^{\top}=\big{(}(\boldsymbol{V}_{n,1})^{\top},(\boldsymbol{U}_{n,1})^{\top},(\boldsymbol{P}_{n,1})^{\top},\ldots,(\boldsymbol{V}_{n,k+1})^{\top},(\boldsymbol{U}_{n,k+1})^{\top},(\boldsymbol{P}_{n,k+1})^{\top}\big{)}\,.$
(4.4)

We note that $\boldsymbol{X}_{n}$ comprises the (spatial) degrees of freedom
for all $(k+1)$ Gauss–Radau nodes, representing the Lagrange interpolation
points in time , of the subinterval $I_{n}$. The approximations at these time
points will be computed simultaneously. Substituting (4.1) and (4.2) into
(B.1) and choosing in (B.1) the test basis
$\\{\chi_{n,m}\boldsymbol{\psi}_{r},\chi_{n,m}\boldsymbol{\psi}_{r},\chi_{n,m}\xi_{s}\\}$,
for $m=1,\ldots,k+1$, $r=1,\ldots,R$ and $s=1,\ldots,S$, built from the trial
basis in (4.1) and (4.2), we obtain the following algebraic system.

######  Problem 4.1 (Algebraic $I_{n}$-problem)

For the vector $\boldsymbol{X}_{n}$, defined in (4.4) along with (4.3), of the
coefficients of the expansions (4.2) solve

$\boldsymbol{A}_{n}\boldsymbol{X}_{n}=\boldsymbol{F}_{n}\,,$ (4.5)

where the matrix $\boldsymbol{A}_{n}$ exhibits the $(k+1)\times(k+1)$ block
structure

$\boldsymbol{A}_{n}=\big{(}\boldsymbol{A}_{a,b}\big{)}_{a,b=1}^{k+1}$ (4.6)

with block submatrices $\boldsymbol{A}_{a,b}$, for $a,b=0,\ldots,k$, defined
by

$\boldsymbol{A}_{a,b}=\left(\begin{array}[]{@{}ccc@{}}-\boldsymbol{M}^{0,\boldsymbol{V}_{h}}_{{a},{b}}&\boldsymbol{M}^{1,\boldsymbol{V}_{h}}_{{a},{b}}&\boldsymbol{0}\\\\[4.30554pt]
\boldsymbol{M}^{1,\boldsymbol{V}_{h}}_{{a},{b}}&\boldsymbol{A}_{{a},{b}}+\boldsymbol{N}^{A}_{{a},{b}}&\boldsymbol{C}^{\top}_{{a},{b}}\\\\[4.30554pt]
\boldsymbol{0}&-\boldsymbol{C}_{{a},{b}}&\boldsymbol{M}^{1,Q_{h}}_{a,b}+\boldsymbol{B}_{a,b}+\boldsymbol{N}^{B}_{a,b}\end{array}\right)\,.$
(4.7)

For the choice $Q_{h}^{l}=Q_{h}^{l,\text{cont}}$ in (2.6), the explicit
representation of the submatrices in (4.7) reads as

$\displaystyle\big{(}\boldsymbol{M}^{1,\boldsymbol{V}_{h}}_{{a},{b}}\big{)}_{i,j}$
$\displaystyle:=Q_{n}\big{(}\langle\rho\partial_{t}\chi_{n,a}\boldsymbol{\psi}_{j},\chi_{n,b}\boldsymbol{\psi}_{i}\rangle\big{)}+\langle\rho\chi_{n,a}(t_{n-1}^{+})\boldsymbol{\psi}_{j},\chi_{n,b}(t_{n-1}^{+})\boldsymbol{\psi}_{i}\rangle\,,$
$\displaystyle\big{(}\boldsymbol{M}^{0,\boldsymbol{V}_{h}}_{{a},{b}}\big{)}_{i,j}$
$\displaystyle:=Q_{n}\big{(}\langle\rho\chi_{n,a}\boldsymbol{\psi}_{j},\chi_{n,b}\boldsymbol{\psi}_{i}\rangle\big{)}\,,$
$\displaystyle\big{(}\boldsymbol{A}_{{a},{b}}\big{)}_{i,j}$
$\displaystyle:=Q_{n}\big{(}\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\chi_{n,a}\boldsymbol{\psi}_{j}),\boldsymbol{\varepsilon}(\chi_{n,b}\boldsymbol{\psi}_{i})\rangle+\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\chi_{n,a}\boldsymbol{\psi}_{j})\boldsymbol{n},\chi_{n,b}\boldsymbol{\psi}_{i}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\big{)}\,,$
$\displaystyle\big{(}\boldsymbol{N}^{A}_{{a},\boldsymbol{\psi}_{b}}\big{)}_{i,j}$
$\displaystyle:=Q_{n}\big{(}a_{\gamma}(\chi_{n,a}\boldsymbol{\psi}_{i},\chi_{n,b}\boldsymbol{\psi}_{i})\big{)}$

as well as

$\displaystyle\big{(}\boldsymbol{C}_{{a},{b}}\big{)}_{r,j}$
$\displaystyle:=Q_{n}\big{(}-\alpha\langle\nabla\cdot(\chi_{n,b}\boldsymbol{\psi}_{j}),\chi_{n,a}\xi_{r}\rangle+\alpha\langle\chi_{n,b}\boldsymbol{\psi}_{j}\cdot\boldsymbol{n},\chi_{n,a}\xi_{r}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\big{)}\,,$
$\displaystyle\big{(}\boldsymbol{M}^{1,Q_{h}}_{p_{a},p_{b}}\big{)}_{r,s}$
$\displaystyle:=Q_{n}\big{(}\langle
c_{0}\partial_{t}\chi_{n,a}\xi_{s},\chi_{n,b}\xi_{r}\rangle\big{)}+\langle
c_{0}\chi_{n,a}(t_{n-1}^{+})\xi_{s},\chi_{n,b}(t_{n-1}^{+})\xi_{r}\rangle\big{)}$

and

$\displaystyle\big{(}\boldsymbol{B}_{{a},{b}}\big{)}_{r,s}$
$\displaystyle:=Q_{n}\big{(}\langle\boldsymbol{K}\nabla(\chi_{n,a}\xi_{s}),\nabla(\chi_{n,b}\xi_{r})\rangle-\langle\boldsymbol{K}\nabla(\chi_{n,a}\boldsymbol{\xi}_{s})\cdot\boldsymbol{n},\chi_{n,b}\xi_{r}\rangle_{\Gamma^{D}_{p}}\big{)}\,,$
$\displaystyle\big{(}\boldsymbol{N}^{B}_{{a},{b}}\big{)}_{r,s}$
$\displaystyle:=Q_{n}\big{(}b_{\gamma}(\chi_{n,a}\xi_{s},\chi_{n,b}\xi_{r})\big{)}\,,$

with $a_{\gamma}(\cdot,\cdot)$ and $b_{\gamma}(\cdot,\cdot)$ being defined in
(3.3) and $i,j=1,\ldots,R$ and $r,s=1,\ldots,S$. The right-hand side vector
$\boldsymbol{F}_{n}$ in (4.5) is defined similarly, according to (B.1) along
with (3.4). Its definition is skipped here for brevity. We note that (3.1a) is
still multiplied by $\rho>0$ for the definition of the first row in (4.7). For
the family of finite element pairs, corresponding to
$Q_{h}^{l}=Q_{h}^{l,\text{disc}}$ in (2.6), the definition of the submatrix
$B_{a,b}$ has to be adjusted to the second of the options in (3.2g). Further,
the contribution $\boldsymbol{N}^{B}_{a,b}$ has to be omitted.

Increasing values of the piecewise polynomial degree in time $k$ enhance the
complexity of the block structure of the system matrix $\boldsymbol{A}_{n}$ in
(4.6) along with (4.7). They impose an additional facet of challenge on the
construction of efficient block preconditioners for (4.5). We solve the linear
system (4.5) for the unknown $\boldsymbol{X}_{n}$ on the subinterval $I_{n}$
by flexible GMRES iterations [56] that are preconditioned by a $V$-cycle
geometric multigrid (GMG) algorithm [62]. Key ingredients of the GMRES–GMG
approach are

  * •

left preconditioning of each GMRES iteration by a single $V$-cycle GMG
iteration,

  * •

prolongation to grids by interpolation and application of the transpose
operator for restriction to grids,

  * •

$J_{\max}$ pre- and post-smoothing steps on each grid level
$\\{\mathcal{T}_{l}\\}$ for $l=1,\ldots,L$ and

  * •

application of a (parallel) direct solver on the coarsest mesh partition
$\mathcal{T}_{0}$.

In our computational studies presented in Sec. 5, we use $J_{\max}=4$ and the
parallel direct solver SuperLU_DIST [44]. For the implementation of the
restriction and prolongation operators the deal.II classes MultiGrid and
MGTransferPrebuilt are used. For the deal.II finite element library we refer
to [6]. For details of the parallel implementation by the message passing
interface (MPI) protocol of our GMG approach we refer to [2].

The choice of the smoothing operator in the GMG method is extremely diverse;
cf., e.g., [23, 39] for a further discussion. We use a collective smoother of
Vanka type that is based on the solution of small local problems. Compared to
the Navier–Stokes system [2], the dynamic Biot problem (1.1) requires
adaptations in the construction of the local Vanka smoother. In inf-sup stable
discretizations of the Navier–Stokes equations by the
$\mathbb{Q}_{r}^{d}/\mathbb{P}_{r-1}^{\text{disc}}$, $r\geq 2$, family of
finite element spaces no coupling between the pressure degrees of freedom is
involved due to the discontinuous approximation of the pressure variable. This
feature leads to excellent performance properties of the Vanka smoother [2].
In contrast to this, the discretization (3.5c) of (1.1b) by
$\mathbb{P}_{r-1}^{\text{disc}}$, $r\geq 2$, elements involves a coupling
between degrees of freedom of the scalar variable $p$ due to the presence of
the face integrals over the average and jump in the second of the options in
(3.2g) for the definition of the bilinear form $B_{\gamma}$. This coupling
reduces the smoothing properties of elementwise Vanka operator. For the
$\mathbb{Q}_{r-1}$, $r\geq 2$, family of elements for the scalar variable $p$,
leading to the first of the options in (3.2g), the coupling of degrees of
freedom of $p$ by its continuous in space approximation and the loss of
smoothing properties arise likewise. As a remedy, for both families of inf-sup
stable approximation in (3.5) the local Vanka smoother is computed on
overlapping patches of adjacent elements. In addition, an averaging of the
patchwise updates is done after the Vanka smoother has been applied on all of
them.

To construct the patchwise Vanka, we let on the mesh partition
$\mathcal{T}_{l}$, for $l=1,\ldots,L$, of the multigrid hierarchy the linear
system, to be solved, be represented by

$\boldsymbol{A}_{l}\boldsymbol{d}_{l}=\boldsymbol{b}_{l}\,,\quad\text{for
}\;\;l=1,\ldots,L\,.$ (4.8)

To each grid node $\boldsymbol{\xi}_{l}^{m}$, for $m=1\,,\ldots\,,M_{l}$,
where $M_{l}$ denotes the total number of grid nodes on the mesh partition
$\mathcal{T}_{l}$, we built a patch of adjacent elements such that

$P_{l}^{m}:=\bigcup\\{K\in\mathcal{T}_{l}\mid\xi_{l}^{m}\in\overline{K}\\}\,,\quad\text{for}\;\;m=1,\ldots,M_{l}\,.$
(4.9)

In two space dimensions $P_{l}^{m}$ is built from four elements, if
$\xi_{l}^{m}\not\in\partial\Omega$. In three space dimensions, $P_{l}^{m}$ has
eight elements, if $\xi_{l}^{m}\not\in\partial\Omega$. If
$\xi_{l}^{m}\in\partial\Omega$, patches of less elements are obtained. On
$\mathcal{T}_{l}$, let $Z_{l}$ denote the index set of all global degrees of
freedom with cardinality $C_{l}$,

$C_{l}:=\operatorname{card}(Z_{l})\,.$

Let $Z_{l}(P_{l}^{m})$ denote the subset of $Z_{l}$ all global degrees of
freedom linked to the patch $P_{l}^{m}$, i.e. degrees of freedom of
$\boldsymbol{u}_{\tau,h}$, $\boldsymbol{v}_{\tau,h}$ and $p_{\tau,h}$ for the
$(k+1)$ Gauss–Radau time points of $I_{n}$. The cardinality of
$Z_{l}(P_{l}^{m})$ is denoted by $C_{l}^{m}$,

$C_{l}^{m}:=\operatorname{card}(Z_{l}(P_{m}^{l}))\,,\quad\text{for}\;\;m=1,\ldots,M_{l}\,.$

Further, we denote the index set of all local degrees of freedom on
$P_{l}^{m}$ by $\hat{Z}_{l}(P_{m}^{l}):=\\{0,\ldots,C_{l}^{m}-1\\}$. For the
notation, we note that the set $\hat{Z}_{l}(P_{m}^{l})$ depends on the
cardinality of the patch $P_{l}^{m}$. For a given patch $P_{l}^{m}$ and a
local degree of freedom with number $\hat{\mu}\in\hat{Z}_{l}(P_{m}^{l})$ let
the mapping

$\operatorname{dof}:\mathcal{T}_{l}\times\hat{Z}_{l}(P_{m}^{l})\rightarrow
Z_{l}\,,\qquad\mu=\mathrm{dof}(P_{l}^{m},\hat{\mu})\in Z_{l}(P_{l}^{m})\,,$
(4.10)

yield the uniquely defined global number $\mu\in Z_{l}$. Finally, we put
$R=\operatorname{dim}\,\boldsymbol{V}_{h}$ and $S=\operatorname{dim}\,Q_{h}$.

We are now in a position to define the local Vanka operator for the patch
$P_{l}^{m}$ in an exact mathematical way; cf., e.g., [2, 32, 37, 66].

######  Definition 4.2 (Patchwise Vanka smoother)

On patch $P_{l}^{m}$, for $m\in\\{1,\ldots,M_{l}\\}$, the local Vanka operator
$\boldsymbol{S}_{P_{l}^{m}}:\mathbb{R}^{(k+1)\,\cdot\,(2R+S)}\rightarrow\mathbb{R}^{C^{l}_{m}}$
is defined by

$\boldsymbol{S}_{P_{l}^{m}}(\boldsymbol{d})=\boldsymbol{R}_{P_{l}^{m}}\boldsymbol{d}+\omega\,\boldsymbol{A}_{P_{l}^{m}}^{-1}\,\boldsymbol{R}_{P_{l}^{m}}\,(\boldsymbol{b}_{l}-\boldsymbol{A}_{l}\boldsymbol{d})\,,$
(4.11)

with underrelaxation factor $\omega>0$ and $P_{l}^{m}$-local restriction
operator
$\boldsymbol{R}_{K}:\mathbb{R}^{(k+1)\cdot(2R+S)}\rightarrow\mathbb{R}^{C_{l}^{m}}$
defined by

$(\boldsymbol{R}_{P_{l}^{m}}\boldsymbol{d})[\hat{\mu}]=\boldsymbol{d}[\operatorname{dof}(P_{l}^{m},\hat{\mu})]\,,\quad\text{for}\;\;\hat{\mu}\in\hat{Z}_{l}(P_{l}^{m})\,.$
(4.12)

Further, the patch system matrix
$\boldsymbol{A}_{P_{l}^{m}}\in\mathbb{R}^{C_{l}^{m},C_{l}^{m}}$ in (4.11) is
defined by

$\boldsymbol{A}_{P_{l}^{m}}[\hat{\nu}][\hat{\mu}]:=\boldsymbol{A}[\operatorname{dof}(P_{l}^{m},\hat{\nu})][\operatorname{dof}(P_{l}^{m},\hat{\mu})]\,,\quad\text{for}\;\;\hat{\nu},\hat{\mu}\in\hat{Z}_{l}(P_{l}^{m})\,.$
(4.13)

In the numerical experiments of Sec. 5 we choose $\omega=0.7$. The
$P_{l}^{m}$-local restriction operator (4.12) assigns to a global defect
vector $\boldsymbol{d}\in\mathbb{R}^{(k+1)\cdot(2R+S)}$ the local block vector
$\boldsymbol{R}_{K}\boldsymbol{d}\in\mathbb{R}^{C_{l}^{m}}$ that contains all
components of $\boldsymbol{d}$ that are associated with all degrees of freedom
(for all $(k+1)$ Gauss–Radau points of $I_{n}$) belonging to the patch
$P_{l}^{m}$. For the computation of the inverse
$(\boldsymbol{A}_{P_{l}^{m}})^{-1}$ of the patch system
matrix$\boldsymbol{A}_{P_{l}^{m}}$ in (4.11) we use LAPACK routines.

The application of the smoother for (4.8) on the mesh partition
$\mathcal{T}_{l}$ is summarized in Algorithm 4.1.

1 Initialize solution $\boldsymbol{d}$ of (4.8) with $\boldsymbol{0}$;

2 for  _$j=1,\ldots,J_{\max}$_ do

3 for  _( $\mu=0$; $\mu<C_{l}$; $\mu\text{++})$_ do

4 $p[\mu]=0$;

5 $z[\mu]=0.0$;

6

7 end for

8 foreach  _patch $P_{l}^{m}$ on $\mathcal{T}_{l}$_ do

9 $\boldsymbol{y}=\boldsymbol{S}_{P_{l}^{m}}\boldsymbol{(}\boldsymbol{d})$;

10 for  _( $\hat{\mu}=0$, $\hat{\mu}<C_{l}^{m}$; $\hat{\mu}\text{++}$)_ do

11
$z[\operatorname{dof}(P_{l}^{m},\hat{\mu})]\text{+=}(\boldsymbol{E}_{P_{l}^{m}}\boldsymbol{y})[{\operatorname{dof}(P_{l}^{m},\hat{\mu})}]$;

12 $p[\operatorname{dof}(P_{l}^{m},\hat{\mu})]\text{++}$ ;

13

14 end for

15

16 end foreach

17 for  _( $\mu=0$; $\mu<C_{l}$; $\mu\text{++})$_ do

18 $d[\mu]=z[\mu]/p[\mu]$;

19

20 end for

21

22 end for

Algorithm 4.1 Smoothing steps for (4.8) on mesh partition $\mathcal{T}_{l}$.

In line 4.1 of Alg. 4.1 the defect and solution vector is pre-initialized with
$\boldsymbol{0}$. In line 4.1 the loop over all $J_{\max}$ smoothing steps
starts. In line 4.1 and 4.1 the counter vector $\boldsymbol{p}$ for the number
of updates of the degrees of freedom and the auxiliary vector $\boldsymbol{z}$
are initialized with $\boldsymbol{0}$. In line 4.1 the loop over all patches
$P_{l}^{m}$, for $m=1,\ldots,M$, starts. In line 4.1 the local Vanka smoother
is applied on patch $P_{l}^{m}$ to the current iterate $\boldsymbol{d}$ of the
defect vector and the image is stored in the local patch vector
$\boldsymbol{y}$. In line 4.1 the local vector $\boldsymbol{y}$ is assigned to
an auxiliary global vector $\boldsymbol{z}$. In line 4.1 the components of the
counter $\boldsymbol{p}$ are incremented for the (global) indices associated
with the degrees of freedom of the patch $P_{l}^{m}$ processed in the loop.
Finally, in line 4.1 the arithmetic mean of the local (patchwise) updates
$\boldsymbol{z}$ is assigned to the update of the global defect vector
$\boldsymbol{d}$. In line 4.1, the $P_{l}^{m}$-dependent extension operator
$\boldsymbol{E}_{P_{l}^{m}}:\mathbb{R}^{C_{l}^{m}}\rightarrow\mathbb{R}^{(k+1)\cdot(2R+S)}$
is defined by

$(\boldsymbol{E}_{P_{l}^{m}}\boldsymbol{y})[\mu]=\left\\{\begin{array}[]{@{}ll}y[{\hat{\mu}}]\,,&\text{if}\;\exists\hat{\mu}\in\hat{Z}_{l}(P_{l}^{m}):\;\mu=\operatorname{dof}(K,\hat{\mu})\,,\\\\[4.30554pt]
0\,,&\text{if}\;\mu\not\in Z_{l}(P_{l}^{m})\,.\end{array}\right.$

Regarding the performance of Alg. 4.1 and the overall GMRES–GMG approach we
note the following.

###### Remark 4.3

  * •

Averaging of the patchwise updates implemented in line 4.1 and 4.1 of Alg.
4.1, that is used instead of overwriting successively the (global) degrees of
freedom within the patch loop starting in line 4.1, is essential and ensures
the convergence and efficiency of the local Vanka smoother and, thereby, the
desired performance of the overall GMRES–GMG linear solver. Without the
averaging step we encountered convergence problems of the GMRES–GMG solver for
the experiments of Sec. 5.

  * •

The same holds for the application of the Vanka smoother on the patches (4.9)
instead of using an elementwise Vanka smoother on the element $K$. The latter
would lead to systems (4.11) of smaller dimension, however causes convergence
problems for the GMRES–GMG solver. This is expected to be due coupling of the
degrees of freedom of the scalar variable $p$ in the spatial discretizations
used here,.

##  5 Numerical studies

In this section we study numerically the proposed space-time finite element
and GMRES–GMG solver approach with respect to its computational and energy
efficiency. Firstly, we demonstrate the accuracy of solutions in terms of
convergence rates. Further, we illustrate the robustness of the GMRES–GMG
solver with respect to the resolution in space and time and the polynomial
approximation order. This is done for a challenging three-dimensional test
problem of practical interest, in order to demonstrate the feasibility of
large scale three-dimensional simulations. The numerical convergence of goal
quantities of the benchmark is addressed. The parallel scaling properties of
our implementation are also investigated. Beyond these studies of classical
performance engineering, the energy efficiency of the approach is considered
further.

The implementation of the numerical scheme and the GMRES–GMG solver was done
in an in-house high-performance frontend solver for the deal.II library [6].
For details of the parallel implementation of the geometric multigrid solver
we refer to [2]. In all numerical experiments, the stopping criterion for the
GMRES iterations is an absolute residual smaller than 1e-8. The computations
were performed on a Linux cluster with 571 nodes, each of them with 2 CPUs and
36 cores per CPU. The CPUs are _Intel Xeon Platinum 8360Y_ with a base
frequency of $2.4\text{\,}\mathrm{GHz}$, a maximum turbo frequency of
$3.5\text{\,}\mathrm{GHz}$ and a level 3 cache of $54\text{\,}\mathrm{MB}$.
Each node has $252\text{\,}\mathrm{GB}$ of main memory.

###  5.1 Accuracy of the discretization: Experimental order of convergence

We study (1.1) for $\Omega=(0,1)^{2}$ and $I=(1,2]$ and the prescribed
solution

$\boldsymbol{u}(\boldsymbol{x},t)=\phi(\boldsymbol{x},t)\boldsymbol{E}_{2}\;\;\text{and}\;\;p(\boldsymbol{x},t)=\phi(\boldsymbol{x},t)\;\;\text{with}\;\;\phi(\boldsymbol{x},t)=\sin(\omega_{1}t^{2})\sin(\omega_{2}x_{1})\sin(\omega_{2}x_{2})$
(5.1)

and $\omega_{1}=\omega_{2}=\pi$. We put $\rho=1.0$, $\alpha=0.9$, $c_{0}=0.01$
and $\boldsymbol{K}=\boldsymbol{E}_{2}$ with the identity
$\boldsymbol{E}_{2}\in\mathbb{R}^{2,2}$. For the fourth order elasticity
tensor $\boldsymbol{C}$, isotropic material properties with Young’s modulus
$E=100$ and Poisson’s ratio $\nu=0.35$, corresponding to the Lamé parameters
$\lambda=86.4$ and $\mu=37.0$, are chosen. For an experiment with larger
values of $\lambda$ and $\mu$ we refer to Table C.2 in the appendix. In our
numerical experiments, the norm of $L^{\infty}(I;L^{2})$ is approximated by
computing the function’s maximum value in the Gauss quadrature nodes $t_{n,m}$
of $I_{n}$, i.e.,

$\|w\|_{L^{\infty}(I;L^{2})}\approx\max\\{\|w_{|I_{n}}(t_{n,m})\|\mid
m=1,\ldots,M\,,\;n=1,\ldots,N\\}\,,\quad\text{with}\;\;M=100\,.$

We study the space-time convergence behavior of the scheme (3.5). For this,
the domain $\Omega$ is decomposed into a sequence of successively refined
meshes of quadrilateral finite elements. The spatial and temporal mesh sizes
are halved in each of the refinement steps. The step sizes of the coarsest
space and time mesh are $h_{0}=1/(2\sqrt{2})$ and $\tau_{0}=0.1$. We choose
the polynomial degree $k=2$ and $r=3$, such that discrete solutions
$\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau}\in
Y_{\tau}^{2}(\boldsymbol{V}_{h})$, $p_{\tau,h}\in Y_{\tau}^{2}(Q_{h})$ with
local spaces $\mathbb{Q}_{3}^{2}/P_{2}^{\text{disc}}$ are obtained, as well as
$k=3$ and $r=4$ with $\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau}\in
Y_{\tau}^{3}(\boldsymbol{V}_{h})$ and $p_{\tau,h}\in Y_{\tau}^{3}(Q_{h})$ and
local spaces $\mathbb{Q}_{4}^{2}/\mathbb{P}_{3}^{\text{disc}}$ ; cf. (2.2),
(2.5) and (2.6). The calculated errors and corresponding experimental orders
of convergence are summarized in Table 5.1 and 5.2, respectively. The error is
measured in the quantities associated with the energy of the system (1.1); cf.
[35, p. 15] and [10]. Table 5.1 and 5.2 nicely confirm the optimal rates of
convergence with respect to the polynomial degrees in space and time of the
overall approach. A notable difference in the convergence behavior between the
pairs $\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$ and
$\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$ of local spaces for the discretization
of the spatial variables is not observed. For completeness, we summarize in
Appendix C the convergence results obtained for the pair
$\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$ of spaces of the Taylor–Hood family. A
minor superiority of the pair
$\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$ over the pair
$\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$ is only seen in the approximation of the
scalar variable $p$. The coincidence of the convergence results also holds for
the application of Problem B.1 instead of Problem 3.1. The two discrete
problems differ from each other by the discretization of the term
$\alpha\nabla\cdot\partial_{t}\boldsymbol{u}$ in (1.1b).

| $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC |  $\|p-p_{\tau,h}\|_{L^{2}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 1.2544218392e-02 | – | 3.4897282317e-02 | – | 2.4070118274e-03 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.5227995262e-03 | 3.04 | 3.9246006564e-03 | 3.15 | 2.8841669021e-04 | 3.06
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 1.8904870171e-04 | 3.01 | 4.8175203148e-04 | 3.03 | 3.5986044195e-05 | 3.00
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 2.3544020929e-05 | 3.01 | 5.9913957453e-05 | 3.01 | 4.4985485037e-06 | 3.00
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.9374127099e-06 | 3.00 | 7.4788750858e-06 | 3.00 | 5.6221925444e-07 | 3.00
| $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{\infty}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 8.1391652440e-02 | – | 1.6415428887e-01 | – | 1.0828068294e-02 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.1474006637e-02 | 2.83 | 2.8042308570e-02 | 2.55 | 1.7292127803e-03 | 2.65
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 1.5142126866e-03 | 2.92 | 3.6968566445e-03 | 2.92 | 2.3686467755e-04 | 2.87
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 1.9337108028e-04 | 2.97 | 4.6572185327e-04 | 2.99 | 3.0596110688e-05 | 2.95
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.4392244385e-05 | 2.99 | 5.8210228994e-05 | 3.00 | 3.8749300434e-06 | 2.98
Table 5.1: $L^{2}(L^{2})$ and $L^{\infty}(L^{2})$ errors and experimental orders of convergence (EOC) for (5.1) with temporal polynomial degree $k=2$ and spatial degree $r=3$ for local spaces $\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$.  | $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{2}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 2.3958455291e-03 | – | 1.7185653242e-02 | – | 7.0604463908e-04 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.0529085600e-04 | 4.51 | 5.4558622568e-04 | 4.98 | 3.4799927094e-05 | 4.34
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 5.9749777345e-06 | 4.14 | 1.9143193313e-05 | 4.83 | 1.2074793168e-06 | 4.85
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 3.6858164230e-07 | 4.02 | 9.6894400389e-07 | 4.30 | 6.1656313807e-08 | 4.29
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.2976498320e-08 | 4.00 | 5.6982903695e-08 | 4.09 | 3.6688337244e-09 | 4.07
| $\tau$ | $h$ |  $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{\infty}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 1.3014842843e-02 | – | 1.1587851244e-01 | – | 4.9827180852e-03 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 8.7944275936e-04 | 3.89 | 4.1645282584e-03 | 4.80 | 2.6180886102e-04 | 4.25
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 5.6086043381e-05 | 3.97 | 2.0121085807e-04 | 4.37 | 9.4019466765e-06 | 4.80
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 3.5171465193e-06 | 4.00 | 1.1873567590e-05 | 4.08 | 5.0100645993e-07 | 4.23
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.1959929657e-07 | 4.00 | 7.3034599295e-07 | 4.02 | 3.0246305328e-08 | 4.05
Table 5.2: $L^{2}(L^{2})$ and $L^{\infty}(L^{2})$ errors and experimental
orders of convergence (EOC) with temporal polynomial degree $k=3$ and spatial
degree $r=4$ for local spaces
$\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$.

Next, we show numerically that the time discretization is even superconvergent
of order $2k+1$ in the discrete time nodes $t_{n}$, for $n=1,\ldots,N$. For
this, we introduce the time mesh dependent norm

$\|w\|_{l^{\infty}(L^{2})}:=\max\\{\|w(t_{n})\|\mid n=1,\ldots,N\\}\,.$ (5.2)

We prescribe the solution

$\boldsymbol{u}(\boldsymbol{x},t)=\begin{pmatrix}-2(x-1)^{2}x^{2}(y-1)y(2y-1)\sin(\omega_{1}t)\\\
\phantom{-}2(x-1)x(2x-1)(y-1)^{2}y^{2}\sin(\omega_{1}t)\end{pmatrix}\;\;\text{and}\;\;p(\boldsymbol{x},t)=-2(x-1)^{2}x^{2}(y-1)y(2y-1)\sin(\omega_{2}t)\;\;$
(5.3)

with $\omega_{1}=40\cdot\pi$ and $\omega_{2}=10\cdot\pi$. We put $\rho=1.0$,
$\alpha=0.9$, $c_{0}=0.01$ and $\boldsymbol{K}=\boldsymbol{E}_{2}$ with the
identity $\boldsymbol{E}_{2}\in\mathbb{R}^{2,2}$. For the elasticity tensor
$\boldsymbol{C}$, isotropic material properties with Young’s modulus $E=100$
and Poisson’s ratio $\nu=0.35$ are used. For the local spaces we choose the
pair $\mathbb{Q}_{5}^{2}/\mathbb{Q}_{4}$ such that, for any $t\in[0,T]$, the
solution (5.3) belongs to the discrete spaces $\boldsymbol{V}_{h}$ and
$Q_{h}$, respectively, and its spatial approximation is exact. This
simplification is done here since we aim to study the convergence of the
temporal discretization only. In the experiment, we choose $k=2$ such that
discrete solutions $\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau}\in
Y_{\tau}^{2}(\boldsymbol{V}_{h})$ and $p_{\tau,h}\in X_{\tau}^{2}(Q_{h})$ are
obtained. We use a spatial mesh that consists of 16 cells with
$h=1/(2\sqrt{2})$ and set $\tau_{0}=0.02$. The calculated errors and
corresponding experimental orders of convergence are summarized in Table 5.3.
Superconvergence of order $2k+1$ in the discrete time nodes is clearly
observed in the second of the tables in Table 5.3.

| $\tau$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{2}(L^{2})}$ | EOC
---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | 3.3600958819e-03 | – | 5.2037881317e-02 | – | 8.5492528946e-05 | –
$\tau_{0}/2^{1}$ | 3.9552579243e-04 | 3.09 | 6.6291344116e-03 | 2.97 | 4.8852611508e-06 | 4.13
$\tau_{0}/2^{2}$ | 4.8913600397e-05 | 3.02 | 8.3404665964e-04 | 2.99 | 1.7945937694e-07 | 4.77
$\tau_{0}/2^{3}$ | 6.1123116212e-06 | 3.00 | 1.0446618892e-04 | 3.00 | 1.0554654134e-08 | 4.09
$\tau_{0}/2^{4}$ | 7.6411713541e-07 | 3.00 | 1.3065222159e-05 | 3.00 | 1.1603768561e-09 | 3.19
$\tau_{0}/2^{5}$ | 9.5518196224e-08 | 3.00 | 1.6333729266e-06 | 3.00 | 1.4370494853e-10 | 3.01
$\tau_{0}/2^{6}$ | 1.1939894351e-08 | 3.00 | 2.0417851944e-07 | 3.00 | 1.7952564965e-11 | 3.00
| $\tau$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{l^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{l^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{l^{\infty}(L^{2})}$ | EOC
---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | 2.1302171198e-03 | – | 1.3451755096e-02 | – | 1.9152238265e-04 | –
$\tau_{0}/2^{1}$ | 1.1003144753e-04 | 4.28 | 6.9871664304e-04 | 4.27 | 1.0401014806e-05 | 4.20
$\tau_{0}/2^{2}$ | 3.6416424630e-06 | 4.92 | 2.3145015656e-05 | 4.92 | 3.5819518738e-07 | 4.86
$\tau_{0}/2^{3}$ | 1.1784707720e-07 | 4.95 | 7.3594357548e-07 | 4.97 | 1.1336271728e-08 | 4.98
$\tau_{0}/2^{4}$ | 3.6966657143e-09 | 4.99 | 2.3096831287e-08 | 4.99 | 3.5849645590e-10 | 4.98
$\tau_{0}/2^{5}$ | 1.1560104141e-10 | 5.00 | 7.2294482894e-10 | 5.00 | 1.1225015472e-11 | 5.00
$\tau_{0}/2^{6}$ | 3.5626018095e-12 | 5.02 | 2.2323871214e-11 | 5.02 | 3.4673336542e-13 | 5.02
Table 5.3: $L^{2}(L^{2})$ and $l^{\infty}(L^{2})$ errors and experimental
orders of convergence (EOC) for (5.3) with temporal polynomial degree $k=2$
and spatial degree $r=5$ for local spaces
$\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$, showing superconvergence of order
$2k+1$ in the discrete time nodes $t_{n}$, i.e., w.r.t. the norm
$\|\cdot\|_{l^{\infty}(L^{2})}$ defined in (5.2).

###  5.2 Computational efficiency: Accuracy of goal quantities and convergence
of the GMRES–GMG solver in a 3d benchmark

For a three-dimensional problem configuration of practical interest (cf. [11,
Sec. 5.1]), we study the potential of the proposed approach to compute
reliably and efficiently goal quantities of physical interest that are defined
in (5.6). Concurrently, we analyze the convergence and performance properties
of the GMRES–GMG solver of Sec. 4 for our families of space-time finite
element methods. More precisely, we consider the L-shaped domain
$\Omega\subset\mathbb{R}^{3}$ sketched in Fig. 5.1 along with the boundary
conditions for $\boldsymbol{u}$ prescribed on the different parts of
$\partial\Omega$. Beyond the boundary conditions (1.1d) and (1.1e), we further
apply (homogeneous) directional boundary conditions on parts of the boundary
for the variable $\boldsymbol{u}$. This is done for the sake of physical
realism. For the portion $\Gamma^{d}_{\boldsymbol{u}}\subset\partial\Omega$
they read as

$\boldsymbol{u}\cdot\boldsymbol{n}=0\qquad\text{and}\qquad(\boldsymbol{\sigma}(\boldsymbol{u})\boldsymbol{n})\cdot\boldsymbol{t}_{i}=0\,,\quad\text{for}\;\;i=1,\ldots,d-1\,,\quad\text{on}\;\;\Gamma^{d}_{\boldsymbol{u}}\times(0,T]\,,$
(5.4)

for the stress tensor
$\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{u})$
and the unit basis vectors $\boldsymbol{t}_{i}$, for $i=1,\ldots,d-1$, of the
tangent space at $\boldsymbol{x}\in\Gamma^{d}_{\boldsymbol{u}}$. In the
definition of $A_{\gamma}$ in (3.2a) and (3.3a), the conditions (5.4) still
need to be implemented properly. For instance, for the left boundary of the
L-shaped domain $\Omega$ the first of the equations in (5.4) implies that

$\boldsymbol{u}_{\tau,h}\cdot\boldsymbol{n}=u_{\tau,h,1}=0\,,\quad\text{for}\;\;\boldsymbol{u}_{\tau,h}=(u_{\tau,h,1},u_{\tau,h,2})^{\top}\,.$

By the second of the conditions in (5.4) we get for the second of the terms on
the right-hand of (3.2a) that

$\langle\boldsymbol{C}\varepsilon(\boldsymbol{u}_{\tau,h})\boldsymbol{n},\boldsymbol{\chi}_{\tau,h}\rangle_{\Gamma_{\boldsymbol{u}}^{l}}=\langle\boldsymbol{C}\varepsilon(\boldsymbol{u}_{\tau,h})\boldsymbol{n}\cdot\boldsymbol{n},\boldsymbol{n}\cdot\boldsymbol{\chi}_{\tau,h}\rangle_{\Gamma_{\boldsymbol{u}}^{d}}\,.$
(5.5)

For the boundary part ${\Gamma^{d}_{\boldsymbol{u}}}$ we then put, similarly
to (3.3a),

$a^{d}_{\gamma}(\boldsymbol{w},\boldsymbol{\chi}_{h}):=-\langle\boldsymbol{w}\cdot\boldsymbol{n},\boldsymbol{C}\varepsilon(\boldsymbol{\chi}_{h})\boldsymbol{n}\cdot\boldsymbol{n}\rangle_{{\Gamma^{d}_{\boldsymbol{u}}}}+\gamma_{a}\,h^{-1}\langle\boldsymbol{w}\cdot\boldsymbol{n},\boldsymbol{\chi}_{h}\cdot\boldsymbol{n}\rangle_{\Gamma^{d}_{\boldsymbol{u}}}\,,$

while leaving $a_{\gamma}$ unmodified for the part
${\Gamma^{D}_{\boldsymbol{u}}}$ where Dirichlet boundary conditions are
prescribed for $\boldsymbol{u}$. In its entirety, we thus have that
$a_{\gamma}(\cdot,\cdot):=a_{\gamma}^{D}(\cdot,\cdot)+a_{\gamma}^{d}(\cdot,\cdot)$
with $a_{\gamma}^{D}(\cdot,\cdot)$ being defined by the right-hand side of
(3.3a). By the arguments of Appendix A, the coercivity of the resulting form
$A_{\gamma}$ is still ensured. On the right boundary we prescribe a
homogeneous Neumann condition for $\boldsymbol{u}$ with
$\boldsymbol{t}_{N}=\boldsymbol{0}$ in (1.1e). On the left part of the upper
boundary, i.e. for
$\Gamma^{N}_{\boldsymbol{u}}:=(0,0.5)\times\\{1\\}\times(0,0.5)$, we impose
the traction force (cf. (1.1e))

$\boldsymbol{t}_{N}=\begin{pmatrix}0\\\ 5\cdot
10^{9}\cdot(32xz-18x-16z+10)\sin(8\pi t)\\\ 0\end{pmatrix}\,.$

For the scalar variable $p$ we prescribe a homogeneous Dirichlet condition
(1.1f) on the left upper part of $\partial\Omega$, i.e., for
$(x_{1},x_{2},x_{3})\in\Gamma_{p}^{D}:=[0,0.5]\times\\{1\\}\times[0,0.5]$. On
$\Gamma_{p}^{N}:=\partial\Omega\backslash\Gamma_{p}^{D}$ a homogeneous Neumann
condition with $p_{N}=0$ in (1.1g) is prescribed. We put $\rho=1.0$,
$\alpha=0.9$, $c_{0}=0.01$ and $\boldsymbol{K}=\boldsymbol{E}_{3}$ with the
identity $\boldsymbol{E}_{3}\in\mathbb{R}^{3,3}$. For the elasticity tensor
$\boldsymbol{C}$, isotropic material properties with Young’s modulus $E=20000$
and Poisson’s ratio $\nu=0.3$, that model the material behavior of human bone
(cf. [67]), are used. We put $T=4$.

(a) Test setting: Cross-section plot of $\Omega$ in $(x_{1},x_{2})$-plane. The
3d domain extends the cross-section to $x_{3}\in[0,0.5]$.

(b) Modulus $|\boldsymbol{u}|$ of $\boldsymbol{u}$ at time t = 3.8.

(c) Variable $p$ at time t = 3.8.

Figure 5.1: Problem setting with boundary conditions for $\boldsymbol{u}$ and
profile of the solution at time t = 3.8.

As goal quantities of this benchmark problem, we measure the displacement in
the normal direction as well as the pressure at the right part of the boundary
by

$b_{\boldsymbol{u}}=\int_{\Gamma_{m}}\boldsymbol{u}\cdot\boldsymbol{n}\,\,\mathrm{d}o\qquad\text{and}\qquad
b_{p}=\int_{\Gamma_{m}}p\,\,\mathrm{d}o\qquad\text{for}\;\;\Gamma_{m}:=\\{1\\}\times(0,0.5)\times(0,0.5)\,.$
(5.6)

The time step size of the temporal discretization is chosen as $\tau=0.01$.
The calculated profile for $t=3.8$ of the modulus of the vectorial variable
$\boldsymbol{u}$ and the scalar variable $p$ are illustrated in Fig. 5.1. In
Table 5.4 we summarize characteristic quantities and results of our
simulations for various spatial resolutions and different temporal and spatial
polynomial degrees of the STFEMs. In Fig. 5.2 we visualize the computed
benchmark quantities (5.6) of some of the simulations over the temporal axis.
For the temporal discretization with polynomial degree $k=0$ and $k=1$ the
computed quantities differ on the spatial refinement level 4 only by a small
shift of the phases of the goal quantifies. This indicates that for the
current benchmark setting (choice of the model parameter and boundary
conditions) the spatial discretization error dominates the temporal one.
Therefore, in order to study numerically the convergence and accuracy of the
goal quantities, we computed the benchmark problem for higher spatial
resolutions (spatial mesh and approximation order; cf. Table 5.4. The lowest
row in Table 5.4 contains the results of the finest simulation that we could
run on our hardware. Fig. 5.2 shows that the solution (i.e., the goal
quantities) is not fully converged yet. This underlines the complexity of the
benchmark. The final column of Table 5.4 shows the convergence statistics of
the proposed GMRES–GMG solver. In terms of the average number of iterations
per time step the solver is (almost) grid independent. This underlines its
capability and robustness for solving efficiently the complex systems arising
from space-time finite element discretizations of partial differential
equations.

Table 5.4: Computed goal quantities (5.6) for different temporal and spatial approximations. Here, $k$ and $r$ are the polynomial degrees (cf. Sec. 2), refn is the number $n$ of refinements of the initial spatial grid (with $(k+1)\cdot 390$ degrees of freedom for $r=2$) and DoF${}_{I}{{}_{n}}$ is the total number of degrees of freedom on the fine grid level $L$ of GMG per subinterval $I_{n}$ (i.e. for all $k+1$ Gauss–Radau quadrature points on $I_{n}$). The coarse grid level of GMG is ref1 (with $(k+1)\cdot 2013$ degrees of freedom for $r=2$). Further, the minimum and maximum of the benchmark quantities $b_{\boldsymbol{u}}$ and $b_{p}$ as well as the average number $\bar{n}_{\text{GMRES}}$ of GMRES iterations per time step are summarized. $k$ |  $r$ | refn | DoF${}_{I_{n}}$ | $b_{\boldsymbol{u}_{\text{min}}}$ | $b_{\boldsymbol{u}_{\text{max}}}$ | $b_{p_{\text{min}}}$ | $b_{p_{\text{max}}}$ | $\bar{n}_{\text{GMRES}}$
---|---|---|---|---|---|---|---|---
0 | 2 | 4 | $682\,950$ | $-7.65210-04$ | $7.62710-04$ | $-33.330$ | $33.125$ | $9$
1 | 2 | 2 | $25\,836$ | $-1.67410-02$ | $1.67610-02$ | $-942.600$ | $945.243$ | $11$
1 | 2 | 3 | $182\,220$ | $-3.26310-03$ | $3.25110-03$ | $-165.600$ | $166.085$ | $12$
1 | 2 | 4 | $1\,365\,900$ | $-7.30510-04$ | $7.33510-04$ | $-34.200$ | $34.292$ | $14$
1 | 2 | 5 | $10\,571\,532$ | $-1.81210-04$ | $1.81210-04$ | $-7.836$ | $7.857$ | $12$
1 | 2 | 6 | $83\,172\,876$ | $-4.78910-05$ | $4.80610-05$ | $-1.892$ | $1.897$ | $12$
1 | 3 | 5 | $34\,596\,492$ | $-4.51810-05$ | $4.51810-05$ | $-1.874$ | $1.879$ | $11$
1 | 3 | 6 | $273\,640\,716$ | $-1.28610-05$ | $1.28010-05$ | $-0.461$ | $0.459$ | $8$

(a) Benchmark quantity $b_{\boldsymbol{u}}$ for refinement levels ref4 to ref6
and $t\in[0,T]$.

(b) Benchmark quantity $b_{p}$ for refinement levels ref4 to ref6 and
$t\in[0,T]$.

(c) Benchmark quantity $b_{\boldsymbol{u}}$ for refinement levels ref5 and
ref6 and $t\in[0,1]$.

(d) Benchmark quantity $b_{p}$ for refinement levels ref5 and ref6 and
$t\in[0,1]$.

Figure 5.2: Goal quantities $b_{\boldsymbol{u}}$ and $b_{p}$ of (5.6) for
different discretizations (spatial mesh resolution and polynomial degrees).

###  5.3 Parallel scaling and energy efficiency

Here, we study briefly the parallel scaling and energy efficiency of our
solver. By studying energy consumption, we’d like to draw attention to this
emerging dimension in the tuning of algorithms. Energy efficiency broadens the
classical hardware-oriented numerics that is applied to enhance the
performance of the current method on the target hardware and/or to find other
numerical methods to improve the numerical efficiency. For the longer term,
energy and power consumption needs to be mapped into a rigorous performance
model. Here, we restrict ourselves to illustrate numerically the parallel
scaling and energy consumption properties of our implementation that uses
Message Passing Interface(MPI) libraries and multi threading parallelism.

We perform a strong scaling benchmark for the test problem of Subsec. 5.2 with
$k=1$, $r=2$ and refn = 5, with $10\,571\,532$ degrees of freedom in each
subinterval $I_{n}$ on the fine level $\mathcal{T}_{L}$ (cf. Table 5.4) and
$98\,304$ mesh cells on $\mathcal{T}_{L}$. Throughout, we assign 18 MPI
processes to each of the nodes used for the computations and vary the number
$n$ of nodes from $n=20$ to $n=200$. For the evaluation of the parallel
scaling properties, we compute the parallel speedup of the code (cf. [1]) that
is approximated by

$S=\frac{t_{\text{wall}}(n=n_{\min})}{t_{\text{wall}}(n)}\,,$ (5.7)

where $t_{\text{wall}}(n)$ denotes the wall time of the simulation of fixed
size on $n$ compute nodes and $t_{\text{wall}}(n=n_{\min})$ is the wall time
of the simulation on the minimum number of nodes involved in the scaling
experiment. Secondly, we compute similarly the energy ratio by means of

$R=\frac{E(n)}{E(n=n_{\min})}\,,$ (5.8)

where $E(n)$ measures the total energy consumption of the simulation on $n$
nodes. The energy consumption is determined by the Linux cluster workload
manager slurm [59]. The energy consumption data is collected from hardware
sensors using Intel’s Running Average Power Limit (RAPL) mechanism. It
measures the energy consumption of the processor and memory. On our system,
the sampling interval of energy consumption is determined by the value of 30
seconds. Fig. 5.3 illustrates the results of the performance test.

(a) Wall time and ideal wall time for increased number of nodes.

(b) Parallel speedup $S$ and energy ratio $R$ for increased number of nodes.

Figure 5.3: Results of the strong scaling and energy consumption benchmark.

To quantify and evaluate the productivity or resource costs of the algorithm
and its implementation, we use a simple model for the $\text{Productivity
}P=\frac{\text{Output }O}{\text{Input }I}$; cf. [21]. In an economic sense,
all outputs should be the desired ones. Therefore, we use the reciprocal of
the wall time $t_{\text{wall}}$ as the output $O=\frac{1}{t_{\text{wall}}}$,
such that a decrease in $t_{\text{wall}}$ represents an increase of the
(abstract) output. As the input $I$ we use the total energy consumption $E$ of
the simulation. We scale the result by multiplying $P$ with the constant
factor $E(n_{\min})\cdot t_{\text{wall}}(n_{\min})$ such that the computation
with $n=n_{\min}$ has a productivity of $P=1.0$:

$P=\frac{\frac{1}{t_{\text{wall}}(n)}}{E(n)}\cdot E(n_{\min})\cdot
t_{\text{wall}}(n_{\min})=\frac{S}{R}\,,$ (5.9)

with $S$ and $R$ being defined in (5.7) and (5.8), respectively. The resulting
productivity curve of our computations is presented in Fig. 5.4. In our model
the simulation on 80 compute nodes is the most productive one, that is, the
ratio of output (low wall-time) to input (energy) is best. The characteristic
quantities of our performance study are also summarized in Table 5.5.

Figure 5.4: Piecewise quadratic interpolation of the productivity function (5.9). Table 5.5: Computed quantities of the scaling and performance benchmark. $n$ |  $t_{\text{wall}}$ [h] | E [kWh] | $S$ | $R$ | $P$
---|---|---|---|---|---
20 | $12.00$ | $98.31$ | 1.00 | 1.00 | 1.00
40 | $7.22$ | $133.49$ | 1.66 | 1.36 | 1.22
80 | $4.87$ | $173.50$ | 2.47 | 1.76 | 1.40
120 | $4.23$ | $220.61$ | 2.83 | 2.24 | 1.26
160 | $3.98$ | $272.48$ | 3.02 | 2.77 | 1.09
200 | $3.72$ | $316.34$ | 3.22 | 3.22 | 1.00

##  6 Summary and outlook

In this work we presented and analyzed families of space-time finite element
discretizations of the coupled hyperbolic-parabolic system (1.1) of poro- or
thermoelasticity. The time discretization uses the discontinuous Galerkin
method. The space discretization is based on inf-sup stable pairs of finite
element spaces with continuous and discontinuous approximation of the scalar
variable $p$. Well-posedness of the discrete problems is proved. For
efficiently solving the arising algebraic systems with complex block structure
in the case of increasing polynomial degrees of the time discretization a
geometric multigrid preconditioner with a local Vanka smoother on patches of
finite elements is proposed and studied. The overall approach is carefully
evaluated for a sequence of numerical experiments, in particular for a
challenging three-dimensional test setting of practical interest. Parallel
scaling and energy consumption is also considered in the experiments. The
benchmark problem of Subsec. 5.2 deserves further studies. In particular, soft
materials with a dominating impact of the hyperbolic equation (1.1a) to the
overall system (1.1) need to be studied as well. The definition of the goal
quantities of the proposed benchmark might be reconsidered and revised with
regard to either an improved choice of such quantities or their computation
with higher rate of convergence. This evaluation remains as a work for the
future. Finally, multi-field formulations of (1.1) with an explicit
approximation of the stress tensor
$\boldsymbol{\sigma}=\boldsymbol{C}\boldsymbol{\varepsilon}$ and the flux
vector $\boldsymbol{q}=-\boldsymbol{K}\nabla p$ might be advantageous for
applications of (1.1) in that these quantities are of interest or desired. By
using the concept of the local (patchwise) Vanka smoother GMG preconditioner,
that was presented here, multi-field formulations of (1.1) are expected to be
become feasible as well. This also remains as a work for the future.

## Acknowledgements

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

## References

  * [1] G. M. Amdahl, Validity of the single processor approach to achieving large scale computing capabilities, in Proceedings of the April 18-20, 1967, Spring Joint Computer Conference on — AFIPS ’67 (Spring), Atlantic City, New Jersey: ACM Press, 1967, p. 483.
  * [2] M. Anselmann, M. Bause, A geometric multigrid method for space-time finite element discretizations of the Navier–Stokes equations and its application to 3d flow simulation, ACM Trans. Math. Softw., accepted (2023), https://doi.org/10.1145/3582492, pp. 1–27; arXiv:2107.10561.
  * [3] M. Anselmann, M. Bause, Efficiency of local Vanka smoother geometric multigrid preconditioning for space-time finite element methods to the Navier-Stokes equations, PAMM Proc. Appl. Math. Mech., accepted (2022), pp. 1–6; arXiv:2210.02690.
  * [4] M. Anselmann, M. Bause, CutFEM and ghost stabilization techniques for higher order space-time discretizations of the Navier–Stokes equations, Int. J. Numer. Meth. Fluids, 94 (2022), pp. 775–802.
  * [5] M. Anselmann, M. Bause, S. Becher, G. Matthies, Galerkin-collocation approximation in time for the wave equation and its post-processing, ESAIM Math. Model. Numer. Anal., 54 (2020), pp. 2099–2123.
  * [6] D. Arndt, W. Bangerth, B. Blais, M. Fehling, R. Gassmöller, T. Heister, L. Heltai, U. Köcher, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, S. Proell, K. Simon, B. Turcksin, D. Wells, J. Zhang, The deal.II Library, Version 9.3, J. Numer. Math., 29 (2021), pp. 171–186.
  * [7] D. N. Arnold, D. Boffi, R. S. Falk, Approximation by quadrilateral finite elements, Math. Comp., 71 (2002), pp. 909–922.
  * [8] W. Bangerth, M. Geiger, R. Rannacher, Adaptive Galerkin finite element methods for the wave equation, Comput. Meth. Appl. Math., 10 (2010), pp. 3–48.
  * [9] W. Bangerth, R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhäuser, Basel, 2003.
  * [10] M. Bause, M. Anselmann, U. Köcher, F. A. Radu, Convergence of a continuous Galerkin method for hyperbolic-parabolic systems, Comput. Math. with Appl., submitted (2022), pp. 1–24; arXiv:2201.12014.
  * [11] M. Bause, R. Radu, U. Köcher, Space-time finite element approximation of the Biot poroelasticity system with iterative coupling, Comput. Methods Appl. Mech. Engrg., 320 (2017), pp. 745–768.
  * [12] R. Becker, Mesh adaptation for Dirirchlet flow control via Nitsche’s method, Commun. Numer. Meth. Engrg., 18 (2002), pp. 669–680.
  * [13] M. Biot, General theory of three-dimensional consolidation, J. Appl. Phys., 12 (1941), pp. 155–164.
  * [14] M. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys., 26 (1955), pp. 182–185.
  * [15] M. Biot, Theory of finite deformations of porous solids, Indiana Univ. Math. J., 21 (1972), pp. 597–620.
  * [16] J. W. Both, N. A. Barnafi, F. A. Radu, P. Zunino, A. Quarteroni, Iterative splitting schemes for a soft material poromechanics model, Comput. Methods Appl. Mech. Engrg., 388 (2022), 114183.
  * [17] A. Brandt, O. E. Livne, Multigrid Techniques—1984 Guide with Applications to Fluid Dynamics, SIAM, Philadelphia, 2011.
  * [18] S. C. Brenner, Korn’s inequalities for piecewise $\boldsymbol{H}^{1}$ vector fields, Math. Comp., 73 (2003), pp. 1067–1087.
  * [19] M. P. Bruchhäuser, U. Köcher, M. Bause, On the implementation of an adaptive multirate framework for coupled transport and flow, J. Sci. Comput., 93:59 (2022), https://doi.org/10.1007/s10915-022-02026-z, pp. 1–29.
  * [20] D. E. Carlson, Linear thermoelasticity, Handbuch der Physik V Ia/2, Springer, Berlin, 1972.
  * [21] W. W. Cooper, L. M. Seiford, K. Tone Data Envelopment Analysis, Kluwer Academic Publishers, Dordrecht, 2000.
  * [22] D. A. Di Pietro, A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Springer, Berlin, 2012.
  * [23] D. Drzisga, L. John, U. Rüde, B. Wohlmuth, W. Zulehner, On the analysis of block smoothers for saddle point problems, SIAM J. Matrix Anal. Appl., 39 (2018), pp. 932–960.
  * [24] G. Fairweather, Finite Element Galerkin Methods for Differential Equations, Lecture Notes in Pure and Applied Mathematics 34, Marcel Dekker Inc., New York, 1978.
  * [25] S. Franz, S. Trostorff, M. Waurick, Numerical methods for changing type systems, IMA J. Numer. Anal., 39 (2019), pp. 1009–1038.
  * [26] B. Gmeiner, M. Huber, L. John, U. Rüde, B. Wohlmuth, A quantitative performance study for Stokes solvers at the extreme scale, J. Comput. Sci., 17 (2016), pp. 509–521.
  * [27] B. Gmeiner, U. Rüde, H. Stengel, C. Waluga, B. Wohlmuth, Performance and scalability of hierarchical hybrid multigrid solvers for Stokes systems, SIAM J. Sci. Comput., 37 (2015), pp. C143–C168.
  * [28] B. Gmeiner, U. Rüde, H. Stengel, C. Waluga, B. Wohlmuth, Towards textbook efficiency for parallel multigrid, Numer. Math. Theory Methods Appl., 8 (2015), pp. 22–46.
  * [29] W. Hackbusch, Multigrid methods and applications, Springer, Berlin, 1985.
  * [30] Q. Hong, J. Kraus, _Parameter-robust stability of classical three-field formulation of Biot’s consolidation model_ , Electron. Trans. Numer. Anal., 48 (2018), pp. 202–226.
  * [31] Q. Hong, J. Kraus, J. Xu, L. Zikatanov, A robust multigrid method for discontinuous Galerkin discretizations of Stokes and linear elasticity equations, Numer. Math., 132 (2016), pp. 23–49.
  * [32] S. Hussain, F. Schieweck, S. Turek, Efficient Newton-multigrid solution techniques for higher order space-time Galerkin discretizations of incompressible flow, Appl. Numer. Math., 83 (2014), pp. 51–71.
  * [33] S. Hussain, F. Schieweck, S. Turek, An efficient and stable finite element solver of higher order in space and time for nonstationary incompressible flow, Internat. J. Numer. Methods Fluids, 73 (2013), pp. 927–952.
  * [34] S. Hussain, F. Schieweck, S. Turek, Higher order Galerkin time discretizations and fast multigrid solvers for the heat equation, J. Numer. Math., 19 (2011), pp. 41–61.
  * [35] S. Jiang, R. Racke, Evolution equations in thermoelasticity, CRC Press, Boca Raton, 2018.
  * [36] V. John, Finite Element Methods for Incompressible Flow Problems, Springer, Cham, 2016.
  * [37] V. John, Higher order finite element methods and multigrid solvers in a benchmark problem for the 3D Navier-Stokes equations, Int. J. Numer. Meth. Fluids, 40 (2002), pp. 775–798.
  * [38] V. John, G. Matthies, Higher-order finite element discretizations in a benchmark problem for incompressible flows, Int. J. Numer. Meth. Fluids, 37 (2001), pp. 885–903.
  * [39] V. John, L. Tobiska, Numerical performance of smoothers in coupled multigrid methods for the parallel solution of the incompressible Navier–Stokes equations, Int. J. Numer. Meth. Fluids, 33 (2000), pp. 453–473.
  * [40] G. Kanschat, B. Riviere, A finite element method with strong mass conservation for Biot’s linear consolidation model, J. Sci. Comput., 77 (2018), pp. 1762–1779.
  * [41] O. Karakashian, C. Makridakis, Convergence of a continuous Galerkin method with mesh modification for nonlinear wave equations, Math. Comp., 74 (2004), pp. 85–102.
  * [42] O. Karakashian, C. Makridakis, A space-time finite element method for the nonlinear Schrödinger equation: The continuous Galerkin method, SIAM J. Numer. Anal., 36 (1999), pp. 1779–1807.
  * [43] R. Leis, Initial boundary value problems in mathematical physics, Teubner, Stuttgart, John Wiley & Sons, Chichester, 1986.
  * [44] X. S. Li, J. W. Demmel, SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems, ACM Trans. Math. Softw., 29 (2003), pp. 110–140.
  * [45] J. Linden, G. Lonsdale, H. Ritzdorf, A. Schüller, Scalability aspects of parallel multigrid, Future Gener. Comput. Syst., 10 (1994), pp. 429–439.
  * [46] S. Manservisi, Numerical analysis of Vanka-type solvers for steady Stokes and Navier–Stokes flows, SIAM. J. Numer. Anal., 44 (2006), pp. 2025–2056.
  * [47] G. Matthies, Mapped finite elements on hexahedra. Necessary and sufficient conditions for optimal interpolation errors, Numer. Algor., 27 (2001), pp. 317–327.
  * [48] G. Matthies, L. Tobiska, The inf-sup condition for the mapped $\mathbb{Q}_{k}^{d}/P_{k-1}^{\operatorname{disc}}$ element in arbitrary space dimensions, Computing, 69 (2002), pp. 119–139.
  * [49] A. Mikelić, M. F. Wheeler, Theory of the dynamic Biot–Allard equations and their link to the quasi-static Biot system, J. Math. Phys., 53 (2012), 123702:1–15.
  * [50] J. Nitsche, _Über ein Variationsprinzip zur Lösung von Dirichlet Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind (in German)_ , Abh. Math. Sem. Univ. Hamburg, 36 (1971), pp. 9–15.
  * [51] P. J. Philips, M. F. Wheeler, A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I, II, Comput. Geosci., 11 (2007), 131–158.
  * [52] P. J. Philips, M. F. Wheeler, A coupling of mixed and discontinuous Galerkin finite element methods for poroelasticity , Comput. Geosci., 12 (2008), 417–435.
  * [53] C. Rodrigo, X. Hu, P. Ohm, J. H. Adler, F. J. Gaspar, L. T. Zikatanov, New stabilized discretizations for poroelasticity and the Stokes’ equations, Comput. Methods Appl. Mech. Engrg., 341 (2018), pp. 467–484.
  * [54] U. Rüde, Algorithmic Efficiency and the Energy Wall, 2017. DOI: 10.13140/RG.2.2.33914.18881.
  * [55] C. Seifert, S. Trostorff, M. Waurick, Evolutionary Equations: Picard’s Theorem for Partial Differential Equations, and Applications, Birkhäuser, Cham, 2022.
  * [56] Y. Saad, Iterative Methods for Sparse Linear Systems. 2nd edition. SIAM, Philadelphia, 2003.
  * [57] R. Showalter, Diffusion in poro-elastic media, J. Math. Anal. Appl., 251 (2000), pp. 310–340.
  * [58] M. Slodička, Application of Rothe’s method to integrodifferential equation, Comment. Math. Univ. Carolinae, 30 (1989), pp. 57–70.
  * [59] Slurm Workload Manager, Version 22.05, https://github.com/SchedMD/slurm/tree/master/src/plugins/acct_gather_energy/rapl, 2023.
  * [60] O. Steinbach, M. Zank, A generalized inf–sup stable variational formulation for the wave equation, J. Math. Anal. Appl., 505 (2022), 125457.
  * [61] V. Thomeé, Galerkin Finite Element Methods for Parabolic Problems, Springer, Berlin, 2006.
  * [62] U. Trottenberg, C. W. Oosterlee, A. Schüller, Multigrid, Academic Press, San Diego, CA, 2001.
  * [63] S. Turek, Efficient Solvers for Incompressible Flow Problems, Springer, Berlin, 1999.
  * [64] S. Turek, C. Becker, D. Kilian, Hardware-oriented numerics and concepts for PDE software, Future Generation Computer Systems, 22 (2006), pp. 217–238.
  * [65] S. Turek, D. Göddecke, C. Becker, S. Buijssen, H. Wobker, FEAST — realisation of hardware-oriented numerics for HPC simulations with finite elements, Concurrency and Computation: Practice and Experience 6 (May), pp. 2247–2265. Special Issue Proceedings of ISC 2008, doi:10.1002/cpe.1584.
  * [66] S. Vanka, Block-implicit multigrid solution of Navier-–Stokes equations in primitive variables, J. Comput. Phys., 65 (1986), pp. 138–158.
  * [67] D. C. Wirtz, N. Schiffers, T. Pandorf, K. Rademacher, D. Weichert, R. Forst, Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur, J. Biomech., 33 (2000), pp. 1325–1330.

## Appendix

##  A Discrete coercivity

Here, we prove the discrete coercivity of the bilinear form $A_{\gamma}$
defined in (3.2a) along with (3.3a). For this we introduce the mesh- and
parameter-dependent norm

$\|\boldsymbol{v}_{h}\|_{h,\widetilde{\gamma}_{a}}:=\Big{(}\|\varepsilon(\boldsymbol{v}_{h})\|^{2}+\widetilde{\gamma}_{a}\,h^{-1}\|\boldsymbol{v}_{h}\|_{\Gamma^{D}_{\boldsymbol{u}}}^{2}\Big{)}^{1/2}\,,\quad\text{for}\;\;\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}\,,$
(A.1)

with some parameter $\widetilde{\gamma}_{a}>0$. Its choice is addressed below.
The norm property of (A.1) is ensured by a variant of Korn’s inequality; cf.
[18, Eq. (1.19)]. We recall the well-known inverse inequality (cf. [61, p.
28])

$\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\boldsymbol{n}\|_{\Gamma^{D}_{\boldsymbol{u}}}\leq
ch^{-1/2}\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|\,,\quad\text{for}\;\;\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}\,.$
(A.2)

From (3.2a) along with (3.3a), it follows by the positive definiteness (1.2a)
of $\boldsymbol{C}$, the inequality of Cauchy–Young with a sufficiently small
constant $\delta>0$ and the trace inequality (A.2) that, for all
$\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$,

$\displaystyle A_{\gamma}(\boldsymbol{v}_{h},\boldsymbol{v}_{h})$
$\displaystyle=\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{v}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\rangle-\langle\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\boldsymbol{n},\boldsymbol{v}_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}-\langle\boldsymbol{v}_{h},\boldsymbol{C}\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\boldsymbol{n}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}+\gamma_{a}\,h^{-1}\langle\boldsymbol{v}_{h},\boldsymbol{v}_{h}\rangle_{\Gamma^{D}_{\boldsymbol{u}}}$
$\displaystyle\geq
c\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|^{2}-c\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\boldsymbol{n}\|_{\Gamma^{D}_{\boldsymbol{u}}}\|\boldsymbol{v}_{h}\|_{\Gamma^{D}_{\boldsymbol{u}}}+\gamma_{a}\,h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}$
$\displaystyle\geq
c\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|^{2}-c\delta\,h\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\boldsymbol{n}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}-c\delta^{-1}h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}+\gamma_{a}\,h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}$
$\displaystyle\geq
c\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|^{2}-c\delta\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|^{2}-c\delta^{-1}h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}+\gamma_{a}\,h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}$
$\displaystyle=c\Big{(}\|\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|^{2}+\widetilde{\gamma}_{a}\,h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{\boldsymbol{u}}}\Big{)}\,,$

with some constant $\widetilde{\gamma}_{a}>0$ for a sufficiently large choice
of the algorithmic parameter $\gamma_{a}$ in (3.3a), such that
$\widetilde{\gamma}_{a}:=\gamma_{a}-c\delta^{-1}>0$. Thus, there holds for
some constant $c>0$ that

$A_{\gamma}(\boldsymbol{v}_{h},\boldsymbol{v}_{h})\geq
c\|\boldsymbol{v}_{h}\|_{h,\widetilde{\gamma}_{a}}\,,\quad\text{for
all}\;\;\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}\,.$ (A.3)

Secondly, we prove the discrete coercivity of the bilinear form $B_{\gamma}$
defined in (3.2g) along with (3.3b). For brevity, we study the case
$Q_{h}=Q_{h}^{l,\text{cont}}$ only, such that $Q_{h}\subset H^{1}(\Omega)$ is
satisfied. The case $Q_{h}=Q_{h}^{l,\text{disc}}$ can be captured similarly.
We introduce the mesh- and parameter-dependent norm

$\|q_{h}\|_{h,\widetilde{\gamma}_{b}}:=\Big{(}\|\nabla
q_{h}\|^{2}+\widetilde{\gamma}_{b}\,h^{-1}\|q_{h}\|_{\Gamma^{D}_{p}}^{2}\Big{)}^{1/2}\,,\quad\text{for}\;\;q_{h}\in
Q_{h}\subset H^{1}(\Omega)\,,$

with some parameter $\widetilde{\gamma}_{b}>0$, where
$Q_{h}=Q_{h}^{l,\text{cont}}$; cf. (2.5). Further, we recall the well-known
inverse inequality

$\|\nabla q_{h}\cdot\boldsymbol{n}\|_{\Gamma^{D}_{p}}\leq ch^{-1/2}\|\nabla
q_{h}\|\,,\quad\text{for}\;\;q_{h}\in Q_{h}\,.$ (A.4)

From (3.2g) along with (3.3b), it follows by the positive definiteness (1.2b)
of $\boldsymbol{K}$, the inequality of Cauchy–Young with a sufficiently small
constant $\delta>0$ and the trace inequality (A.4) that, for all $q_{h}\in
Q_{h}$,

$\displaystyle B_{\gamma}(q_{h},q_{h})$
$\displaystyle=\langle\boldsymbol{K}\nabla q_{h},\nabla
q_{h}\rangle-\langle\boldsymbol{K}\boldsymbol{\nabla}q_{h}\cdot\boldsymbol{n},q_{h}\rangle_{\Gamma^{D}_{p}}-\langle
q_{h},\boldsymbol{K}\nabla
q_{h}\cdot\boldsymbol{n}\rangle_{\Gamma^{D}_{p}}+\gamma_{b}\,h^{-1}\langle
q_{h},q_{h}\rangle_{\Gamma^{D}_{p}}$ $\displaystyle\geq c\|\nabla
q_{h})\|^{2}-c\|\nabla
q_{h}\cdot\boldsymbol{n}\|_{\Gamma^{D}_{p}}\|q_{h}\|_{\Gamma^{D}_{p}}+\gamma_{b}\,h^{-1}\|q_{h}\|^{2}_{\Gamma^{D}_{p}}$
$\displaystyle\geq c\|\nabla q_{h}\|^{2}-c\delta\,h\|\nabla
q_{h}\cdot\boldsymbol{n}\|^{2}_{\Gamma^{D}_{p}}-c\delta^{-1}h^{-1}\|q_{h}\|^{2}_{\Gamma^{D}_{p}}+\gamma_{b}\,h^{-1}\|q_{h}\|^{2}_{\Gamma^{D}_{p}}$
$\displaystyle\geq c\|\nabla q_{h}\|^{2}-c\delta\|\nabla
q_{h}\|^{2}-c\delta^{-1}h^{-1}\|\boldsymbol{v}_{h}\|^{2}_{\Gamma^{D}_{p}}+\gamma_{b}\,h^{-1}\|q_{h}\|^{2}_{\Gamma^{D}_{p}}$
$\displaystyle=c\Big{(}\|\nabla
q_{h}\|^{2}+\widetilde{\gamma}_{b}\,h^{-1}\|q_{h}\|^{2}_{\Gamma^{D}_{p}}\Big{)}\,,$

with some constant $\widetilde{\gamma}_{b}>0$ for a sufficiently large choice
of the algorithmic parameter $\gamma_{b}$ in (3.3b), such that
$\widetilde{\gamma}_{b}:=\gamma_{a}-c\delta^{-1}>0$. Thus, there holds for
some constant $c>0$ that

$B_{\gamma}(q_{h},q_{h})\geq
c\|q_{h}\|_{h,\widetilde{\gamma}_{b}}\,,\quad\text{for all}\;\;q_{h}\in
Q_{h}\,.$ (A.5)

##  B Alternative formulation of the fully discrete problem

Here, we present an alternative formulation for the fully discrete system to
(1.1). The difference to Problem 3.1 comes through using
$\partial_{t}\boldsymbol{u}_{\tau,h}$ instead of $\boldsymbol{v}_{\tau,h}$ in
(3.5c) by means of (3.1a) and (3.5a), respectively. However, this modification
requires the inclusion of an additional boundary integral in the resulting
equation (B.1c). This non-obvious adaptation results from the proof of well-
posedness of the discrete problem. The following problem is thus considered.

######  Problem B.1 (Numerically integrated $I_{n}$-problem with
$\partial_{t}\boldsymbol{u}_{\tau,h}$)

For given
$\boldsymbol{u}_{h}^{n-1}:=\boldsymbol{u}_{\tau,h}(t_{n-1})\in\boldsymbol{V}_{h}$,
$\boldsymbol{v}_{h}^{n-1}:=\boldsymbol{v}_{\tau,h}(t_{n-1})\in\boldsymbol{V}_{h}$,
and $p_{h}^{n-1}:=p_{\tau,h}(t_{n-1})\in Q_{h}$ with
$\boldsymbol{u}_{\tau,h}(t_{0}):=\boldsymbol{u}_{0,h}$,
$\boldsymbol{v}_{\tau,h}(t_{0}):=\boldsymbol{u}_{1,h}$ and
$p_{\tau,h}(t_{0}):=p_{0,h}$, find
$(\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau,h},p_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$
such that

$\displaystyle\begin{aligned}
&Q_{n}\big{(}\langle\partial_{t}\boldsymbol{u}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle-\langle\boldsymbol{v}_{\tau,h},\boldsymbol{\phi}_{\tau,h}\rangle\big{)}+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle=\langle\boldsymbol{u}_{h}^{n-1},\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle\,,\\\\[4.30554pt]
\end{aligned}$ (B.1a) $\displaystyle\begin{aligned}
&Q_{n}\Big{(}\langle\rho\partial_{t}\boldsymbol{v}_{\tau,h},\boldsymbol{\chi}_{\tau,h}\rangle+A_{\gamma}(\boldsymbol{u}_{\tau,h},\boldsymbol{\chi}_{\tau,h})+C(\boldsymbol{\chi}_{\tau,h},p_{\tau,h})\Big{)}+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\chi}_{\tau,h}^{+}(t_{n-1})\rangle\\\\[4.30554pt]
&\qquad=Q_{n}\Big{(}F_{\gamma}(\boldsymbol{\chi}_{\tau,h})\Big{)}+\langle\rho\boldsymbol{v}_{h}^{n-1},\boldsymbol{\chi}_{\tau,h}^{+}(t_{n-1})\rangle\,,\end{aligned}$
(B.1b) $\displaystyle\begin{aligned} &Q_{n}\Big{(}\langle
c_{0}\partial_{t}p_{\tau,h},\psi_{\tau,h}\rangle-C(\partial_{t}\boldsymbol{u}_{\tau,h},\psi_{\tau,h})+B_{\gamma}(p_{\tau,h},\psi_{\tau,h})\Big{)}\\\\[4.30554pt]
&\qquad\quad+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1})+\alpha\nabla\cdot\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\psi_{\tau,h}^{+}(t_{n-1})\rangle-\alpha\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1})\cdot\boldsymbol{n},\psi_{\tau,h}^{+}(t_{n-1})\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\\\\[4.30554pt]
&\qquad=Q_{n}\Big{(}G_{\gamma}(\psi_{\tau,h})\Big{)}+\langle
c_{0}p_{h}^{n-1}+\alpha\nabla\cdot\boldsymbol{u}_{h}^{n-1},\psi_{\tau,h}^{+}(t_{n-1})\rangle-\alpha\langle\boldsymbol{u}_{D}(t_{n-1})\cdot\boldsymbol{n},\psi_{\tau,h}^{+}(t_{n-1})\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\end{aligned}$
(B.1c)

for all
$(\boldsymbol{\phi}_{\tau,h},\boldsymbol{\chi}_{\tau,h},\psi_{\tau,h})\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})\times\mathbb{P}_{k}(I_{n};Q_{h})$.

######  Lemma B.2 (Existence and uniqueness of solutions to Problem B.1)

Problem B.1 admits a unique solution.

#### Proof

The proof follows basically the lines of the proof of Lem. 3.2 and is kept
short. Only differences to Lem. 3.2 are depicted. For the differences
$(\boldsymbol{u}_{\tau,h},\boldsymbol{v}_{\tau,h},p_{\tau,h})$ of two triples
satisfying (B.1) and the test functions
$\boldsymbol{\phi}_{\tau,h}=A_{\gamma}\boldsymbol{u}_{\tau,h}$,
$\boldsymbol{\chi}_{\tau,h}=\boldsymbol{v}_{\tau,h}$ and
$\psi_{\textacutedbl\tau,h}=p_{\tau}$ there holds that

$\displaystyle\frac{1}{2}\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}\big{(}A_{\gamma}(\boldsymbol{u}_{\tau,h},\boldsymbol{u}_{\tau,h})+\langle\rho\boldsymbol{v}_{\tau,h},\boldsymbol{v}_{\tau,h}\rangle+\langle
c_{0}p_{\tau,h},p_{\tau,h}\rangle\big{)}\,\mathrm{d}t+Q_{n}\big{(}B_{\gamma}(p_{\tau,h},p_{\tau,h})\big{)}$
(B.2)
$\displaystyle\quad+Q_{n}\big{(}C_{\gamma}(\boldsymbol{v}_{\tau,h}-\partial_{t}\boldsymbol{u}_{\tau,h},p_{\tau,h})\big{)}+\alpha\langle\nabla\cdot\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle{\color[rgb]{0,0,0}-\alpha\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1})\cdot\boldsymbol{n},p_{\tau,h}^{+}(t_{n-1})\rangle_{\Gamma^{D}_{\boldsymbol{u}}}}$
$\displaystyle\quad+A_{\gamma}(\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),u_{\tau,h}^{+}(t_{n-1}))+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\chi_{\tau,h}^{+}(t_{n-1})\rangle+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle=0\,.$

Now, let $l\in\\{1,\ldots,k+1\\}$ be arbitrary but fixed and
$\boldsymbol{\phi}_{\tau,h}\in\mathbb{P}_{k}(I_{n};\boldsymbol{V}_{h})$ be
chosen as

$\boldsymbol{\phi}_{\tau,h}(t):=\xi_{n,l}(t)\boldsymbol{\phi}_{h}\quad\text{with}\quad\xi_{n,l}(t):=\left(\prod_{i=1\atop
i\neq l}^{k+1}\big{(}t-t_{n,i}^{\text{GR}}\big{)}\right)\left(\prod_{i=1\atop
i\neq
l}^{k+1}\big{(}t_{n,l}^{\text{GR}}-t_{n,i}^{\text{GR}}\big{)}\right)^{-1}\in\mathbb{P}_{k}(I_{n};\mathbb{R})\,,\;\;\boldsymbol{\phi}_{h}\in\boldsymbol{V}_{h}\,,$

and the Gauss–Radau quadrature nodes $t_{n,\mu}^{\text{G}R}$, for
$\mu=1,\ldots,k+1$; cf. (2.4). By the exactness of the Gauss–Radau quadrature
formula (2.4) for all polynomials in $\mathbb{P}_{2k}(I_{n};\mathbb{R})$ we
deduce from (B.1a) that

$\displaystyle 0$
$\displaystyle=\frac{\tau_{n}}{2}\sum_{\mu=1}^{k+1}\hat{\omega}_{\mu}^{\text{GR}}(\langle\partial_{t}\boldsymbol{u}_{\tau,h}(t_{n,\mu}^{\text{GR}}),\boldsymbol{\phi}_{\tau,h}(t_{n,\mu}^{\text{GR}})\rangle-\langle\boldsymbol{v}_{\tau,h}(t_{n,\mu}^{\text{GR}}),\boldsymbol{\phi}_{\tau,h}(t_{n,\mu}^{\text{GR}})\rangle)+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{\phi}_{\tau,h}^{+}(t_{n-1})\rangle$
$\displaystyle=\frac{\tau_{n}}{2}\hat{\omega}_{l}^{\text{GR}}\langle\partial_{t}\boldsymbol{u}_{\tau,h}(t_{n,l}^{\text{GR}})-\boldsymbol{v}_{\tau,h}(t_{n,l}^{\text{GR}}),\boldsymbol{\phi}_{h}\rangle+\langle\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\xi^{+}_{n,l}(t_{n-1})\boldsymbol{\phi}_{h}\rangle\,.$

Thus, we have that ($l=1,\ldots,k+1$)

$\boldsymbol{v}_{\tau,h}(t_{n,l}^{\text{GR}})-\partial_{t}\boldsymbol{u}_{\tau,h}(t_{n,l}^{\text{GR}})=c_{n,l}\,\boldsymbol{u}^{+}_{\tau,h}(t_{n-1})\quad\text{with}\;\;c_{n,l}=2\tau_{n}^{-1}\,\left(\hat{\omega}_{l}^{\text{GR}}\right)^{-1}\,\xi^{+}_{n,l}(t_{n-1})\,.$
(B.3)

Substituting (B.3) into the third term on the left-hand side of (B.2), we get
that

$\displaystyle
Q_{n}\big{(}C_{\gamma}(\boldsymbol{v}_{\tau,h}-\partial_{t}\boldsymbol{u}_{\tau,h},p_{\tau,h})\big{)}$
$\displaystyle=\sum_{\mu=1}^{k+1}C_{\gamma}(\boldsymbol{u}_{\tau_{h}}^{+}(t_{n-1}),\xi^{+}_{n,\mu}(t_{n-1})p_{\tau,h}(t_{n,\mu}^{\text{GR}}))=C_{\gamma}(\boldsymbol{u}_{\tau_{h}}^{+}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1}))$
$\displaystyle=-\alpha\langle\nabla\cdot\boldsymbol{u}_{\tau_{h}}^{+}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle+\alpha\langle\boldsymbol{u}_{\tau_{h}}^{+}(t_{n-1})\cdot\boldsymbol{n},p_{\tau,h}^{+}(t_{n-1})\rangle_{\Gamma^{D}_{\boldsymbol{u}}}\,.$
(B.4)

Combining (B.2) with (B.4) then implies that

$\displaystyle
A_{\gamma}(\boldsymbol{u}_{\tau,h}(t_{n}),\boldsymbol{u}_{\tau,h}(t_{n}))+\langle\rho\boldsymbol{v}_{\tau,h}(t_{n}),\boldsymbol{v}_{\tau,h}(t_{n})\rangle+\langle
c_{0}p_{\tau,h}(t_{n}),p_{\tau,h}(t_{n})\rangle+2Q_{n}\big{(}B_{\gamma}(p_{\tau,h},p_{\tau,h})\big{)}$
(B.5)
$\displaystyle\quad+A_{\gamma}(\boldsymbol{u}^{+}_{\tau,h}(t_{n-1}),\boldsymbol{u}_{\tau,h}^{+}(t_{n-1}))+\langle\rho\boldsymbol{v}^{+}_{\tau,h}(t_{n-1}),\chi_{\tau,h}^{+}(t_{n-1})\rangle+\langle
c_{0}p^{+}_{\tau,h}(t_{n-1}),p_{\tau,h}^{+}(t_{n-1})\rangle=0\,.$

From (B.5) along with the discrete coercivity properties (A.3) and (A.5) we
deduce that

$\boldsymbol{u}_{\tau,h}(t_{n})=\boldsymbol{u}_{\tau,h}^{+}(t_{n-1})=\boldsymbol{0}\,,\quad\boldsymbol{v}_{\tau,h}(t_{n})=\boldsymbol{v}_{\tau,h}^{+}(t_{n-1})=\boldsymbol{0}\quad\text{and}\quad
p_{\tau,h}(t_{n})=p_{\tau,h}^{+}(t_{n-1})=0$

as well as

$p_{\tau,h}\big{(}t_{n,\mu}^{\text{GR}}\big{)}=0\,,\quad\text{for}\;\;\mu=1,\ldots,k+1\,.$

The rest then follows as in the proof of Lem. 3.2. $\blacksquare$

##  C Additional numerical experiments

Here, we present some additional results of our numerical experiments. Table
C.1 shows the results of the convergence study introduced in Subsec. 5.1 for
the solution (5.1) of (1.1). For the computations of Table C.1 the Taylor-Hood
pair of finite element spaces $\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$, with
$r=4$, is used instead of the pair
$\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$ chosen for the results of
Table 5.2. For the pressure approximation marginally smaller errors are
observed. In Table C.2 the computed errors of the approximation of (5.1) for
larger values of the Lamé parameters $\lambda$ and $\mu$ in the elasticity
tensor $\boldsymbol{C}$ are summarized. No significant increase of the errors
is observed, indicating the independence of the error constant on the
magnitude of $\boldsymbol{C}$.

| $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{2}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 2.3958497242e-03 | – | 1.7185669625e-02 | – | 7.6634974254e-04 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.0529091363e-04 | 4.51 | 5.4558642696e-04 | 4.98 | 3.8993048591e-05 | 4.30
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 5.9749805593e-06 | 4.14 | 1.9143193940e-05 | 4.83 | 1.6140682324e-06 | 4.59
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 3.6858181764e-07 | 4.02 | 9.6894400439e-07 | 4.30 | 9.0805047399e-08 | 4.15
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.2976509245e-08 | 4.00 | 5.6982904446e-08 | 4.09 | 5.5524640600e-09 | 4.03
| $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{\infty}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 1.3014883491e-02 | – | 1.1587851233e-01 | – | 4.9980604461e-03 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 8.7944378641e-04 | 3.89 | 4.1645287352e-03 | 4.80 | 2.6289118267e-04 | 4.25
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 5.6086076018e-05 | 3.97 | 2.0121086363e-04 | 4.37 | 9.4136464267e-06 | 4.80
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 3.5171476997e-06 | 4.00 | 1.1873567381e-05 | 4.08 | 5.0122396011e-07 | 4.23
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.1959938532e-07 | 4.00 | 7.3034502924e-07 | 4.02 | 2.8669553969e-08 | 4.13
Table C.1: $L^{2}(L^{2})$ and $L^{\infty}(L^{2})$ errors and experimental orders of convergence (EOC) with temporal polynomial degree $k=3$ and spatial degree $r=4$ for local spaces $\mathbb{Q}_{r}^{2}/\mathbb{Q}_{r-1}$.  | $\tau$ | $h$ | $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{2}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{2}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 1.1835824122e-02 | – | 2.8040747896e-02 | – | 2.9610298753e-03 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.5714797742e-03 | 2.91 | 5.7609393975e-03 | 2.28 | 3.9787786632e-04 | 2.90
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 1.8935735637e-04 | 3.05 | 5.1606189537e-04 | 3.48 | 3.7075698634e-05 | 3.42
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 2.3548556912e-05 | 3.01 | 6.0883693658e-05 | 3.08 | 4.5049814701e-06 | 3.04
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.9374871960e-06 | 3.00 | 7.5088512252e-06 | 3.02 | 5.6225700500e-07 | 3.00
| $\tau$ | $h$ |  $\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{\tau,h})\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|\boldsymbol{v}-\boldsymbol{v}_{\tau,h}\|_{L^{\infty}(\boldsymbol{L}^{2})}$ | EOC | $\|p-p_{\tau,h}\|_{L^{\infty}(L^{2})}$ | EOC
---|---|---|---|---|---|---|---
$\tau_{0}/2^{0}$ | $h_{0}/2^{0}$ | 6.5064235342e-02 | – | 1.3580389516e-01 | – | 1.4838995481e-02 | –
$\tau_{0}/2^{1}$ | $h_{0}/2^{1}$ | 1.1842822460e-02 | 2.46 | 1.8023863764e-02 | 2.91 | 2.1492360068e-03 | 2.79
$\tau_{0}/2^{2}$ | $h_{0}/2^{2}$ | 1.5191266076e-03 | 2.96 | 3.4010640997e-03 | 2.41 | 2.3825320454e-04 | 3.17
$\tau_{0}/2^{3}$ | $h_{0}/2^{3}$ | 1.9340613477e-04 | 2.97 | 4.6154471881e-04 | 2.88 | 3.0597316621e-05 | 2.96
$\tau_{0}/2^{4}$ | $h_{0}/2^{4}$ | 2.4393096482e-05 | 2.99 | 5.8122744877e-05 | 2.99 | 3.8749301950e-06 | 2.98
Table C.2: $L^{2}(L^{2})$ and $L^{\infty}(L^{2})$ errors and experimental
orders of convergence (EOC) with temporal polynomial degree $k=2$ and spatial
degree $r=3$ for local spaces
$\mathbb{Q}_{r}^{2}/\mathbb{P}_{r-1}^{\text{disc}}$ for Young’s modulus
$E=10000$ and Poisson’s ratio $\nu=0.35$, corresponding to the Lamé parameters
$\lambda=8642$ and $\mu=3704$.