A Cartesian grid based tailored finite point method for reaction-diffusion equation on complex domains

Abstract

This paper presents a Cartesian grid based tailored finite point method (TFPM) for singularly perturbed reaction-diffusion equation on complex domains. The method is incorporated with the kernel-free boundary integral algorithm, where the semi-discrete boundary value problems after time integration are reformulated into corresponding Fredholm boundary integral equations (BIEs) of the second kind, however with no algorithmic dependence on the exact analytical expression for the kernels of integrals. The BIEs are iteratively solved by the GMRES method while integral evaluation during each iteration resorts to solving an equivalent interface problem, which in practice is achieved by a series of manipulations in the framework of TFPM including discretization, correction, solution, and interpolation. The proposed method has second-order accuracy for the reaction-diffusion equation as demonstrated by the numerical examples.

Introduction

This paper considers the initial boundary value problems (IBVPs) for the reaction-diffusion equation

$$\frac{\partial \mathbf{\text{u}}}{\partial t}=\mathit{D}\mathrm{\Delta }\mathbf{\text{u}}+\mathit{R}(\mathbf{\text{u}})\text{in}\mathrm{\Omega },$$

with the initial condition

$$\mathbf{\text{u}}(\mathbf{\text{x}},0)={\mathbf{\text{u}}}^0(\mathbf{\text{x}})\text{for}t=0\text{and}\mathbf{\text{x}}\in \mathrm{\Omega },$$

and either the Dirichlet boundary condition

$$\mathbf{\text{u}}(\mathbf{\text{x}},t)={\mathbf{\text{g}}}_D(\mathbf{\text{x}},t)\text{for}t>0\text{and}\mathbf{\text{x}}\in \mathrm{\Gamma },$$

or the Neumann boundary condition

$${\partial }_{\mathbf{\text{n}}}\mathbf{\text{u}}(\mathbf{\text{x}},t)={\mathbf{\text{g}}}_N(\mathbf{\text{x}},t)\text{for}t>0\text{and}\mathbf{\text{x}}\in \mathrm{\Gamma }.$$

Here, Ω is a bounded two-dimensional complex domain whose boundary Γ is at least twice continuously differentiable; $\mathbf{\text{u}}=\mathbf{\text{u}}(\mathbf{\text{x}},t)$ is a vector of concentration variables with $\mathbf{\text{x}}\in {\mathbb{R}}^2$ being a space point; n is the unit outward normal on Γ; ${\mathbf{\text{u}}}^0$, ${\mathbf{\text{g}}}_D$, ${\mathbf{\text{g}}}_N$ are given smooth data; $\mathit{R}(\mathbf{\text{u}})$ describes a local reaction kinetics and D denotes a diagonal diffusion coefficient matrix. The diffusion parameters for each equation are assumed as small positive constants. When a certain ${\mathit{D}}_{ii}\ll 1$ (${\mathit{D}}_{ij}$ denotes the element of D in the i-th row and j-th column), the IBVP becomes a singular perturbation problem whose solution may have rapid transitions in some narrow regions [1]. Problems of this type and its steady state version arise in many areas of physical mechanism [2], [3], chemical model [4] and biological systems [5].

It is well known that classical numerical methods for singular perturbation problem may produce wild oscillations throughout the whole domain [6]. Traditional finite difference method (FDM) is no longer appropriate with uniform meshes [7], [3]. Instead, two kinds of general purpose techniques are proposed as improvement. The one is fitted operator method on either uniform or quasi-uniform meshes [8], [9], [10], [11], [12], [13], where numerical schemes are designed for a modified governing equation with the effect of perturbation parameter ε taken account. In the same framework of operator fitting, an optional approach proceeds by constructing a finite difference scheme with local calculation for each interior mesh point, so that the approximation is exact on some specific linear spaces. The space spanned by polynomials with finite degrees generates the original High-Order Difference approximation with Identity Expansion (HODIE) [14], [15], [16], based on which the exponentially fitted scheme [17], [11] and the operator compact implicit (OCI) scheme [18] also come into being. The HODIE kind methods can give high-order approximation for a broad class of singular perturbation problems in space dimensions from one to three [19], [15]. In [20], a Cartesian grid-based immersed finite element method (IFEM) with piecewise bilinear shape functions on multi-interface elements is used to enhance accuracy and stability for elliptic interface problems with multi-domain and triple-junction points.

Another effective treatment is fitted mesh method, i.e., method defined on special meshes condensing the grid points in layer regions. In [21], a non-uniform graded mesh was first set up with the optimum node distribution minimizing the error of the classical finite difference solution. Beyond that, graded mesh has limited application in conjunction with the FDM for reaction-diffusion equations. Relatively, the Shishkin mesh [22], which may be the simplest piecewise uniform mesh, is well received from its beginning [7], [23], [24], [25], [26], [27], [28], [12]. However, most of these works utilize standard or reformed central difference scheme, despite ε-uniform convergence, the results often remain low accuracy. In contrast, high-order compact finite difference (HOCFD) scheme for singularly perturbed reaction-diffusion problems, when well-fitted Shishkin-type meshes are utilized, gives rise to high-order ε-uniformly convergent numerical method [29], [30], [6], [31]. The HOCFD scheme can be constructed in different ways like polynomial interpolation based direct discretization [32], [33] or employment of HODIE fitted operator method [34]. Notice that the latter can be considered as a combination of the fitted mesh method and the fitted operator method, thus benefits from both with strict description of the layer and accurate approximation of the equation.

Throughout the achievements above, FDM for singular perturbed reaction-diffusion equations is confined to those defined in square problem domain or even in one space dimension, and Dirichlet boundary conditions (BCs) are more often encountered than others. More possibilities for domain geometry and BCs are explored by finite element and finite volume method [35], [36], [37], [38], [39], [40]. In this paper, we propose a numerical scheme from another perspective, which is compatible for both small perturbation parameter and complex problem domain but performed on uniform Cartesian grid. The scheme is constructed in a local sense to ensure certain basis functions exactly satisfied, thus somewhat resembles the HODIE. Application of uniform meshes allows the convenience in grid generation and needless of a priori knowledge for the layer. In fact, the essential idea originates from the mesh-free method or the finite point method [41], [42], [43]. However, with the particular properties attached, the method is termed with tailored finite point method (TFPM).

The TFPM was proposed by Han, Huang and Kellogg a decade ago for solving singular perturbation problems of the second-order elliptic type [44]. Subsequent studies on this topic include [45], [46], [47], [48], [49], [50], which have been comprehensively reviewed in [51]. Motivation of this work is to enable the TFPM adapted for the problems in more general complex domain whose boundary is at least twice continuously differentiable. To this end, we follow the process of kernel-free boundary integral (KFBI) method [52], [53], where the governing equation with its BC is integrated into a BIE by Green's third identity. The boundary or volume integral involved is regarded as the solution to corresponding equivalent interface problem, which can be solved by a Cartesian grid based FDM if no singular perturbation exists [54], [55], [56], [57], [58], [59]. Otherwise, interface problem is singularly perturbed so that standard FDM on uniform mesh will lose its efficacy as mentioned before. In [46], a TFPM for singularly perturbed interface problem is given under a simple interface circumstance, whereas the present work proceeds with TFPM but starts from a new point for more complicated cases.

The rest of this paper is organized as follows. The first section gives the time integration method to get the semi-discretized scheme, which is in the form of modified Helmholtz equation with some Dirichlet or Neumann boundary conditions satisfied. The two BVPs are reformulated into the corresponding BIEs in section 3. The ensuing section 4 gives the algorithmic implementation details for integral evaluation with the jump calculation supplemented in Appendix A. Numerical results are presented in section 5 with two elliptic BVP examples and one time-dependent reaction-diffusion model in addition. We make a short summary and brief discussion in section 6. Meanwhile, the TFPM in another form is concisely brought up as an alternative in Appendix B.