Categories

## Using Gray-Scott Reaction-Diffusion Model to Create Spots (and Stripes) on 3D Spot

By SGI Fellow Natalie Patten

The Gray-Scott reaction-diffusion model consists of two partial differential equations: $$\frac{\partial U}{\partial t}=D_u\nabla^2U-UV^2+F(1-U)$$ $$\frac{\partial V}{\partial t}=D_v\nabla^2V+UV^2-(F+k)V.$$Using a semi-implicit approach, each $$U^{i+1}$$ can be calculated as follows: $$(I -{\Delta}tD_uL)U^{i+1}=U^i + {\Delta}t (F(1-U^i) – U^iV^iV^i)$$. Here, $$I$$ is the identity matrix, $${\Delta}t$$ is change in time, $$D_u$$ is the diffusion rate, $$L$$ is the Laplacian of the mesh, $$F$$ is the feed rate, and $$k$$ is the degrading rate. Similarly, $$V^{i+1}$$ can be calculated as follows: $$(I -{\Delta}tD_vL)V^{i+1}=V^i + {\Delta}t (U^iV^iV^i-(F+k)V^i)$$. Using Eigen‘s SimplicialLDLT solver, which uses Cholesky decomposition, we can solve for each:

Eigen::SparseMatrix<double> temp1 = (I - timeStep * du * L).eval();
Eigen::SparseMatrix<double> temp2 = (I - timeStep * dv * L).eval();

Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver1;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver2;

solver1.compute(temp1);
solver2.compute(temp2);

for (int i = 0; i < numSteps; i++) {
Eigen::VectorXd UVV = U.array() * V.array() * V.array();
Eigen::VectorXd FU = F.array() * U.array();
U = solver1.solve((U + timeStep * (F - FU - UVV)).eval());
Eigen::VectorXd kV = k.array() * V.array();
Eigen::VectorXd FV = F.array() * V.array();
V = solver2.solve((V + timeStep * (UVV - kV - FV)).eval());
}


U was initialized to be 1.0 everywhere, and V was initialized to be 0.0 everywhere except in one small area where it was initialized to be 1.0. Du and Dv were set to be 1.0 and 0.5, respectively. I then modelled this onto a 3D mesh of Spot the cow using Polyscope with a time step of 1 and various F and k values, based on parameters found in Pearson’s “Complex Patterns in a Simple System” (1993). Below is a summary of the initial results. You can see various patterns formed by the Gray-Scott model.

## F = 0.044 and k = 0.063 at 10,000 steps:

Acknowledgments:

• Thank you to my project partners Mariem Khlifi and Gimin Nam for help with (lots of!) debugging.
• Thank you to my project mentor Etienne Vouga and project TA Erick Jimenez Berumen for guidance, debugging, and resources.