做了这么多年聚变等离子体模拟电脑里装过的输运代码一只手数不过来t3code是其中资历最老、却最不该被忽略的一个。如果你盯着“t3code”这个名字发愁多半是在查托卡马克等离子体的一维径向输运演化怎么做或者看到某篇老文献里引用它作为基准工具想拿它复现一个密度剖面或温度剖面却不知道从哪下手。T3不是那种开箱即用、文档齐全的现代框架它更像一台做工扎实的老式机床——没有花哨界面但物理建模思路非常通透跑一个典型场景只需要几秒到几分钟特别适合做参数扫描、概念设计、以及快速验证“如果我把加热功率抬高到多少台基会不会塌”这类问题。这篇文章我按自己的使用经验把t3code从物理模型、编译运行、算例设置到常见坑点完整梳理了一遍给准备入门或被迫接手老代码的你省点时间。1. t3code到底解决什么问题一维径向输运的取舍1.1 为什么聚变模拟偏要选一维托卡马克等离子体真实状态是三维的有环向、极向和径向三个空间方向再加上时间是个四维系统。但你真想研究“密度分布怎么随时间演化”“温度台基怎么形成”直接上三维湍流代码是杀鸡用牛刀——成本高、周期长且大量细节对你要回答的问题并没有帮助。t3code做的是把问题简化成一条径向线。它假设等离子体在磁面上是均匀的所有物理量只随小半径r变化再用通量面平均的方式把三维效应塞进一维系数里。这样做的核心逻辑是输运过程主要由径向梯度驱动环向和极向的快速输运会把不均匀抹平所以描述长时间尺度演化时“径向一维时间演化”是一个非常好的近似。这个取舍的收益非常直接运行速度快、参数空间容易扫、稳定性和可复现性都远好于大代码。用T3做一个基准算例从给定初始剖面到接近稳态通常就是几分钟的事情。它不会给你湍流的细节但能给你整个等离子体的“骨架演化”这在做设计参数扫描和趋势判断时足够用了。1.2 时间尺度分离欧姆定律是准静态的T3背后有个关键物理思想时间尺度分离。等离子体里有极快的磁流体阿尔芬过程有微秒量级的湍流涨落也有毫秒到秒量级的能量约束时间还有更慢的电流扩散和杂质积累。如果我们把全部过程放在一起算数值步长必须追到最快过程计算成本会爆炸。所以t3code只追踪慢过程等离子体密度、电子温度、离子温度、电流密度剖面、杂质浓度和极向磁场的演化走“输运方程组”而快过程如磁流体平衡、粒子源和加热沉积则用准静态平衡和经验模型即时更新。实际实现中这意味着每个时间步内T3会先调用一个平衡计算基于当前的电流剖面和压强剖面当场解出磁面几何更新通量坐标然后才推进输运方程。这是一个典型的“算子分裂”思路快物理作为背景慢物理作为演化主体两步交替前进。如果你读T3源码会发现主循环里就是“求平衡-推进输运-更新源项”这个结构在不断重复逻辑清楚到即使Fortran老代码也不难看懂。2. 核心物理模型与代码结构速览2.1 输运方程组粒子、热、电流一条龙T3求解的核心方程组可以分成三块都源自标准通量面平均输运理论。密度演化方程写为粒子平衡∂n/∂t -(1/V)(∂/∂r)(V Γ) S_n其中V是体积对径向坐标的导数Γ是径向粒子通量S_n是粒子源项包括中性束注入、壁解吸和再循环等。Γ由扩散系数D和对流速度V convective构成Γ -D ∂n/∂r V_conv n。电子和离子的能量方程则是热平衡(3/2)∂(nT_e)/∂t -(1/V)(∂/∂r)(V q_e) P_heat_e - P_eq - P_rad (3/2)∂(nT_i)/∂t -(1/V)(∂/∂r)(V q_i) P_heat_i P_eqq是径向热流P_heat是外加加热功率密度P_eq是电子离子碰撞能量交换项P_rad是辐射损失。真正精细的地方在于P_rad通常包含轫致辐射、线辐射和杂质辐射模型其中杂质含量一旦涨起来辐射会把边缘温度压得很低。电流密度演化则是一个磁扩散方程∂B_θ/∂t (η/μ₀)(1/r)(∂/∂r)(r ∂B_θ/∂r) 驱动项这里η是电阻率由Spitzer公式给出驱动项包括欧姆电流、束驱动电流和自举电流。自举电流的模型直接决定了剖面是否“空心”也影响最终是否进入稳态。这三套方程用径向网格离散然后时间推进。T3默认推进方式是显式的优点是简单、不需要大型矩阵求解缺点是时间步长会被稳定性条件卡住。所以跑长时间演化时注意时间步长设到满足库朗条件具体经验我在后文会讲。2.2 径向网格与边界条件近轴点最容易出问题T3的径向网格一般是等间距或略微加密的近轴分布。网格数不用太多一二百个点已经能描述剖面形状多了反而会拖慢速度。反正它是一维代码单步计算成本本来就低。边界条件分两个端点。中心r0处所有通量归零所有物理量的径向梯度为零——这是对称性条件不能乱改。边界ra处也就是等离子体边缘温度密度要被设定为有限的边界值。在实践中T3的边界条件通常由刮削层模型或直接给定边界温度和密度表格来提供。踩坑提示多数运行发散问题都出在边界条件不自洽。你设的边界温度低于辐射和输运决定的自然边界值边缘立刻出现陡梯度导致热流数值振荡。如果你看到T_e剖面在最后一个网格点上跳来跳去十有八九是边界值不匹配先调边界比调时间步长更管用。2.3 输运系数新经典基础加反常提升输运系数是整个模拟的“灵魂”T3里通常拆成两部分写。新经典输运部分基于Hinton-Hazeltine模型处理碰撞导致的不连续轨道效应包括香蕉区、平台区和箍缩区的输运。这一项在高温低碰撞率区域较弱但它是自洽的由等离子体参数直接推导出来不需要额外调节。反常输运部分则是“调参的重头戏”。T3默认提供几种参数化模型比如定常D、定常χ以及更复杂的缩放模型。实际操作中绝大多数人都会在D和χ上乘一个系数去匹配实验约束时间。我的做法一般是先跑一组基准给定密度和加热调χ_e和χ_i让温度剖面匹配实验点然后再用这个标定过的输运系数去扫描其他参数。需要警惕的是参数过拟合。T3这种参数化输运代码总是能通过调系数拟合出一两组实验剖面但到了预测新场景时系数不一定可迁移。所以做预测前至少要用两三个不同放电的基准点做交叉验证不要只靠一个放电定系数。3. 从拿到t3code到跑出第一张剖面图的完整流程3.1 获取源码与编译环境准备t3code是老牌Fortran代码很多研究所内部都有本地拷贝如果你能访问项目仓库或代码分发渠道直接拿到源码包最省事。拿到后先看压缩包里的文件结构一般包含一个主程序文件、若干物理模块子程序和一个示例输入文件还可能附带一个用于绘图的后处理脚本。编译前确认你装了Fortran编译器。这里我推荐gfortran免费、跨平台、对老代码兼容性好。Intel的ifort也行有些老代码针对ifort做过向量化优化但如果没有licence就没什么好纠结的。在Linux环境下一条命令就能完成编译gfortran -O2 -o t3code t3main.f90 t3transport.f90 t3balance.f90 ...如果你的源码包里有Makefile就直接make编译报错一般是两类一类是Fortran语法版本差异老代码用的固定格式.f后缀和自由格式.f90后缀混在一起gfortran默认可能按后缀区分格式注意把固定格式文件加上 -ffixed-form另一类是缺少某些工具库或数学函数库引用用 -lm 或添加对应链接路径即可解决。一个务实的建议先用小网格、最少物理模块编译通过再逐步打开完整物理选项。一次全开报起错来就要排查很久。3.2 输入文件结构namelist风格要耐心t3code的输入控制通常走NamellistFortran的命名列表方式你会在示例输入文件里看到类似下面的结构name_list nr 128 rmax 0.85 dt 0.01 ntstep 10000 ne_bc 1.0e19 te_bc 100.0 ti_bc 100.0 ipedestal 1 chi_e 1.2 chi_i 0.8 d_coef 0.5 /每个参数都有明确含义。nr是径向网格数rmax是等离子体半径米dt是时间步长秒ntstep是总时间步数ne_bc、te_bc、ti_bc是边界电子密度、电子温度、离子温度ipedestal开关决定是否启动台基高度模型chi_e、chi_i、d_coef分别控制电子热扩散系数、离子热扩散系数和粒子扩散系数的基值。第一次跑请务必先保留所有参数默认值只改输出文件名。T3内置的基准算例大概率能直接跑通给你一个已知正确的参照。先建立“能跑”的基线再做物理修改这个顺序能帮你避免“改了十个参数出了问题却不知道是哪一个改坏的”的尴尬局面。3.3 正式运行与标准输出解读编译和输入文件准备完毕后运行非常简单./t3code t3code_input.dat终端会实时打印当前时间步、等离子体电流、约束时间、中心温度等关键指标。如果你看到这些量稳定增长并趋于饱和恭喜这个算例走上了正轨。T3的每个时间步都会把剖面数据写入输出文件一般是纯文本格式多列分别为径向位置、电子密度、电子温度、离子温度、电流密度等方便你用任何工具处理。我一般会把终端输出重定向到log文件方便逐段回看./t3code t3code_input.dat run.log 21如果中途跑飞数值发散run.log里通常能看到时间步推进失败或剖面出现NaN的痕迹这是第一手排查资料丢掉很可惜。3.4 把输出文件画成剖面演化图T3的历史输出文件记录的是“每一小间隔时间步的全剖面快照”。文件如果很大画图前要先抽样否则gnuplot或matplotlib都能被卡死。用Python清洗和绘图可以参考如下脚本思路import numpy as np import matplotlib.pyplot as plt data np.loadtxt(t3code_profiles.dat) # 假设第1列是时间轴第2列是r第3列是Te tvals np.unique(data[:, 0]) fig, ax plt.subplots(figsize(8, 5)) for t in tvals[:: max(1, len(tvals)//10)]: # 抽样10个时刻 mask data[:, 0] t ax.plot(data[mask, 1], data[mask, 2], labelft{t:.2f}s) ax.set_xlabel(r (m)) ax.set_ylabel(Te (eV)) ax.legend() plt.tight_layout() plt.savefig(te_evolution.png, dpi150)如果你拿到的输出格式是“每个时间步一个文件”的形式那就更简单了循环读取文件绘制即可。画图时注意横轴是“小半径”还是“归一化半径”T3默认常见输出是物理半径画图前记得看单位避免出现一根看起来像是被压缩过的奇怪曲线。3.5 参数扫描实例加热功率对温度剖面和约束时间的影响用t3code做参数扫描是它最擅长的事情。现在演示一个典型任务扫描总加热功率从5MW到25MW看中心电子温度和能量约束时间怎么变化。做法很简单写一个小循环每次改输入文件里的加热功率参数调t3code抓取输出里的中心温度和约束时间。我习惯用Python脚本控制这个过程import subprocess import numpy as np powers np.linspace(5.0, 25.0, 9) results [] for p in powers: with open(t3code_input.dat, r) as f: text f.read() text text.replace(p_heat 10.0, fp_heat {p:.2f}) with open(t3code_input.dat, w) as f: f.write(text) subprocess.run([./t3code t3code_input.dat scan_run.log 21], shellTrue) # 从输出文件读取中心温度和约束时间 data np.loadtxt(t3code_profiles_last.dat) results.append((p, data[-1, 2], data[-1, 0]))跑完之后把结果画出来你会看到中心温度随功率上升、能量约束时间却往下掉——这正是典型L模约束定标率的体现。如果你给代码开了台基模型可能还会看到功率超过某个阈值后台基高度跳变、约束时间突然变化的“相变”行为。这种参数空间里的拐点用大代码扫描是奢侈的用t3code可以轻松探索。4. 实战中经常翻车的几个点问题定位与修正4.1 初始剖面的选择冷启动还是热启动T3的初始剖面直接决定前几百步的稳定性。最常见的翻车是初始剖面与世界平衡差别太大第一大步的梯度就超出了显式时间步的稳定性极限直接振荡。解决方式有两条路。第一使用更平滑的初始剖面比如从抛物线形状开头不要上来就塞一个很尖锐的台基状剖面。第二减小初始时间步长让代码自己慢慢适应通常跑过前几百步后就可以适当增大步长。T3允许你在输入文件里做时间步长的分段控制可以设定前500步用dt0.001之后切换到dt0.01这样的“预热”能有效避免冷启动发散。如果你的研究对象是台基演化那就必须从高约束模开始做热启动。给一个事先跑好的稳态剖面作为初值更稳妥。T3支持读入上次运行的最终剖面作为新运行的初始条件这个功能强烈建议用起来别每回都从同一个粗糙初值算起。4.2 显式时间步长、网格数与源项增量的三角约束T3采用显式推进时间步长受限于网格分辨率与最大输运系数的乘积。公式上大概满足dt ≤ Δr² / (2 max(D, χ))如果你把网格数从128加到256时间步长理论上要缩小4倍才能维持稳定。实操中我遇到过有人只加密网格、忘改时间步长然后五分钟前还在正常跑的算例突然就爆了。排查这种问题先看run.log如果发散点固定出现在同一个时间步附近并且剖面在某个径向位置开始锯齿状振荡基本可以判定是稳定性条件被破坏。对策是缩小dt或者用输运系数更保守的参数值。源项增量也得注意加热功率密度在网格边界处跳变过大会让对应网格点的温度更新过猛局部产生虚假热流。更平滑的沉积剖面能缓解这个问题。4.3 台基高度与边界温度两个容易混淆的调参旋钮很多T3的中文用户上来就调“台基参数”但台基模型效果不明显就开始怀疑代码。这里有个容易混淆的点T3这类代码里的“台基”模型一般说的是控制台基高度或宽度的经验关系式它在漂移动理方程里不是自觉出现的而是被“钉”在某个值附近。真正的边界温度则直接影响台基底部的位置。边界温度给得太低台基难形成给得太高台基会异常高。我的经验是把边界温度作为一个“锚点”配合台基标度关系一起调。先固定边界温度不变量扫描台基模型参数找到匹配实验高度的数值再回头微调边界温度让台基底面位置也match上。两步分开做不要一起动否则你永远不知道是谁在起主导作用。4.4 杂质含量与辐射损失隐藏的剖面“刺客”如果你跑参数扫描时发现温度剖面整体偏低特别是边缘温度上不去而密度剖面看起来正常先别急着调热输运系数——看一眼杂质浓度参数。T3支持输入有效电荷数Z_eff剖面或固定杂质浓度辐射功率密度跟Z_eff近似成正比Z_eff从1涨到3辐射损失能吃掉一大块加热功率边缘温度应声而落。这是个非常容易“误诊”的问题。我见过有人为了抬高边缘温度费劲半天把chi_e调了五倍结果发现只是辐射模型开关忘了关。建议每次换工况时都先确认杂质相关的参数和辐射模型开关是否符合你的物理预期再谈其他调整。5. 快速排查速查表症状可能原因排查方向运行几步后出现NaN时间步长过大缩小dt检查初始剖面梯度边缘温度振荡边界温度与输运不平衡调边界温度或输运系数密度剖面不变化粒子源项或扩散系数过小检查S_n和D的取值量级温度剖面整体偏低杂质辐射过大查Z_eff、辐射模型开关台基高度异常边界温度、台基参数不匹配固定一个变量单独扫描另一个长时间跑不到稳态电流扩散时间尺度远大于热输运时间尺度跑足够长时间的模拟或干脆用更小的初始dt预热编译报格式错误老代码固定格式/自由格式混乱用-ffixed-form或调整后缀名查问题时我习惯只改一个变量跑短时间步对比再改下一个。T3跑得快试错成本低不要一上来就大改参数组合那是把简单问题复杂化。6. 从t3code往后走扩展思路与个人体会t3code这类代码的价值在框架而不在某个参数的精确值。跑通它、看懂它等于把一个完整的一维输运思路装进脑子里。后续如果要做更复杂的分析完全可以以T3的输运方程和推进框架为蓝本加入自己的模型比如新的反常输运系数表达式、更精细的边界模型、或者耦合一个简单的刮削层模型。我在实际项目里就用T3做过沙箱测试把一种新的台基高度标度式写进去对比基准算例确认趋势合理后才决定要不要移植到大代码里跑高成本模拟。还有个很容易被忽略的小技巧T3速度快非常适合作不确定性量化。你可以在输入参数上加随机扰动跑几十次ensemble看输出剖面的散布范围。这种事在TRANSP这种级别的大代码上几乎不敢想但在t3code上只是多花几分钟的事。如果你是被课题或项目“强制”需要弄懂这个代码的建议你从基准算例开始逐行过一遍主循环搞清楚“平衡—输运—源项”三个模块之间的调用顺序。这个顺序搞懂了T3就已经是你手里的工具而不是一团黑箱。后续想往里面加物理或者换数值格式都有明确的下手点。