Energy-corrected FEM and explicit time-stepping for parabolic problems

TL;DR

This paper introduces an energy-corrected finite element method for parabolic problems with corners, enabling stable explicit time-stepping on quasi-uniform meshes and demonstrating effectiveness in complex geometries.

Contribution

It develops a novel energy-corrected FEM approach that handles corner singularities without mesh refinement, allowing efficient explicit time integration in 2D and 3D.

Findings

Method achieves optimal convergence despite corner singularities.

Explicit time-stepping scheme is stable and efficient on quasi-uniform meshes.

Numerical results confirm theoretical predictions and method flexibility.

Abstract

The presence of corners in the computational domain, in general, reduces the regularity of solutions of parabolic problems and diminishes the convergence properties of the finite element approximation introducing a so-called "pollution effect". Standard remedies based on mesh refinement around the singular corner result in very restrictive stability requirements on the time-step size when explicit time integration is applied. In this article, we introduce and analyse the energy-corrected finite element method for parabolic problems, which works on quasi-uniform meshes, and, based on it, create fast explicit time discretisation. We illustrate these results with extensive numerical investigations not only confirming the theoretical results but also showing the flexibility of the method, which can be applied in the presence of multiple singular corners and a three-dimensional setting. We also propose a fast explicit time-stepping scheme based on a piecewise cubic energy-corrected discretisation in space completed with mass-lumping techniques and numerically verify its efficiency.

Figures

Tables

Table 1. Table 1 . Summary of convergence rates obtained using two different approximations of the heat equation on the L-shape domain

L	${\Vert u-u_h\Vert }_0$	rate	${\Vert u-u_h\Vert }_{0,\alpha }$	rate	${\Vert u-u_h^m\Vert }_0$	rate	${\Vert u-u_h^m\Vert }_{0,\alpha }$	rate	${\Vert u-{\tilde{u}}_h^m\Vert }_0$	rate
1	9.9471e-02		7.7723e-02		1.0172e-01		7.9832e-02		6.8808e-02
2	3.3940e-02	1.55	2.3263e-02	1.74	3.2843e-02	1.63	2.2848e-02	1.80	2.2432e-02	1.62
3	1.2351e-02	1.46	7.6646e-03	1.60	9.5573e-03	1.78	5.6370e-03	2.02	5.7524e-03	1.96
4	4.6633e-03	1.41	2.7192e-03	1.50	2.7736e-03	1.78	1.3492e-03	2.06	1.3562e-03	2.08
5	1.7942e-03	1.38	1.0130e-03	1.42	8.2665e-04	1.74	3.2694e-03	2.04	3.1236e-04	2.12
6	6.9729e-04	1.36	3.8825e-04	1.38	2.5229e-04	1.71	8.0793e-05	2.02	7.1894e-05	2.12

Equations

$$u_t-\Delta u$$

$$u$$

$$⟨u_t(t),v⟩+a(u(t),v)$$

$$\displaystyle H^{m}_{\beta}(\Omega)\coloneqq\big{\{}u\text{ - measurable }:\;r^{\beta+|\mu|-m}D^{\mu}u\in L^{2}(\Omega),\;\;0\leq|\mu|\leq m\big{\}},$$

$$\displaystyle\|u\|_{m,\beta}\coloneqq\bigg{(}\sum_{|\mu|\leq m}\big{\|}r^{\beta+|\mu|-m}D^{\mu}u\big{\|}^{2}_{L^{2}(\Omega)}\bigg{)}^{1/2}.$$

$$H_{\beta +l}^{m+l}(\Omega )\hookrightarrow H_{\beta }^m(\Omega ),\text{for any }l\in {\mathbb{Z}}_{+}.$$

$$\Vert v{\Vert }_{-\alpha }\le c_{\alpha }\Vert \nabla v{\Vert }_0,\text{for any }v\in H_0^1(\Omega ).$$

$$-\Delta w=f\text{in }\Omega ,w=0\text{on }\partial \Omega .$$

$$w=\sum _{0<n<(1+\alpha )\Theta /\pi }{k}_n{s}_n+W.$$

$$\displaystyle s_{n}(r,\phi)=\eta(r)r^{n\pi/\Theta}\sin\Big{(}\frac{n\pi}{\Theta}\phi\Big{)},$$

$$k_n=-\frac{1}{n\pi }{\int }_{\Omega }f{s}_{-n}+u\Delta s_{-n}.$$

$$u(t,\mathbf{x})=U(t,\mathbf{x})+k_1(t)s_1(\mathbf{x}).$$

$$\displaystyle\max_{0\leq t\leq T}\|\Delta u\|_{-\alpha}<\infty,\quad\text{and }\int_{0}^{T}\|\Delta u_{t}\|_{-\alpha}^{2}<\infty,\quad\text{and }\int_{0}^{T}\big{\|}u_{tt}\big{\|}_{\alpha}^{2}<\infty.$$

$$S_h:=\{v\in \mathcal{C}(\overline{\Omega }):v{|}_T\in P_1(T),\text{for all }T\in {\mathcal{T}}_h\},$$

$$V_h:=S_h\cap H_0^1(\Omega )\subset H_0^1(\Omega ).$$

$$|v_h{|}_{H^m(\Omega )}\le c_i{h}_{\min }^{l-m}|v_h{|}_{H^l(\Omega )},\text{for all }{v}_h\in S_h,$$

$$|v_h{|}_{H^m(\Omega )}\le c_i{h}^{l-m}|v_h{|}_{H^l(\Omega )},\text{for all }{v}_h\in S_h,$$

$$⟨u_{h,t}(t),v_h⟩+a(u_h(t),v_h)$$

$$u_h(0)$$

$$\displaystyle\big{\|}\nabla(u-u_{h})\big{\|}_{0}\lesssim h,\quad\|u-u_{h}\|_{0}\lesssim h^{2}.$$

$$a(w_h,v_h)=⟨f,v_h⟩\text{for all }{v}_h\in V_h.$$

$$\Vert w-w_h{\Vert }_{\alpha }\gtrsim \Vert \nabla (w-w_h){\Vert }_0^2\gtrsim h^{2\pi /\Theta }.$$

$$\displaystyle\big{\|}\nabla(u-u_{h})\big{\|}_{0}\lesssim h^{\pi/\Theta},\quad\|u-u_{h}\|_{0}\lesssim h^{2\pi/\Theta}.$$

$$a_h(w_h^{\text{m}},v_h)=⟨f,v_h⟩\text{for all }{v}_h\in V_h,$$

$$a_h(v_h,w_h)\le c_b\Vert \nabla v_h{\Vert }_{H_0^1(\Omega )}\Vert \nabla w_h{\Vert }_{H_0^1(\Omega )},\text{and }{a}_h(v_h,v_h)\ge c_c\Vert \nabla v_h{\Vert }_{H_0^1(\Omega )}^2.$$

$$c_h(w,v)=\gamma {\int }_{{\omega }_h}\nabla w\cdot \nabla v,$$

$$a_h(R_h^{\text{m}}w,v_h)=a(w,v_h)\text{for all }{v}_h\in V_h.$$

$$a(s_1-s_{1,h}^{\mathrm{m}},s_1-s_{1,h}^{\mathrm{m}})-c_h(s_{1,h}^{\mathrm{m}},s_{1,h}^{\mathrm{m}})=\mathcal{O}(h^2).$$

$$\Vert w-w_h^{\mathrm{m}}{\Vert }_{\alpha }\lesssim h^2\Vert f{\Vert }_{-\alpha },\Vert w-w_h^{\mathrm{m}}{\Vert }_{1,\alpha }\lesssim h\Vert f{\Vert }_{-\alpha }.$$

$$w_h^{\text{m}}=k_1{s}_{1,h}^{\text{m}}+R_h^{\text{m}}W,$$

$$k_1^h=-\frac{1}{\pi }{\int }_{\Omega }f{s}_{-1}+w_h^{\text{m}}\Delta s_{-1}.$$

Full text

1. Introduction

Numerical approximations of parabolic problems have been extensively studied in multiple settings and are of great interest due to the many technical applications, in which they appear [15, 31, 34]. Standard approximation methods involve finite differences [30, 41], but to allow for computations on more complicated domains, the theory of finite element methods has been developed, among many others, in [19, 22, 43]. For a more exhaustive discussion of the standard finite element approximations of parabolic problems and an extensive list of references, we refer the reader to [42].

The presence of corners in the computational domain negatively influences regularity properties of the solutions of elliptic [23] and parabolic problems [7, 23], due to the appearance of certain known singular functions. For corners with angles $\Theta>\pi$, in general, $H^{2}$ regularity in space cannot be guaranteed. This introduces the so-called pollution effect, diminishing the convergence order of the standard finite element schemes both for elliptic [9, 10] and parabolic problems [16]. Standard methods for improving the approximation properties of the finite element method in the elliptic setting include refinement [4, 5] and grading [1, 2, 38] of the mesh around the singular corner. These results can also be extended to parabolic problems on non-convex polygonal domains [16]. However, due to the very small mesh size in the vicinity of the re-entrant corner, the CFL condition [30], which guarantees the stability of an explicit time integration scheme, becomes very restrictive. This, in turn, means that explicit time-stepping schemes cannot be efficiently used, as they would require prohibitively small time-steps.

In this article, we follow a different approach, based on the so-called energy corrected finite element method. The idea was originally proposed for finite difference schemes in [35, 37, 44]. Recently, it has been extended to cover finite element methods for the Poisson equation in [20, 25, 36, 40], the Stokes equation in [26] and optimal control problems in [27]. The energy-correction method eliminates the pollution effect from the finite element approximation by a scaling of a fixed number of entries in the stiffness matrix. As opposed to adaptivity and grading, it can be successfully applied on quasi-uniform meshes, and hence, does not suffer from a too severe CFL condition. This permits the construction of fast explicit time-stepping schemes combined with the energy-corrected finite element in space.

This article is structured as follows: in Section 2 we discuss the regularity properties of the parabolic equations and briefly review the main results concerning the energy-correction method. In Section 3, we generalize the energy-corrected finite element to parabolic problems and provide a complete error analysis. We illustrate the analysis with numerical investigations in Section 4. In Section 5, we present several potential extensions of the scheme. We show that it can be applied to the advection-diffusion problem with a moderate advection term. Furthermore, we introduce higher-order energy-corrected discretisation combined with mass-lumping technique and a post-processing approach improving convergence properties of the scheme also in the vicinity of the singular corner. We complete the discretisation with an explicit Runge-Kutta time-stepping and show that the energy-corrected scheme exhibits superior performance compared with other commonly used discretisation methods. Finally, we present a potential application of the method, showing the flexibility of the energy-corrected finite element, which can be applied to problems with multiple re-entrant corners, and in a three-dimensional setting.

2. Parabolic problem

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded, non-convex, polygonal domain, i.e., a domain containing a re-entrant corner of size $\pi<\Theta<2\pi$. For the sake of simplicity we shall assume that this domain contains only one such corner. Note however that the analysis presented here also applies in a more general setting of domains with an arbitrary number of re-entrant corners.

Consider a standard heat equation defined on $\Omega$ in a time interval $[0,T]$ with $T>0$

[TABLE]

We define a corresponding weak solution as a continuous function $u:[0,T]\rightarrow H^{1}_{0}(\Omega)$ such that $u(0)=u_{0}\;a.e.$ and for all $v\in H^{1}_{0}(\Omega)\text{ and }a.e.\;0<t\leq T$

[TABLE]

where a(u(t),v)\coloneqq\big{\langle}\nabla u(t),\nabla v\big{\rangle}.

The unique solution to this problem exists and its regularity on smooth and convex domains has been studied extensively, see [21]. Having smooth initial conditions $u_{0}\in C^{\infty}(\overline{\Omega})$ and a smooth source term $f\in C^{\infty}(\overline{\Omega}\times[0,T])$ the solution of Equation (2.1) is also smooth, see [21, Chapter 7, Theorem 7]. However, this does not hold anymore, if polygonal domains are concerned.

2.1. Definitions and auxiliary results

We present our analysis in weighted Sobolev spaces as they constitute a convenient framework for the description of elliptic and parabolic problems defined on domains with re-entrant corners. Let $r(\mathbf{x})$ be the Euclidean distance of $\mathbf{x}\in\mathbb{R}^{2}$ from the re-entrant corner. For any $\beta\in\mathbb{R}$ we define a weighted Sobolev space as the following vector space

[TABLE]

where $\mu$ is a multiindex with non-negative entries. Moreover, we equip this space with a norm

[TABLE]

For convenience we denote $L^{2}_{\beta}(\Omega)\coloneqq H^{0}_{\beta}(\Omega)$. We also use the notational conventions $\|u\|_{0,\beta}=\|u\|_{\beta}$ and $\|u\|_{L^{2}(\Omega)}=\|u\|_{0}$. Similarly as standard Sobolev spaces, weighted spaces form a natural hierarchy, see [29]

[TABLE]

In particular for $m=0$ and $l=1$ we have

[TABLE]

for some $c_{\alpha}>0$ depending only on $\alpha$ and the domain $\Omega$.

2.2. Regularity results

In order to present regularity results for the parabolic problem on non-convex domains, we would like to first summarise regularity properties of a related elliptic problem. We choose this approach, as elliptic and parabolic problems exhibit the same type of singular behaviour in the vicinity of a re-entrant corner.

Consider

[TABLE]

Moreover, let $1-\alpha<\pi/\Theta$. It is well known, see, e.g., [28], that in the case of domains with re-entrant corners, the $H^{2}(\Omega)$ regularity of the solution $w$ usually cannot be obtained regardless of the regularity of the forcing term $f$. However, the solution $w$ can be split into regular and singular parts, namely

[TABLE]

Here, $W\in H^{2}_{-\alpha}(\Omega)$ denotes the smooth remainder, and $s_{n}$ are the singular functions defined as

[TABLE]

where $(r,\phi)$ are the polar coordinates in the vicinity of a re-entrant corner corresponding to the angle $\Theta$. Moreover, $\eta$ is a smooth cut-off function equal to $1$ in an arbitrary neighbourhood with a fixed distance from the re-entrant corner and equal to [math] far from it. Notice that $w\in H^{2}_{\alpha}(\Omega)$. Furthermore, a precise formula for the stress-intensity factors $k_{n}$ is known [23, Chapter 8]

[TABLE]

The following regularity result for the parabolic system (2.1) was proposed in [7].

Theorem 1**.**

Let $u(t)$ be a weak solution of the heat equation (2.4) on a polygonal domain $\Omega\subset\mathbb{R}^{2}$ with a re-entrant corner of size $\Theta$. Let also f\in C^{\sigma}\big{(}[0,T],L^{2}(\Omega)\big{)}, $\sigma>0$. Then there exists U\in C\big{(}[0,T];H^{2}(\Omega)\big{)} and k_{1}(t)\in C\big{(}[0,T)\big{)}\cap C^{1}\big{(}(0,T)\big{)} such that

[TABLE]

Moreover, $u(t,\cdot)\in H^{2}_{\alpha}(\Omega)$ for all $t>0$.

The regularity of parabolic equations on domains with conical points was further studied in [32, Theorem 3.2] with the results presented in the framework of weighted Sobolev spaces. The following theorem is an important consequence of the analysis presented there

Theorem 2**.**

Suppose that the assumptions of Theorem 1 are satisfied. Furthermore, let f\in L^{2}\big{(}0,T;H^{4}_{-\alpha}(\Omega)\big{)}, \frac{\mathrm{d}f}{\mathrm{d}t}\in L^{2}\big{(}0,T;H^{2}_{-\alpha}(\Omega)\big{)}, \frac{\mathrm{d}^{2}f}{\mathrm{d}t^{2}}\in L^{2}\big{(}0,T;L^{2}_{-\alpha}(\Omega)\big{)}, and $u_{0}\in H^{3}_{-\alpha}(\Omega)$. Assume also that the standard compatibility condition $f(0)+\Delta u_{0}\in H^{1}(\Omega)$ is satisfied. Then,

[TABLE]

We would like to stress out that the same singular functions arise in the solution of both elliptic and parabolic problems on non-convex polygonal domains. Note also that the results of Theorem 1 and Theorem 2 can be easily extended to multiple singular functions.

2.3. Finite element discretisation

Suppose that $\mathcal{T}_{H}$ is an admissible and shape-regular coarse initial triangulation of the domain $\overline{\Omega}$. We introduce the space of piecewise linear and globally continuous functions by

[TABLE]

and define the conforming finite element space by

[TABLE]

The following inverse inequality is a fundamental property of the discrete spaces and follows from the equivalence of norms in the finite dimentional spaces, see [12, Section 6.8] for more details.

Lemma 3** (Global inverse inequality).**

Let $0\leq l\leq m\leq k$. There exists a constant $c>0$ such that

[TABLE]

where $h_{\min}=\min_{T\in\mathcal{T}_{h}}h_{T}$. In particular, for uniformly refined meshes, we have

[TABLE]

where $c_{1}h\leq h_{\min}\leq c_{2}h$ for some $c_{1},c_{2}>0$.

From now on, we assume that $\mathcal{T}_{h}$ be a mesh obtained by uniform refinement of $\mathcal{T}_{H}$.

We define the semi-discretisation of the problem (2.4) in space as finding $u_{h},u_{h,t}[0,T]\rightarrow V_{h}$ such that

[TABLE]

The operator $\mathcal{P}_{h}:L^{2}(\Omega)\rightarrow S_{h}$ is linear and will be precisely specified later. For now, we assume that $\mathcal{P}_{h}$ is the $L^{2}$-projection. It is well known [11, 42] that then the finite element approximation in a sufficiently regular setting exhibits optimal second-order convergence in the $L^{2}$-norm and first-order convergence in the energy $H^{1}$-norm, namely for $a.e.\;t\in[0,T]$

[TABLE]

The situation is significantly different when re-entrant corners of maximum angle $\Theta>\pi$ are present in the computational domain.

Let us now consider the model elliptic boundary value problem (2.7). Remember, $\Omega$ is a polygonal domain with a re-entrant corner of angle $\pi<\Theta<2\pi$. In case of the standard piecewise linear finite element approximation of the model problem (2.7), we find $w_{h}\in V_{h}$ such that

[TABLE]

Due to the reduced regularity of the solution of (2.7), as summarised in Section 2.2, the convergence order of the finite element approximation (2.14) is also not optimal, when measured in the standard and weighted $L^{2}(\Omega)$-norms. This behaviour is known as the so-called pollution effect, see, e.g., [9, 10, 14, 17, 39].

Theorem 4** (Pollution effect).**

Let $w$ be the solution of (2.7) with $f\in L^{2}_{-\alpha}(\Omega)$ for some $1-\alpha<\pi/\Theta$. Further, assume that $k_{1}\neq 0$, then

[TABLE]

The proof of Theorem 4 can be found in [20]. Notice that the suboptimal approximation order is also obtained far from the re-entrant corner, so even in the case of elliptic equations, standard piecewise polynomial finite element approximation yields suboptimal convergence order.

This translates directly to parabolic problems, since in the presence of non-convex corners in the polygonal domain, the following convergence rates can be observed [16]

[TABLE]

These rates can be improved using suitable mesh-grading techniques, so that the optimal convergence in space is regained [42]. However, the corresponding CFL condition for explicit time-stepping schemes, meaning that the time-step needs to be scaled like a square of the size of the smallest element [30], gets very prohibitive and makes the use of explicit time-stepping schemes less attractive.

2.4. Energy-corrected finite element for elliptic equations

Here, we give a brief overview of the energy-correction techniques used for improving the convergence order in the finite element approximations of elliptic problems on polygonal domains. The idea was originally proposed for finite difference schemes in [35, 37, 44] and has been further developed recently in the finite element setting in [20, 36, 40] and is based on a local modification of the bilinear form governing the problem. It was later extended to piecewise polynomial approximation spaces in [24]. One of the advantages of this method, which we will exploit later, is the possibility of using quasi-uniform meshes.

The pollution effect in Theorem 4 is a result of an insufficient approximation of the energy $|a(u,u)-a(u_{h},u_{h})|$ by standard finite element techniques on uniform meshes. In order to remove this effect in the finite element approximation (2.14), we introduce a modification of the bilinear form $a(\cdot,\cdot)$, which shall mitigate the stiffness of the problem in the vicinity of the singularity. The modified finite element approximation of (2.7) reads then: find $w_{h}^{\text{m}}\in V_{h}$ such that

[TABLE]

where the bilinear form is defined as $a_{h}(w,v)\coloneqq a(w,v)-c_{h}(w,v)$. We assume that $a_{h}(\cdot,\cdot)$ is bilinear, continuous and elliptic, namely there exist $c_{b},c_{b}>0$ such that for all $v_{h},w_{h}\in V_{h}$

[TABLE]

Furthermore, we assume that $c_{h}(\cdot,\cdot)$ is symmetric. One possible choice of the modification is

[TABLE]

where $\omega_{h}$ is a one element patch around the re-entrant corner and $0<\gamma<1/2$. Due to the choice of the modification $c_{h}(\cdot,\cdot)$, we preserve the sparsity structure of the stiffness matrix, as only a small, fixed number of its entries needs to be suitably scaled. An additional assumption of symmetry of the nodal patch $\omega_{h}$ is necessary for angles $\Theta\geq 3\pi/2$.

Consequently, taking into account the energy-corrected formulation introduced in (2.16), we define a modified Ritz projection $R^{\text{m}}_{h}:H^{1}_{0}(\Omega)\rightarrow V_{h}$ by

[TABLE]

Let $s_{1,h}^{\text{m}}\in V_{h}$ denote the modified Ritz projection of the singular function $s_{1}$. The following theorem, providing sufficient conditions for the optimal convergence of the energy corrected method, was proposed in [20].

Theorem 5**.**

Let $1-\alpha<\pi/\Theta$ and $f\in L^{2}_{-\alpha}(\Omega)$. Let the modification $c_{h}(\cdot,\cdot)$ be defined as above and satisfy

[TABLE]

Then, for the energy-corrected finite element solution we obtain the following optimal error estimates

[TABLE]

Asymptotically, as $h\rightarrow 0$, there exists a unique, independent of $h$ optimal parameter $\gamma^{*}$ on the correction patch $\omega_{h}$. Note that the correction parameter depends on the number and shape of the elements $T$ of the correction patch $\omega_{h}$ and on the angle $\Theta$ of the re-entrant corner. Several effective procedures for finding it, based on nested Newton strategies, were proposed in [36].

Since weighted norms are not as commonly used as standard $L_{2}$-norms, we would like to obtain similar results to Theorem 5 for the latter. This can be done by exploiting some prior knowledge about the form of the singularity in the solution.

According to the decomposition (2.8) into singular and regular parts of the exact solution of (2.7), and by the linearity of the problem, we can represent the energy-corrected finite element approximation of $w$ as

[TABLE]

and the stress-intensity factor (2.10) can be efficiently approximated by

[TABLE]

Then, due to Theorem 5, we immediately obtain

[TABLE]

We define the post-processed approximation by

[TABLE]

Theorem 6**.**

The post-processed solution $\widetilde{w}_{h}^{\text{m}}$ defined in (2.23) converges with an optimal rate in standard norms, namely

[TABLE]

For a more detailed analysis of the post-processing we refer the reader to [24, Section 2].

3. Energy-corrected finite elements for parabolic problem

In this section, we investigate the energy-corrected finite element approximation of the parabolic problem (2.1). We begin by analysing the semi-discretisation, where the time variable is continuous. The reasoning is then further extended to a fully discrete case with the explicit Euler time-stepping.

3.1. Energy-corrected semi-discrete scheme

We define a modified semi-discrete finite element approximation as

[TABLE]

where $a_{h}(\cdot,\cdot)$ is an energy-corrected bilinear form introduced in Section 2.4. To impose the initial conditions we also use the modified Ritz projection (2.19) in place of the operator $\mathcal{P}_{h}$ in (2.12).

Lemma 7** (Stability of the semi-discrete scheme).**

Let $f\in L^{2}_{\alpha}(\Omega)$ for some $\alpha<1$. The semi-discrete solution of problem (3.1) satisfies for some $C^{*}>0$

[TABLE]

Proof.

An application of the ellipticity of the bilinear form $a_{h}(\cdot,\cdot)$, see (2.17), the Cauchy-Schwarz and Young inequalities, together with the choice $v_{h}=u_{h}^{\text{m}}$ in Equation (3.1), yields for any $\epsilon>0$

[TABLE]

Using the embedding (2.6) $H^{1}(\Omega)\hookrightarrow L^{2}_{-\alpha}(\Omega)$ for $1-\alpha>0$, we get

[TABLE]

where $c_{\alpha}$ is a constant coming from Inequality (2.6). Choosing $\epsilon\leq 2c_{c}/c^{2}_{\alpha}$ we arrive at

[TABLE]

where $C^{*}\geq c^{2}_{\alpha}/4c_{c}$. Integrating both sides over the time interval $[0,T]$ completes the proof of the lemma. ∎

Remark 8**.**

In the case of the standard choice of the modified bilinear form (2.18) in the energy-corrected scheme we have $c_{c}=1$ and the stability constant $C^{*}$ above can be reduced to $C^{*}=c^{2}_{\alpha}/4$. The embedding constant $c_{\alpha}$ from (2.6) is an equivalent of the Poincaré-Friedrichs constant in weighted spaces and depends only on the domain $\Omega$.

Theorem 9**.**

Suppose that functions $u_{0}$ and $f$ satisfy the regularity requirements stated in Theorem 1 and Theorem 2 and let $0<1-\alpha<\pi/\Theta$. The energy corrected semi-discretisation (3.1) of Problem (2.1) yields optimal convergence rate in the weighted $L^{2}$-norm, namely for some $c>0$ independent of $u$

[TABLE]

Proof.

We proceed in a standard manner by splitting the discretisation error into two independent parts

[TABLE]

where $R_{h}^{\text{m}}$ denotes the energy-corrected Ritz projection defined in Equation (2.19). Hence, due to Theorem 5

[TABLE]

and

[TABLE]

Using the definition of the modified Ritz projection (2.19), definition of the continuous solution (2.4) and the energy-corrected discretisation (3.1) we arrive at

[TABLE]

Finally, due to Lemma 7, we get

[TABLE]

Note that due to Theorem 2, the right-hand side of the inequality above is well-defined. Moreover, the discrete initial conditions were chosen in a way that $\eta(0)=0$.

Combining this with Equation (3.3) gives

[TABLE]

Finally, application of the results stated in (3.4) and (3.5) completes the proof.

∎

The right-hand side of Equation (3.2) is finite, see Theorem 2. The above theorem shows that the application of the energy-corrected finite element scheme to the parabolic equations results in the optimal accuracy of the scheme, when compared to the interpolation error.

3.2. Energy-corrected fully discrete scheme

Now, we move to the fully discrete setting, where also the temporal dimention is discretised. We consider only explicit Euler time-stepping, which later will serve as a foundation for building fast numerical schemes. The extension to a more general case of $\theta$-scheme in time is straightforward.

The fully discrete energy-corrected finite element approximation of the model problem (2.4) reads as follows: Find $U_{h}^{\text{m},n}\in V_{h}$ for $0\leq n\leq N$ such that

[TABLE]

The initial condition, as before, is imposed using the modified Ritz projection (2.19)

[TABLE]

In order to investigate the stability of the explicit Euler scheme, we need to introduce the so-called CFL (Courant–Friedrichs–Lewy) condition, which holds for all admissible triangulations.

Definition 10**.**

Consider the explicit Euler time-stepping introduced above. We define the CFL condition as

[TABLE]

where $h_{\min}=\min_{T\in\mathcal{T}_{h}}h$ and $c_{s}$ is a constant independent of the triangulation and the time-step $\Delta t$.

This condition was first introduced in [30] in the context of finite difference methods. The extension concerning the finite element methods for time-dependent problems can be found in [6, 13]. Note that in the case of uniform meshes $h_{\min}$ can be replaced with the mesh size $h$ in (3.7).

We begin the convergence analysis of the scheme by showing an auxiliary result bouding the finite difference in the formulation (4.1). Let $c_{i}>0$ denote the constant appearing in the inverse inequality from Lemma 3, when $l=0$ and $m=1$. Furthermore, let $c_{\alpha}>0$ and $c_{b}>0$ be respectively the constants in inequalities (2.6) and (2.17).

Lemma 11**.**

Suppose that f\in C\big{(}0,T;L^{2}_{\alpha}(\Omega)\big{)} for some $0\leq\alpha<1$. Then for all $0\leq n\leq N-1$

[TABLE]

Proof.

Let us set $v_{h}=\frac{U^{\text{m},n+1}_{h}-U^{\text{m},n}_{h}}{\Delta t}$ in (4.1). Then, applying the Cauchy-Schwarz inequality and using the boundedness of the bilinear form $a_{h}(\cdot,\cdot)$, we get

[TABLE]

Due to (2.6), we obtain

[TABLE]

Finally, application of the inverse inequality from Lemma 3 yields the desired result. ∎

Now, we can state the stability result in weighted Sobolev spaces, which will prove crucial for showing the error estimates for the fully discrete scheme. Similarly as in the case of the standard norms, the explicit time-stepping scheme is stable only under an additional assumption that the CFL condition (3.7) is satisfied. We provide the precise value of the stability constant.

Theorem 12** (Stability of the fully discrete scheme).**

Suppose that for some $0\leq\alpha<1$ we have f\in C\big{(}0,T;L^{2}_{\alpha}(\Omega)\big{)} and let $0<\epsilon<1/2$, $0<\delta<\frac{c_{c}}{c_{\alpha}^{2}}$. Suppose also that the CFL condition proposed in Definition 10 is satisfied with the constant $c_{s}=2\frac{c_{c}-c_{\alpha}^{2}\delta}{c_{i}^{2}c_{b}^{2}(1+\epsilon)}$. Then for some $c_{\epsilon,\delta}>0$ independent of $h$ and $\Delta t$, we have

[TABLE]

Proof.

We set $v_{h}=U^{\text{m},n}_{h}$ in (4.1). Notice that

[TABLE]

Hence

[TABLE]

Note that for any numbers $\epsilon,a,b>0$ we have

[TABLE]

Therefore, Lemma 11 gives us

[TABLE]

Furthermore, for any $\delta>0$ we get due to the Cauchy-Schwarz inequality and Equation (2.6)

[TABLE]

Using the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$ (2.17) and applying (3.8)–(3.10), we obtain

[TABLE]

The CFL condition (3.7) states that

[TABLE]

and therefore

[TABLE]

Setting

[TABLE]

and applying induction we finally obtain

[TABLE]

∎

Upon the right choice of the values $\epsilon,\delta$, we see that any $c_{s}<2\frac{c_{c}}{c_{i}^{2}c_{b}^{2}}$ is a feasible stability constant. Note however that when $\delta\rightarrow 0$ or $\epsilon\rightarrow 0$, then $c_{\epsilon,\delta}\rightarrow\infty$.

Finally, we are in a position to state the convergence result for the fully discrete scheme.

Theorem 13**.**

Suppose that functions $u_{0}$ and $f$ satisfy the regularity requirements stated in Theorem 9 and let $1-\lambda_{1}<\alpha<1$. Suppose also that the CFL condition stated in Definition 10 holds with the constant $c_{s}=2\frac{c_{c}-c_{\alpha}^{2}\delta}{c_{i}^{2}c_{b}^{2}(1+\epsilon)}$ for some $0<\epsilon<1/2$, $0<\delta<\frac{c_{c}}{c_{\alpha}^{2}}$. Then, the following error estimate for the energy-corrected discretisation $U^{\text{m},n}_{h}$, see (4.1), of Problem (2.1) holds for some $c>0$ independent of $u$

[TABLE]

Proof.

Similarly as in the proof of Theorem 9, we begin the proof by splitting the error into two components.

[TABLE]

where $R_{h}^{m}$ denotes the energy-corrected Ritz projection defined in Equation (2.19). Due to Theorem 5

[TABLE]

We focus now on estimating the remaining $\eta^{n}$ component of the error. Due to the definition of the energy-corrected Ritz projection and the problem formulation (2.4), we have

[TABLE]

Therefore, equation (4.1) yields for all $v_{h}\in V_{h}$

[TABLE]

Thus, we can write

[TABLE]

where

[TABLE]

Thanks to the stability estimate stated in Theorem 12 we obtain

[TABLE]

We now estimate $\psi_{1}^{n}$ and $\psi_{2}^{n}$ separately. Note that

[TABLE]

and hence

[TABLE]

Further, due to the linearity of the modified Ritz projection, we have

[TABLE]

and thus, see Theorem 5, we also get

[TABLE]

Further, applying Cauchy-Schwarz inequality, we obtain

[TABLE]

Since the initial conditions in the discretisation are imposed using the modified Ritz projection, see (4.1), we automatically have $\eta^{0}=0$. Using this and combining (3.14) with (3.15) and (3.16), we arrive at

[TABLE]

Finally, combining this result with (3.13), and applying to the error splitting (3.12) we get

[TABLE]

Hence, the proof is completed upon taking the square root of both sides of the inequality. The boundedness of the right-hand side is ensured by Theorem 2.

∎

As opposed to the mesh grading strategy, the energy-correction works on uniform meshes with less restrictive CFL stability condition (3.7). We shall exploit this fact further in the next section, when creating fast time-stepping schemes.

4. Numerical results

In this section, we propose and numerically investigate a fast solver for parabolic problems based on energy-corrected finite element (3.1). We show that, as opposed to the algorithms involving mesh grading and adaptivity, explicit time-stepping schemes are a feasible choice in the proposed setting.

We discretize Eq. (3.1) using the Explicit Euler time-stepping scheme. Let us divide the time interval $[0,T]$ into $N\in\mathbb{Z}_{+}$ time steps of equal lengths $\Delta t$, so $t_{n}=n\Delta t$. The fully discrete approximation of the model problem (2.4) reads as follows: Find $u_{h}^{m,n}\in V_{h}$ for $0\leq n\leq N$ such that

[TABLE]

Let \big{\{}\varphi_{i}\big{\}}_{i=1}^{K} be a standard nodal finite element basis. The discrete systems (4.1) can be rewritten in a matrix-vextor formulation as

[TABLE]

where \mathbf{M}=\Big{[}\big{(}\varphi_{i},\varphi_{j}\big{)}\Big{]}_{i,j=1}^{K} and \mathbf{S}^{m}=\Big{[}a_{h}(\varphi_{i},\varphi_{j}\big{)}\Big{]}_{i,j=1}^{K} denote the standard mass and energy-corrected stiffness matrices respectively, and \mathbf{F}^{n}=\Big{[}\big{(}f(t_{n},\cdot),\varphi_{i}\big{)}\Big{]}_{i=1}^{K}. The choice of a nodal, vertex-based quadrature rule for assembling the mass matrix leads to a lumped diagonal matrix $\widetilde{\mathbf{M}}$, which can be used in place of $\mathbf{M}$, see [8] for more details. This results in a fast time-stepping scheme, where at each time step, multiplication by a diagonal matrix $\widetilde{\mathbf{M}}$ and a sparse matrix $\mathbf{S}^{m}$ needs to be performed.

Stability of the Explicit Euler scheme is guaranteed by the CFL condition [42], meaning that the size of the time step needs to scale like the square of the mesh size, i.e. $\Delta t\sim h^{2}$. This is very prohibitive when mesh grading or adaptivity is concerned. However, this is not an issue in the case of the energy-corrected FEM, which works on uniform meshes. Then, balancing the error of order $\mathcal{O}(\Delta t)$ coming from the time-stepping discretization with $\mathcal{O}(h^{2})$ order of error measured in the weighted $L^{2}$-norm, see Thm. 9, exactly the same relationship needs to be kept.

In order to improve the convergence of the scheme at a fixed point in time $T$, we complete the algorithm with a post-processing strategy, following the post-processing strategy as in (2.23). As stated in Eq. (2.10), the stress-intensity factor, see Thm. 1, can be computed by

[TABLE]

We define its discrete approximation using (2.21) as

[TABLE]

This leads to the post-processed solution of the form

[TABLE]

Note that the additional cost of performing the post-processing is equal to the cost of solving one additional elliptic equation and evaluating one integral.

In Table 1, we summarise the errors and the convergence rates of the proposed scheme. We choose the L-shape domain \Omega=(-1,1)^{2}\setminus\big{(}[0,1]\times[-1,0]\big{)} with the largest interior angle of size $\Theta=3\pi/2$. We also choose a known exact solution $u=\sin(t)s_{1}+\sin(2t)s_{2}-\sin(3t)s_{3}$, being a linear combination of singular functions (2.9) with smooth time-dependent coefficients.

The parameter $\gamma$ in the modification (2.18) is computed using a version of the Newton algorithm described in [36] and in the numerical experiments we choose the weight $\alpha=1-\pi/\Theta$. This choice of the weight induces a slightly stronger norm than assumed in Sec. 3 but the optimal convergence order of the energy-corrected scheme can be observed regardless of this.

We consider uniform refinement of the initial mesh and together with refining the mesh, we also divide the time-step $\Delta t$ by $4$, initially set to be equal to $0.1$. For the purpose of comparison, in the first two columns of Table 1, we summarise the results obtained using the uncorrected scheme. The suboptimal convergence rates in the sense of the best approximation error, when measured both in standard and weighted $L^{2}(\Omega)$-norms, are in line with the results in (2.15). For the energy-corrected scheme (4.1), we see that no pollution in the $L^{2}(\Omega)$-norm appears. Moreover, second-order convergence in the weighted norm means that the error is relatively large only in the vicinity of the re-entrant corner, so the pollution effect from Thm. 4 has been removed. Finally, the post-processing approach yields second-order convergence in the standard $L^{2}(\Omega)$-norm. Numerical tests confirm the theoretical results of Theorem 9.

5. Extensions

In this section, we present extensions of the methods introduced above. We show that the energy-corrected finite element can be applied to domains with multiple re-entrant corners, also in the presence of a moderate advection in the problem (2.1). Furthermore, we show a possible extension to the higher-order piecewise polynomial finite elements and propose a fast explicit time-stepping scheme based on cubic elements combined with mass-lumping techniques. Finally, in order to show the flexibility of the energy-correction method, we present numerical experiments involving multiple re-entrant corners in three dimensions.

5.1. Advection-diffusion equation

In [40], pointwise error estimates for the energy-corrected finite element method for the elliptic problems on polygonal domains were studied. It was shown that the energy-corrected discretisation (2.16) of (2.7) yields, up to a logarithmic factor, optimal convergence in the sense of the optimal approximation property. In the following, we show that the improved pointwise convergence of the energy-corrected scheme can be also expected in the case of parabolic problems.

We consider the following advection-diffusion problem

[TABLE]

In the numerical example, we consider $T=1$, $u_{0}=0$ and $b=(1,1)$, so the problem is equipped with moderate advection. The computational domain $\Omega$, together with its triangulation, is presented in Figure 1 and consists of a rectangle $(0,4)\times(0,3)$ with a right, isosceles triangle cut out. There are three re-entrant corners in the domain $\Omega$, two of sizes $7\pi/4$ and one of size $3\pi/2$. We use a computational grid with one-element patches around the singular corners consisting of the identical isosceles triangles. For the right-hand side we choose f=\sin(\pi t)\Big{(}(x-2)^{2}+(y-3/2)^{2}\Big{)}^{-1}, which has a singularity in the middle of the cut-out triangle. For the time discretisation, we choose Explicit Euler time-stepping as described in Section 4 with an initial step-size $\Delta t=0.02$, which is small enough to guarantee the stability of the scheme.

We investigate the behaviour of a quantity of interest $\text{QoI}=\|u_{h}(T)\|_{L^{\infty}(\Omega)}$ for standard finite element method and the energy-corrected finite element method on $5$ consecutive refinement levels. We use the same modification of the bilinear form as before, namely (2.18). The results of the simulations are summarised in the plot on the right-hand side of Figure 1. For completeness, we also include the extrapolated approximation $\|u(T)\|_{\infty}^{ex}$ of the real value in the plot. The estimated order of convergence of $|\|u(T)\|_{\infty}^{ex}-\|u_{h}(T)\|_{\infty}|$ is equal to $1.85$ and $1.55$ in the case of the energy-corrected scheme and standard finite element respectively. The energy-corrected finite element can be successfully applied also in the cases of several different re-entrant corners in the domain and the presence of a moderate advection in the problem.

5.2. Piecewise polynomial FEM

In [24] ideas presented in Section 2.4 were extended to a more general setting of piecewise polynomial finite elements. Optimal convergence of the approximation (2.16) of the model Poisson problem (2.7), when using $k$-th order polynomial basis functions, is then obtained upon the choice of $f\in H^{k-1}_{-\alpha}(\Omega)$, and one of the modifying bilinear functions

[TABLE]

Here, $\omega^{i}_{h}$ denotes the $i$-th layer of elements, counting from the considered corner and $\hat{r}$ is a distance from the corner measured on the reference triangle. We also assume that the patch of elements around the corner consists of identical isosceles triangles. Asymptotically, unique optimal sequences of parameters \underline{\gamma}^{R*}=\big{(}\gamma^{R*}_{i}\big{)}_{i=1}^{K} and \underline{\gamma}^{F*}=\big{(}\gamma^{F*}_{i}\big{)}_{i=1}^{K} on a given correction patch exist.

We now focus our attention on cubic basis functions, which means that for $f\in H^{2}_{-\alpha}(\Omega)$ we have

[TABLE]

We apply the method to the parabolic problem (2.4) in a straightforward manner, suitably modifying the formulation (3.1).

We are aiming to construct a fast, explicit solver and this means that, in order to guarantee the stability of the method, the CFL condition $\Delta t\sim h^{2}$ needs to be satisfied. The use of cubic finite element basis yields fourth-order convergence in weighted $L^{2}(\Omega)$-norm, and we would like to balance it with a second-order time-stepping scheme. Similarly to the piecewise linear case discussed in Section 4, optimal balancing of the errors stemming from the time and space discretisations means that the CFL condition is automatically satisfied and makes the use of explicit time-stepping scheme feasible. We use the second-order Runge-Kutta scheme, also known as the Heun’s method. The fully discrete scheme can be written as

[TABLE]

Note that the application of the mass-lumping strategy is not as straightforward as in the piecewise linear case, since in general it results in a singular mass matrix for the higher-order finite element. To overcome this, we follow the method proposed in [18] in the context of the wave equation. It is based on the enrichment of the cubic finite element space with three fourth-order polynomial bubble functions, which are uniformly equal to [math] at the elements’ edges. Moreover, the mass matrix is assembled using a positive quadrature rule, which is exact for seventh-order polynomials, with quadrature points located in the nodal points of the enriched space. Such a construction yields a diagonal matrix $\widetilde{\mathbf{M}}$ used in place of $\mathbf{M}$ in Scheme (5.5).

In Figure 2, we compare the accuracy and normalised computational times of the energy-corrected schemes with several commonly used methods for discretising the heat equation. We use a known exact solution $u=\sin(t)s_{1}+\sin(2t)s_{2}-\sin(3t)s_{3}$ on the L-shaped domain \Omega=(-1,1)^{2}\setminus\big{(}[0,1]\times[-1,0]\big{)} with the largest interior angle of size $\Theta=3\pi/2$ and compare the $L^{2}(\Omega)$ and $L^{2}(\Omega^{\prime})$ errors of the schemes at the last time step $T=1$, where $\Omega^{\prime}=\Omega\cap\{|x|>0.25\}$.

On the left-hand side of Figure 2, standard $L^{2}(\Omega)$ errors and normalised computational times are shown. The standard piecewise linear finite element combined with an Explicit Euler time-stepping and mass-lumping provides the worst results among the ones studied since its performance is limited by (2.15). Initially $\Delta t=0.1$ is chosen and with each space refinement the time step is divided by $4$.

Application of mesh grading techniques improves the performance of the solver [16]. However, the use of explicit time-stepping schemes is infeasible due to the very restrictive CFL condition. In order to recover the optimal convergence order in the $L^{2}(\Omega)$ norm, it is necessary to grade the mesh towards the singular corner, introducing elements of the size $h^{1/\mu}$, where $\mu<\pi/\Theta$. This in particular means that in the case of the L-shape domain with $\Theta=3\pi/2$, time-steps $\Delta t$ smaller than $\mathcal{O}(h^{3})$ need to be used. Therefore, we use an unconditionally stable Crank-Nicolson scheme in time allowing for a coarser time discretisation. In order to keep the right balance between space and time discretisation errors, with each mesh refinement we divide the time-step by $2$, beginning with $\Delta t=0.1$.

Piecewise linear energy-corrected finite element scheme with Explicit Euler time-stepping (4.2) yields significantly better results than the standard piecewise linear discretisations. It also gives comparable results with the mesh grading scheme completed with Crank-Nicolson time-stepping. An application of the post-processing additionally improves the accuracy of the method resulting in a better error-to-time ratio than the mesh grading method equipped with Crank-Nicolson time-stepping. We use in-built MATLAB linear system solvers. Note that the application of fast iterative solvers, such as multigrid methods, could additionally improve the performance of the relevant implicit methods.

Application of the cubic energy-corrected finite element scheme with second-order Runge-Kutta scheme in time gives similar results to the piecewise linear energy-corrected scheme. However, additional application of the post-processing yields the best results in terms of the balance between the computational time and the accuracy of the scheme out of all tested methods. This can be attributed to the use of the scheme eliminating the pollution effect in the solution, completion with the post-processing strategy yielding optimal convergence in the standard norms, and the use of mass-lumping strategy. Note that the use of $c_{h}^{F}(\cdot,\cdot)$ modification gives quantitatively better results than $c_{h}^{R}(\cdot,\cdot)$. This phenomenon was previously observed in [24, Section 6.3] and can be attributed to the smaller modification subregion in the computational domain.

As shown in Section 3, energy-correction method gives optimal convergence rates in terms of the best-approximation property, however, in weighted norms. This, in particular, means that the method converges optimally when measured far from the re-entrant corner. Therefore, no additional post-processing needs to be applied, when one is interested in the solution far from the singular corner.

On the right-hand side of Figure 2, a comparison of $L^{2}(\Omega^{\prime})$ errors and normalised computational times of the previously described methods are shown. Again, due to the pollution effect, the standard finite element discretisation results in the worst error-to-time ratio. It can be improved by the application of the mesh grading together with Crank-Nicolson time-stepping, which yields only slightly worse results than the piecewise-linear energy-corrected scheme. The cubic energy-corrected finite element, together with Heun’s time-stepping and mass-lumping strategy, results in by far the best method when the $L^{2}$-error far from the re-entrant corner is concerned. Small variations in the convergence rates in the cubic finite element scheme, when using $c^{R}_{h}(\cdot,\cdot)$ modification, appear because of insufficient initial resolution of the mesh. Again, modification $c^{F}_{h}(\cdot,\cdot)$ yields a better performance than $c^{R}_{h}(\cdot,\cdot)$.

5.3. Application

In this section, we apply the piecewise linear energy-corrected finite element method to a real 3D geometry of a graphite moderator brick of a nuclear power plant. Such a moderation type is commonly used in Advanced Gas-cooled Reactors (AGR) [33]. Efficient simulations of heat distribution in moderator bricks play an important role in the analysis of the material properties of the whole nuclear core, and accurate computations of temperature distribution can help determine the lifetime of nuclear materials, which often suffer from large temperature gradients and fast neutron fluxes [3].

The 3D geometry $\Omega$ of the graphite brick is presented in Figure 3(a). The shape of the brick has a tensorial structure, with identical 2D horizontal cross-sections containing $16$ equally-sized re-entrant corners with angle $\Theta=3\pi/2$ at the external boundary. Heat transfer in a graphite moderator brick, in its simplest form, can be described by

[TABLE]

Here, $\partial\Omega_{1}$ is the interior, cylindrical boundary of the domain and $g$ is the heating produced due to the nuclear reaction occurring in the fuel assembly. The system is thermally isolated from below on the $\partial\Omega_{2}$ part of the boundary, which is reflected by the uniform Neumann boundary conditions. Finally, the remaining part of the domain’s boundary - $\partial\Omega_{3}$, is subject to a circulating coolant of a constant temperature. Note that the solution is rescaled, so that the temperature there is uniformly distributed.

Exploiting the tensorial structure of the domain $\Omega$, we divide it into prismatic elements of equal length $h_{z}$ in the vertical dimension. Moreover, each cross-section is triangulated as shown in Figure 3(b). It is worth noting that around each of the $16$ re-entrant corners in the cross-sections, we use identical one-element patches consisting of congruent isosceles triangles. This, together with the tensorial structure of the mesh, allows us to reuse the parameter $\gamma$ in the energy-corrected scheme (2.18) once computed in the two-dimensional setting. Following (4.1), we complete the finite element discretisation in space with the Explicit Euler time-stepping scheme.

In the simulations, we choose a homogenuous initial temperature distribution $u_{0}=0$ and the heating on the internal boundary to be given by g=10\cos\big{(}\pi z/(2L)\big{)}\big{(}1+\sin(\phi+4\pi t)\big{)}\sin(\pi t), where $L=3$ is the height of the domain $\Omega$ and $(r,\phi,z)$, $z\in[0,L]$ are cylindrical coordinates associated with the domain. The simulation is performed using $N=5000$ time steps with the final time at $T=1$ and at the end of the simulation the solution is rescaled to a physically meaningful value range by $U=20u+400$. The solution at two intermediate time-steps is presented in the vertical cross-sections of the domain taken in the middle at the height $z=1.5$ at times $t=0.25$ and $t=0.5$.

To investigate the convergence of the scheme, we measure two different quantities of interest, namely the average temperature in the body at the final time-step and the average temperature in the whole space-time cylinder $\Omega\times[0,T]$

[TABLE]

In order to investigate the convergence properties of the quantities of interest, we perform the computations on four different refinement levels. Initially, we use the mesh described above, which is then uniformly refined and the time-step size is divided by $4$. The results of the computations are presented in Figure 5. We compare the accuracy of the standard finite element with the energy-corrected finite element method. Together with the quantities of interest $\text{QoI}_{1}$ and $\text{QoI}_{2}$, we include also the extrapolated values $\text{QoI}^{ex}_{1}$ and $\text{QoI}^{ex}_{2}$ in the plots. As presented in Figure 5, the application of energy-correction improves the approximation properties of the quantities of interest. Moreover, in the case of the standard finite element approximation the respective estimated orders of convergence of $|\text{QoI}^{ex}_{1}-\text{QoI}_{1}|$ and $|\text{QoI}^{ex}_{2}-\text{QoI}_{2}|$ are $1.37$ and $1.23$. The application of the energy-correction improves these orders and yields estimated values of $2.1$ and $2.12$, respectively.

6. Concluding remarks

In this article, we proposed an energy-corrected finite element discretisation for parabolic problems on non-convex polygonal domains. We showed rigorously that the pollution effect inherited from elliptic problems can also be eliminated in the case of time-dependent problems, resulting in the optimal order of convergence in the sense of the best approximation property. The use of uniform meshes in the energy-correction method leads to a less restrictive CFL condition compared to standard methods based on mesh grading and adaptivity. Further, this allows for the application of explicit time-stepping schemes and creation of fast numerical schemes based on them.

Moreover, we proposed a post-processing approach and showed how the higher-order energy-corrected scheme can be combined with explicit Runge-Kutta type discretisation in time. This, together with the mass-lumping techniques, results in efficient solvers of parabolic problems for both piecewise linear and piecewise polynomial finite elements. Finally, we showed that the proposed algorithm can be successfully applied to 3D geometries with multiple re-entrant corners.

Acknowledgements

We gratefully acknowledge the support of the German Research Fundation (DFG) through the grant WO 671/11-1 and, together with the Austrian Science Fund, through the IGDK1754 Training Group.