4 天 215 次提交 91 张缺陷单：把 LAMMPS 改成能跑 PFC 脚本的内核
先说结果。
过去 4 天，我们让开源离散元引擎 LAMMPS 能原样跑通 Itasca PFC2D 6.0 的脚本，并用 6 个真实工程算例做了两侧对照验证。为此在 LAMMPS 的 fork 上提交了 215 次代码（涉及 157 个文件），开了 91 张缺陷单，其中 79 张已修复并通过独立复验。同一份 PFC 脚本，分别喂给 LAMMPS 和 PFC 6.0，两侧都能跑完、结果对得上。
这一篇不讲技术细节，只讲这件事到底花了多少功夫、功夫花在哪。
一句话背景：PFC 是商业软件，脚本生态成熟但闭源、按 license 收费；LAMMPS 开源免费、算力可扩展，但它不认识 PFC 的脚本。我们的目标不是"用 LAMMPS 仿一个类似的东西"，而是让一份没改过的 PFC 脚本，原封不动喂进 LAMMPS 就能跑。
1. 为什么要做这件事
用户手上的工程积累，全是一份份 PFC2D 的 .dat 脚本——成样、加胶结、换边界、加载，一条多步链跑下来。这些脚本是多年沉淀的资产。
问题在于 PFC 本身：license 按核数收费，大规模并行的成本随算例涨；闭源也意味着遇到不符合预期的地方，只能等原厂。
LAMMPS 恰好相反——开源、免费、并行扩展性好，但它原生是分子动力学的词汇体系，看不懂 ball distribute、contact cmat、wall create vertices 这些 PFC 命令。
于是有了这个项目：把 LAMMPS 改造成一个能执行 PFC2D 脚本的内核。改造完成度不能用"我说改好了"来证明，只能用真实算例来证明。
2. 验证方法：同一份脚本，两边各跑一遍
验证的口径很简单，但执行起来很严：
从真实工程算例里挑脚本，按 PFC 脚本的原始顺序跑； 同一份脚本， lmp_2d（我们的 fork）与 Itasca PFC 6.0 各跑一遍；母本两侧逐字节相同——脚本改一个字都不行，改了就等于把"内核差异"和"脚本差异"混在一起，结论作废； 比对宏观量：加载曲线形状、峰值强度、应力水平、接触数、能量、收敛步数。
判据有个前提要先说清楚：两侧的颗粒是各自随机生成的，分布必然不同，所以逐颗粒比对没有意义，跨侧的绝对值差只能作量级参照。真正立得住的，是曲线形状、几何形态，以及同侧内部的统计量。
这条口径是我们踩了坑才定下来的：早期我们拿"两侧峰值强度相差几倍"当缺陷报，后来发现换一个 model random 种子，同一侧自己的峰值就能散布 2.4 倍——那个判据根本判别不了内核差异。差多少算"对得上"，得先把判据本身测到站得住。
3. 六个算例
下面每一个算例，都是两侧各自成样、各自加载，跑完后把末态颗粒场画出来对照。
3.1 1. 巴西劈裂（平行胶结）
圆盘试样，径向加载，测间接抗拉强度。
图1 巴西劈裂 —— 末态颗粒场
球数 321 / 306。劈裂面形态一致。
⚠️ 这个算例的峰值之比同样不能当判据——换种子后巴西峰值 σt 在 1.52e7~3.67e7 之间散布，还出现过一个"试样压根没被压到"的废点。但图仍然有用：形状与回弹形态可比。
3.2 2. 直剪（平行胶结）
剪切盒里施加法向压力，再推剪。
图2 直剪 —— 末态颗粒场
球数 160 / 159。剪切的位移场与破坏位置两侧对应。
3.3 3. 地基开挖
分层开挖、地表沉降。这个算例的几何最大——球心分布横跨约 39 m。
图3 地基开挖 —— 末态颗粒场
球数 327 / 319。两侧挖方量相近、坑几何相同。
⚠️ 这个算例有一条还没查清的差异：LAMMPS 侧持续下沉到 -0.0407 m，PFC 侧在第 6 个开挖块后冻结在 -0.0163 m，差约 2.5 倍。成因未定，如实记在案。
3.4 4. 边坡滑动（课程第四章）
先重力沉降成样，再用 range plane 切出 40° 坡面。
图4 边坡削坡后 —— 末态颗粒场（40° 坡面）
球数 8242 / 8015。削坡的几何高度一致：左高台 → 40° 斜坡 → 右低台，坡脚坡顶位置几乎重合。切坡删除占比 17.70%（LAMMPS）vs 17.59%（PFC）。
3.5 5. CPT 贯入（第三章）
CPT 锥（wall create vertices 建的两条 facet）以恒定速度贯入 2 m，记录贯入阻力。
图5 CPT 贯入末态颗粒场（锥已贯入 2 m）
球数 9948 / 9683，锥的几何两侧逐位吻合（差 < 1e-4）。
这个算例曾经完全跑不动：model restore 之后 wall create vertices 被拒，报 Walls can only be created before cycle 0，而 PFC 正常建墙。根因是建墙守卫读的是全局步数，而 restore 会把快照里的步数带进来（这一步本身是对的），于是守卫把"快照带进来的步数"误读成"本次运行已经 cycle 过"。修法是改判"此刻是否正在 cycling"。
修完之后的贯入阻力曲线：
图6 CPT 贯入阻力曲线（两侧对照）
而且我们发现曲线起点不一样——LAMMPS 侧要到 0.6 m 才开始受力，PFC 侧 0.1 m 就起了。追下去是成样随机性：脚本用 ball delete range position-y 按硬阈值 5.175 削顶，锥正上方那一列有没有球恰好高于阈值全凭随机。实测两侧锥正上方那列试样顶面差 0.36 m，于是 LAMMPS 侧的锥起始悬空 0.29 m、PFC 侧已嵌入 0.08 m——接触时刻天然不同。
3.6 6. 双轴压缩（第七章）
最完整的算例：成样 → 伺服围压 → 加平行胶结 → 换承压板 → 加载到 3%。
图7 双轴压缩 —— 破坏后颗粒场 + 裂纹（红点）
球数 10020 / 10006。两侧都形成了斜向剪切带，但形态有别：PFC 侧裂纹聚得更紧、呈明显的交叉带状；LAMMPS 侧更弥散铺开（断键数 5211 vs 2695）。
⚠️ 断键数是跨侧量，只作量级参照。
这个算例还暴露了一类最难查的缺陷：跑得完、也不报错，exit code 全是 0，但数字是错的。比如——
平行键的断裂模式恒报"拉伸"，剪切永远读不到（因为 svector[13]在断裂被观察到时读到 0，落回默认值）；读平行键刚度 pb_kn时返回 0，而 PFC 返回1e9（根因是 deform 参数化下从没设过这个成员）；建键当场 bonded=1797，一存盘就掉到354，而 PFC 一个不掉。
这类"静默读错值"比崩溃危险得多——崩溃会告诉你它坏了，静默错误不会。
4. 91 张单长什么样
修复不是"改到不报错"就算完，每一条都要走完这套流程：
从真实算例里撞出问题，做出最小复现（常常缩到 1 个球）； 两侧对照，确认是缺陷而非脚本问题； 写成缺陷单，含现象、复现脚本、两侧日志； 派给开发、修复后由测试侧独立复验——自己克隆、自己编译、跑反向对照（把修复撤掉，该红的面必须红）。
最后一条是关键。"我不能把它测成红的"和"它确实是绿的"是两回事，而复验这件事必须由没有写这段代码的人做。
91 张单里，问题大致分几类：
wall create vertices restore 后被拒、brick import 整条命令缺失 | ||
pb_state pb_kn 读 0 | ||
wall.facet.* 名字，fork 一个都没有 | ||
不是所有"对不上"都是 bug。 91 张里有 1 张最终判定为"不是 bug"、8 张经用户裁定搁置（比如需要新增接触本构模型这类难度较大的），目前还有 3 张在办。
5. 钱花在哪了
这 4 天，光 token 就烧了 350 块。
听起来不少。但要记住这是一套并行体系在跑——一个测试工作区加三个开发工作区，全部围绕同一个缺陷单信箱协作。摊到每张已结案的缺陷上，大约 4 块 4。
更重要的是这些钱换来了什么：如果靠人工，光是"把 91 个问题各自缩成最小复现 + 两侧对照 + 独立复验"这一套流程，就远不止这个数，而且人的耐心会在第 20 个"静默读错值"上耗尽。
成本结构上值得记一笔：单次复验大约 17 万 token，大头不在"跑"，而在重建产物与冷启动。所以我们后来立了规矩——复验前先查现成的快照表，避免重复编译；同一提交的两臂对照能复用快照就复用。
6. 结语
这次工作的产出不是一个"能跑"的内核，而是一套可验证的内核：91 张单是它的体检报告，79 张已结案是它的合格证明，6 个真实算例是它的实测成绩。
还开着的口子也如实列在这里：
地基开挖沉降差 2.5 倍，成因未定； 双轴弹性段模量跨侧差 23.3%，已排除"单个接触刚度"这一层，嫌疑在"接触网络怎么组织"； 平行键的断裂模式仍报不准（拉断/剪断分不开）； wall.facet家族还缺 8 个名字。
诚实声明：本文所有跨侧数值差都只作量级参照。两侧颗粒分布本就不同，逐颗粒对比无意义；能立得住的判据是曲线形状、几何形态与同侧内部统计量。上面每一个结论，都有对应的日志与复现脚本可查。
VibePFC 社区 · 一起把 AI 玩进 PFC
我成立了一个 VibePFC 社群，专门聊 AI 怎么落地到 PFC / 数值模拟——让大模型读懂模型、自己改脚本、自动排查"看不懂"的问题。收一个 50 元的门槛费。
想进群的同学，加 QQ 763388012，备注"vibePFC"即可。
