4.2 A Problem and its Decomposition

We will model a Belusov-Zhabotinsky oscillator in which we've crafted the initial conditions to prevent oscillations, resulting in a travelling wave solution instead.

The system can be thought of as a reactant ( $\phi_{{1}}^{}$ ) turning into a product ( $\phi_{{2}}^{}$ ) in a temperature field ( $\phi_{{0}}^{}$ ). The temperature, reactants, and products diffuse in space, from regions of higher concentration to regions of lower concentration. The initial conditions are such that the reactant is unevenly distributed in space, so it starts diffusing while simultaneously producing products. These products are produced only where the reactant is present, so they too start diffusing the moment they are produced.

The reaction term for the reactant is such that it is always positive (think of the reactant being produced from an infinite substrate), so it does not run out. It also diffuses out, thus forming a traveling wave. The same thing is reflected in the product; however its production gets shut off when it reaches a certain level, vis-a-vis the reactant concentration.

This system can be described by the following PDE on domain $\mathcal {D}$

$${\frac{{\partial \mathbf{\Phi}}}{{\partial t}}} \Phi \nabla \Phi \nabla^{2}_{} \Phi \Phi$$ =

$$\Phi \Phi$$ =

$$\Phi \mbox{on $\partial {\mathcal{D}}$}$$ = (4.1)

This equation consists of two terms (on the right hand side): D , which includes spatial derivatives and R , which is a local term. Each of these terms are dependent on parameters k_D and k_R . Initial and Dirichlet boundary conditions are also prescribed. We further define D and R to be

$$\left(\vphantom{ \begin{array}{c} \epsilon \nabla^2 \phi_0 \ \epsilon \nabla^2 \phi_1 \ \epsilon \nabla^2 \phi_2 \end{array} }\right. \begin{array}{c} \epsilon \nabla^2 \phi_0 \ \epsilon \nabla^2 \phi_1 \ \epsilon \nabla^2 \phi_2 \end{array} \left.\vphantom{ \begin{array}{c} \epsilon \nabla^2 \phi_0 \ \epsilon \nabla^2 \phi_1 \ \epsilon \nabla^2 \phi_2 \end{array} }\right) \left(\vphantom{ \begin{array}{c} 0 \ \frac{1}{\epsilon} \left(... ...\phi_1) (\phi_1 - \phi_{th}) \right) \ (\phi_1 - \phi_2) \end{array} }\right. \begin{array}{c} 0 \ \frac{1}{\epsilon} \left( \phi_1 ( 1 - \phi_1) (\phi_1 - \phi_{th}) \right) \ (\phi_1 - \phi_2) \end{array} \left.\vphantom{ \begin{array}{c} 0 \ \frac{1}{\epsilon} \left(... ...\phi_1) (\phi_1 - \phi_{th}) \right) \ (\phi_1 - \phi_2) \end{array} }\right)$$ (4.2)

where k_D = $\epsilon$ , k_R = {0,$\epsilon$, 1} and $\phi_{{th}}^{}$ = ($\phi_{2}^{}$ + q)/f , q = 0.01, f = 0.3 . $\nabla^{2}_{}$ is the Laplacian operator.

We wish to solve this equation on the 2D unit square defined on 0 $\leq$ x $\leq$ 1, 0 $\leq$ y $\leq$ 1 . q(x) = {0.1, 0.1, 0.1},$\epsilon$ = 10^-3 . The initial condition for $\phi_{0}^{}$,$\phi_{2}^{}$ are 0.1, while for $\phi_{1}^{}$ ,

$$\phi_{1}^{} \left(\vphantom{ \exp( - a(\mathbf{x}- \mathbf{x}_1)^2 ) + \exp( - a (\mathbf{x}- \mathbf{x}_2)^2) }\right. \left.\vphantom{ \exp( - a(\mathbf{x}- \mathbf{x}_1)^2 ) + \exp( - a (\mathbf{x}- \mathbf{x}_2)^2) }\right)$$ (4.3)

where a = (0.03)^-2,x₁ = {0.35, 0.35},x₂ = {0.65, 0.65}

While this equation can be solved in many ways, we will use second-order centered finite-difference stencils for $\nabla^{2}_{}$ , i.e.

$$\nabla^{2}_{} \phi {\frac{{ \phi_{i+1,j} - 2\phi_{ij} + \phi_{i-1,j} }}{{ (\Delta x)^2}}} {\frac{{ \phi_{i,j+1} - 2\phi_{ij} + \phi_{i,j-1} }}{{ (\Delta y)^2}}}$$ (4.4)

where $\phi_{{ij}}^{}$ = $\phi$(i$\Delta$x, j$\Delta$y) . We will use a conventional RK2 time-integrator to solve this set of equations. We will use a N_x x N_y mesh (i.e. there will be N_x x N_y grid cells) to discretize $\mathcal {D}$ and a time step of $\Delta$t ; N_x , N_y and $\Delta$t will be inputs into the code. Note, that N_x , N_y and $\Delta$t are not independent; $\Delta$t $\propto$ min{N_x^-2, N_y^-2} for stability. $\Phi$ will be stored at cell-centers, and each cell will be $\Delta$x x $\Delta$y in size, $\Delta$x = 1/N_x and $\Delta$y = 1/N_y . Fig. 4.1 depicts the mesh.

Figure 4.1: 2D central differences stencil.

Solving this equation requires the following functionalities:

• a Mesh, which discretizes the domain;
• a data object, that contains the dependent variables on the cell-centers of the mesh;
• initial conditions, to impose $\Phi$⁽⁰⁾ ;
• boundary conditions, to impose q(x) ;
• an RK2 integrator, to advance in time;
• physics components, to evaluate D and R ;
• data output components, to store the data and a visualizer, to visualize the data; and
• a driver, to orchestrate the creation of variables, and the time-stepping.

These functionalities, implemented as components, can be used to solve the PDE and form a starting point for a PDE-solving toolkit.