模拟 Gray–Scott 反应扩散并生成图灵斑图动画

原文标题: Gray-Scott Model: Turing Patterns from a Coupled Reaction-Diffusion System
作者: FiniteVolumeMethod.jl 文档维护者(页面未署个人名)
来源: FiniteVolumeMethod.jl 官方教程
原文日期: 页面未注明。页面页脚记载其由 Documenter.jl 1.8.0 于 2025 年 1 月 1 日生成,使用 Julia 1.11.2;这是文档构建信息,不是文章首发日期。

Gray–Scott 模型描述两种化学物质在空间中的反应与扩散。FiniteVolumeMethod.jl 的教程把两个浓度场写成耦合偏微分方程,在二维正方形区域上用有限体积网格离散,设置零通量边界条件并积分到较长时间。随后,教程取第二个浓度场的数值结果,重塑为二维数组,用 CairoMakie 绘制热图并录制动画。

模型和初始条件

令 u、v 表示两种化学物质的浓度。Gray–Scott 系统为:

∂u/∂t = ε₁ ∇²u + b(1 − u) − uv²
∂v/∂t = ε₂ ∇²v − dv + uv²

其中,ε₁ 和 ε₂ 分别控制两种物质的扩散,b 和 d 是反应项参数。教程使用以原点为中心的高斯初始分布:

u(x, y, 0) = 1 − exp(−80(x² + y²))
v(x, y, 0) =     exp(−80(x² + y²))

计算区域为 [-1, 1]²,两个场的边界均采用零通量条件。

创建网格与设置边界条件

教程使用 DelaunayTriangulation 生成正方形网格,并把它包装为 FiniteVolumeMethod 的 FVMGeometry。网格参数是每个方向 200 个点,并启用单一边界:

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

文档显示,该几何对象包含 40,000 个控制体积、79,202 个三角形和 119,201 条边。两个变量使用同一个零通量函数,并分别创建 Neumann 边界条件:

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

定义通量、反应项和初始场

扩散通量使用各变量的梯度分量乘以相应的扩散系数;源项则实现 Gray–Scott 方程中的反应项。原文代码如下:

ε₁ = 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、v 建立 FVMProblem。两者共用同一几何网格与边界类型,最终时间为 6000;接着把两个问题组合为一个耦合系统:

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)

教程页面显示生成了包含两个方程、时间范围为 0.0 到 6000.0 的 FVMSystem。

积分并保存解

教程选择 OrdinaryDiffEq 的 TRBDF2 时间积分方法,并通过 LinearSolve 的 KLUFactorization 求解线性系统。每隔 10 个时间单位保存一次结果,关闭并行执行:

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

页面中嵌入的示例输出显示求解状态为 Success,保存时间点共有 601 个,从 0 到 6000,间隔为 10。教程还展示了各保存时刻的二维解矩阵。这里引用的是原文页面中的示例输出,并非本文重新运行所得。

把第二个浓度场录制成动画

动画只显示 v 变量。解向量的第二行对应耦合系统的第二个场;原文将这一行按 200×200 重塑为二维网格数据,通过 Observable 更新时间索引,再交给 CairoMakie 的热图和 record 函数生成 MP4:

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

原文把输出写到当前脚本目录上一级的 figures 子目录,文件名为 gray_scott_patterns.mp4,帧率为每秒 60 帧。该目录需要可写;若把示例移出文档构建环境,应根据本机项目结构调整输出路径。教程使用了 40,000 个控制体积,并保存 601 个时刻的解,因此执行时间和内存开销都应纳入本地资源规划。

复现时的版本与许可信息

该页面属于 FiniteVolumeMethod.jl 的 stable 文档,页脚明确记录了文档生成器和 Julia 构建版本,但没有为这段教程锁定 FiniteVolumeMethod.jl、OrdinaryDiffEq、LinearSolve、CairoMakie 等依赖包的具体版本。因此,这些代码是否与当前安装环境完全兼容,需要由实际项目环境决定。

页面没有署个人作者,也未注明文章首发日期;页脚中的生成日期不应当作发布日期。该文档页面未列出单篇文章的转载许可说明。本文保留项目文档身份、构建信息和原始教程链接;没有运行代码或下载动画。

来源说明: 本文依据 SciML / FiniteVolumeMethod.jl 指定教程页面整理与中文化。原文页面由 Literate.jl 生成.

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容