简介SPH流体仿真项目资源包面向流体力学模拟、计算机图形学与物理仿真方向的开发者。项目基于光滑粒子动力学SPH方法用C与OpenGL实现流体运动过程通过大量粒子求解动力学方程并实时渲染可帮助学习者理解粒子法流体模拟的完整流程。压缩包共75个文件包含25个头文件、13个C源文件另有若干CUDA源文件与头文件用于GPU加速计算以及项目工程文件、静态库和动态链接库能直接支撑编译与调试整体大小仅5.64MB。目前已有1057人学习下载适合具备一定C和图形学基础的读者进行源码研读与实践。资源中提供了完整的工程配置与可执行程序便于快速运行观察效果也可在此基础上修改粒子数量、边界条件或外力场进一步掌握SPH核函数估算、邻居搜索、压强求解等关键环节。1. 粒子代替网格的流体模拟SPH 能算什么、不能算什么第一次接触 sph流体仿真大多数人都是冲着“水花四溅、飞沫乱飞”来的真正跑起来之后才会意识到这个算法的价值不在画面而在它把流体方程换成了一种完全不同的离散方式没有网格、没有拓扑重建每个采样点自带质量和速度流体的大变形、飞溅和破碎都变成粒子位置的重新排列。SPH 的基本思路是用一组“光滑核”把邻居粒子的物理量加权平均把连续场方程变成一组可迭代的求和。它特别适合自由表面剧烈变化的问题比如溃坝、泼溅、液滴撞击液面缺点也写在明面上——粒子数上去后算力陡增参数没调好会出现粒子炸飞或扎堆。这套方法适合三类人需要自己掌握流体内核的仿真工程师、做物理动画的技术美术、以及研究粒子方法的算法人员。如果你只是想要一个“一键出效果”的流体工具SPH 不是首选如果你要的是把物理过程完整握在手里它值得投入。下文按我自己的落地路径讲先立住物理再给最小可运行实现然后谈参数和排坑最后说怎么验证结果。2. 从 Navier-Stokes 到粒子求和SPH 的物理内核与核函数2.1 拉格朗日视角粒子真正承载的是速度场而不是“水分子”先弄清一个常被误解的点SPH 粒子不是水分子它是携带物质属性的体积单元。每个粒子在初始化时继承了质量、位置、速度、密度后续所有物理量都跟着粒子走。这和网格法的欧拉视角正好相反——欧拉法盯着固定的空间点流体流过时该点的速度、压力随时间变化SPH 盯着采样点本身点移动到哪里属性就带到哪里。这个视角带来的第一个好处是质量守恒变得极其自然。每个粒子质量在初始化后保持不变总质量就是所有粒子质量之和不存在网格法中“单元间通量计算误差”导致的质量流失。第二个好处是自由表面不需要额外追踪自由表面就是“有粒子和没粒子”的边界溃坝里的波浪翻卷、液滴脱离水面本质上都是粒子集合的重新分布没有网格畸变问题。代价在于粒子的邻居是动态变化的任何物理量的计算都要靠自己周围的邻居来插值。这就引出了 SPH 的核心工具——核函数。你可以把核函数理解为一种“带半径的加权平均模板”距离越近贡献越大超过一定半径直接为零。所有连续方程到这里都被改写成“对邻居求和”的形式这正是 SPH 名字里“光滑粒子”的由来。2.2 核函数与光滑长度 h把连续积分写成粒子求和SPH 的数学起点是一个恒等式任意连续场 A(r) 可以用核函数 W 表示成积分再用粒子离散化成求和A(r) ≈ Σ_j (m_j / ρ_j) A_j W(|r - r_j|, h)其中 W 就是核函数h 是光滑长度决定 W 的衰减半径。工程中最常用的是 Wendland C2 核因为它比传统三次样条更不容易产生粒子成团现象是近年多数 SPH 实现的首选。它的二维形式如下W(q) α (1 - q/2)^4 (2q 1), 0 ≤ q 2 q r / h α 7 / (4π h²) # 2D 归一化系数核函数必须满足三个性质归一化权重在支持域内积分等于 1对称性W 只取决于距离紧支撑距离超过 2h 权重直接为零。紧支撑决定了近邻搜索的半径是 2h也就是每个粒子只需要跟半径 2h 内的粒子通信这直接决定了性能模型。光滑长度 h 是一个全局敏感参数。h 太小支持域内邻居粒子太少权重分布抖动大密度场和压力场都会出现高频噪声h 太大所有权重都趋向均匀流体细节被抹平看起来像一团“缓慢膨胀的果冻”。常见做法是取初始粒子间距 dp 的 1.2 到 2.0 倍我一般从 1.5 起步具体影响在第 4 章展开。2.3 密度、压力与粘性最简 SPH 需要哪三项一套能稳定跑起来的最小 SPH至少包含三项计算。第一项是密度直接对邻居求和ρ_i Σ_j m_j W(|r_i - r_j|, h)注意求和要包含粒子自身因为 r0 时 W 取最大值这是密度公式里不可省略的自贡献项。第二项是压力弱可压缩 SPH 使用状态方程而不是求解压力泊松方程p_i (c0² ρ0 / γ) [ (ρ_i / ρ0)^γ - 1 ]这是 Tait 方程γ 通常取 7ρ0 是参考密度c0 是人工声速量级取预估最大流速的 10 倍。c0 越大流体越“硬”可压缩性越小但时间步长会被迫缩小计算成本上升。第三项是动量方程也就是加速度怎么算dv_i/dt -Σ_j m_j (p_i/ρ_i² p_j/ρ_j²) ∇W_ij Σ_j m_j Π_ij ∇W_ij g第一项是压力梯度力压力高的粒子把压力低的推开这是流体流动和波浪产生的主要驱动第二项是人工粘性项抑制粒子之间的高频相对运动防止数值震荡g 是重力。人工粘性 Π_ij 只在粒子相互靠近时激活作用类似物理粘性但它属于数值手段调整幅度要克制调大了会让流动显得“稠”。我不建议一上来就上投影法。虽然不可压缩 SPH 在精度上更优但每一步要解压力方程组实现复杂度上升一个量级。先用弱可压缩形式跑通物理流程确认参数和边界处理没问题再考虑换更高级的压力求解器工程节奏更合理。3. 最小可运行 SPH用 Python 写一个 2D 溃坝模拟选择溃坝作为首个场景原因很直接初始条件简单一个大水柱放在矩形区域一侧撤去支撑后水流自由坍塌物理现象丰富又是公开文献里的标准基准测试。下面这套基于 NumPy 的实现是我的教学惯用版本性能不是目标但结构完整能直接跑出肉眼可见的水流效果。3.1 参数与粒子初始化先把“水柱子”摆进计算域先用一组常量把全局参数固定下来再写初始化函数。粒子参数集中定义后续调参会很顺手import numpy as np # ---------- 模型参数 ---------- DP 0.02 # 初始粒子间距 0.02 m H 1.5 * DP # 光滑长度h/dp1.5 是常用起点 RHO0 1000.0 # 参考密度 kg/m^3 C0 20.0 # 人工声速按最大流速 2m/s 的 10 倍取 GAMMA 7.0 # 水的 Tait 状态方程指数 ALPHA 0.1 # 人工粘性系数 G np.array([0.0, -9.8]) # 重力2D 场景只有 y 分量 DT 2e-4 # 时间步长先保守取小 BOUNDARY_MARGIN 3 * DP # 边界斥力作用距离 K_BOUNDARY 500.0 # 边界斥力强度 def init_dambreak(water_w, water_h, domain_w, domain_h): # 在计算域左侧生成一块矩形水体 xs np.arange(0.5 * DP, water_w, DP) ys np.arange(0.5 * DP, water_h, DP) gx, gy np.meshgrid(xs, ys) pos np.stack([gx.ravel(), gy.ravel()], axis1).astype(float) # 加微小随机抖动打破规则排列 pos np.random.uniform(-0.05 * DP, 0.05 * DP, pos.shape) vel np.zeros_like(pos) mass np.full(len(pos), RHO0 * DP * DP) return pos, vel, mass参数里值得重点说明的是人工声速 C0。很多人误以为它应该取真实水中的声速 1500 m/s那样时间步长会被压到不可用的程度。SPH 中 C0 是数值参数取预估最大流速的 10 倍就能把密度变化控制在 1% 以内同时保证时间步长在合理范围。溃坝实验的最大流速常在两三米每秒取 20 是个稳妥起点。质量赋值 m ρ0 · dp² 对应的是 2D 场景每个粒子代表一个面积元而不是体积元。如果你改成 3D质量要换成 ρ0 · dp³。最后一行的 np.full 确保所有粒子质量一致这是 SPH 里“每个粒子质量相等”的常见简化——非均匀质量的 SPH 也能做但没必要在起步阶段给自己加复杂度。3.2 核函数与近邻搜索用 KDTree 拿到邻居点对接下来定义核函数和梯度。我用 Wendland C2它在粒子分布不均匀时的表现优于三次样条def wendland(q): # q r / h返回核函数数值 alpha 7.0 / (4.0 * np.pi * H * H) out np.zeros_like(q) mask q 2.0 qm q[mask] out[mask] alpha * (1.0 - qm / 2.0) ** 4 * (2.0 * qm 1.0) return out def grad_wendland(dx, q): # dx 为从粒子 i 指向粒子 j 的位移向量返回梯度向量 alpha 7.0 / (4.0 * np.pi * H * H) g np.zeros_like(dx) mask q 2.0 qm q[mask, None] r np.linalg.norm(dx[mask], axis1, keepdimsTrue) 1e-12 # dW/dr 的解析式 dq -5.0 * alpha * qm * (1.0 - qm / 2.0) ** 3 / H g[mask] dq * dx[mask] / r return g梯度公式里的符号要格外小心。dx 从 i 指向 jdW/dr 为负值距离越远权重越小两者相乘后梯度方向实际指向 j。这个方向约定直接关系到压力梯度的正负后面主循环里所有力的方向都以它为准。近邻搜索用 SciPy 的 cKDTree 实现。它一次性返回所有距离小于 2h 的点对省去自己写空间哈希的麻烦from scipy.spatial import cKDTree def build_neighbors(pos): tree cKDTree(pos) # output_typearray 直接返回 Nx2 的点对数组 return tree.query_pairs(2.0 * H, output_typearray)query_pairs 返回的点对不含“自己和自己的配对”这是符合预期的因为密度求和里粒子的自贡献需要单独补上。用 KDTree 的教学版本在几千到几万粒子时够用到了十万级我会换成均匀网格配合空间哈希原因在第 5 章的近邻搜索瓶颈里展开。3.3 密度与加速度计算主循环的每次迭代密度求和要特别注意自贡献项。每个粒子对自身密度的贡献是 m_i · W(0)W(0) 在二维 Wendland 下等于 7/(4πh²)必须手动加上否则密度整体偏低压力也随之失真def compute_density(pos, mass, pairs): n len(pos) rho np.full(n, 0.0) dx pos[pairs[:, 0]] - pos[pairs[:, 1]] r np.linalg.norm(dx, axis1) q r / H w wendland(q) # 对称累加每对邻居的贡献同时加到双方 np.add.at(rho, pairs[:, 0], mass[pairs[:, 1]] * w) np.add.at(rho, pairs[:, 1], mass[pairs[:, 0]] * w) # 自贡献项最容易漏 rho wendland(np.zeros(1))[0] * mass return rho密度求出后压力由状态方程直接计算然后组装加速度。压力梯度、人工粘性和重力依次加上def compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h): n len(pos) rho compute_density(pos, mass, pairs) p C0 * C0 * RHO0 / GAMMA * ((rho / RHO0) ** GAMMA - 1.0) acc np.tile(G, (n, 1)) # 重力项 dx pos[pairs[:, 0]] - pos[pairs[:, 1]] r np.linalg.norm(dx, axis1, keepdimsTrue) q r / H grad grad_wendland(dx, q) # 梯度指向 j # 压力梯度力 coef mass[pairs[:, 1]] * ( p[pairs[:, 0]] / rho[pairs[:, 0]] ** 2 p[pairs[:, 1]] / rho[pairs[:, 1]] ** 2 ) np.add.at(acc, pairs[:, 0], -coef[:, None] * grad) np.add.at(acc, pairs[:, 1], coef[:, None] * grad) # 人工粘性只在粒子逼近时激活 dv vel[pairs[:, 0]] - vel[pairs[:, 1]] prod (dv * dx[:, 0]).sum(axis1) if False else (dv * dx).sum(axis1) r2 r.squeeze() ** 2 0.01 * H * H mu H * prod / r2 rho_bar 0.5 * (rho[pairs[:, 0]] rho[pairs[:, 1]]) Pi np.where(prod 0.0, -ALPHA * C0 * mu / rho_bar, 0.0) np.add.at(acc, pairs[:, 0], mass[pairs[:, 1]] * Pi[:, None] * grad) np.add.at(acc, pairs[:, 1], -mass[pairs[:, 1]] * Pi[:, None] * grad) acc compute_boundary_force(pos, domain_w, domain_h) return acc人工粘性里的 prod 是 (v_i − v_j)·(r_i − r_j)当它小于 0 表示粒子正在靠近此时叠加一个与逼近速度成正比的阻力。mu 表达式里 0.01h² 是为了避免分母为零的数值保护也防止距离极近时粘性力发散。ALPHA 取 0.1 是经验值调大可以让模拟更稳但水花会显得偏黏。3.4 时间积分与边界斥力稳住的最后两环时间积分采用 kick-drift-kick 的 leapfrog 格式一次力计算可以推进半步稳定性和二阶精度都够用def run(pos, vel, mass, steps, domain_w, domain_h, save_every50): pairs build_neighbors(pos) acc compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h) for step in range(steps): vel 0.5 * DT * acc pos DT * vel pairs build_neighbors(pos) acc compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h) vel 0.5 * DT * acc if step % save_every 0: np.save(fpos_{step:06d}.npy, pos)边界处理是新手最容易翻车的位置。纯斥力函数虽然简单但只适合在低压场景用冲击壁面时压力可能把粒子从斥力区间里硬挤出去。我常用配合式内壁外侧放两层固定边界粒子让它们参与密度和压力计算流体粒子靠近时自然感受到压力升高而被推开而不是依赖解析斥力def compute_boundary_force(pos, domain_w, domain_h): f np.zeros_like(pos) margin BOUNDARY_MARGIN # 左墙粒子离左侧内壁越近向右的推力越强 dist pos[:, 0] - 0.5 * DP m dist margin f[m, 0] K_BOUNDARY * (1.0 - dist[m] / margin) # 右墙 dist (domain_w - 0.5 * DP) - pos[:, 0] m dist margin f[m, 0] - K_BOUNDARY * (1.0 - dist[m] / margin) # 底墙与顶墙同理沿 y 方向 dist pos[:, 1] - 0.5 * DP m dist margin f[m, 1] K_BOUNDARY * (1.0 - dist[m] / margin) dist (domain_h - 0.5 * DP) - pos[:, 1] m dist margin f[m, 1] - K_BOUNDARY * (1.0 - dist[m] / margin) return f把上面两个函数拼进 run一个最小的二维溃坝就能跑起来。初期你会看到水柱先微微松弛然后顶端开始翻卷底部水流加速冲向右侧壁面——这就是 SPH 对强自由表面流动的直接响应。若出现粒子乱飞先不要动核函数按第 5 章的排查顺序走一遍。4. 参数怎么定光滑长度、粒子间距、声速与时间步长的匹配SPH 调参没有“一套通吃”的黄金组合但参数之间有强耦合关系。我把最影响结果的四个参数单独拎出来讲清楚它们各自的职责以及调一个参数连带影响了什么。4.1 h/dp 比值从 1.2 到 2.0细节和稳定性的交换光滑长度 h 与粒子间距 dp 的比值是所有参数里最核心的一个它决定每个粒子支持域内的平均邻居数量。二维里邻居数约等于 (2h/dp)² 再乘一个密度因子h/dp1.5 时每个粒子大约有几十个有效邻居这个数量足以让密度和压力场平滑。h/dp 取 1.2 时支持域内粒子偏少压力场容易在局部出现空洞或尖峰好处是细节更锐利液滴脱离、薄水舌这类小尺度结构保留得更好。取 2.0 时数值更稳粒子成团倾向减弱但薄片结构会被过度平滑看起来偏“黏”。我的习惯是先取 1.5 跑基准再把 1.2 和 2.0 各跑一遍对比液面前沿位置曲线视差值决定往哪边靠。这个对比实验成本很低但能避免拍脑袋定 h。还有一个容易忽略的点h 是全局固定值还是可以随密度自适应自适应 h 能改善粒子稀疏区域的表现但会让核函数失去对称性动量守恒变差。除非粒子分布跨度极大否则我不会在第一个项目里用自适应 h。4.2 粒子间距与规模感从 1 万到 100 万粒子怎么换算粒子间距 dp 决定分辨率也直接决定总粒子数。二维中粒子数 N 与面积成正比三维中与体积成正比dp 减半时三维粒子数变成 8 倍。拿一个 1m×0.6m 的二维水槽估算dp0.02m 时大约 1500 个粒子演示没问题dp0.005m 时直接跳到 24000 个计算量和内存同步上涨。我的经验是先把规模分成三个档位来规划。教学和原理验证在几万粒子内用 CPU 单线程配 NumPy 足够十到几十万粒子的工程测试需要上 C/Rust 配合均匀网格近邻搜索百万级场景基本绕不开 GPU 或大规模并行。下表给出我自己常用的规模参照粒子规模常见场景计算侧重点内存压力1 万以下算法验证、参数初调近邻搜索逻辑正确性极低1 万到 10 万小溃坝、液滴撞击密度与压力场稳定性中等10 万到 100 万工程仿真、动画生产近邻搜索结构、并行化较高100 万以上大场景工业仿真GPU 显存调度、多卡分解极高这里的“内存压力”不是指粒子数组本身——一个二维粒子含位置、速度、密度、压力约 5 个 float64 数组每个粒子 40 字节三维会更多压力主要来自近邻点对数组它的长度能到粒子数的几十倍才是内存消耗的大头。4.3 人工声速与 CFL 条件压缩性和时间步长的拉扯弱可压缩 SPH 里人工声速 c0 决定了流的“軟硬程度”。c0 取 10 倍最大流速对应马赫数 0.1密度振荡控制在 1% 量级。提高 c0 可以进一步逼近不可压缩但声速越高信息传播越快显式时间积分被迫用更小的 dt计算量随 c0 近似线性上涨这笔账要把模拟时长一起算进去。时间步长受 CFL 条件约束常用判据是dt ≤ 0.25 * h / (c0 v_max)v_max 是当前最大流速。这个式子能解释大部分“模拟爆炸”事故水柱刚倒下时流速快速增大如果 dt 按初始流速算好就不再调整一旦瞬时流速超过设计值CFL 条件被破坏局部分离会在一两个时间步内指数放大。我把常用的起步参数整理成一张速查表方便照抄参数起步值调大影响调小影响h/dp1.5更稳、细节模糊细节保留、易震荡c010 倍 vmax更不可压缩、dt 变小更可压缩、压力偏软γ7.0压力响应更硬压力响应更软ALPHA0.1耗散强、流动发黏易出现高频抖动DT0.25h/(c0vmax)可能触发爆炸计算变慢参数之间存在联动比如 c0 调大一倍CFL 条件里的 vmax 贡献相对变小但绝对限制更严DT 往往要同步减半ALPHA 调大后高频噪声被压住看似稳定了但真实的细小涡旋也被抹掉。因此调参时要有“一次只动一个”的纪律并且永远用同一个基准场景对比。5. SPH 模拟的 4 个常见坑从粒子爆炸到边界穿透SPH 的坑很有辨识度物理引擎报错会告诉你“穿模了”“飞出去了”SPH 则是直接给出一堆 NaN 或者一帧比一帧离谱的粒子云。下面四条是我自己踩过也帮别人排查过的高频问题按现象、原因、解决三步写。5.1 粒子炸飞密度尖峰后的压力雪崩现象模拟跑到几十步某几个粒子速度突然变成几百甚至几千随后整个粒子场在几步内全面发散数值变成 NaN。断点调试时能看到爆炸前一个时间步的密度数组里出现异常大值。原因压力状态方程里 (ρ/ρ0)^γ 的 γ7密度一旦超过参考密度压力以七次方放大。密度尖峰最常见的来源是粒子初始抖动过大导致支持域内粒子重叠或者是 dt 太大导致相邻粒子在一步内穿过彼此。自由表面处粒子密度偏低产生负压又把周围粒子吸附过来形成局部坍缩。解决先把 dt 减半看是否还炸如果还炸检查初始化抖动幅度把 0.05dp 改成 0.01dp再不行就检查核函数支持域内是否有粒子缺失比如水柱角点处粒子天然稀疏此时需要稍微加大 h。最后的手段是给压力做一个上下限截断但那是治标只能让模拟“不炸”不能让你理解物理不推荐新手优先用。5.2 粒子扎堆拉伸不稳定性在捣乱现象模拟中后期粒子在流体内部或自由表面附近聚成小团形成规则晶格状的分布密度场出现与初始网格结构一致的花纹。这不是真实物理是 SPH 的经典数值病——拉伸不稳定性。原因在受拉伸区域负压环境SPH 的力对粒子位置扰动呈正反馈规则排列被不动点吸引粒子趋于成团。早期三次样条核函数更容易触发这个现象这也是我改用 Wendland 核的原因之一。粒子间距相对 h 过密时也更明显。解决优先把 h/dp 从 1.5 提到 1.8 或 2.0让权重更平滑再把人工粘性略微提高抑制局部相对速度。如果还不行检查 c0 是否过高——c0 越大压力对密度误差越敏感负压抖动越剧烈。换成更高阶的核函数或加密度修正项也能缓解但代码复杂度会上去建议作为第二步。5.3 边界穿透斥力挡不住高速冲击现象水柱坍塌后高速冲击底部或侧壁一部分流体粒子直接穿出计算域卡在边界外。用固定斥力函数的模拟尤其容易出现因为斥力只在粒子进入边界的某个区间内生效速度足够大时一个时间步内就能穿过整个区间。原因纯斥力的本质是“事后补救”它不参与流体的密度和压力计算相当于一道没有物理支撑的墙。粒子速度接近或超过 margin/dt 时斥力还没来得及起作用粒子已经越过边界。解决最可靠的做法是改用固定边界粒子。在域边界外侧布置二到三层不参与积分的虚拟粒子让它们正常参与流体粒子的密度和压力求和。流体粒子接近边界时密度真实升高压力真实增大自然被推回——这比任何解析排斥力都符合 SPH 的物理逻辑。边界粒子层的厚度至少要覆盖核函数支持域三层是稳妥起步。5.4 近邻搜索拖垮性能O(N²) 的死胡同现象粒子数从 1 万翻到 4 万运行时间不是涨 4 倍而是涨将近 10 倍。用 KDTree 的教学代码在粒子数上升后每步重建树的开销越来越大如果还为每个粒子单独查球就更慢。原因近邻查询的复杂度与支持域半径、粒子密度相关但实现不当很容易退化成全局搜索。每步重建 KDTree 本身也有常数开销在十万粒子量级会明显拖后腿。解决换均匀网格空间哈希网格尺寸直接取 h每个粒子只需要查询自己所在网格及周围 3×3二维或 3×3×3三维的网格。网格大小和近邻搜索只遍历邻近网格复杂度降到 O(N)这是 SPH 工程实现的标配。静态边界粒子可以预先建好网格每步只重建流体粒子部分。6. 往工程走GPU 加速、验证实验与参数标定的习惯6.1 三个验证实验静水压、溃坝、能量曲线跑通模拟只是第一步真正判断“做对了”要依赖验证手段。我习惯先做三个实验静水压测试、溃坝对照、能量曲线检查。静水压测试很简单把水体放在容器里静置几千步检查底部压力是否符合 ρgh 的理论值且粒子速度长期保持接近零。这个测试能暴露核函数、状态方程和边界处理的系统性问题比一上来就跑动态场景省心得多。溃坝对照是用公开实验数据里水柱前沿位置随时间的变化跟自己的模拟结果画在同一张图上看趋势是否一致这是判断整体参数是否合理的金标准。能量曲线则是记录系统动能加势能E_kin 0.5 * (mass * np.sum(vel ** 2, axis1)).sum() E_pot (mass * (-G[1]) * pos[:, 1]).sum() E_total E_kin E_potSPH 的总能量不严格守恒因为人工粘性会耗散能量但曲线应该平滑下降而不是上下剧跳。如果总能量出现明显回升多半是时间步长或边界处理有问题值得回头排查。6.2 从 CPU 到 GPU性能上台阶的正确姿势粒子规模到了十万级CPU 版本的速度会让人失去耐心。常见的加速路径是把近邻搜索和密度、压力、力的累加写成 GPU kernel同一份数据在 GPU 上的计算速度通常比 CPU 高一个数量级。现有开源 SPH 框架大多提供 CPU 和 GPU 双后端但我不建议直接跳到麻烦的 framework——先用第 3 章的最小实现理解每个阶段的数据流动再迁移到并行版本遇到瓶颈时你能清楚知道是近邻写入冲突还是内存带宽不足。6.3 参数标定的第一个习惯我把压箱底的习惯写在这里每一个需要调试的项目永远先建立一个“最小复现场景”它要足够小到几十秒内能跑完又足够典型到能复现你关心的物理现象。参数全部放进配置字典每次只改一个记录下它前后的能量曲线和液面前沿别靠肉眼判断“看起来差不多”。我吃过不少亏最典型的一次是同时调了 h 和 ALPHA结果模拟从炸飞到死水都找不到原因最后只能重跑分项排查。血泪经验换来的教训就是参数实验是科学实验要有对照和记录。希望帮到你。本文还有配套的精品资源点击获取