Analysis of Gray Scott’s Model Numerically

Ahmed Abdulrahim Ahmed Amin, Daoud Suleiman Mashat^*

Department of Mathematics, Numerical Analysis, King Abdulaziz University, Jeddah, Saudi Arabia.

DOI: 10.4236/ajcm.2021.114018                

Abstract

In this paper, a two-dimensional nonlinear coupled Gray Scott system is simulated with a finite difference scheme and a finite volume technique. Pre and post-processing lead to a new solution called GSmFoam by understanding geometry settings and mesh information. The concentration profile changes over time, as does the intensity of the contour patterns. The OpenFoam solver gives you the confidence to compare the pattern result with efficient numerical algorithms on the Gray Scott model.

Keywords

Fourth Order Compact Scheme, Finite Volume Method, Fully Implicit Scheme, Alternating Direction Implicit (ADI) Scheme,
Gray Scott Solver, OpenFoam

1. Introduction

Chemical physics typically uses reaction-diffusion equations to explain temperature and concentration distributions with some pattern formulations [1] . The diffusion term describes the rate of heat and mass transport, while the reaction term describes the rate of heat and mass creation [1] [2] . Nonlinear concepts are frequently referred to as mass action laws because of their broad application [1] [2] . Similar models can also be employed in other applications where transport phenomena are determined by a random movement and where production or consumption conditions should be taken into consideration [1] [2] [3] . A system of reacting chemicals serves as a classic illustration [3] [4] . Population dynamics uses diffusion and reaction terms to characterize the density of a population. Diffusion describes the random movement of individuals, while reaction describes the reproduction of those individuals. The kind of reactions used may be different depending on the research. For a long time, they were thought of as a result of reacting components or interacting populations, according to Gray Scott [3] [4] [5] [6] . There are still some applications today that use the mass-action law, although more intricate functional forms like the Monod functional response in biochemistry or the Holling types II and III in population dynamics are frequently considered more realistic alternatives. It is common for reaction-diffusion systems to have equilibrium conditions where the reaction terms vanish. If there are multiple such equilibrium points, then a transition between them is likely to occur. Reaction-diffusion waves are responsible for these transitions between states. Examples of such occurrences include the spread of flames, the migration of biological species, and the growth of tumours.

A survey of reaction-diffusion wave applications in biology will be presented here [3] [4] [5] [6] . We shall begin with a brief introduction to the mathematical theory of reaction-diffusion patterns. Algorithms using nonlinear reaction-diffusion (NRD) equations are addressed in depth by Alan Turing in 1952, where linear stability and pattern creation are critical in biochemistry. Alan Turing suggested that a chemical be diffused at a uniform initial condition, and the reaction that results is monitored and explained to explain biological growth patterns. There are a lot of easy-to-study response diffusion systems out there. However, the mathematical modelling that goes along with these models is notoriously difficult to understand. In a liquid medium, reaction-diffusion creates spatial patterns that arise and dissipate quickly. To solve this challenge, a viscous permeable support liquid medium, such as gels, is utilized to slow the formation of patterns, which is essential for vision. The structure would arise if a chemical system in a stable state were disturbed. These pattern-forming systems are interesting because they remain out of steady-state over lengthy periods of time, which is a feature of biological systems [3] [4] [5] [6] . Systems for analyzing reaction-diffusion measure the impact of chemical reactions distributed over space.

In well-stirred environments, such changes produce non-spatial periodic oscillations. Transport phenomena like diffusion, which spread these oscillations across space, make them significantly more fascinating [6] [7] [8] . Structure development in the activator-inhibitor Gierer-Meinhardt model is thought to need the interaction of these antagonistic feedback processes (the synthesis of new molecules and their dispersion) [8] [9] [10] . Filip Buric [7] investigated the patterns generation in different species at different parameters in term of behavior which emerges population composition and formation of their structure [10] [11] [12] . M. Cronhjart observed replicators possess continuous variability but behave as small molecules emphasis being placed on their interactions which are explained in such that the Gray Scott model and its extensions presented here specify a limited availability of substrate used in replication and therefore introduce competition between species in fuel consumption but unlike completely generalized models elaborated by H. Takagi [11] [12] [13] [14] [15] . Also, H. Takagi explained spots patterns in a reaction diffusion system with many chemicals. Jeff S McGough invested cubic auto-catalytic reactions while studying Gray Scott model. Numerically Gray Scott system is solved by Ascher in 1995 by using implicit explicit numerical technique [16] while W. Chen investigated numerical behavior of nonlinear convection diffusion equation by multistep finite element methods in 2001 [15] [16] [17] [18] [19] . Also, Kai Zhang investigation to solve such a complex pattern-oriented system by using second order implicit explicit schemes on Gray Scott model reported in 2008 [15] [16] [17] [18] [19] . Seyed Ali Madani studied the significant effects on two dimensional pattern formations of chemical reactions concerned with diffusion of species like Gray Scott model by using the explicit finite difference method in 2014 [15] [16] [17] [18] [19] .

Dejene Gizaw used FEM to Gray Scott reaction diffusion problem by applying FEniCs Software in 2016 [20] [21] [22] [23] . S. Hasnain and Muhammad Saqib studied the reaction diffusion phenomena in three dimension in order to apply implicit finite difference higher order scheme in 2017. While high efficiency techniques to solve convection diffusion problems is observed in the article by S. Hasnain and Muhammad Saqib. Also, S. Hasnain and Muhammad Saqib explained three dimensional coupled nonlinear reaction diffusion system by compact schemes to get high efficiency of the numerical algorithms which constitute the experimental and simulation oriented study of Gray Scott model [22] - [27] .

2. Gray Scott Model Governing Equations

The Gray Scott model is a two dimensional model of diffusion reaction that can produce interesting patterns. Two species ${\xi }_i$ and ${\xi }_j$ interact in the following ways:

${\xi }_i+2{\xi }_j\to 3{\xi }_j$

and

$\Phi \to {\xi }_i\to \Phi$ .

While ${\xi }_j$ can only be produced by the above reactions involving two values of ${\xi }_j$ and one ${\xi }_i$ , ${\xi }_j$ decays on its own with a rate $\left(f+g\right)$ larger than the one controlling the decay of ${\xi }_i$ .

${\xi }_j\to \Phi$ .

Introduction leads us to study Gray Scott model by simulations oriented environment attached with simple geometry and meshing. Let us consider the mathematical equation oriented Gray Scott model,

$\frac{\partial {\xi }_i}{\partial t}={\mathbb{D}}_{{\xi }_i}\left(\frac{{\partial }^2{\xi }_i}{\partial x^2}+\frac{{\partial }^2{\xi }_i}{\partial y^2}\right)-{\xi }_i{\xi }_j^2+\chi \left(1-{\xi }_i\right),$ (1)

$\frac{\partial {\xi }_j}{\partial t}={\mathbb{D}}_{{\xi }_j}\left(\frac{{\partial }^2{\xi }_j}{\partial x^2}+\frac{{\partial }^2{\xi }_j}{\partial y^2}\right)+{\xi }_i{\xi }_j^2+{\xi }_j\left({\chi }_0+\mathbb{K}\right),$ (2)

Where,

· ${\mathbb{D}}_{{\xi }_i}$ diffusion coefficient for first specie.

· ${\mathbb{D}}_{{\xi }_j}$ diffusion coefficient for second specie.

· The first term ${\mathbb{D}}_{{\xi }_i}\left(\frac{{\partial }^2{\xi }_i}{\partial x^2}+\frac{{\partial }^2{\xi }_i}{\partial y^2}\right)$ , we called as diffusive term which is

concentration coefficient (to be constant). Such term shows that first concentration increase in proportion to operator, we say Laplacian (multidimensional, second order derivative explains variation in the gradient). If the quantity of species is higher in neighboring vicinity, the value of ${\xi }_i$ will increase which constitutes diffusion system known as heat equation [9] [13] [19] [23] [24] [28] .

· The second term represents nonlinear reaction rate which is $-{\xi }_i{\xi }_j^2$ . As there is no constant involved in reaction term but strength of the reaction can be adjustified by concentration constants.

· The third term is $\chi \left(1-{\xi }_i\right)$ , known as replenishment term, completely uses of ${\xi }_i$ give rise to generation of ${\xi }_j$ such that all the chemical utilize eventually to replenish it. The third term in Equation (2) shows that ${\xi }_i$ will be increased by the rate proportional to the difference between its initial level and 1. If no term make interference between two concentrations, 1 will be the highest value for the ${\xi }_i$ [9] [13] [19] [23] [24] [28] .

· $\mathbb{K}$ represents the rate of conversion of ${\xi }_j$ to new.

· ${\chi }_0$ represents the feeds rate.

3. Numerical Methods

Next, we will go over numerical solutions for the two-dimensional non-linear coupled Gray Scott problem using finite difference approximation. To demonstrate the difference and similarity between the finite difference approximation and the finite volume method (FVM), and to compare the outcomes of the two methods, the coupled system approximation uses finite difference and finite volume. To begin, the domain must be discretized (split) into equal pieces or control volumes (it is not necessary for the length of each segment to be equal, but for simplicity, we use equal segments). Consider a finite domain $\omega =\left\{\left(x,y\right)|a<x<b,c<y<d\right\}$ with step sizes in x-direction $h_x=\left(b-a\right)/\left(N-1\right)$ and y-direction $h_y=\left(d-c\right)/M-1$ respectively. Also $\left(x_i\mathrm{,}{y}_j\right)$ , are grid points or Cartesian mesh in x & y directions for which $i=1,2,\cdots ,N-1$ and $j=1,2,\cdots ,M-1$ represent the interior nodes or interior points of mesh. We assume $t_{final}=T_f$ which discretized into finite steps in time t [19] [23] [24] [28] [29] [30] . Methods are categorized into two parts:

· Implicit Finite Difference coupled with ADI scheme.

· Replace ${\xi }_i\to \xi$ and ${\xi }_j\to \eta$ .

· Fully Implicit Finite Volume scheme.

· Replace ${\xi }_i\to \varphi$ and ${\xi }_j\to \phi$ .

3.1. Scheme 01

Scheme procedure can be as follows:

Interior Boundary Points:

$\begin{array}{l}{\xi }_{\left(xx\right)}=\gamma {{\xi }^{″}}_{\left(i-\mathrm{1,}:\right)}+{{\xi }^{″}}_{\left(i\mathrm{,}:\right)}+\gamma {{\xi }^{″}}_{\left(i+\mathrm{1,}:\right)}\\ {\xi }_{\left(yy\right)}=\gamma {{\xi }^{″}}_{\left(:\mathrm{,}j-1\right)}+{{\xi }^{″}}_{\left(:\mathrm{,}j\right)}+\gamma {{\xi }^{″}}_{\left(:\mathrm{,}j+1\right)}\end{array}$ (3)

$\begin{array}{l}{\eta }_{\left(xx\right)}=\gamma {{\eta }^{″}}_{\left(i-\mathrm{1,}:\right)}+{{\eta }^{″}}_{\left(i\mathrm{,}:\right)}+\gamma {{\eta }^{″}}_{\left(i+\mathrm{1,}:\right)}\\ {\eta }_{\left(yy\right)}=\gamma {{\eta }^{″}}_{\left(:\mathrm{,}j-1\right)}+{{\eta }^{″}}_{\left(:\mathrm{,}j\right)}+\gamma {{\eta }^{″}}_{\left(:\mathrm{,}j+1\right)}\end{array}$ (4)

$\begin{array}{l}{\xi }_{\left(xx\right)}=\frac{a_0}{4h^2}\left({\xi }_{\left(i+\mathrm{2,}:\right)}-2{\xi }_{\left(i\mathrm{,}:\right)}+{\xi }_{\left(i-\mathrm{2,}:\right)}\right)+\frac{a_1}{h^2}\left({\xi }_{\left(i+\mathrm{1,}:\right)}-2{\xi }_{\left(i\mathrm{,}:\right)}+{\xi }_{\left(i-\mathrm{1,}:\right)}\right)\\ {\xi }_{\left(yy\right)}=\frac{a_0}{4h^2}\left({\xi }_{\left(:\mathrm{,}j+2\right)}-2{\xi }_{\left(:\mathrm{,}j\right)}+{\xi }_{\left(:\mathrm{,}j-2\right)}\right)+\frac{a_1}{h^2}\left({\xi }_{\left(:\mathrm{,}j+1\right)}-2{\xi }_{\left(:\mathrm{,}j\right)}+{\xi }_{\left(:\mathrm{,}j-1\right)}\right)\end{array}$ (5)

$\begin{array}{l}{\eta }_{\left(xx\right)}=\frac{a_0}{4h^2}\left({\eta }_{\left(i+\mathrm{2,}:\right)}-2{\eta }_{\left(i\mathrm{,}:\right)}+{\eta }_{\left(i-\mathrm{2,}:\right)}\right)+\frac{a_1}{h^2}\left({\eta }_{\left(i+\mathrm{1,}:\right)}-2{\eta }_{\left(i\mathrm{,}:\right)}+{\eta }_{\left(i-\mathrm{1,}:\right)}\right)\\ {\eta }_{\left(yy\right)}=\frac{a_0}{4h^2}\left({\eta }_{\left(:\mathrm{,}j+2\right)}-2{\eta }_{\left(:\mathrm{,}j\right)}+{\eta }_{\left(:\mathrm{,}j-2\right)}\right)+\frac{a_1}{h^2}\left({\eta }_{\left(:\mathrm{,}j+1\right)}-2{\eta }_{\left(:\mathrm{,}j\right)}+{\eta }_{\left(:\mathrm{,}j-1\right)}\right)\end{array}$ (6)

· Such scheme constitutes tridiagonal structure with $a_1=\frac{4}{3}\left(1-\gamma \right)$ and $a_0=\frac{1}{3}\left(-1+10\gamma \right)$ .

· For $\gamma =0$ , scheme can be recognized as fourth order.

1st Point at Left Boundary:

$\begin{array}{l}{{\xi }^{″}}_{\left(\mathrm{1,}:\right)}+\gamma {{\xi }^{″}}_{\left(\mathrm{2,}:\right)}=\frac{c_I{\xi }_{\left(I\mathrm{,}:\right)}}{h^2}\mathrm{,}{{\xi }^{″}}_{\left(:\mathrm{,1}\right)}+\gamma {{\xi }^{″}}_{\left(:\mathrm{,2}\right)}=\frac{c_I{\xi }_{\left(:\mathrm{,}I\right)}}{h^2}\\ {{\eta }^{″}}_{\left(\mathrm{1,}:\right)}+\gamma {{\eta }^{″}}_{\left(\mathrm{2,}:\right)}=\frac{c_I{\eta }_{\left(I\mathrm{,}:\right)}}{h^2}\mathrm{,}{{\eta }^{″}}_{\left(:\mathrm{,1}\right)}+\gamma {{\eta }^{″}}_{\left(:\mathrm{,2}\right)}=\frac{c_I{\eta }_{\left(:\mathrm{,}I\right)}}{h^2}\end{array}$ where $I=\mathrm{1,}\cdots .7$ . (7)

For $\gamma =0$ , the linear system in constant c’s can be written as:

· $c_1=\frac{2077}{157},c_2=\frac{-2943}{110},c_3=\frac{573}{44},c_4=\frac{167}{99},c_5=\frac{-18}{11},c_6=\frac{57}{110},c_7=\frac{-131}{1980}$

2nd Point at Right Boundar:

$\begin{array}{l}\gamma {{\xi }^{″}}_{\left(\mathrm{1,}:\right)}+{{\xi }^{″}}_{\left(\mathrm{2,}:\right)}+\gamma {{\xi }^{″}}_{\left(\mathrm{3,}:\right)}=\frac{c_I{\xi }_{\left(:\mathrm{,}I\right)}}{h^2}\mathrm{,}\gamma {{\xi }^{″}}_{\left(:\mathrm{,1}\right)}+{{\xi }^{″}}_{\left(:\mathrm{,2}\right)}+\gamma {{\xi }^{″}}_{\left(:\mathrm{,3}\right)}=\frac{c_I{\xi }_{\left(:\mathrm{,}I\right)}}{h^2}\mathrm{,}\\ \gamma {{\eta }^{″}}_{\left(\mathrm{1,}:\right)}+{{\eta }^{″}}_{\left(\mathrm{2,}:\right)}+\gamma {{\eta }^{″}}_{\left(\mathrm{3,}:\right)}=\frac{c_I{\eta }_{\left(:\mathrm{,}I\right)}}{h^2}\mathrm{,}\gamma {{\eta }^{″}}_{\left(:\mathrm{,1}\right)}+{{\eta }^{″}}_{\left(:\mathrm{,2}\right)}+\gamma {{\eta }^{″}}_{\left(:\mathrm{,3}\right)}=\frac{c_I{\eta }_{\left(:\mathrm{,}I\right)}}{h^2}\mathrm{,}\end{array}$ where $I=\mathrm{1,}\cdots .7$ .

(8)

For $\gamma =0$ , the linear system in constant c’s can be written as:

· $c_1=\frac{585}{512},c_2=\frac{-141}{64},c_3=\frac{459}{512},c_4=\frac{9}{32},c_5=\frac{-81}{512},c_6=\frac{3}{64},c_7=\frac{-3}{512}$

Top Left Point at Boundary:

$\begin{array}{l}\gamma {{\xi }^{″}}_{\left(N-\mathrm{2,}:\right)}+{{\xi }^{″}}_{\left(N-\mathrm{1,}:\right)}+\gamma {{\xi }^{″}}_{\left(N\mathrm{,}:\right)}=\frac{c_I{\xi }_{\left(N-I\mathrm{,}:\right)}}{h^2}\mathrm{,}\gamma {{\xi }^{″}}_{\left(:\mathrm{,}N-2\right)}+{{\xi }^{″}}_{\left(:\mathrm{,}N-1\right)}+\gamma {{\xi }^{″}}_{\left(:\mathrm{,}N\right)}=\frac{c_I{\xi }_{\left(:\mathrm{,}N-I\right)}}{h^2}\mathrm{,}\\ \gamma {{\eta }^{″}}_{\left(N-\mathrm{2,}:\right)}+{{\eta }^{″}}_{\left(N-\mathrm{1,}:\right)}+\gamma {{\eta }^{″}}_{\left(N\mathrm{,}:\right)}=\frac{c_I{\eta }_{\left(N-I\mathrm{,}:\right)}}{h^2}\mathrm{,}\gamma {{\eta }^{″}}_{\left(:\mathrm{,}N-2\right)}+{{\eta }^{″}}_{\left(:\mathrm{,}N-1\right)}+\gamma {{\eta }^{″}}_{\left(:\mathrm{,}N\right)}=\frac{c_I{\eta }_{\left(:\mathrm{,}N-I\right)}}{h^2}\mathrm{,}\end{array}$ (9)

For $\gamma =0$ , the linear system in constant c’s can be written as:

· $c_1=\frac{585}{512},c_2=\frac{-141}{64},c_3=\frac{459}{512},c_4=\frac{9}{32},c_5=\frac{-81}{512},c_6=\frac{3}{64},c_7=\frac{-3}{512}$

Top Right Point at Boundary:

$\begin{array}{l}\gamma {{\xi }^{″}}_{\left(N-\mathrm{1,}:\right)}+{{\xi }^{″}}_{\left(N\mathrm{,}:\right)}=\frac{c_I{\xi }_{\left(N+1-I\mathrm{,}:\right)}}{h^2}\mathrm{,}\gamma {{\xi }^{″}}_{\left(:\mathrm{,}N-1\right)}+{{\xi }^{″}}_{\left(:\mathrm{,}N\right)}=\frac{c_I{\xi }_{\left(:\mathrm{,}N+1-I\right)}}{h^2}\\ \gamma {{\eta }^{″}}_{\left(N-\mathrm{1,}:\right)}+{{\eta }^{″}}_{\left(N\mathrm{,}:\right)}=\frac{c_I{\eta }_{\left(N+1-I\mathrm{,}:\right)}}{h^2}\mathrm{,}\gamma {{\eta }^{″}}_{\left(:\mathrm{,}N-1\right)}+{{\eta }^{″}}_{\left(:\mathrm{,}N\right)}=\frac{c_I{\eta }_{\left(:\mathrm{,}N+1-I\right)}}{h^2}\end{array}$ where $I=\mathrm{1,}\cdots .7$ . (10)

For $\gamma =0$ , the linear system in constant c’s can be written as:

· $c_1=\frac{2077}{157},c_2=\frac{-2943}{110},c_3=\frac{573}{44},c_4=\frac{167}{99},c_5=\frac{-18}{11},c_6=\frac{57}{110},c_7=\frac{-131}{1980}$

ADI Algorithm Setting:

Combining Equations (3)-(10) according to the following procedure,

$\begin{array}{l}{\xi }_{\left(xx\right)\mathrm{,}j}=\frac{1}{{\left(\Delta x\right)}^2}{C}^{-1}D{\xi }_{\left(\mathrm{,}j\right)}\mathrm{,}\\ {\eta }_{\left(xx\right)\mathrm{,}j}=\frac{1}{{\left(\Delta x\right)}^2}{C}^{-1}D{\eta }_{\left(\mathrm{,}j\right)}\mathrm{,}\end{array}\}$ (11)

$\begin{array}{l}\left(I_x-\frac{1}{2}\frac{\Delta t}{{\left(\Delta x\right)}^2}{C}^{-1}D\right){\xi }_{\left(\mathrm{,}j\right)}^{\star }={\xi }_{\left(\mathrm{,}j\right)}^n+\frac{1}{2}\Delta t\left[{\xi }_{\left(yy\right)\left(\mathrm{,}j\right)}^n+F\right]\mathrm{,}\\ \left(I_x-\frac{1}{2}\frac{\Delta t}{{\left(\Delta x\right)}^2}{A}^{-1}B\right){\xi }_{\left(i\mathrm{,}\right)}^{n+1}={\xi }_{\left(i\mathrm{,}\right)}^{\star }+\frac{1}{2}\Delta t\left[{\xi }_{\left(xx\right)\left(i\mathrm{,}\right)}^{\star }+F\right]\mathrm{,}\end{array}\}\text{where}F=F\left({\xi }_{\left(i\mathrm{,}j\right)}^n\mathrm{,}{\eta }_{\left(i\mathrm{,}j\right)}^n\right)$ (12)

$\begin{array}{l}\left(I_x-\frac{1}{2}\frac{\Delta t}{{\left(\Delta x\right)}^2}{C}^{-1}D\right){\eta }_{\left(\mathrm{,}j\right)}^{\star }={\eta }_{\left(\mathrm{,}j\right)}^n+\frac{1}{2}\Delta t\left[{\eta }_{\left(yy\right)\left(\mathrm{,}j\right)}^n+G\right]\mathrm{,}\\ \left(I_x-\frac{1}{2}\frac{\Delta t}{{\left(\Delta x\right)}^2}{A}^{-1}B\right){\eta }_{\left(i\mathrm{,}\right)}^{n+1}={\eta }_{\left(i\mathrm{,}\right)}^{\star }+\frac{1}{2}\Delta t\left[{\eta }_{\left(xx\right)\left(i\mathrm{,}\right)}^{\star }+G\right]\mathrm{,}\end{array}\}\text{where}G=G\left({\xi }_{\left(i\mathrm{,}j\right)}^n\mathrm{,}{\eta }_{\left(i\mathrm{,}j\right)}^n\right)$ (13)

where F and G are nonlinear terms in Equations (11)-(13). Also, the matrices A, B, C and D are (N × N) sparse with triangular nature.

3.2. Scheme 02

It is possible to separate the scheme implementation process into two components, as follows:

1) Steady state two dimensional system with source.

2) Transient two dimensional system with source.

3.2.1. Algorithm Strategy for Steady Phenomena

· Divide the domain into the finite sized subdomains (finite control volumes) and each subdomain is represented by a finite number of grid points (like Nodes).

· Integrate the governing differential equation (GDE) over each subdomain.

· Consider a profile assumption for the dependent variable (like, interpolation function) to evaluate the above integral which expresses the result as an algebraic quantity at the grid points.

${\displaystyle {\int }_V}{\xi }_{\left(xx\right)}\text{d}x\cdot \text{d}y+{\displaystyle {\int }_V}{\xi }_{\left(yy\right)}\text{d}x\cdot \text{d}y+{\displaystyle {\int }_V}Sourc{e}_{\xi }\text{d}x\cdot \text{d}y=0$ (14)

${\displaystyle {\int }_V}{\eta }_{\left(xx\right)}\text{d}x\cdot \text{d}y+{\displaystyle {\int }_V}{\eta }_{\left(yy\right)}\text{d}x\cdot \text{d}y+{\displaystyle {\int }_V}Sourc{e}_{\eta }\text{d}x\cdot \text{d}y=0$ (15)

${\left[e{\xi }_{\left(x\right)}\right]}_e-{\left[A_w{\xi }_{\left(x\right)}\right]}_w+{\left[A_n{\xi }_{\left(y\right)}\right]}_n-{\left[A_s{\xi }_{\left(y\right)}\right]}_s+\bar{S}\Delta V=\mathrm{0,}$ (16)

${\left[e{\eta }_{\left(x\right)}\right]}_e-{\left[A_w{\eta }_{\left(x\right)}\right]}_w+{\left[A_n{\eta }_{\left(y\right)}\right]}_n-{\left[A_s{\eta }_{\left(y\right)}\right]}_s+\bar{S}\Delta V=\mathrm{0,}$ (17)

Using the assumptions,

$\begin{array}{l}WestFaceFlu{x}_{\varphi }=\frac{A_w\left({\varphi }_P-{\varphi }_W\right)}{{\delta }_{WP}}\mathrm{,}\\ WestFaceFlu{x}_{\phi }=\frac{A_w\left({\phi }_P-{\phi }_W\right)}{{\delta }_{WP}}\mathrm{,}\end{array}$ (18)

$\begin{array}{l}EastFaceFlu{x}_{\varphi }=\frac{A_e\left({\varphi }_E-{\varphi }_P\right)}{{\delta }_{PE}}\mathrm{,}\\ EastFaceFlu{x}_{\phi }=\frac{A_e\left({\phi }_E-{\phi }_P\right)}{{\delta }_{PE}}\mathrm{,}\end{array}$ (19)

$\begin{array}{l}SouthFaceFlu{x}_{\varphi }=\frac{A_s\left({\varphi }_P-{\varphi }_S\right)}{{\delta }_{SP}}\mathrm{,}\\ SouthFaceFlu{x}_{\phi }=\frac{A_s\left({\phi }_P-{\phi }_S\right)}{{\delta }_{SP}}\mathrm{,}\end{array}$ (20)

$\begin{array}{l}NorthFaceFlu{x}_{\varphi }=\frac{A_n\left({\varphi }_N-{\varphi }_P\right)}{{\delta }_{PN}}\mathrm{,}\\ NorthFaceFlu{x}_{\phi }=\frac{A_n\left({\phi }_N-{\phi }_P\right)}{{\delta }_{PN}}\mathrm{.}\end{array}$ (21)

Combining Equations (14)-(21) to get the following system,

$\left(\frac{A_w}{{\delta }_{WP}}+\frac{A_e}{{\delta }_{PE}}+\frac{A_s}{{\delta }_{SP}}+\frac{A_n}{{\delta }_{PN}}-S_P\right){\varphi }_P=\frac{A_w}{{\delta }_{WP}}{\varphi }_W+\frac{A_e}{{\delta }_{PE}}{\varphi }_E+\frac{A_s}{{\delta }_{SP}}{\varphi }_S+\frac{A_n}{{\delta }_{PN}}{\varphi }_N+S_u$ (22)

$\left(\frac{A_w}{{\delta }_{WP}}+\frac{A_e}{{\delta }_{PE}}+\frac{A_s}{{\delta }_{SP}}+\frac{A_n}{{\delta }_{PN}}-S_P\right){\phi }_P=\frac{A_w}{{\delta }_{WP}}{\phi }_W+\frac{A_e}{{\delta }_{PE}}{\phi }_E+\frac{A_s}{{\delta }_{SP}}{\phi }_S+\frac{A_n}{{\delta }_{PN}}{\phi }_N+S_u$ (23)

In a two-dimensional case, the face areas are assumed to be constant and are treated as $A=1$ . The distribution of the $\varphi$ & $\phi$ in a particular two-dimensional scenario is obtained by formulating discretized equations at each grid node of the subdivided domain. To account for boundary conditions, discretized equations must be changed where flux information is available. The boundary-side coefficient is set to zero, and any flux crossing the boundary is added as a new source to any existing $S_u$ and $S_p$ components. The resulting equations are then solved to obtain the $\varphi$ & $\phi$ two-dimensional distribution.

3.2.2. Algorithm Strategy for Transient Phenomena

To determine the right-hand side of the aforementioned Equation (22) and Equation (23), we must make an assumption regarding the change in ${\varphi }_P$ & ${\phi }_P$ , ${\varphi }_E$ & ${\phi }_E$ and ${\varphi }_W$ & ${\phi }_W$ with time. When calculating the time integral, we can utilize the values from the previous step as well as the values from the current step plus the values from the step plus the step after that. The approach can be generalized by using a weighting value between 0 and 1 and writing the integral in the following way:

${\mathbb{I}}_{\varphi }={\displaystyle {\int }_t^{t+\Delta t}}{\varphi }_P\text{d}t$ (24)

$\begin{array}{l}{\mathbb{I}}_{\varphi }=\left[\theta {\varphi }_P+\left(1-\theta \right){\varphi }_P^0\right]\text{d}t\hfill \end{array}$ (25)

${\mathbb{I}}_{\phi }=\left[\theta {\phi }_P+\left(1-\theta \right){\phi }_P^0\right]\text{d}t$ (26)

For example, if $\theta$ is equal to 0, then scheme is explicit while $\theta$ is equal to 1, scheme is fully implicit [17] [23] [26] [31] [32] .

4. Test Problems

We can organize the problem into two parts one dimensional coupled nonlinear system and two dimensional. To demonstrate the effectiveness of the approaches described, we numerically solved the following example.

4.1. One Dimensional Coupled Problem

$\begin{array}{l}{\varphi }_t={\mathbb{D}}_{\varphi }\Delta \varphi -\varphi {\phi }^2+\chi \left(1-\varphi \right)\mathrm{,}t>0\mathrm{<mo>\; \&\; </mo>}x\in \omega \hfill \\ {\phi }_t={\mathbb{D}}_{\phi }\Delta \phi +\varphi {\phi }^2+\phi \left(\chi +\mathbb{K}\right)\mathrm{,,}t>0\mathrm{<mo>\; \&\; </mo>}x\in \omega \hfill \end{array}\}\text{where}{\mathbb{D}}_{\varphi }={\mathbb{D}}_{\phi }$ (27)

with initial and boundary conditions,

$\begin{array}{l}\varphi \left(x,0\right)={\varphi }_{SteadyState}+0.01\sin \left(\frac{\pi x}{L}\right),0\le x\le L\\ \phi \left(x,0\right)={\phi }_{SteadyState}-0.12\sin \left(\frac{\pi x}{L}\right),0\le x\le L\\ \varphi \left(0,t\right)={\varphi }_{SteadyState},\varphi \left(L,t\right)={\varphi }_{SteadyState},t>0,\\ \phi \left(0,t\right)={\phi }_{SteadyState},\phi \left(L,t\right)={\phi }_{SteadyState},t>0,\\ {\varphi }_{SteadyState}=0,{\phi }_{SteadyState}=1,\\ \chi =0.09\&\mathbb{K}=-\mathrm{0.004.}\end{array}\}\text{where}{\mathbb{D}}_{\varphi }={\mathbb{D}}_{\phi }=0.01,$ (28)

4.2. Two Dimensional Coupled Problem

Two dimensional coupled non-linear Gray Scott model system can be written as:

$\begin{array}{l}{\varphi }_t={\mathbb{D}}_{\varphi }\Delta \varphi -\varphi {\phi }^2+\chi \left(1-\varphi \right),t>0\&\left(x,y\right)\in \omega \\ {\phi }_t={\mathbb{D}}_{\phi }\Delta \phi +\varphi {\phi }^2+\phi \left(\chi +\mathbb{K}\right),t>0\&\left(x,y\right)\in \omega \\ {\varphi }_{SteadyState}=0,{\phi }_{SteadyState}=1,\\ \chi =1\&\mathbb{K}=0.\end{array}\}\text{where}{\mathbb{D}}_{\varphi }={\mathbb{D}}_{\phi }=1,$ (29)

4.3. Simulation Setup

The one of the main purposes of this research study is to get good understandings of the simulation of the reaction diffusion model like Gray Scott system including geometry construction, related mesh, boundary and initial conditions, convergence criterion to solver selection, pre & post processing. Cases are settled under block meshes for pre processing with some control parameters as setup simulations in Open Foam (OpenFOAM solvers). Post processing is viewed by using ParaView which is embedded in OpenFoam software [20] [27] [30] [31] [32] .

4.4. Geometry

The geometry consists of the square block with all the boundaries of the square are walls & stationary. Initially, the concentration will be assumed stationary and will be solved on a uniform mesh using the GSmFoam which is designed to solve nonlinear reaction diffusion system oriented problems. The effect of increased mesh resolution and mesh grading at the center of the walls will be investigated. Convert the units to meters as along eight vertices (0 0 0), (1 0 0), (1 1 0), (0 1 0), (0 0 0.1), (1 0 0.1), (1 1 0.1) & (0 1 0.1) with blocks and edges. Boundaries are as top, bottom, left & right with type cyclic while front & back as empty due to two dimensional case. Faces are along (3 7 6 2) as top, (0 4 7 3) bottom, (1 5 4 0) left, (2 6 5 1) right, (0 3 2 1) front and (4 5 6 7) back [20] [27] [30] [31] [32] .

4.5. Mesh

A system in three dimensional cartesian coordinate form can easily solved by the use of OpenFoam solver with all type of geometries as a default parameters

Figure 1. Block type geometry for GS model in two dimensions.

selection. It can be instructed to solve in two dimensions by specifying a special empty boundary condition on boundaries normal to the (3rd) dimension for which no solution is required. Such domain consists of a square of side length $d=0.1\text{m}$ in the xy plane. A uniform mesh of 200 by 200 cells will be used initially. The mesh generator supplied with OpenFOAM (blockMesh) generates meshes from a description specified in an input dictionary located in the system directory for a given case. For the sake of clarity, the file first specifies coordinates of the block vertices. The mesh is generated by running blockMesh on this blockMeshDict file by simplifying in terminal writing blockMesh (Figure 1). The running status of blockMesh is reported in the terminal window. Any mistakes in the blockMeshDict file are picked up by blockMesh and the resulting error message directs the user to the line in the file where the problem occurred [27] [30] [33] [34] .

4.6. Boundary and Initial Conditions

The case is set up for geometry with mesh generation is completed in which initial fields are designed in 0 folder. Such case is set up to start at time $t=0\text{s}$ which contains two sub files named as ${\xi }_i$ and ${\xi }_j$ whose initial and boundary conditions must be specified by considering dimensions $\left[0000000\right]$ with internal field uniform along top, bottom, left & right boundaries with type as cyclic (pattern formations). Also consider front & back boundaries type as empty because of two dimensional coupled system for both components ${\xi }_i$ and ${\xi }_j$ . Selection to boundary conditions are under three principals as internal field is the data which can be uniformly described by a single cell (or at point) or nonuniform for all the other points of the field specified [27] [30] [33] [34] . Boundary conditions are stored in boundary field file which includes data for all the boundary patches which consist of walls. Zero gradient boundary conditions for both components ${\xi }_i$ and ${\xi }_j$ are observed due to normal gradient of the concentration which is zero. Front & back walls are considered empty due to two dimensional case with initial fields being uniform. No slip boundary conditions on fixed stationary walls of the two dimensional case is set up [27] [30] [33] [34] .

4.7. Space Time Discretisation & Solver Selection

We applied finite volume discretisation schemes in the fv (finite volume) Schemes dictionary in the system directory. The specification of the linear equation solvers, tolerances and other algorithm controls is made in the fv (finite volume) Solution dictionary.

4.8. Numerical Schemes

The finite volume schemes (fvSchemes) dictionary in the system directory sets the numerical schemes for terms such as derivatives in equations that appear in applications being run. As standard Gaussian finite volume integration is selected which is the general choice for every user in OpenFoam. Gaussian integration is based on summing values on cell faces which must be interpolated from cell centers. Diffusion term can be solved under Laplacian schemes ( Δ 2 )

while time $\frac{\partial }{\partial t}$ is discretized under time oriented scheme known as Euler

scheme which is a default choice along with Gauss linear for $div\left(\xi \mathrm{,}U\right)$ [27] [30] [33] [34] [35] .

5. Results

Numerical computations have been performed using the uniform grid. For solving problem 01, we set some parameters which can be seen from data in Tables 1-3 by using implicit schemes. Figures 2-5 give results for component ${\xi }_i$ at different time levels which show pattern formulation for the initial level to the final level. Such results indicate for judgement of pattern, the time scale must be larger enough to see in a clear glance. While Figures 6-9 indicate the different settings of parameters that output in different patterns.

Table 1. Shows test problem 01 at time $T_{final}=100$ at different locations. Components are ${\xi }_i$ and ${\xi }_j$ .

Table 2. Shows test problem 01 at time $T_{final}=1000$ for components are ${\xi }_i$ and ${\xi }_j$ for different grid sizes and error during calculations.

Table 3. Shows test problem 01 at different time $T_{final}$ for components are ${\xi }_i$ and ${\xi }_j$ for different grids and CPU usage.

Figure 2. Shows test problem 02 at time $T_{final}=100$ , $\xi =0.022$ and $K=0.051$ .

Figure 3. Shows test problem 02 at time $T_{final}=500$ , $\xi =0.022$ and $K=0.051$ .

Figure 4. Shows test problem 02 at time $T_{final}=1000$ , $\xi =0.022$ and $K=0.051$ .

Figure 5. Shows test problem 02 at time $T_{final}=10000$ , $\xi =0.022$ and $K=0.051$ .

Figure 6. Shows test problem 02 at time $T_{final}=100$ , $\xi =0.01$ and $K=0.041$ .

Figure 7. Shows test problem 02 at time $T_{final}=1000$ , $\xi =0.01$ and $K=0.041$ .

Figure 8. Shows test problem 02 for component $\eta$ at time $T_{final}=1000$ , $\xi =0.01$ and $K=0.45$ .

Figure 9. Shows test problem 02 for component $\eta$ at time $T_{final}=10000$ , $\xi =0.01$ and $K=0.45$ .

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Turing, A.M. (1952) The Chemical Basis of Morphogenesis. Philosophical Transactions of the Royal Society B, 237, 37-72. https://doi.org/10.1098/rstb.1952.0012
[2]	Grzybowski, B.A. (2009) Chemistry in Motion, Reaction-Diffusion Systems for Micro- and Nanotechnology, Chapter 3 & 4. John Wiley & Sons, West Sussex, 72-96.
[3]	Gierer, A. and Meinhardt, H. (1972) A Theory of Biological Pattern Formation. Kybernetik, 12, 30-39. https://doi.org/10.1007/BF00289234
[4]	Chasnov, J.R. (2016) Mathematical Biology, Lecture Notes for MATH 4333, Chapter 4. The Hong Kong University of Science and Technology, Hong Kong, 49-59.
[5]	Meron, E. (1992) Pattern Formation in Excitable Media. Physics Reports, 218, 1-66. https://doi.org/10.1016/0370-1573(92)90098-K
[6]	Murray, J.D. (2003) Mathematical Biology, Volume I, II. 3rd Edition, Springer, New York.
[7]	Buric, F. (2014) Pattern Formation and Chemical Evolution in Extended Gray Scott Models. Master’s Thesis, Chalmers University of Technology, Sweden.
[8]	Cronhjort, M.B. and Blomberg, C. (1994) Hyper-Cycles versus Parasites in a Two Dimensional Partial Differential Equations Model. Journal of Theoretical Biology, 169, 31-49. https://doi.org/10.1006/jtbi.1994.1128
[9]	Cronhjort, M.B. and Blomberg, C. (1997) Cluster Compartmentalization May Provide Resistance to Parasites for Catalytic Networks. Physica D: Nonlinear Phenomena, 101, 289-298. https://doi.org/10.1016/S0167-2789(97)87469-6
[10]	Grzybowski, B.A. (2009) Chemistry in Motion, Reaction-Diffusion Systems for Micro- and Nanotechnology, Chapter 8. John Wiley & Sons, West Sussex, 149-176.
[11]	Takagi, H. and Kaneko, K. (2001) Differentiation and Replication of Spots in a Reaction Diffusion System with Many Chemicals. Europhysics Letters, 56, 145-151. https://doi.org/10.1209/epl/i2001-00500-3
[12]	McGough, J.S. and Riley, K. (2004) Pattern Formation in the Gray-Scott Model. Nonlinear Analysis: Real World Applications, 5, 105-121. https://doi.org/10.1016/S1468-1218(03)00020-8
[13]	Epstein, I.R. (1984) Complex Dynamical Behavior in Simple Chemical Systems. The Journal of Physical Chemistry, 88, 187-198. https://doi.org/10.1021/j150646a007
[14]	Gray, P. and Scott, S.K. (1985) Sustained Oscillations and Other Exotic Patterns of Behavior in Isothermal Reactions. The Journal of Physical Chemistry, 89, 22-32. https://doi.org/10.1021/j100247a009
[15]	Pearson, J.E. (1993) Complex Patterns in a Simple System. Science, 261, 189-192. https://doi.org/10.1126/science.261.5118.189
[16]	Ascher, U.M. and Ruuth, S.J. (1995) Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. SIAM Journal on Numerical Analysis, 32, 797-823. https://doi.org/10.1137/0732037
[17]	Chen, W. (2001) Implicit Explicit Multistep Finite Element Methods for Nonlinear Convection Diffusion Problems. Numerical Methods for Partial Differential Equations, 17, 93-104. https://doi.org/10.1002/1098-2426(200103)17:2<93::AID-NUM1>3.0.CO;2-B
[18]	Zhang, K., Jeff, C., Wong, F. and Zhang, R. (2008) Second-Order Implicit-Explicit Scheme on Gray-Scott Model. Journal of Computational and Applied Mathematics, 213, 559-581. https://doi.org/10.1016/j.cam.2007.01.038
[19]	Madani, S.A. and Shademani, A. (2014) Mathematical Investigation of Two Dimensional Pattern Formation. International Journal of Biological Research, 2, Article ID: 1502. https://doi.org/10.14419/ijbr.v2i1.1502
[20]	Gizaw, D. (2016) Finite Element Method Applied to Gray-Scott Reaction-Diffusion Problem. Using FEniCs Software. Advances in Physics Theories and Applications, 58, 5-8.
[21]	Hasnain, S., Saqib, M. and Mashat, D.S. (2017) Fourth Order Douglas Implicit Scheme for Solving Three Dimension Reaction Diffusion Equation with Non-Linear Source Term. AIP Advances, 7, Article ID: 075011. https://doi.org/10.1063/1.4986322
[22]	Hasnain, S., Saqib, M. and Mashat, D.S. (2017) Highly Efficient Computational Methods for Two Dimensional Coupled Nonlinear Unsteady Convection-Diffusion Problems. IEEE Access, 5, 7139-7148. https://doi.org/10.1109/ACCESS.2017.2699320
[23]	Hasnain, S. and Saqib, M. (2019) Numerical Study to Coupled Three Dimensional Reaction Diffusion System. IEEE Access, 7, 46695-46705. https://doi.org/10.1109/ACCESS.2019.2903977
[24]	Tompkinsa, N., Lia, N., Girabawea, C., Heymanna, M., Bard Ermentrout, G., Epsteind, I.R. and Fraden, S. (2014) Testing Turing’s Theory of Morphogenesis in Chemical Cells. Proceedings of the National Academy of Sciences of the United States of America, 111, 4397-4402. https://doi.org/10.1073/pnas.1322005111
[25]	Toiya, M., Vanag, V.K. and Epstein, I.R. (2008) Diffusively Coupled Chemical Oscillators in a Microfluidic Assembly. Angewandte Chemie International Edition, 47, 7753-7755. https://doi.org/10.1002/anie.200802339
[26]	Vanag, V.K. and Epstein, I.R. (2003) Diffusive Instabilities in Heterogeneous Systems. The Journal of Chemical Physics, 119, 7297-7307. https://doi.org/10.1063/1.1606677
[27]	Greenshields, C.J. (2017) OpenFOAM: The OpenFOAM Foundation User Guide Version 5.0. OpenFOAM Foundation Ltd., U20-U72.
[28]	Litschel, T., Norton, M.M., Tserunyan, V. and Fraden, S. (2018) Engineering Reaction-Diffusion Networks with Properties of Neural Tissue. Lab on a Chip, 5, 2-46. https://doi.org/10.1039/C7LC01187C
[29]	Noyes, R.M., Field, R. and Koros, E. (1972) Oscillations in Chemical Systems. I. Detailed Mechanism in a System Showing Temporal Oscillations. Journal of the American Chemical Society, 94, 1394-1395. https://doi.org/10.1021/ja00759a080
[30]	Cross, M. and Greenside, H. (2009) Pattern Formation and Dynamics in Non Equilibrium Systems. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511627200
[31]	Strogatz, S.H. (1994) Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry, and Engineering. Addison-Wesley Pub., Reading, MA.
[32]	Versteeg, H.K. and Malalasekera, W. (2007) An Introduction to Computational Fluid Dynamics. 2nd Edition, Pearson, London.
[33]	Löhner, R. (2008) Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods. 2nd Edition, John Wiley & Sons, Hoboken. https://doi.org/10.1002/9780470989746
[34]	Li, J.C. (2008) Computational Partial Differential Equations Using MATLAB (Chapman and Hall/CRC Applied Mathematics and Nonlinear Sciences). Chapman and Hall/CRC, Boca Raton.
[35]	Ferziger, J.H. and Peric, M. (2002) Computational Methods for Fluid Dynamics. 3rd Edition, Springer, Berlin.