convergence-study

收敛性研究

目标

通过脚本驱动的收敛性分析，验证数值解在网格或时间步长细化时是否以预期速率收敛。

要求

• Python 3.10+
• 无第三方软件包——脚本仅使用 Python 标准库（math）。

需收集的输入

输入	说明	示例
网格间距	网格尺寸序列（从粗到细）	0.4,0.2,0.1,0.05
时间步长	dt 值序列	0.04,0.02,0.01
解值	各细化级别的 QoI	1.16,1.04,1.01,1.0025
预期阶数	数值格式的形式阶数	2.0
安全系数	GCI 安全系数（默认为 1.25）	1.25

脚本输出（JSON 字段）

脚本	关键输出
scripts/h_refinement.py	results.observed_orders、results.mean_order、results.richardson_extrapolated_value、results.convergence_assessment
scripts/dt_refinement.py	与 h_refinement 相同，但用于时间收敛性分析
scripts/richardson_extrapolation.py	results.extrapolated_value、results.error_estimate、results.observed_order
scripts/gci_calculator.py	results.observed_order、results.gci_fine、results.gci_coarse、results.asymptotic_ratio、results.in_asymptotic_range、results.extrapolated_value、results.notes

工作流

1. 开展网格/时间步长细化研究，至少使用 3 个级别
2. 使用 h_refinement.py 或 dt_refinement.py 计算观测收敛阶数
3. 将观测阶数与格式的预期阶数进行比较
4. 通过 Richardson 外推估算离散误差
5. 使用 gci_calculator.py 报告 GCI，用于形式解验证
6. 记录收敛结果及任何异常

决策指南

Do you have 3+ refinement levels?
+-- YES --> Run h_refinement.py or dt_refinement.py
|           +-- Observed order matches expected? --> Solution verified
|           +-- Order too low? --> Check: pre-asymptotic, coding error, insufficient resolution
|           +-- Order too high? --> Check: superconvergence or cancellation effects
+-- NO (only 2 levels) --> Use richardson_extrapolation.py with assumed order
                           (less reliable without order verification)

CLI 示例

# Spatial convergence with 4 grid levels
python3 scripts/h_refinement.py --spacings 0.4,0.2,0.1,0.05 --values 1.16,1.04,1.01,1.0025 --expected-order 2.0 --json

# Temporal convergence with 3 timestep levels
python3 scripts/dt_refinement.py --timesteps 0.04,0.02,0.01 --values 2.12,2.03,2.0075 --expected-order 2.0 --json

# Richardson extrapolation with assumed 2nd-order
python3 scripts/richardson_extrapolation.py --spacings 0.02,0.01 --values 1.0032,1.0008 --order 2.0 --json

# GCI for 3-mesh verification
python3 scripts/gci_calculator.py --spacings 0.04,0.02,0.01 --values 1.0128,1.0032,1.0008 --json

错误处理

错误	原因	解决方法
spacings and values must have the same length	输入数组长度不匹配	提供长度相同的列表
At least 2 refinement levels required	数据点过少	增加细化级别
Exactly 3 refinement levels required	GCI 需要 3 个级别	提供细、 中、粗网格
Oscillatory convergence detected	收敛非单调	检查网格质量或数值格式

结果解读指南

情况	含义	操作
观测阶数与预期阶数一致	处于渐近区间的最有力证据	报告 GCI 并进行外推
观测阶数 < 预期阶数	尚未进入渐近区间或存在代码错误	进一步细化或调试
观测阶数为负	解正在发散	检查实现
GCI 渐近比接近 1.0	见下方注意事项	通过阶数比较进行确认
GCI 渐近比远离 1.0	不在渐近区间内	进一步细化

渐近比注意事项（等比例细化）。 当细化比相等（r21 == r32，这是常见情况）时，渐近比 AR = GCI_coarse / (r^p * GCI_fine) 在代数上可化简为 f1/f2。此时它仅衡量两个最细 QoI 值之间的相对差距，并不反映数据是否符合假设的幂律。因此，即使观测阶数与预期相差甚远，AR 接近 1.0 也可能造成虚假的安心。gci_calculator.py 的 JSON 输出会在 notes 条目中对此进行标记。要真正判定是否处于渐近区间：

1. 将观测阶数 p 与数值格式的理论/预期阶数进行比较——两者一致才是处于渐近区间的有意义证据。
2. 如需更有力的验证，请使用 4 个或更多按规律细化的网格，并检查连续网格三元组上的观测阶数是否一致（h_refinement.py 会报告每个网格三元组的阶数及 mean_order）。

验证清单

• 使用了 >= 3 个系统性细化的网格/时间步，使 h_refinement.py / dt_refinement.py 能够报告 results.mean_order；记录了每组三元组的 results.observed_orders（单对 Richardson 运行无法验证阶数）。
• 记录了 results.mean_order，并确认 results.convergence_assessment 显示 PASS（观测阶数与方案的预期阶数相差在 10% 以内）；FAIL 或 unknown 表示结果尚未验证。
• 确认 results.in_asymptotic_range 为 true，并且在引用任何 GCI 或外推值之前，notes 中没有报告预渐近（阶数变化 >50%）、负/非正阶数或误差差值为零的条目。
• 检查了 gci_calculator.py 的细化比 r21/r32 是否 >= 1.3（舍入误差噪声下限），并确认未触发 Oscillatory convergence detected 错误。
• 记录了 results.gci_fine（以及所用的安全系数：经过验证的阶数且网格数 >= 3 时为 1.25；网格数为 2 且假设阶数时为 3.0），作为报告的离散化不确定性；同时将 results.extrapolated_value 作为最佳估计值。
• 对于恒定细化比，未将接近 1.0 的 asymptotic_ratio 视为渐近范围的证明（它会退化为 f1/f2）；通过比较观测阶数与预期阶数确认渐近范围，并留意 gci_calculator.py 的 notes 中的注意事项。

常见误区与错误理由

看似诱人的捷径	错误原因 / 应采取的做法
“两个网格的结果很接近，所以已经收敛”	两个层级无法估算观测阶数。使用 >= 3 个层级运行 h_refinement.py/dt_refinement.py；双网格 Richardson 运行使用的是假设阶数，且需要安全系数 3.0，而非 1.25。
“渐近比约为 1.0，所以我们处于渐近范围内”	在细化比恒定时，AR = GCI_coarse/(r^p*GCI_fine) 会简化为 f1/f2，只衡量两个最细网格 QoI 值之间的差距。应通过比较观测阶数与预期阶数（并检查 4 个以上网格的一致性）来验证渐近范围。
“GCI_fine 很小，所以解与网格无关”	GCI 接近零也可能意味着 QoI 对细化不敏感，或者差值处于舍入误差噪声范围内（比值 < 1.3）。应先确认 in_asymptotic_range 为 true，且细化比 >= 1.3。
“观测阶数高于预期，这样更好”	阶数远高于形式阶数通常意味着超收敛或误差抵消，而非精度更高。应将其视为警示（偏差 >10% 时，convergence_assessment 为 FAIL），并通过增加网格层级进行验证。
“脚本输出了外推值，所以就用它”	当观测阶数为非正值时，解正在发散，h_refinement.py 会返回 richardson_extrapolated_value = null，并附带发散说明。此时不要引用外推值或 GCI。
“求解器是隐式/稳定的，所以研究中时间步可以任意取值”	时间上的稳定性不等于时间上的精度。dt_refinement.py 仍需要 >= 3 个系统性缩小的时间步，才能恢复该方案的阶数；时间步序列过粗会一直处于预渐近状态。

安全性

输入验证

• 所有数值参数（spacings、timesteps、values、expected-order、order）都会被验证为有限正数
• 逗号分隔的值列表会进行长度匹配（spacings 和 values 的长度必须相同），且最多包含 10,000 个条目
• GCI 计算器强制要求恰好 3 个细化层级；Richardson 外推至少需要 2 个
• 安全系数会被验证为不小于 1.0 的有限数值（Roache 使用集合 {1.25, 3.0} 中的 Fs）

文件访问

• 脚本不读取外部文件；所有输入均通过 CLI 参数提供
• 脚本只向 stdout 写入内容（JSON 输出）；除非 Agent 显式使用 Write 工具，否则不会创建文件

工具限制

• Bash：用于使用明确的参数列表执行四个 Python 分析脚本（h_refinement.py、dt_refinement.py、richardson_extrapolation.py、gci_calculator.py）
• Read：用于检查脚本源代码和参考文档

安全措施

• 不使用 eval()、exec() 或动态代码生成
• 所有子进程调用均使用明确的参数列表（不使用 shell=True）
• 脚本仅使用 Python 标准库（math 模块）；不加载 pickle 或反序列化不可信数据
• 最小化工具范围（仅限 Bash 和 Read），以限制 Agent 修改文件系统的能力

参考资料

• references/convergence_theory.md - 正式的收敛阶、双对数分析、渐近范围
• references/gci_guidelines.md - Roache 的 GCI 方法、ASME V&V 20、安全系数