Space-Time Finite Element Methods for Parabolic Evolution Problems with Non-smooth Solutions

TL;DR

This paper introduces space-time finite element methods with local stabilization and adaptive refinement for solving non-autonomous parabolic problems with low-regularity solutions, achieving improved convergence and efficient solutions.

Contribution

It develops new stabilized finite element schemes on unstructured meshes, provides a priori estimates for low-regularity solutions, and implements adaptive refinement with efficient GMRES solvers.

Findings

New a priori estimates for low-regularity solutions.

Adaptive refinement improves convergence rates.

Efficient GMRES solver with algebraic multigrid preconditioning.

Abstract

We propose consistent locally stabilized, conforming finite element schemes on completely unstructured simplicial space-time meshes for the numerical solution of non-autonomous parabolic evolution problems under the assumption of maximal parabolic regularity. We present new a priori estimates for low-regularity solutions. In order to avoid reduced convergence rates appearing in the case of uniform mesh refinement, we also consider adaptive refinement procedures based on residual a posteriori error indicators. The huge system of space-time finite element equations is then solved by means of GMRES preconditioned by algebraic multigrid.

Equations

$${\partial }_{t}u-{\mathrm{div}}_x(\nu {\nabla }_{x}u)=f\text{mbox}inQ,u=0\text{mbox}on\Sigma ,u=u_0\text{mbox}on{\Sigma }_0$$

$$\int_{Q}\bigl{(}-u\,\partial_{t}v+\nu\,\nabla_{x}u\cdot\nabla_{x}v\bigr{)}\;\mathrm{d}(x,t)=\int_{Q}\!f\,v\;\mathrm{d}(x,t)+\int_{\Omega}\!u_{0}\,v|_{t=0}\;\mathrm{d}x$$

$$a_h(u_h,v_h)=l_h(v_h),\forall v_h\in V_{0h},$$

$$a_h(u_h,v_h)=$$

$$\displaystyle\qquad\qquad+\nu\nabla_{x}u_{h}\cdot\nabla_{x}v_{h}-\theta_{K}h_{K}\mathrm{div}_{x}(\nu\nabla_{x}u_{h})\partial_{t}v_{h}\big{]}\mathrm{d}(x,t),$$

$$l_h(v_h)=$$

$$\|v\|_{h}^{2}=\frac{1}{2}\|v\|_{L_{2}(\Sigma_{T})}^{2}+\sum_{K\in\mathcal{T}_{h}}\Bigl{[}\theta_{K}h_{K}\|\partial_{t}v\|_{L_{2}(K)}^{2}+\|\nabla_{x}v\|_{L_{2}^{\nu}(K)}^{2}\Bigr{]},$$

$$\Vert v{\Vert }_{h,\ast }^2$$

$$a_h(u-u_h,v_h)=0,\forall v_h\in V_{0h}.$$

$${\mathbf{K}}_h{\mathbf{u}}_h={\mathbf{f}}_h,$$

$$\|u-u_{h}\|_{h}\leq\biggl{(}1+\frac{\mu_{b}}{\mu_{c}}\biggr{)}\inf_{v_{h}\in V_{0h}}\|u-v_{h}\|_{h,*},$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{H^k(K)}\le C_{{\mathfrak{I}}_{SZ}}{h}_K^{l-k}|v{|}_{H^l(S_K)},k=0,1.$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{L_2({\Sigma }_T)}$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_h$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{h,\ast }$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{L_2({\Sigma }_T)}^2$$

$$\displaystyle\leq\sum_{\mathclap{\begin{subarray}{c}K\in\mathcal{T}_{h}\\ \partial K\cap\Sigma_{T}\neq\emptyset\end{subarray}}}\bigl{[}c_{Tr}^{2}h_{K}^{-1}(\|v-\mathfrak{I}_{S\!Z}v\|_{L_{2}(K)}^{2}+h_{K}^{2}\|\nabla(v-\mathfrak{I}_{S\!Z}v)\|_{L_{2}(K)}^{2})\bigr{]}$$

$$\displaystyle\leq c_{Tr}^{2}\sum_{\mathclap{\begin{subarray}{c}K\in\mathcal{T}_{h}\\ \partial K\cap\Sigma_{T}\neq\emptyset\end{subarray}}}\bigl{[}C_{\mathfrak{I}_{S\!Z}}^{2}\ h_{K}^{-1}h_{K}^{2l}|v|_{H^{l}(S_{K})}^{2}+C_{\mathfrak{I}_{S\!Z}}^{2}\ h_{K}\,h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}^{2}\bigr{]}$$

$$\displaystyle\leq\max_{K\in\mathcal{T}_{h}}\bigl{(}2\,c_{Tr}^{2}C_{\mathfrak{I}_{S\!Z}}^{2}\bigr{)}\sum_{\mathclap{\begin{subarray}{c}K\in\mathcal{T}_{h}\\ \partial K\cap\Sigma_{T}\neq\emptyset\end{subarray}}}\bigl{[}h_{K}^{2l-1}|v|_{H^{l}(S_{K})}^{2}\bigr{]},$$

$$\Vert v-$$

$$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\Bigl{[}\theta_{K}h_{K}\,\|\partial_{t}(v-\mathfrak{I}_{S\!Z}v)\|_{L_{2}(K)}^{2}+\|\nu^{1/2}\nabla_{x}(v-\mathfrak{I}_{S\!Z}v)\|_{L_{2}(K)}^{2}\Bigr{]}$$

$$+\frac{1}{2}\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{L_2({\Sigma }_T)}^2$$

$$\displaystyle\leq\sum_{K\in\mathcal{T}_{h}}\Bigl{[}\theta_{K}h_{K}C_{\mathfrak{I}_{S\!Z}}^{2}h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}^{2}+{\nu}_{K}C_{\mathfrak{I}_{S\!Z}}^{2}h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}\Bigr{]}$$

$$+c_1\sum _{\begin{array}{c}K\in {\mathcal{T}}_h\\ \partial K\cap {\Sigma }_T\ne \varnothing \end{array}}{h}_K^{2l-1}|v{|}_{H^l(S_K)}^2$$

$$\displaystyle\leq\sum_{K\in\mathcal{T}_{h}}\bigl{(}C_{\mathfrak{I}_{S\!Z}}^{2}(\theta_{K}h_{K}+{\nu}_{K})+c_{1}h_{K}\bigr{)}h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}^{2}\leq c_{2}\sum_{K\in\mathcal{T}_{h}}h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}^{2},$$

$$\sum_{K\in\mathcal{T}_{h}}\Bigl{[}(\theta_{K}h_{K})^{-1}\|v-\mathfrak{I}_{S\!Z}v\|_{L_{2}(K)}^{2}+\theta_{K}h_{K}\|\mathrm{div}_{x}(\nu\nabla_{x}(v-\mathfrak{I}_{S\!Z}v))\|_{L_{2}(K)}^{2}\Bigr{]}.$$

$$\sum _{K\in {\mathcal{T}}_h}({\theta }_K{h}_K{)}^{-1}\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{L_2(K)}^2$$

$$\displaystyle\leq\sum_{K\in\mathcal{T}_{h}}C_{\mathfrak{I}_{S\!Z}}^{2}\Bigl{(}\frac{h_{K}}{\theta_{K}}\Bigr{)}h_{K}^{2(l-1)}|v|_{H^{l}(S_{K})}^{2}.$$

$$\Vert {\mathrm{div}}_x(\nu {\nabla }_x(v-{\mathfrak{I}}_{SZ}v)){\Vert }_{L_2(K)}^2=\Vert {\mathrm{div}}_x(\nu {\nabla }_{x}v)-{\mathrm{div}}_x(\nu {\nabla }_x({\mathfrak{I}}_{SZ}v)){\Vert }_{L_2(K)}^2$$

$${\mathrm{div}}_x(\nu {\nabla }_x({\mathfrak{I}}_{SZ}v))={\nu }_K{\mathrm{div}}_x({\nabla }_x({\mathfrak{I}}_{SZ}v))=0$$

$$\sum _{K\in {\mathcal{T}}_h}{\theta }_K{h}_K\Vert {\mathrm{div}}_x(\nu {\nabla }_{x}v){\Vert }_{L_2}^2=\sum _{K\in {\mathcal{T}}_h}{h}_K^{2(l-1)}{\theta }_K{h}_K^{1-2(l-1)}\Vert {\mathrm{div}}_x(\nu {\nabla }_{x}v){\Vert }_{L_2(K)}^2,$$

$$\Vert v-{\mathfrak{I}}_{SZ}v{\Vert }_{h,\ast }^2\le$$

$$+\sum _{K\in {\mathcal{T}}_h}({\theta }_K{h}_K^{3-2l})h_K^{2(l-1)}\Vert {\mathrm{div}}_x(\nu {\nabla }_{x}v){\Vert }_{L_2(K)}^2$$

$$\le$$

$$\sum _{K\in {\mathcal{T}}_h}{\theta }_K{h}_K$$

$$\le \sum _{K\in {\mathcal{T}}_h}2{\theta }_K{h}_K\left(\Vert {\mathrm{div}}_x(\nu {\nabla }_{x}v){\Vert }_{L_2}^2+\Vert {\mathrm{div}}_x(\nu {\nabla }_x({\mathfrak{I}}_{SZ}v)){\Vert }_{L_2}^2\right)$$

$$\sum _{K\in {\mathcal{T}}_h}2{\theta }_K{h}_K$$

$$\le \sum _{K\in {\mathcal{T}}_h}2{\theta }_K{h}_K{h}_K^{-2}{c}_{I,3}^2\Vert \nu {\nabla }_x({\mathfrak{I}}_{SZ}^{1}v-{\mathfrak{I}}_{SZ}v){\Vert }_{L_2(K)}^2$$

$$\le \sum _{K\in {\mathcal{T}}_h}2{\theta }_K{h}_K^{-1}{c}_{I,3}^2{\nu }_K^2\Vert {\nabla }_x({\mathfrak{I}}_{SZ}^{1}v-{\mathfrak{I}}_{SZ}v){\Vert }_{L_2(K)}^2$$

Full text

Keywords: Parabolic initial-boundary-value problems; Space-time finite element methods; Unstructured meshes; Adaptivity

1 Introduction

Parabolic initial-boundary value problems of the form

[TABLE]

describe not only heat conduction and diffusion processes but also 2D eddy current problems in electromagnetics and many other evolution processes, where $Q=\Omega\times(0,T)$, $\Sigma=\partial\Omega\times(0,T)$, and $\Sigma_{0}=\Omega\times\{0\}$ denote the space-time cylinder, its lateral boundary, and the bottom face, respectively. The spatial computational domain $\Omega\subset\mathbb{R}^{d}$, $d=1,2,3$, is supposed to be bounded and Lipschitz. The final time is denoted by $T$. The right-hand side $f$ is a given source function from $L_{2}(Q)$. The given coefficient $\nu$ may depend on the spatial variable $x$ as well as the time variable $t$. In the latter case, the problem is called non-autonomous. We suppose at least that $\nu$ is uniformly positive and bounded almost everywhere. We here consider homogeneous Dirichlet boundary conditions for the sake of simplicity. In practice, we often meet mixed boundary conditions. Discontinuous coefficients, non-smooth boundaries, changing boundary conditions, non-smooth or incompatible initial conditions, and non-smooth right-hand sides can lead to non-smooth solutions.

In contrast to the conventional time-stepping methods in combination with some spatial discretization method, or the more advanced, but closely related discontinuous Galerkin (dG) methods based on time slices, we here consider space-time finite element discretizations treating time as just another variable and the term $\partial_{t}u$ in (1) as convection term in time. Following [8], we derive consistent, locally stabilized, conforming finite element schemes on completely unstructured simplicial space-time meshes under the assumption of maximal parabolic regularity; see, e.g., [6]. Unstructured space-time schemes have clear advantages with respect to adaptivity, parallelization, and the numerical treatment of moving interfaces or special domains. We refer the reader to the survey paper [10] that provides an excellent overview of completely unstructured space-time methods and simultaneous space-time adaptivity. In particular, we would like to mention the papers [9] that is based on an inf-sup-condition, [4] that uses mesh-grading in time, and [1] that also uses stabilization techniques. All three papers treat the autonomous case.

We here present new a priori discretization error estimates for low-regularity solutions. In order to avoid reduced convergence rates appearing in the case of uniform mesh refinement, we also consider adaptive refinement procedures in the numerical experiments presented in Section 5. The adaptive refinement procedures are based on residual a posteriori error indicators. The huge system of space-time finite element equations is then solved by means of Generalized Minimal Residual Method (GMRES) preconditioned by an algebraic multigrid cycle. In particular, in the 4D space-time case that is 3D in space, simultaneous space-time adaptivity and parallelization can considerably reduce the computational time. The space-time finite element solver was implemented in the framework of MFEM. The numerical results nicely confirm our theoretical findings. The parallel version of the code shows an excellent parallel performance.

2 Weak formulation and maximal parabolic regularity

The weak formulation of the model problem (1) reads as follows: find $u\in H^{1,0}_{0}(Q):=\{u\in L_{2}(Q):\nabla_{x}u\in[L_{2}(Q)]^{d},\,u=0\ \text{on}\ \Sigma\}$ such that (s.t.)

[TABLE]

for all $v\in\hat{H}^{1}_{0}(Q)=\{u\in H^{1}(Q):u=0\ \text{on}\ \Sigma\cup\Sigma_{T}\}$, where $\Sigma_{T}:=\Omega\times\{T\}$. The existence and uniqueness of weak solutions is well understood; see, e.g., [7]. It was already shown in [7] that $\partial_{t}u\in L^{2}(Q)$ and $\Delta u\in L^{2}(Q)$ provided that $\nu=1$, $f\in L^{2}(Q)$, and $u_{0}=0$. This case is called maximal parabolic regularity. Similar results can be obtained under more general assumptions imposed on the data; see, e.g., [6] for some very recent results on the non-autonomous case.

3 Locally stabilized space-time finite element

methods

In order to derive the space-time finite element scheme, we need an admissible, shape regular decomposition $\mathcal{T}_{h}=\{K\}$ of the space-time cylinder $Q=\bigcup_{K\in\mathcal{T}_{h}}\overline{K}$ into finite elements $K$. On $\mathcal{T}_{h}$, we define a $H^{1}$ conforming finite element space $V_{h}$ by means of polynomial simplicial finite elements of the degree $p$ in the usual way; see, e.g., [2]. Let us assume that the solution $u$ of (2) belongs to the space $V_{0}=H^{\mathcal{L},1}_{0,\underline{0}}(\mathcal{T}_{h}):=\{u\in L_{2}(Q):\partial_{t}u\in L_{2}(K),\ \mathcal{L}u:=\mathrm{div}_{x}(\nu\nabla_{x}u)\in L_{2}(K)\ \forall K\in\mathcal{T}_{h},\ \text{and}\ u|_{\Sigma\cup\Sigma_{0}}=0\}$, i.e., we only need some local version of maximal parabolic regularity, and, for simplicity, we assume homogeneous initial conditions, i.e., $u_{0}=0$. Multiplying the PDE (1) on $K$ by a local time-upwind test function $v_{h}+\theta_{K}h_{K}\partial_{t}v_{h}$, with $v_{h}\in V_{0h}=\{v_{h}\in V_{h}:v_{h}=0\,\mbox{on}\,\Sigma\cup\Sigma_{0}\}$, $h_{K}=\mathrm{diam}(K)$, and a parameter $\theta_{K}>0$ which we will specify later, integrating over $K$, integrating by parts, and summing up over all elements $K\in\mathcal{T}_{h}$, we arrive at the following consistent space-time finite element scheme: find $u_{h}\in V_{0h}$ s.t.

[TABLE]

with

[TABLE]

The bilinearform $a_{h}(\,.\,,\,.\,)$ is coercive on $V_{0h}\times V_{0h}$ wrt to the norm

[TABLE]

i.e., $a_{h}(v_{h},v_{h})\geq\mu_{c}\|v_{h}\|_{h}^{2}$, $\forall v_{h}\in V_{0h}$, and bounded on $V_{0h,*}\times V_{0h}$ wrt to the norm

[TABLE]

i.e., $a_{h}(u_{h},v_{h})\leq\mu_{b}\|u_{h}\|_{h,*}\|v_{h}\|_{h}$, $\forall u_{h}\in V_{0h,*},\forall v_{h}\in V_{0h}$, where $V_{0h,*}:=H^{\mathcal{L},1}_{0,\underline{0}}(Q)+V_{0h}$; see [8, Lemma 3.8] and [8, Remark 3.13], respectively. The coercivity constant $\mu_{c}$ is robust in $h_{K}$ provided that we choose $\theta_{K}=\mathcal{O}(h_{K})$; see Section 5 or [8, Lemma 3.8] for the explicit choice. From the above derivation of the scheme, we get consistency $a_{h}(u,v_{h})=l_{h}(v_{h}),\;\forall v_{h}\in V_{0h}$, provided that the solution $u$ belongs to $H^{\mathcal{L},1}_{0,\underline{0}}(Q)$ that is ensured in the case of maximal parabolic regularity. The space-time finite element scheme (3) and the consistency relation immediately yield Galerkin orthogonality

[TABLE]

We deduce that (3) is nothing but a huge linear system of algebraic equations. Indeed, let $V_{0h}=\mathrm{span}\{p^{(j)},j=1,\dots,N_{h}\}$, where $\{p^{(j)},j=1,\dots,N_{h}\}$ is the nodal finite element basis and $N_{h}$ is the total number of space-time degrees of freedom (dofs). Then we can express each function in $V_{0h}$ in terms of this basis, i.e., we can identify each finite element function $v_{h}\in V_{0h}$ with its coefficient vector $\mathbf{v}_{h}\in\mathbb{R}^{N_{h}}$. Moreover, each basis function $p^{(j)}$ is also a valid test function. Hence, we obtain $N_{h}$ equations from (3), which we rewrite as a system of linear algebraic equations, i.e.

[TABLE]

with the solution vector $\mathbf{u}_{h}=(u_{j})_{j=1,\ldots,N_{h}}$, the vector \mathbf{f}_{h}=\bigl{(}l_{h}(p^{(i)})\bigr{)}_{i=1,\ldots,N_{h}}, and system matrix \mathbf{K}_{h}=\bigl{(}a_{h}(p^{(j)},p^{(i)})\bigr{)}_{i,j=1,\ldots,N_{h}} that is non-symmetric, but positive definite.

4 A priori discretization error estimates

Galerkin orthogonality (8), together with coercivity and boundedness, enables us to prove a Cèa-like estimate, where we bound the discretization error in the $\|\,.\,\|_{h}$-norm by the best-approximation error in the $\|\,.\,\|_{h,*}$-norm.

Lemma 1**.**

Let the bilinearform (4) be coercive [8, Lemma 3.8] with constant $\mu_{c}$ and bounded [8, Lemma 3.11, Remark 3.13] with constant $\mu_{b}$, and let $u\in H^{\mathcal{L},1}_{0,\underline{0}}(\mathcal{T}_{h})$ be the solution of the space-time variational problem (2). Then there holds

[TABLE]

where $u_{h}\in V_{0h}$ is the solution to the space-time finite element scheme (3).

Proof.

Estimate (9) easily follows from triangle inequality and Galerkin-orthogonality; see [8, Lemma 3.15, Remark 3.16] for details. ∎

Next, we estimate the best approximation error by the interpolation error, where we have to choose a proper interpolation operator $\mathfrak{I}_{*}$. For smooth solutions, i.e., $u\in H^{l}(Q)$ with $l>(d+1)/2$, we obtained a localized a priori error estimate, see [8, Theorem 3.17], where we used the standard Lagrange interpolation operator $\mathcal{I}_{h}$; see e.g. [2]. In this paper, we are interested in non-smooth solutions, which means that we only require $u\in H^{l}(Q)$, with some real $l>1$. Hence, we cannot use the Lagrange interpolator. We can, however, use so-called quasi-interpolators, e.g. Clément [3] or Scott-Zhang [2]. For this kind of operators, we need a neighborhood $S_{K}$ of an element $K\in\mathcal{T}_{h}$ which is defined as $S_{K}:=\{K^{\prime}\in\mathcal{T}_{h}:\overline{K}\cap\overline{K}^{\prime}\neq\emptyset\}.$ Let $u\in H^{l}(Q)$, with some real $l>1$, then, for the Scott-Zhang quasi-interpolation operator $\mathfrak{I}_{S\!Z}:L_{2}(Q)\rightarrow V_{0h}$, we have the local estimate (see e.g. [2, (4.8.10)])

[TABLE]

For details on the construction of such a quasi-interpolator, we refer to [2] and the references therein. For simplicity, we now assume that the diffusion coefficient $\nu$ is piecewise constant, i.e., $\nu|_{K}=\nu_{K}$, for all $K\in\mathcal{T}_{h}$. Then we can show the following lemma.

Lemma 2**.**

Let $l>1$ and $v\in V_{0}\cap H^{l}(\mathcal{T}_{h})$. Then the following interpolation error estimates are valid:

[TABLE]

with $s=\min\{l,p+1\}$ and positive constants $c_{1},c_{2}$ and $c_{3}$, that do not depend on $v$ or $h_{K}$ provided that $\theta_{K}=\mathcal{O}(h_{K})$ for all $K\in\mathcal{T}_{h}$. Here, $p$ denotes the polynomial degree of the finite element shape functions on the reference element.

Proof.

For the first estimate, we use the scaled trace inequality and the quasi-interpolation estimate (10) with $k=0,1$

[TABLE]

which corresponds to (11) with c_{1}=\max_{K\in\mathcal{T}_{h}}\bigl{(}2\,c_{Tr}^{2}C_{\mathfrak{I}_{S\!Z}}^{2}\bigr{)}. To show the second estimate (12), we use definition (6) and that $\nu$ is piecewise constant, the quasi-interpolation error estimate (10) with $k=1$, and the above estimate (11), and obtain

[TABLE]

with c_{2}=\max_{K\in\mathcal{T}_{h}}\bigl{(}C_{\mathfrak{I}_{S\!Z}}^{2}(\theta_{K}h_{K}+{\nu}_{K})+c_{1}h_{K}\bigr{)}. For the third estimate, we deduce from (7) that we only have to estimate the additional sum

[TABLE]

We start with the first $L_{2}$-term. We apply the quasi-interpolation estimate (10) with $k=0$ and obtain

[TABLE]

Note that the term $(h_{K}/\theta_{K})$ is bounded for $\theta_{K}=\mathcal{O}(h_{K})$. For the $L_{2}$-norm of the spatial divergence, we can distinguish between two cases: linear basis functions ($p=1$) and higher order basis functions ($p\geq 2$). For linear basis functions, we split the divergence of the gradient, obtaining

[TABLE]

for each element $K\in\mathcal{T}_{h}$. Since $\mathfrak{I}_{S\!Z}v$ is a linear polynomial and $\nu$ is piecewise constant, we deduce

[TABLE]

for each element $K\in\mathcal{T}_{h}$. Hence, we get

[TABLE]

where the norm is bounded since $v\in V_{0}$. Moreover, for $\theta_{K}=\mathcal{O}(h_{K})$, the term $\theta_{K}h_{K}^{1-2(l-1)}$ is bounded for $1\leq l\leq 2$, and the case $l>2$ is already treated in [8]. Combining all above estimates, we obtain

[TABLE]

with c_{3}=c_{2}+\max_{K\in\mathcal{T}_{h}}\bigl{\{}C_{\mathfrak{I}_{S\!Z}}^{2}(h_{K}/\theta_{K}),\theta_{K}h_{K}^{3-2l}\bigr{\}}. For the general case of higher order basis functions, i.e., $p\geq 2$, the divergence of the gradient does not vanish. First we split the divergence of the gradient and also the norms, obtaining

[TABLE]

The first term in the sum we have already estimated above, where we now replace $l$ by $s=\min\{l,p+1\}$, which is in our case ($l\leq 2$) again just $l$. For the second term, we insert $\nu_{K}\mathrm{div}_{x}(\nabla_{x}(\mathfrak{I}_{S\!Z}^{1}v))$ into the norm, where $\mathfrak{I}_{S\!Z}^{1}v$ is the linear quasi-interpolation of $v$. We observe that $\mathfrak{I}_{S\!Z}^{1}v-\mathfrak{I}_{S\!Z}v$ is a finite element function, i.e., we can apply the inverse equality for the $H(\mathrm{div}_{x})$-norm [8, Lemma 3.5]. Since the diffusion coefficient is piecewise constant, we obtain

[TABLE]

Now we insert and subtract $\nabla_{x}v$ and use the triangle inequality, which yields

[TABLE]

For both terms we can apply (10) (quasi-interpolation) and obtain

[TABLE]

Combining all of the above,

[TABLE]

where c_{3,p}=c_{2}+\max_{K\in\mathcal{T}_{h}}\bigl{\{}{8\,\theta_{K}h_{K}^{-1}c_{I,3}^{2}{\nu}_{K}^{2}C_{\mathfrak{I}_{S\!Z}}^{2}},2(\theta_{K}h_{K}^{3-2l})\bigr{\}} . ∎

Now we are in the position to prove our main theorem.

Theorem 3**.**

Let $p$ be the polynomial degree used, and let $u\in H^{l}(Q)\cap V_{0}$, with $l>1$, be the exact solution, and $u_{h}\in V_{0h}$ be the approximate solution of the finite element scheme (3). Furthermore, let the assumptions of Lemma 1 (Céa-like estimate) and 2 (quasi-interpolation estimates) hold. Then the a priori discretization error estimate

[TABLE]

holds, with $s=\min\{l,p+1\}$ and a positive generic constant $C$.

Proof.

Choosing the quasi-interpolant $v_{h}=\mathfrak{I}_{S\!Z}u$ in (9), we can apply the quasi-interpolation estimate (13) to obtain

[TABLE]

We set $C=c_{3}(1+\mu_{b}/\mu_{c})$ to obtain (14), which closes the proof. ∎

5 Numerical results

We implemented the space-time finite element scheme (3) in C++, where we used the finite element library MFEM111http://mfem.org/ and the solver library hypre222https://www.llnl.gov/casc/hypre/. The linear system was solved by means of the GMRES method, preconditioned by the algebraic multigrid BoomerAMG. We stopped the iterative procedure if the initial residual was reduced by a factor of $10^{-8}$. Both libraries are already fully parallelized with the Message Passing Interface (MPI). Therefore, we performed all numerical tests on the distributed memory cluster RADON1333https://www.ricam.oeaw.ac.at/hpc/ in Linz. For each element $K\in\mathcal{T}_{h}$, we choose $\theta_{K}=h_{K}/(\tilde{c}^{2}\overline{\nu}_{K})$, where $\tilde{c}$ is computed by solving a local generalized eigenvalue problem which comes from an inverse inequality, see [8] for further details.

5.1 Example: Highly oscillatory solution

As first test example, we consider the unit (hyper-)cube $Q=(0,1)^{d+1}$, with $d=2,3$, as space-time cylinder, and $\nu\equiv 1$. The manufactured function

[TABLE]

serves as the exact solution, where we compute the right hand side $f$ accordingly. This solution is highly oscillatory. Hence, we do not expect optimal rates for uniform refinement in the pre-asymptotic range. However, using adaptive refinement, we may be able to recover the optimal rates. We used the residual based error indicator proposed by Steinbach and Yang in [10]. For each element $K\in\mathcal{T}_{h}$, we compute

[TABLE]

where $u_{h}$ is the solution of (3), $R_{h}(u_{h}):=f+\mathrm{div}_{x}(\nu\nabla_{x}u_{h})-\partial_{t}u_{h}$ in $K$, and $J_{h}(u_{h}):=[\nu\nabla_{x}u_{h}]_{e}$ on $e\subset\partial K$, with $[\,.\,]_{e}$ denoting the jump across one face $e\subset\partial K$. We mark each element where the condition $\eta_{K}\geq\sigma\max_{K\in\mathcal{T}_{h}}\eta_{K}$ is fulfilled, with $\sigma$ an a priori chosen threshold, e.g., $\sigma=0.5$. Note that $\sigma=0$ results in uniform refinement. In Figure 1, we observe indeed reduced convergence rates for all polynomial degrees and dimensions tested. However, using an adaptive procedure, we are able to recover the optimal rates. Moreover, we significantly reduce the number of dofs required to reach a certain error threshold. For instance, in the case $d=3$ and $p=2$, we need $276\,922\,881$ dofs to get an error in the $\|\,.\,\|_{h}$-norm of $\sim 10^{-1}$, whereas we only need $26\,359$ dofs with adaptive refinement. In terms of runtime, the uniform refinement needed $478.57s$ for assembling and solving, while the complete adaptive procedure took $156.5s$ only. The parallel performance is also shown in Figure 1, where we obtain a nice strong scaling up to 64 cores. Then the local problems are too small (only $\sim 10\,000$ dofs for 128 cores) and the communication overhead becomes too large.

5.2 Example: Moving peak

For the second example, we consider the unit-cube $Q=(0,1)^{3}$, i.e. $d=2$. As diffusion coefficient, the choice $\nu\equiv 1$ is made. We choose the function

[TABLE]

as exact solution and compute all data accordingly. This function is smooth, and almost zero everywhere, except in a small region around the line from the origin $(0,0,0)$ to $(1,1,1)$. This motivates again the use of an adaptive method. We use the residual based indicator $\eta_{K}$ introduced in Example 5.1. In Figure 2, we can observe that we indeed obtain optimal rates for both uniform and adaptive refinement. However, using the a posteriori error indicator, we reduce the number of dofs needed to reach a certain threshold by one order of magnitude. For instance, in the case $p=2$, we need $16\,974\,593$ dofs to obtain an error of $\sim 7\times 10^{-5}$ with uniform refinement. Using adaptive refinement, we need $1\,066\,777$ dofs only.

6 Conclusions

Following [8], we introduced a space-time finite element solver for non-autonomous parabolic evolution problems on completely unstructured simplicial meshes. We only assumed that we have so-called maximal parabolic regularity, i.e., the PDE is well posed in $L_{2}$. We note that this property is only required locally in order to derive a consistent space-time finite element scheme. We extended the a priori error estimate in the mesh-dependent energy norm to the case of non-smooth solutions, i.e. $u\in H^{1+\epsilon}(Q)$, with $0<\epsilon\leq 1$. This is necessary, since we cannot expect a smooth solution, especially when dealing with discontinuous diffusion coefficients. For simplicity, we considered only piecewise constant diffusion coefficients. The extension to piecewise smooth coefficients is straight-forward, but more technical. In comparison to the previous result for sufficiently smooth solutions, we no longer have a completely localized estimate. We also have to include the neighborhood of an element $K\in\mathcal{T}_{h}$. This may not be sufficient if we have to deal with solutions that have different regularity in different subdomains of the whole space-time cylinder. In applications, these subdomains correspond, for instance, to different materials. However, Duan et al. [5] have shown a quasi-interpolation estimate that fits into this setting.

We performed two numerical examples with known solutions. The first example had a highly oscillatory solution, and the second one was almost zero everywhere except along a line through the space-time cylinder. Using a high-performance cluster, we solved both problems on a sequence of uniformly refined meshes, where we also obtained good strong scaling rates. In order to reduce the computational cost, we also applied an adaptive procedure, using a residual based error indicator. Indeed, using a simultaneously in space and time adaptive procedure reduced the computational as well as the memory cost by a large factor. Moreover, we could observe that, especially for $d=3$, the AMG preconditioned GMRES method solves the problem quite efficiently.