Toys

Reaction Diffusion

Click the canvas to start/pause the simulation.

These simulations are based on reaction-diffusion models of the form

$$\frac{\partial \mathbf{u}}{\partial t} = \gamma \Delta \mathbf{u} + g(\mathbf{u})$$

where $\mathbf{u}$ represents one or more chemical concentrations or phase fields evolving in a spatial domain. Different choices of the reaction term $g$ and the number of fields give rise to qualitatively different patterns.

The Gray-Scott model tracks two chemicals $u$ and $v$ with the stoichiometry

$$\begin{array}{rcl}U + 2V &\rightarrow& 3V,\\ V &\rightarrow& P.\end{array}$$

Species $U$ is fed into the system at rate $f$ and $V$ is removed at kill rate $k$, giving the PDE

$$\left\{\begin{array}{rcl}\dfrac{\partial u}{\partial t} &=& \gamma_u \Delta u - uv^2 + f(1-u),\\[6pt]\dfrac{\partial v}{\partial t} &=& \gamma_v \Delta v + uv^2 - (f+k)v.\end{array}\right.$$

Different $(\gamma_u,\gamma_v,f,k)$ produce spots, stripes, spirals, and coral-like branching -- all called Turing patterns. The display shows $|u - v|$.

Numerical scheme: Forward Euler, 9-point Laplacian, periodic BCs.

References: MIT GrayScott, Karl Sims, Pearson (1993).

Cahn-Hilliard (Cahn and Hilliard, 1958) describes phase separation in a binary mixture through a conserved order parameter $\phi \in [-1,1]$ that labels the two phases:

$$\left\{\begin{array}{rcl} \dfrac{\partial \phi}{\partial t} &=& \Delta\mu,\\[8pt] \mu &=& -\varepsilon^2\Delta\phi + W'(\phi), \end{array}\right. \quad W(\phi) = \frac{(1-\phi^2)^2}{4}.$$

The chemical potential $\mu$ is the variational derivative of the Ginzburg-Landau energy, so Cahn-Hilliard is the $H^{-1}$ gradient flow of $$E[\phi] = \int_\Omega \frac{\varepsilon^2}{2}|\nabla\phi|^2 + W(\phi)\,d\Omega.$$ Unlike Allen-Cahn, the total mass $\int \phi$ is conserved -- the two phases can only rearrange, not appear or disappear. Starting from near-uniform noise (Spinodal preset), the mixture spontaneously separates into interconnected bicontinuous domains, then slowly coarsens by Ostwald ripening: smaller features dissolve and feed the growth of larger ones. In the Off-Critical preset the minority phase instead forms isolated droplets. The Two Bubbles preset uses the same large-and-small tangent-circle geometry as the Allen-Cahn Two Bubbles preset -- compare the two: in Allen-Cahn both circles shrink to nothing (volume is not conserved); in Cahn-Hilliard the small bubble transfers its mass to the large one and the total volume stays fixed. Dark = phase $-1$, bright = phase $+1$.

Numerical scheme: Two-pass forward Euler split -- pass 1 computes $\mu$ from $\phi$ via the chemical potential equation and stores $(\phi, \mu)$ in an intermediate texture; pass 2 samples that texture to advance $\phi$ via $\phi_t = \Delta\mu$. Periodic BCs. Only conditionally stable -- large $\Delta t$ or small $\varepsilon^2$ can cause blow-up.

Research connection: Cahn-Hilliard is a central model in Justin's dissertation on energy-stable schemes for gradient flow problems. The scheme here is forward Euler -- only conditionally stable and not guaranteed to decrease $E$. The dissertation develops semi-implicit schemes that are unconditionally energy-stable. See /p/research/.

References: Cahn and Hilliard (1958), J. W. Cahn (1961), L.-Q. Chen (2002).

The Oseen-Frank model describes a nematic liquid crystal using a director field $\mathbf{n}$ -- a unit vector giving the local average molecular orientation (with $\mathbf{n} \equiv -\mathbf{n}$ since molecules have no head or tail). In the one-constant approximation all elastic modes cost the same energy, giving

$$E[\mathbf{n}] = \frac{K}{2}\int_\Omega |\nabla\mathbf{n}|^2\,d\Omega.$$

The overdamped (no-flow) gradient flow of this energy subject to $|\mathbf{n}|=1$ is

$$\frac{\partial \mathbf{n}}{\partial t} = K\bigl(\nabla^2\mathbf{n} - (\mathbf{n}\cdot\nabla^2\mathbf{n})\,\mathbf{n}\bigr).$$

The projection term $-(\mathbf{n}\cdot\nabla^2\mathbf{n})\,\mathbf{n}$ enforces the unit-norm constraint. The display maps director angle $\theta$ to hue with period $\pi$ (reflecting $\mathbf{n}\equiv-\mathbf{n}$); defect cores appear as points where the hue cycles rapidly.

Topological defects are points where $\mathbf{n}$ cannot be continuously defined. In 2D nematics the fundamental defects have half-integer winding numbers $\pm\tfrac{1}{2}$. Integer ($\pm 1$) defects also appear; the $+1$ radial ("hedgehog" or "aster") and the $-1$ hyperbolic ("saddle") are shown in the Hedgehog preset as an attracting pair -- watch them collide and annihilate. Annihilation shows two $+\tfrac{1}{2}$ and two $-\tfrac{1}{2}$ defects in a symmetric cross pattern; each $+\tfrac{1}{2}$ is nearest to a $-\tfrac{1}{2}$, so two simultaneous annihilation events occur. Coarsening starts from random directors -- many defect pairs nucleate and annihilate until only a few survive.

Numerical scheme: Forward Euler, 9-point Laplacian on each director component, with projection and renormalization each step. Stable for $K \cdot \Delta t \le 0.5$ (stencil max eigenvalue $-2$); current parameters give margin $1 - K\Delta t \cdot 2 = 0.2$.

Research connection: The Oseen-Frank model breaks down inside defect cores where $|\mathbf{n}|\to 0$. The full Landau-de Gennes Q-tensor model regularizes the core by allowing the order parameter to melt there -- that is the model in Justin's dissertation work. See /p/research/.

References: Oseen (1933), Frank (1958), de Gennes and Prost "The Physics of Liquid Crystals" (1993).

Several models here connect directly to Justin's dissertation on energy-stable schemes for gradient flow problems. Allen-Cahn and (upcoming) Cahn-Hilliard are canonical gradient flows; the versions in this toy use forward Euler for simplicity. The energy-stable semi-implicit schemes that unconditionally preserve the free energy appear in the dissertation -- see /p/research/.

Source code: github