An unstructured mesh control volume method for two-dimensional space fractional diffusion equations with variable coefficients on convex domains

TL;DR

This paper introduces a novel unstructured mesh control volume method for two-dimensional space fractional diffusion equations on convex domains, offering improved efficiency and applicability over traditional finite element methods.

Contribution

The paper presents a new finite volume scheme for space fractional diffusion equations on arbitrary convex domains, including implementation details and an efficient solver for the resulting linear system.

Findings

The method reduces CPU time significantly compared to finite element methods.

The stiffness matrix is sparse and irregular, requiring specialized storage and solvers.

Numerical experiments confirm the method's accuracy, efficiency, and applicability to complex domains.

Abstract

In this paper, we propose a novel unstructured mesh control volume method to deal with the space fractional derivative on arbitrarily shaped convex domains, which to the best of our knowledge is a new contribution to the literature. Firstly, we present the finite volume scheme for the two-dimensional space fractional diffusion equation with variable coefficients and provide the full implementation details for the case where the background interpolation mesh is based on triangular elements. Secondly, we explore the property of the stiffness matrix generated by the integral of space fractional derivative. We find that the stiffness matrix is sparse and not regular. Therefore, we choose a suitable sparse storage format for the stiffness matrix and develop a fast iterative method to solve the linear system, which is more efficient than using the Gaussian elimination method. Finally, we present several examples to verify our method, in which we make a comparison of our method with the finite element method for solving a Riesz space fractional diffusion equation on a circular domain. The numerical results demonstrate that our method can reduce CPU time significantly while retaining the same accuracy and approximation property as the finite element method. The numerical results also illustrate that our method is effective and reliable and can be applied to problems on arbitrarily shaped convex domains.

Tables

Table 1. Table 1: The size and density of matrix 𝐌 𝐌 \mathbf{M} for different h ℎ h on a square domain [ 0 , 1 ] × [ 0 , 1 ] 0 1 0 1 [0,1]\times[0,1]

$h$	Size	Density
5.2693E-01	4$\times$4	100%
3.1123E-01	15$\times$15	86.667%
1.6759E-01	64$\times$64	57.715%
8.6682E-02	258$\times$258	34.002%
4.3719E-02	1115$\times$1115	17.705%
2.3063E-02	5255$\times$5255	8.517%

Table 2. Table 2: The L 2 subscript 𝐿 2 L_{2} error, L ∞ subscript 𝐿 L_{\infty} error, convergence order and CPU time of h ℎ h with τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} for the linear coefficient case at t = 1 𝑡 1 t=1

	$h$	$L_2$ error	Order	$L_{\infty }$ error	Order	Time
	3.1123E-01	3.5684E-04	–	1.4774E-03	–	4.90s
$\alpha =0.3$	1.6759E-01	1.0880E-04	1.92	4.3735E-04	1.97	19.50s
$\beta =0.5$	8.6682E-02	2.2391E-05	2.40	1.3895E-04	1.74	2.30min
	4.3719E-02	6.9379E-06	1.71	3.7632E-05	1.91	28.42min
	3.1123E-01	3.7935E-04	–	1.4827E-03	–	4.91s
$\alpha =0.4$	1.6759E-01	1.2435E-04	1.80	4.2971E-04	2.00	19.98s
$\beta =0.8$	8.6682E-02	2.5152E-05	2.42	1.3725E-04	1.73	2.36min
	4.3719E-02	7.2675E-06	1.81	3.5722E-05	1.97	28.56min
	3.1123E-01	3.9259E-04	–	1.3844E-03	–	4.91s
$\alpha =0.7$	1.6759E-01	1.4100E-04	1.65	4.1957E-04	1.93	19.87s
$\beta =0.9$	8.6682E-02	2.8670E-05	2.42	1.4117E-04	1.65	2.37min
	4.3719E-02	7.5385E-06	1.95	3.3666E-05	2.09	28.47min

Table 3. Table 3: The L 2 subscript 𝐿 2 L_{2} error, L ∞ subscript 𝐿 L_{\infty} error, convergence order and CPU time of h ℎ h with τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} for the quadratic coefficient case at t = 1 𝑡 1 t=1

	$h$	$L_2$ error	Order	$L_{\infty }$ error	Order	Time
	3.1123E-01	3.1608E-04	–	1.3430E-03	–	4.97s
$\alpha =0.3$	1.6759E-01	1.0064E-04	1.85	4.0906E-04	1.92	20.48s
$\beta =0.5$	8.6682E-02	2.0661E-05	2.40	1.3852E-04	1.64	2.45min
	4.3719E-02	6.2709E-06	1.74	3.7584E-05	1.91	28.69min
	3.1123E-01	3.6299E-04	–	1.4108E-03	–	4.88s
$\alpha =0.4$	1.6759E-01	1.2145E-04	1.77	4.1614E-04	1.97	20.51s
$\beta =0.8$	8.6682E-02	2.4646E-05	2.42	1.3823E-04	1.67	2.46min
	4.3719E-02	6.7517E-06	1.89	3.3858E-05	2.06	28.78min
	3.1123E-01	3.8524E-04	–	1.3424E-03	–	4.97s
$\alpha =0.7$	1.6759E-01	1.3952E-04	1.64	4.0669E-04	1.93	20.56s
$\beta =0.9$	8.6682E-02	2.8522E-05	2.41	1.4126E-04	1.60	2.44min
	4.3719E-02	7.1520E-06	2.02	3.1880E-05	2.17	28.68min

Table 4. Table 4: The L 2 subscript 𝐿 2 L_{2} error, L ∞ subscript 𝐿 L_{\infty} error, convergence order and CPU time of h ℎ h with τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} for the exponential coefficient case at t = 1 𝑡 1 t=1

	$h$	$L_2$ error	Order	$L_{\infty }$ error	Order	Time
	3.1123E-01	5.1809E-04	–	1.9033E-03	–	4.97s
$\alpha =0.3$	1.6759E-01	1.6296E-04	1.87	5.3973E-04	2.04	20.62s
$\beta =0.5$	8.6682E-02	3.8817E-05	2.18	1.6032E-04	1.84	2.45min
	4.3719E-02	1.1574E-05	1.77	4.8226E-05	1.76	28.46min
	3.1123E-01	4.5022E-04	–	1.6750E-03	–	4.93s
$\alpha =0.4$	1.6759E-01	1.4896E-04	1.79	1.0117E-04	2.01	20.52s
$\beta =0.8$	8.6682E-02	3.4126E-05	2.24	4.8309E-04	1.84	2.45min
	4.3719E-02	1.1238E-05	1.62	4.3016E-05	1.76	28.66min
	3.1123E-01	4.2412E-04	–	1.4994E-03	–	4.93s
$\alpha =0.7$	1.6759E-01	1.5286E-04	1.65	4.6520E-04	1.89	20.50s
$\beta =0.9$	8.6682E-02	3.3401E-05	2.31	1.4533E-04	1.76	2.45min
	4.3719E-02	1.0565E-05	1.68	4.0322E-05	1.87	28.56min

Table 5. Table 5: Comparison of the consumed CPU time of Gaussian elimination versus Bi-CGSTAB

$N_e$	$h$	Gauss elimination	Bi-CGSTAB
44	3.1123E-01	4.90s	4.90s
158	1.6759E-01	22.57s	19.50s
578	8.6682E-02	5.39min	2.30min
2356	4.3719E-02	5.48h	28.42min

Table 6. Table 6: The comparison of the density of stiffness matrix generated by FEM and CVM for different h ℎ h

$N_e$	$h$	Size	FEM	CVM
174	2.8917E-01	$74\times 74$	65.413 %	55.332%
570	1.6444E-01	$260\times 260$	41.814 %	33.521%
2310	8.6550E-02	$1104\times 1104$	22.233 %	17.469%
8744	4.5873E-02	$4271\times 4271$	11.712 %	9.107%

Table 7. Table 7: The L 2 subscript 𝐿 2 L_{2} error, L ∞ subscript 𝐿 L_{\infty} error and convergence order of h ℎ h for FEM with τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} at t = 1 𝑡 1 t=1

FEM	$h$	$L_2$ error	Order	$L_{\infty }$ error	Order
	2.8917E-01	6.7022E-03	–	5.8841E-03	–
$\alpha =0.80$	1.6444E-01	2.0787E-03	2.07	2.8557E-03	1.28
$\beta =0.80$	8.6550E-02	5.2077E-04	2.16	8.1791E-04	1.95
	4.5873E-02	1.3554E-04	2.12	2.3520E-04	1.96
	2.8917E-01	6.9018E-03	–	5.5925E-03	–
$\alpha =0.70$	1.6444E-01	2.1713E-03	2.05	2.7718E-03	1.24
$\beta =0.90$	8.6550E-02	5.4452E-04	2.16	7.9048E-04	1.95
	4.5873E-02	1.4147E-04	2.12	2.2242E-04	2.00

Table 8. Table 8: The L 2 subscript 𝐿 2 L_{2} error, L ∞ subscript 𝐿 L_{\infty} error and convergence order of h ℎ h for CVM with τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} at t = 1 𝑡 1 t=1

CVM	$h$	$L_2$ error	Order	$L_{\infty }$ error	Order
	2.8917E-01	1.4782E-02	–	2.1786E-02	–
$\alpha =0.80$	1.6444E-01	4.5014E-03	2.11	7.5230E-03	1.88
$\beta =0.80$	8.6550E-02	1.2275E-03	2.02	1.8279E-03	2.20
	4.5873E-02	3.4069E-04	2.02	5.4557E-04	1.90
	2.8917E-01	1.4950E-02	–	2.1864E-02	–
$\alpha =0.70$	1.6444E-01	4.5530E-03	2.11	7.6462E-03	1.86
$\beta =0.90$	8.6550E-02	1.2566E-03	2.01	1.8659E-03	2.20
	4.5873E-02	3.4898E-04	2.02	5.4606E-04	1.94

Table 9. Table 9: The comparison of running time between FEM and CVM for different h ℎ h with α = β = 0.80 𝛼 𝛽 0.80 \alpha=\beta=0.80 , τ = 10 − 3 𝜏 superscript 10 3 \tau=10^{-3} at t = 1 𝑡 1 t=1

$N_e$	$h$	FEM	CVM
174	2.8917E-01	3.49 min	35.01 s
570	1.6444E-01	12.90 min	2.63min
2310	8.6550E-02	1.38 h	28.41min
8744	4.5873E-02	17.89h	6.59h

Equations

$$\frac{\partial u(x,y,t)}{\partial t}=$$

$$+$$

$$u(x,y,0)=\varphi (x,y),(x,y)\in \overline{\Omega },$$

$$u(x,y,t)=0,(x,y,t)\in \partial \Omega \times [0,T],$$

$$\frac{{\partial }^{\alpha }u(x,y,t)}{\partial x^{\alpha }}$$

$$\frac{{\partial }^{\alpha }u(x,y,t)}{\partial (-x{)}^{\alpha }}$$

$$\frac{{\partial }^{\beta }u(x,y,t)}{\partial y^{\beta }}$$

$$\frac{{\partial }^{\beta }u(x,y,t)}{\partial (-y{)}^{\beta }}$$

$$\frac{{\partial }^{\alpha }u(x,y,t)}{\partial x^{\alpha }}$$

$$\frac{{\partial }^{\alpha }u(x,y,t)}{\partial (-x{)}^{\alpha }}$$

$$\frac{{\partial }^{\beta }u(x,y,t)}{\partial y^{\beta }}$$

$$\frac{{\partial }^{\beta }u(x,y,t)}{\partial (-y{)}^{\beta }}$$

$$K_1(x,y,t)$$

$$K_3(x,y,t)$$

$$\frac{\partial u(x,y,t)}{\partial t}=K_x\frac{{\partial }^{1+\alpha }u(x,y,t)}{\partial |x{|}^{1+\alpha }}+K_y\frac{{\partial }^{1+\beta }u(x,y,t)}{\partial |y{|}^{1+\beta }}+f(x,y,t),$$

$$\frac{{\partial }^{1+\alpha }u(x,y,t)}{\partial |x{|}^{1+\alpha }}$$

$$\frac{{\partial }^{1+\beta }u(x,y,t)}{\partial |y{|}^{1+\beta }}$$

$$A(u,v)$$

$$\displaystyle+\frac{{K_{y}}}{2\cos\frac{\pi(1+\beta)}{2}}\Big{\{}\Big{(}{{}_{c(x)}{D}_{y}^{\frac{(1+\beta)}{2}}u},{{}_{y}{D}_{d(x)}^{\frac{(1+\beta)}{2}}v}\Big{)}+\Big{(}{{}_{y}{D}_{d(x)}^{\frac{(1+\beta)}{2}}u},{{}_{c(x)}{D}_{y}^{\frac{(1+\beta)}{2}}v}\Big{)}\Big{\}}.$$

$$\displaystyle\frac{K_{x}}{2\cos\frac{\pi(1+\alpha)}{2}}\oint_{\Gamma_{i}}\bigg{[}\frac{\partial^{\alpha}u(x,y,t)}{\partial x^{{\alpha}}}-\frac{\partial^{\alpha}u(x,y,t)}{\partial(-x)^{{\alpha}}}\bigg{]}{d}{y}$$

$$-$$

$$M_h=\{P_i:P_i\text{is a vertex of the element}K\in {\mathcal{T}}_h\text{and}{P}_i\in \Omega \},$$

$${\int }_{V_i}\frac{\partial u(x,y,t)}{\partial t}d{V}_i=$$

$$+$$

$$\displaystyle\Delta{V_{i}}\frac{\partial u(x,y,t)}{\partial t}\bigg{|}_{(x_{i},y_{i})}=$$

$$-$$

$$+$$

$$\Delta V_i=\sum _{j=1}^{m_i}\Delta V_{ij},$$

$$\displaystyle\oint_{\Gamma_{i}}\bigg{[}K_{1}(x,y,t)\frac{\partial^{\alpha}u(x,y,t)}{\partial x^{{\alpha}}}-K_{2}(x,y,t)\frac{\partial^{\alpha}u(x,y,t)}{\partial(-x)^{{\alpha}}}\bigg{]}d{y}$$

$$=$$

$$\displaystyle\oint_{\Gamma_{i}}\bigg{[}K_{3}(x,y,t)\frac{\partial^{\beta}u(x,y,t)}{\partial y^{{\beta}}}-K_{4}(x,y,t)\frac{\partial^{\beta}u(x,y,t)}{\partial(-y)^{{\beta}}}\bigg{]}d{x}$$

$$=$$

$$\displaystyle\Delta{V_{i}}\frac{\partial u(x,y,t)}{\partial t}\bigg{|}_{(x_{i},y_{i})}$$

$$=$$

$$-$$

$$+$$

$$\frac{\partial u(x,y,t_n)}{\partial t}=\frac{u(x,y,t_n)-u(x,y,t_{n-1})}{\tau }+O(\tau ).$$

$$u^p(x,y)=\sum _{j=1}^3{u}_j{\phi }_j(x,y)+O(h^2),$$

$$\displaystyle\varphi_{j}(x,y)\Big{|}_{(x,y)\in e_{p}}=\frac{1}{2\Delta_{e_{p}}}(a_{j}\,x+b_{j}~{}y+c_{j}),\quad\varphi_{j}(x,y)\Big{|}_{(x,y)\notin e_{p}}=0,$$

Full text

††journal: arXiv.org

1 Introduction

In the past two decades, fractional differential equations have been applied in many fields of science [1, 2, 3, 4, 5, 6, 7], in which space fractional diffusion equations are used to model the anomalous transport of solute in groundwater hydrology [8, 9]. For space fractional diffusion equations with constant coefficients, analytical solutions can be obtained by utilising the Fourier transform methods. However, many practical problems involve variable coefficients [10, 11], in which the diffusion velocity can vary over the solution domain. The work involving space fractional diffusion equations with variable coefficients is numerous. Meerschaert et al. [8, 12] considered the finite difference method for the one-dimensional one-sided and two-sided space fractional diffusion equations with variable coefficients, respectively. Zhang et al. [13] explored the homogeneous space-fractional advection-dispersion equation with space-dependent coefficients. Ding et al. [14] presented the weighted finite difference methods for a class of space fractional partial differential equations with variable coefficients. Moroney and Yang [15, 16] proposed some fast preconditioners for the numerical solution of a class of two-sided nonlinear space-fractional diffusion equations with variable coefficients. Chen and Deng [17] discussed the alternating direction implicit method to solve a two-dimensional two-sided space fractional convection-diffusion equation on a finite domain. Wang and Zhang [18] developed a high-accuracy preserving spectral Galerkin method for the Dirichlet boundary-value problem of a one-sided variable-coefficient conservative fractional diffusion equation. Feng et al. [19] proposed the finite volume method for a two-sided space-fractional diffusion equation with variable coefficients. Chen et al. [20] considered an inverse problem for identifying the fractional derivative indices in a two-dimensional space-fractional nonlocal model with variable diffusivity coefficients. Jia and Wang [21] presented a fast finite volume method for conservative space-fractional diffusion equations with variable coefficients. In [22], Feng et al. presented a new second order finite difference scheme for a two-sided space-fractional diffusion equation with variable coefficients.

In fact, many mathematical models and problems from science and engineering must be computed on irregular domains and therefore seeking effective numerical methods to solve these problems on such domains is important. Although existing numerical methods for fractional diffusion equations are numerous [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34], most of them are limited to regular domains and uniform meshes. Research involving unstructured meshes and irregular domains is sparse. Yang et al. [35] proposed the finite volume scheme for a two-dimensional space-fractional reaction-diffusion equation based on the fractional Laplacian operator $-(-\nabla^{2})^{\frac{\alpha}{2}}$, which was computed using unstructured triangular meshes on a unit disk. Burrage et al. [36] developed some techniques for solving fractional-in-space reaction diffusion equations using the finite element method on both structured and unstructured grids. Qiu et al. [37] developed the nodal discontinuous Galerkin method for fractional diffusion equations on a two-dimensional domain with triangular meshes. Liu et al. [38] presented the semi-alternating direction method for a two-dimensional fractional FitzHugh-Nagumo monodomain model on an approximate irregular domain. Qin et al. [39] also used the implicit alternating direction method to solve a two-dimensional fractional Bloch-Torrey equation using an approximate irregular domain. Karaa et al. [40] proposed a finite volume element method implemented on an unstructured mesh for approximating the anomalous subdiffusion equations with a temporal fractional derivative. Yang et al. [41] established the unstructured mesh finite element method for the nonlinear Riesz space fractional diffusion equations on irregular convex domains. Fan et al. [42] extended the unstructured mesh finite element method developed by Yang et al. [41] to the time-space fractional wave equation. Feng et al. [43] investigated the unstructured mesh finite element method for a two-dimensional time-space Riesz fractional diffusion equation on irregular arbitrarily shaped convex domains and a multiply-connected domain. Le et al. [44] studied the finite element approximation for a time-fractional diffusion problem on a domain with a re-entrant corner. To the best of our knowledge, the control volume finite element method (see Carr et al. [45] for an illustration of the method applied to wood drying) has not been generalised to allow the solution of space fractional diffusion equations with variable coefficients.

In this paper, we will consider the unstructured mesh control volume method for the following two-dimensional space fractional diffusion equation with variable coefficients (2D SFDE-VC) [20] on an arbitrarily shaped convex domain:

[TABLE]

subject to the initial condition

[TABLE]

and boundary conditions

[TABLE]

where $0<\alpha,~{}\beta<1$, $K_{i}(x,y,t)\geq 0$, $i=1,2,3,4$, $f(x,y,t)$ and $\phi(x,y)$ are assumed to be two known smooth functions. When the solution domain is rectangular $\Omega=(a,b)\times(c,d)$, we define the Riemman-Liouville fractional derivative as [46]:

[TABLE]

When the boundary of the solution domain is nonconstant or curved, for example a convex domain shown in Figure 1 with left boundary $a(y)$, right boundary $b(y)$, lower boundary $c(x)$ and upper boundary $d(x)$, we define the Riemman-Liouville fractional derivative as [43]:

[TABLE]

Remark 1.1**.**

When $K_{i}(x,y,t)~{}i=1,2,3,4$ take the special form

[TABLE]

equation (1) can be written as the following Riesz space fractional diffusion equation [38, 41]

[TABLE]

where

[TABLE]

One important application of equation (4) is in the study of cardiac arrhythmias. In two dimensions, the fractional FitzHugh-Nagumo monodomain model can be rewritten as a two-dimensional Riesz space fractional reaction-diffusion model, which can be used to describe the propagation of the electrical potential in heterogeneous cardiac tissue [38, 47]. This electrophysiological model of the heart can describe how electrical currents flow through the heart controlling its contraction and can be used to ascertain the effects of certain drugs designed to treat heart problems.

The major contribution of this paper is as follows.

• $\bullet$

Different from [35] and [40], we consider the control volume method for the two-dimensional space fractional diffusion equation with variable coefficients, in which the space fractional operator is either the Riemman-Liouville fractional derivative or Riesz space fractional derivative. To the best of our knowledge, this is a new contribution to the literature.

• $\bullet$

We propose a novel technique utilizing the control volume method implemented with an unstructured triangular mesh to deal with the space fractional derivative on an irregular convex domain, which we believe provides a very flexible solution strategy because our considered solution domain can be arbitrarily convex. Compared to the finite difference method in [38, 39], our method requires fewer grid nodes to generate the meshes in the solution domain partition.

• $\bullet$

For the methods considered in this paper, we construct the control volumes using triangular meshes and transform the problem (1) from the solution domain to a single control volume. Then we integrate problem (1) over an arbitrary control volume and change the control volume integral to a line integral over the control volume faces, which is approximated by the midpoint approximation. Moreover, we utilise the linear basis function to approximate the fractional derivatives at the midpoints of the control volume faces, in which some numerical techniques are used to handle the non-locality of the fractional derivative of the basis function.

• $\bullet$

We explore the property of the stiffness matrix generated by the integral of space fractional derivative. We find that the stiffness matrix is sparse and not regular. Specially, the more small the maximum edge of the triangulation is, the more sparse of the stiffness matrix becomes. Therefore, we choose a suitable sparse storage format for the stiffness matrix and utilise the bi-conjugate gradient stabilized method (Bi-CGSTAB) iterative method to solve the linear system, which is more efficient than using the Gaussian elimination method.

• $\bullet$

We present several examples to verify our method, in which we make a comparison of our method with the finite element method proposed in [41] for solving the Riesz space fractional diffusion equation (4) on a circular domain. In [41], the authors develop an algorithm to form the stiffness matrix and derive the bilinear operator as

[TABLE]

The bilinear form involves 8 fractional derivative terms and the approximation of two-fold multiple integrals, which are approximated by Gauss quadrature. While for the control volume method, we use the following form to generate the stiffness matrix form,

[TABLE]

in which we only need to calculate 4 fractional derivative terms and the approximation of line integrals. The numerical results demonstrate that our method can reduce CPU time significantly while retaining the same accuracy and approximation property as the finite element method. The numerical results also illustrate that our method is effective and reliable and can be applied to problems on arbitrarily convex domains.

The outline of this paper is as follows. In Section 2, the unstructured mesh control volume method for problem (1) is proposed and the full implementation details are provided. Then the property of the stiffness matrix is explored and a fast iterative solver is developed for the linear system. In Section 3, several numerical examples are presented to verify the effectiveness of the method and comparisons are made with existing methods to highlight its computational performance. Finally, some conclusions of the work are drawn.

2 Control volume finite element method

In this section, we will generalise the control volume method to solve equation (1), placing particular emphasis on the way the Riemman-Liouville fractional derivatives are discretised in space. Firstly, we divide the solution domain $\Omega$ into a number of regular triangular regions. Let $\mathcal{T}_{h}$ denote this triangulation and $h$ be the maximum diameter of the triangular elements. Then we introduce the control volumes, which are constructed as follows. Let $M_{h}$ be a set of vertice,

[TABLE]

and $M_{h}^{0}$ be the set of interior nodes in $\mathcal{T}_{h}$. We denote $P_{0}$ as the interior node of the triangulation $\mathcal{T}_{h}$ and $P_{i}~{}(i=1,2,\cdots,m)$ as its adjacent nodes (see Figure 2 with $m=6$). Let $S_{i}~{}(i=1,2,\cdots,m)$ be the midpoints of the line segments $\overline{P_{0}P_{i}}$ and $Q_{i}~{}(i=1,2,\cdots,m)$ the barycenters of the triangle $\Delta P_{0}P_{i}P_{i+1}$ with $P_{m+1}=P_{1}$. The control volume $K_{P_{0}}^{*}$ is constructed by joining successively $S_{1},~{}Q_{1},~{}\cdots,~{}S_{m},~{}Q_{m},~{}S_{1}$ (see Figure 2). We call the line segments $\overline{S_{i}Q_{i}}$ and $\overline{Q_{i}S_{i+1}}$ ($i=1,2,\cdots,m$ and $S_{m+1}=S_{1})$ control volume faces. Consequently, each of the triangular elements is divided into three sub-domains by these control surfaces. These quadrilateral shapes are called sub-control volumes and are illustrated in Figure 2 (for example, the quadrilateral $S_{1}Q_{1}S_{2}P_{0}$). Thus, a control volume consists of the sum of all neighbouring sub-control volumes that surround the given node $P_{0}$. The control volume is polygonal in shape and can be assembled in a straightforward and efficient manner at the element level. The flow across each control surface must be determined by an integral. Therefore, the finite volume method discretization process is initiated by utilising the integrated form of equation (1).

Integrating (1) over an arbitrary control volume $V_{i}$ ($i=1,2,\cdots,N_{p}$), yields

[TABLE]

Utilising a lumped mass approach for the time derivative and source term and applying Green's theorem to the other two integrals terms, gives

[TABLE]

where $\Gamma_{i}$ is the boundary of control volume $V_{i}$. We assume the finite volume integration is an anticlockwise traversal and the outward unit normal surface vector to the control surface with $\Delta x=x_{b}-x_{a}$ and $\Delta y=y_{b}-y_{a}$. Denote $\Delta V_{i}$ and $\Delta V_{ij}$ the area of the control volume and the sub-control volume surrounding the point $(x_{i},y_{i})$, then we have

[TABLE]

where $m_{i}$ is the total number of sub-control volumes that make up the control volume associated with the node $i$. The integral term on the right-hand side of equation (1) is a line integral, which can be approximated by the midpoint approximation for each control surface. Hence, the first integral term in equation (6) can be rewritten as

[TABLE]

where $(x_{r},y_{r})$ is the mid-point of the control face (CF). Similarly, for the second integral term in equation (6), we have

[TABLE]

Substituting equations (7) and (8) into (6), we obtain

[TABLE]

To discretise the time derivative in equation (9) at $t=t_{n}$, we use the backward Euler difference scheme

[TABLE]

In the following, we discuss the spatial discretisation of $u(x,y,t_{n})$. We consider the computation process for piecewise linear polynomials on the triangular element $e_{p}$, $p=1,2,...,N_{e}$, where $N_{e}$ is the total number of triangles. Then, within element $e_{p}$, the field function $u^{p}(x,y)$ can be written as

[TABLE]

where the triangle vertices are numbered in a counter-clockwise order as $1,~{}2,~{}3$ and the basis function $\varphi_{j}(x,y)$ is defined as

[TABLE]

where $\Delta_{e_{p}}$ is the area of triangle element $p$. It is well-known that

[TABLE]

where $\delta$ is the Kronecker function. With these local field functions and basis functions, we can obtain a global approximation of $u(x,y)$ for the whole triangulation:

[TABLE]

where $l_{k}(x,y)$ is the new basis function whose support domain is $\Omega_{e_{k}}$ (see Figure 3 the green polygonal domain) and $N_{p}$ is the total number of vertices on the convex domain $\Omega$.

Now, we denote $u_{h}(x,y,t_{n})$ as the approximation solution of $u(x,y,t_{n})$ and write $u_{h}(x,y,t_{n})$ in the form

[TABLE]

where $u_{k}^{n}$ are the coefficients that are to be solved for. Substituting equations (10) and (11) into equation (9), we discretise equation (9) at $t=t_{n}$ as follows:

[TABLE]

Using the fact that

[TABLE]

we obtain

[TABLE]

Equation (13) can be written in the following matrix form

[TABLE]

where $\mathbf{A}=$diag$~{}[\Delta V_{1},\Delta V_{2},\ldots,\Delta V_{N_{p}}]$, $\mathbf{U}^{n}=[u_{1}^{n},u_{2}^{n},\ldots,u_{N_{p}}^{n}]^{T}$, $\mathbf{F}^{n}=[f(x_{1},y_{1},t_{n}),f(x_{2},$ $y_{2},t_{n}),\ldots,f(x_{N_{p}},y_{N_{p}},t_{n})]^{T}$. Rearranging we obtain

[TABLE]

To form matrix $\mathbf{M}$, we need to calculate the fractional derivative of the basis function $l_{k}(x,y)$. In the following, we focus on the calculation of $\frac{\partial^{\alpha}l_{k}(x,y)}{\partial x^{{\alpha}}}$, $\frac{\partial^{\alpha}l_{k}(x,y)}{\partial(-x)^{{\alpha}}}$, $\frac{\partial^{\beta}l_{k}(x,y)}{\partial y^{{\beta}}}$ and $\frac{\partial^{\beta}l_{k}(x,y)}{\partial(-y)^{{\beta}}}$ at $(x_{r},y_{r})$. To evaluate \frac{\partial^{\alpha}l_{k}(x,y)}{\partial x^{{\alpha}}}\big{|}_{(x_{r},y_{r})} and \frac{\partial^{\alpha}l_{k}(x,y)}{\partial(-x)^{{\alpha}}}\big{|}_{(x_{r},y_{r})}, suppose that line $y=y_{r}$ intersects $n_{q}$ points with the support domain $\Omega_{e_{k}}$ of $l_{k}(x,y)$ (see Figure 3 with $n_{q}=5$).

Then we have

[TABLE]

Using the important observation that

[TABLE]

where $\varphi_{kp}(x,y)$ is the basis function of node $k$ on the triangular element $e_{p}$, we obtain

[TABLE]

As $l_{k}(x,{y_{r}})$ is a linear function on each sub integral interval, equation (16) can be evaluated using integration by parts over each sub integral interval. For the right fractional derivative of $l_{k}(x,{y_{r}})$ at $(x_{r},y_{r})$, we obtain

[TABLE]

Now we consider the case that point $(x_{r},y_{r})$ is in the support domain $\Omega_{e_{k}}$ of $l_{k}(x,y)$. Suppose that line $y=y_{r}$ intersects $n_{q}$ points with the support domain $\Omega_{e_{k}}$ (see Figure 4 with $n_{q}=4$). In this case, we have

[TABLE]

Then

[TABLE]

and

[TABLE]

If line $y=y_{r}$ intersects zero points with the support domain $\Omega_{e_{k}}$, then we have

[TABLE]

The calculation of $\frac{\partial^{\beta}l_{k}(x,y)}{\partial y^{{\beta}}}$ and $\frac{\partial^{\beta}l_{k}(x,y)}{\partial(-y)^{{\beta}}}$ at $(x_{r},y_{r})$ can be derived in a similar manner for the $y$ direction. Finally, we summarise the whole computation process in the following algorithm (see Algorithm 1).

Remark 2.1**.**

When the boundary of the solution domain is nonconstant or curved, all of the above calculation is applicable as well.

Here, we discuss the structure of matrix $\mathbf{M}$. Firstly, the matrix $\mathbf{M}$ generated by scheme (13) is sparse and not regular. Then we explore the sparsity of matrix $\mathbf{M}$ for different $h$. Table 1 shows the size and density (nonzero entries percentage) of matrix $\mathbf{M}$ for different $h$ where we can observe that as $h$ decreases the density of matrix $\mathbf{M}$ reduces significantly. We can infer that when $h$ is small enough, matrix $\mathbf{M}$ is extremely sparse and this facilitates the use of a sparse matrix storage format to reduce the memory usage of our computational method. Furthermore, we employ an efficient sparse iterative solver Bi-CGSTAB [48] to solve the linear system (15) (see Algorithm 2), which is more efficient than using Gaussian elimination method. The CPU time comparison of the two methods is studied numerically in Example 3.1.

3 Discussion of Numerical Results

In this section, we provide some numerical examples to verify the effectiveness of our method presented in Section 2. We adopt linear polynomials on triangles and define $h$ as the maximum length of the triangle edges. $N_{e}$ is taken as the number of triangles in $\mathcal{T}_{h}$. Here, the numerical computations were carried out using MATLAB R2014b on a Dell desktop with configuration: Intel(R) Core(TM) i7-4790, 3.60 GHz and 16.0 GB RAM. We use the following formula to calculate the convergence order:

[TABLE]

Example 3.1**.**

Firstly, we consider the following 2D SFDE-VC on a rectangular domain

[TABLE]

subject to

[TABLE]

where $\Omega=(0,1)\times(0,1)$, $T=1$,

[TABLE]

This is a two-dimensional anomalous diffusion model, which can describe anomalous transport in heterogeneous porous media and can be used to explain the region-scale anomalous dispersion with heavy tails [20].

The exact solution of this problem is given by $u(x,y,t)=(t^{2}+1)x^{2}(1-x)^{2}y^{2}(1-y)^{2}$. Here, we consider three different coefficient cases [22]: linear coefficients $K_{1}(x,y,t)=2-x$, $K_{2}(x,y,t)=2+x$, $K_{3}(x,y,t)=2-y$, $K_{4}(x,y,t)=2+y$, quadratic coefficients $K_{1}(x,y,t)=2-x^{2}$, $K_{2}(x,y,t)=2+x^{2}$, $K_{3}(x,y,t)=2-y^{2}$, $K_{4}(x,y,t)=2+y^{2}$ and exponential coefficients $K_{1}(x,y,t)=3-e^{x}$, $K_{2}(x,y,t)=3+e^{x}$, $K_{3}(x,y,t)=3-e^{y}$, $K_{4}(x,y,t)=3+e^{y}$. The numerical results are given in Tables 2 to 4. Table 2 illustrates the $L_{2}$ error, $L_{\infty}$ error and corresponding convergence order of $h$ for the linear coefficient case for different $\alpha$, $\beta$ with $\tau=10^{-3}$ at $t=1$. Tables 3 and 4 show the $L_{2}$ error, $L_{\infty}$ error and corresponding convergence order of $h$ for the quadratic coefficient case and exponential coefficient case, respectively. From these tables we can see that the convergence order of both the $L_{2}$ error and $L_{\infty}$ error is $2-\max\{\alpha,\beta\}$ order [19] and the numerical results are in excellent agreement with the exact solution, which demonstrates the effectiveness of the numerical method. We can also observe that with $h$ deceasing, the CPU time grows considerably, which we believe is mainly due to the non-locality of the fractional derivative of the basis function and the computational cost to generate the matrix $\mathbf{M}$. In addition, we give a comparison between the Bi-CGSTAB and Gaussian elimination. In the Bi-CGSTAB solver, we set $10^{-10}$ as the stopping criterion and the maximum iteration number is $10^{2}$. Table 5 displays the consumed CPU time of these two algorithms at $t=1$ with $\tau=10^{-3}$, $\alpha=0.3$, $\beta=0.5$, $K_{1}(x,y,t)=2-x$, $K_{2}(x,y,t)=2+x$, $K_{3}(x,y,t)=2-y$, $K_{4}(x,y,t)=2+y$ for different $h$. Compared to the Gaussian elimination, Bi-CGSTAB has significantly reduced 90% of the computational time for $h=4.3719E-02$. Another advantage of Bi-CGSTAB to be mentioned is that the average iteration number does not appear to increase significantly as $h$ decreases. Here, the average iteration number is approximately $10$ regardless of the model dimensions. We conclude that the Bi-CGSTAB solver is more efficient than Gaussian elimination for solving this problem.

Example 3.2**.**

Next, we consider the following two-dimensional Riesz space fractional diffusion equation on a circular domain, which can be used to describe the propagation of the electrical potential in heterogeneous cardiac tissue [38, 41, 47].

[TABLE]

where $\Omega=\{(x,y)|{x^{2}}+{y^{2}}<1\}$, $K_{x}=1$, $K_{y}=1$, $T=1$,

[TABLE]

The exact solution is given by $u(x,y,t)=e^{-t}({x^{2}}+{y^{2}}-1)^{2}$. Figure 5 shows the circular domain partitioned by unstructured triangular meshes and control volumes for different $h$. In [41], Yang et al. applied the Galerkin finite element method for solving the two-dimensional Riesz space fractional diffusion equation with a nonlinear source term on convex domains. They developed an algorithm to form the stiffness matrix on triangular meshes, which can deal with space fractional derivatives on any convex domain. Here, we will make a comparison between our method (CVM) and Yang's method (FEM) for solving the two-dimensional Riesz space fractional diffusion equation (21) on a circular domain using the same triangular meshes. Firstly, we present a comparison of the density of the two stiffness matrices generated by FEM and CVM for different $h$ in Table 6. We can see that with $h$ decreasing the density of the two stiffness matrices reduces significantly. Compared to the stiffness matrix generated by FEM, the stiffness matrix generated by CVM is slightly more sparse. Next, we present a comparison of the error and convergence. Table 7 displays the $L_{2}$ error, $L_{\infty}$ error and corresponding convergence order of $h$ for different $\alpha$, $\beta$ with $\tau=10^{-3}$ at $t=1$ by applying FEM. Table 8 highlights the error and convergence order by using FVM. We can see that the accuracy of our method is similar to FEM, both of which are second order. Then, we present a comparison of CPU time for the two methods in Table 9 both using the Bi-CGSTAB solver. We choose $\alpha=\beta=0.8$ and $\tau=10^{-3}$ at $t=1$ to observe the running time for different $h$. We observe that compared to the running time of FEM, CVM can reduce the running time significantly, which illustrates that CVM is more effective for solving the two-dimensional Riesz space fractional diffusion equation on convex domains. This is mainly due to the bilinear form in [41] that involves 8 fractional derivative terms and the approximation of two-fold multiple integrals, which are approximated by Gauss quadrature, while for CVM we only need to calculate 4 fractional derivative terms and the approximation of line integrals. We can see that the numerical solution is in excellent agreement with the exact solution, which demonstrates the effectiveness of our numerical method again.

4 Conclusions

In this paper, we considered the unstructured mesh control volume method for the two-dimensional space fractional diffusion equation with variable coefficients on convex domains. We partitioned the irregular convex domain using triangular meshes. Then we constructed the control volumes and solved the space fractional diffusion equation by utilising the finite volume method. Finally, numerical examples on irregular convex domains were studied, which verified the effectiveness and reliability of the method. We concluded that the numerical method can be extended to other arbitrarily shaped convex domains. Furthermore, according to the property of the stiffness matrix generated by the finite volume method, we chose a suitable sparse matrix format for the stiffness and utilised the Bi-CGSTAB iterative method to solve the linear system, which is more efficient than using Gauss elimination method. In addition, we made a comparison of our method with the finite element method proposed in [41], which demonstrated that our method can reduce CPU time significantly while retaining the same accuracy and approximation property as the finite element method. In future work, we shall investigate the unstructured mesh control volume method applied to other fractional problems on irregular convex domains, such as the two-dimensional multi-term time-space fractional diffusion equation with variable coefficients or three-dimensional space fractional diffusion equations with variable coefficients.