fatcuda

English · 中文

English

fatcuda is an intent-first optical tweezers force, torque, and field engine. Users declare the beam, optical system, particle, observable task, and execution policy. The package compiles those declarations into an auditable intermediate representation and dispatches verified CPU, CUDA, or Metal operators.

The current package lives in fatcuda/. The historical direct-port implementation is isolated in fatcuda_old/ and is used only as an oracle, comparison source, and history archive.

Project status

The project is pre-1.0. Its current verified particle surface includes homogeneous and layered spheres, axisymmetric spheroids, and dielectric cap-coated Janus spheres. Version 0.2.0 freezes this verified dielectric-Janus baseline before the next solver/preconditioner workstream. Public PyPI releases exist; APIs may still change before 1.0.

The package is designed for vectorial, non-paraxial optical calculations. It does not treat a visualization as a numerical validation result, and it does not silently replace an unsupported GPU operation with a different physical model.

Design

beam + optical system + particle + task + execution policy
                            |
                            v
                 compile_problem(...)
                            |
                 auditable IR + route plan
                            |
             response construction or store lookup
                            |
            CPU / CUDA / Metal numerical operators
                            |
             force / torque / field + audit trace

The main design rules are:

1. Physical intent is the public boundary. Users declare a problem through solve(...) or compile_problem(...); kernel factories remain available for validation and performance work, not as the recommended application API.
2. The particle response is task-independent. A body-frame T-matrix is built from the optical system and particle, then consumed by force, field, or orientation tasks.
3. Device is execution policy. CPU, CUDA, and Metal share the same problem declarations. The audit records where each operator actually ran and at which precision.
4. Every expensive decision is inspectable. Lowering steps, selected routes, response provenance, cache behavior, numerical order, and backend diagnostics are returned in the audit trace.
5. Oracles remain explicit. MATLAB OTS fixtures, the current NumPy path, and closed-form identities define the verification hierarchy.

Current capability

Area	Current surface
Beam intent	Gaussian, Laguerre-Gauss, Hermite-Gauss, Airy, vortex, calibrated SLM phase, and illumination-driven SLM beams
Focusing	Vectorial non-paraxial Debye-Wolf focusing; dense and FFT multipole coefficient routes
Particles	Sphere, LayeredSphere, Spheroid, JanusSphere, plus the Particle protocol
Responses	Homogeneous/layered Mie, axisymmetric EBCM, and C4-reduced DDA-to-T
Tasks	ForceAtPositions, FieldSlice, and OrientationSweep
Observables	Optical force, optical torque, incident/scattered/total electric fields, orientation-dependent force and torque
Reuse	Content-addressed c128 response store with integrity checks, cross-process locking, and least-recently-used retention
Inspection	Structured progress events, compiled IR, operator plan, route plan, response provenance, and backend diagnostics

Force and torque are evaluated from the multipole-coefficient Maxwell stress tensor. The plane-wave Mie shortcut is used as a validation identity, not as the focused-beam force path.

Backend boundaries

Backend	Response construction	Accelerated execution
CPU	Reference Mie, layered Mie, EBCM, and DDA-to-T paths	Full fp64/c128 reference execution
CUDA	Janus DDA-to-T uses c128 dense cuSOLVER while it fits, then switches to matrix-free FFT-iterative construction	CuPy focusing, coefficient evaluation, force/torque sweeps, field reconstruction, and Debye field paths
Metal	Responses are built on CPU in c128	MLX focusing, coefficient application, sweeps, and field reconstruction; cancellation-sensitive torque work retains its verified c128 CPU boundary

Mie and EBCM responses are device-independent c128 artifacts. A CUDA or Metal run may therefore build a response on CPU and consume it on the requested device; the audit distinguishes build_device from load_device.

Installation

Once a distribution is published:

python -m pip install fatcuda
python -m pip install "fatcuda[metal]"   # Apple Silicon / MLX
python -m pip install "fatcuda[cuda12]"  # CUDA 12 / CuPy

For development from this checkout:

python -m pip install -e ".[test,demo]"
python -m pip install -e ".[test,demo,metal]"
python -m pip install -e ".[test,demo,cuda12]"

The portable default is CPU-only and depends on NumPy and SciPy. Python 3.10 or newer is required.

Quick start

import numpy as np

from fatcuda import (
    ForceAtPositions,
    GaussianBeam,
    OpticalSystem,
    PupilGrid,
    Sphere,
    solve,
)

system = OpticalSystem(
    lambda0=1.064e-6,
    NA=1.20,
    nm=1.33,
    grid=PupilGrid(Nphi=32, Nr=16),
    power=1.0e-3,
)
beam = GaussianBeam(Ex0=1.0, Ey0=0.0)
particle = Sphere(radius=0.2e-6, n_p=1.59, L=8)
task = ForceAtPositions(
    positions=np.array(
        [
            [0.0, 0.0, 0.0],
            [0.1e-6, 0.0, 0.0],
        ]
    )
)

result = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    strategy="auto",
    device="cpu",
)

print(result.force)                    # shape: (2, 3), N
print(result.torque)                   # shape: (2, 3), N m
print(result.audit.operator_plan)
print(result.audit.route_plan)

Use compile_problem(...) when the compiled IR and operator plan should be inspected before execution:

from fatcuda import compile_problem

compiled = compile_problem(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    device="cpu",
)
print(compiled.ir)
print(compiled.audit.operator_plan)
result = compiled.run()

Incident and scattered fields

An incident-only Debye slice needs no particle or T-matrix:

from fatcuda import FieldSlice

incident = solve(
    beam=beam,
    system=system,
    task=FieldSlice(
        plane="xz",
        extent=1.0e-6,
        n_pix=64,
        kind="incident",
    ),
    strategy="debye",
    device="cpu",  # also "cuda" or "metal"
)
print(incident.field.shape)      # (3, 64, 64)
print(incident.intensity.shape)  # (64, 64)

strategy="auto" and "debye" use the faithful Debye-Wolf reconstruction for incident fields. Explicit "direct" and "fft" incident routes are finite-order multipole diagnostics controlled by FieldSlice.L_field. Scattered and total fields require particle= and use the particle response order; the total field combines a Debye incident field with a multipole scattered field.

Progress and audit

Progress is opt-in and synchronous:

def on_progress(event):
    print(event.stage, event.state, event.completed, event.total, event.unit)

result = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    device="cuda",
    progress_callback=on_progress,
    progress_granularity="detailed",  # or "stage"
)

Callbacks do not become part of the task, IR, audit, or result. CUDA detailed events are emitted only after the selected stream events complete.

Expensive responses and caching

solve() defaults to cache="off". Enable the content-addressed response store explicitly when a response will be reused:

first = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    cache="prefer_hit",  # miss -> build -> store
)
again = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    cache="prefer_hit",  # same response -> hit
)

print(again.audit.artifacts["response_cache"]["cache_hit"])

The key includes the quantities that determine the body-frame response, but not particle orientation, requested device, or requested precision. Stored responses are exact c128 arrays with manifest and artifact integrity checks. Inspect the store with:

python -m fatcuda.cache.inventory
python -m fatcuda.cache.inventory --sort build --json

CUDA DDA capacity

The CUDA Janus builder keeps the dense C4/cuSOLVER path while its memory ledger fits. Larger responses automatically use a C4-reduced, matrix-free FFT/BiCGSTAB path. The FFT route solves one RHS at a time: live RTX 4090 measurements found that wider RHS batches consumed more VRAM without a product-relevant speed benefit. It performs no RHS memory-ledger admission; an actual CUDA OOM defines that route's capacity. Projection still accumulates directly from C4 representative moments in bounded site blocks, so its device phase is bounded by Q * projection_site_chunk rather than Q * N.

The largest live-measured RTX 4090 capacity point is:

Quantity	Measured workload
Vacuum wavelength / medium index	1064 nm / 1.33
Janus body / cap used in the run	500 nm radius / 30 nm hemispherical cap
Multipole order	12
Mesh density	135 dipoles per medium wavelength
Lattice spacing	5.926 nm
Dipoles	2,757,332
RHS execution	fixed single RHS
Projection site chunk	16,384
Build time	1,980 s (33.0 min)
Maximum iterative residual	9.997e-13
Measured incremental peak VRAM	22.902 GiB

For this optical system, mesh=135 means (1064 nm / 1.33) / 135 = 5.926 nm. At that lattice spacing, a 15 nm length scale spans about 2.53 cells; the 30 nm cap actually used in the run spans 5.06 cells. This is a qualitative capacity statement, not a thin-film convergence result: quantitative thin-film claims require multiple cells through the coating and an observable-convergence study.

The streamed route has passed its live mesh-40 parity gate and the formerly failing mesh-80 capacity gate. Mesh 135 is the largest measured successful scratch build, not an OOM bracket; its 22.902 GiB incremental peak left about 0.195 GiB of the run's initial free memory. Separately, the mesh-120 c128 response is committed with manifest/NPZ integrity checks and has passed a fresh Mac/CPU CachePolicy.require_hit() without importing CuPy.

Physics and numerical scope

• SI units are used throughout.
• Focusing is vectorial and non-paraxial, with the verified Debye-Wolf pupil mapping and apodization.
• The forward convention is exp(+i k z) with implicit exp(-i omega t).
• Force and torque are optical Maxwell-stress-tensor observables.
• Photophoretic and self-thermophoretic forces are not modeled. They may rival or dominate optical forces for absorbing or metal-capped Janus particles.
• The validated Janus target in 0.2.0 is dielectric-dielectric. For the declared real-Au cap, the current diagonal-preconditioned FFT-iterative DDA route does not reach its production residual at practical meshes. Metal caps remain unsupported until a stronger solver/preconditioner and a separate mesh-convergence campaign land; thermal physics remains out of scope.
• DDA responses are mesh-limited. Use Sphere or LayeredSphere for analytic spherical limits instead of approximating them with DDA.

Demos, numerical QA, and verification

Usage demos and numerical QA are deliberately separate:

python -m fatcuda.demos.run_specific --list
python -m fatcuda.demos.run_specific A1 --quick --device cpu \
  --outdir /tmp/fatcuda-demo

python -m performance_n_qa.run_all --list
python -m performance_n_qa.run_all Q1 --quick --device cpu \
  --outdir /tmp/fatcuda-qa

• fatcuda/demos/ answers “How do I declare and solve this problem?” through the public intent API.
• performance_n_qa/ answers “Is this kernel correct?” and “How fast is it?” and may call lower-level kernels directly.

Run the default test suite with:

python -m pytest -q -p no:cacheprovider

CUDA numerical validation requires a CUDA machine. Metal validation requires Apple Silicon, the metal extra, and an available MLX Metal device.

Repository map and documentation

Path	Role
fatcuda/	Current intent-first package
tests/	Numerical, physics, backend, and architecture gates
fatcuda/demos/	Standalone public-API usage figures
performance_n_qa/	Numerical QA, validation reports, and performance tools
doc/	Project, physics, API, and implementation documentation
ots/	Read-only MATLAB OTS/OTT oracle sources
fatcuda_old/	Self-contained historical archive and comparison oracle

• Documentation index
• Quickstart
• Forward plan
• Detailed code-layer notes

---

中文

fatcuda 是一个 intent-first 的光镊力、力矩与场计算引擎。用户声明光束、光学 系统、粒子、观测任务和执行策略；软件将这些声明编译成可审计的中间表示，再交给 经过验证的 CPU、CUDA 或 Metal 算子执行。

当前生产包位于 fatcuda/。历史 direct-port 实现已隔离到 fatcuda_old/，只作为 oracle、对照来源和历史归档使用。

项目状态

项目尚未进入 1.0。当前已验证的粒子表面包括均匀球、同心分层球、轴对称椭球，以及 介电 cap-coated Janus 球。0.2.0 在下一轮 solver/preconditioner 工作开始前冻结这套 已经验证的介电 Janus 基线。项目已有公开 PyPI release；API 在 1.0 前仍可能调整。

本项目面向矢量、非傍轴光学计算。可视化结果不替代数值验证；当 GPU 不支持某个 算子时，软件也不会静默切换到不同的物理模型。

设计

光束 + 光学系统 + 粒子 + 任务 + 执行策略
                    |
                    v
          compile_problem(...)
                    |
            可审计 IR + 路由计划
                    |
         构建 response 或查询 response store
                    |
          CPU / CUDA / Metal 数值算子
                    |
          力 / 力矩 / 场 + audit trace

核心设计规则：

1. 公开边界是物理意图。 应用代码通过 solve(...) 或 compile_problem(...) 声明问题；kernel factory 主要服务验证与性能工作，不是推荐 的应用入口。
2. 粒子 response 与任务解耦。 body-frame T-matrix 由光学系统和粒子决定，再由 力、场或取向任务消费。
3. device 是执行策略。 CPU、CUDA、Metal 共用同一组问题声明；audit 记录每个 算子的真实执行设备与精度。
4. 昂贵决策都可检查。 lowering、算子路由、response 来源、cache 行为、数值阶数 和 backend diagnostics 都进入 audit trace。
5. oracle 显式存在。 MATLAB OTS fixtures、当前 NumPy 路径和闭式恒等式构成验证 层级。

当前能力

范围	当前表面
光束意图	Gaussian、Laguerre-Gauss、Hermite-Gauss、Airy、vortex、经过标定的 SLM phase，以及 illumination-driven SLM
聚焦	矢量非傍轴 Debye-Wolf；dense 与 FFT 多极系数路线
粒子	Sphere、LayeredSphere、Spheroid、JanusSphere 和 Particle protocol
Response	均匀/分层 Mie、轴对称 EBCM、C4-reduced DDA-to-T
任务	ForceAtPositions、FieldSlice、OrientationSweep
观测量	光学力、光学力矩、incident/scattered/total 电场、取向相关力与力矩
复用	内容寻址的 c128 response store；带完整性检查、跨进程锁和 least-recently-used 保留策略
审计	结构化 progress、compiled IR、operator/route plan、response provenance 和 backend diagnostics

力和力矩由多极系数空间的 Maxwell stress tensor 计算。平面波 Mie shortcut 只作为 验证恒等式，不是 focused-beam force 路径。

Backend 边界

Backend	Response 构建	加速执行
CPU	Mie、layered Mie、EBCM 与 DDA-to-T 的参考实现	完整 fp64/c128 参考路径
CUDA	Janus DDA-to-T 在显存允许时使用 c128 dense cuSOLVER，超过容量墙后自动切换 matrix-free FFT-iterative	CuPy focusing、coefficient evaluation、force/torque sweep、field reconstruction 与 Debye field
Metal	response 在 CPU 以 c128 构建	MLX focusing、coefficient apply、sweep 与 field reconstruction；对 cancellation 敏感的 torque 保留已验证的 CPU c128 边界

Mie 和 EBCM response 是与 device 无关的 c128 artifact。因此 CUDA 或 Metal 任务可能 在 CPU 构建 response，再在指定设备消费；audit 会区分 build_device 与 load_device。

安装

发布 distribution 后：

python -m pip install fatcuda
python -m pip install "fatcuda[metal]"   # Apple Silicon / MLX
python -m pip install "fatcuda[cuda12]"  # CUDA 12 / CuPy

从当前 checkout 开发：

python -m pip install -e ".[test,demo]"
python -m pip install -e ".[test,demo,metal]"
python -m pip install -e ".[test,demo,cuda12]"

默认 portable 路径只依赖 NumPy 和 SciPy，使用 CPU。要求 Python 3.10 或更高版本。

快速开始

下面的完整示例与英文部分相同，可直接运行：

import numpy as np

from fatcuda import (
    ForceAtPositions,
    GaussianBeam,
    OpticalSystem,
    PupilGrid,
    Sphere,
    solve,
)

system = OpticalSystem(
    lambda0=1.064e-6,
    NA=1.20,
    nm=1.33,
    grid=PupilGrid(Nphi=32, Nr=16),
    power=1.0e-3,
)
beam = GaussianBeam(Ex0=1.0, Ey0=0.0)
particle = Sphere(radius=0.2e-6, n_p=1.59, L=8)
task = ForceAtPositions(
    positions=np.array(
        [
            [0.0, 0.0, 0.0],
            [0.1e-6, 0.0, 0.0],
        ]
    )
)

result = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    strategy="auto",
    device="cpu",
)

print(result.force)
print(result.torque)
print(result.audit.operator_plan)
print(result.audit.route_plan)

需要在执行前检查 IR 和 operator plan 时，使用两步式接口：

from fatcuda import compile_problem

compiled = compile_problem(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    device="cpu",
)
print(compiled.ir)
print(compiled.audit.operator_plan)
result = compiled.run()

Incident 与 scattered field

incident-only Debye slice 不需要粒子或 T-matrix：

from fatcuda import FieldSlice

incident = solve(
    beam=beam,
    system=system,
    task=FieldSlice(
        plane="xz",
        extent=1.0e-6,
        n_pix=64,
        kind="incident",
    ),
    strategy="debye",
    device="cpu",  # 也可以是 "cuda" 或 "metal"
)
print(incident.field.shape)
print(incident.intensity.shape)

对 incident field，strategy="auto" 和 "debye" 使用 faithful Debye-Wolf reconstruction；显式 "direct" / "fft" 是由 FieldSlice.L_field 控制的有限阶多极 诊断。scattered/total field 必须传入 particle=，并使用粒子 response order；total field 是 Debye incident 与 multipole scattered field 的和。

Progress 与 audit

progress 是 opt-in、同步执行的：

def on_progress(event):
    print(event.stage, event.state, event.completed, event.total, event.unit)

result = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    device="cuda",
    progress_callback=on_progress,
    progress_granularity="detailed",  # 或 "stage"
)

callback 不进入 task、IR、audit 或 result。CUDA detailed event 只在选定的 stream event 完成后发出。

昂贵 response 与缓存

solve() 默认 cache="off"。需要重复使用 response 时，显式启用内容寻址 store：

first = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    cache="prefer_hit",  # miss -> build -> store
)
again = solve(
    beam=beam,
    system=system,
    particle=particle,
    task=task,
    cache="prefer_hit",  # 相同 response -> hit
)

print(again.audit.artifacts["response_cache"]["cache_hit"])

key 包含决定 body-frame response 的量，但不包含粒子 orientation、requested device 或 requested precision。store 中的 response 是 exact c128 array，并带 manifest 与 artifact 完整性检查。盘点命令：

python -m fatcuda.cache.inventory
python -m fatcuda.cache.inventory --sort build --json

CUDA DDA 容量

CUDA Janus builder 在显存账本允许时保留 dense C4/cuSOLVER；更大的 response 自动使用 C4-reduced matrix-free FFT/BiCGSTAB。RTX 4090 实测表明 wider RHS batch 只增加显存，没有 product-relevant 性能收益，因此 FFT production 固定逐个 RHS 求解，不再用 RHS 显存账本做 admission；真实 CUDA OOM 直接定义容量边界。 projection 仍从 C4 representative moments 按 bounded site block 累加，因此 device phase 的上界是 Q * projection_site_chunk，不再是 Q * N。

RTX 4090 最大的 live 实测成功容量点为：

量	实测 workload
真空波长 / 介质折射率	1064 nm / 1.33
实测 Janus 本体 / cap	500 nm radius / 30 nm hemispherical cap
多极阶数	12
Mesh density	每个介质波长 135 个 dipoles
Lattice spacing	5.926 nm
Dipoles	2,757,332
RHS execution	fixed single RHS
Projection site chunk	16,384
Build time	1,980 s（33.0 min）
最大 iterative residual	9.997e-13
Measured incremental peak VRAM	22.902 GiB

在该光学系统中，mesh=135 表示 (1064 nm / 1.33) / 135 = 5.926 nm。在这一网格间距下，15 nm 长度尺度对应 约 2.53 个 cells；本次实际计算的 30 nm cap 对应 5.06 个 cells。这是 定性容量判断，不是薄膜收敛结论：定量薄膜结论仍需 coating 内有多个网格， 并完成 observable convergence。

streamed route 已通过 live mesh-40 parity gate，以及此前失败的 mesh-80 capacity gate。mesh 135 是最大实测成功的 scratch build，并非 OOM bracket；22.902 GiB incremental peak 后，本次运行开始时的 free memory 只剩约 0.195 GiB。另行持久化 的 mesh-120 c128 response 已通过 manifest/NPZ 完整性检查；新启动的 Mac/CPU 进程可通过 CachePolicy.require_hit() 命中，且不 import CuPy。

物理与数值范围

• 全部使用 SI 单位。
• 聚焦采用经过验证的矢量非傍轴 Debye-Wolf pupil mapping 与 apodization。
• 前向传播约定为 exp(+i k z)，隐含时间因子为 exp(-i omega t)。
• 力和力矩是 optical Maxwell-stress-tensor observables。
• 未建模 photophoresis 和 self-thermophoresis；对吸收或 metal-capped Janus 粒子， 它们可能与光学力相当或占主导。
• 0.2.0 当前已验证的 Janus 目标是 dielectric-dielectric。对于已经声明的 real-Au cap，现有 diagonal-preconditioned FFT-iterative DDA route 在实用 mesh 上无法达到 production residual；在更强的 solver/preconditioner 和独立 mesh-convergence campaign 落地前，metal cap 不属于受支持能力，thermal physics 仍在范围外。
• DDA response 受 mesh 限制。球形解析极限应使用 Sphere 或 LayeredSphere，而不是 用 DDA 近似。

Demo、numerical QA 与验证

usage demo 与 numerical QA 刻意分开：

python -m fatcuda.demos.run_specific --list
python -m fatcuda.demos.run_specific A1 --quick --device cpu \
  --outdir /tmp/fatcuda-demo

python -m performance_n_qa.run_all --list
python -m performance_n_qa.run_all Q1 --quick --device cpu \
  --outdir /tmp/fatcuda-qa

• fatcuda/demos/ 回答“如何用公开 intent API 声明并求解问题”。
• performance_n_qa/ 回答“kernel 是否正确”和“性能如何”，因此可以直接调用底层 kernel。

运行默认测试：

python -m pytest -q -p no:cacheprovider

CUDA 数值验证需要 CUDA 机器。Metal 验证需要 Apple Silicon、metal extra，以及 MLX 能识别的 Metal device。

仓库结构与文档

路径	角色
fatcuda/	当前 intent-first package
tests/	数值、物理、backend 与架构 gates
fatcuda/demos/	独立、使用公开 API 的 usage figures
performance_n_qa/	numerical QA、validation report 与 performance tools
doc/	项目、物理、API 与实现文档
ots/	只读 MATLAB OTS/OTT oracle source
fatcuda_old/	self-contained 历史归档与 comparison oracle

• 文档索引
• Quickstart
• Forward plan
• 详细 code-layer notes