Non-unified Compact Residual-Distribution Methods for Scalar Advection–Diffusion Problems

Abstract

This paper solves the advection–diffusion equation by treating both advection and diffusion residuals in a separate (non-unified) manner. An alternative residual distribution (RD) method combined with the Galerkin method is proposed to solve the advection–diffusion problem. This Flux-Difference RD method maintains a compact-stencil and the whole process of solving advection–diffusion does not require additional equations to be solved. A general mathematical analysis reveals that the new RD method is linearity preserving on arbitrary grids for the steady-state advection–diffusion equation. The numerical results show that the flux difference RD method preserves second-order accuracy on various unstructured grids including highly randomized anisotropic grids on both the linear and nonlinear scalar advection–diffusion cases.

Appendix by Rémi Abgrall

Institut für Mathematik Universität Zürich, Zurich, Switzerland. E-mail: remi.abgrall@math.uzh.ch

The problem is defined in ${\varOmega }\subset \mathbb {R}^2$. The case $\mathbb {R}^3$ could be done the same by replacing the 1/2 factor to 1/3. The boundary conditions could also be dealt similarly.

The residual is

$$\begin{aligned} \phi _i^T=\frac{1}{2} (\mathbf {f}-\mathbf {f}^\star _T)\cdot \mathbf {n}_i+\nu \int _T\mathbf {\nabla }\varphi _i\cdot \mathbf {\nabla } u \;d\mathbf {x}+\phi _i^{\text {art},T} \end{aligned}$$

(43)

where $\varphi _i$ is the basis function, $\mathbf {n}_i=\int _T\nabla \varphi _i \;d \mathbf {x}$. The as in [12,4], consider a test function v and its interpolant $v^h=\sum \limits _iv_i\varphi _i$. The scheme is

$$\begin{aligned} \sum \limits _{T, i\in T}\phi _i^T= -\sum \limits _{T, i\in T} \mathbf {f}^\star _T\cdot \mathbf {n}_i+\nu \int _{{\varOmega }}\nabla \varphi _i\cdot \nabla u d\mathbf {x}+\phi ^{\text {art}}_i \end{aligned}$$

(44)

where

$$\begin{aligned} \phi ^{\text {art}}_i=\sum \limits _{T, i\in T}\phi _i^{\text {art},T}. \end{aligned}$$

Then we multiply by $v_i$, sum over all degrees of freedom and get:

$$\begin{aligned} \begin{aligned} 0&=\sum \limits _{i}v_i \left( -\sum \limits _{T, i\in T} \mathbf {f}^\star _T\cdot \mathbf {n}_i+\nu \int _{{\varOmega }}\nabla \varphi _i\cdot \nabla u d\mathbf {x}+\phi ^{\text {art}}_i\right) \\&=-\sum _{T}\left( \sum _{j\in T} v_j\mathbf {n}_j\right) \mathbf {f}^\star _T +\nu \int _{\varOmega }\nabla v^h\cdot \nabla u d\mathbf {x}+\sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T} \end{aligned} \end{aligned}$$

(45)

Similarly, if $u^h$ is the piecewise linear interpolant of the exact solution, we can define the truncation error as,

$$\begin{aligned} \text {TE}=-\sum _{T}\left( \sum _{j\in T} v_j\mathbf {n}_j\right) \mathbf {f}^\star _T(u^h) +\nu \int _{\varOmega }\nabla v^h\cdot \nabla u^h d\mathbf {x}+\sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}(u^h). \end{aligned}$$

Now we also know that the exact solution $u^{ex}$ satisfies

$$\begin{aligned} -\sum _T \int _T\nabla v^h f( u^{ex}) d\mathbf {x}+\nu \int _T\nabla v^h\cdot \nabla u^{ex}\; d\mathbf {x}=0. \end{aligned}$$

So taking the difference, the truncation error is,

$$\begin{aligned} \begin{aligned} \text {TE}&=-\sum _{T}\left( \sum _{j\in T} v_j\mathbf {n}_j\cdot \mathbf {f}^\star _T\left( u^{ex,h}\right) \right) +\nu \int _{\varOmega }\nabla v^h\cdot \nabla u^{ex,h} d\mathbf {x}+\sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}( u^{ex,h})\\&=-\sum _{T}\left( \sum _{j\in T} \int _T \nabla v^h\cdot \mathbf {f}^\star _T\left( u^{ex,h}\right) \; d\mathbf {x}\right) +\nu \int _{\varOmega }\nabla v^h\cdot \nabla u^{ex,h} d\mathbf {x}\\&\quad +\sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}\left( u^{ex,h}\right) \\&=-\sum _{T}\left( \sum _{j\in T} \int _T \nabla v^h\cdot \left( \mathbf {f}^\star _T\left( u^{ex,h}\right) -\mathbf {f}( u^{ex})\right) \; d\mathbf {x}\right) {+}\nu \int _{\varOmega }\nabla v^h\cdot \nabla \left( u^{ex,h}- u^{ex}\right) \; d\mathbf {x}\\&\quad +\sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}\left( u^{ex,h}\right) . \end{aligned} \end{aligned}$$

(46)

Then, clearly if v is sufficiently regular, $\mathbf {f}^\star _T( u^{ex,h})-\mathbf {f}( u^{ex})=O(h^2)$ (at least for the choices of the paper), and $\nabla \big ( u^{ex,h}- u^{ex}\big )=O(h)$. This indicates that

$$\begin{aligned} -\sum _{T}\left( \sum _{j\in T} \int _T \nabla v^h\cdot \left( \mathbf {f}^\star _T\left( u^{ex,h}\right) -\mathbf {f}\left( u^{ex}\right) \right) \right. \; d\mathbf {x}= & {} O(h^2), \\ \nu \int _{\varOmega }\nabla v^h\cdot \nabla \left( u^{ex,h}- u^{ex}\right) \; d\mathbf {x}= & {} O(h) \end{aligned}$$

The second inequality shows it is O(h) but using the Aubin–Nitsche approach [4], we have that

$$\begin{aligned} \sup \limits _{v\in L^2} \frac{ \int _{\varOmega }\nabla v^h\cdot \nabla \left( u^{ex,h}- u^{ex}\right) \; d\mathbf {x}}{||v||_{L^2}} =O(h^2), \end{aligned}$$

and using Poincaré inequality (possible since v has a compact support), we get that $||v||_{L^2}\le C ||\nabla v||_2$ for some $C>0$ that only depends on ${\varOmega }$. Collecting the two together, we see that the viscous term

$$\begin{aligned} \nu \int _{\varOmega }\nabla v^h\cdot \nabla \left( u^{ex,h}- u^{ex}\right) \; d\mathbf {x}\end{aligned}$$

behaves like $O(h^2)$ in reality.

So in generality the scheme has a truncation error $O(h^2)$ provided that

$$\begin{aligned} \sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}( u^{ex,h})=O(h^2) \end{aligned}$$

This is also true under the assumptions of the paper since on each element,

$$\begin{aligned} \sum _{i\in T} \phi _i^{\text {art},T}\left( u^{ex,h}\right) =0, \end{aligned}$$

and then,

$$\begin{aligned} \sum _{T}\sum _{j\in T} v_j \phi _i^{\text {art},T}\left( u^{ex,h}\right) =\sum _T\sum _{j\in T} (v_i-v_T)v_j \phi _i^{\text {art},T}\left( u^{ex,h}\right) , \end{aligned}$$

where $v_T$ is the value at one arbitrarily chosen degree of freedom. In this paper, $\phi _i^{\text {art},T}( u^{ex,h})$ is a sum of differences of $u^{ex}$ multiplied by coefficient that are $O(h^2)$ (before the normalisation, since the normalisation removes an h and the tilded coefficients are O(h), so in summary $O(h^2)$ is preserved. Hence the Flux-Difference approach is LP in full generality and taking into account the diffusion. This however is a necessary but not a sufficient condition to preserve the order of accuracy on advection–diffusion problems.