Keller-Segel Model of Chemotaxis

This tutorial considers the following Keller-Segel chemotaxis model:

\[\begin{equation*} \begin{aligned} \pdv{u}{t} &= \grad^2u - \div \left(\frac{cu}{1+u^2}\grad v\right) + u(1-u), \\ \pdv{v}{t} &= D\grad^2 v + u - av, \end{aligned} \end{equation*}\]

on the square $[0, 100]^2$ with homogeneous Neumann boundary conditions.

This example is intentionally marked as documentation-only metadata in CI: it remains part of the published tutorial set, but because of its runtime and memory cost it is not executed in pull-request CI. The validation manifest records that policy explicitly so the page still has source ownership.

using FiniteVolumeMethod, DelaunayTriangulation

tri = triangulate_rectangle(0, 100, 0, 100, 250, 250, single_boundary = true)
mesh = FVMGeometry(tri)
bc_u = (x, y, t, (u, v), p) -> zero(u)
bc_v = (x, y, t, (u, v), p) -> zero(v)
BCs_u = BoundaryConditions(mesh, bc_u, Neumann)
BCs_v = BoundaryConditions(mesh, bc_v, Neumann)
q_u = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    u = αu * x + βu * y + γu
    ∇u = (αu, βu)
    ∇v = (αv, βv)
    χu = p.c * u / (1 + u^2)
    return χu .* ∇v .- ∇u
end
q_v = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    ∇v = (αv, βv)
    return -p.D .* ∇v
end
S_u = (x, y, t, (u, v), p) -> u * (1 - u)
S_v = (x, y, t, (u, v), p) -> u - p.a * v
q_u_parameters = (c = 4.0,)
q_v_parameters = (D = 1.0,)
S_v_parameters = (a = 0.1,)
u_initial_condition = 0.01rand(DelaunayTriangulation.num_points(tri))
v_initial_condition = zeros(DelaunayTriangulation.num_points(tri))
final_time = 1000.0
u_prob = FVMProblem(
    mesh,
    BCs_u;
    flux_function = q_u,
    flux_parameters = q_u_parameters,
    source_function = S_u,
    initial_condition = u_initial_condition,
    final_time = final_time,
)
v_prob = FVMProblem(
    mesh,
    BCs_v;
    flux_function = q_v,
    flux_parameters = q_v_parameters,
    source_function = S_v,
    source_parameters = S_v_parameters,
    initial_condition = v_initial_condition,
    final_time = final_time,
)
prob = FVMSystem(u_prob, v_prob)

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

To solve and animate the problem:

using OrdinaryDiffEq, Sundials, CairoMakie
sol = solve(prob, CVODE_BDF(linear_solver = :GMRES), saveat = 1.0, parallel = Val(false))
fig = Figure(fontsize = 44)
x = LinRange(0, 100, 250)
y = LinRange(0, 100, 250)
i = Observable(1)
axu = Axis(fig[1, 1], width = 600, height = 600,
    title = map(i -> L"u(x,~ y,~ %$(sol.t[i]))", i), xlabel = L"x", ylabel = L"y")
axv = Axis(fig[1, 2], width = 600, height = 600,
    title = map(i -> L"v(x,~ y,~ %$(sol.t[i]))", i), xlabel = L"x", ylabel = L"y")
u = map(i -> reshape(sol.u[i][1, :], 250, 250), i)
v = map(i -> reshape(sol.u[i][2, :], 250, 250), i)
heatmap!(axu, x, y, u, colorrange = (0.0, 2.5), colormap = :turbo)
heatmap!(axv, x, y, v, colorrange = (0.0, 10.0), colormap = :turbo)
resize_to_layout!(fig)
record(
    fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol);
    framerate = 60) do _i
    i[] = _i
end

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(0, 100, 0, 100, 250, 250, single_boundary = true)
mesh = FVMGeometry(tri)
bc_u = (x, y, t, (u, v), p) -> zero(u)
bc_v = (x, y, t, (u, v), p) -> zero(v)
BCs_u = BoundaryConditions(mesh, bc_u, Neumann)
BCs_v = BoundaryConditions(mesh, bc_v, Neumann)
q_u = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    u = αu * x + βu * y + γu
    ∇u = (αu, βu)
    ∇v = (αv, βv)
    χu = p.c * u / (1 + u^2)
    return χu .* ∇v .- ∇u
end
q_v = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    ∇v = (αv, βv)
    return -p.D .* ∇v
end
S_u = (x, y, t, (u, v), p) -> u * (1 - u)
S_v = (x, y, t, (u, v), p) -> u - p.a * v
q_u_parameters = (c = 4.0,)
q_v_parameters = (D = 1.0,)
S_v_parameters = (a = 0.1,)
u_initial_condition = 0.01rand(DelaunayTriangulation.num_points(tri))
v_initial_condition = zeros(DelaunayTriangulation.num_points(tri))
final_time = 1000.0
u_prob = FVMProblem(
    mesh,
    BCs_u;
    flux_function = q_u,
    flux_parameters = q_u_parameters,
    source_function = S_u,
    initial_condition = u_initial_condition,
    final_time = final_time,
)
v_prob = FVMProblem(
    mesh,
    BCs_v;
    flux_function = q_v,
    flux_parameters = q_v_parameters,
    source_function = S_v,
    source_parameters = S_v_parameters,
    initial_condition = v_initial_condition,
    final_time = final_time,
)
prob = FVMSystem(u_prob, v_prob)

This page was generated using Literate.jl.