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.4325702738964372e-67 4.731756865128093e-67 … 4.73175686512834e-67 2.4325702738939466e-67]
 [1.0 1.0 … 1.0 1.0; 2.6649626397622582e-65 4.8760160587993155e-65 … 4.876016058799581e-65 2.6649626397597336e-65]
 [1.0 1.0 … 1.0 1.0; 1.741316725097496e-63 3.0344727420643513e-63 … 3.034472742064537e-63 1.7413167250959649e-63]
 [1.0 1.0 … 1.0 1.0; 7.3509129131461305e-62 1.2283640430811412e-61 … 1.2283640430812101e-61 7.350912913140033e-62]
 ⋮
 [0.8918491025657457 0.8881162523381297 … 0.8881162523380723 0.8918491025657383; 0.003955601948747031 0.004614252937495002 … 0.0046142529375096325 0.0039556019487484665]
 [0.8906346525997196 0.8868642481189674 … 0.8868642481189092 0.8906346525997128; 0.0040511366616753 0.0047244002668829855 … 0.0047244002668981174 0.004051136661676758]
 [0.8894526278633508 0.8856454772509151 … 0.8856454772508565 0.8894526278633446; 0.004144588572213786 0.004832137901849384 … 0.004832137901865002 0.004144588572215255]
 [0.8883011126227521 0.8844579558940665 … 0.8844579558940076 0.8883011126227468; 0.004236038956508563 0.0049375630113348504 … 0.004937563011350942 0.004236038956510029]
 [0.8871747464537841 0.8832961228782099 … 0.8832961228781507 0.8871747464537798; 0.0043257406427064694 0.005040979593428648 … 0.005040979593445194 0.004325740642707921]

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.