求解对流扩散后比较线性与自然邻域插值

作者： FiniteVolumeMethod.jl 文档维护者（原页面未署个人名）
发表日期： 原页面未注明
原文： Piecewise Linear and Natural Neighbour Interpolation for an Advection-Diffusion Equation
来源： SciML / FiniteVolumeMethod.jl 官方文档

目标与方程

本教程依次演示三件事：求解一个二维对流扩散方程、在每个时间点对有限体积解做分片线性插值，以及使用 NaturalNeighbours.jl 比较其他自然邻域插值方式。

考虑方程：

[ = D + D – , ]

初始条件为 (u(,0)=())，并采用齐次 Dirichlet 边界。理论区域是 (^2)，示例把它截为 (=[-L,L]^2)，取 (L=30)。

建立在原点附近加密的三角网格

示例先建立覆盖方形区域的初始三角网格，再通过自定义面积约束细化网格，让原点附近得到更多三角形。约束根据每个三角形质心到原点的距离计算允许面积，再调用 refine!：

using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie

L = 30
tri = triangulate_rectangle(-L, L, -L, L, 2, 2, single_boundary=true)
tot_area = get_area(tri)

max_area_function = (A, r) -> 1e-6tot_area * r^2 / A

area_constraint = (_tri, T) -> begin
    u, v, w = triangle_vertices(T)
    p, q, r = get_point(_tri, u, v, w)
    c = (p .+ q .+ r) ./ 3
    dist_to_origin = norm(c)
    A = DelaunayTriangulation.triangle_area(p, q, r)
    flag = A ≥ max_area_function(A, dist_to_origin)
    return flag
end

refine!(tri; min_angle=33.0, custom_constraint=area_constraint)
triplot(tri)

mesh = FVMGeometry(tri)

原文的示例输出将该几何对象描述为包含 8,541 个控制体、16,947 个三角形和 25,487 条边。这是文档中的展示输出，并非本稿运行得到的结果。

把扩散与平流写成通量

边界条件在整条边界上取零：

BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> zero(u), Dirichlet)

有限体积问题需要把方程改写为守恒形式 (_t u+=0)。将原方程中的平流项与扩散项合并，可得

[ =u-Du, =(,0)^. ]

对有限体积方法中的局部线性近似 (u(x,y)=x+y+)，梯度为 ((_xu,_yu)=(,))。相应的通量函数和参数如下：

ε = 1 / 10
f = (x, y) -> 1 / (ε^2 * π) * exp(-(x^2 + y^2) / ε^2)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]

flux_function = (x, y, t, α, β, γ, p) -> begin
    ∂x = α
    ∂y = β
    u = α * x + β * y + γ
    qx = p.ν * u - p.D * ∂x
    qy = -p.D * ∂y
    return (qx, qy)
end

flux_parameters = (D=0.02, ν=0.05)
final_time = 250.0

prob = FVMProblem(mesh, BCs;
    initial_condition,
    flux_function,
    flux_parameters,
    final_time)

点源 (()) 用高斯函数近似，()。本例选择扩散系数 (D=0.02)、平流参数 ()，并积分到 (t=250)。

求解并查看时间演化

示例使用 OrdinaryDiffEq 的 TRBDF2，并让 KLUFactorization 处理线性系统；保存时刻为 0、10、25、50、100、200 和 250：

using OrdinaryDiffEq, LinearSolve

times = [0, 10, 25, 50, 100, 200, 250]
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=times)

原文用 CairoMakie 的 tricontourf! 在三角网格上绘制多个保存时刻的解。可视化展示了随后插值的对象：非均匀三角网格上的有限体积解。

将分片线性解映射到规则网格

有限体积方法假定每个三角单元内的 (u) 是分片线性的，因此可以先找出规则网格每个点所在的三角形，再用 pl_interpolate 求值。示例把目标网格设为两个方向各 250 个点：

x = LinRange(-L, L, 250)
y = LinRange(-L, L, 250)

triangles = Matrix{NTuple{3,Int}}(undef, length(x), length(y))
for j in eachindex(y)
    for i in eachindex(x)
        triangles[i, j] = jump_and_march(tri, (x[i], y[j]))
    end
end

interpolated_vals = zeros(length(x), length(y), length(sol))
for k in eachindex(sol)
    for j in eachindex(y)
        for i in eachindex(x)
            interpolated_vals[i, j, k] =
                pl_interpolate(prob, triangles[i, j], sol.u[k], x[i], y[j])
        end
    end
end

这里把每个目标点所属的三角形先存起来，再对每个保存时刻进行插值。如果每个目标点只需评估一个时刻，也可以先计算对应时刻的插值，而不保存整张三角形索引矩阵。原文还说明，先批量定位所需三角形再插值，通常比对每个点重复查找更有效。

用 NaturalNeighbours.jl 比较不同插值器

在三角剖分上，除分片线性方法外，自然邻域插值也是一种自然的选择。先为 (t=50) 的解构造插值器；因为后面会用到高阶方法，需要启用导数：

using NaturalNeighbours

itp = interpolate(tri, sol.u[4], derivatives=true) # sol.t[4] == 50

若只用普通一阶插值，derivatives=true 并非必需；如果之后想用 differentiate 对插值结果求导，也需要导数信息。

把规则网格坐标收集成两个向量，再一次传入插值器。原文指出，批量传点比逐点广播更有效，因为该调用形式可以使用多线程：

_x = [x for x in x, _ in y] |> vec
_y = [y for _ in x, y in y] |> vec

教程将 NaturalNeighbours.jl 提供的方法与 PDE 解作图比较，方法包括：

sibson_vals = itp(_x, _y; method=Sibson())
triangle_vals = itp(_x, _y; method=Triangle()) # 与 pl_interpolate 相同
laplace_vals = itp(_x, _y; method=Laplace())
sibson_1_vals = itp(_x, _y; method=Sibson(1))
nearest_vals = itp(_x, _y; method=Nearest())
farin_vals = itp(_x, _y; method=Farin())
hiyoshi_vals = itp(_x, _y; method=Hiyoshi(2))
pde_vals = sol.u[4]

比较时包括 Sibson、Triangle、Laplace、Sibson(1)、Nearest、Farin 与 Hiyoshi(2)，并把 (t=50) 的 PDE 解作为参照。原文示例通过二维等高线和三维表面展示各结果；这里保留方法与比较流程，不附加图片，也不宣称已执行这些计算或获得新的数值结论。

文档在 interpolate 调用后的展示块仍出现 ODE 解对象的 retcode、t 和 u 输出，而不是插值器的单独结果。因此，不应把该展示块误读为自然邻域插值计算的验证输出。

受约束网格上的边界限制

自然邻域插值在受约束三角剖分上并非严格总是有定义。文档指出，本例使用时可以工作；但如果区域带孔或边界非凸，可能出现问题。对于这类情形，原文建议通常尝试在调用插值器时使用 project=false，并指出 identify_exterior_points 也可能有帮助。不要据此假设对任意带孔或非凸区域都能直接得到良好定义的插值。

环境与复现边界

Notion 核验记录标注所读的是 stable 文档：页面设置区记录 Documenter.jl 1.8.0 于 2025-01-01 生成，并注明使用 Julia 1.11.2。构建日期不是文章发布日期；该教程没有锁定 FiniteVolumeMethod 或其他依赖的具体版本，因此不能由此保证示例与当前最新版兼容。本稿没有运行数值代码，也没有下载或附加原文图片。

许可说明： 原页面未列出明确的文章转载许可；Notion 核验记录也说明没有检查仓库许可证。