2.9 Time Step and Mesh Control

A time step and mesh should always be chosen in a such a way that the error made during the numerical solving procedure is acceptably small. In this chapter we consider a general system of nonlinear PDEs and there is no rigorous theory available which covers the error behavior of such systems. However for certain special classes of problems there are well-known results which can be used as basis for the study of more complex systems.

For our purpose the case $\mathcal{L}(\mathbf{c})=\triangle c$ and $\mathbf{f}=0$ in (2.1) is of special interest, since it represents the problems discussed in this work in its simplest form.

Using a backward Euler scheme and a finite element method with linear basic nodal functions we have the estimate [13,10,9],

$$\Vert c_h^n-c(t_n)\Vert _{L_2(\Omega)}\leq \gamma \Bigl(1+ \frac{t_n}{\Delta t_n}\Bigr)^{1/2} \Bigl(\underset{m\leq n}{\text... ...a)}dt+h^2\underset{t\leq t_n}{\text{max}} \Vert c(t)\Vert _{H^2(\Omega)} \Bigr)$$ (2.51)

for $I_n=[t_{n-1},t_{n}]$ and $t_{n}=t_{n-1}+\Delta t_n$, where $\Vert f\Vert _{L_2(\Omega)} = \sqrt{(f,f)}$, and

$$\Vert f\Vert _{H^2(\Omega)}= \Bigl(\int_{\Omega}f^2 +\sum_{i}\Big... ..._{i,j}\Bigl(\frac{\partial^2 f}{\partial x_i x_j}\Bigr)^2 d\Omega\Bigr)^{1/2},$$ (2.52)

$c^n_h$ is the finite element approximation at the discrete time step $t_n$, and $h$ is the the maximal edge length of the discretization. The estimate (2.51) is valid independently of the dimension of a given problem [13,10,9].

Taking into account that the logarithmic quantity in ( $\ref{eq:time-step-and-mesh-control:1}$) changes slowly, it can be absorbed in the constant $\gamma$ and considering that,

$$\int_{I_n}\Bigl\Vert\frac{\partial c}{\partial t}\Bigr\Vert _{L_2... ...}{\text{sup}} \Bigl\Vert\frac{\partial c}{\partial t}\Bigr\Vert _{L_2(\Omega)},$$ (2.53)

we can write (2.51) as,

$$\underset{t\leq t_N}{\text{max}}\Vert c(t)-C(t)\Vert\leq \gamma \... ...Bigr)+h^2\underset{t\leq t_N}{\text{max}} \Vert c(t)\Vert _{H^2(\Omega)}\Bigr).$$ (2.54)

The first term on the right side of (2.54) corresponds to the time discretization error and the second to the spatial discretization error [13].

If we want to achieve the time discretization error to be bounded by $\delta$ for $t\leq t_N$, we have to ensure that,

$$\Delta t_n \underset{t\in I_n}{\text{sup}} \Bigl\Vert\frac{\partial c}{\partial t} \Bigl\Vert \leq \frac{\delta}{\gamma},\quad n=1,\dots,N.$$ (2.55)

Making the plausible assumption that,

$$\Delta t_n \underset{t\in I_n}{\text{sup}} \Bigl\Vert\frac{\partial c}{\partial t} \Bigl\Vert\thicksim \Vert c_h^n-c_h^{n-1}\Vert,$$ (2.56)

we obtain a rule for choosing the local time step $\Delta t_n$. Starting with $c_h^{n-1}$ we calculate $c_h^{n}$ using the time step $\Delta t_n= \Delta t_{n-1}$. If $\Vert c_h^n-c_h^{n-1}\Vert\leq \delta/\gamma$, the value $c_h^n$ is accepted and otherwise $\Delta t_n$ is decreased. Variations of this approach are also possible. For example, the initial $\Delta t_n$ can be chosen as [13],

$$\Delta t_n = \frac{\delta \Delta t_{n-1}}{\gamma\Vert c_h^{n-1}-c_h^{n-2}\Vert}.$$ (2.57)

To estimate the error of the spatial finite element discretization several methods are known today. These methods are based on the estimation of the second term on the right side of (2.54) using the solution obtained by means of the finite element method (a posteriori error estimators)[13,10,9].

First we will present the flux error estimator based on the work by Babuska, Bieterman, and Rheinboldt [9].

We consider a tethraedron $T_0\in T_h(\Omega)$ and define a tehraedron-patch $P(T_0)$ as the subset of $T_h(\Omega)$ of all tethraedrons which have one common face with $T_0$. Each point $p^i$ of the discretization holds a concentration value $C_i\geq 0$ and inside of the each element the concentration is linearly interpolated as $C(\mathbf{x},t_{n})=\sum_{i=1}^4 C_i^n\varphi_i(\mathbf{x})$. In the following we omit the discrete time index $n$, because we consider only a single time step.

The concentration gradient, calculated by $\nabla C(\mathbf{x})=\sum_{i=1}^4 C_i\nabla\varphi_i(\mathbf{x})$, is constant on each element $T$ and discontinuously changes between two neighboring elements. We write $\nabla C(\mathbf{x}) = \nabla_T C$ for $\mathbf{x}\in T$. This discontinuity expresses the error of the approximation and can be quantified on the tehraedron-patch $P(T_0)$ as the flux error estimator,

$$e_{F}^2(T_0)= h^2\sum_{T\in P(T_0)}\int_{T\cap T_0}\Bigl[\frac{\partial C}{\partial\mathbf{n}} \Bigr]^2_{T,T_0} df,$$ (2.58)

$$\Bigl[\frac{\partial C}{\partial\mathbf{n}} \Bigr]_{T,T_0}=\mathbf{n}\cdot(\nabla_T C - \nabla_{T_0} C),$$ (2.59)

where $\mathbf{n}$ is the outside normal on the face $f=T\cap T_0$ and $h$ is the length of the longest edge of the tethraedron $T$.

The second spatial discretization error estimator presented here is the gradient recovery error estimator which is widely used today in the finite element community [14].

For the derivation of the gradient recovery error estimator we first need the recovered gradient. The recovered gradient is calculated as the Lagrange interpolate at the nodes of the mesh. Taking a node $a\in T_h(\Omega)$ and defining $P(a)$ as a point-patch of all elements which contain node $a$ as one of its nodes. The recovered gradient $\nabla^RC(a)$ at the node $a$, is then calculated as,

$$\nabla^RC(a) = \sum_{T\in P(a)} \theta_{a,T} \nabla_T C.$$ (2.60)

We chose the weights $\theta_{a,T}$ of the Lagrange interpolate as,

$$\theta_{a,T}=\frac{\vert T\vert}{\sum_{T'\in P(a)}\vert T'\vert}.$$ (2.61)

The values of the concentration gradient on the element nodes $p^j, (j=1,2,3,4)$ enable to define a linear approximation of the gradient inside the element as,

$$\nabla^RC(\mathbf{x})=\sum_{i=1}^4 \nabla^RC(p^i)\varphi_i(\mathbf{x}).$$ (2.62)

The deviation of $\nabla_T C$ from $\nabla^RC(\mathbf{x})$ on the element $T$ gives the gradient recovery error estimator,

$$e_{R}^2(T)=\int_T \vert\nabla^RC(\mathbf{x})-\nabla_T C\vert^2 d\Omega.$$ (2.63)

Although the gradient recovery error estimator and the flux error estimator are apparantly crude and derived on the basis of a simple parabolic model problem, they result in a suprisingly good error estimate [14]. These error estimators are used as indicators for the local mesh adaptation in several FEM based software tools [15,16,17].

For the applications discussed in this thesis the physically motivated dosis error estimator (Section 3.9) has even more importance than the error estimators presented above.

Independent of the actual mesh adaptation strategy and error estimator, the following algorithm for the finite element adaptive computation can be applied [15],

Here $E$ denotes the total error, $e(T)$ is an error estimator, $X=L_2(\Omega),H^1(\Omega),$ or $H^2(\Omega)$, and $\eta$ is a given tolerance for the relative error.