Local and parallel finite element algorithms based on two-grid discretizations
By Jinchao Xu and Aihui Zhou
Abstract
A number of new local and parallel discretization and adaptive finite element algorithms are proposed and analyzed in this paper for elliptic boundary value problems. These algorithms are motivated by the observation that, for a solution to some elliptic problems, low frequency components can be approximated well by a relatively coarse grid and high frequency components can be computed on a fine grid by some local and parallel procedure. The theoretical tools for analyzing these methods are some local a priori and a posteriori estimates that are also obtained in this paper for finite element solutions on general shape-regular grids. Some numerical experiments are also presented to support the theory.
1. Introduction
In this paper, we will propose some new parallel techniques for finite element computation. These techniques are based on our understanding of the local and global properties of a finite element solution to some elliptic problems. Simply speaking, the global behavior of a solution is mostly governed by low frequency components while the local behavior is mostly governed by high frequency components. The main idea of our new algorithms is to use a coarse grid to approximate the low frequencies and then to use a fine grid to correct the resulted residue (which contains mostly high frequencies) by some local/parallel procedures.
Let us now give a somewhat more detailed but informal (and hopefully informative) description of the main ideas and results in this paper. We consider the following very simple model problem posed on a convex polygonal domain $\Omega \subset R^2$:
The main philosophy behind this paper is that we should treat phenomena of different scales by different tools. In multigrid and domain decomposition methods, this kind of idea is used to devise iterative methods for solving a given discretization scheme (see e.g. Bank Reference 10, Bramble Reference 19, Chan and Mathew Reference 22, Hackbusch Reference 31, Xu Reference 49 and Yserentant Reference 53); while in our approach, we try to use this type of idea for designing discretization schemes. The two-grid method proposed by the first author Reference 48Reference 50Reference 51 (and later further investigated by many others such as Reference 3Reference 16Reference 24Reference 25Reference 33Reference 34Reference 42) is a result of such a consideration. The two-grid method is based on the observation that, for an equation like Equation 1.1, the symmetric positive definite leading term $-\Delta u$ dominates the equation on high frequencies, while the low frequencies, as in a multigrid method, can be well approximated by a relatively coarse grid. Therefore, by first approximating the equation on a coarse grid, say $T^H(\Omega )$, we can then correct the residue (in which high frequencies dominate) on a finer grid, say $T^h(\Omega )$, by ignoring the lower order term and solving the resulting symmetric positive definite system.
For elliptic problems, the low frequencies are more global, while the high frequencies are more local. This fact is crucial in the multigrid methodology, in which high frequency errors are damped out by local relaxation techniques while low frequencies are handled by coarse grids. If we consider this fact more carefully, we can then imagine that if we first approximate the equation Equation 1.1 on a coarse grid $T^H(\Omega )$, the residue which is dominated by high frequencies can then be resolved locally. This is precisely the central idea of the new algorithms in this paper, and is based on the local behavior of finite element approximations presented in Section 3.
One technical tool for motivating this idea is the local error estimate for finite element approximations. Let $u_h$ be a finite element approximation to Equation 1.1 on a quasi-uniform grid $T^h(\Omega )$. Then the following kind of local error estimate holds (see Theorem 3.4):
where $S^h_0(\Omega )$ is the finite element space associated with $T^h(\Omega )$ and $D\subset \subset \Omega _0\subset \subset \Omega$ (here $D\subset \subset \Omega _0$ means that $\mathrm{dist} (\partial D\setminus \partial \Omega , \partial \Omega _0\setminus \partial \Omega )>0$).
At first glance, this type of estimate does not appear to be clearly related to what we said above, but we shall soon explain the connection. The above kind of estimate is available in the literature for quasi-uniform grids (see Nitsche and Schatz Reference 35, Schatz Reference 38, Schatz and Wahlbin Reference 39Reference 40 and Wahlbin Reference 46Reference 47), but we need them on locally refined grids with different mesh scales, which were also discussed in Reference 54 for local quasi-uniform grids. We are indeed able to extend these estimates to very general grids, which is one technical part of this paper—see Section 3.
Now we consider a very special grid that is obtained by refining a given coarse grid $T^H(\Omega )$ for the region $\Omega _0$, and obtain a locally refined grid $T^h_H(\Omega )$ with mesh size $h$ in $\Omega _0$ and size $H$ away from $\Omega _0$ (see Figure 1). Let us for example consider linear finite element discretizations on this grid. Then, by Equation 1.2 and some well-known finite element error estimates, we obtain (see Corollary 3.5)
This estimate means that we can obtain an asymptotically optimal error in the $H^1$ norm locally by taking $H=O(\sqrt h)$.
With this basic result, it is then not difficult to devise a parallel algorithm on a fine grid by using a collection of overlapped subdomains. See Section 4.
In the above approach, all the “local” solvers need to be coupled with the global coarse grid $T^H(\Omega )$ in some way. We can actually improve the above procedure by using a residue-correction technique used in Xu Reference 48Reference 50Reference 51. One prototype algorithm is as follows (see Section 4):
(1)
Solve on a coarse grid: Find $u_H\in S^H_0(\Omega )$ such that$$\begin{eqnarray*} \int _{\Omega }\nabla u_H\cdot \nabla v+\mathbf{b}\cdot \nabla u_H v=\int _{\Omega }fv, \quad \forall v\in S^H_0(\Omega ). \end{eqnarray*}$$
(2)
Correct the residue (with SPD part only) on a fine grid: Find $e_h\in S^h_0(\Omega _0)$ such that$$\begin{eqnarray*} \int _{\Omega }\nabla e_h\cdot \nabla v=\int _{\Omega }fv -\int _\Omega \nabla u_H\nabla v -\int _{\Omega }\mathbf{ b}\cdot \nabla u_H v, \quad \forall v\in S^h_0(\Omega _0). \end{eqnarray*}$$
In this algorithm, a coarse grid problem only needs to be solved once and it does not have to be coupled with the subsequence of parallel local solvers. For the above algorithm, we can still establish the following result (see Theorem 4.3):
This is a very satisfying result in many ways. As a consequence, for example, we can then design the following type of parallel algorithms: first solve the problem on a coarse grid, and then correct the residue in parallel on a collection of overlapped subdomains on a fine grid.
In practical finite element computations, it is desirable to carry out the finite element computations in an adaptive fashion, cf. Ainsworth and Oden Reference 2, Babuška, Duran and Rodriguez Reference 4, Babuška and Rheinboldt Reference 5, Babuška, Zienkiewick, Gago and Oliveira Reference 6, Bank and Weiser Reference 15, Johnson Reference 32 and Verfürth Reference 45. A typical procedure is first to start with a coarse grid and then use some a posteriori estimates as a guidance to properly refine the mesh to achieve the desired accuracy. In the existing literature, a posteriori error estimates are often obtained globally, but in practical applications, they are often used locally (see e.g. Reference 2Reference 4Reference 5Reference 14Reference 26Reference 27Reference 28Reference 36Reference 43Reference 44Reference 45 and references cited therein). In this paper, we shall also present some local a posteriori error estimates that would give a certain justification of the local application of the a posteriori estimates. For example, we shall prove the following type of local a posteriori estimate on a general grid (see Theorem 3.9):
where $E_h(\Omega _0)$ is the usual a posteriori estimator on the domain $\Omega _0$. Again we notice that the global term $\|u-u_h\|_{0,\Omega }$ in the above estimate is a high order term.
More important, based on a posteriori estimates like Equation 1.5 and a priori estimates like Equation 1.2, we are able to design an adaptive procedure that can be carried locally in a given subdomain and hence in parallel. We believe this type of parallel adaptive techniques will have great implications in practical parallel computations.
The rest of the paper is organized as follows. In Section 2, some preliminary materials are provided. In Section 3, a number of local a priori and a posteriori error estimates are obtained for finite element discretizations on general shape regular grids. Based upon these local error estimates, several new local/parallel algorithms are devised and analyzed in Section 4, and local and parallel adaptive processes are discussed in Section 5. In Section 6, some numerical experiments, which support our theory, are reported. Finally in Section 7, some further remarks are presented.
2. Preliminaries
In this section, we shall describe some basic notation and basic assumptions on the finite element spaces, and then study properties of the finite element approximation to a general linear second order elliptic boundary value problem.
Let $\Omega$ be a bounded domain in $R^d\ (d\ge 1)$. We shall use the standard notation for Sobolev spaces $W^{s,p}(\Omega )$ and their associated norms and seminorms, see e.g. Reference 1Reference 23. For $p=2$, we denote $H^s(\Omega )=W^{s,2}(\Omega )$ and $H^1_0(\Omega )=\{v\in H^1(\Omega ): v\mid _{\partial \Omega }=0\}$, where $v\mid _{\partial \Omega }=0$ is in the sense of trace, $\|\cdot \|_{s,\Omega }= \|\cdot \|_{s,2,\Omega }$ and $\|\cdot \|_{\Omega }=\|\cdot \|_{0,2,\Omega }$. (In some places of this paper, $\|\cdot \|_{s,2,\Omega }$ should be viewed as piecewise defined, if necessary.) The space $H^{-1}(\Omega )$, the dual of $H^1_0(\Omega )$, will also be used.
For $D\subset G\subset \Omega$, we write $D\subset \subset G$ to mean that $\mathrm{dist}( \partial D\setminus \partial \Omega , \partial G\setminus \partial \Omega )>0$, see Figure 2.
Throughout this paper, we shall use the letter C (with or without subscripts) to denote a generic positive constant which may stand for different values at its different occurrences. For convenience, following Reference 49, the symbols $\lesssim , \gtrsim$ and $\cong$ will be used in this paper. Specifically, $x_1\lesssim y_1, x_2\gtrsim y_2$ and $x_3\cong y_3,$ mean that $x_1\le C_1 y_1,\quad x_2\ge c_2 y_2$ and $c_3 x_3\le y_3 \le C_3x_3$ for some constants $C_1, c_2, c_3$ and $C_3$ that are independent of the mesh parameters.
Note that any $w\in H^1_0(\Omega _0)$ can be naturally extended to be a function in $H^1_0(\Omega )$ with zero outside of $\Omega _0$. We shall state this fact by the slightly abused notation $H^1_0(\Omega _0)\subset H^1_0(\Omega )$.
2.1. Finite element spaces
For generality, we will not concentrate on any specific finite element space, rather we shall study a class of finite element spaces that satisfy certain assumptions. We now describe those assumptions.
Assume that $T^h(\Omega )=\{\tau \}$ is a mesh of $\Omega$ with mesh-size function $h(x)$ whose value is the diameter $h_{\tau }$ of the element $\tau$ containing $x$. One basic assumption on the mesh is that it is not exceedingly over-refined locally, namely
where $h_{_\Omega }=\max _{x\in \Omega }h(x)$ is the (largest) mesh size of $T^h(\Omega )$.
This is apparently a very mild assumption, and most practical meshes should satisfy it. Sometimes, we will drop the subscript in $h_{_\Omega }$, writing $h$ for the mesh size on a domain that is clear from the context.
Associated with a mesh $T^h(\Omega )$, let $S^h(\Omega )\subset H^1(\Omega )$ be a finite dimensional subspace on $\Omega$ and $S^h_0(\Omega ) =S^h(\Omega )\cap H^1_0(\Omega )$. Given $G\subset \Omega$, we define $S^h(G)$ and $T^h(G)$ to be the restriction of $S^h(\Omega )$ and $T^h(\Omega )$ to $G$, and
A.3.Superapproximation. For $G\subset \Omega _0$, let $\omega \in C^{\infty }_0(\Omega )$ with supp $\omega \subset \subset G$. Then for any $w\in S^h(G)$, there is $v\in S^0_h(G)$ such that
Next we shall give an example of finite element spaces. The assumptions mentioned above are satisfied by most of the finite element spaces used in practice. We shall now briefly describe a special but popular finite element space for illustration. For simplicity, let us assume that $\Omega$ is a polygonal domain. Let $T^h(\Omega )$ consist of shape-regular simplices and define $S^h(\Omega )$ to be a space of continuous functions on $\Omega$ such that for $v\in S^h(\Omega ), v$ restricted to each $\tau$ is a polynomial of total degree $\le r$, namely
where $P^r_{\tau }$ is the space of polynomials of degree not greater than a positive integer $r$. These are the Lagrange finite element spaces, and they satisfy all the above assumptions.
The approximation assumptions A.1 and A.1$'$ are well-known for the Lagrange finite element spaces. A.2 and A.4 are well-known and they can be easily proven by a standard scaling argument. The superapproximation assumption has been discussed in many papers, cf. Reference 35Reference 39Reference 40Reference 46Reference 47. This superapproximation assumption can be easily verified for the Lagrange finite element spaces Equation 2.8, using a locally defined interpolation operator $I_h$ satisfying
where $|\cdot |_{r, \infty ,\Omega }$ denotes the $r^{th}$ semi-norm involving only $r^{th}$ derivatives. Setting $\phi =\omega w$ in Equation 2.9 and noting that $r^{th}$ derivatives of $w$ vanish, and using the inverse estimates, one obtains Assumption A.3.
The verification of A.5 can go as follows. For $v\in S^h(G)$, let $\chi \in S^h_0(G)$ be the unique function satisfying
Then $v-\chi$ is discrete harmonic. The desired result then follows from the following well-known (cf. Reference 52) estimate for discrete harmonic functions:
In this subsection, we shall study some basic properties of general second order elliptic boundary value problems and their finite element approximations. We consider the homogeneous boundary value problem
$$\begin{equation} \left\{\begin{array}{rl} L u &= f,\quad \text{in} \quad \Omega ,\\u &= 0,\quad \text{on}\nobreakspace\nobreakspace \partial \Omega .\end{array}\right. \cssId{problem}{\tag{2.10}} \end{equation}$$
Here $L$ is a general linear second order elliptic operator:
Our basic assumption is that (Equation 2.11) is well-posed, namely (Equation 2.11) is uniquely solvable for any $f\in H^{-1}(\Omega )$. (A simple sufficient condition for this assumption to be satisfied is that $c\ge 0$.) An application of the open-mapping theorem yields
Throughout this paper, we will assume that $h_{_\Omega }\ll 1$ (when $N(\cdot ,\cdot )\not =0$) and Assumption A.1 holds, so that the above two estimates hold. From the above two estimates, we can then define Galerkin-projections $P_h\ (\equiv P_h^{\Omega }): H^1_0(\Omega )\mapsto S^h_0(\Omega )$ and $P^*_h\ (\equiv (P^{\Omega }_h)^*): H^1_0(\Omega ) \mapsto S^h_0(\Omega )$ by
From (Equation 2.17), various a global priori error estimates can be obtained from the approximate properties of the finite element subspaces $S^h(\Omega )$ (cf. Reference 23). Particularly, by Assumption A.1, if $u\in H^1_0(\Omega )$, then
Similarly, if Assumption R(G) holds, we can define $\rho _{_G}(h)$ well.
3. Local a priori and a posteriori error estimates
In this section, we shall present a number of local a priori and a posteriori error estimates for finite element discretizations on general shape regular grids. Novelties of our estimates lie in, for example, the weak assumption on the underlying grids as well as the generality of model continuous problems. Although these general estimates are theoretically interesting in their own right, our main motivation is to use them to devise and analyze some new local/parallel methods to be presented in the following sections.
We shall now present a local a priori estimate for finite element approximation that will play a crucial role in our analysis. This type of estimates can be found in Reference 35Reference 39Reference 40Reference 46Reference 47; a new feature here is the generality of the underlying finite element grid for which this estimate is proven valid.
The local stability of the Galerkin-projection is stated as follows.
3.2. Local a posteriori error estimates
In this section, we shall present local a posteriori error estimate in energy-norm. First, we need the following technical result.
To state the a posteriori estimates, we need some more notation. Let $\partial T^h (\Omega )$ be the set of all the interior faces of the mesh $T^h(\Omega )$, and $\partial T^h(\tau )=\{F\in \partial T^h(\Omega ): F\subset {\bar{\tau }}\}.$ For each $F\in \partial T^h(\Omega )$, let $\mathbf{n}_{F}$ be a unit vector normal to $F$, and define for $v\in S^h(\Omega )$
One sees that $\eta (u_h), h\eta (u_h)$ and $h^2\eta (u_h)$ are computable in terms of the finite element solution $u_h(\equiv P_hu)$.
Remark 3.11.
The argument here easily extends to other boundary conditions, provided the problems are well-posed.
4. New local and parallel algorithms
In this section we shall present some new local and parallel finite element algorithms. These algorithms are motivated by the local error estimates studied in the previous section. We shall first discuss the local algorithms. The generalization of the local algorithms to parallel algorithms is straightforward.
For clarity, let $\Omega$ be a polygonal domain and $S_0^h(\Omega )\subset H^1_0(\Omega )$ be a finite element subspace of degree $r$ associated with a grid $T^h(\Omega )$. Let $u_h\in S^h_0(\Omega )$ be the solution of the standard finite element scheme for solving (Equation 2.11):
Either locally or globally, with proper regularity assumption, we have the following error estimate:
$$\begin{equation*} \|u-u_h\|_{1}\lesssim h^s, \quad 1\le s \le r. \end{equation*}$$
With this type of error estimates in mind, in the rest of this section, we will only compare the approximate solutions from our new algorithms with $u_h$ instead of the exact solution $u$.
4.1. Local algorithms
The local algorithms we shall now present can be used to obtain approximate solution on a given subdomain, mostly by local computation. The main idea is that the more global component of a finite element solution may be obtained by a relatively coarser grid, and the rest of the computation can then be localized.
Roughly speaking, our new algorithms will be based on sometimes one coarse grid of size $H$ and one fine grid of size $h\ll H$, and sometimes on a grid that is fine in a subdomain and coarse on the rest of the domain. The fine grid may be only defined locally. In our analysis, we shall use an auxiliary fine grid, say $T^h(\Omega )$, that is globally defined. One basic assumption for this auxiliary fine grid is that it should coincide with the local fine grid in the subdomain of interest.
Let $T^H(\Omega )$ be a shape-regular coarse grid, of size $H\gg h$, so that the highly locally refined mesh $T^h(\Omega _0)$ can be obtained, where $\Omega _0$ is a slightly larger subdomain containing a subdomain $D\subset \Omega$ (namely $D\subset \subset \Omega _0$), see Figure 3. More precisely, we let $T_H^h(\Omega )$ denote a locally refined shape-regular mesh that may be viewed as being obtained by refining $T^H(\Omega )$ locally around the subdomain $D$ in such a way that $T^h_H(\Omega _0)=T^h(\Omega _0)$. We are interested in obtaining the approximation solution in the given subdomain $D$ with an accuracy comparable to that from $T^h(\Omega )$. We shall propose two different gridding strategies for obtaining finite element approximations on the subdomain $D$ (see Figure 4). We denote the corresponding finite element space by $S_0^{H,h}(\Omega )\subset H_0^1(\Omega )$, which consists of piecewise polynomial of degree less than or equal to $r$ in this section.
4.1.1. First approach
The first strategy is simply to solve a standard finite element solution in $S_0^{H,h}(\Omega )$.
Algorithm A0.Find $u_H^h\in S_0^{H,h}(\Omega )$ such that
Strictly speaking, this algorithm is still a global algorithm as a global problem is solved. But it is designed to obtain a local approximation in the subdomain $D$ and it makes use of a mesh that is much coarser away from $D$.
Theorem 4.1.
Assume that $u_H^h\in S_0^{H,h}(\Omega )$ is obtained by Algorithm A0. Then
We would like to remark here that similar locally refined grids have been used for different purposes in the literature, cf. Bramble Reference 19, Bramble, Ewing, Parashkevov and Pasciak Reference 20, and Bramble, Ewing, Pasciak and Schatz Reference 21.
4.1.2. Second approach
Our second strategy is in a way an improvement of the first strategy. In this strategy, we first solve a global problem only on the given coarse grid $T^H(\Omega )$, and we then correct the residue locally on the fine mesh $T^h(\Omega _0)\ (=T^h_H(\Omega _0))$. Let $S^H_0(\Omega )\subset H^1_0(\Omega )$ be the finite element space of degree $r$ defined on $T^H(\Omega )$.
A prototype of our new local algorithms is as follows.
Algorithm B0. 1. Find a global coarse grid solution $u_H\in S_0^H(\Omega )$:
To estimate $\|e_h\|_{0,\Omega _0}$, we use the Aubin-Nitsche duality argument. Given any $\phi \in L^2(\Omega _0)$, there exists $w\in H^1_0(\Omega _0)$ such that
Following the basic idea in Xu Reference 48Reference 50Reference 51, for non-SPD problems Algorithm B0 may be modified in such a way that the local fine grid correction in Algorithm B0 is only carried out for the symmetric positive definite leading part of the equation.
Algorithm C0. 1. Find a global coarse grid solution $u_H\in S_0^H(\Omega )$:
The rest of the proof is similar to that of Theorem 4.2, and we leave the details to the interested readers.
■
4.2. New parallel algorithms based on local algorithms
The parallel algorithms we shall present here are naturally obtained from the local algorithms that we studied above. Given an initial coarse triangulation $T^H(\Omega )$, let us first divide $\Omega$ into a number of disjoint subdomains $D_1, \ldots , D_m$ (see Figure 5), then enlarge each $D_j$ to obtain $\Omega _j$ that align with $T^H(\Omega )$. The basic idea of our parallel algorithm is very simple: we just apply the local algorithms in parallel in all $\Omega _j$’s.
4.2.1. Basic parallel algorithms
Let us first discuss the parallel version of Algorithm A0. For each $j$, we use some adaptive process to obtain a shape-regular mesh $T_j(\Omega )$ and the corresponding finite element solution denoted by $u_j$. We note that each $T_j(\Omega )$ looks like the mesh depicted in Figure 1, and it has a substantially finer mesh inside $\Omega _j$. We note that all $T_j(\Omega )$ are different triangulations for $\Omega$ and they can be very arbitrary; but for simplicity of exposition, we assume each $T_j(\Omega )$ has the same size $h$ in $\Omega _j$ (more precisely, $T_j(\Omega _j)=T^h(\Omega _j)$) and has the size $H$ away from $\Omega _j$. Let $S_0^{h_j}(\Omega )\subset H_0^1(\Omega )$ be the corresponding finite element spaces.
We now discuss the parallel versions of Algorithms B0 and C0. Although there are many possibilities for the generalization, for clarity of exposition, it appears to be most convenient to discuss these using two globally defined grids: an initial coarse grid $T^H(\Omega )$ and a refined (from $T^H(\Omega )$) grid $T^h(\Omega )$ that satisfies $h\ll H$.
Algorithm B1. 1. Find a global coarse grid solution $u_H\in S_0^H(\Omega )$:
We note that the approximations $u^h$ obtained by Algorithms A1, B1 and C1 are piecewise defined and they are in general discontinuous. We also point out that $\|u_h-u^h\|_{0,\Omega }$ does not in general have higher order than $\Vvert u_h-u^h\Vvert _{1,\Omega }$. In this subsection, we shall propose some further modifications for these algorithms to achieve the following two goals:
(1)
smooth $u^h$ to obtain a global $H^1(\Omega )$ approximation;
(2)
improve the error $\|u_h-u^h\|_{0,\Omega }$.
The first goal will be achieved by using some more local fine grid problems; the second will be achieved by carrying out a further global coarse grid correction. We note that the second goal is realized after the first goal has been achieved.
We now proceed to present a modified algorithm that addresses both of the aforementioned two issues. To do this, we pick another sequence of subdomains $G_j\subset \subset D_j$ and $G_{m+1}=\Omega \setminus (\bigcup _{j=1}^m{\bar{G}_j})$ (see Figure 6).
Algorithm C2. 1. Find a global coarse grid solution $u_H\in S^{H}_0(\Omega )$:
3. Set $u^h=u_H+e_{h}^{j}$, in $G_j\ (j=1,\dots ,m)$ and $u^h$ on $\bar{G}_{m+1}$ is defined by $\displaystyle u^h\mid _{\partial G_j\cap \partial G_{m+1}}=u_H+e^{j}_h \ (j=1,\dots ,m)$ and satisfying
The local error estimates and local/parallel algorithms presented in previous sections make it possible to design many new local and parallel adaptive algorithms for finite element computations. In this section, we shall give some examples.
5.1. On the traditional adaptive process
Before we present our new local and parallel adaptive method, we would like to recall some traditional finite element adaptive process based on a posteriori error estimates. We would like to illustrate that our local a posteriori error estimates give a certain theoretical justification of some traditional finite element adaptive techniques.
The basic idea of an adaptive algorithm is to use a given computed finite element solution to detect the behavior of the exact solution so that the underlying finite element meshes can get properly refined or de-refined in certain regions of the domain according to the behavior of the solution. The behavior of the solution is detected by using certain a posteriori error estimates like the ones presented in §3.2. In the existing literature, these a posteriori error estimates are often given and analyzed in a global form. For example, in view of Equation 3.19, one can often use the following kind of global a posteriori error estimate:
The constant $C_0$ only depends on the shape of the grids, and it may be properly estimated, or, for simplicity, one may take $C_0=1$. In practice, we wish to find a mesh $T^h(\Omega )$ (with least possible number of nodes) such that the corresponding finite element approximation $u_h$ satisfies
If the above estimate is satisfied for the given mesh, then we have achieved our goal. Otherwise, we need to further refine the mesh locally. In order to use the global error estimate for local mesh refinement, one often uses the principle of equi-distribution for the error. Let $N_h$ be the total number of elements. For each element, we then check if the following is satisfied:
Here $\delta ^2/N_h$ is the averaged error on $\tau$ obtained by the aforementioned equi-distribution principle.
One natural question to ask is why a global a posteriori error estimate like Equation 5.1 can be used locally as in Equation 5.3. We would like to argue that our local estimates would give a theoretical justification of the aforementioned local application of a posteriori error estimates. One argument we can make is that the a posteriori error estimate itself is essentially local, according to Theorem 3.9. Another related argument can be based on the locality of a priori error estimates for $u-u_h$ according to Theorem 3.4.
5.2. A local adaptive process
The locality of both a priori and a posteriori error estimates can be used to devise some local/parallel adaptive algorithms. As an example, we shall now propose a local adaptive algorithm that can be used to obtain an approximate solution that has a desired accuracy in a given subdomain locally.
Let $D\subset \Omega$ be a given subdomain and $\delta$ a given tolerance. Suppose we want to obtain an approximate solution $u_h$ satisfying
According to the local error estimates in previous sections, this error tolerance can be achieved by only local mesh refinement around the domain $D$. Let $\Omega _0$ be a subdomain that is slightly larger than $D$. By Theorem 3.9, we have the estimate
$$\begin{equation} \|u-u_h\|_{1,D}\le C_0\|h\eta (u_h)\|_{0,\Omega _0} +\text{higher order global error term}. \cssId{LocalEh}{\tag{5.4}} \end{equation}$$
If our initial grid is reasonably fine, the higher order global error term is negligible. This means we may only need to refine the mesh in the subdomain $\Omega _0$ so that
This adaptive process corresponds precisely to the local algorithm described in §4.1.1, where a priori error estimates are discussed.
In view of the local algorithm described in §4.1.2, apparently, a corresponding local adaptive process can also be obtained. We omit the details here.
5.3. A parallel adaptive process
As before, a simultaneous application of a local algorithm on a number of subdomains naturally leads to a parallel algorithm. In this subsection, we shall give some details for a parallel adaptive algorithm corresponding to the local adaptive algorithm described in the previous subsection.
Our goal is to design a parallel procedure to find a finite element approximation $u_h$ (which may be piecewise defined) satisfying (see Equation 4.1)
Based on a reasonable good initial mesh, denoted by $T^{H}(\Omega )$, and its corresponding solution, denoted by $u_{H}$, an adaptive process is to make use of some a posteriori estimates based on information from $T^{H}(\Omega )$ and $u_0$ to adaptively come up with better and better meshes until a desired error tolerance is reached. Traditionally, after a stage of refinement/de-refinement, the a posteriori estimates are evaluated on the whole domain. Thanks to the local estimate Equation 3.21, we propose that a posteriori estimates can be evaluated concurrently on a number of proper subdomains and a parallel adaptive process can then be brought about.
As in §4.2, given an initial coarse triangulation $T^H(\Omega )$, we divide $\Omega$ into a number of disjoint subdomains $D_1, \ldots , D_m$, then enlarge each $D_j$ to obtain $\Omega _j$’s that align with $T^H(\Omega )$.
We aim to reach Equation 5.5 by refining the mesh $T^{H}(\Omega )$. Note that $\|u-u_{H}\|_{0,\Omega }$ is of higher order compared with $\|u-u_{H}\|_{1,\Omega }$; for convenience of exposition, we may assume that our initial mesh is fine enough so that $\|u-u_{H}\|_{0,\Omega }$ is much smaller than $\delta$. (This assumption is not crucial in practice, as $T^{H}(\Omega )$ can get updated by finer meshes in an adaptive process.) Thanks to the estimate Equation 3.21, we have (with $h=H$)
We use a standard principle of equi-distribution for the error control, in which we equalize the contribution from each subdomain. More precisely, the finite element approximation computed on the targeted mesh $T^h(\Omega )$ in terms of computational work satisfies
$$\begin{eqnarray} \|u-u_h\|^2_{1,D_j}\le \frac{\delta ^2}{m},\quad j=1,2,\cdots , m. \tag{5.6} \end{eqnarray}$$
Therefore we can just carry out the mesh refinement locally on $\Omega _j$ until the following estimates are satisfied:
We note that the refinement process on each $\Omega _j$ is independent of those from other subdomains. Associated with each $\Omega _j$, we get a locally refined mesh as in Figure 1 and then find the corresponding finite element solution, denoted by $u_j$, on such a mesh. After all these are done, we then take final solutions that are defined piecewise on each $D_j$ restricted from $\Omega _j$.
The above exposition contains the main ideas of a parallel adaptive process, but for its application, there are many practical issues that need to be addressed. We refer to Bank and Holst Reference 12 for a similar approach and implementation details.
6. Some numerical experiments
In this paper, we have presented many new estimates and new algorithms. It is perhaps a little too much of an undertaking to carry out and report numerical experiments for all these results in this single work. For illustration, we choose to report some simple numerical experiments only for Algorithms B1, C1 and C2.
We consider the simple unit square domain $\Omega =(0, 1)\times (0, 1)$ and a uniform triangulation $T^h(\Omega )=\{\tau \}$ (see Figure 7) and piecewise linear finite element spaces:
Now we apply Algorithm B1, Algorithm C1 and Algorithm C2 with coarse mesh size $H = \sqrt {h}$ to solve two partial differential equations of second order, respectively. For the exact solver of all the nonsymmetric and/or indefinite systems of coarse spaces, the Gaussian elimination is used. On the other hand, a standard V-cycle multigrid algorithm is used to solve all the SPD systems on fine spaces.
We first consider the following simple Poisson equation:
We shall apply Algorithm B1 and Algorithm C2 to solve this problem, using fine meshes of sizes $h=2^{-j}\ (j=6, 8, 10)$ and corresponding coarse meshes of size $H=\sqrt {h}$.
Let $u_h$ be the standard finite element solution, let $u^h$ be obtained by Algorithm B1, and let $\tilde{u}^h$ be obtained by Algorithm C2. Then, by Theorems 4.4 and 4.6, we obtain
$$\begin{equation} \Vvert u_h-u^h\Vvert _{1,\Omega }=\mathcal{O}(H^2)\approx c h,\quad \Vvert u_h-\tilde{u}^h\Vvert _{0,\Omega }=\mathcal{O}(H^3)\approx c h^{3/2}. \cssId{texmlid2}{\tag{6.2}} \end{equation}$$
The results shown in Table 1 support the above estimate. Actually, for Algorithm C2, this numerical example shows a better order of convergence than our theory predicted.
We next consider a simple example of convection-diffusion problems:
where $\mathbf{b}=(2x-e^y, 3y\cos (\pi x)),\nobreakspace f= 70\log ((x+1/10)(\sin (\pi y)+1))$.
As for the convection-diffusion equation above, we again apply Algorithm C1 and Algorithm C2 to solve this problem. The corresponding computational results are shown in Table 2, and again support our theory.
7. Some further remarks
In this last section, we shall make a few technical comments and a concluding remark.
7.1. On the dependence of subdomains
We would like to point out that most of the error estimates presented in this paper depend on the distance between the boundaries of the underlying subdomains (cf. Schatz and Wahlbin Reference 39Reference 40, and Wahlbin Reference 46Reference 47). To avoid notational complication, we chose not to explicitly spell out this kind of dependence.
7.2. Estimates in terms of different norms
Most of the local estimates in this paper can be generalized to other norms such as the $W^{1,\infty }$ and $L^\infty$ norms. As an illustration, let us discuss briefly the possible maximum norm estimates in two dimensions.
For any $z\in \Omega$, let $G_z$ be the Green function with respect to the singular point $z$:
In our local error estimates, the global errors are all measured in $L^2$ norms. As in the existing literature on local a priori error estimates (cf. Schatz and Wahlbin Reference 39Reference 40, and Wahlbin Reference 46Reference 47), it is possible to replace the global $L^2$ norm by some negative Sobolev norms. For example, the following estimate may be obtained:
For simplicity and generality, we did not get into details in this paper when the above improved estimates may be obtained. We will report this kind of results in our future work.
7.4. Conclusion
In this paper, we have used a simple second oder elliptic model problem and a class of finite element discretization methods to demonstrate how to use a coarse grid to capture the global component of the approximate solution and then parallelize the major computation in a much finer grid. We believe this is a general and powerful parallel-computing technique that can be used for a variety of partial differential equations with different types of discretization methods.
Acknowledgments
The authors wish to thank H. Kim and L. Zikatanov for their assistance with numerical experiments and helpful discussions. This paper was completed while the first author was visiting ETH Zurich in Switzerland during the summer of 1998, and the author wishes to thank Professors C. Schwab and R. Jeltsch for their hospitality.
Adams R.A.(1975): Sobolev Spaces, Academic Press, New York. MR 56:9247
Reference [2]
Ainsworth, M. and Oden, J.T.(1993): A unified approach to a posteriori error estimation using element residual methods, Numer. Math., 65, 23-50. MR 95a:65185
Reference [3]
Axelsson, O. and Layton, W.(1996): A two-level discretization of nonlinear boundary value problems, SIAM J. Numer. Anal., 33, 2359-2374. MR 98c:65181
Reference [4]
Babuska, I., Duran, R. and Rodriguez, R.(1992): Analysis of the efficiency of an posteriori error estimator for linear triangular finite elements, SIAM J. Numer. Anal., 29, 947-946. MR 93d:65096
Reference [5]
Babuska, I. and Rheinboldt, C.(1978): Error estimates for adaptive finite element computations, SIAM J. Numer. Anal., 15, 736-754. MR 58:3400
Reference [6]
Babuska, I., Zienkiewicz, O.C., Gago, J. and Oliveira, E.R. de A. (eds.)(1986): Accuracy Estimates and Adaptive Refinements in Finite Element Computations, Wiley, New York. MR 87j:65004
[7]
Babuska, I., Strouboulis, T. and Gangaraj, S.K.(1997): A posteriori estimation of the error in the recovered derivatives of the finite element solution, Comput. Methods Appl. Mech. Engrg., 150, 369-396. CMP 98:06
[8]
Babuska, I., Strouboulis, T., Gangaraj, S.K. and Upadhyay, C.S.(1997): Pollution error in the $h-$version of the finite element method and the local quality of the recovered derivatives, Comput. Methods Appl. Mech. Engrg., 140, 1-37. MR 97i:73092
[9]
Babuska, I., Strouboulis, T. and Upadhyay, C.S.(1994): A model study of the quality of a posteriori error estimators for linear elliptic problems. Error estimation in the interior of patchwise uniform grids of triangles, Comput. Methods Appl.Mech. Engrg., 114, 307-378. MR 95d:65093
Reference [10]
Bank, R.E.(1996): Hierarchical bases and the finite element method, Acta Numerica, 5, 1-43. CMP 98:14
[11]
Bank, R.E.(1998): A simple analysis of some a posteriori error estimates, Appl. Numer. Math., 26, 153-164. CMP 98:08
Reference [12]
Bank, R.E. and Holst, M.(1998): A new paradigm for parallel adaptive meshing algorithms (manuscript).
[13]
Bank, R.E. and Smith, R.K.(1993): A posteriori error estimates based on hierarchical bases, SIAM J. Numer. Anal., 30, 921-935. MR 95f:65212
Reference [14]
Bank, R.E. and Smith, R.K.(1997): Mesh smoothing using a posteriori error estimates, SIAM J. Numer. Anal., 34, 979-997. MR 98M:65162
Reference [15]
Bank, R.E. and Weiser, A.(1985): Some a posteriori error estimates for elliptic partial differential equations, Math. Comp., 44, 283-301. MR 86g:65207
Reference [16]
Bedivan, D.M.(1995): A two-grid method for solving elliptic problems with inhomogeneous boundary conditions, Comput. Math. Appl., 29, 59-66. MR 95k:65103
Reference [17]
Blum, H., Lin, Q. and Rannacher, R.(1986): Asymptotic error expansion and Richardson extrapolation for linear finite elements, Numer. Math., 49, 11-38. MR 87m:65172
[18]
Bornemann, F.A., Erdmann, B. and Kornhuber, R.(1996): A posteriori error estimates for elliptic problems in two and three space dimensions, SIAM J. Numer. Anal., 33, 1188-1204. MR 98a:65161
Reference [19]
Bramble, J.H.(1993): Multigrid Methods, Pitman Research Notes in Mathematics, 294, London Co-published in the USA with Wiley, New York. MR 95b:65002
Reference [20]
Bramble, J.H., Ewing, R.E., Parashkevov, R.R. and Pasciak, J.E.(1992): Domain decomposition methods for problems with partial refinement, SIAM J. Sci. Stat. Comp., 13, 397-410. MR 92i:65179
Reference [21]
Bramble, J.H., Ewing, R.E., Pasciak, J.E. and Schatz, A.H.(1988): A preconditioning technique for the efficient solution of problems with local grid refinement, Comp. Meth. Appl. Mech. Eng., 67, 149-159.
Reference [22]
Chan, T. and Mathew, T.(1994): Domain decomposition algorithms, Acta Numerica, 3, 61-143. MR 95f:65214
Reference [23]
Ciarlet, P.G. and Lions J.L.(1991): Handbook of Numerical Analysis, Vol.II, Finite Element Methods (Part I), North-Holland. MR 92f:65001
Reference [24]
Dawson, C.N. and Wheeler, M.F.(1994): Two-grid methods for mixed finite element approximations of nonlinear parabolic equations, Contemp. Math., 180, 191-203. MR 95j:65117
Reference [25]
Dawson, C.N., Wheeler, M.F. and Woodward, C.S.(1998): A two-grid finite difference scheme for nonlinear parabolic equations, SIAM J. Numer. Anal., 35, 435-452. MR 99b:65097
Reference [26]
Eriksson, K., Estep, D., Hansbo, P. and Johnson, C.(1995): Introduction to adaptive methods for differential equations, Acta Numerica, 105-158. MR 96k:65057
Reference [27]
Eriksson, K., Estep, D., Hansbo, P. and Johnson, C.(1996): Computational Differential Equations, Cambridge University Press. MR 97m:65006
Reference [28]
Eriksson, K. and Johnson, C.(1991): Adaptive finite element methods for parabolic problems I: a linear model problem, SIAM J. Numer. Anal., 28, 43-77. MR 91m:65274
[29]
Eriksson, K. and Johnson, C.(1995): Adaptive finite element methods for parabolic problems IV: Nonlinear problems, SIAM J. Numer. Anal., 32, 1729-1749. MR 96i:65081
Hackbusch, W.(1985): Multigrid Methods and Applications, Springer, New York. MR 87e:65082
Reference [32]
Johnson, C.(1990): Adaptive finite element methods for diffusion and convection problems, Comp. Methods Appl. Mech. Engrg., 82, 301-322. MR 91k:65134
Reference [33]
Layton, W. and Lenferink, W.(1995): Two-level Picard and modified Picard methods for the Navier-Stokes equations, Appl. Math. Comp., 69, 263-274. MR 95m:65191
Reference [34]
Marion, M. and Xu, J.(1995): Error estimates on a new nonlinear Galerkin method based on two-grid finite elements, SIAM J. Numer. Anal., 32, 1170-1184. MR 96f:65136
Reference [35]
Nitsche, J. and Schatz, A.H.(1974): Interior estimates for Ritz-Galerkin methods, Math. Comp., 28, 937-955. MR 51:9525
Reference [36]
Nochetto, R. H.(1995): Pointwise a posteriori error estimates for elliptic problems on highly graded meshes, Math. Comp., 64, 1-22. MR 95c:65172
Reference [37]
Rannacher, R. and Scott, R.(1982): Some optimal error estimates for piecewise linear finite element approximations, Math. Comp., 38, 437-445. MR 83e:65180
Reference [38]
Schatz, A.H.(1998): Pointwise error estimates and asymptotic error expansion inequalities for the finite element method on irregular grids: Part I. Global estimates, Math. Comp., 67, 877-899. MR 98j:65082
Reference [39]
Schatz, A.H. and Wahlbin, L.B.(1977): Interior maximum-norm estimates for finite element methods, Math. Comp., 31, 414-442. MR 55:4748
Reference [40]
Schatz, A.H. and Wahlbin, L.B.(1995): Interior maximum-norm estimates for finite element methods, Part II, Math. Comp., 64, 907-928. MR 95j:65143
[41]
Schatz, A.H. and Wang, J.(1996): Some new error estimates for Ritz-Galerkin methods with minimal regularity assumptions, Math. Comp., 65, 19-27. MR 96d:65190
Reference [42]
Utnes, T.(1997): Two-grid finite element formulations of the incompressible Navier-Stokes equations, Comm. Numer. Methods Engrg., 34, 675-684. MR 98d:76110
Reference [43]
Verfürth, R.(1994): A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations, Math. Comp., 62, 445-475. MR 94j:65136
Reference [44]
Verfürth, R.(1995): A posteriori error estimates for nonlinear problems. Finite element discretizations of parabolic equations, Bericht Nr. 180, Fakultät für Mathematik, Ruhr-Universität Bochum.
Reference [45]
Verfürth, R.(1996): A Review of A-Posteriori Error Estimation and Adaptive Mesh Refinement, Wiley-Teubner.
Reference [46]
Wahlbin, L.B.(1991): Local behavior in finite element methods, in Reference 23, pp. 355-522. MR 92f:65001
Reference [47]
Wahlbin, L.B.(1995): Superconvergence in Galerkin Finite Element Methods, Vol. 1605, Lecture Notes in Math., Springer. MR 98j:65083
Reference [48]
Xu, J.(1992): A new class of iterative methods for nonselfadjoint or indefinite problems, SIAM J. Numer. Anal., 29, 303-319. MR 92k:65063
Reference [49]
Xu, J.(1992): Iterative methods by space decomposition and subspace correction, SIAM Review, 34, 4, 581-613. MR 93k:65029
Reference [50]
Xu, J.(1994): A novel two-grid method for semilinear equations, SIAM J. Sci. Comput., 15, 231-237. MR 94m:65178
Reference [51]
Xu, J.(1996): Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal., 33, 1759-1777. MR 97i:65169
Reference [52]
Xu, J. and Zou, J.(1998): Some non-overlapping domain decomposition methods, SIAM Review 40, 4, 857-914.
Reference [53]
Yserentant, H.(1993): Old and new proofs for multigrid algorithms, Acta Numerica, 2, 285-326. MR 94i:65128
Reference [54]
Zhou, A., Liem, C.L., Shih, T.M. and Lü, T.(1998): Error analysis on bi-parameter finite elements, Comput. Methods Appl. Mech. Engrg., 158, 329-339. CMP 98:14
This work was partially supported by NSF DMS-9706949, NSF ACI-9800244 and NASA NAG2-1236 through Penn State and Center for Computational Mathematics and Applications, The Pennsylvania State University.
Journal Information
Mathematics of Computation, Volume 69, Issue 231, ISSN 1088-6842, published by the American Mathematical Society, Providence, Rhode Island.
Show rawAMSref\bib{1654026}{article}{
author={Xu, Jinchao},
author={Zhou, Aihui},
title={Local and parallel finite element algorithms based on two-grid discretizations},
journal={Math. Comp.},
volume={69},
number={231},
date={2000-07},
pages={881-909},
issn={0025-5718},
review={1654026},
doi={10.1090/S0025-5718-99-01149-7},
}
Settings
Change font size
Resize article panel
Enable equation enrichment
(Not available in this browser)
Note. To explore an equation, focus it (e.g., by clicking on it) and use the arrow keys to navigate its structure. Screenreader users should be advised that enabling speech synthesis will lead to duplicate aural rendering.