Uniform Convergent Tailored Finite Point Method for Advection–Diffusion Equation with Discontinuous, Anisotropic and Vanishing Diffusivity

Abstract

We propose two tailored finite point methods for the advection–diffusion equation with anisotropic tensor diffusivity. The diffusion coefficient can be very small in one direction in some part of the domain and be discontinuous at the interfaces. When flows advect from the vanishing-diffusivity region towards the non-vanishing diffusivity region, standard numerical schemes tend to cause spurious oscillations or negative values. Our proposed schemes have uniform convergence in the vanishing diffusivity limit, even when the solution exhibits interface and boundary layers. When the diffusivity is along the coordinates, the positivity and maximum principle can be proved. We use the value as well as their derivatives at the grid points to construct the scheme for nonaligned case, which makes it can achieve good accuracy and convergence for the derivatives as well, even when there exhibit boundary or interface layers. Numerical experiments are presented to show the performance of the proposed scheme.

Appendix

Here, we show that the limiting scheme for the degenerate equation is consist and explain heuristically the monotonicity. We only consider the case when $\beta _1=\beta _2=0$ and $\theta \ne 0$. In the limit that $\epsilon _{i,j}\rightarrow 0$ in (25),

$$\begin{aligned} r_1\rightarrow \gamma _1=-\alpha \cos ^2\theta ,\quad r_2\rightarrow \gamma _2=-\alpha \sin ^2\theta ,\quad r_3\rightarrow \gamma _3=2\alpha \sin \theta \cos \theta . \end{aligned}$$

Let the coefficient matrix $D_k$ in (34) be $(D_k^1,~D_k^2),~k=1,2,3,4$, where $D_k^i$ ($i=1,2$) are $3\times 1$ columns. We need to show that (34) becomes a consistent and monotone discretization for

$$\begin{aligned} \displaystyle \gamma _1\frac{\partial ^2 v}{\partial x^2}+\gamma _2\frac{\partial ^2 v}{\partial y^2}+\gamma _3\frac{\partial ^2 v}{\partial x\partial y}+v=0. \end{aligned}$$

(49)

Due to the form of $\gamma _1$, $\gamma _2$, $\gamma _3$, (49) can be rewritten as

$$\begin{aligned} -\alpha \frac{\partial ^2v}{\partial \xi ^2}+v=0,\qquad \text{ with } \xi =x\cos \theta +y\sin \theta . \end{aligned}$$

(50)

If we are able to show that when $\epsilon _{i,j}\rightarrow 0$, the limit of (34) is a discretization to (50), the proof completes.

For the degenerate case we consider here, $\beta _{1,i,j}=\beta _{2,i,j}=\epsilon _{i,j}=0$, so that $(\lambda _{i,j}, \mu _{i,j})$ in (31) become $(\pm \sqrt{\frac{1}{\alpha \cos ^2\theta }},0)$, $(0,\pm \sqrt{\frac{1}{\alpha \sin ^2\theta }})$. Let

$$\begin{aligned} \lambda =\sqrt{\frac{1}{\alpha \cos ^2\theta }}, \qquad \mu =\sqrt{\frac{1}{\alpha \sin ^2\theta }}, \end{aligned}$$

(51)

the linear space $H_{i,j}$ becomes

$$\begin{aligned}&\{e^{\lambda x},e^{-\lambda x},e^{\mu y},e^{-\mu y},(x\sin \theta -y\cos \theta )e^{\lambda x},(-x\sin \theta +y\cos \theta )e^{-\lambda x},\\&(x\sin \theta -y\cos \theta )e^{\mu y}, (-x\sin \theta +y\cos \theta )e^{-\mu y}\} \end{aligned}$$

When $(\lambda _{k,i,j},\mu _{k,i,j})=(\pm \lambda ,0)$, $(\lambda _{k,i,j},\mu _{k,i,j})=(0,\pm \mu )$, (35) are respectively

$$\begin{aligned}&\begin{aligned} \left( \begin{array}{l}1\\ \pm \lambda \\ 0\end{array}\right)&= D_1\left( \begin{array}{l}1\\ 0\end{array}\right) e^{\pm \lambda h_x} +D_2\left( \begin{array}{l}1\\ 0\end{array}\right) e^{\mp \lambda h_x} +D_3\left( \begin{array}{l}1\\ \pm \lambda \end{array}\right) +D_4\left( \begin{array}{l}1\\ \pm \lambda \end{array}\right) ,\\&=D_1^1e^{\pm \lambda hx}+D_2^1e^{\mp \lambda hx}+D_3^1+D_4^1\pm \lambda (D_3^2+D_4^2). \end{aligned} \end{aligned}$$

(52)

$$\begin{aligned}&\begin{aligned} \left( \begin{array}{l}1\\ 0\\ \pm \mu \end{array}\right)&= D_1\left( \begin{array}{l}1\\ \pm \mu \end{array}\right) +D_2\left( \begin{array}{l}1\\ \pm \mu \end{array}\right) +D_3\left( \begin{array}{l}1\\ 0\end{array}\right) e^{\pm \mu h_y} +D_4\left( \begin{array}{l}1\\ 0\end{array}\right) e^{\mp \mu h_y},\\&=D_1^1+D_2^1\pm \mu (D_1^2+D_2^2)+D_3^1e^{\pm \mu h_y}+D_4^1e^{\mp \mu h_y}. \end{aligned} \end{aligned}$$

(53)

and (36) are respectively (after dividing some constant)

$$\begin{aligned}&\begin{aligned} \left( \begin{array}{l}0\\ \sin \theta \\ -\cos \theta \end{array}\right)&= D_1\left( \begin{array}{l}h_x\sin \theta \\ -\cos \theta \end{array}\right) e^{\pm \lambda h_x} +D_2\left( \begin{array}{l} -h_x\sin \theta \\ -\cos \theta \end{array}\right) e^{\mp \lambda h_x}\\&\quad +D_3\left( \begin{array}{l} h_y\cos \theta \\ \sin \theta \mp h_y\lambda \cos \theta \end{array}\right) +D_4\left( \begin{array}{l}-h_y\cos \theta \\ \sin \theta \pm h_y\lambda \cos \theta \end{array}\right) .\\&=h_x\sin \theta (D_1^1e^{\pm \lambda h_x}-D_2^1e^{\mp \lambda h_x})-\cos \theta (D_1^2e^{\pm \lambda h_x}+D_2^2e^{\mp \lambda h_x})\\&\quad -h_y\cos \theta (D_3^1-D_4^1)+\sin \theta (D_3^2+D_4^2)\mp h_y\lambda \cos \theta (D_3^2-D_4^2) \end{aligned} \end{aligned}$$

(54)

$$\begin{aligned}&\begin{aligned} \left( \begin{array}{l}0\\ \sin \theta \\ -\cos \theta \end{array}\right)&= D_1\left( \begin{array}{l}h_x\sin \theta \\ -\cos \theta \pm h_x\mu \sin \theta \end{array}\right) +D_2\left( \begin{array}{l}-h_x\sin \theta \\ -\cos \theta \mp h_x\mu \sin \theta \end{array}\right) \\&\quad +D_3\left( \begin{array}{l} -h_y\cos \theta \\ \sin \theta \end{array}\right) e^{\pm \mu h_y} +D_4\left( \begin{array}{l}h_y\cos \theta \\ \sin \theta \end{array}\right) e^{\mp \mu h_y}.\\&=h_x\sin \theta (D_1^1-D_2^1)-\cos \theta (D_1^2+D_2^2)\pm h_x\mu \sin \theta (D_1^2-D_2^2)\\&\quad -h_y\cos \theta (D_3^1e^{\pm \mu h_y}-D_4^1e^{\mp \mu h_y})+\sin \theta (D_3^2e^{\pm \mu h_y}+D_4^2e^{\mp \mu h_y}) \end{aligned} \end{aligned}$$

(55)

Let

$$\begin{aligned} D_k^i=\big (D_{k,1}^i, D_{k,2}^i,D_{k,3}^i\big )^\intercal , \qquad k=1,2,3,4,\quad i=1,2. \end{aligned}$$

(56)

The subtraction of the two equations in (52), the subtraction of the two equations in (53) and the summation of the two equations in (54), the summation of the two equations in (54) yield a linear system with scalar coefficnets for

$$\begin{aligned} D_1^1-D_2^1,\quad D_3^1-D_4^1,\quad D_1^2+D_2^2,\quad D_3^2+D_4^2. \end{aligned}$$

Since the right hand side of the subsystem for $D_{1,1}^1-D_{2,1}^1$, $D_{3,1}^1-D_{4,1}^1$, $D_{1,1}^2+D_{2,1}^2$ and $D_{3,1}^3+D_{4,1}^2$ is zero, we get

$$\begin{aligned} D_{1,1}^1=D_{2,1}^1=g_1,\quad D_{3,1}^1=D_{4,1}^1=g_2,\quad D_{1,1}^2=-D_{2,1}^2=g_3,\quad D_{3,1}^3=-D_{4,1}^2=g_4. \end{aligned}$$

The summation of the two equations in (52), the summation of the two equations in (53) and the subtraction of the two equations in (54), the subtraction of the two equations in (54) yield a linear system with scalar coefficients for

$$\begin{aligned} D_1^1+D_2^1,\quad D_3^1+D_4^1,\quad D_1^2-D_2^2,\quad D_3^2-D_4^2, \end{aligned}$$

which provides the following system for $g_1$, $g_2$, $g_3$, $g_4$

$$\begin{aligned} \left\{ \begin{array}{ll} \displaystyle g_1\cosh (\lambda h_x)+g_2=\frac{1}{2}, \\ \displaystyle g_1+g_2\cosh (\mu h_y)=\frac{1}{2},\\ \displaystyle g_1h_x\sin \theta \sinh (\lambda h_x)+g_3\cos \theta \sinh (\lambda h_x)+g_4\lambda h_y\cos \theta =0,\\ \displaystyle g_3\mu h_x\sin \theta +g_2h_y\cos \theta \sinh (\mu h_y)+g_4\sin \theta \sinh (\mu h_y)=0.\end{array} \right. \end{aligned}$$

(57)

The first line of the scheme (34) is

$$\begin{aligned} \displaystyle u_{i,j}= & {} g_1 u_{i+1,j}+g_3\partial _y u_{i+1,j}+g_1 u_{i-1,j}-g_3\partial _y u_{i-1,j}+g_2 u_{i,j+1}+g_4\partial _x u_{i,j+1}\nonumber \\&+g_2 u_{i,j-1}-g_4\partial _x u_{i,j-1}. \end{aligned}$$

(58)

The system (57) can be solved exactly such that

$$\begin{aligned} \left\{ \begin{array}{ll} \displaystyle g_1=\frac{1}{2}\frac{\cosh (\mu h_y)-1}{\cosh (\lambda h_x)\cos (\mu h_y)-1}, \\ \displaystyle g_2=\frac{1}{2}\frac{\cosh (\lambda h_x)-1}{\cosh (\lambda h_x)\cos (\mu h_y)-1},\\ \displaystyle g_3=\frac{g_1h_x\sin ^2\theta \sinh (\lambda h_x)\sinh (\mu h_y)-g_2\lambda h_y^2\cos ^2\theta \sinh (\mu h_y)}{\cos \theta \sin \theta \sinh (\lambda h_x)\sinh (\mu h_y)-\lambda \mu \sin \theta \cos \theta h_x h_y},\\ \displaystyle g_4=\frac{-g_2h_y\cos ^2\theta \sinh (\lambda h_x)\sinh (\mu h_y)+g_1\mu h_x^2\sin ^2\theta \sinh (\lambda h_x)}{\cos \theta \sin \theta \sinh (\lambda h_x)\sinh (\mu h_y)-\lambda \mu \sin \theta \cos \theta h_x h_y}. \end{array} \right. \end{aligned}$$

Then by Taylor expansion, when $\theta$ is away from $n\pi /2$ ($n=0,1,2,3$) and the mesh is fine enough,

$$\begin{aligned} \begin{aligned} \displaystyle g_1&\approx 1/2\frac{ (\mu h_y)^2}{(\lambda h_x)^2+(\mu h_y)^2}, \quad \displaystyle g_2\approx 1/2\frac{ (\lambda h_x)^2}{(\lambda h_x)^2+(\mu h_y)^2},\\ \displaystyle \frac{g_3}{\frac{\sin \theta }{\cos \theta }h_x}&\approx \frac{ (\lambda h_x)^2}{(\lambda h_x)^2+(\mu h_y)^2}g_1, \quad \displaystyle \frac{g_4}{\frac{\cos \theta }{\sin \theta }h_y}\approx \frac{ (\mu h_y)^2}{(\lambda h_x)^2+(\mu h_y)^2}g_2. \end{aligned} \end{aligned}$$

(59)

• The second order convergence. Let

  $$\begin{aligned} g=\frac{\cosh (\lambda h_x)\cos (\mu h_y)-1}{(\cosh (\lambda h_x)-1)(\cosh (\mu h_y)-1)} \end{aligned}$$

  Multiplying both sides of (58) by g and using (51), (59), we can obtain

  $$\begin{aligned} \begin{aligned}&\displaystyle -\alpha \cos ^2\theta \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2}-\alpha \sin ^2\theta \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_y^2}\\&\displaystyle -2\alpha \sin ^3\theta \cos \theta \frac{\partial _y u_{i+1,j}-\partial _y u_{i-1,j}}{2h_x}-2\alpha \cos ^3\theta \sin \theta \frac{\partial _x u_{i,j+1}-\partial _x u_{i,j-1}}{2h_y}\\&+u_{i,j}=O(h_x^2+h_y^2). \end{aligned} \end{aligned}$$

  (60)

  Therefore the limiting scheme is a second order discretization for (50).

• Heuristic argument about the monotonicity. When $\theta$ is away from $n\pi /2$ with $n=0,1,2,3$, i.e. $\sin \theta ,\cos \theta \ne 0$, let

  $$\begin{aligned} \begin{aligned} \digamma _1&= u_{i+1,j}+\frac{\sin \theta }{\cos \theta }h_x\partial _y u_{i+1,j}-2u_{i,j}+u_{i-1,j}-\frac{\sin \theta }{\cos \theta }h_x\partial _y u_{i-1,j},\\ \digamma _2&= u_{i,j+1}+\frac{\cos \theta }{\sin \theta }h_y\partial _x u_{i,j+1}-2u_{i,j}+u_{i,j-1}-\frac{\cos \theta }{\sin \theta }h_y\partial _x u_{i,j-1},\\ \digamma _3&=\left( g_1-\frac{g_3}{\frac{\sin \theta }{\cos \theta }h_x}\right) \left( u_{i+1,j}+u_{i-1,j}\right) +\left( g_2-\frac{g_4}{\frac{\cos \theta }{\sin \theta }h_y}\right) (u_{i,j+1}+u_{i,j-1})-u_{i,j}. \end{aligned} \end{aligned}$$

  (61)

  We can reformulate (58) into

  $$\begin{aligned} \displaystyle \frac{g_3}{\frac{\sin \theta }{\cos \theta }h_x}\digamma _1+\frac{g_4}{\frac{\cos \theta }{\sin \theta }h_y} \digamma _2+\digamma _3+\left( 2\frac{g_3}{\frac{\sin \theta }{\cos \theta }h_x}+2\frac{g_4}{\frac{\cos \theta }{\sin \theta }h_y}\right) u_{i,j}=0. \end{aligned}$$

  (62)

  According to the definition of monotone schemes for degenerate elliptic equation in [16], in order to have a consist and convergent scheme for the degenerate elliptic equation, we have to show that the coefficients in front of all neighboring grids are positive while the coefficient in front of $u_{i,j}$ are negative. For simplicity we only consider the case when $\theta \in (0,\pi /2)$. First of all, $g_1,g_2,g_3,g_4>0$ in this case.

  $$\begin{aligned} \begin{aligned} u_{i+1,j}+\frac{\sin \theta }{\cos \theta }h_x\partial _y u_{i+1,j},\qquad&u_{i-1,j}-\frac{\sin \theta }{\cos \theta }h_x\partial _y u_{i-1,j},\\ u_{i,j+1}+\frac{\cos \theta }{\sin \theta }h_y\partial _x u_{i,j+1},\qquad&u_{i,j-1}-\frac{\cos \theta }{\sin \theta }h_y\partial _x u_{i,j-1} \end{aligned} \end{aligned}$$

  can be respectively considered as approximations to the function values at $Z_{i,j}^1$, $Z_{i,j}^2$, $Z_{i,j}^3$, $Z_{i,j}^4$ in Fig. 12. Besides, from (59) and $g_1,g_2>0$, the coefficients in front of $u_{i\pm 1,j}$, $u_{i,j\pm 1}$ in $\digamma _3$ are positive. Finally, the coefficient in front of $u_{i,j}$ is $-1$ which is negative, therefore the limiting scheme is monotone, which gives the convergence of the scheme for the degenerate equation according to [16]. When $\theta$ is in other quadrant, the discussion is the same.

  The approximations in (59) can only be valid the mesh have to be fine enough, especially when $\theta$ is really close to $n\pi /2$. One may argue that when $\theta$ is close to $n\pi /2$ with $n=0,1,2,3$, the locations of $Z_{i,j}^1$, $Z_{i,j}^2$, $Z_{i,j}^3$, $Z_{i,j}^4$ can be far away from the cell. This is not a proof but a heuristic argument. However, the numerical results seem quite promising.

Fig. 12

The limiting scheme