CelerisLab/tests/specs/streakline_notes.md

Streakline postprocessing design

这份说明面向一个最小可实现的离线后处理程序。输入是按时间保存的二维速度场序列，输出是接近实验染色图的 streakline 图像。核心路线是连续释粒，因为 streakline 正是“固定位置持续释放染料后，在某个时刻所有已释放粒子的位置集合” [Lan96]。对非定常尾迹，这比单帧 streamline 更贴近实验图像 [Lan96]。

搜索和全文阅读给出的信息已经足够支撑第一版实现。最直接的方案是按保存时刻连续注入粒子，用时空插值后的速度场推进粒子，再将粒子云渲染为图像 [Lan96, Ken96]。如果纯粒子结果过于尖锐，再在推进后加入一个很小的随机扩散项，作为对真实染色液扩散的近似。更完整的被动标量法当然更物理，但实现成本更高，不适合作为第一版 [Kim04]。

Core recommendation

方法	与实验染色图的对应	实现复杂度	适合作为第一版
单帧 streamline	弱	很低	不推荐
连续释粒 streakline	强	低	最推荐
streakline 加少量扩散	很强	低到中	推荐作为第二步
被动标量对流扩散	最强	中到高	暂不作为第一版

最小可实现路线如下：

• 从 CelerisLab 导出一串二维速度场快照 u_x(x,y,t_k), u_y(x,y,t_k)
• 选定一个固定释放点或一小段释放线
• 在每个保存时刻都注入新粒子
• 在两个相邻快照之间，用时间插值和空间插值计算粒子速度
• 用二阶或四阶时间积分推进粒子
• 删去出域粒子与进入固体的粒子
• 将当前时刻存活粒子按位置和年龄渲染成图像

这个程序逻辑直接对应 [Lan96] 的 streakline 算法，插值与推进细节由 [Ken96] 和 [Dar96] 支撑。

Data contract

默认输入是规则 Cartesian 网格上的二维时间序列速度场。对当前项目，最自然的输入约定是：

名称	含义
t_k	第 k 个保存时刻
ux[k, j, i]	时刻 t_k 的 x 方向速度
uy[k, j, i]	时刻 t_k 的 y 方向速度
mask[j, i]	可选，固体与流体标记
x_i, y_j	网格坐标

第一版最好满足下面三条：

• 时间快照间隔固定或至少已知
• 网格坐标固定不动
• 固体几何位置已知，至少能判断粒子是否进入圆柱内部

如果输出来自 CelerisLab，本质上只需要把每个保存时刻的 ux, uy 和对应时间写出来即可。算法本身是通用的，并不依赖 LBM 本身。

Mathematical model for the base streakline

Particle motion

粒子位置满足拉格朗日运动方程 [Lan96]：

\frac{d \boldsymbol{x}}{dt} = \boldsymbol{u}(\boldsymbol{x}, t)

其中

\boldsymbol{x} = (x, y), \qquad \boldsymbol{u} = (u_x, u_y)

积分形式为 [Lan96]：

\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \int_t^{t+\Delta t} \boldsymbol{u}(\boldsymbol{x}(\tau), \tau) \, d\tau

streakline 在观察时刻 t_n 的定义是：从固定释放位置 \boldsymbol{x}_s 在过去各时刻持续注入的粒子，在 t_n 时刻的全部位置集合 [Lan96]。

Time interpolation

速度场只在离散时刻存储，所以必须做时间插值。对位于 t_k \le t \le t_{k+1} 的任意子步，最简单的做法是线性时间插值 [Ken96]：

\delta = \frac{t - t_k}{t_{k+1} - t_k}

\boldsymbol{u}(\boldsymbol{x}, t) = (1 - \delta) \, \boldsymbol{u}_k(\boldsymbol{x}) + \delta \, \boldsymbol{u}_{k+1}(\boldsymbol{x})

这里 \boldsymbol{u}_k(\boldsymbol{x}) 和 \boldsymbol{u}_{k+1}(\boldsymbol{x}) 仍需通过空间插值得到。

Spatial interpolation

在规则网格上，第一版直接用双线性插值即可。虽然 [Ken96] 讨论的是非结构网格上的点定位与线性插值，但其一般原则完全适用于规则网格：

• 先定位粒子所在单元
• 再用单元顶点速度插值得到粒子点速度 [Ken96]

若粒子位于单元局部坐标 ((\xi, \eta) \in [0,1]^2)，四个角点速度为 (\boldsymbol{u}{00}, \boldsymbol{u}{10}, \boldsymbol{u}{01}, \boldsymbol{u}{11})，则

\boldsymbol{u}(\xi, \eta) = (1-\xi)(1-\eta) \boldsymbol{u}_{00}
+ \xi (1-\eta) \boldsymbol{u}_{10}
+ (1-\xi)\eta \boldsymbol{u}_{01}
+ \xi\eta \boldsymbol{u}_{11}

Time integration choice

[Lan96] 和 [Ken96] 都直接使用四阶 Runge Kutta。对你的第一版程序，推荐两个选项：

积分器	优点	缺点	建议
RK2	简洁，容易调试	精度一般	最小原型可用
RK4	文献最一致，精度更稳	每步插值次数更多	默认推荐

RK4 的更新式为 [Lan96, Ken96]：

\boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \frac{1}{6}(\boldsymbol{a} + 2\boldsymbol{b} + 2\boldsymbol{c} + \boldsymbol{d})

其中

\boldsymbol{a} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n, t_n)

\boldsymbol{b} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \tfrac{1}{2}\boldsymbol{a}, t_n + \tfrac{1}{2}\Delta t)

\boldsymbol{c} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \tfrac{1}{2}\boldsymbol{b}, t_n + \tfrac{1}{2}\Delta t)

\boldsymbol{d} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \boldsymbol{c}, t_n + \Delta t)

[Ken96] 特别强调，RK4 的每一个子步都要重新做点定位与时空插值。这在规则网格上不难实现，只是意味着一次完整步进需要四次速度查询 [Ken96]。

Practical timestep guidance

第一版不必上自适应步长，但步长不能随意设大。搜索里最重要的精度提醒来自 [Dar96]：

• 非定常粒子积分的误差常常由时间离散控制，而不是概念本身
• 合理步长必须不大于流动主要非定常时间尺度的同量级 [Dar96]
• 若时间步过大，粒子轨迹会在拓扑上都出错，而不仅是位置略偏 [Dar96]

对当前项目，最简单的工程准则是：

\Delta t_{trace} \le \min \bigl(\alpha_t \, \Delta t_{save},\; \alpha_x \, \frac{\Delta x}{\max |u|} \bigr)

其中 \Delta t_{save} 是两个保存快照之间的时间间隔，\Delta x 是网格尺度。第一版可取：

• \alpha_t = 0.1 \sim 0.25
• \alpha_x = 0.25 \sim 0.5

这不是文献中的严格上界，而是结合 [Dar96] 的结论给出的实现准则：粒子推进子步既要显著小于保存间隔，也不要一子步跨过太多网格。

[Ken96] 还给了一个很实用的自适应思路：根据相邻速度方向夹角调节步长。如果速度方向变化太快就减半，变化很小就加倍 [Ken96]。这很适合后续增强版，但不属于第一版必需项。

Base streakline algorithm

Release strategy

最符合染色实验的是连续释粒 [Lan96]。第一版可用两种释放方式：

方式	适用场景	建议
单点释放	针头式染料入口	最简单
短线段释放	更像细缝或细带注入	更稳健

如果实验是在圆柱上游某一点持续释放染色液，第一版就直接使用单点释放。

Particle state

每个粒子只需保存：

字段	含义
x, y	当前坐标
t_birth	释放时刻
age	当前年龄
alive	是否仍在域内

第一版不需要保存整条历史轨迹。因为你要的是当前时刻的 streakline 图，而不是每个粒子的 pathline 曲线。

常见误区（务必区分）

概念	画什么	是否等于水洞染色迹线
Streakline（迹线）	固定源点持续释粒，在观察时刻 t_n 画出所有仍存活粒子的当前位置（可带年龄衰减权重）	是
Pathline（迹线/轨道）	单个粒子从释放到当前的整条历史轨迹	否
错误实现	每隔一段时间放一个粒子，再把所有粒子走过的路径段全部叠加到图上	否（这是 pathline 叠加，不是 streakline）

水洞实验：针头在固定点连续注 dye → 某一时刻拍照 → 看到的是“此刻染料粒子在流场里的分布”，一条下游色带由不同释放时刻、当前仍在本流场中的粒子共同构成，而不是把每个粒子从出生到现在的轨迹都画出来。

Main loop pseudocode

given velocity snapshots U[k] = {ux[k], uy[k]} at times t[k]
given seeding point or seeding segment S
initialize empty particle list P

for k = 0 to N-2:

    inject new particles at source S at time t[k]

    set t_local = t[k]
    while t_local < t[k+1]:
        dt = min(dt_trace, t[k+1] - t_local)

        for each particle p in P with p.alive:
            v = interpolate_velocity(p.position, t_local)
            advance p by one integration step using RK2 or RK4
            update p.age
            if p leaves domain:
                p.alive = false
            if p enters solid body:
                p.alive = false

        t_local = t_local + dt

    optionally remove very old particles
    optionally render current particle cloud

这个结构与 [Lan96] 的离散 streakline 算法一致，只是把“从 t_k 到 t_{k+1} 的一次推进”细化成多个更小的粒子子步，以满足 [Dar96] 对精度的要求。

Velocity query pseudocode

function interpolate_velocity(position x, time t):
    find k such that t[k] <= t <= t[k+1]
    compute delta = (t - t[k]) / (t[k+1] - t[k])

    u_k   = bilinear_interpolation(U[k],   x)
    u_k1  = bilinear_interpolation(U[k+1], x)

    return (1 - delta) * u_k + delta * u_k1

这是第一版最核心的数值部件。只要这部分实现正确，streakline 程序主体就很直接。

Small diffusion particle model

当纯 streakline 太细、太锐利、不像染色液图像时，可以给每个粒子加一个很小的随机扩散项。这不是严格的被动标量求解，但能以很低代价增加染色带宽度。

Guiding idea

[Kim04] 表明，染料图像更接近一个被动标量浓度场，其基本控制方程是对流扩散方程：

\frac{\partial c}{\partial t} + \boldsymbol{u} \cdot \nabla c = D \nabla^2 c

其中 c 是染料浓度，D 是扩散系数。对第一版粒子法，一个常见近似是将每个粒子的位置更新写成“对流加随机扩散”：

\boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \Delta \boldsymbol{x}_{adv} + \Delta \boldsymbol{x}_{diff}

其中 \Delta \boldsymbol{x}_{adv} 由 RK2 或 RK4 给出，扩散项取二维各向同性随机增量：

\Delta \boldsymbol{x}_{diff} = \sqrt{2 D \Delta t}
\begin{bmatrix}
\eta_x \\
\eta_y
\end{bmatrix}

这里 (\eta_x, \eta_y \sim \mathcal{N}(0,1))。

这个形式本身是对扩散过程的标准随机游走近似。搜索结果提醒，随机游走模型如果处理不当会产生假漂移与错误浓度偏置 [Mac92]。因此第一版的使用原则应当很克制：

• 只用于加入少量模糊和厚度
• 不把粒子密度当成严格浓度
• 不在强近壁统计上过度解读结果

Diffusive particle update pseudocode

for each particle p in P with p.alive:
    x_adv = RK4_step(p.position, t_local, dt)

    sigma = sqrt(2 * D * dt)
    dx_rand = sigma * normal(0, 1)
    dy_rand = sigma * normal(0, 1)

    x_new = x_adv + [dx_rand, dy_rand]

    if x_new leaves domain:
        p.alive = false
    else if x_new enters solid:
        p.alive = false
    else:
        p.position = x_new
        p.age += dt

Choosing the diffusion level

第一版不必试图从真实染料物性严格标定 (D)。更实用的做法是把 D 当作视觉匹配参数，并保持它足够小，使图像结构仍主要由对流控制。

[Kim04] 用 Schmidt 数控制扩散强弱。若已有参考速度 U 和长度尺度 (L)，则

Sc = \frac{\nu}{D}

也可写成

D = \frac{\nu}{Sc}

对第一版，可以用下面的思路选扩散强度：

目标效果	建议
只想让线条略微变厚	取较大 Sc，即很小的 D
想模拟明显洗开和模糊	取较小 Sc，即较大的 D
不确定	先从几乎看不出的弱扩散开始

因为当前目标是“少量扩散”，所以推荐先把扩散当成弱修饰，而不是主导机制。

Rendering logic

最终图像不需要把每个粒子的完整轨迹都画出来。更像实验染色图的做法是把当前存活粒子投影到图像网格上，生成粒子密度图或带年龄权重的强度图。

Minimal rendering choices

方法	图像风格	实现难度
直接散点	最简陋	很低
网格计数直方图	像浓度图	低
高斯核累积	更平滑	低到中

第一版推荐：

• 把每个粒子投到像素网格
• 对像素做计数或加权累积
• 最后做一次轻微 Gaussian blur

这与先做完整被动标量相比便宜很多，但视觉上已经会很接近实验染色图。

Optional particle weighting

可以给粒子一个简单权重：

I = \sum_p w_p K(\boldsymbol{x} - \boldsymbol{x}_p)

其中 K 是像素核或 Gaussian 核。第一版里，权重 w_p 可直接取 1，也可对年龄做衰减，例如：

w_p = \exp\left(- \frac{\mathrm{age}_p}{\tau_f} \right)

这样旧粒子会逐渐淡出，图像不会无限堆积。

Boundary handling

这部分先按最小原则处理即可。

情况	第一版处理
粒子出计算域	直接删除
粒子进入圆柱内部	直接删除
粒子贴近边界滑动	暂不专门处理
扩散后跨入固体	直接删除

这样做的优点是简单稳妥。后续若发现近壁 streakline 误差明显，再考虑反射、投影或更物理的壁面处理。

Minimal implementation plan

Version 1

目标是尽快得到可信的 streakline 图：

• 读取 ux, uy, t
• 单点连续释粒
• 双线性空间插值
• 线性时间插值
• RK4 粒子推进
• 出域删除与入固体删除
• 粒子计数成图

Version 2

在不改变主体架构的前提下增强视觉效果：

• 增加短线段释放
• 增加年龄衰减
• 增加弱随机扩散
• 用 Gaussian 核代替简单计数

Version 3

若后续发现实验图像明显受扩散与混合主导，再转向更重的模型：

• 被动标量对流扩散
• 用 Sc 或 D 做参数标定 [Kim04]

Design decisions that are already justified by the literature

下面这些设计已经有足够文献支撑，可以直接采用：

• 用 streakline 而不是 streamline 对应持续染色实验 [Lan96]
• 用连续释放粒子离线重构 streakline [Lan96]
• 用时空插值从离散速度快照查询粒子速度 [Ken96]
• 用 RK4 作为默认推进器 [Lan96, Ken96]
• 把步长取得明显小于保存时间间隔，并避免一次跨过太多网格 [Dar96]
• 若需要更像染料图，可以加入弱扩散，或后续上升到被动标量 [Kim04]

Recommended default choices

项目	默认选择
释放方式	单点连续释粒
空间插值	双线性
时间插值	线性
积分器	RK4
粒子子步	保存步长的 0.1 到 0.25
固体处理	入固体即删除
输出图像	粒子密度图加轻微模糊
扩散	默认关闭，作为第二步

Bottom line

对当前项目，最简洁且足够准确的程序逻辑不是去解新的染料场，而是先做一个离线连续释粒 streakline 后处理器。它只依赖时间序列速度场，程序结构清楚，数值风险也集中在可控的几个环节：时空插值、积分步长和边界删除 [Lan96, Dar96, Ken96]。在此基础上，再增加一个弱随机扩散项，就能以很小代价把图像从“几何上正确的粒子线”推进到“更像实验染色照片”的粒子云 [Kim04]。

CelerisLab implementation mapping

Current implementation lives in:

• library: src/CelerisLab/common/streakline.py
• Kan99b demo CLI: tests/run_kan99b_streakline.py
• experiment notebook: tests/experiment.ipynb (CLEAN 5/6)

Available runtime modes:

• online: sample Simulation.get_macroscopic() directly in memory (no velocity snapshot files).
• offline: replay from an existing snapshot directory and run the same streakline integrator.

Dense release support:

• base upstream points are expanded by release_mode (point|line|strip).
• strip mode densifies both cross-stream (line_count, line_span) and downstream (downstream_count, downstream_spacing) directions.

Performance notes (why it can feel slow)

Streakline post-processing is CPU-side in the current stack:

1. GPU→host velocity copies dominate for large grids (nx=6000, ny=1200). Each get_macroscopic() pulls full ux/uy arrays.
2. Particle RK4 + bilinear interpolation runs in NumPy on CPU (vectorized over particles, not GPU).
3. Rendering accumulates polylines into an image (optionally multi-threaded via ThreadPoolExecutor).

Practical tuning:

knob	effect
sample_every	larger → fewer host copies, faster, coarser streaks
max_particle_age	None → no hard cutoff; trails survive until boundary/solid
num_threads	0 = auto; sets OpenBLAS/OMP threads for blur/host ops
blur_sigma	0 disables Gaussian blur pass

For experiment.ipynb triangle case, use render_streakline_density(..., minimal_axes=True) on the final alive particle cloud (positions, ages). Do not accumulate and draw full path history (render_snapshot_trails is pathline-style debug only).

Clean render style (experiment)

render_snapshot_trails defaults:

• white background
• red streaks, brighter toward the downstream end (fade_along_trail=True)
• black filled cylinders
• no axes, labels, colorbar, or release-point markers

Example commands:

# Online in-memory streakline (no snapshot write)
conda run -n pycuda_3_10 python tests/run_kan99b_streakline.py \
  --mode online --domain M --re 100 --alpha 1.0 \
  --sample-every 300 --n-snapshots 20 \
  --release-mode strip --line-count 7 --downstream-count 6

# Offline replay from existing velocity snapshots
conda run -n pycuda_3_10 python tests/run_kan99b_streakline.py \
  --mode offline \
  --snapshot-dir tests/output/final_validation_round/streakline_kan99b_k2/velocity_snapshots \
  --release-mode strip --line-count 7 --downstream-count 6