Solving 3-D Gray–Scott Systems with Variable Diffusion Coefﬁcients on Surfaces by Closest Point Method with RBF-FD

: The Gray–Scott (GS) model is a non-linear system of equations generally adopted to describe reaction–diffusion dynamics. In this paper, we discuss a numerical scheme for solving the GS system. The diffusion coefﬁcients of the model are on surfaces and they depend on space and time. In this regard, we ﬁrst adopt an implicit difference stepping method to semi-discretize the model in the time direction. Then, we implement a hybrid advanced meshless method for model discretization. In this way, we solve the GS problem with a radial basis function–ﬁnite difference (RBF-FD) algorithm combined with the closest point method (CPM). Moreover, we design a predictor–corrector algorithm to deal with the non-linear terms of this dynamic. In a practical example, we show the spot and stripe patterns with a given initial condition. Finally, we experimentally prove that the presented method provides beneﬁts in terms of accuracy and performance for the GS system’s numerical solution.


Introduction
Reaction-diffusion systems model many physical phenomena. In fact, such a system is frequently encountered into many branches of physics, such as electromagnetism [1], heat transfer [2], and diffusion [3]. For the problem of fluid mechanics in porous formations (aquifers and petroleum reservoirs), reaction-diffusion equations lend themselves as a powerful diagnostic tool to determine the hydraulic properties of formations. Briefly, from a well (or a battery of wells), the water is pumped/injected over a measured interval. The effect of such a stimulation upon the pressure distribution in the flow domain is recorded, and the hydraulic properties are detected by matching the measured pressure values against the theoretical ones.
In the context of hydrological applications, several analytical studies are focused on solving reaction-diffusion equations (see [4] for a comprehensive collection of the classical solutions). These analytical studies commonly model the formation as homogeneous (or at most as a sequence of homogeneous layers). However, geological formations are de facto heterogeneous with hydraulic conductivity, in particular, varying in the space somewhat largely (see, e.g., [5]). These irregular changes have a tremendous impact upon transport processes taking place in porous formations [6]. The common (and widely accepted) approach to tackle these erratic spatial variations is a stochastic framework that regards the conductivity as a random space function [7], therefore rendering the reaction-diffusion equation(s) stochastic [7]. Consequently, mean values are insufficient to adequately characterize the spatial distribution of the aquifer's parameters (that now become random fields) since one is required to quantify the uncertainty. Hence, uncertainty quantification has witnessed developments in response to the promise of achieving reliable predictions. In this research field, the integration of ideas from multiple disciplines to effectively design actions is being adopted by engineers, mathematicians, physicists. Uncertainty of any random field Σ can be quantified by deriving deterministic equations governing the probability density function of Σ or, alternatively, by computing moments of Σ [8]. Although quantifying the uncertainty of such random fields is a mathematically challenging problem, it has been scarcely considered in the past (partly due to technical difficulties).
In this paper, we consider one of the most interesting models in engineering, physics, and chemistry, the Gray-Scott (GS) system, which describes the following chemical reaction: where U, V, and P are chemicals. The corresponding GS surface model has the form where u(x, t) and v(x, t) are the concentrations of chemical species; D u (x, t) and D v (x, t) are the variable diffusion coefficients of the chemicals species U and V, respectively; K is the conversion rate from U to P; and F is the in-flow rate of u from the outside; therefore (F + K), is the removal rate of v from the reaction field; ∇ Γ denotes the curve gradient operator. We aim to design and implement a numerical method to solve (2) using a radial basis function combined with finite differences method (RBF-FD). This represents the first (but most important) step toward uncertainty quantification of the random fields uv. However, in this work, a more general novel case is considered in which the system coefficients depend on space and time. For this reason, first, an implicit difference stepping method is adopted to semi-discretize the model in the time direction. Then, a hybrid advanced meshless approach is considered to fully discretize the model. In this way, an effective sparse meshless method based on radial basis function-finite difference (RBF-FD) technique combined with the closest point method (CPM), as an efficient embedding procedure for solving PDEs on surfaces, is proposed. Here, we suggest a predictor-corrector algorithm to deal with the non-linear terms of the dynamics. Finally, we prove through experiments that RFB-FD provides benefits in terms of accuracy and performance. The paper is organized as follows: Section 2 provides a review of related works; Section 3 describes the numerical discretization procedure. In Section 4, some numerical experiments are considered.

Related Numerical Works
The so-called Turing models are coupled partial differential equations describing the reaction and diffusion behaviour of chemicals. The fundamental concept introduced by Turing was that diffusion favours model instability; as a consequence, these models describe the arising of spatially heterogeneous configurations [9]. These formations are called Turing patterns and were firstly observed in chemical experiments; then, they attracted the attention of mathematical biologists due to their ability to imitate biological both regular patterns, such as the stripes of a zebra or the spots of a cheetah, and more irregular patterns, such as those of leopards and giraffes. Moreover, branching in plants and spot patterns of vegetation in the deserted area are also described. As for vegetation dynamics, results concerning the numerical counterpart of the Turing instability have been investigated [10]. In the Gray-Scott model, the reactions create patterns that successfully describe formations in living things. Data on the intriguing history of the Gray-Scott reaction-diffusion model and its applications to the real world are freely available at Robert Munafo's web page [11]. The main results concerning evidence of the Gray-Scott equations documented in the following: Turing 1952 (embryo gastrulation; multiple spots in a 1-D system), Bard and Lauder 1974 (leaves, hair follicles), Bard, 1981 (spots on deer and giraffe), Murray 1981 (butterfly wings), Meinhardt 1982 (stripes, veins on leaf), Young 1984 (development of eye), Meinhardt and Klinger 1987 (mollusk shells), Turk 1991 (leopard, jaguar, zebra). Reaction and diffusion of chemical species can produce a variety of patterns. These PDE systems model pattern formation on curved surfaces. Here, we present a numerical study for studying Turing patterns [12]. In our approach, we adopt a meshless scheme based on radial basis function-finite difference (RBF-FD). Moreover, we use the closest point discretization to deal with the computational complexity of the problem. Radial basis functions (RBFs) are a fundamental mesh-free method to numerically solve PDEs on irregular domains (see [13][14][15] for the global collocation approach). This is due to the versatility of scattered data interpolation techniques used in several applications, e.g., surface reconstruction, image restoration, and in painting, meshless/Lagrangian methods for fluid dynamics [16]. The numerical solution of elliptic partial differential equations by a global collocation approach based on RBF refers to a strong-form solution in the PDE literature [13].
The main drawback for the global approach, although spectrally accurate, consists of solving large, ill-conditioned, dense linear systems, and many attempts are known to deal with it [17,18]. In some cases, local methods are preferred, giving up spectral accuracy for a sparse, better-conditioned linear system and more flexibility for handling non-linearities. The comprehensive literature concerns local RBF schemes by partitioning the domain, referred to as partition of unity (PU) [19][20][21].
An open issue in RBF theory concerns the choice of the shape parameter. The study of this parameter is a crucial topic to guarantee the stability of numerical methods. To overcome the problem, hybrid approaches called variably scaled kernels (VSKs) have been recently introduced [22].
In [23], the authors propose a different local approach based on a generalization of the classical finite difference (FD) method. The FD weights are computed using polynomial interpolation [24,25]. Finite difference methods are applied to solve PDEs on triangulated surfaces [26]. Some numerical issues are involved in the geometric quantities calculation, including values of the normal vector or the curvature of a surface [27]. For solving PDEs on surfaces with complex geometries [28], a level set representation of surfaces and the use of smart projection operators are generally adopted. RBF approaches are also discussed in the literature for the sphere [29] and to deal with embedded surfaces [30,31].
Regarding the level set representation of a surface, the closest point represents a surface, whereby grid nodes store the closest point in Euclidean distance to the surface, successfully used in the computation of diffusion-generated motions of curves on surfaces [32]. The closest point operator is generally designed to solve PDEs on surfaces in embedding space to extend the problem from the surface to the surrounding space. This method allows for boundary conditions at surface boundaries and immediately generalizes beyond surfaces embedded in R 3 to objects of any dimension embedded in any R n [33].
The CPM is generally designed for problems posed on static surfaces [33], using standard finite difference schemes. The classical formulation adopts explicit time stepping, and the stability results are given in the surface case (see [34] for details). Finally, the solution of PDEs on moving surfaces is also of considerable interest. In this paper, we solve the GS system on surfaces by proposing an implicit formulation in the time direction. Moreover, the RBF-FD [35,36] is considered to discretize the time-independent problem in the spatial direction.

Numerical Discretization Scheme
Here, an accurate implicit formulation is proposed to discretize the problem (2) in the time direction. For this purpose, firstly, the time interval [0, T] is uniformly decomposed , where t j = jτ, j = 0, · · · , M and τ = T/M is time step size. Then, by substituting t = t n+1 in the system (2), the following relations are obtained: with x ∈ Γ ⊂ R 3 and t ∈ [0, T]. The time integer derivative can be discretized at two sequential time levels, n + 1 and n, as follows: where u n = u(x, t n ) and v n = v(x, t n ) . In addition, R n+1 denotes the truncation error that is bounded by |R n+1 | < Cτ. By substituting Equation (4) in (3), the following relations at the (n + 1)−th time level for n = 0, 1, 2, · · · , M − 1 are obtained: By rearranging the above relations, the following equations result in the following: Now, we consider the closest point method for solving PDEs on surfaces using the RBF-FD to discretize the time-independent problems (6) in the spatial direction. In the following subsection, we discuss space decomposition schemes.

Closest Point Method
CPM is an embedding approach for numerically solving partial differential equations on surfaces. It considers a limited, narrow band surrounding the surface Γ. Then, the closest point function maps each point of the embedding band to its closest point on the surface. CPM applies different standard numerical techniques to solve the embedding PDEs. Since it preserves the solution constant along with normal directions to the surface, the gradient and Laplacian over the embedding space all retrieve their inherent surface properties when restricted to the surface. The obtained embedding solutions are equal to the original PDEs' solutions on the surface. In this work, an efficient local meshless method based on RBFs is combined with CPM.
We consider the surface Γ ⊂ R d and ξ to be a point in the embedding space Ω ⊂ R d . The closest point function for Γ is defined by which identifies the closest point of ξ on the surface Γ.
In a neighborhood of the surface, CP Γ is C l -smooth for a C l+1 -smooth surface Γ. Moreover, the embedding space Ω is considered a tubular neighborhood of surface Γ.
To prevent the entry discontinuities into CP Γ , it should be assumed that the radius of the computational tube, Ω, has been chosen such that it consists only of points within a distance κ −1 ∞ of the surface Γ, where κ ∞ is an upper bound on the curvatures of Γ [37]. Regarding the CPM, solving the surface PDE on the surface is converted to solving a proper embedding PDE on the computational domain; the surface derivatives should be replaced with Cartesian derivative and closest points operators according to two principles [33]:

1.
Suppose V to be a function defined on R d that is constant along normal directions of Γ. Then, at the surface, ∇ Γ V = ∇V.

2.
Suppose V to be a vector field on R d that is tangent to Γ and also tangent at all surfaces displaced by a fixed distance from Γ. Then, at the surface, The first principle states that if a function is constant in the normal direction, it only varies along the surface. Moreover, the second principle explains that a flux that is directed everywhere along the surface can only spread out within the surface directions.
If u is a function defined on the surface, u(CP) is constant along the directions normal to the surface; thus, the first principle implies Since u(CP) is tangent to the distance function's level-sets, applying the second principle yields the surface Laplacian, which is known as the Laplace-Beltrami operator in the following form: More generally, we consider the surface diffusion operator ∇ Γ · (A(x)∇ Γ ), where the diffusion coefficient is not constant and depends on position. In this case, A(CP)∇(u(CP)) is tangent to the level surfaces, which implies By combining the two principles, all surface differential operators in the governing problem can be replaced with the corresponding Cartesian differential operators. Therefore, by using the relation (10), the surface Equations (6) are changed to the embedding equations, which depend only on Cartesian derivatives as follows: where the constants along the normal direction extension of u and v are denoted bŷ u : Ω → R andv : Ω → R, respectively. The above relations are equivalent to where G n+1 =û n+1 (v n+1 ) 2 . Finally, the points are on the tubular computational domain, and the corresponding closest point on the surface forms the closest point representation. Therefore, a suitable platform was provided for performing a suitable numerical method in order to solve the embedding problem.

CPM Based on RBF-FD Technique
The radial basis function-finite difference (RBF-FD) method is a hybrid and advanced procedure in numerical analysis constructed by combining the beneficent characteristics of the radial basis function (RBF) and easy implementation of finite difference [35,36]. The main factor behind the RBF-FD method's extension is to reduce the computational cost of global methods. Particularly, reducing computational cost in the face of the big problems is essential.
In this section, the RBF-FD method on the closest point representation is introduced. For this purpose, corresponding to each point of the surface Γ, a local domain containing n c nearest points of the computational domain Ω is considered. The group of points in the local domain is called a stencil. Although the problem raised in this work is threedimensional, a schematic view in two dimensions is shown to better understand the RBF-FD method combined with CPM. Figure 1 is a schematic demonstration of the circular curve in the two-dimensional perspective. This figure shows an example of an RBF-FD stencil that consists of the m = 14 closest grid points to a particular surface point. The aim is to approximate a linear operator L at the point x j on the surface Γ using a linear combination of the function value u j in the corresponding stencil. In this work, the linear Cartesian operator L can be ∇ 2 , ∇, or the identity operator such that where x j = CP Γ (ξ j ), ξ j ∈ Ω and the surface operator L Γ can be demonstrate the ∇ 2 Γ , ∇ Γ , or identity operator. Here, to solve the time discretization embedding Equation (12), we could apply CPM based on RBF-FD. For this purpose, for each scattered point In the RBF-FD procedure, the operator L is approximated at the evaluation point x e by a linear weighted combination of the functionsû j =û(ξ j ) at n c nearest neighboring points Ξ j . Therefore, we aim to identify weights {ω j } n c j=1 as follows: The combined RBF and polynomial interpolant is assumed to be a linear combination of the radial and polynomial functions; thus ,it takes the form where φ(.) is a radial basis function, p i (.) is the i−th monomial in the polynomial basis function in the space coordinates x T = [x, y, z], and n p is the number of monomials. To solve for the unknowns, we force the interpolant to match the data at the point locationŝ where φ(.) is a radial basis function, p i (.) is the i−th monomial in the polynomial basis function in the space coordinates x T = [x, y, z], and n p is the number of monomials. Moreover, to find an unique result for unknown coefficients, the following n c constraint orthogonality conditions are imposed Equations (15) and (16) are arranged in a linear system of n c + n p equations for the as follows: where the coefficient matrix M is a non-singular (n c + n p ) × (n c + n p ) matrix [16,38], in which and 0 is the n p × n p zero matrix. Therefore, solving the linear system (18) yields the unknown coefficients as follows: Evaluating the operator L of the interpolant at x e and using relation (20) give = Lφ x e − ξ 1 · · · Lφ x e − ξ n c Lp 1 (x e ) · · · Lp n p (x e ) Once multiplied, B T M −1 gives the weights vectorW eL = [ω 1 , · · · ω n c , ωn c + 1, · · · ω n c +n p ] for approximating Lû(x c ). Therefore, the solution vectorW jL in relation (14) could be determined in the following form: Note that in the vectorW eL , the entries ω n c +1 , · · · , ω n c +n p should be scrapped. Hence, the first n c components of vectorW eL are constructed of the weight vector We L = [ω 1 , ω 2 , · · · , ω n c ] T corresponding to stencil points Ξ j .
Therefore, for evaluation point x e on the surface, the local approximation (14) in the related stencil can be written in novel notation as follows: Then, by computing W jL for all closest points {x j } N j=1 , the global sparse N × N matrix W L = {W 1L W 2L · · · W NL } T could be assembled such that For computing the numerical solutions of the problem (12), each differential operator is replaced by the relation (26) in the following form: Therefore, relation (27) can be rewritten as follows: where the N × N matrices A u , A v , B u , and B v are defined in the following forms: In numerical implementation to face the non-linear term G, a predictor-corrector Algorithm 1 is used.

end if end while
Moreover, for implementing the procedure, the thin plate spline (TPS) RBFs φ(r) = r 2p ln(r), p = 1, 2, · · · , and generalized multiquadric (GMQ) RBFs are applied. These radial basis functions depend on scatter points on the computational domain, such as x T = [x, y, z] through the radial variable r = (x − x i ) 2 + (y − y i ) 2 + (z − z i ) 2 , which is the Euclidean distance between the point of interest x and a computational point x i in three-dimensional space. The TPS belongs to C 2p−1 and has strictly conditionally positive definite radial functions of order p + 1 [16,38]. Moreover, The GMQ belongs to C ∞ , and for 0 < q < 1, it has a strictly, conditionally positive definite order of one [38]. Moreover, for q = 0.5 we have the standard multiquadric (MQ) RBFs. The accuracy and stability of the methods based on the GMQ radial basis functions depend on the shape parameter c. Furthermore, a set of three-dimensional quadratic monomial basis functions {1, x, y, z, x 2 , y 2 , z 2 , xy, xz, yz, xyz} is used to implement the procedure. It is noteworthy that adding polynomials and constructing the augmented RBFs, which guarantee the non-singularity and uniquely solvability of the linear system when using conditionally positive RBFs, could improve the accuracy of the numerical results [39]. Therefore, adding polynomials could improve the stability of the procedure using RBFs. Furthermore, the shape parameters' sensitivity is reduced by augmenting the polynomials into RBFs, and the range of selection of the shape parameters is extended [40].
To verify the numerical stability of the proposed method, the noise effects on error estimates are investigated. For this purpose, we suppose that the initial solutions u 0 and v 0 in the numerical process are perturbed toũ 0 = (1 + δ)u 0 andṽ 0 = (1 + δ)v 0 , respectively. Therefore, the influence of noise δ on error estimates is studied for perturbation solutions u n andṽ n at T = 1.

Numerical Results
In this section, the proposed technique's performance and accuracy in dealing with the GS system with variable coefficients (2) are verified. The suggested method is implemented to solve four test problems. In the numerical implementation, a tubular-embedding computational domain with a radius of r Ω and N uniformly distributed nodes surrounding the intended surface is considered. Furthermore, the shape parameters' sensitivity is reduced by augmenting the polynomials into RBFs, and the range of selection of the shape parameters is extended [40]. Furthermore, The TPS radial basis function (30) with p = 2 and the MQ radial basis function combined with three-dimensional quadratic monomials basis functions are used to construct shape functions. The root mean square error ( r ) and maximum absolute error ( ∞ ) are considered to measure the accuracy of the following method: which is applied to make comparisons, where w e (x, t) and w a (x, t) demonstrate the exact and approximate solutions, respectively. Furthermore, to verify convergence rates of the presented time discretization scheme, the following rate is calculated:

Example 1.
As the first example, we consider the following problem on a sphere with radius of one and a center of zero: The exact solutions of Equation (1) are u(x, y, t) = 1 2 sin(t) cos(πx) cos(πy) cos(πz), v(x, y, t) = 1 4

sin(t) cos(πx) cos(πy) cos(πz).
The proposed numerical technique is applied to solve the problem. The resulting error estimates and CPU times for various values of N by letting τ = 0.5(dx) 2 for two different RBFs are reported in Table 1. As observed, by increasing the number of computational points, the suggested numerical method's accuracy is improved. Table 2 demonstrates the estimated errors and computational temporal convergence rate R in terms of time step sizes τ. The results are obtained by introducing stencil size n s = 33 and N = 2670 data points in the TPS and MQ cases. The shape parameter for the MQ radial basis function is c = 0.35. The given results show the first-order convergence of the proposed temporal discretization scheme. Furthermore, the suggested method's stability is reported in Table 3. This table shows the impression of the noise (δ) on absolute computational errors and indicates that the proposed method has a reasonable and stable behaviour against the enforcing noise. Moreover, the graphs of exact and computed solutions related to u and v by taking N = 4416 , n s = 33, and τ = 0.5(dx) 2 for TPS radial basis functions are shown in Figure 2. Furthermore, the coefficient matrix's sparsity pattern for two types of RBFs is plotted in Figure 3.
The presented computational technique is used to solve this example. Several reports are given to verify the accuracy and efficiency of the suggested procedure. Table 4 reports the absolute and RMS errors and CPU times of the various values of N by taking τ = 0.5(dx) 2 for the two kinds of RBFs tat are presented. This table shows the accuracy and convergence of the method by increasing the number of computational points. In Table 5, the absolute and RMS errors and computational temporal convergence rate R in terms of time step sizes τ are demonstrated. The results are given by taking stencil size n s = 35 and N = 4509 data points for the two kinds of RBFs. The obtained numerical results illustrate the first-order convergence of the proposed temporal discretization scheme. Moreover, the proposed numerical method's stability affected by the noise (δ) on absolute computational errors is reported in Table 6. In addition, the graphs of exact and computed solutions related to u and v are shown in Figure 4. These results are achieved by letting N = 4509, n s = 35, and τ = 0.5(dx) 2 for TPS radial basis function. The sparsity pattern of the coefficient matrix for two types of RBFs is plotted in Figure 5.     The suggested numerical method is applied to solve this example. The computational results show the accuracy and efficiency of the proposed technique. The accuracy and convergence of the numerical technique in terms of the number of computational points N are demonstrated in Table 7. In this table, the error estimates and CPU time are obtained by letting τ = 0.5(dx) 2 for both TPS and MQ radial basis functions. The computational error estimates and convergence rate in terms of time step size are reported in Table 8. The achieved numerical results are computed by letting N = 2670 and n s = 33. The results show that the computational convergence rate R follows the analytical convergence rate O(τ). Furthermore, the computational method's stability is investigated by considering the effect of noise (δ) on absolute computational errors, as shown in Table 9. The graphs of exact and computed solutions related to u and v are demonstrates in Figure 6. All figures are achieved by letting N = 3316 , n s = 33, and τ = 0.5(dx) 2 for the TPS radial basis function. Furthermore, the sparsity patterns of the coefficient matrix related to two kinds of RBFs are plotted in Figure 7.

Example 4.
In this example, we examine several practical problems to show that the GS model's pattern reconstruction property is preserved with variable diffusion coefficients. Considering that in real-world problems, only some conditions such as initial conditions can be determined, the proposed method's effect and performance in solving the GS model were investigated according to the initial conditions. In these cases, only the initial condition is given on the desired surface. The results show the efficiency of the suggested numerical method in detecting spot and stripe patterns based on the initial condition on the surfaces.
In the first case, we consider a spherical surface with the following initial condition: Furthermore, the coefficients are D u (x, y, z, t) = 10 −4 t(x + y + z). 2 and D v (x, y, z, t) = 5 × 10 −5 t(x + y + z) 2 . The parameters are F = 20 and K = 10 −3 . Figure 8 shows the numerical solutions v and u on the spherical surface.
In the second case, we consider a bumpy spherical surface with the following initial condition: In addition, the coefficients are D u (x, y, z, t) = 10 −4 cos(πx) cos(πy) cos(πz) sin(t) and D v (x, y, z, t) = D u (x, y, z, t), and the constant parameters are F = 5 and K = 10 −3 . Figure 9 shows the numerical solutions v and u on the irregular bumpy sphere surface.
In the third case, we consider an ellipsoidal surface with the following initial condition: where v 0 = 1 − u 0 . Additionally, the coefficients are D u (x, y, z, t) = 10 −4 sin(πx) sin(πy) sin(πz) exp(t) and D v (x, y, z, t) = D u (x, y, z, t), and the parameters are F = 20 and K = 10 −3 . Figure 10 shows the approximate solutions v and u on the surface. As can be seen, in Figures 8-10, each pair of shapes corresponds to the variables u and v. Therefore, the inverse colours in these figures are due to different initial conditions for the two variables u and v. That is, each shape represents the pattern recognition corresponding to each variable by different initial conditions.

Conclusions
In this work, an advanced and powerful computational technique is used to numerically investigate the three-dimensional Gray-Scott system with variable diffusion coefficients on surfaces. Firstly, the appearing time derivatives are discretized by employing an implicit difference stepping approach. The time-independent discretization problem is solved using a meshfree method based on the combination of RBF-FD and CPM. A meshless local RBF-FD method is used to discretize the governing time-independent problem in the spatial direction. In our implementation, two radial basis functions, thin plate splines and multiquadrics radial basis functions, are used. The presented method is applied to four numerical examples. The presented results, through the tables and figures, show the performance and accuracy of the proposed numerical technique.