Gray-Scott Model: Turing Patterns from a Coupled Reaction-Diffusion System

In this tutorial, we explore some pattern formation from the Gray-Scott model:

\[\begin{equation} \begin{aligned} \pdv{u}{t} &= \varepsilon_1\grad^2u+b(1-u)-uv^2, \\ \pdv{v}{t} &= \varepsilon_2\grad^2 v - dv+uv^2, \end{aligned} \end{equation}\]

where $u$ and $v$ are the concentrations of two chemical species. The initial conditions we use are:

\[\begin{align*} u(x, y, 0) &= 1 -\exp\left[-80\left(x^2 + y^2\right)\right], \\ v(x, y, 0) &= \exp\left[-80\left(x^2+y^2\right)\right]. \end{align*}\]

The domain we use is $[-1, 1]^2$, and we use zero flux boundary conditions.

using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(-1, 1, -1, 1, 200, 200, single_boundary=true)
mesh = FVMGeometry(tri)

FVMGeometry with 40000 control volumes, 79202 triangles, and 119201 edges

bc = (x, y, t, (u, v), p) -> zero(u) * zero(v)
u_BCs = BoundaryConditions(mesh, bc, Neumann)
v_BCs = BoundaryConditions(mesh, bc, Neumann)

BoundaryConditions with 1 boundary condition with type Neumann

ε₁ = 0.00002
ε₂ = 0.00001
b = 0.04
d = 0.1
u_q = (x, y, t, α, β, γ, _ε₁) -> (-α[1] * _ε₁, -β[1] * _ε₁)
v_q = (x, y, t, α, β, γ, _ε₂) -> (-α[2] * _ε₂, -β[2] * _ε₂)
u_S = (x, y, t, (u, v), _b) -> _b * (1 - u) - u * v^2
v_S = (x, y, t, (u, v), _d) -> -_d * v + u * v^2
u_qp = ε₁
v_qp = ε₂
u_Sp = b
v_Sp = d
u_icf = (x, y) -> 1 - exp(-80 * (x^2 + y^2))
v_icf = (x, y) -> exp(-80 * (x^ 2 + y^2))
u_ic = [u_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
v_ic = [v_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
u_prob = FVMProblem(mesh, u_BCs;
    flux_function=u_q, flux_parameters=u_qp,
    source_function=u_S, source_parameters=u_Sp,
    initial_condition=u_ic, final_time=6000.0)
v_prob = FVMProblem(mesh, v_BCs;
    flux_function=v_q, flux_parameters=v_qp,
    source_function=v_S, source_parameters=v_Sp,
    initial_condition=v_ic, final_time=6000.0)
prob = FVMSystem(u_prob, v_prob)

FVMSystem with 2 equations and time span (0.0, 6000.0)

Now that we have our system, we can solve.

using OrdinaryDiffEq, LinearSolve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=10.0, parallel=Val(false))

retcode: Success
Interpolation: 1st order linear
t: 601-element Vector{Float64}:
    0.0
   10.0
   20.0
   30.0
   40.0
    ⋮
 5960.0
 5970.0
 5980.0
 5990.0
 6000.0
u: 601-element Vector{Matrix{Float64}}:
 [1.0 1.0 … 1.0 1.0; 3.257488532207521e-70 1.6133794454928614e-69 … 1.6133794454928614e-69 3.257488532207521e-70]
 [1.0 1.0 … 1.0 1.0; 2.4209733463794378e-67 4.710199528842436e-67 … 4.710199528842678e-67 2.4209733463769578e-67]
 [1.0 1.0 … 1.0 1.0; 2.586382429956201e-65 4.737055348542031e-65 … 4.737055348542284e-65 2.5863824299537447e-65]
 [1.0 1.0 … 1.0 1.0; 1.89453206542487e-63 3.291853021046029e-63 … 3.291853021046226e-63 1.8945320654232093e-63]
 [1.0 1.0 … 1.0 1.0; 7.595315302554674e-62 1.2679633107409087e-61 … 1.2679633107409795e-61 7.595315302548382e-62]
 ⋮
 [0.8918649868335269 0.8881326202595985 … 0.8881326202595577 0.8918649868335151; 0.003954731905333725 0.004613268056650112 … 0.004613268056659179 0.003954731905335481]
 [0.8906490236269553 0.8868790470372127 … 0.8868790470371722 0.8906490236269442; 0.004050427796292943 0.004723604858117508 … 0.0047236048581268415 0.004050427796294711]
 [0.8894650411492144 0.8856582493134005 … 0.8856582493133602 0.8894650411492041; 0.004144073264221219 0.00483156833529601 … 0.004831568335305598 0.00414407326422299]
 [0.8883113209628068 0.884468449874595 … 0.884468449874555 0.8883113209627975; 0.004235739107344492 0.004937242868783922 … 0.00493724286879375 0.004235739107346254]
 [0.8871844913663948 0.8833061550545532 … 0.8833061550545136 0.8871844913663864; 0.00432557149610877 0.005040804341900593 … 0.0050408043419106485 0.004325571496110511]

Here is an animation of the solution, looking only at the $v$ variable.

using CairoMakie
fig = Figure(fontsize=33)
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
tightlimits!(ax)
i = Observable(1)
u = map(i -> reshape(sol.u[i][2, :], 200, 200), i)
x = LinRange(-1, 1, 200)
y = LinRange(-1, 1, 200)
heatmap!(ax, x, y, u, colorrange=(0.0, 0.4))
hidedecorations!(ax)
record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
    framerate=60) do _i
    i[] = _i
end

"/home/runner/work/FiniteVolumeMethod.jl/FiniteVolumeMethod.jl/docs/build/tutorials/../figures/gray_scott_patterns.mp4"

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(-1, 1, -1, 1, 200, 200, single_boundary=true)
mesh = FVMGeometry(tri)

bc = (x, y, t, (u, v), p) -> zero(u) * zero(v)
u_BCs = BoundaryConditions(mesh, bc, Neumann)
v_BCs = BoundaryConditions(mesh, bc, Neumann)

ε₁ = 0.00002
ε₂ = 0.00001
b = 0.04
d = 0.1
u_q = (x, y, t, α, β, γ, _ε₁) -> (-α[1] * _ε₁, -β[1] * _ε₁)
v_q = (x, y, t, α, β, γ, _ε₂) -> (-α[2] * _ε₂, -β[2] * _ε₂)
u_S = (x, y, t, (u, v), _b) -> _b * (1 - u) - u * v^2
v_S = (x, y, t, (u, v), _d) -> -_d * v + u * v^2
u_qp = ε₁
v_qp = ε₂
u_Sp = b
v_Sp = d
u_icf = (x, y) -> 1 - exp(-80 * (x^2 + y^2))
v_icf = (x, y) -> exp(-80 * (x^ 2 + y^2))
u_ic = [u_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
v_ic = [v_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
u_prob = FVMProblem(mesh, u_BCs;
    flux_function=u_q, flux_parameters=u_qp,
    source_function=u_S, source_parameters=u_Sp,
    initial_condition=u_ic, final_time=6000.0)
v_prob = FVMProblem(mesh, v_BCs;
    flux_function=v_q, flux_parameters=v_qp,
    source_function=v_S, source_parameters=v_Sp,
    initial_condition=v_ic, final_time=6000.0)
prob = FVMSystem(u_prob, v_prob)

using OrdinaryDiffEq, LinearSolve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=10.0, parallel=Val(false))

using CairoMakie
fig = Figure(fontsize=33)
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
tightlimits!(ax)
i = Observable(1)
u = map(i -> reshape(sol.u[i][2, :], 200, 200), i)
x = LinRange(-1, 1, 200)
y = LinRange(-1, 1, 200)
heatmap!(ax, x, y, u, colorrange=(0.0, 0.4))
hidedecorations!(ax)
record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
    framerate=60) do _i
    i[] = _i
end

This page was generated using Literate.jl.