0 Fork 0

Files

T

wyj 265d7b95d5 Fix: reject invalid inverse lens-map weights

Check spherical area weight sums against a fixed 1e-8 tolerance and normalize accepted weights before interpolation. Skip invalid images in both catalog paths.

Add a thin-triangle regression for both parities, document the measured tolerance margin, and refresh HDR fixtures with strict comparisons restored.

2026-09-06 17:47:09 -04:00

32 KiB

Raw Blame History

Numerical-Relativity Spacetime Movie Renderer

1. Repo 目标

本项目的目标是开发一个离线、尽可能物理正确的数值相对论时空视频渲染器，最终用于渲染双黑洞（BBH）并合等动态数值时空的 4K 视频。

核心目标不是实时渲染，也不是构造视觉上“像双黑洞”的近似度规，而是：

直接消费数值相对论演化输出的真实 4D 时空；
对每一帧的相机光线做完整的 time-dependent backward null ray tracing；
支持真实巡天恒星数据作为无穷远天球背景；
正确处理恒星点源的多像、放大、频移和最终 PSF；
允许相机沿一般 4D worldline 运动，并使用预先生成的 tetrad/标架轨迹；
让平直时空、解析时空、数值时空在同一渲染框架中作为可替换 backend；
最终能够“看到”每次 NR 代码实际跑出来的时空，而不是只看 waveform 或标量诊断。

第一阶段不考虑物质辐射、吸积盘、流体、等离子体等局域发射源。每条 ray 的终点暂时只有两类：

被黑洞捕获；
到达无穷远天球。

2. 明确的非目标

当前 repo 不以以下内容为主要目标：

实时渲染；
ShaderToy / 游戏式视觉近似；
frozen-snapshot 作为最终物理方案；
直接渲染吸积盘、喷流或其他 matter emission；
GPU 加速；
通过 OO/class hierarchy 构建复杂抽象层；
将恒星背景预烘焙为 RGB 天球纹理；
依赖单个固定坐标位置的“相机”。

可以保留 frozen-snapshot 模式作为 nmesh 内部的廉价在线 diagnostic，但它与本 repo 的最终 4D renderer 是不同层次的工具。

3. 语言与总体实现风格

首选 C。

C++ 并不是当前项目的必要条件。需要的抽象主要可以通过：

struct
函数指针表
显式 ownership
大块连续数组
thread-local workspace

完成。

性能敏感部分尽量采用：

bulk array processing；
SoA（structure of arrays）；
避免 per-ray malloc；
避免大量虚调用/函数指针 dispatch；
避免细粒度锁；
避免每个 triangle/ray 建 task；
尽量顺序读取 metric 数据；
OpenMP 做 CPU 并行。

核心实现原则：

bulk arrays + time-slab streaming + coarse-grained parallelism

4. 总体物理与数据流

整个程序不是“逐帧、逐像素独立 trace”，而是：

预先生成相机 worldline 与 tetrad 轨迹；
建立所有视频帧的初始 image-plane adaptive mesh；
从所有帧收集当前 refinement level 需要的新 ray samples；
将这些 rays 组成一个全局 RayPool；
从视频结束时刻向过去，按 time slab 顺序加载数值时空；
在每个 slab 内，把所有 active rays 一起推进到 slab 左边界；
ray 若到达无穷远或进入黑洞，则立即终止；
一轮 ray tracing 完成后，把 endpoint 数据回填到各帧 image mesh；
根据局部 lens mapping 误差判断哪些 image-plane triangles 需要进一步细分；
生成下一批新增 rays；
重复若干 refinement passes；
最终用局部 inverse lens maps 查询恒星 catalog；
计算每个恒星像的位置、放大率和频移；
将对应黑体 PSF 加到 HDR framebuffer；
tone mapping / 编码输出视频。

总结构：

Observer Track
z^μ(t), e_(a)^μ(t)
        │
        ▼
Movie / Frame Set
        │
        ▼
Adaptive image-plane meshes
        │
        ▼
Refinement analyzer
        │
        ▼
new ray samples
        │
        ▼
RayPool
        │
        ▼
Spacetime time-slab stream
(newest → oldest)
        │
        ▼
ray endpoint:
n∞, frequency shift, captured/escaped
        │
        ├──────────► next refinement pass
        │
        ▼
Final local lens maps
        │
        ▼
Star catalog
(n∞, T, amplitude)
        │
        ▼
inverse mapping + μ + g + PSF
        │
        ▼
HDR video frames

5. 为什么必须按 time slab 组织，而不是按帧组织

最终 BBH 4D 数据很可能超过单机内存。

因此不能对每帧独立：

for frame:
    load all required 4D metric data
    trace all rays

否则同一段时空会被反复从硬盘读取。

更合理的是：

for time_slab = latest -> earliest:
    load metric data for this slab

    activate all rays whose camera time enters this slab

    evolve every active ray backward through this slab

    remove terminated rays

    discard slab

这样：

metric I/O 主要取决于 NR 数据量；
不随 4K 分辨率、帧数、supersampling 线性重复；
数据访问可以是大块顺序读；
适合本地 NVMe；
适合未来超过 128 GB 的数据集。

6. Ray 生命周期与跨帧整合

不应一开始生成整个视频所有 rays 并常驻内存。

随着倒序时间扫描，在经过某个 frame 的 t_camera 时才激活该帧当前 refinement pass 的新 rays。

因此 active-ray pool 大小由以下因素决定：

单位相机时间产生的 ray 数
×
ray 在 NR domain 内的典型 coordinate-time 停留长度

而不是整个视频总 ray 数。

对于长时间停留在 strong-field region 的 rays，可以一直跨多个 slabs 保留。

7. Adaptive image-plane mesh

最终点源渲染不适合简单的“每 pixel backward lookup sky color”。

恒星在无穷远天球上是点源：

[ I(\hat n) \sim \sum_s F_s \delta(\hat n-\hat n_s) ]

因此需要先通过 backward ray tracing 得到局部 lens mapping：

[ F:\ (x,y)\text{image} \to \hat n\infty ]

然后在每个局部可逆 patch 上建立 inverse：

[ F^{-1}:\ \hat n_\infty \to (x,y)_\text{image} ]

基本单元：triangle

每个 image-plane triangle 的三个顶点都保存：

image-plane 坐标 (x,y)；
ray 是否 escaped/captured；
若 escaped：无穷远方向 n_inf；
frequency shift / redshift accumulator。

示意：

typedef struct {
    double x, y;

    double n_inf[3];
    double log_g;

    uint8_t ray_status;
} LensVertex;

typedef struct {
    uint32_t v[3];

    uint8_t level;
    uint8_t flags;
} LensTriangle;

8. 局部 adaptive refinement

局部 inverse map 的球面面积权重必须满足凸组合约束。当前实现以无符号子面积除以总面积得到权重，在插值位置和频移前检查权重和有限、为正，且 abs(sum(weights) - 1) <= 1e-8；不满足时只跳过当前恒星在当前 triangle 中的像。通过检查后，将权重除以其总和，消除小的归一化误差。整 tile 包含的 catalog 快速路径也必须执行这项检查。

这是反演有效性保护，不是 refinement 阈值：细长 source triangle 的带容差边测试可能误接纳外部点，无符号子面积之和便大于总面积，直接插值会把像放到 image triangle 外。不能先归一化明显无效的权重来掩盖误接纳。固定的无量纲容差 1e-8 依据 Schwarzschild 3840×2160、45° FOV、R=100、 coarse cell=32、level=4、J=0.2 网格上的诊断选定：正对中心及偏离 0.00161°、方位为 1°/44°/23° 的单星测试中，严格包含样本的 double 权重和误差不超过约 2.5e-13，正对中心的越界样本最小误差约 5.9e-4。这不保证临界曲线附近的离散映射已经收敛，也不保证数值亮环无缺口。

每个 triangle 只关心自身局部映射。

不需要先把整个像平面映射到天球、再全局分类“第几阶像”。

对一个 triangle：

已有三个 corner rays；
计算或已有 edge midpoint / center rays；
比较真实 mapping 与低阶插值；
检查 orientation consistency；
检查 Jacobian 是否接近奇异；
若误差过大则 subdivide。

这意味着：

Schwarzschild 下无需显式按 winding number 分类；
BBH 下也不需要定义全局“第几阶像”；
每个局部 patch 只需保持局部单值可逆；
多像会自然表现为同一个 source direction 被多个 image patches 覆盖。

9. 为什么 refinement 需要多 pass

adaptive refinement 与 out-of-core time streaming 存在天然冲突。

只有当一条 ray 已经完整追到无穷远后，才知道某个 triangle 是否需要细分。

但此时早期 metric slab 已经释放。

因此不能在单次 time sweep 中随时创建新 ray 并从相机时刻重新追。

采用：

Pass 0:
    trace coarse samples through whole spacetime
    analyze maps

Pass 1:
    trace only newly requested midpoint/center rays
    analyze maps

Pass 2:
    trace only newly requested rays
    ...

until converged

每个 pass：

metric data 仍然只顺序扫描；
只追踪新增样本；
refinement 通常若干轮即可完成。

10. 相机不是固定坐标点

相机不应由：

double x[3];

定义。

否则结果强依赖某个 NR 坐标系中的“固定点”。

应由独立 trajectory generator 预先生成：

[ z^\mu(\tau), \qquad e_{(a)}{}^\mu(\tau) ]

其中：

(z^\mu)：observer worldline；
(e_{(0)}{}^\mu=u^\mu)：observer 四速度；
(e_{(i)}{}^\mu)：相机空间 tetrad。

trajectory generator 可以采用不同物理规则：

timelike geodesic；
Fermi-Walker transport；
parallel transport；
人为给定加速度；
人为叠加 camera pointing/rotation。

电影 renderer 本身只读取轨迹；单张入口可以用下述参数直接构造一个事件处的完整 tetrad。

示意：

typedef struct {
    double t;
    double tau;

    double x[3];
    double tetrad[4][4];
} ObserverSample;

typedef struct {
    ObserverSample *samples;
    size_t n;
} ObserverTrack;

每帧根据 t_camera 插值得到完整 observer state。

10.1 单张的瞬时相机

单张与电影共享 ObserverState 和 ray 初始化。单张不是只指定三维位置：由事件、坐标速度、指向与 roll 生成完整四速度和 tetrad。目前解析单张事件取 $t=0$；此构造不要求静态时空或静止观测者。

--observer-position X Y Z 与 --look-ra-deg / --look-dec-deg 独立。只有位置时取朝原点的坐标方向；只有指向时令 $\mathbf x=-R\mathbf d$。 --observer-radius R 仅用于缺省位置的推导（默认 30），与显式位置互斥。几何参数全部省略时保留 Minkowski 原点、Schwarzschild 半径 30 的默认相机。原点无法推导指向，要求显式 look；由位置推导方向且在极轴上时固定 RA=0。 RA 或 Dec 只给一个时，另一个取默认值（90°、−90°）。

--observer-velocity VX VY VZ 表示 $v^i=dx^i/dt$，默认零。由事件处的 \alpha,\beta^i,\gamma_{ij} 计算


Q=-\alpha^2+\gamma_{ij}(v^i+\beta^i)(v^j+\beta^j),\qquad
u^\mu=\frac{(1,v^i)}{\sqrt{-Q}}.

要求 $Q<0$，未来指向选择 $u^t>0$。不能用欧氏 |\mathbf v|<1 或 r>2M 代替类时性检查；删除旧的 --observer-inward-speed 参数。构造失败直接报错，不修改速度或退回静止相机。

采用标准 ICRS 坐标种子


\mathbf d=(\cos\delta\cos\alpha_{\rm RA},\cos\delta\sin\alpha_{\rm RA},\sin\delta),

将 D^\mu=(0,\mathbf d) 投影为 $F^\mu=D^\mu+u^\mu(u_\nu D^\nu)$，归一化得到 forward。随后把天球北向和西向的坐标种子投影、正交化为 up、right，保留标准画面定向。因此 look 表示坐标方向在相机静止空间内的投影，不保证特定远方物体位于画面中心。 --camera-roll-deg 默认零，正角按 $e'_2=\cos\rho,e_2+\sin\rho,e_3$、 e'_3=-\sin\rho\,e_2+\cos\rho\,e_3 定义。

observer 构造只接收当地 metric 和已补全的参数，不加载 slab、不分类 ray。调用方在昂贵的 catalog/PSF 初始化前验证相机及 backend 数据域。现有 Schwarzschild cutoff 为 $r=1.5M$；相机必须在 cutoff 外，但允许在视界内。此功能不改变捕获 cutoff 或向过去追踪的高红移终止条件。单张相机参数与轨迹输入、lens-map 导入互斥；导入仍跳过 metric 与 observer 初始化。

验证包括 tetrad 正交归一和 null 初始化、平直时空平移不变性与解析光行差/多普勒、 KS 视界处和视界内的向过去径向逃逸及频移，以及相同 observer 的单张/电影单帧一致性。

11. Ray 初始化

给定相机 film coordinate (x,y)：

camera projection 得到 observer local frame 中的 photon direction；
与 observer tetrad 相乘得到 spacetime 中的 (k^\mu)；
转换到 geodesic integrator 所使用的 3+1 ray variables。

例如：

[ k^\mu = e_{(a)}{}^\mu k^{(a)} ]

film model 只负责：

[ (x,y) \to \hat n^{(i)}_\text{camera} ]

它与 gauge、metric backend、observer worldline 都解耦。

12. Ray 状态

采用 coordinate time (t) 作为 ODE 自变量，以便天然和 time slabs 对齐。

单 ray 概念结构：

typedef struct {
    double t;

    double x[3];
    double Pi[3];

    double log_E;
    double h;

    uint32_t sample_id;
    uint32_t frame_id;

    uint8_t status;
} Ray;

真正批量计算时建议 SoA：

typedef struct {
    double *restrict t;

    double *restrict x0;
    double *restrict x1;
    double *restrict x2;

    double *restrict p0;
    double *restrict p1;
    double *restrict p2;

    double *restrict log_E;
    double *restrict h;

    uint32_t *sample_id;
    uint32_t *frame_id;

    uint8_t *status;

    size_t n;
    size_t capacity;
} RayPool;

目的：

cache locality；
SIMD；
OpenMP；
stream compaction；
避免 per-ray allocation。

13. Spacetime backend

所有上层渲染逻辑都不应该知道 metric 来自哪里。

需要支持至少：

Minkowski；
analytic Schwarzschild；
未来可加 analytic Kerr；
numerical nmesh 4D spacetime。

核心不是“单点 metric 对象”，而是当前已加载 time slab 的 evaluator。

概念接口：

typedef struct SpacetimeSource SpacetimeSource;
typedef struct MetricSlab MetricSlab;

typedef struct {
    int (*load_slab)(
        SpacetimeSource *src,
        double t_hi,
        size_t memory_budget,
        MetricSlab **out);

    void (*free_slab)(MetricSlab *slab);

    int (*eval)(
        const MetricSlab *slab,
        double t,
        const double x[3],
        void *workspace,
        MetricData *out);

    int (*classify)(
        const MetricSlab *slab,
        double t,
        const double x[3]);
} SpacetimeOps;

解析 backend 的 load_slab() 可以近似 no-op。

nmesh backend 则真正加载一段时间的 3D outputs。

14. 统一的 3+1 metric 数据

初步采用：

[ \alpha,\quad \beta^i,\quad \gamma_{ij},\quad K_{ij} ]

作为主要 3+1 变量。

空间导数：

[ \partial_i \alpha,\quad \partial_i\beta^j,\quad \partial_i\gamma_{jk} ]

由 backend evaluator 提供。

示意：

typedef struct {
    double alpha;

    double beta[3];
    double gamma[6];
    double K[6];

    double d_alpha[3];
    double d_beta[3][3];
    double d_gamma[3][6];
} MetricData;

15. 时间导数

时间导数不能简单假设全部从文件直接输出。

应按量的性质选择来源。

例如：

[ \partial_t \gamma_{ij}

-2\alpha K_{ij} + \mathcal L_\beta \gamma_{ij} ]

可以直接由 3+1 内部量构造。

对于 lapse/shift：

若 gauge evolution equations 已知，可从 gauge RHS 构造；
或从 temporal interpolant 解析求导。

原则：

若已经建立 (u(t,\mathbf x)) 的时间插值，就从同一个 interpolation polynomial 解析求 (\partial_tu)，而不是另做 finite difference。

且不要每个 geodesic RHS 都计算用不到的时间导数。

16. nmesh 数值时空 backend

nmesh 当前 3D 输出格式已经非常适合 renderer：

对每个需要输出的 timestep：

每个 AMR leaf node 一个文件；
文件名按 node 的树状 refinement 命名；
从文件名可恢复 refinement tree；
文件内保存每个 nodal point 的 global coordinates；
保存所选物理量。

nmesh：

DG；
node 内用 local coordinates 上 Legendre nodes 对应的 Lagrange basis；
production 中每边一般最多约 12 个 nodal points；
AMR 为 h-refinement；
典型网格为中心六面体 + 外层 cubic-sphere shells。

因此读入一个 node 文件后，可以直接构造 node 内：

[ u(\xi,\eta,\zeta)

\sum_{ijk} u_{ijk} L_i(\xi)L_j(\eta)L_k(\zeta) ]

空间导数直接对 Lagrange polynomial 求导。

不需要：

resample 到 Cartesian uniform grid；
额外输出 derivative fields。

17. nmesh time interpolation

h-AMR 可能使相邻 timestep 的树结构不同。

但不要求时间方向上 node 一一对应。

对任意固定 global point (\mathbf x)：

[ u_n(\mathbf x)

I_{\mathcal M_n}[u](\mathbf x) ]

[ u_{n+1}(\mathbf x)

I_{\mathcal M_{n+1}}[u](\mathbf x) ]

先分别在各自空间 mesh 上求值，再做：

[ u(t,\mathbf x)

T[u_n(\mathbf x),u_{n+1}(\mathbf x),\ldots] ]

因此：

每个 time slice 可独立恢复空间函数；
AMR topology 不需要跨时间匹配；
time slab 只负责保存足够多邻近 slices 供时间插值。

slab 边界需要少量 overlapping temporal ghost slices。

18. 黑洞捕获判据

目标使用 moving-puncture BBH，而不是 excision。

production renderer 不希望依赖每次 NR run 都开启昂贵的 AH finder。

计划：

用低分辨率 single-BH / BBH calibration run 开 AH finder；
测量 horizon 相对于 puncture 的最小 coordinate radius；
选择明显保守、始终位于 AH 内部的 puncture-centered cutoff；
正式 renderer 只根据 puncture trajectory 做判断。

形式：

[ |\mathbf x-\mathbf x_p(t)|<r_\text{cut} \Rightarrow \text{captured} ]

未来 BBH merger 后若仅用两个 puncture-centered 小球导致大量 doomed rays 继续积分，可以再加入 common-horizon-derived termination 优化。

19. 恒星 catalog 的内部表示

不要预存 RGB。

恒星应保留：

无穷远天球方向；
温度 (T)；
辐射 normalization / amplitude。

概念结构：

typedef struct {
    float n[3];

    float temperature;
    float amplitude;
} Star;

原因：

ray tracing 会给每个 image 一个频率比：

[ g= \frac{\nu_\text{obs}} {\nu_\text{emit}} ]

对于黑体谱：

[ T_\text{obs}=gT_\text{emit} ]

因此温度是更自然的源数据。

最终根据：

intrinsic temperature；
frequency shift (g)；
lensing magnification (\mu)；
source amplitude；

得到 observed blackbody spectrum，再转换为 linear RGB。

20. 点源渲染

不将恒星 catalog 插值为天球亮度纹理。

最终 frame mesh 已知局部：

[ (x,y) \to \hat n_\infty ]

对每个局部可逆 triangle，反过来得到：

[ \hat n_\infty \to (x,y) ]

并查询 source triangle 内的 catalog stars。

固定 2MASS tile 的基准结论（2026-08-27）

在 5d80bb1 上，以 54,160-source 的 2MASS 1 deg tile、完全位于该 tile 内的 0.6 deg Minkowski 视场、1920 x 1920 输出和 6 个 catalog render worker 测得：每个 source triangle 扫描完整 tile 的有效墙钟成本约为 105.9 us，并且 wall time 对 triangle 数的线性拟合为 11.855 s + 105.88 us * T（R^2 = 0.99955）。完整记录、命令和原始时间见 benchmarks/2mass_fixed_tile_triangle_scan_2026-08-27.md。

这说明局部稠密 tile 的全扫本身可以接受；catalog spatial index 的首要职责是让每个 source triangle 只枚举相交的固定 1-degree tile，而不是扫描全天。不要仅为把叶子降至数百颗星而预先激进细分；是否需要对极稠密 tile 再做浅层细分，应由实际 Schwarzschild/BBH source-triangle 覆盖和负载 profile 决定。

每颗星获得：

image-plane 连续位置 (x,y)；
frequency shift (g)；
magnification (\mu)。

然后：

[ T_\text{obs}=gT ]

最终用一个 PSF：

[ P(x-x_s,y-y_s) ]

直接 splat 到 HDR framebuffer。

不要先把星压进 pixel 再卷积，因为那会丢失亚像素位置。

PSF 第一版可用 Gaussian；以后可换成 Airy 或其他相机模型。

21. Lens magnification

局部 mapping：

[ (x,y) \leftrightarrow \hat n_\infty ]

自然给出 Jacobian。

点源每阶像的总 flux 应乘：

[ \mu

\left| \frac{d\Omega_\text{image}} {d\Omega_\text{source}} \right| ]

可由：

triangle area ratio；
或局部 interpolation Jacobian；

计算。

critical curve / caustic 附近需要更多 adaptive refinement。

22. 并行策略

22.1 ray evolution

使用 OpenMP 对当前 active ray pool 做 coarse-grained parallel loop。

优先：

#pragma omp parallel for schedule(static, CHUNK)
for (...) {
    evolve_ray_through_current_slab(...);
}

若 strong-field rays 工作量差异明显，再测试：

schedule(dynamic, CHUNK)

当前 movie RayPool 使用 schedule(dynamic, 32)：samples 按 frame 排列，连续静态分块会把同一 time slab 的 active rays 集中到少数 worker。 32 条一块的动态分配让空闲 worker 继续领取工作，也允许不同光线的积分开销不一致。每条 ray 仍由一个 worker 推进，pool slot 与 frame/sample ID 不移动；没有新增共享可变缓存或 per-ray allocation。

这一 chunk 大小来自 393 帧 Schwarzschild 对照实验，并非所有 backend 的最优参数。实验同时比较连续静态、轮转静态和动态分块，用完整 lens-map 文件校验调度前后的结果一致性。当前仍遍历整个 pool； active-index compaction、activation 和 mesh analysis 的并行化属于后续优化。

不要：

一 ray 一个 OpenMP task；
每 ray malloc；
细粒度 mutex。

22.2 adaptive mesh analysis

recursive refinement 逻辑用迭代 work queue实现，不用真正的递归 task。

概念：

while (!queue_empty()) {
    tri = queue_pop();

    if (triangle_is_good(tri)) {
        accept(tri);
    } else {
        request_missing_samples(tri);
    }
}

并行粒度应是：

(frame, root_tile)

或至少按 frame。

每线程使用自己的 local queue 和 local request buffer。

最后 merge：

sort；
deduplicate；
assign vertex IDs；
bulk append。

避免多个 threads 同时修改一个 global mesh。

23. Metric evaluator 的 thread-local workspace

nmesh 单点查询需要：

global x
→ locate AMR leaf
→ global/local coordinate inversion
→ Lagrange interpolation
→ derivatives

应给每线程一个 workspace：

typedef struct {
    uint32_t last_node[...];

    double lx[12];
    double ly[12];
    double lz[12];

    /* scratch arrays */
} MetricWorkspace;

连续 RK stage 中 ray 大概率仍在同一 node，因此优先：

is x still inside last_node?
    yes: reuse
    no: tree / neighbor lookup

不要在 evaluator 内维护共享 mutable cache。

24. 顶层主要数据结构

typedef struct ObserverTrack ObserverTrack;
typedef struct SpacetimeSource SpacetimeSource;
typedef struct MetricSlab MetricSlab;
typedef struct RayPool RayPool;
typedef struct StarCatalog StarCatalog;

typedef struct {
    double t_camera;

    ObserverState observer;

    FrameMesh mesh;
    HDRImage image;
} Frame;

typedef struct {
    Frame *frames;
    size_t nframes;
} Movie;

25. 顶层执行逻辑

概念 main()：

int main(void)
{
    ObserverTrack observer;
    SpacetimeSource spacetime;
    StarCatalog stars;
    Movie movie;

    load_observer_track(&observer, ...);
    open_spacetime(&spacetime, ...);
    load_star_catalog(&stars, ...);

    init_movie(&movie, &observer, ...);
    movie_init_coarse_meshes(&movie);

    for (;;) {
        SampleRequest *req = NULL;

        size_t nreq =
            movie_collect_refinement_requests(
                &movie,
                &req
            );

        if (nreq == 0)
            break;

        RayPool rays;

        rays_build_generation(
            &rays,
            &movie,
            req,
            nreq
        );

        rays_trace_generation(
            &rays,
            &spacetime
        );

        movie_install_ray_results(
            &movie,
            &rays
        );

        free_ray_pool(&rays);
        free(req);
    }

    movie_render_star_catalog(
        &movie,
        &stars
    );

    movie_write_output(&movie);

    return 0;
}

26. `rays_trace_generation()` 的核心流程

void rays_trace_generation(
    RayPool *rays,
    SpacetimeSource *st)
{
    while (spacetime_has_previous_slab(st)) {

        MetricSlab *slab =
            spacetime_load_previous_slab(st);

        activate_rays_for_slab(
            rays,
            slab
        );

        #pragma omp parallel
        {
            MetricWorkspace ws;
            metric_workspace_init(&ws);

            #pragma omp for schedule(static)
            for (size_t i = 0; i < rays->n; ++i) {
                if (!ray_is_active(rays, i))
                    continue;

                evolve_ray_to_slab_left_boundary(
                    rays,
                    i,
                    slab,
                    &ws
                );
            }

            metric_workspace_free(&ws);
        }

        compact_ray_pool(rays);

        spacetime_free_slab(slab);
    }
}

27. repo 的几个主要模块

这里讨论的是代码职责，不强调具体文件树。

`observer`

负责：

读取 worldline；
插值 observer event；
插值/正交化 tetrad；
给每个 frame 提供 observer state。

不负责：

geodesic；
metric；
image rendering。

`movie / frame`

负责：

frame 时间表；
每帧 adaptive image mesh；
vertex / triangle；
refinement passes；
ray sample request；
安装 ray endpoint results。

`ray`

负责：

批量 ray state；
ray 初始化；
active/inactive 状态；
stream compaction；
endpoint 保存。

`geodesic`

负责：

3+1 null geodesic RHS；
frequency shift evolution；
ODE stepper；
在给定 metric slab 中推进 ray。

不负责：

disk I/O；
frame mesh；
star rendering。

`spacetime`

负责统一 backend：

Minkowski；
analytic Schwarzschild；
nmesh 4D numerical spacetime。

负责：

time slab；
metric interpolation；
spatial derivatives；
必要的 temporal derivatives；
capture/infinity classification。

`nmesh_backend`

负责：

读取 nmesh node 文件；
从 filename 恢复 AMR tree；
node spatial lookup；
global/local coordinate conversion；
Lagrange interpolation；
Lagrange derivatives；
temporal interpolation；
puncture trajectory。

`catalog`

负责：

读取 Gaia / 2MASS 等 catalog；
转换为内部 (direction, temperature, amplitude)；
天球 spatial index；
triangle/source-region 查询。

`optics / psf`

负责：

blackbody frequency shift；
temperature → linear RGB；
flux normalization；
PSF splatting；
HDR accumulation。

`output`

负责：

HDR frame；
exposure；
tone mapping；
PNG/EXR；
最终视频编码接口。

28. 开发顺序

Phase 0 — 平直时空 + 星表

目的：

验证 observer/camera projection；
验证天球方向；
验证 Gaia/2MASS catalog；
验证 temperature/amplitude；
验证 PSF；
验证 adaptive image mesh 数据结构。

平直时空：

[ (x,y)\leftrightarrow \hat n_\infty ]

可以作为严格 benchmark。

Phase 1 — analytic Schwarzschild

建议使用 horizon-penetrating coordinates。

验证：

capture；
Einstein ring；
multiple images；
adaptive refinement；
local inverse maps；
magnification；
frequency shift。

Phase 2 — numerical single Schwarzschild

使用当前已能跑的 FOCCZ4 moving-puncture Schwarzschild。

特点：

开始存在 gauge relaxation；
随后 metric 趋于稳定；
可真实测试 4D metric streaming；
可和 analytic Schwarzschild 做 regression comparison。

比较：

[ \Delta \hat n_\infty, \quad g, \quad \text{captured/escaped classification} ]

Phase 3 — FOZ4c BBH

等 FOZ4c：

双曲性；
constraint damping；
数值实现；
BBH 稳定演化；

成熟后接入。

renderer 顶层架构原则上不应为 BBH 重新设计。

29. 尚未决定的问题

需要以后通过 prototype / convergence test 决定：

时间输出 cadence；
time interpolation order；
time slab 最佳内存大小；
BBH production node 数；
4D 输出总数据量；
AH calibration 后 puncture cutoff；
adaptive triangle refinement criterion；
critical curve 附近最大 refinement level；
Gaia 与 2MASS 的最终组合；
PSF 模型；
ODE integrator；
observer trajectory 标准；
tone mapping / exposure 规则；
是否需要 diffuse Milky Way background；
是否将 ray redshift 变量定义为 log(alpha p^0) 或其他更方便的量。

30. 当前最核心的设计原则

整个 repo 可以压缩成以下几条：

最终目标是动态 4D NR 视频渲染，不是单帧 ray tracer。
所有帧的 rays 按 coordinate time 整合，通过 time slabs 倒序推进。
metric backend 可替换，但顶层 scheduler 与具体时空无关。
相机是 worldline + tetrad track，不是固定三维坐标。
恒星是点源 catalog，不是 RGB 天球纹理。
恒星存温度与 amplitude，频移通过温度缩放自然处理。
image-plane 使用 adaptive local triangles 构造局部 inverse lens maps。
adaptive refinement 采用多 pass，避免和 out-of-core metric streaming 冲突。
ray evolution 用 SoA + OpenMP bulk loops。
recursive refinement 用 thread-local iterative queues，不用细粒度 OpenMP tasks。
nmesh 原生 DG node 输出直接作为空间插值数据，不做 Cartesian resampling。
时间导数优先由 3+1 identities、gauge RHS 或 temporal interpolant 解析获得。
所有高开销 mutable cache 都 thread-local。
第一版 CPU-only，先把物理与数据流做正确，再谈 GPU。
从 Minkowski → analytic Schwarzschild → numerical Schwarzschild → BBH 逐级验证。

Schwarzschild 自由落体轨迹辅助脚本

scripts/schwarzschild_camera_track.py 在入射 Cartesian Kerr–Schild 坐标中使用 $g_{\mu\nu}=\eta_{\mu\nu}+(2/r)\ell_\mu\ell_\nu$、$\ell_\mu=(1,x_i/r)$，固定 G=c=M=1 与解析 backend 一致。解析微分 metric 构造四维 Christoffel，以本征时联合积分 dz^\mu/d\tau=u^\mu 与 $de_{(a)}^\mu/d\tau=-\Gamma^\mu_{\alpha\beta}u^\alpha e_{(a)}^\beta$。 e_{(0)}=u 同时满足自由落体方程；四加速度为零时费米–沃克输运就是平行输运。初始坐标速度及指向构造沿用 10.1 的约定，也接受经验证的完整标架。不在积分期间重新正交化以隐藏误差；输出正交归一误差超过 10^{-6} 时失败。

DOP853 默认 rtol=$10^{-10}$、atol=$10^{-12}$，可配置并通过解析径向自由落体、圆轨道和圆轨道平行输运的收敛回归验证。视界不终止相机；默认 r=10^{-3}M 只是可配置的奇点数值保护边界，不等于精确撞击奇点。它独立于光线的 1.5M 捕获 cutoff；该 cutoff 内的轨迹可输出，但当前 renderer 的光线会立即被捕获。

CSV 保持 21 列不变，记录 \tau=k/\mathrm{fps} 与积分得到的真实坐标时间。只保留不超过请求持续本征时或提前终止时刻的规则采样，包含 $\tau=0$。 --movie-track-samples 选择每个 CSV 样本直接生成一帧，绕过坐标时间均匀插值，忽略 renderer 的 start-time/duration/fps。默认 movie 路径不变；两条路径随后共用按真实 coordinate time 的 time-slab 调度，不把本征时冒充坐标时间。

32 KiB Raw Blame History Unescape Escape