A-posteriori error estimation and adaptivity for nonlinear parabolic equations using IMEX-Galerkin discretization of primal and dual equations

TL;DR

This paper develops a duality-based a posteriori error estimation method for nonlinear parabolic PDEs discretized with IMEX-Galerkin schemes, enabling adaptive mesh refinement and time-stepping.

Contribution

It introduces a novel error decomposition technique that separates spatial and temporal errors for nonlinear parabolic PDEs using IMEX discretization.

Findings

Effective error estimates for heat and Allen-Cahn equations

Adaptive algorithms improve accuracy and efficiency

Error decomposition guides adaptive refinement

Abstract

While many methods exist to discretize nonlinear time-dependent partial differential equations (PDEs), the rigorous estimation and adaptive control of their discretization errors remains challenging. In this paper, we present a methodology for duality-based a posteriori error estimation for nonlinear parabolic PDEs, where the full discretization of the PDE relies on the use of an implicit-explicit (IMEX) time-stepping scheme and the finite element method in space. The main result in our work is a decomposition of the error estimate that allows to separate the effects of spatial and temporal discretization error, and which can be used to drive adaptive mesh refinement and adaptive time-step selection. The decomposition hinges on a specially-tailored IMEX discretization of the dual problem. The performance of the error estimates and the proposed adaptive algorithm is demonstrated on two canonical applications: the elementary heat equation and the nonlinear Allen-Cahn phase-field model.

Tables

Table 1. Table 1: Effectivity of error estimate Eq.   54 under spatial refinement

$M$	$𝒬\left(u\right)-𝒬\left(I{u}_{\tau h}\right)$	${ℰ}_{\mathrm{st}}$	Effectivity
16	0.0003046	0.0002871	0.943
64	1.711e-05	1.606e-05	0.939
256	1.041e-06	9.756e-07	0.938
1024	6.458e-08	6.054e-08	0.937
4096	4.028e-09	3.776e-09	0.937

Table 2. Table 2: Effectivity of error estimate Eq.   54 under temporal refinement with uniform initial time steps and nonuniform initial time steps

	Uniform initial time steps			Nonuniform initial time steps
$N$	$𝒬\left(u\right)-𝒬\left(I{u}_{\tau h}\right)$	${ℰ}_{\mathrm{st}}$	Eff	$𝒬\left(u\right)-𝒬\left(I{u}_{\tau h}\right)$	${ℰ}_{\mathrm{st}}$	Eff
4	0.0001378	5.893e-05	0.428	0.0001058	4.773e-05	0.451
8	3.584e-05	1.659e-05	0.463	2.612e-05	1.253e-05	0.480
16	8.918e-06	4.272e-06	0.479	6.247e-06	3.054e-06	0.489
32	2.130e-06	1.032e-06	0.484	1.444e-06	7.063e-07	0.489
64	4.778e-07	2.303e-07	0.482	3.122e-07	1.507e-07	0.483

Table 3. Table 3: Effectivity of error estimate Eq.   54 under space-time refinement with uniform initial time steps and nonuniform initial time steps

	Uniform initial time steps			Nonuniform initial time steps
$M\times N$	$𝒬\left(u\right)-𝒬\left(I{u}_{\tau h}\right)$	${ℰ}_{\mathrm{st}}$	Eff	$𝒬\left(u\right)-𝒬\left(I{u}_{\tau h}\right)$	${ℰ}_{\mathrm{st}}$	Eff
256	6.061e-05	2.583e-05	0.426	4.184e-05	2.013e-05	0.481
2048	2.471e-05	1.086e-05	0.440	1.651e-05	7.581e-06	0.459
16384	7.926e-06	3.717e-06	0.469	5.295e-06	2.534e-06	0.479
131072	2.233e-06	1.082e-06	0.485	1.491e-06	7.299e-07	0.489

Equations

$$\int_{0}^{T}\Big{(}\langle\partial_{t}u,v\rangle+\mathcal{B}(u,v)+\mathcal{N}(u;v)\Big{)}\,\mathrm{d}t=\int_{0}^{T}\langle f,v\rangle\,\mathrm{d}t$$

$$\mathcal{Q}(u):=\left(\bar{q},u(T)\right)+{\int }_0^T(q,u)\mathrm{d}t,$$

$$\displaystyle\mathcal{N}^{s}(u,\hat{u};w,v)=\int_{0}^{1}\mathcal{N}^{\prime}\big{(}su+(1-s)\hat{u}\big{)}(w,v)\mathrm{d}s,\qquad\forall w,v\in\mathcal{V}$$

$${\mathcal{N}}^{\prime }(\hat{w})(w,v)=\underset{s\to 0}{\lim }\frac{\mathcal{N}\left(\hat{w}+sw;v\right)-\mathcal{N}\left(\hat{w};v\right)}{s}\mathrm{d}s,\forall w,v\in \mathcal{V}.$$

$${\mathcal{N}}^s(u,\hat{u};u-\hat{u},v)=\mathcal{N}(u;v)-\mathcal{N}(\hat{u};v).$$

$$\int_{0}^{T}\Big{(}\langle-\partial_{t}z,w\rangle+\mathcal{B}(z,w)+\mathcal{N}^{s}(u,\hat{u};w,z)\Big{)}\,\mathrm{d}t=\int_{0}^{T}(q,w)\,\mathrm{d}t$$

$${\mathcal{R}}_{\mathrm{PDE}}(\hat{u}(t);v)$$

$${\mathcal{R}}_0(\hat{u}(0);w)$$

$$\mathcal{Q}(u)-\mathcal{Q}(\hat{u})=\mathcal{R}_{0}\big{(}\hat{u}(0);z(0)\big{)}+\int_{0}^{T}\mathcal{R}_{\mathrm{PDE}}\left(\hat{u}(t);z(t)\right)\,\mathrm{d}t$$

$$\mathcal{N}(u;v)={\mathcal{N}}_c(u;v)-{\mathcal{N}}_e(u;v),$$

$$\left(\frac{u_{k+1}-u_k}{{\tau }_{k+1}},v\right)+\mathcal{B}(u_{k+1},v)+{\mathcal{N}}_c(u_{k+1};v)-{\mathcal{N}}_e(u_k;v)=(f_{k+1},v)$$

$$(u_0,v)=(u^0,v)\forall v\in L^2(\Omega ).$$

$${\mathcal{S}}_k^{h,p}:=\left\{v\in \mathcal{V}:v(x){|}_K\in {\mathbb{P}}^p(K),\forall K\in {\mathcal{K}}_k\right\}$$

$$\left(\frac{u_{k+1}^h-u_k^h}{{\tau }_{k+1}},v^h\right)+\mathcal{B}(u_{k+1}^h,v^h)+{\mathcal{N}}_c(u_{k+1}^h;v^h)-{\mathcal{N}}_e(u_k^h;v^h)=(f_{k+1},v^h)$$

$$(u_0^h,v^h)=(u^0,v^h)\forall v^h\in {\mathcal{S}}_0^{h,p}.$$

$$I{u}_{\tau h}(t)=\frac{t_{k+1}-t}{{\tau }_{k+1}}{u}_k^h+\frac{t-t_k}{{\tau }_{k+1}}{u}_{k+1}^{h}t\in [t_k,t_{k+1}].$$

$$⟨r_h^{k+1},v⟩:=$$

$$⟨r_{\tau }^{k+1}(t),v⟩:=$$

$$⟨r_f^{k+1}(t),v⟩:=$$

$${\mathcal{R}}_{\mathrm{PDE}}(I{u}_{\tau h}(t);v)$$

$$=⟨r_h^{k+1},v⟩+⟨r_{\tau }^{k+1}(t),v⟩+⟨r_f^{k+1}(t),v⟩$$

$$I{w}_{\tau }(\mathbf{x},t)=\sum _{k=0}^N{w}_k(\mathbf{x}){\mathbf{N}}_k(t)$$

$$\displaystyle\mathbf{N}_{k}(t):=\left\{\begin{array}[]{cl}{\displaystyle{\frac{t_{k+1}-t}{\tau_{k+1}}}}&\text{if}\;t\in\mathcal{I}_{k+1},k\leq N-1\\[8.0pt] {\displaystyle{\frac{t-t_{k-1}}{\tau_{k}}}}&\text{if}\;t\in\mathcal{I}_{k},k\geq 1\\ 0&\text{otherwise}\\ \end{array}\right.$$

$$q_k=\frac{1}{{\tau }_k}{\int }_0^{T}q{\mathbf{N}}_k(t)\mathrm{d}t\text{for }k=1,2,\dots ,N$$

$$q_0=\frac{1}{{\tau }_1}{\int }_0^{T}q{\mathbf{N}}_0(t)\mathrm{d}t.$$

$$\mathcal{Q}(I{w}_{\tau })={\tau }_1(q_0,w_0)+\sum _{k=1}^N{\tau }_k(q_k,w_k)+\left(\bar{q},w_N\right).$$

$$\mathcal{Q}(I{w}_{\tau })=\left(\bar{q},w_N\right)+{\int }_0^T\left(q,\sum _{k=0}^N{w}_k{\mathbf{N}}_k(t)\right)\mathrm{d}t=\left(\bar{q},w_N\right)+\sum _{k=0}^N\left({\int }_0^{T}q{\mathbf{N}}_k(t)\mathrm{d}t,w_k\right).$$

$$\sum _{k=0}^N\left({\int }_0^{T}q{\mathbf{N}}_k(t)\mathrm{d}t,w_k\right)=\left({\int }_0^{T}q{\mathbf{N}}_0(t),w_0\right)+\sum _{k=1}^N\left({\int }_0^{T}q{\mathbf{N}}_k(t)\mathrm{d}t,w_k\right),$$

$$-\left(\frac{z_1-z_0}{{\tau }_1},w\right)-{\mathcal{N}}_e^s(u_0,u_0^h;w,z_1)=(q_0,w)\forall w\in \mathcal{V}$$

$$\begin{gathered} -\left(\frac{z_{k+1}-z_k}{{\tau }_k},w\right)+\mathcal{B}(z_k,w)+{\mathcal{N}}_c^s(u_k,u_k^h;w,z_k) \\ -\frac{{\tau }_{k+1}}{{\tau }_k}{\mathcal{N}}_e^s(u_k,u_k^h;w,z_{k+1})=(q_k,w)\forall w\in \mathcal{V} \end{gathered}$$

$$\left(z_N,w\right)+{\tau }_N\mathcal{B}(z_N,w)+{\tau }_N{\mathcal{N}}_c^s(u_N,u_N^h;w,z_N)={\tau }_N(q_N,w)+(\bar{q},w)\forall w\in \mathcal{V}.$$

$$\mathcal{Q}(I{u}_{\tau })-\mathcal{Q}(I{u}_{\tau h})=\left(u^0-u_0^h,z_0-v_0^h\right)+\sum _{k=1}^N{\tau }_k⟨r_h^k,z_k-v_k^h⟩,$$

Full text

{AMS}

65M15, 65M20, 65M50

1 Introduction

Nonlinear parabolic PDEs are ubiquitous in science, however, their efficient numerical solution remains challenging. Implicit-explicit (IMEX) methods have been widely used for the time integration of complex time-dependent PDEs with terms of different type [1, 8]. Recently, a number of IMEX time-stepping schemes, paired with spatial Galerkin finite-element discretizations, have been proposed for phase-field models [35, 18, 37, 30, 34], which are currently a much-studied class of nonlinear parabolic problems [21, 24, 20, 28, 19]. When the PDE solution displays alternating fast and slow variations, the numerical discretization can, obviously, benefit significantly from adaptivity in both space and time.

This paper is devoted to the development of a posteriori error estimates and corresponding adaptive algorithms for these popular discretizations. In particular, we consider dual-based error estimates that assess the discretization error with respect to user specified quantities of interest describing the goal of the analyses. The quantities of interest might, for instance, be physical quantities or some appropriate norms of the error of the solution (e.g. energy norm, $L^{2}$ norm). To efficiently drive adaptive mesh refinement and adaptive time-step selection, the error estimates need to address the temporal and the spatial discretization errors separately.

There have been several studies on goal-oriented adaptive techniques for parabolic equations during the last decade, but mostly in the context of space-time (discontinuous) Galerkin finite element discretization, see for instance Eriksson and Johnson [13, 14, 15, 16], Schmich and Vexler [27], Carey et al. [10], Bermejo and Carpio [6], Braack et al. [9], Besier and Rannacher[7], and Asner et al. [2]. Very little progress has been made for parabolic equations discretized using IMEX time-stepping schemes.

Recently, Chaudhry et al. [11, 12] proposed a posteriori error estimates for various IMEX schemes, based on an equivalence relation between IMEX schemes and time-Galerkin finite element methods. They rewrite the time-Galerkin method using special numerical quadrature rules and carry out a standard duality-based analysis for the resultant approximations. The splitting of the temporal and the spatial error contributions in these error estimates are commonly achieved by inserting and subtracting suitable projections of the dual solution.

The objective of this paper is to present an alternative approach to duality-based a posteriori error estimates for fully discretized semi-linear parabolic PDEs using conforming finite elements in space and first-order IMEX schemes in time. Contrary to Chaudhry et al, in our approach we directly obtain a posteriori error estimates without resorting to an interpretation of IMEX as a Galerkin-in-time method. This paper is a follow-up to our recent paper [29], where we only considered errors due to spatial discretization. The focus of this work is on the total discretization error which contains both the spatial and temporal parts.

The starting point of our analysis is the exact duality-based error representation, which is a duality pairing of the global space-time residual with the solution of the mean-value-linearized (backward-in-time) dual problem. This error representation can be decomposed into various distinct residuals weighted by the same dual solution. A fundamental framework for successfully decomposing the residuals for (non)linear parabolic PDEs, discretized by a classical A-stable $\theta$-scheme in time, has been developed by Verfürth [33] in the context of energy-based a posteriori error analysis.

By extending Verfürth’s framework to IMEX schemes and a duality-based error analysis, we will decompose the error representation into three contributions which can be associated to the temporal and spatial discretization error, and additionally data oscillation. This novel decomposition hinges on a special nonstandard IMEX discretization of the dual problem. We then propose a general space-time adaptive algorithm for an efficient distribution of the discretization parameters: a set of time steps and the refined mesh at each time step.

This work is structured as follows. In Section 2, we introduce the abstract setting for a general (non)linear parabolic PDE and its IMEX-Galerkin discretization. Section 3 is devoted to the methodology for a space-time decomposition of a duality-based a posteriori error estimate. After having established computable error estimates in Section 4, we propose the associated adaptive algorithm in Section 5. The application to the elementary heat equation and the nonlinear Allen–Cahn equation (an elementary phase-field model), together with numerical results are presented in Section 6, after which we present our conclusions.

2 Abstract setting

In this section, we start by introducing an abstract setting of nonlinear parabolic PDEs and the corresponding dual problem in weak formulations. Then we present the discretization of the primal problem using IMEX time-stepping schemes and conforming finite elements in space.

2.1 Weak formulation and error representation

As a model problem we consider a general semi-linear parabolic equation in a bounded space-time domain $\Omega\times(0,T]$ with $\Omega\subset\mathbb{R}^{d}$, $d=1,2,3$, having natural boundary conditions. To provide a setting for the weak formulation, we denote by $\mathcal{V}$ a suitable Hilbert space and by $\mathcal{V}^{*}$ its dual space, such that $\mathcal{V}\subset L^{2}(\Omega)\subset\mathcal{V}^{*}$ with continuous embeddings. We denote the inner product in $L^{2}(\Omega)$ by $(\cdot,\cdot)$, and duality pairings between $\mathcal{V}^{*}$ and $\mathcal{V}$ by $\left\langle\cdot,\cdot\right\rangle$. By defining $\mathcal{W}:=\left\{v\in L^{2}(0,T;\mathcal{V}),\partial_{t}v\in L^{2}(0,T;\mathcal{V}^{*})\right\}$ as a suitable space for $u$, the weak form reads: find $u\in\mathcal{W}_{u^{0}}:=\left\{v\in\mathcal{W}:v(0)=u^{0}\right\}$ such that $\forall v\in L^{2}(0,T;\mathcal{V})$

[TABLE]

where $f\in L^{2}(0,T;\mathcal{V}^{*})$, $u^{0}\in L^{2}(\Omega)$, the semi-linear form $\mathcal{N}(\cdot;\cdot)$ of a sufficiently smooth nonlinear operator represents the nonlinear components which is linear with respect to arguments on the right of the semicolon, and $\mathcal{B}(\cdot,\cdot)$ is the bilinear form of a elliptic self-adjoint operator. A prime example of the abstract setting is the Allen–Cahn equation $\partial_{t}u-\Delta u+\frac{1}{\varepsilon^{2}}\psi^{\prime}(u)=0$, which will be discussed later in Section 6.1.

Given the solution $u$, we consider the quantity

[TABLE]

with $\bar{q}\in L^{2}(\Omega)$ and $q\in L^{2}(0,T;L^{2}(\Omega))$ so that $\mathcal{Q}:\mathcal{W}\to\mathbb{R}$ is a continuous linear functional.111One can more generally consider $q\in L^{2}(\tau_{1},T;\mathcal{V}^{*})\cup L^{2}(0,T;L^{2}(\Omega))$ where $\tau_{1}>0$ is the size of the first time step. For technical reasons later on (i.e., $\langle q_{0},u_{0}\rangle$ must be well-defined, with $q_{0}$ defined in Eq. 21), we can not take $\tau_{1}=0$. One example of $\mathcal{Q}(u)$ would be the value of the solution at the final time $t=T$ at a critical area of the domain centered at $\mathbf{x}_{0}$, $\mathcal{Q}(u)=\int_{\Omega}\rho_{\epsilon}(\mathbf{x}_{0}-\mathbf{x})u(\mathbf{x},T)\,\mathrm{d}\mathbf{x},$ where $\rho_{\epsilon}\in C^{\infty}$ is a kernel function with radius and center of $\epsilon$ and $\mathbf{x}_{0}$. Alternatively, one might wish to estimate the error in the $L^{2}$ norm at the final time $T$. To achieve this, we set $\bar{q}=u(T)-\hat{u}(T)$ and $q=0$ where $\hat{u}$ is an approximation of the solution $u$. Then we have $\mathcal{Q}(u)-\mathcal{Q}(\hat{u})={\|u(T)-\hat{u}(T)\|}^{2}_{L^{2}(\Omega)}$.

For any $u,\hat{u}\in\mathcal{V}$, we denote by $\mathcal{N}^{s}(u,\hat{u};\cdot,\cdot)$ the mean-value linearization of $\mathcal{N}(\cdot,\cdot)$ performed at a value in between $u$ and $\hat{u}$, namely,

[TABLE]

where $\mathcal{N}^{\prime}$ is the Gâteaux derivative of $\mathcal{N}$, i.e.

[TABLE]

Note that if we set $w=u-\hat{u}$, the chain rule gives

[TABLE]

The mean-value-linearized (backward-in-time) dual problem takes the form: find $z\in\mathcal{W}^{\bar{q}}:=\left\{v\in\mathcal{W}:v(T)=\bar{q}\right\}$ such that $\forall w\in L^{2}(0,T;\mathcal{V})$

[TABLE]

Let $\hat{u}\in\mathcal{W}$ denote any approximation of the solution $u$ in Eq. 1. We define the residual of the primal PDE, $\mathcal{R}_{\mathrm{PDE}}$, and the residual of the initial condition, $\mathcal{R}_{0}$, as

[TABLE]

for all $v\in\mathcal{V}$ and $w\in L^{2}(\Omega)$. Following the general framework of goal-oriented error analysis (see, e.g. [5, 25]), we obtain an exact error representation assessing the error in $\mathcal{Q}$, which can generally be represented by a global space-time residual weighted by the solution of the dual problem Eq. 5.

Theorem 2.1** (Global space-time error representation).**

Given any approximation $\hat{u}\in\mathcal{W}$ of the solution $u$ of the primal problem Eq. 1, we have the following a posteriori error representation:

[TABLE]

*where $z\in\mathcal{W}^{\bar{q}}$ is the solution of the dual problem Eq. 5. *

Proof 2.2**.**

*The proof is standard, see, e.g. [32] or [29, Theorem 2.A]. *

Note that errors in norm are also included in Theorem 2.1 by suitably changing $\bar{q}$ and $q$, e.g., as in the example above.

2.2 IMEX - FEM Discretization

We next describe a full discretization of problem Eq. 1 by partitioning $[0,T]$ as $0=t_{0}<t_{1}<t_{2}<\cdots<t_{k}<\cdots<t_{N}=T$ into $N$ subintervals $\mathcal{I}_{k+1}=[t_{k},t_{k+1}]$ of length $\tau_{k+1}=t_{k+1}-t_{k}$, $k=0,1,\cdots,N-1$. Because of the nonlinearity in the system Eq. 1, one has to be careful in choosing a time discretization to avoid prohibitive stability restrictions and high computational complexity. In this paper, we focus on first-order IMEX time-stepping schemes, which employ a splitting of the nonlinear term $\mathcal{N}$ according to

[TABLE]

The notation $\mathcal{N}_{c}$ and $\mathcal{N}_{e}$ comes from the phase-field modeling community, and refers to the contractive and expansive part, respectively, which can also refer to the stiff and non-stiff term. The fundamental idea is to treat the contractive part implicitly and the expansive part explicitly. Such a time scheme for problem Eq. 1 is defined recursively by: find $u_{k+1}\in\mathcal{V}$ such that $\forall v\in\mathcal{V}$

[TABLE]

for $k=0,1,\cdots,N-1$, where the initial condition is

[TABLE]

Here, $f_{k+1}=f(\cdot,t_{k+1})$, which is well-defined upon assuming that the function $f$ is sufficiently regular, e.g., $f\in C^{0}((0,T];L^{2}(\Omega))$. We remark that instead of the time approximation $f_{k+1}$, a time-averaged approximation $\bar{f}=\frac{1}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}f(\cdot,t)\,\mathrm{d}t$ can be used, provided that $f\in L^{2}(0,T;L^{2}(\Omega))$. We also implicitly assume in Eq. 9 that $\mathcal{N}_{e}(u_{0};v)$ is bounded for $u_{0}\in L^{2}(\Omega)$. If $\mathcal{N}_{e}(u_{0};v)$ is not well-defined, one can remove this term from Eq. 9 for the first time step. For simplicity, we continue our analysis assuming that $\mathcal{N}_{e}(\cdot;\cdot)$ is bounded on $L^{2}(\Omega)\times L^{2}(\Omega)$.

To fully discretize the primal problem Eq. 1, we consider a standard shape-regular mesh $\mathcal{K}_{k}$ of $\Omega$ and an associated conforming finite element space $\mathcal{S}^{h,p}_{k}$ defined by

[TABLE]

for $k=0,1,\ldots,N$, where $\mathbb{P}^{p}(K)$ is the space of polynomials up to order $p$ on element $K$ and $h$ denotes the mesh parameter. The fully discrete approximation is then formulated as: find $u_{k+1}^{h}\in\mathcal{S}^{h,p}_{k+1}$ such that $\forall v^{h}\in\mathcal{S}^{h,p}_{k+1}$

[TABLE]

for $k=0,1,\cdots,N-1$, where the initial condition is

[TABLE]

We assume that the solutions $u_{\tau}:=\{u_{k}\}_{k=0}^{N}$ and $u_{\tau h}:=\{u_{k}^{h}\}_{k=0}^{N}$ exist for the time-discrete primal problem Eq. 9 and the fully-discrete primal problem Eq. 11, respectively.

3 Space-time decomposed a posteriori error estimate

Space-time adaptivity is heavily dependent on an appropriate decomposition of error estimates, which will be derived in this section. Our approach to isolate error contributions from different sources is inspired by the work of Verfürth in [33, Chapter 6], which contains a general framework for deriving residual-based a posteriori error estimates for nonlinear parabolic problems with the $\theta$-scheme. In the following Lemma, we adapt Verfürth’s residual decomposition to our fully discrete primal problem Eq. 11.

Lemma 3.1** (Residual decomposition).**

Let $u_{\tau h}:=\{u_{k}^{h}\}_{k=0}^{N}$ denote the solution of the fully discrete problem Eq. 11, and $Iu_{\tau h}$ denote the piecewise-linear time reconstruction of $u_{\tau h}$ on time intervals $[t_{k},t_{k+1}]$, $k=0,1,\ldots,N-1$, i.e.,

[TABLE]

Let the spatial residual $r_{h}^{k}$, the temporal residual $r_{\tau}^{k}$ and the data-oscillation contribution $r_{f}^{k}$ be defined, for each $k=0,1,\cdots,N-1$, by

[TABLE]

for all $v\in\mathcal{V}$ and $t\in(t_{k},t_{k+1}]$. Then, for each $k=0,1,\cdots,N-1$, the following decomposition of the space-time residual Eq. 6 holds:

[TABLE]

*where $t\in(t_{k},t_{k+1}]$. *

Proof 3.2**.**

*Since $\partial_{t}Iu_{\tau h}=\frac{u^{h}_{k+1}-u^{h}_{k}}{\tau_{k+1}}$ on $(t_{k},t_{k+1}]$, the identities in Eq. 18 follow from a straightforward substitution in Eq. 17 using the definition Eq. 14, Eq. 15 and Eq. 16. *

Remark 3.3**.**

*We note that the spatial residuals Eq. 14 are independent of time, and due to Galerkin orthogonality, the spatial residuals will be equal to zero if $v\in\mathcal{S}^{h,p}_{k+1}$. Furthermore, upon convergence $u_{k}^{h}\to u_{k}$ as $h\to 0^{+}$, for all $k$, we also have $r_{h}^{k+1}\to 0$ (see Eq. 9). Similarly, assuming sufficient smoothness in time, then $u_{k}^{h},u_{k+1}^{h}\to Iu_{\tau h}(t)$ for $t\in[t_{k},t_{k+1}]$ as $\tau_{k+1}\to 0^{+}$, which implies $r_{\tau}^{k+1}(t)\to 0$ as $\tau_{k+1}\to 0^{+}$. This is the motivation for calling $r_{\tau}^{k+1}$ and $r_{h}^{k+1}$ the temporal residual and the spatial residual, respectively. *

3.1 Time-discrete error representation

The first step toward a decomposition of duality-based error estimates is to introduce a time-discrete error representation identifying only the spatial discretization error. To this end, we introduce a novel and specially-tailored IMEX time-discrete dual problem. This time-discrete problem is driven by the following discrete representation of $\mathcal{Q}$.

Let us rewrite the piecewise-linear time reconstruction $Iw_{\tau}\in\mathcal{W}$ of any sequence $w_{\tau}:=\{w_{k}\}_{k=0}^{N}$, $w_{k}\in\mathcal{V}$, as

[TABLE]

where

[TABLE]

We consider the following discrete representation of $\mathcal{Q}:\mathcal{W}\to\mathbb{R}$ when applied to $Iw_{\tau}$.

Lemma 3.4**.**

Let us define

[TABLE]

and

[TABLE]

Then, the following time-discrete representation of $\mathcal{Q}:\mathcal{W}\to\mathbb{R}$ holds

[TABLE]

Proof 3.5**.**

For $Iw_{\tau}$ defined in Eq. 19, we observe that, according to Eq. 2,

[TABLE]

By virtue of

[TABLE]

*we obtain Eq. 22 by substituting the definition Eq. 20 and Eq. 21. *

We now state the novel IMEX time-stepping scheme to discretize the dual problem backwards in time: Find $z_{k}\in\mathcal{V}$, $k=0,1,\ldots,N,$ such that

[TABLE]

and for $k=1,2,\ldots,N-1$:

[TABLE]

where the terminal condition is

[TABLE]

The time-discrete dual Eq. 23-Eq. 25 has been defined so as to provide an exact error representation for $\mathcal{Q}(Iu_{\tau h})$ with respect to $\mathcal{Q}(Iu_{\tau})$.

Theorem 3.6** (Time-discrete error representation).**

Let $u_{\tau}=\{u_{k}\}_{k=0}^{N}$ denote the solution of the time discrete system Eq. 9, and $u_{\tau h}=\{u_{k}^{h}\}_{k=0}^{N}$ denote the solution of the fully discrete system Eq. 11. Let $z_{\tau}=\{z_{k}\}_{k=0}^{N}$ denote the time discrete approximation of the dual problem obtained from Eq. 23-Eq. 25. Then the following error representation holds:

[TABLE]

*for any $v_{k}^{h}\in\mathcal{S}^{h,p}_{k}$, $k=0,1,\ldots,N$. *

Proof 3.7**.**

From Eq. 22, it follows that $\mathcal{Q}(Iu_{\tau})-\mathcal{Q}(Iu_{\tau h})$ can be formulated as

[TABLE]

Substituting the time-discrete dual problem Eq. 23–Eq. 25 into Eq. 27, we get

[TABLE]

After applying summation by parts on $(z_{k+1}-z_{k},u_{k}-u_{k}^{h})$, i.e.,

[TABLE]

it follows that

[TABLE]

Then, by shifting the indices of the arguments of $\mathcal{B}$ and $\mathcal{N}^{s}_{c}$:

[TABLE]

and employing the mean-value linearization property Eq. 4 on $\mathcal{N}_{c}^{s}$ and $\mathcal{N}_{e}^{s}$, we arrive at

[TABLE]

After substituting the time-discrete primal problem Eq. 9 weighted by dual solution $z_{k+1}$, we finally obtain

[TABLE]

*This is Eq. 26 by the definition in Eq. 14, and noting that $\forall v_{k+1}^{h}\in\mathcal{S}^{h,p}_{k+1}$ it holds that $\langle r_{k+1}^{h},v_{k+1}^{h}\rangle=0$ by Eq. 11 for $k=0,1,\ldots,N-1$ and $\left(u^{0}-u_{0}^{h},v_{0}^{h}\right)=0$ $\forall v_{0}^{h}\in\mathcal{S}^{h,p}_{0}$ by Eq. 12. *

Remark 3.8**.**

Note that there is an alternative discrete representation of $\mathcal{Q}$ to the one in Eq. 22:

[TABLE]

where

[TABLE]

Instead of Eq. 23-Eq. 25, one then would expect an alternative time-stepping scheme for the dual problem in Eq. 5 (solved backwards in time) as: find $z_{k}\in\mathcal{V}$, $k=0,1,\ldots,N$, such that

[TABLE]

and for $k=1,2,\ldots,N-1$:

[TABLE]

where the terminal condition is

[TABLE]

*However, this alternative time scheme is equivalent to Eq. 23-Eq. 25, and therefore leads to exactly the same time-discrete error representation Eq. 26. *

3.2 Spatial and temporal error representation

Building on Verfürth’s residual decomposition Eq. 18 and the time-discrete error representation Eq. 26, we are now ready to state our main result: A suitable decomposition of the dual-weighted residual Eq. 8.

Theorem 3.9** (Decomposed error representation).**

Let the assumptions of Theorem 3.6 hold. Let $u$ denote the solution of the primal problem Eq. 1 and $z$ denote the solution of the dual problem Eq. 5. Then the following error representation holds:

[TABLE]

for any $v_{\tau}^{h}:=\{v_{k}^{h}\}_{k=0}^{N}$, $v_{k}^{h}\in\mathcal{S}^{h,p}_{k}$, where $\mathcal{R}_{\mathrm{s}}(u_{\tau h};z_{\tau}-v_{\tau}^{h})$ is the spatial error representation

[TABLE]

$\mathcal{R}_{\mathrm{t}}(u_{\tau h},u_{\tau},z,z_{\tau})$* is the temporal error representation*

[TABLE]

and $\operatorname{Osc}$ denotes the data-oscillation contribution

[TABLE]

Remark 3.10**.**

*By virtue of Galerkin orthogonality, $\mathcal{R}_{\mathrm{s}}(u_{\tau h};z_{\tau})$ will vanish if $z_{k}\in\mathcal{S}^{h,p}_{k}$, $k=0,1,\ldots,N$. In addition, as $h\to 0^{+}$ it holds that $u_{k}^{h}\to u_{k}$ for all $k$ and, accordingly, $\mathcal{R}_{\mathrm{s}}(u_{\tau h};\cdot)\to 0$ (see Eq. 9 and Eq. 10). Similarly, assuming sufficient smoothness in time, then $z_{k},z_{k+1}\to z(t)$ for $t\in[t_{k},t_{k+1}]$ as $\tau_{k+1}\to 0^{+}$, which implies $\left(u^{0}-u_{0}^{h},z(0)-z_{0}\right)\to 0$ and $\mathcal{R}_{\mathrm{PDE}}(Iu_{\tau h}(t);z(t)-z_{k+1})\to 0$ as $\tau_{k+1}\to 0^{+}$. And since $r_{\tau}^{k+1}(t)\to 0$ as $\tau_{k+1}\to 0^{+}$ (see Remark 3.3), we conclude that $\mathcal{R}_{\mathrm{t}}(u_{\tau h},u_{\tau},z,z_{\tau})\to 0$ as $\tau_{k+1}\to 0^{+}$. This is the motivation for calling $\mathcal{R}_{\mathrm{t}}(u_{\tau h},u_{\tau},z,z_{\tau})$ and $\mathcal{R}_{\mathrm{s}}(u_{\tau h};z_{\tau})$ the temporal error representation and the spatial error representation, respectively. *

Remark 3.11**.**

*If we choose $f_{k+1}=\frac{1}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}f(t)\,\mathrm{d}t$, then the data-oscillation contribution Eq. 31 will vanish. *

Proof 3.12**.**

(of Theorem 3.9)  The global space-time error representation is Theorem 2.1 with $\hat{u}=Iu_{\tau h}$:

[TABLE]

The spatial error representation Eq. 29 satisfies the representation in Theorem 3.6.

The temporal error representation is obtained by subtracting the spatial error representation Eq. 29 and the data-oscillation contribution Eq. 31 from the space-time error representation Eq. 32, i.e.,

[TABLE]

Adding and subtracting $\mathcal{R}_{\mathrm{PDE}}(Iu_{\tau h};z_{k+1})$ yields

[TABLE]

Since \partial_{t}Iu_{\tau h}=\big{(}u^{h}_{k+1}-u^{h}_{k}\big{)}/\tau_{k+1} on $(t_{k},t_{k+1}]$ according to the definition of $Iu_{\tau h}$ Eq. 19, we employ the definition of the residuals in Eq. 6, Eq. 7 and Eq. 14, and obtain

[TABLE]

*Finally, substituting the definitions in Eq. 15, Eq. 16 and Eq. 31 gives the result Eq. 30. *

A useful interpretation of the spatial and the temporal error representation can be obtained by writing the global space-time error as:

[TABLE]

The following Corollary holds:

Corollary 3.13**.**

Under the assumptions of Theorem 3.6 and Theorem 3.9, we have

[TABLE]

and

[TABLE]

*for any $v_{\tau}^{h}:=\{v_{k}^{h}\}_{k=0}^{N}$, $v_{k}^{h}\in\mathcal{S}^{h,p}_{k}$. *

Proof 3.14**.**

*The identity in Eq. 35 is a direct consequence of Eq. 26 and Eq. 29. Equation Eq. 34 then follows from Eq. 28 and Eq. 29. *

4 Computable error estimate

There are two approximations commonly involved in evaluating the exact error representations Eq. 32, Eq. 29 and Eq. 30:

• •

one for approximating the exact primal solutions $u$ and $u_{\tau}$ in the mean-value-linearized dual problem Eq. 5 and its time-discrete system Eq. 23-Eq. 25,

• •

the other for approximating the exact dual solutions $z$ and $z_{\tau}$ in the error representation formulas Eq. 32, Eq. 29 and Eq. 30.

The resulting error estimate can only be accurate if the approximations are sufficiently close to the true solutions.

Here, to obtain computable and asymptotically effective error estimates, we consider a hierarchical two-level methodology developed in [29] where the estimate is directly evaluated with an enriched dual approximation that is computed with help of an additional primal approximation at an enriched discretization level for the mean-value-linearization. Since the focus of [29] is on the spatial discretization error, we now extend this methodology to our problem: We need two additional discretization levels: one which is spatially-enriched for evaluating the spatial error representation Eq. 29 and the other which is space-time enriched for evaluating Eq. 30 and Eq. 32.

We first introduce the following notations: {remunerate}

$\mathcal{S}^{h,p}_{k}$: the original FE space with spatial mesh of size $h=h_{k}$ at time $t_{k}$ for $k=0,1,\ldots,N$.

$\mathcal{S}^{h/2,p}_{k}$: an enriched FE space with finer spatial mesh of size ${\displaystyle{\frac{h_{k}}{2}}}$ at time $t_{k}$ for $k=0,1,\ldots,N$. $\mathcal{S}^{h/2,p}_{k}$ is obtained by global refinement of all the element in $\mathcal{S}^{h,p}_{k}$.

$\mathcal{S}^{h/2,p}_{k+1/2}$: an enriched FE space with finer spatial mesh of size ${\displaystyle{\frac{h_{k}}{2}}}$ at the intermediate time level ${\displaystyle{t_{k+1/2}=\frac{t_{k+1}+t_{k}}{2}}}$ for $k=0,1,\ldots,N-1$. We set $\mathcal{S}^{h/2,p}_{k+1/2}=\mathcal{S}^{h/2,p}_{k+1}$.

$u_{\tau,h}=\{u_{k}^{h}\}_{k=0}^{N}$: the solution of Eq. 11 and Eq. 12 using time-step sizes ${\displaystyle{\{\tau_{k}\}_{k=1}^{N}}}$ and FE spaces ${\displaystyle{\{\mathcal{S}^{h,p}_{k}\}_{k=0}^{N}}}$.

$u_{\tau,h/2}=\{u_{k}^{h/2}\}_{k=0}^{N}$: the solution of Eq. 11 and Eq. 12 using time-step sizes ${\displaystyle{\{\tau_{k}\}_{k=1}^{N}}}$ and enriched FE spaces ${\displaystyle{\{\mathcal{S}^{h/2,p}_{k}\}_{k=0}^{N}}}$; $u_{\tau,h/2}$ represents an approximation of the time-discrete primal solution $u_{\tau}$.

$u_{\tau/2,h/2}=\{\tilde{u}_{l}^{h/2}\}_{l=0,1/2,\ldots,N}$: the solution of Eq. 11 and Eq. 12 using half time-step sizes ${\displaystyle{\left\{\frac{\tau_{1}}{2},\frac{\tau_{1}}{2},\frac{\tau_{2}}{2},\frac{\tau_{2}}{2},\ldots,\frac{\tau_{N}}{2},\frac{\tau_{N}}{2}\right\}}}$ and enriched FE spaces $\{\mathcal{S}^{h/2,p}_{l}\}_{l=0,1/2,\ldots,N}$; $u_{\tau/2,h/2}$ represents an approximation of the exact primal solution $u$.

$z_{\tau,h/2}=\{z_{k}^{h/2}\}_{k=0}^{N}$: the solution of the approximate dual problem obtained by replacing $u_{\tau}$ with $u_{\tau,h/2}$ in Eq. 23-Eq. 25, using time-step sizes ${\displaystyle{\{\tau_{k}\}_{k=1}^{N}}}$ and enriched FE spaces ${\displaystyle{\{\mathcal{S}^{h/2,p}_{k}\}_{k=0}^{N}}}$; $z_{\tau,h/2}$ represents an approximation of the time-discrete dual $z_{\tau}$.

$z_{\tau/2,h/2}=\{\tilde{z}_{l}^{h/2}\}_{l=0,1/2,\ldots,N}$: the solution of the approximate dual problem obtained by replacing $u_{\tau}$ with $u_{\tau/2,h/2}$ in Eq. 23-Eq. 25, using half time-step sizes ${\displaystyle{\left\{\frac{\tau_{1}}{2},\frac{\tau_{1}}{2},\frac{\tau_{2}}{2},\frac{\tau_{2}}{2},\ldots,\frac{\tau_{N}}{2},\frac{\tau_{N}}{2}\right\}}}$ and enriched FE spaces $\{\mathcal{S}^{h/2,p}_{l}\}_{l=0,1/2,\ldots,N}$; $z_{\tau/2,h/2}$ represents the approximation of the exact dual solution $z$.

The strategy for computing the primal and dual solutions is illustrated in Figure 1. For evaluating the error representations Eq. 32, Eq. 29 and Eq. 30, we compute two enriched dual solutions $z_{\tau,h/2}$ and $z_{\tau/2,h/2}$ solved backwards in time to approximate $z_{\tau}$ and $z$, respectively. In order to make $z_{\tau,h/2}$ computable, an additional primal approximation $u_{\tau,h/2}$ is computed forwards in time using FE spaces $\{\mathcal{S}^{h/2,p}_{k}\}_{k=0}^{N}$ to approximate the mean-value-linearization of the dual problem Eq. 23-Eq. 25. Similarly, another additional primal approximation $u_{\tau/2,h/2}$ is computed using spaces $\{\mathcal{S}^{h/2,p}_{l}\}_{l=0,1/2,\ldots,N}$ for obtaining $z_{\tau/2,h/2}$.

Now let us denote by $\hat{z}\in\mathcal{V}$ a time-reconstruction of the dual solution $z_{\tau/2,h/2}$ (e.g. a piecewise-constant time-reconstruction will be used in numerical applications; see Section 6). By replacing $z$ with $\hat{z}$ in Eq. 32, the estimate of the space-time error in $\mathcal{Q}$ can then be computed as:

[TABLE]

Replacing $z_{\tau}$ with the computable $z_{\tau,h/2}$ in Eq. 29, we compute the spatial error estimate $\mathcal{E}_{\mathrm{s}}$ as:

[TABLE]

Finally, replacing $z$ and $z_{\tau}$ with the computable $\hat{z}$ and $z_{\tau,h/2}$ in Eq. 30, respectively, we compute the temporal error estimate $\mathcal{E}_{\mathrm{t}}$ as:

[TABLE]

Remark 4.1** (Reduced-cost implementation).**

*If one wants to reduce the number of distinct approximations in the error estimates Eq. 36–Eq. 38, the most straightforward strategy is to simply take $u_{\tau,h/2}=u_{\tau/2,h/2}$ and $z_{\tau,h/2}=z_{\tau/2,h/2}$ at concurrent time steps (i.e. $u_{k}^{h/2}=\tilde{u}_{k}^{h/2}$ and $z_{k}^{h/2}=\tilde{z}_{k}^{h/2}$ for $k=0,1,\ldots,N$). In this manner, one only needs to compute $u_{\tau h}$, $u_{\tau/2,h/2}$ and $z_{\tau/2,h/2}$, without the gray columns in Figure 1. More detailed description and examples are given in numerical applications; see Section 6. An even cheaper alternative is to compute higher-order reconstructions using only $u_{\tau h}$ and $z_{\tau h}$; see, Becker and Rannacher [4] for an overview, or coarse-scale adjoints [10]. *

5 Adaptive algorithm

Our goal is now to design an adaptive algorithm to iteratively increase the accuracy of the numerical solution by using the error estimates. In this section, we first derive error indicators of local contributions that serve as the basis to control adaptive mesh refinement and adaptive time-step selection, and then present the space-time adaptive algorithm.

5.1 Error indicators

To drive space-time adaptivity, the information of the global error estimates has to be localized to time-intervals and spatially-local contributions. To this end, we rewrite the computable error estimates $\mathcal{E}_{\mathrm{st}}$, $\mathcal{E}_{s}$ and $\mathcal{E}_{\mathrm{t}}$ in Eq. 36–Eq. 38 as a sum of their local contributions on each time intervals $[t_{k},t_{k+1}]$, $k=0,1,\ldots,N-1$, respectively. The absolute values of these local contributions are identified as the local indicators, which can directly be used for adaptive time-step selection. For adaptive mesh refinement, the local contributions associated to the spatial discretization error have to be localized further in space. We summarize the result in the following propositions.

Proposition 5.1**.**

The error estimate $\mathcal{E}_{\mathrm{st}}$, $\mathcal{E}_{\mathrm{t}}$ and $\mathcal{E}_{\mathrm{s}}$ can be bounded from above by

[TABLE]

where the local space-time error indicators $\mathcal{E}^{0}_{h\tau}$ and $\mathcal{E}^{k+1}_{h\tau}$ are defined by

[TABLE]

the temporal error indicators $\mathcal{E}^{0}_{\tau}$ and $\mathcal{E}^{k+1}_{\tau}$ are defined by

[TABLE]

and the spatial error indicators $\mathcal{E}^{0}_{h}$ and $\mathcal{E}^{k+1}_{h}$ are defined by

[TABLE]

Proof 5.2**.**

We split the error estimator $\mathcal{E}_{\mathrm{st}}$ Eq. 36 into local space-time error indicators Eq. 39 by

[TABLE]

*Following the same procedure as above, the temporal error estimate $\mathcal{E}_{\mathrm{t}}$ Eq. 38 is localized to the temporal error indicators LABEL:eq:est_time on each time intervals $[t_{k},t_{k+1}]$, $k=0,1,\ldots,N-1$, and the spatial error estimate $\mathcal{E}_{\mathrm{s}}$ Eq. 37 is localized in time to the spatial error indicators Eq. 42. *

For spatial adaptivity we consider hierarchical mesh-refinement indicators as explained in [29, Section 4.2]. In addition, let us note that instead of a traditional element-wise marking strategy, we use the function-support marking strategy introduced in [23] (see also [26, 31]).

Remark 5.3**.**

*In general in Eq. 42, one needs to subtract an interpolant of the dual solution, $\Pi z_{\tau h/2}$ from $z_{\tau,h/2}$ to get a sharp spatial indicator. However, this is not needed for hierarchical indicators; see, [29, Section 4.2] for more details. *

5.2 The space-time adaptive algorithm

In Algorithm 1, we propose a global space-time adaptive procedure using the above duality-based indicators. The pseudocode consists of three parts: (1) The computation of the primal and dual approximations (comprised of lines 3-8, 9-11 and 13-14), (2) the evaluation of the error estimates (given in lines 12, 15, 18-19 and 23), and (3) the error control (comprised of the remaining lines of Algorithm 1).

Within the adaptive procedure, the error control is built on a two-step approach. First, in lines 16-17 of the pseudocode we apply the maximum marking strategy (following Babuška and Vogelius [3]) with fraction $\theta\in[0,1]$ on space-time error indicators $\{\mathcal{E}^{k}_{h\tau}\}_{k=0}^{N}$ to globally select time steps $\{k^{\prime}\}$ (i.e. $|\mathcal{E}^{k^{\prime}}_{h\tau}|\geq\theta\,\max\{|\mathcal{E}^{0}_{h\tau}|,|\mathcal{E}^{1}_{h\tau}|,\ldots,|\mathcal{E}^{N}_{h\tau}|\}$), which contain the largest error contributions throughout the time period. Second, in line 20 we locally check the leading causes of the error at the targeted time steps $\{k^{\prime}\}$, whether from the spatial error or from the temporal error. An example of this is shown and explained in Figure 2. The figure on the left indicates the time steps $\{k^{\prime}\}$, where the major error contributions are located. Then, if the spatial indicator $\mathcal{E}^{k^{\prime}}_{h}$ is larger than the temporal indicator $\mathcal{E}^{k^{\prime}}_{\tau}$, the spatial mesh is targeted for refinement according to the mesh indicators $\mathcal{E}_{i}^{k^{\prime}}$; see Figure 2 (center). Otherwise, the time step size $\tau_{k^{\prime}}$ is marked and reduced by half ; see Figure 2 (right).

The adaptive spatial mesh refinement is also based on a maximum marking strategy. The nodes $\{i^{\prime}\}$ are marked for which their mesh-refinement indicators are at least a fraction $\lambda\in[0,1]$ of the maximal mesh indicator (i.e. $|\mathcal{E}_{i^{\prime}}^{k}|\geq\lambda\max\{|\mathcal{E}_{0}^{k}|,\ldots,|\mathcal{E}_{M}^{k}|\}$). The addition of the basis function on selected nodes is performed using hierarchical refinement for finite element methods [22, 23, 26, 31]. Moreover, instead of projection, we introduce a common refinement to transfer the solution from one mesh to another without loss of accuracy in any quadrature approximations.

Remark 5.4**.**

*A standard adaptive algorithm for time-dependent problems starts with an initial coarse mesh, and proceeds sequentially. Based on the mesh for the current time step, a new space mesh is generated for each new time step. Such a sequential procedure commonly uses residual-based error estimates in space which only contain information at the current time step. Duality-based error estimates, on the contrary, contain the entire evolution history of the error dependence implicitly via the dual solution. *

6 Applications

In this section, we give two examples of problems which fit into the abstract framework introduced in Section 2: the nonlinear Allen–Cahn equation and the linear heat equation (as a special case of the Allen–Cahn equation). We numerically investigate the performance of the duality-based error estimates and the proposed adaptive algorithm.

Let us point out that the abstract framework easily accommodates other applications, for example, systems of parabolic equations; see [36, Section 6.3] for the application to a phase-field tumor-growth system.

6.1 Allen–Cahn equation

We subject the (forced) Allen–Cahn equation, $\partial_{t}u-\Delta u+\epsilon^{-2}\psi^{\prime}(u)=f(t)$, to homogeneous Neumann boundary conditions. We choose the function spaces as $\mathcal{V}=H^{1}(\Omega)$, $\mathcal{V}^{*}=[H^{1}(\Omega)]^{*}$ and set $\mathcal{B}(u,v)=(\nabla u,\nabla v)$, $\mathcal{N}(u;v)=\frac{1}{\epsilon^{2}}\left(\psi^{\prime}(u),v\right)$ in Eq. 1, where $\epsilon$ is a parameter that controls the thickness of the diffuse interface (typical in phase-field models), and the nonlinear double-well function $\psi(u)$ is defined as (a standard truncated quartic polynomial)222Note that by choosing $\psi(u)=0$, one obtains the linear heat equation.

[TABLE]

Then, we obtain the weak form of the Allen–Cahn equation is: Find $u\in\mathcal{W}_{u^{0}}:=\left\{v\in L^{2}(0,T;\mathcal{V}),\partial_{t}v\in L^{2}\left(0,T;\mathcal{V}^{*}\right):v(0)=u^{0}\right\}$ such that $\forall v\in L^{2}(0,T;H^{1}(\Omega))$

[TABLE]

In our setting the IMEX scheme for Eq. 44 leads to the energy-stable time-stepping scheme introduced in [17]: find $u_{k+1}\in H^{1}(\Omega)$ such that $\forall v\in H^{1}(\Omega)$

[TABLE]

for $k=0,1,\ldots,N-1$, where the initial condition is $(u_{0},v)=(u^{0},v)$, $\forall v\in L^{2}(\Omega)$, and where we choose $\bar{f}_{k+1}=\frac{1}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}f(\cdot,t)\,\mathrm{d}t$. In particular, for a splitting of $\psi$ with a quadratic convex part, the resulting system is linear, for example:

[TABLE]

We then have the full discretization: find $u_{k+1}^{h}\in\mathcal{S}^{h,1}_{k+1}$ such that $\forall v^{h}\in\mathcal{S}^{h,1}_{k+1}$

[TABLE]

for $k=0,1,\ldots,N-1$, where the initial condition is $(u^{h}_{0},v^{h})=(u^{0},v^{h})$, $\forall v^{h}\in\mathcal{S}^{h,1}_{0}$.

According to the definition of $\mathcal{N}^{s}$ in Eq. 3, we can explicitly write the mean-value linearization of $\psi^{\prime}(u)$ in terms of $u$ and $\hat{u}$:

[TABLE]

although, because of its piecewise definition Eq. 43, this is an elaborate expression. For example, for $u,\hat{u}>1$ or $u,\hat{u}<-1$, we have $\psi^{\prime s}(u,\hat{u})=2$, and for $u,\hat{u}\in[-1,1]$, we have $\psi^{\prime s}(u,\hat{u})=u^{2}+\hat{u}^{2}+u\hat{u}-1$.

Then, by setting \mathcal{N}^{s}(u,\hat{u};w,z)=\frac{1}{\epsilon^{2}}\big{(}\psi^{\prime s}(u,\hat{u})z,w\big{)} in Eq. 5, the dual problem reads: find $z\in\mathcal{W}^{\bar{q}}:=\left\{v\in L^{2}(0,T;\mathcal{V}),\partial_{t}v\in L^{2}\left(0,T;\mathcal{V}^{*}\right):v(T)=\bar{q}\right\}$ such that $\forall w\in L^{2}(0,T;H^{1}(\Omega))$

[TABLE]

And the IMEX time-discrete dual problem, based on Eq. 23-Eq. 25, is defined by: find $z_{k}\in H^{1}(\Omega)$, $k=0,1,\ldots,N$, such that

[TABLE]

and for $k=1,2,\ldots,N-1$:

[TABLE]

where the terminal condition is

[TABLE]

with $q_{k}$, $k=0,1,\ldots,N$ defined in Lemma 3.4. Note that for our choice of $\psi_{c}$, the derivative $\psi^{\prime s}_{c}$ reduces to a constant. The main results of Section 3 hold, as shown in the following corollary.

Corollary 6.1** (Decomposed error representation for Allen–Cahn equation).**

The following error representation holds

[TABLE]

where the spatial error representation reduces to

[TABLE]

and the temporal error representation reduces to

[TABLE]

Proof 6.2**.**

*The result simply follows from Theorem 3.9 applied to Eq. 44–Eq. 50. *

6.2 Computable error indicator

Let $u_{\tau h}=\{u_{k}^{h/2}\}_{k=0}^{N}$ denote the solution of the fully-discrete primal problem Eq. 46 using time-step sizes $\{\tau_{k}\}_{k=1}^{N}$ and FE spaces $\{\mathcal{S}^{h,1}_{k}\}_{k=0}^{N}$, and let $u_{\tau,h/2}=\{u_{k}^{h/2}\}_{k=0}^{N}$ denote the solution of Eq. 46 using time-step sizes $\{\tau_{k}\}_{k=1}^{N}$ and enriched FE spaces $\{\mathcal{S}^{h/2,1}_{l}\}_{l=0,1/2,\ldots,N}$. Replacing $u_{\tau}$ with the computable $u_{\tau,h/2}$ in Eq. 48-Eq. 50, we obtain the full discretization of the dual problem using enriched FE spaces: find $z_{k}^{h/2}\in\mathcal{S}^{h/2,1}_{k}$, $k=0,1,\ldots,N$, such that

[TABLE]

and for $k=1,2,\ldots,N-1$:

[TABLE]

where the terminal condition is

[TABLE]

We denote by $z_{\tau,h/2}=\{z_{k}^{h/2}\}_{k=0}^{N}$ the solution of Eq. 51-Eq. 53. To get $z_{\tau/2,h/2}$, we first compute the space-time enriched approximation $u_{\tau/2,h/2}=\{\tilde{u}_{l}^{h/2}\}_{l=0,1/2,\ldots,N}$ of the primal problem using half time-step sizes $\{\tau_{1}/2,\tau_{1}/2,\tau_{2}/2,\tau_{2}/2,\ldots,\tau_{N}/2,\tau_{N}/2\}$ and enriched FE spaces $\{\mathcal{S}^{h/2,1}_{l}\}_{l=0,1/2,\ldots,N}$. Then, replacing $u_{\tau}$ with $u_{\tau/2,h/2}$ in Eq. 48-Eq. 50, we compute the space-time enriched approximation $z_{\tau/2,h/2}=\{\tilde{z}_{l}^{h/2}\}_{l=0,1/2,\ldots,N}$ using the same time-step sizes and FE spaces as $u_{\tau/2,h/2}$.

In the numerical examples for the Allen–Cahn equation in Section 6.3 we compute $u_{\tau,h/2}$ and $z_{\tau,h/2}$ directly by taking $u_{\tau,h/2}=u_{\tau/2,h/2}$ and $z_{\tau,h/2}=z_{\tau/2,h/2}$ at concurrent time steps as in Remark 4.1, and consider a piecewise-constant time-reconstruction $\hat{z}$ of $z_{\tau/2,h/2}$ for each time interval $[t_{k},t_{k+1})$, $k=0,1,\ldots,N-1$, i.e.

[TABLE]

According to Eq. 36 and Eq. 39-Eq. 42, we then get the following global space-time error estimate:

[TABLE]

and the local error indicators

[TABLE]

for $k=0,1,\ldots,N-1$ and $i=0,1,\ldots,M_{k}$. And since $z_{0}^{h/2}=\hat{z}(0)=\tilde{z}_{0}^{h/2}$,

[TABLE]

Remark 6.3**.**

*For the linear heat equation, $\partial_{t}u-\Delta u=f$, note that the dual problem and the computable error indicators are easily obtained from the above results by neglecting the nonlinear terms $\psi^{\prime}_{e}$ and $\psi^{\prime}_{c}$. *

6.3 Numerical results

In the following numerical experiments, we investigate the efficiency of the duality-based error estimates and the performance of the proposed adaptive algorithm. The results will be demonstrated in three parts. In the first part, we illustrate the consistency of the dual time scheme Eq. 48–Eq. 50 since we introduced a non-standard IMEX time-discrete dual problem Eq. 23–Eq. 25 which contains a nonstandard coefficient $\tau_{k+1}/\tau_{k}$. The second part is on the convergence of error estimate $\mathcal{E}_{\mathrm{st}}$ under uniform refinements for the nonlinear Allen–Cahn equation. In the third part, we apply the proposed duality-based adaptive algorithm to the linear heat equation and the nonlinear Allen–Cahn equation. We compare our adaptive result for the heat equation with the sequential-in-time adaptive algorithm (specifically for the heat equation) of Verfürth stated in [33, Sec. 6.8]. For the Allen–Cahn equation, we compare our adaptive results with uniform space-time refinements.

Consistency test

We begin with verifying numerically that the special IMEX time scheme of the dual problem Eq. 23-Eq. 25 is first-order accurate in time with respect to refinements of uniform initial time steps (i.e. $\tau_{k+1}/\tau_{k}=1$) and, in particular, nonuniform initial time steps (i.e. $\tau_{k+1}/\tau_{k}\neq 1$). Here, the Allen–Cahn equation is considered in 1D on the domain $\Omega=(-3,3)$ with parameter $\epsilon=1$. The spatial mesh is composed of $256$ elements along the axis. We consider a manufactured solution which oscillates in time:

[TABLE]

The convergence results are presented on a double logarithmic scale in Figure 3. Figure 3a is the convergence of the error in $L^{2}$ norm at time $T$ based on time-step refinements using uniform initial time steps $\{\tau_{k}\}_{k=1}^{5}=\{0.1,0.1,0.1,0.1,0.1\}$, and Figure 3b uses nonuniform initial time steps $\{\tau_{k}\}_{k=1}^{5}=\{2\cdot 10^{-2}(3k+1)\}_{k=4}^{0}=\{0.26,0.2,0.14,0.08,0.02\}$. For both uniform and nonuniform initial time steps, the observed rates are close to 1, which demonstrates that Eq. 23-Eq. 25 is a first-order time-accurate scheme.

Effectivity test

In this numerical experiment, we consider the Allen–Cahn equation with an exact solution

[TABLE]

on the domain $\Omega=(0,1)\times(0,1)$ where $\epsilon=1$. We suppose that we are interested in the error at final time $T$, i.e., in (2) we take $q=0$ and $\bar{q}=u(T)-Iu_{\tau h}(T)$ (which is approximated in the computations by $\bar{q}\approx\tilde{u}_{N}^{h/2}-u_{N}^{h}$, see Section 4). The error of interest is thus $\mathcal{Q}(u)-\mathcal{Q}(Iu_{\tau h})=\|u(T)-Iu_{\tau h}(T)\|_{L^{2}(\Omega)}^{2}=\|u(T)-u_{N}^{h}\|_{L^{2}(\Omega)}^{2}$.

To investigate the error estimate with respect to the spatial discretization, temporal discretization and space-time discretization, we compute $\mathcal{E}_{\mathrm{st}}$ according to Eq. 54 under uniform spatial refinement, uniform temporal refinement and uniform space-time refinement, respectively. In view of the nonstandard time discretization scheme for the dual problem, we also investigate the accuracy of the error estimate $\mathcal{E}_{\mathrm{st}}$ with respect to uniform and nonuniform initial time step sizes. Table 1 presents the convergence of $\mathcal{E}_{\mathrm{st}}$ under uniform spatial refinement for a sufficiently small time step size ($\tau_{k}=1e-4$). For uniform initial time steps $\{\tau_{k}\}_{k=1}^{4}=\{0.05,0.05,0.05,0.05\}$, the left side of Table 2 shows the convergence of $\mathcal{E}_{\mathrm{st}}$ under uniform temporal refinement for a sufficiently fine spatial mesh ($128\times 128$ elements) and the left side of Table 3 shows the convergence under uniform space-time refinements. For nonuniform initial time steps $\{\tau_{k}\}_{k=1}^{4}=\{0.08,0.06,0.04,0.02\}$, the convergence result for uniform temporal refinement with a sufficiently fine spatial mesh is presented on the right side of Table 2, and the convergence result for uniform space-time refinements is presented on the right side of Table 3. The effectivity results for different refinements are also presented in the two tables.

From all the results of Table 1–Table 3, we observe that the error and the estimate converge with the same order under various refinements. The effectivity indices are always between $0.4$ and $1$, and seem to converge to a constant. This indicates the asymptotic effectivity of the error estimate.

Adaptivity test for the linear heat equation: Clockwise moving source

Let $\Omega=(0,1)\times(0,2)$. We take the example from Asner, Tavener and Kay [2] by choosing the right hand side and the initial and the boundary condition so that the exact solution is given by

[TABLE]

which is an exponential peak moving clockwise inside the domain, see Figure 4. We set up the final time $T=0.5$ and $\theta=\lambda=0.8$, and start the adaptive procedure with a uniform mesh containing $4\times 8$ elements and $10$ equally distributed time steps (i.e. $\tau=0.05$). We aim to minimize the error at final time $T$, i.e., as before, $\mathcal{Q}(u)-\mathcal{Q}(Iu_{\tau h})=\|u(T)-Iu_{\tau h}(T)\|_{L^{2}(\Omega)}^{2}=\|u(T)-u_{N}^{h}\|_{L^{2}(\Omega)}^{2}$. For the sequential adaptive algorithm, we choose an initial coarse mesh with $4\times 8$ elements and an initial $\tau=0.01$.

In Figure 4, we show a comparison of the results obtained by the sequential adaptive algorithm and our duality-based adaptive algorithm. The first row corresponds to the exact solution. The second and third row are the snapshots of adaptively refined meshes with corresponding approximations by using the two adaptive algorithms. Figure 6 illustrates the various time steps over time for the two adaptive algorithms. A comparison of the convergence of the error is shown in Figure 6 where ‘Total dof’ refers to total number of degrees of freedom, $\sum_{k=0}^{N}M_{k}$.

From these results, it can be seen that the spatial and temporal refinements of the sequential adaptive algorithm focus on tracing the movement of the Gaussian shaped peak of the solution in the space-time domain. In contrast, since our duality-based adaptive algorithm targets the $L^{2}$ norm of the final error, the refinements of our algorithm are only concentrated at the final moment of the space-time domain. Away from this final moment, where residuals contribute much less to the final error, the original mesh and time step already provide sufficient resolution and need not be refined. This leads to fewer degrees of freedom and time steps for reaching the same accuracy of the final solution than for the sequential adaptive algorithm.

In Figure 6, we also show the convergence of the error and the error estimate $\mathcal{E}_{\mathrm{st}}$ for the duality-based adaptive algorithm, which demonstrates the accuracy of the error estimate under adaptive refinement.

Adaptivity test for the Allen–Cahn equation: Shrinking ring

The Allen–Cahn dynamics of this test case is a shrinking ring (with diffuse interfaces) in the middle of the domain $\Omega=(-1,1)^{2}$. The initial condition is set as

[TABLE]

where $\epsilon=0.0625$, see Figure 7. The inner circle has a small radius of $0.15$ which is expected to vanish much earlier than the outer circle. We are interested in the final error (i.e. $\mathcal{Q}(u)-\mathcal{Q}(Iu_{\tau h})=\|u-Iu_{\tau h}\|_{L^{2}}^{2}(T)$) when the inner circle has disappeared, i.e., the final time $T=0.02$. To have a reference value for the error, we compute an approximation to Eq. 46 on a uniform mesh with $512^{2}$ elements and a uniform time step size $\tau=1e-5$. For the adaptive algorithm, we take a coarse initial time-step size $\tau_{k}=5e-3$ and a coarse initial mesh with $16^{2}$ elements. The fractions in the adaptive algorithm are selected as $\theta=\lambda=0.8$.

Snapshots of the results are presented in Figure 7. The first row shows snapshots of the reference solution, while the second and fourth row show the primal approximation and the corresponding adaptive mesh obtained by the proposed adaptive algorithm. The computed dual solution $z$ is displayed in the third row. It can be observed that the dual solution grows as time progresses. This growth is localized at the interface of the outer ring. Let us mention that, for better visualization, in the plots of the dual solution the range of the color bars is adapted to each plot. Figure 9 displays the various time steps over time for the shrinking ring. As expected, smaller time steps are need at the vanishing moment of the inner ring and towards the final time $T$. Furthermore, the interfaces in the solution are well-resolved throughout the simulation time, with a significant increase of resolution towards the final time $T$.

The convergence of the error estimate and the error for the duality-based adaptive algorithm in comparison to uniform space-time refinement is shown in Figure 9 where ‘Total dof’ refer to total number of degrees of freedom $\sum_{k=0}^{N}M_{k}$. Note that the error exhibits a plateau for the most refined approximations because the accuracy of the adaptively-refined approximations surpasses that of the reference approximation.

7 Conclusion

In this work we carried out a comprehensive study of duality-based a posteriori error estimates for semi-linear parabolic problems, with a special focus on discretizations using the finite element method in space combined with IMEX time stepping. We introduced a decomposition of the error estimates to identify the separate error contributions due to temporal and spatial approximation. The key idea is to adapt the residual decomposition by Verfürth to our duality-based error representation and propose a specially-tailored time-discrete dual problem. The resultant error indicators quantify the spatial and temporal discretization errors and provide information to drive adaptive mesh refinement and adaptive time-step selection.

To illustrate the performance of the duality-based error estimates and the proposed adaptive algorithm, we presented numerical experiments for the heat equation and Allen–Cahn equation. We refer to [36, Section 6.3] for the application to systems. The numerical results verified the accuracy and the effectivity of the error estimate in test problems. We also observed the overall good quality of the adaptive algorithm.

The proposed methodology can be further extended to other finite difference time-stepping schemes which do not fit in our considered abstract setting, e.g. higher-order multi-stage Runge-Kutta schemes. The key challenge in any such extension is the derivation of a specially-tailored time-discrete dual problem. Our analysis indicates that these dual problems can be derived systematically by means of “backwards” summation-by-parts on the time-discrete system.