fmq/TLUSTY
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
TLUSTY/cno_grid/PIPELINE.md
fmq 92d9c515c9 2
2026-07-23 19:06:45 +08:00

20 KiB
Raw Permalink Blame History

CNO 网格完整计算流程

本文档说明完整理论光谱网格的计算流程：每个网格点的计算阶段、每阶段的配置、 配置原理、CPU/并行机制、以及如何统计每个阶段的信息（时间、收敛等）。

1. 总体架构

config.yaml (网格点 + 收敛链配置)
       │
       ▼
run_grid.py  ── 生成 6 维笛卡尔积参数点
       │        断点续算(跳过已成功) / 种子复用(最近邻) /
       │        失败隔离 / 冷启动失败→种子步进回退
       │
       ├── worker 1 ── run_one.py ── 点 A
       ├── worker 2 ── run_one.py ── 点 B        每个 worker 独立工作目录
       ├── ...                                    互不干扰，并行(默认16核)
       └── worker N ── run_one.py ── 点 X
                          │
                          ▼
            冷启动链(lte→nc→nl) + synspec
                          │  冷启动发散且有干净邻居种子？
                          ▼  → 移失败结果到 <model>.coldfail/，
            种子步进链(seed_nc→nl) + synspec   改用 LTGRAY=F 热启动重试
                          │
                          ▼
                   results/<模型名>/
                     conv.json    ← 阶段信息(收敛/迭代/时间/seed_step_used)
                     *.spec/.cont ← 光谱
                     *.7          ← 各阶段大气

2. 每个网格点的计算阶段

每个网格点（一组 Teff/logg/logHe/logC/logN/logO 参数）经过 4 个阶段。

默认走"冷启动链"（DEFAULT_CHAIN，适用 20-40K 及大部分易收敛点）：

阶段 程序 做什么 NITER 典型耗时
1. LTE 灰大气 tlusty T T 模式，解析求灰色 T(τ) 结构 0 1-3 秒
2. nc（NLTE 连续谱） tlusty F F + ilvlin=0，收敛电离平衡（无线跃迁） 10 1-5 分钟
3. nl（NLTE 含线） tlusty F F + ilvlin=100，加全部谱线跃迁（要求收敛） 100 5-20 分钟
4. synspec synspec 用 nl 大气合成可观测光谱 — 3-10 秒

阶段间依赖：1→2→3→4 严格顺序。每阶段用上一阶段的 .7 大气作种子（fort.8）。

冷启动链在高温/He-poor/富金属区会发散，此时 run_grid 自动改走下文 §2.1 的种子步进链（seed_nc→nl，跳过 LTE grey 直接热启动）。

为什么是这 4 个阶段（原理）

Tlusty 的 NLTE 求解用迭代线性化（complete linearization）。线性化的收敛半径 有限——当初猜离真解太远时迭代发散。四个阶段逐步缩小初猜与真解的差距：

  • 阶段1（LTE 灰大气）：LTE + 灰色不透明度假设下解析求温度结构。提供物理 合理的起点，不需要种子（从零开始）。
  • 阶段2（nc 连续谱）：切换到 NLTE，但 ilvlin=0 不含束缚-束缚线跃迁。 线跃迁是统计平衡方程中最敏感的非线性项；先不加线，只收敛电离平衡（光致电离 +复合），得到稳定的 NLTE 布居数结构。跳过此步直接 grey→含线 NLTE 必发散。
  • 阶段3（nl 含线）：加入全部线跃迁（ilvlin=100）。从已收敛的 nc 种子起步， 线扰动小，快速收敛（典型 ~15 次迭代）。
  • 阶段4（synspec）：用 nl 阶段收敛的大气模型，计算指定波长范围的合成光谱。

2.1 种子步进链（seed_step）—— 高温区破局（NEW）

当冷启动链发散时（典型：60000-80000K + He-poor + 富金属），run_grid 自动 切换到种子步进链（seed_step.SEED_STEP_CHAIN），跳过 LTE grey 冷启动， 直接用邻居已收敛的 .7 热启动：

阶段 做什么 关键标志 NITER
seed_nc 从种子热启动 NLTE 连续谱 LTGRAY=F(读 fort.8) + ICHANG=0 80
nl 含谱线完整 NLTE（要求收敛） ilvlin=100 100

触发流程（run_grid.py 的 _worker）：

  1. 先跑冷启动链（DEFAULT_CHAIN）。
  2. 若 converged=false 且 seed_step_fallback=true 且找到了干净邻居种子：
    • 把失败结果整体移到 results/<model>.coldfail/（避免污染种子库）；
    • 用 SEED_STEP_CHAIN（seed=<邻居>.7）重算；
    • conv.json 里记 seed_step_used=true、coldfail_backup=<路径>。

原理（tlusty208.f 源码确认，详见 EXPERIENCE.md §5Y）：

  • LTGREY=T → CALL LTEGR 生成灰大气（忽略 fort.8，冷启动）；
  • LTGREY=F → CALL INPMOD 读 fort.8 作初猜（热启动）。
  • ICHANG=0：种子与目标模型原子完全一致（同为 H/He/CNO 设置），只改 Teff/logg/丰度，不需要重映射布居数。

实测突破（80K + He-poor + logCNO=-4，冷启动必败点）：种子步进 126s 收敛 （冷启动 970s 发散到 NaN）。详见 EXPERIENCE.md §5Y 验证表。

物理极限（如实标注）：80000K + He-poor + logCNO=-1 即使种子步进也发散 （CNO 高价离子主导不透明度，金属 ×10 跳跃线性化无法阻尼）。网格如实标 converged=false，不强制成功——这些极端参数组合观测上本就罕见。

ICRSW 是死代码：源码中 SUBROUTINE SWITCH 定义但从未被 CALL， 默认 CRSW≡1.0 无阻尼效果。不要依赖 ICRSW 稳定化（见 §7.8）。

3. 每阶段的配置信息与原理

3.1 .5 文件（每阶段一份，三阶段相同 NATOMS/ions，只改 3 处）

第1行: TEFF  GRAV              （三阶段相同：目标参数）
第2行: LTE  LTGRAY             （阶段1=T T，阶段2/3=F F；种子步进链全 F）
第3行: nst 文件名               （固定写 'nst'，内容每阶段由 write_nst 生成）
第4行: NFREAD                   （=2000 → 展开约 75443 个频率点；不要用 50）
第5行: NATOMS                   （=8: H,He,空×3,C,N,O）
第6+行: atoms (mode abn modpf)  （C/N/O 的 mode=2 显式NLTE, abn=10^logX）
ions段: iat iz nlevs ilast ilvlin nonstd typion filei
       （阶段1/2/seed_nc: ilvlin=0; 阶段3/nl: ilvlin=100 ← 关键区别）

为什么 NATOMS/ions 三阶段必须相同：每阶段的 .7 大气记录了每个能级的 布居数。种子与目标的能级结构必须一一对应，否则读取时索引错位 → NaN。

3.2 nst 文件（非标准参数，每阶段不同）

阶段1（LTE 灰大气）—— 保持干净，不加稳定化参数：

ND=50,VTB=2.,NITER=0
  • NITER=0：灰大气不迭代，只做一次形式解。

阶段2（nc）和阶段3（nl）—— 频率细化（用户验证配方，不设 CHMAX/ITEK）：

ND=50,NLAMBD=3,VTB=2.,ISPODF=1,DDNU=50.,CNU1=6.,NITER=<阶段>
IELCOR=-1
  • nc: NITER=10, nl: NITER=100
  • 不设 CHMAX（用默认 0.001，强迫 nc 真正收敛）
  • 不设 ITEK（用默认 4）

每个参数的作用与原理：

参数 nc值 nl值 作用 为什么这样设
ND 50 50 大气深度点数 sdB 标准配置
NLAMBD 3 3 lambda 迭代频率点数 频率网格细化（用户验证配方）
VTB 2. 2. 微湍流速度 km/s sdB 典型值
ISPODF 1 1 频率网格开关 启用细化频率网格
DDNU 50. 50. 频率间隔因子 频率网格细化参数
CNU1 6. 6. 频率网格起点 频率网格细化参数
NITER 10 100 最大迭代数 nc 给 10 次足够（实测最优）；nl 给 100 次
IELCOR -1 -1 电子密度修正 关闭

关键：不设 CHMAX（用默认 0.001）、不设 ITEK（用默认 4）、不设 IDLTE/ORELAX。 之前版本设了 CHMAX=0.1 导致 nc 没真正收敛，是大部分失败的根本原因。 详见 EXPERIENCE.md §4 的错误记录。

nc 的 NITER=10 是实测最优（GUIDE.md NITER 扫描结论）：

  • nc（纯连续谱）缺少谱线约束，外层温度永不真正收敛，只在外层漂移；
  • NITER=10 vs NITER=50 的最终光谱差异仅 6e-6（完全等价）；
  • nl（含谱线）会自修正到正确解，无论 nc 给什么初值；
  • NITER=10 总耗时 ~12 分钟（35000K CNO），NITER=50 浪费 2.2× 时间。

3.3 synspec 配置（fort.55.lin + 谱线表）

fort.55.lin 第6行: WLMIN  WLMAX  WLSTEP  ...  CUTOFF  ...
谱线表 fort.19:     data/gfVIS99.dat (含 C 1412 / N 2396 / O 1885 条线)
  • 当前用 3000-7000Å（光学波段，覆盖 C II 4267、C III 4647 等）。
  • 大气来自 nl 阶段的 .7（复制为 fort.8）。

4. CPU 与并行机制

4.1 每个网格点只用一个 CPU 核

是的。 tlusty.exe 和 synspec.exe 是 Fortran 编译的单线程程序，每个实例只用 1 个 CPU 核。网格点的并行不是靠程序内部的多线程，而是靠同时启动多个程序实例。

4.2 如何做到并行

run_grid.py 用 Python 的 multiprocessing.Pool（run_grid.py 的 _worker）：

with Pool(nworkers) as pool:
    for res in pool.imap_unordered(_worker, worker_args):
        ...
  • nworkers（config.yaml，当前=16）：同时运行的 worker 进程数。
  • 每个 worker 是一个独立的 Python 子进程，调用 run_one.py 跑一个网格点 （在独立的工作目录里，互不干扰）。
  • imap_unordered：哪个点先完成就先回收，立即分配下一个点（动态负载均衡）。
  • 16 核机器跑 16 个 worker = 16 个 tlusty 实例同时跑 = 满载利用。 （按机器核数调整；每个 tlusty 运行是单线程的，调大 nworkers 即可吃更多核。）

关键：每个 worker 用独立工作目录（results/<模型名>/），避免 fort.* 文件 冲突。这是并行安全的基础。

4.3 吞吐量估算

模型类型 单点耗时 16核并行吞吐 收敛性
20000-40000K（标准 sdB 区） ~12-25 分钟 ~48-80 点/小时 冷启动全区间可靠
60000K + He-rich/低金属 ~10-15 分钟 ~64-96 点/小时 冷启动或一步种子步进
80000K + He-rich/低金属 ~3-5 分钟（种子步进） 约同上 冷启动失败→种子步进成功
80000K + He-poor + logCNO=-1 — — 真实物理极限，发散（如实标注）

当前 config.yaml 共 432 点（4×2×2×3×3*3）；中等参数区冷启动为主， 高温区走种子步进回退，整体约需数小时到一天。

5. 如何统计每阶段信息

5.1 当前已记录的信息（conv.json）

每个网格点完成后，results/<模型名>/conv.json 记录：

{
  "name": "t40000_g6.0_he0_c-1_n-1_o-1",
  "params": {"teff":40000, "logg":6.0, "loghe":0, "logc":-1, "logn":-1, "logo":-1},
  "converged": true,
  "final_max_relc": 0.0069,
  "atmosphere_has_nan": false,
  "synspec_rc": 0,
  "elapsed_sec": 715.0,            ← 总耗时（所有阶段+synspec之和）
  "seed": null,                    ← 冷启动为 null；走种子步进时为邻居 .7 路径
  "seed_step_used": false,         ← true=种子步进链跑成功的（含 coldfail_backup 路径）
  "stages": [
    {
      "label": "lte",
      "converged": true,
      "final": {"itek":null, "rc":0, "max_relc":0.0,
                "note":"NITER=0 grey start (no iterations)"}
    },
    {
      "label": "nc",
      "converged": false,          ← nc 不要求收敛，作种子即可（NITER=10）
      "final": {"itek":null, "rc":0, "max_relc":0.957,
                "worst_depth":1, "last_iter":10, "n_depths":50}
    },
    {
      "label": "nl",
      "converged": true,
      "final": {"itek":null, "rc":0, "max_relc":0.0069,
                "worst_depth":1, "last_iter":17, "n_depths":50}
    }
  ]
}

走种子步进链时：seed 指向邻居 .7，seed_step_used=true， coldfail_backup 指向 <model>.coldfail/，stages 里没有 lte，而是 seed_nc（LTGRAY=F 热启动，NITER=80）→ nl。

每阶段记录：converged（是否收敛）、max_relc（最大相对变化）、 worst_depth（最差深度点）、last_iter（迭代次数）、n_depths（深度点数）、 elapsed_sec（本阶段耗时）。

5.2 每阶段时间记录（已实现）

run_one.py 现在在每个阶段的循环开始/结束处计时，conv.json 里每个 stage 有 elapsed_sec，synspec 也有单独的 synspec_sec：

"stages": [
  {"label":"lte", "elapsed_sec": 2.1,  "converged":true,  ...},
  {"label":"nc",  "elapsed_sec": 62.4, "converged":false, ...},
  {"label":"nl",  "elapsed_sec": 7.8,  "converged":true,  ...}
],
"synspec_sec": 3.1,
"elapsed_sec": 75.4

统计所有模型的阶段时间分布：

python3 -c "
import json,glob
for f in sorted(glob.glob('results/*/conv.json')):
    j=json.load(open(f))
    times = {s['label']:s.get('elapsed_sec',0) for s in j['stages']}
    print('%-30s lte=%5.0fs nc=%5.0fs nl=%5.0fs syn=%4.0fs total=%5.0fs' % (
        j['name'], times.get('lte',0), times.get('nc',0), times.get('nl',0),
        j.get('synspec_sec',0), j['elapsed_sec']))
"

5.3 统计整个网格的信息

run_grid.py 完成后写 results/grid_status.json：

{
  "total": 432,
  "elapsed_sec": 36000,
  "counts": {"converged": 400, "unfinished": 18, "error": 3, "skipped": 11},
  "seed_step_retries": 47,
  "models": [
    {"name":"t20000_...", "status":"converged", "max_relc":0.0065,
     "seed_step_used": false},
    {"name":"t80000_...", "status":"converged", "max_relc":0.00091,
     "seed_step_used": true},
    ...
  ]
}

汇总统计命令：

# 成功率 + 种子步进命中数
python3 -c "import json; j=json.load(open('results/grid_status.json')); print(j['counts'], 'seed_step_retries=', j['seed_step_retries'])"
# 所有收敛模型的 max_relc 分布（标注是否走了种子步进）
python3 -c "
import json,glob
for f in sorted(glob.glob('results/*/conv.json')):
    j=json.load(open(f))
    if j['converged']:
        tag='SEED' if j.get('seed_step_used') else 'cold'
        print(j['name'], tag, j['final_max_relc'], str(j['elapsed_sec'])+'s')
"
# 失败/未收敛的模型
python3 -c "
import json,glob
for f in sorted(glob.glob('results/*/conv.json')):
    j=json.load(open(f))
    if not j['converged']:
        print(j['name'], 'FAILED', j.get('note',''))
"

5.4 单个网格点的详细收敛诊断

# 看某阶段的迭代收敛趋势（fort.9）
python3 src/check_conv.py results/<模型>/<模型>.nl.9 --chmax 0.01

# 画光谱（标出 CNO 诊断线位置）
python3 src/plot_spec.py results/<模型>

6. 完整操作步骤

第1步：配置网格密度（config.yaml 的 grid 段）

grid:
  teff:  [20000, 30000, 40000, 60000]   # 各维采样点列表
  logg:  [5.0, 6.0]
  loghe: [-2, 0]
  logc:  [-4, -2, -1]                    # 亚太阳范围（已修正；旧版 -1..1 会发散）
  logn:  [-4, -2, -1]
  logo:  [-4, -2, -1]
  # 共 4*2*2*3*3*3 = 432 个点

第2步：设置环境变量

export TLUSTY=/home/dckj/program/tlusty/tl208-s54

第3步：预览（dry-run）

cd $TLUSTY/cno_grid
python3 src/run_grid.py config.yaml --dry-run
# 输出：grid: 432 points total, N already done, M to compute

第4步：启动批量计算（后台并行）

nohup python3 src/run_grid.py config.yaml > results/grid_run.log 2>&1 &
# 冷启动失败的点会自动尝试种子步进回退（seed_step_fallback: true）。
# 想关闭回退：在 config.yaml 设 seed_step_fallback: false。

第5步：监控

tail -f results/grid_run.log          # 实时进度
grep seed_step results/grid_run.log   # 看哪些点走了种子步进
cat results/grid_status.json          # 汇总（完成后才有）

第6步：断点续算（中断后恢复，自动跳过已成功的）

python3 src/run_grid.py config.yaml   # 重跑同一命令即可

第7步：检查结果 + 画图

# 成功率 + 种子步进命中数
python3 -c "import json;j=json.load(open('results/grid_status.json'));print(j['counts'],'seed_step=',j['seed_step_retries'])"
# 画某个模型光谱
python3 src/plot_spec.py results/<模型名>

单点调试（不走批量）

# 冷启动单点（适用 20-40K 大部分点）
python3 src/run_one.py --teff 40000 --logg 6.0 --loghe 0 --logc -1 --logn -1 --logo -1
# 种子步进单点（高温 He-poor 等冷启动失败点）
python3 src/seed_step.py --teff 80000 --logg 6.5 --loghe -4 \
    --logc -4 --logn -4 --logo -4 \
    --seed results/<seed-model>/<seed-model>.7

第8步：加密网格（Phase 2）

在 config.yaml 的各维列表里加更多点，重跑 run_grid.py（自动跳过已完成的， 只算新点）。高温区加密时建议先冷启动算 He-rich/低金属的"桥头堡"模型， 建立种子库再扩散到难收敛点（见 §7.1）。

7. 注意事项

  1. 种子步进回退（核心机制）：seed_step_fallback: true（默认开）时， 冷启动失败的点会自动改走种子步进链（seed_step.SEED_STEP_CHAIN）： 把失败结果备份到 <model>.coldfail/，用 find_seed 找最近邻已收敛的 .7 作种子，LTGRAY=F + ICHANG=0 热启动重跑。这是高温/He-poor/富金属区的 决定性破局手段（详见 §2.1）。 种子跨度限制：种子与目标的参数差不能太大（logg ≤0.5/步， Teff ≤5000K/步；logCNO 一步 ≤100×）。跨度过大（如 cno-4 直接跳 cno-1， 1000× 金属跳跃）即便种子步进也发散 → 这是真实物理极限，网格如实标注。 种子库建立策略：高温区建议先冷启动算 He-rich 或低金属的"桥头堡" 模型，建立种子库再扩散到难收敛点（如 cno-4 → cno-2 → cno-1 多步跳板）。

  2. 断点续算判定：conv.json 里 converged=true 的点会被跳过。假收敛 （atmosphere_has_nan=true）的点会被重算。

  3. 失败隔离：单点失败（发散/崩溃）不中断整个网格，记入 grid_status.json 的 error/unfinished 列表。高温难点的失败大多是真实物理极限（见 §2.1）， 不强制成功；如确需重试，可用 seed_step.py 手动从更近的种子起跳。

  4. 磁盘空间：每个模型约 50-100MB（含中间文件，走种子步进的点还会留 <model>.coldfail/ 备份）。432 点约需 20-40GB。如不够可定期清理中间文件 （保留 .spec/.cont/.7/conv.json）。

  5. 并行安全：每个 worker 用独立工作目录（results/<模型名>/），fort.* 文件 不冲突。可安全并行。

  6. nst 文件行长限制（已修复）：tlusty 的 nst 解析器有 ~72 字符的行宽限制。 如果参数太多写在一行（如加了 IDLTE/IACC 后 >72 字符），行尾的参数会被 静默截断。write_nst() 现在把参数分两行写（line1≤64c, line2 余下参数）。 这是个隐蔽 bug——截断后 tlusty 不报错而是用默认值，导致"看似成功实则参数 没生效"。

  7. fort.84 残留（已修复）：tlusty 运行时会在工作目录写 fort.84（nst 参数的 内部表示）。如果下一次 tlusty 运行（不同 NATOMS）读到旧的 fort.84，会报 "Bad integer for item 48" 崩溃。run_tlusty() 现在每次运行前删除 fort.84。

  8. 收敛可靠性（如实）：用正确配方（无 CHMAX/ITEK, NFREAD=2000, nc NITER=10） 后，20000-40000K 全区间冷启动可靠收敛；60000-80000K + He-rich/低金属用 种子步进可解决；80000K + He-poor + logCNO=-1 是真实物理极限 （CNO 高价离子主导不透明度，金属 ×10 跳跃线性化无法阻尼），网格如实标记 未收敛。之前版本的"8/8 边界全部成功"不准确（基于错误的 CHMAX=0.1）。 ICRSW（Hummer & Voels 切换）在 tlusty208 原版中是死代码 （SUBROUTINE SWITCH 从未被 CALL，CRSW≡1.0），不要依赖它稳定化； 即便源码修复启用后实测对难点也无帮助。详见 EXPERIENCE.md §5Y。