Simple Positivity Preserving Nonlinear Finite Volume Scheme for Subdiffusion Equations on General Non-conforming Distorted Meshes

We propose a positivity preserving ﬁnite volume scheme on non-conforming quadrilateral distorted meshes with hanging nodes for subdiffusion equations, where the differential equations have a sum of time-fractional derivatives of different orders, and the typi-cal solutions of the problem have a weak singularity at the initial time t = 0 for given smooth data. In this paper, a positivity-preserving nonlinear method with centered unknowns is obtained by the two-point ﬂux technique, where a new method to handling vertex-unknown including hanging nodes is the highlight of our paper. For each time derivative, we apply the L1 scheme on a temporal graded mesh. Especially, the existence of a solution is strictly proved for the nonlinear system by applying the Brouwer’s ﬁxed point theorem. Numerical results show that the proposed positivity-preserving method is effective for strongly anisotropic and heterogeneous full tensor subdiffusion coefﬁcient problems.

such a property can lead to nonphysical oscillations or nonphysical numerical solutions, especially on distorted meshes. For example, the temperatures in Kelvin, temperature in thermal conduction problem and other physical variables governed by some subdiffusion equations are generally positive quantities. For Lagrangian computation of hydrodynamic problems and reservoir simulations, meshes are usually large deformation and the subdiffusion coefficients are heterogeneous and anisotropic, which may cause the discrete schemes more readily to produce non-physical negative solutions and oscillations. Therefore, effective numerical methods with positivity-preserving property is regarded as an indispensable requirement in constructing discrete schemes in order to avoid these problem. In the present paper we consider a more complex model, where several time fractional derivatives present in the differential equation.
We shall consider the multi-term subdiffusion problem where (a) l k is the given positive constants, u = u(x,t) is the solution of the multi-term subdiffusion equation. Without loss of generality we can consider l k = 1. (c) κ(x,t) is 2 × 2 diffusion tensor (possibly anisotropic and discontinuous on some interfaces) such that it is piecewise C 1 on Ω and is uniformly bounded above and below on Ω . That is to say, there exist two positive constants λ min > 0 and λ max > 0 such that the set of eigenvalues of κ(x,t) is included in [λ min , λ max ] with λ min > 0.

(d)
The source term f , the initial data v, and the boundary data ψ are piecewise smooth on Ω × (0, T ], respectively. It is well known that large relative displacement tends to occur near multi-material interface or sliding interface for high temperature and high pressure Lagrangian radiation hydrodynamic problems. Because the computational grid moves with fluid flow, which causes general non-conforming meshes (see Figure 1) to occur naturally, especially result in the appearance of hanging nodes. Therefore, efficient numerical methods with positivitypreserving property on general non-conforming meshes are necessary for solving these problems. Positivity is an important feature of the sub-diffusion equation. Recently, Luchko [1][2][3] established a series of interesting works, proved that the solution u(t) of problem (1.2) satisfied the property of positivity-preserving. Hence, it is natural to expect that numerical solutions could maintain the discrete analogy of this property of exact solution. So far, few numerical methods have considered the positivity preservation for fractional differential equation, especially on various distorted meshes. Ye et al. [20] used a fractional predictorcorrector method combining the L1 and L2 discrete schemes for the 1D multi-term timespace Riesz-Caputo fractional differential equations and proved a maximum principle on a uniform grid. Brunner et al. [21] provided a finite difference maximum principle preserving scheme for the fractional diffusion equation by introducing an equivalent definition of the Caputo fractional derivative, but there is no numerical implementation. Jin et al. [10] studied three positivity preservation numerical methods with some special meshes and geometric restriction for the problem (1.2)-(1.4), including the standard Galerkin (SG) scheme, the lumped mass (LM) scheme and the finite volume element (FVE) scheme, however, positivity is not preserved in small time or time-step for all three schemes, and positivity may reappear after a positivity threshold. Especially, for lumped mass method positivity is preserved if and only if the triangulation is of Delaunay type. Based on energy stable, Liao et al. [15,16] considered the second-order and nonuniform adaptive time-stepping maximum-principle preserving scheme for time-fractional Allen-Cahn equations, where the convolution structure of consistency error is used and sharp maximum-norm error estimates with the temporal regularity is proved. Ji et al. [17] provided the fast L1 formula preserving the discrete maximum principle for the time-fractional Allen-Cahn equation with Caputo's derivative, then extended to volume constraint problem [18]. Based on the piecewise linear interpolation, Ji et al. [19] considered an adaptive second-order Crank-Nicolson time-stepping schemes for time-fractional molecular beam epitaxial growth models, and proved the provided schemes preserves the positive semidefinite property of the integral kernel. However, all of these existing preserving positivity schemes are designed on conforming meshes, which cannot be applied directly to non-conforming grids.
It is very necessary to construct monotone diffusion scheme on general non-conforming meshes. For some applications such as inertial confinement fusion and nuclear reactor, it is important to make sure that the heat propagation travels in the right direction, that is heat should flow from higher to lower temperature on distorted, especially on non-conforming meshes. However, to the best of our knowledge, there is no numerical method to consider the positivity preservation for subdiffusion equation on non-conforming distorted meshes in the appearance of hanging nodes. In this paper, the nonlinear preserving positivity finite volume scheme will be used to non-conforming meshes, moreover, a new method will be introduced to eliminate auxiliary vertex unknowns, in particular those defined vertex unknowns at hanging nodes. We find that our new method is adapted to almost arbitrary mesh geometry including non-conforming cells, and is a cell-centered conservative scheme.
The rest of this paper is organized as follows. In Section 2, we first present some definitions and notations for the meshes, and then we consider the construction of the nonlinear finite volume scheme, and introduce a formula of calculating hanging-node unknowns. In section 3, positivity and nonlinear iteration are studied. Existence of a solution for the positivity-preserving scheme is analysed in section 4. Numerical examples are presented to demonstrate the performance of our new methods in section 5.

Construction of a fully positivity-preserving scheme
Firstly, we present some definitions and notations for the meshes. Set (M , E , P) be a mesh partition of Ω , where -M is a set of polygonal cells K,Ω = K∈MK , and ∂ K is the boundary of K.
-E is the set of edges of all the cells, denoted by σ . We denote by E K ⊂ E the set of edges of K ∈ M . Set E = E ext E int , E ext and E int the set of boundary and interior edges, respectively.
The point x K is the so-called collocation point or cell center which associated with K.
We assume that -Each cell K is a star-shaped polygon with respect to its collocation point x K ; if σ ∈ E ext , the perpendicular line of σ , which starts from x K , intersects on the boundary of Ω , not on any interior edge.
Next, we present more notations (see Figure 2). Let m(K) be the area of cell K, h = (sup K∈M m(K)) 1/2 . Denote the cell vertexes by P 1 , P 2 , P 3 , · · · and A, B. Set n Kσ (resp. n Lσ ) be the unit outer normal on the edge σ of cell K (resp. L), and n Kσ = −n Lσ for σ = K|L = AB. Let t KP i and t LP i be the unit tangential vectors on the line KP j and LP j , j = 1, 2, · · · . The distance between A and B is denoted by |AB|. The unknown u defined at K is denoted by u K . Now we describe the construction of a fully positivity-preserving scheme. Let N be a positive integer. Firstly, by integrating (1.2) on the cell K, we have by using the divergence formula to equation above, we have that is and let be the continuous subdiffusion flux going out of K through on the edge σ .
To get a fully discrete scheme, we first discretise (2.5) in time. Taking into account this initial singularity of the solution of problem (1.2)-(1.4), we use a graded temporal meshes. Set t n = (n/N) r T for n = 0, 1 · · · , N, where r ≥ 1 is the constant mesh grading. Let τ n = t n − t n−1 , I n = (t n−1 ,t n ) for n = 1 · · · , N, and τ = max 1≤k≤N τ k .
For each α k and n ≥ 1, we approximate Caputo time fractional derivative D α k t u(x,t n ) of (2.5) by the classical L1 formula on graded meshes Also, note thatb , by using the mean value theorem, we easily obtaiñ Based on the proposed graded meshes above, for each k = 1, 2, · · · m, there exists a constant C so that The proof of the above results can be found in Lemma 3.1 of Huang and Martin [11].
In order to attain a fully positivity-preserving finite volume method, we also need to approximate the continuous subdiffusion flux (2.6) of the edge σ .
By using Figure 2, the Gauss theorem, and the following well-known formula, where κ T (x,t n ) is the transpose of subdiffusion matrix κ(x,t n ). By the mesh stencil in Figure  2, we can obtain herein, R n K,σ and R n L,σ are the truncation error, which are estimated in Yuan et al. [12,13].
which are discrete normal flux on σ of cell K (resp. L). To obtain a finite volume scheme satisfying the positivity-preserving property and conservation, we use the continuity of the normal flux component on edge σ and introduce the weighted combination ofF n 1 andF n 2 to get: where µ 1 ≥ 0, µ 2 ≥ 0 and µ n 1 + µ n 2 = 1.
By combining the expressionF n 1 andF n 2 , we have To get the two-point flux approximation, we take For the sake of writing, we let Therefore, µ 1 and µ 2 satisfy µ n 1 + µ n 2 = 1, −a n 1 µ n 1 + a n 2 µ n 2 = 0.
By the definition of a n 1 and a n 2 , there is a constant C such that |a n 1 (u) − a n 1 (w)| ≤ C u − w , |a n 2 (u) − a n 2 (w)| ≤ C u − w . Therefore, that is to sayμ n 1 (u) is a continuous function of u. Similarly, we can prove thatμ n 2 (u) is also a continuous function of u.
If σ ∈ E ext (see Fig 2), let Then, we obtain Therefore, by using (2.9), (2.17) and (2.18), we can obtain the positivity-preserving nonlinear finite volume scheme for problem (1.2)-(1.4) as follow where f n K = f (K,t n ). Remark: The reason the scheme is nonlinear is that these coefficients A n K,σ and A n L,σ depend on cell-vertex unknowns.
To obtain a central type scheme, next, we will discuss how to eliminate vertex unknowns in the flux expression (2.17). For cell-vertex unknowns at non-hanging nodes, we use the methods in [13] to approximate. However, for hanging nodes the method in reference [13] fails, especially for "T" style cell-vertex, that is, the hanging node lie a straight segment of a cell-edge. Now, we focus on how to eliminate the auxiliary unknowns at hanging nodes. For convenience, we first introduce a local discrete stencil, which cell-vertex A is a hanging node (see Figure 3).
(2.32) Denote U n P = (u n P 1 , u n P 2 , u n P 3 ), U n K = (u n K 1 , u n K 2 , u n K 3 ), and Then, by solving the system the vertex unknowns u n p 2 can be written as the convex combination of cell-central unknowns. That is u n p 2 = ν 1 u n K 1 + ν 2 u n K 2 + ν 3 u n K 3 . (2.33)

Positivity and nonlinear iteration
Denote U n = (U n K ) K∈M be the vector discrete cell-centered unknowns at the n-th time level. A(U n ) be the matrix corresponding to the spacial discretizaiton. f n = ( f n K ) K∈M , U 0 = (v(K)) K∈M and ψ n = (ψ(K,t n )) K∈M are three vectors corresponding to the source term, initial condition and boundary condition. Then, the positivity-preserving nonlinear finite volume scheme (2.19) can be written in the following matrix form, Let S be the diagonal matrix of the areas of control volumes, L be the diagonal matrix of with a diagonal element of n,α k j , F n is a vectors corresponding to the right-hand-side vector of (3.34). Finally, (3.34) can be further rewritten as For any vector U ∈ R n and U ≥ 0, the matrix A(U) has the following properties: 1. All diagonal entries (A(U)) ii of matrix A(U) are positive; 2. All off-diagonal entries (A(U)) i j (i = j) of A(U) are non-positive; 3. Column sum corresponding to the boundary nodes is positive, and column sum corresponding to the interior nodes is 0.
For the nonlinear algebraic system (3.35), we will use the relaxed version of Picard nonlinear iterations to solve. By choosing a small value E non > 0 and initial vector U 0 0 and repeat for k = 1, 2, . . . , The value of U n is then given by the last U k .
The linear system with non-symmetric matrix (L+A(U k )) is solved by the Bi-Conjugate Gradient Stabilized (BiCGStab) method, and the linear iterations are stop when relative norm of the initial residual becomes smaller then E lin .
Based on the definition of M-matrix [14], for any vector U ≥ 0, A(U) T is an M-matrix. Note that b n,α k n = τ −α k n > 0, then, is also an M-matrix. By using the properties of M-matrix [14], it holds that the matrix (L + A(U n )) −T has non-negative elements. Moreover, from (3.36), we have b n,α k j < 0 for j = 0, 1, · · · , n − 1.
Therefore, under the assumptions f n ≥ 0, ψ n ≥ 0, and U 0 ≥ 0, the linear systems in the above Picard iterations are exactly solved. Moreover, it is clear that all iterative solutions are non-negative vectors.

Existence of a solution for the positivity-preserving scheme
In this section, we prove the existence of a solution for the positivity-preserving nonlinear finite volume scheme (2.19).
Theorem 1 If f n ≥ 0, ψ n ≥ 0 and U 0 ≥ 0, then the nonlinear scheme (3.35) has at least one solution U n .
Proof Define a compact set in R N where C 1 be a non-negative constant.
Define a map φ : C → R N such that Now we need to prove that φ has a fixed point in order to prove that system (4.36) has a solution. In view of the discussion in the above section, we prove that the matrix (L + A(w)) −1 has non-negative elements. Since F n = S f n + ψ n + n−1 ∑ j=0 L j U j ≥ 0, we have φ (w) ≥ 0, f or ∀w ∈ C .
The righthand side of (4.39) is a non-negative constant. We take Then, φ maps C into itself. The set C is a convex compact subset of R N . Since each coefficient of A(w) is a continuous function of w and A(w) → (L + A(w)) −1 is continuous from the set M-matrices to the set of matrices, φ is continuous. Hence, we may apply Brouwer's theorem, which implies that φ has a fixed point in C , i.e., the scheme (3.35) has a solution. The proof of theorem 1 is finished.

Numerical experiments
In this section, we examine numerical performance of the positivity-preserving finite volume method proposed in this paper. The positivity preservation and convergence study of the discrete solution are presented for four subdiffusion problems on both uniform and random distorted quadrilateral meshes. Our numerical experiments cover mainly single-term (ℓ = 1). We use the graded meshes in time with N intervals. The optimal value of the mesh grading parameter r = 2−α 1 α 1 is used in all five examples. By choosing the time step τ and the number of cell so that the error stemming from the temporal and spatial is consistent. In this paper we take E non = 1.0e − 15 and E lin = 1.0e − 20. Moreover, to illustrate the efficiency of our scheme, we apply it to solve an equilibrium radiation sub-diffusion equation.
Here for the Picard nonlinear iterations we take E non = 1.0e − 15, and E lin = 1.0e − 20 for BiCGStab iterations. We use following discrete norms to evaluate approximation errors.
For the solution u, we use the following L 2 -norm and L ∞ -norm For the flux F, we use the following L 2 -norm and L ∞ -norm

A problem with anisotropic diffusion coefficient
where θ = 5π 12 , κ 1 (x, y) = 1 + 2x 2 + y 2 , κ 2 (x, y) = 1 + x 2 + 2y 2 . The solution is chosen to be u(x, y) = t α 1 sin(πx) sin πy.   We take α 1 = 0.5 and test the proposed monotone finite volume scheme on non-conforming rectangular mesh and random meshes shown in Figure 4 and 7 for four levels of grid refinement. This random mesh is generated from a perturbation of mesh in Figure 4. For a vertex P, we denote h P = 1 2 min {σ |P∈σ } |σ |, and the random mesh is defined by x PR = x P + β (R − 0.5)h P , where β ∈ [0, 1] is a parameter, and R is a normalized random variable. In this test, we take β = 0.5. The numerical solutions obtained by our proposed monotone finite volume scheme are shown in Figure 5 and 8 for four levels of grid refinement, and the numerical solution counterparts at the grid point are shown in Figures 6 and 9. In these figures, we observe that the monotone finite volume scheme preserves the positivity of the solution. Numerical errors and convergence rates about solution and flux are obtained by running cases for four levels of grid refinement, staring with the initial grid and refining the grids uniformly on each successive level in Tables 1 and 2.      In the Example 2, the exact solution u(x, y,t) is unknown. We test our scheme on nonconforming rectangular mesh and non-conforming random quadrilateral meshes (see Figure  10-13). We take α 1 = 0.6 and the perturbation parameter β = 0.2. The numerical solutions are shown in Figures 11-12 for non-conforming rectangular meshes and Figures 14-15 for non-conforming random quadrilateral meshes, which demonstrates that our scheme preserves positivity of the solution in this case. Although our method is positive on the nonconforming rectangular meshes and non-conforming random meshes, however, we find that our method does not preserve the maximum at several distorted grid points, such as, the grid We will discuss a finite volume scheme preserving maximum principle in our subsequent work.

A problem with scalar diffusion coefficient
In the following Example 3, we compare preserving positivity obtained by our proposed nonlinear finite volume method with the methods proposed by Jin et al. [10], in which time was discretized by convolution quadrature generated by the backward Euler method and space was approximated by the three methods on different meshes, including the standard FVE method, SG method and LM method. We take the perturbation parameter β = 0.5. The numerical solutions obtained by our proposed scheme on non-conforming random meshes (see Figure 16) with α 1 = 0.8. From the Figures 17-18, we find that our scheme preserves the positivity of the continuous solution on non-conforming random meshes. However, the numerical solutions obtained by the Jin et al. [10] can produce negative values when the time step is less than some positivity thresholds.   Fig. 18 The projection of the solution at the grid point in the y direction on non-conforming random meshes for Example 3.

Conclusion
In this paper, we construct a nonlinear positivity-preserving finite volume scheme for subdiffusion equations on dirtorted non-conforming quadrilateral meshes, where a new method of eliminating vertex unknowns is constructed. Our method is adapted to subdiffusion equations with anisotropic coefficients on general non-conforming meshes with hanging nodes. Besides, the proofs of the positivity and the existence of a solution for this nonlinear finite volume scheme are also presented.

Figure 1
Non-conforming mesh   Non-conforming rectangular meshes for Example 1 Figure 5 The projection of the solution in the x direction on non-conforming rectangular meshes for Example 1. The projection of the solution at the grid point in the x direction on non-conforming rectangular meshes for Example 1. Non-conforming random meshes for Example 1 Figure 8 The projection of the solution in the x direction on non-conforming random meshes for Example 1. The projection of the solution at the grid point in the x direction on non-conforming random meshes for Example 1.

Figure 10
Non-conforming random meshes for Example 2 Figure 11 The projection of the solution in the x direction on non-conforming rectangular meshes for Example 2.

Figure 12
The projection of the solution at the grid point in the x direction on non-conforming rectangular meshes for Example 2. Non-conforming random meshes for Example 2 Figure 14 The projection of the solution in the y direction on non-conforming random meshes for Example 2.

Figure 15
The projection of the solution at the grid point in the y direction on non-conforming random meshes for Example 2. Non-conforming random meshes for Example 3 Figure 17 The projection of the solution in the y direction on non-conforming random meshes for Example 3.

Figure 18
The projection of the solution at the grid point in the y direction on non-conforming random meshes for Example 3.