# 第一梯队补算脚本

配合 `修改报告_FEM分析审查.md` / https://fenxi.aiimage.icu 使用。全部脚本**只驱动你现有的代码，不修改任何上游源文件**——两个没有命令行开关的选项（均质核、稳定性模态）用运行时 patch `_options` 实现，这样上游代码树保持干净、可 diff。

## 先决条件：脚本必须跑在真实工程树上，不能跑在投稿包上

审计时发现投稿包 `05_Source_Data` 里的脚本**无法独立运行**——它们 import 的这些东西一个都不在包里：

- `src/nuclear_envelope_fem/`（`nonlinear_solver`、`chromatin_heterogeneity`、`tags`）
- `run_extreme_finite_deformation.py`（提供 `BASE_CONFIG`、`MESHES`、`_run_case`）
- `run_vaziri_finan_two_group.py`、`validate_literature_reproduction.py`
- gmsh 网格 `ellipsoid_core_shell.msh` 及**生成它的脚本**
- 染色质体数据 `reynolds_2021_figure2_chromatin_50cube.npz`

这与 §10 "custom code ... included in this local submission package" 的声明冲突，**投稿前必须补齐**，否则任何人都无法复现。下面所有 `--project-root` 都指向你本地的真实工程树（含 `src/` 和 `results/` 的那个），不是解压出来的投稿包。

先用 `--dry-run` 确认计划，再去掉它真跑。

## 运行顺序

### 1. 机制分解（最高优先级，6 次求解）

```bash
python3 02_run_mechanism_decomposition.py \
    --project-root /path/to/project --out-root /path/to/out/decomposition --dry-run
```

三个配置：`probe_only`（预张力=0）、`pretension_only`（探针=0）、`full`（复现现状）。产出 `mechanism_decomposition.json`，直接给出总差值中预张力松弛、探针响应、耦合项各占多少。

为什么必须做：预张力 0.05 mN/m 占 Control 端点 0.066 的 **76%**；代码确认它是残差化初应力（`PK1 += (F−I)·S₀`，u=0 时贡献为零），且被推前进入 proxy 的 Cauchy 应力，所以头条数字**确实包含全部预张力贡献**。而 UQ 的 pretension 范围 0.025–0.1 从不含 0，这个分解从未做过。

### 2. 均质核对照（第二优先级，先算等效模量再跑 4 次求解）

```bash
python3 04_equivalent_core_modulus.py \
    --project-root /path/to/project --json-out equivalent_modulus.json

python3 03_run_homogeneous_control.py \
    --project-root /path/to/project --out-root /path/to/out/homogeneous \
    --core-young-pa <REUSS> <VOIGT> --label reuss voigt --dry-run
```

`04` 调用**上游同一个 loader、同一套网格**算体积加权的 Voigt/Reuss 等效剪切模量（并按 K/G=30 换成 E），保证对照公平。不要用 Table S1 里的 800 Pa——那只是上游的均质回退值（Table S1 里 Core E=800/ν=0.45 的来源之谜就是它），与异质材料无关。Voigt/Reuss 是任何合理均匀化的上下界，两端都跑就把答案夹住了。

Main Fig 3c 里 chromatin scale 的 Spearman ρ 只有 **−0.08**，已经预告了这个对照大概率显示"全局端点对异质性不敏感"。

### 3. production 换 P2（2 次求解，`--include-main-p2` 上游已支持）

```bash
./01_run_main_p2.sh /path/to/project /path/to/out/p2
```

P1 族还在漂（−37.9 / −40.9 / −42.3%），P2 族已收敛（−44.64 / −44.51%，差 0.13 pp）。用 `projected_minres`——之前 augmented direct 在大 P2 上显存失败的记录还在。

### 4. 稳定性检查（第二梯队但近乎免费）

```bash
python3 05_run_stability_check.py \
    --project-root /path/to/project --out-root /path/to/out/stability --dry-run
```

`stability_modes` 在求解器里默认 6、`_solve_tangent_stability_modes` 已完整实现，但运行脚本显式设成了 **0**。这个脚本把它打开。必要性：高渗组压缩面积达 62.6%（Fig S4d），且残差在载荷因子 ≈0.6 处飙升两个量级（Fig S1）——都是接近临界点的样子。

### 5. 汇总出表

```bash
python3 06_collect_results.py \
    --runs /path/to/out/p2/p2_production /path/to/out/decomposition/full \
    --decomposition /path/to/out/decomposition/mechanism_decomposition.json \
    --homogeneous /path/to/out/homogeneous/homogeneous_control.json \
    --markdown-out revised_tables.md
```

产出可直接贴进手稿的修订版 Table S2（**新增组间变化列**，P1/P2 分族），以及分解表和对照表。

收敛判定是内置的、不迁就：单调才算 Richardson/GCI，非单调直接拒绝出数并说明理由。已用已发表的 Table S2 数据回归验证，输出与手工复算完全一致：

```
Control/P1        观测阶 2.89, Richardson 极限 0.066378, GCI 0.59%
Hyperosmotic/P1   NOT convergent（增量变号），拒绝给 GCI
跨离散化散布       −37.91% ~ −44.64%（6.7 pp）
```

## 这批脚本改不了的四件事

以下四条无法靠重跑解决，只能改代码或改写正文，见修改报告：

1. **材料场随网格变**（逐四面体赋组）——P1/P2 外推到不同极限的根源，需要网格无关的 L² 投影；
2. **壳是 ν=0.35 的可压缩实体薄层**，最小 J_e≈0.69 意味着 31% 体积压缩，真实 lamina 靠减薄不靠丢体积，应改膜/近不可压公式；
3. **proxy 公式正文缺失**——代码里是：外表面切向主 Cauchy 应力（含推前初应力）× 当前厚度（0.12 µm × 由 J/面积伸缩运动学估计的法向伸缩），当前面积加权。照抄进方法学即可；
4. **探针定义**——代码是均匀对称名义应力张量 σ = τ(eₓ⊗e_y + e_y⊗eₓ)、t₀ = σ·N 施加于参考构型（dead load，两组外载严格相同）。实现没问题，是正文没写。我此前据文档判断"名实不符"的指控据此**撤回**。