reaction diffusion

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. D_u and D_v 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.024 and k = 0.06 at 10,000 steps:

F = 0.04 and k = 0.06 at 10,000 steps:

F = 0.05 and k = 0.06 at 10,000 steps:

F = 0.03 and k = 0.056 at 10,000 steps:

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.