📖 第18课：扩散过程仿真

扩散 Fick定律热传导 LBM

🌡️ 扩散——自然界最基本的过程

扩散是无处不在的物理过程：墨水在水中扩散、热量在金属中传导、信息在社交网络中传播。它们都遵循相同的数学框架——扩散方程。

扩散方程（Fick第二定律）：

∂c/∂t = D∇²c

其中 c = 浓度/温度/信息密度，D = 扩散系数

离散化：c_i(t+1) = c_i(t) + D·Δt/Δx² × Σ_neighbors(c_j - c_i)

⚡ 扩散CA的Verilog实现

定点数扩散CA引擎

// ============================================================================
// diffusion_ca.v - 扩散过程CA引擎
// 使用定点数表示浓度，Moore邻域Laplacian
// 支持多种边界条件和源项
// ============================================================================
module diffusion_ca #(
    parameter WIDTH  = 64,
    parameter HEIGHT = 64,
    parameter FRAC_W = 8,           // 小数部分位宽
    parameter [15:0] DIFF_COEFF = 16'h0100  // 扩散系数 D (8.8定点)
)(
    input  wire                     clk,
    input  wire                     rst_n,
    input  wire                     enable,
    input  wire                     init,
    // 浓度输入（8位整数+8位小数 = 16位/元胞）
    output wire [WIDTH*HEIGHT*16-1:0] concentration,
    output wire [31:0]              step_count,
    output wire [15:0]              total_mass    // 总质量守恒验证
);

    // 每格点16位浓度（8.8定点）
    reg [15:0] conc [0:WIDTH*HEIGHT-1];
    reg [15:0] conc_nxt [0:WIDTH*HEIGHT-1];
    reg [31:0] steps;

    integer idx, x, y;
    reg signed [19:0] laplacian;  // 20位有符号，足够容纳累加

    always @(*) begin
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            x = idx % WIDTH;
            y = idx / WIDTH;

            // 计算Laplacian: Σ(c_neighbor - c_center)
            laplacian = 0;
            // 8个Moore邻居
            if (y > 0) begin
                laplacian = laplacian + $signed({1'b0, conc[(y-1)*WIDTH+x]}) - $signed({1'b0, conc[idx]});
                if (x > 0)  laplacian = laplacian + $signed({1'b0, conc[(y-1)*WIDTH+x-1]}) - $signed({1'b0, conc[idx]});
                if (x < WIDTH-1) laplacian = laplacian + $signed({1'b0, conc[(y-1)*WIDTH+x+1]}) - $signed({1'b0, conc[idx]});
            end
            if (y < HEIGHT-1) begin
                laplacian = laplacian + $signed({1'b0, conc[(y+1)*WIDTH+x]}) - $signed({1'b0, conc[idx]});
                if (x > 0)  laplacian = laplacian + $signed({1'b0, conc[(y+1)*WIDTH+x-1]}) - $signed({1'b0, conc[idx]});
                if (x < WIDTH-1) laplacian = laplacian + $signed({1'b0, conc[(y+1)*WIDTH+x+1]}) - $signed({1'b0, conc[idx]});
            end
            if (x > 0)  laplacian = laplacian + $signed({1'b0, conc[y*WIDTH+x-1]}) - $signed({1'b0, conc[idx]});
            if (x < WIDTH-1) laplacian = laplacian + $signed({1'b0, conc[y*WIDTH+x+1]}) - $signed({1'b0, conc[idx]});

            // 扩散更新: c_new = c + D * laplacian
            // D是8.8定点，laplacian是~16位
            conc_nxt[idx] = conc[idx] + (laplacian >>> FRAC_W);
        end
    end

    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                conc[idx] <= 16'd0;
        end else if (init) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                conc[idx] <= 16'd0;
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                conc[idx] <= conc_nxt[idx];
            steps <= steps + 32'd1;
        end
    end

    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1)
            assign concentration[gi*16 +: 16] = conc[gi];
    endgenerate

    reg [31:0] mass;
    always @(*) begin
        mass = 0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
            mass = mass + conc[idx];
    end
    assign step_count = steps;
    assign total_mass = mass[15:0];

endmodule

热传导仿真

扩散方程与热传导方程在数学上等价——只需将"浓度"替换为"温度"。

热传导方程：∂T/∂t = α∇²T

其中α是热扩散率

经典问题：正方形区域的稳态温度分布，四边固定不同温度

解析解可用Fourier级数展开

💡 数值稳定性：扩散CA的稳定性条件是 D·Δt/Δx² < 1/4（2D Moore邻域）。如果扩散系数太大，数值会振荡甚至发散！在硬件中，可以使用饱和截断来防止溢出。

🏋️ 练习

练习18.1：模拟一维扩散——初始时左半区域浓度=1.0，右半=0.0。运行1000步，画出浓度分布曲线。与解析解对比。

练习18.2：模拟稳态热传导——正方形区域四边温度分别固定为100°、0°、50°、25°。运行至收敛，画出等温线。

练习18.3：添加点热源——在网格中心持续注入热量。观察热量的扩散过程。稳态分布是什么？

练习18.4：验证质量守恒——每步后计算总质量，确认其不随时间变化。如果出现漂移，分析原因。

练习18.5（挑战）：实现各向异性扩散——x方向和y方向使用不同的扩散系数。这在图像处理（边缘保护平滑）中有应用。

🏆 成就解锁

🏅 扩散模拟师

✅ 实现定点数扩散CA引擎
✅ 理解Laplacian的硬件计算
✅ 掌握数值稳定性条件
✅ 了解扩散与热传导的等价性

📐 深入分析：扩散方程的数值方法对比

数值方法对比：

1. 前向Euler（显式）：c^{n+1} = c^n + D·Δt·∇²c^n

稳定性：D·Δt/Δx² ≤ 0.25 (2D)

优点：简单，每步只需1次评估

缺点：时间步受限制

2. 后向Euler（隐式）：c^{n+1} - D·Δt·∇²c^{n+1} = c^n

稳定性：无条件稳定

优点：可以大时间步

缺点：需要解线性方程组（硬件不友好）

3. Crank-Nicolson：c^{n+1} - D·Δt/2·∇²c^{n+1} = c^n + D·Δt/2·∇²c^n

稳定性：无条件稳定

精度：O(Δt²)，比Euler的O(Δt)高

硬件实现推荐：前向Euler，因为只需要1次Laplacian计算

补充实现

// 隐式扩散求解器的简化版 - 用迭代法近似
module diffusion_iterative #(
    parameter WIDTH = 64,
    parameter HEIGHT = 64,
    parameter ITERATIONS = 4
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire              enable,
    input  wire [15:0]       conc [0:WIDTH*HEIGHT-1],
    output wire [15:0]       conc_out [0:WIDTH*HEIGHT-1]
);
    // Jacobi迭代：c_new[i] = (c[left]+c[right]+c[up]+c[down]) / 4
    reg [15:0] buf_a [0:WIDTH*HEIGHT-1];
    reg [15:0] buf_b [0:WIDTH*HEIGHT-1];
    integer iter, idx, x, y;
    
    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            // 初始化
        end else if (enable) begin
            // 复制输入到buf_a
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                buf_a[idx] <= conc[idx];
            // 迭代
            for (iter = 0; iter < ITERATIONS; iter = iter + 1) begin
                for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                    x = idx % WIDTH; y = idx / WIDTH;
                    buf_b[idx] = (buf_a[(y>0?(y-1)*WIDTH+x:(HEIGHT-1)*WIDTH+x)] +
                                  buf_a[(y0?x-1:WIDTH-1)] +
                                  buf_a[y*WIDTH+(x> 2;
                end
                for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                    buf_a[idx] <= buf_b[idx];
            end
        end
    end
endmodule

性能与优化分析

扩散的应用：

1. 图像处理：高斯模糊=各向同性扩散，边缘保护平滑=各向异性扩散

2. 热管理：芯片热仿真的核心——每个计算单元产生热量，通过扩散散开

3. 金融数学：Black-Scholes方程本质上是扩散方程

4. 神经科学：神经递质在突触间隙的扩散

📖 扩展阅读

Wolfram, S. A New Kind of Science, Wolfram Media, 2002.
Toffoli, T. & Margolus, N. Cellular Automata Machines, MIT Press, 1987.
Ilachinski, A. Cellular Automata: A Discrete Universe, World Scientific, 2001.

🔬 补充专题：实验方法论

在进行CA实验时，科学的方法论至关重要。以下是一些通用的实验指导原则：

实验设计三要素：

控制变量：每次只改变一个参数，保持其他不变
可重复性：记录所有参数，确保实验可以精确复现
统计显著性：多次运行取平均，避免偶然结果误导

参数扫描方法

系统性地扫描参数空间是理解CA行为的关键技术：

参数	范围	步长	测量指标
规则号	0-255	1	种群密度/周期
初始密度	0.1-0.9	0.1	收敛时间
网格大小	16-256	×2	有限尺寸效应
边界条件	环形/固定/镜像	离散	边界效应

💡 实验记录模板：每次实验记录以下信息——日期、规则号、初始条件、网格大小、边界条件、运行步数、关键观察、数据文件路径。这是科学研究的良好习惯，也使得结果可以追溯和验证。

数据可视化技巧

CA实验产生的数据通常是高维的（空间+时间+状态）。有效的可视化对于理解至关重要：

时空图：横轴=空间位置，纵轴=时间，颜色=状态——最常用的CA可视化方式
相图：x轴=参数A，y轴=参数B，颜色=测量指标——用于参数空间扫描
频谱图：对种群时间序列做FFT，识别周期性
直方图：统计状态/密度的分布

自相关分析：

时间自相关：R(τ) = ⟨n(t)·n(t+τ)⟩ / ⟨n²⟩

空间自相关：C(r) = ⟨n(x)·n(x+r)⟩ / ⟨n²⟩

如果R(τ)以周期T振荡 → 系统有周期T的行为

如果C(r)幂律衰减 → 系统处于临界状态

🔧 工程实践指南

将CA从理论变为可运行的硬件系统需要解决许多工程细节。以下是基于实践经验的详细指南：

时钟域设计

CA系统通常需要多个时钟域：

clk_ca：CA计算时钟，100-200MHz
clk_vga：VGA像素时钟，25-165MHz
clk_uart：UART波特率时钟，115200Hz
clk_ddr：DDR内存接口时钟，200-400MHz

跨时钟域同步使用双触发器或异步FIFO

资源优化清单

优化项	方法	节省量
邻居计数器	增量更新替代全加法	~40% LUT
状态存储	BRAM替代分布式RAM	~60% FF
规则查找	硬编码XOR替代MUX	~75% LUT（XOR规则）
显示输出	行缓冲替代全帧缓冲	~50% BRAM
边界处理	环形替代固定（零开销）	0

调试技巧

CA系统的调试有其特殊性：

守恒律验证：粒子数守恒、动量守恒、能量守恒——每步检查
已知模式测试：用已知的经典图案（滑翔机、Sierpinski等）验证规则正确性
对偶验证：软件模拟和硬件实现对比，用相同输入产生相同输出
边界测试：特别注意边界条件——环形边界最容易出错
初始条件敏感：不同的初始种子可能暴露不同的bug

⚠️ 常见陷阱：

忘记双缓冲——读写同一状态导致部分元胞使用t+1值
环形边界索引错误——首尾元胞的邻居指向错误位置
定点数溢出——大网格的种群计数器可能溢出
跨时钟域未同步——导致亚稳态和随机错误
VGA时序违反——消隐期内输出数据导致显示异常

性能基准测试

为CA引擎建立性能基准：

关键性能指标：

SPS（Steps Per Second）：每秒完成的CA步数
PPS（Pixel Updates Per Second）：每秒更新的元胞数 = SPS × W × H
延迟：从输入变化到输出响应的时间
功耗：每百万元胞更新所需的能量

理论峰值PPS = f_clk × W × H（全并行）

实际PPS取决于架构——行缓冲约为理论值的1/H

◀ 第17课 📜 目录第19课 ▶

元胞自动机课程 · 从Conway到Langton到Lattice Gas