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.5776838957747154e-67 5.00362424408089e-67 … 5.003624244081145e-67 2.5776838957720818e-67]
 [1.0 1.0 … 1.0 1.0; 2.86995351512663e-65 5.238449482197094e-65 … 5.238449482197378e-65 2.8699535151239194e-65]
 [1.0 1.0 … 1.0 1.0; 1.770669240540597e-63 3.0844243515241024e-63 … 3.084424351524289e-63 1.7706692405390405e-63]
 [1.0 1.0 … 1.0 1.0; 7.38362131075508e-62 1.2337176370957813e-61 … 1.2337176370958503e-61 7.383621310748956e-62]
 ⋮
 [0.8918439814248814 0.8881110650845869 … 0.8881110650845422 0.8918439814248684; 0.003955466259312813 0.004614079644990367 … 0.004614079645000409 0.003955466259314679]
 [0.8906300037351358 0.8868595112490326 … 0.8868595112489878 0.8906300037351236; 0.004050987598363976 0.00472421924470387 … 0.004724219244714215 0.004050987598365852]
 [0.8894462412294697 0.8856388952186759 … 0.8856388952186312 0.8894462412294583; 0.004144547157895905 0.004832092877430623 … 0.004832092877441258 0.004144547157897783]
 [0.8882909600685053 0.8844474315258799 … 0.8844474315258354 0.8882909600684946; 0.004236222175381564 0.004937790244989291 … 0.004937790245000205 0.004236222175383431]
 [0.8871624264128646 0.8832833347030079 … 0.8832833347029637 0.8871624264128551; 0.004326089888293916 0.005041401049198541 … 0.00504140104920972 0.004326089888295759]

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.