arXiv CSopen access

GPUPhysBench: Benchmarking Coding Agents for Correct and Efficient GPU Physics Simulation

Yuchen Sun et al. · arxiv_cs

arXiv CS · Papers · License: Open Access

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

GPUP HYS B ENCH : B ENCHMARKING C ODING AGENTS FOR C ORRECT AND E FFICIENT GPU P HYSICS S IMU LATION Yuchen Sun, Jinjin He, Sinan Wang, Bo Zhu Georgia Institute of Technology

arXiv:2609.35639v1 [cs.DC] 28 Sep 2026

A BSTRACT Writing fast GPU code for physical simulation is difficult: implementations must preserve numerical accuracy while handling irregular data access, synchronization, and iterative solvers. We introduce GPUPhysBench, a benchmark of 50 tasks testing whether coding agents can meet these demands. Tasks cover fluids, deformable solids, and granular materials, from individual simulation operators to complete simulators. Agents write, compile, test, and optimize GPU code with access to a NVIDIA GPU under fixed time budgets. We report pass rates and runtime performance relative to expert-optimized reference implementations. In a single-attempt evaluation of six frontier model-harness pairs, the two strongest pass all 50 tasks, but even the fastest reaches at least 0.9× the reference speed on only 22% of them, and no submission is more than 5% faster than the reference. The largest gaps arise in collision detection, constraint solving, and iterative solvers. GPUPhysBench brings physical simulation workloads to coding-agent evaluation, testing both the ability to implement numerical methods correctly and the ability to make them run efficiently.

I NTRODUCTION

Modern physical simulation and machine learning rely on GPUs for large-scale workloads, but realizing this performance requires carefully optimized CUDA implementations that demand expertise in numerical algorithms, parallel programming, and hardware-aware optimization. LLM-based coding agents (Chen et al., 2026; Wei et al., 2025) offer an opportunity to automate this labor-intensive process, but evaluating them requires benchmarks that measure both correctness and speed. Physical simulation is an important domain. Beyond long-standing uses in video games (Macklin et al., 2014), aerospace engineering, and manufacturing, it offers a scalable and controllable source of physically plausible interaction data for world models and embodied intelligence (Makoviychuk et al., 2021), where real-world data collection is costly and slow. This demand calls for simulators that are both accurate and efficient, increasing the value of automating their GPU implementation. Existing benchmarks do not cover this setting. Simulation-oriented evaluations emphasize numerical correctness or rely on CPU-based scientific software and established solver libraries (Somasekharan et al., 2025; Hang et al., 2026), while benchmarks for efficient GPU code are dominated by machinelearning operators (Ouyang et al., 2025; Guan et al., 2026); even CUDABench (Zhu et al., 2026) includes only a handful of elementary numerical examples. Unlike regular tensor operators, simulators combine structured grids with particles, meshes, and dynamically evolving neighborhoods, so their bottlenecks include irregular memory access, atomic contention (Gao et al., 2018), neighbor search, and load imbalance. Optimizing them also requires numerical reasoning: agents must preserve the prescribed discretization and stability properties, and in solver-dominated tasks the choice of algorithm and preconditioner determines runtime through convergence. Existing benchmarks are also limited in how they evaluate generated code. Most rely on singleturn generation or short, fixed loops of generation, verification, and profiling. These assess local kernel synthesis but not the ability of modern coding agents to plan, modify files, compile and run programs, diagnose failures, and improve an implementation over an extended trajectory. Evaluating 1

such agents requires a controlled environment that preserves their autonomy while standardizing resource budgets and hardware access. To address both gaps, we introduce GPUPhysBench, a benchmark of 50 coding tasks spanning classical methods for fluid dynamics, deformable solids, and granular materials. Every task includes standardized test data and a reference implementation that human experts have numerically validated and optimized. We evaluate frontier models in full-featured coding-agent harnesses, such as Codex CLI and Claude Code, inside a sandbox that isolates references and evaluators and fixes time budgets and hardware. This design measures pass rates and runtime relative to expert implementations while capturing the long-horizon, tool-using optimization of complete coding-agent systems.

R ELATED W ORK

GPU Kernel Generation A line of benchmarks studies how well LLMs can write highperformance GPU kernels for ML workloads. KernelBench (Ouyang et al., 2025) and MultiKernelBench (Wen et al., 2025) score generated kernels on correctness and speedup over PyTorch baselines. TritonBench (Li et al., 2025), Geak (Wang et al., 2025), and TritonGym (Guan et al., 2026) target the Triton DSL (Tillet et al., 2019). Beyond one-shot generation, agentic systems iteratively optimize kernels using compilation, testing, and profiling feedback. They explore the optimization space through multi-agent refinement loops (Wei et al., 2025; Sun et al., 2026; Zhang et al., 2025), tree search (Dong et al., 2026b; Cao et al., 2026), or evolutionary search with LLM-based variation operators (Chen et al., 2026; Liao et al., 2025; Yoo et al., 2026). LLMs for Scientific Computing Several benchmarks test LLM-written scientific code: SciCode (Tian et al., 2024) curates research coding problems across the natural sciences, while CFDLLMBench (Somasekharan et al., 2025) and PDEAgent-Bench (Hang et al., 2026) check generated CFD and PDE solvers for accuracy and efficiency. Beyond benchmarks, agentic systems apply LLMs to scientific computing workflows: they generate PDE solver code and refine it with execution feedback (Li et al., 2026a; Dong et al., 2026a), automate end-to-end OpenFOAM workflows (Yue et al., 2026; 2025), or tackle neighboring tasks such as PDE control (Soroco et al., 2025), discovery (Luo et al., 2025), and reduced-order modeling (Wang et al., 2026).

B ENCHMARK

3.1

TASK D ESIGN AND C OVERAGE

GPUP HYS B ENCH comprises 50 coding tasks in fluid dynamics, deformable solids, and granular materials, drawn from well-established methods in physical simulation. We select tasks that represent widely used algorithms, pose nontrivial numerical and GPU optimization challenges, and support automatic correctness verification and reliable performance measurement. Using established methods provides clear mathematical specifications and interpretable failure modes while leaving agents substantial freedom in algorithmic and low-level implementation. The fluid tasks use Eulerian grid solvers (Zehnder et al., 2018), particle-in-cell/fluid-implicit-particle (PIC/FLIP) and affine particle-in-cell (APIC) methods (Jiang et al., 2015), position-based fluids (PBF) (Macklin & Müller, 2013), and the lattice Boltzmann method (LBM) (Li et al., 2026b). The solid and granular tasks use mass-spring systems (Liu et al., 2013), the finite element method (FEM) (Sifakis & Barbic, 2012), the material point method (MPM) (Stomakhin et al., 2013), extended position-based dynamics (XPBD) (Macklin et al., 2016), continuous collision detection (CCD) (Brochu et al., 2012), and the discrete element method (DEM) (Lu et al., 2022). Appendix E describes the covered methods and lists all tasks with their categories. The tasks operate on structured grids, particles, lattices, and meshes, covering GPU computation patterns such as stencils, iterative linear solves, atomic scatter and gather, irregular neighborhood interactions, constraint projection, and collision processing. Operator tasks isolate performancecritical operations for fine-grained analysis, while full-simulator tasks require agents to integrate multiple stages into a complete, efficient simulator that preserves numerical behavior. Each task comes with an expert-optimized CUDA reference implementation that defines both the correctness baseline and the performance target. We build the references from open-source, 2

(a) Vortex-ring collision (mc_r)

(b) Water drop (pic_flip)

(d) Jelly cube (fem_explicit)

(e) Snowball (mpm_explicit)

(f) Cloth on a sphere (xpbd)

Figure 1: Multi-step simulations driven by GPUPhysBench reference implementations. high-performance simulation code released with SIGGRAPH papers and courses, such as Fast UAAMG (Shao et al., 2022), GPUMPM (Gao et al., 2018), and an LBM course (Li et al., 2026b); this code is written in C++, earlier versions of CUDA, and NVIDIA Warp (Macklin, 2022). Claude Fable translates it into CUDA, and human experts then rewrite and optimize the result for modern GPUs. Although the upstream code may appear in the pretraining data of the evaluated models, the expert-optimized references are never visible to agents, and the performance gaps in our experiments (Section 4.2) indicate that recalling the upstream code does not by itself reach reference performance. We validate each reference numerically against a serial baseline and run complete multi-step simulations to confirm that it produces physically valid behavior (Figure 1; Appendix G). The references are also competitive with established GPU libraries: on the tasks’ public inputs, they are 9.0× faster than AMGX (Naumov et al., 2015) on the Poisson solve and 2.6–7.0× faster than ports of Warp and Taichi (Hu et al., 2019) examples (Appendix F). 3.2

TASK SPECIFICATION

In each task, an agent implements a specified simulation computation in CUDA and minimizes its GPU execution time subject to numerical correctness requirements (Figure 2). It receives a naturallanguage prompt and a workspace with starter code and a Python driver, while the references and correctness evaluators are withheld. Coding-agent harness. We evaluate agents in production-grade coding harnesses, such as Claude Code, that support file editing and tool use. Prior evaluations may sample or refine over multiple API calls, but each code-generating call emits a complete implementation, both in KernelBench (Ouyang et al., 2025) and in the scripted multi-call workflows evaluated by TritonGym (Guan et al., 2026). In our tests, generating a complete implementation in one call breaks down for large simulators, such as semi-implicit MPM, because the model’s reasoning and code together can exceed the per-call token limit. Persistent harnesses instead let agents build a simulator incrementally, interleaving edits with compilation, testing, and optimization in their own order within the time budget. Task input and environment. The prompt specifies the numerical method, including the governing equations or update rules, discretization, boundary conditions, and convergence criteria where applicable, together with the input and output arrays, scalar parameters, performance objective, and execution constraints. The workspace contains a CUDA source file, an xmake build file, and a 3

(a) Physics tasks

(b) Agent development

GPU workspace Submit </>

Write CUDA

Compile

Correctness

Profile & optimize

Compile

Public + hidden inputs

Coding agent Run & inspect

Task specification CUDA starter + public inputs

Audits for violations

GPU timing Limited execution time Hidden inputs & reference withheld

Hidden inputs

Correctness

Pass rate

fastp

Figure 2: Overview of GPUPhysBench. Agents implement and optimize CUDA physics code under a time budget, without access to hidden inputs, references, or evaluators. Submissions are rebuilt, checked on public and hidden inputs, audited, and timed against expert references. Python driver, and the environment provides a GPU and the tooling to build and run the extension. The driver constructs the public inputs and runs the compiled module, so the agent can inspect inputs, outputs, and timing, but it provides no reference outputs or correctness verdicts. Online access is prohibited. Starter code and required output. The starter code defines a task-specific class whose pybind11 bindings and input()/exec()/output() methods are fixed: they upload the inputs to device buffers, time the computation, and copy the output back. The agent implements allocate(), compute(), and release(), and may add CUDA kernels, internal buffers, and its own data layout; any layout conversion or preprocessing runs inside the timed compute(), while the untimed allocate() may only allocate memory from shapes and scalar parameters. The fixed code and module interface must be preserved, and build changes are limited to compilation flags. Appendix A gives the complete prompt and starter code for a Neumann Poisson task. Numerical requirements. All floating-point computation and storage must use FP32. A submission must pass task-specific checks, such as relative ℓ2 error or solver convergence, on the outputs of its timed calls for both public and hidden inputs. The hidden inputs are withheld during development. They keep the problem size but change the data through different random seeds and initial conditions and, for some tasks, different geometry or solver coefficients (Appendix E.2). Any computational strategy is allowed as long as it meets these accuracy requirements. Correctness tolerances. Most checks compare the submission’s output with the reference output by relative ℓ2 error on the quantity the step changes; for example, particle positions are compared by their displacement over the step, so that large absolute coordinates cannot hide an error in the update. Linear and Newton solves instead require the true residual of the returned solution to fall below the prescribed solver tolerance. Each tolerance is set above the variation between correct FP32 implementations, which differ in rounding, reduction order, and the order of atomic accumulation, and well below the effect of plausible implementation errors, such as a missing term or a flipped sign. Appendices C.5 and C.6 confirm both bounds: correct outputs stay well below their tolerances and erroneous ones far above them, so moderately different thresholds would not change any outcome. Performance objective. The objective is to minimize GPU execution time while satisfying the numerical specification. CUDA events measure all GPU work performed during exec(), including auxiliary computation and solver iterations, while input upload and output download are excluded through the separate interface methods. The generated and reference implementations are evaluated on identical inputs under the same timing protocol. Time budget. The agent receives a wall-clock budget of 30 minutes for operator tasks and 60 minutes for full-simulator tasks, covering implementation, compilation, testing, and optimization. 4

The prompt states the budget and absolute deadline. The agent may independently check the current time using the date command and compare it with the deadline to decide whether to continue optimizing or finalize its submission. At most one agent, including the main agent and any subagents, may be active at a time. Any delegated execution shares the same task-level wall-clock budget. The run is terminated when the budget expires, and the submitted source and build configuration present when the agent finishes or reaches the deadline are used for evaluation. 3.3

A NTI - CHEATING

File-system isolation. Each task runs in a Docker container, in a separate working directory that contains only the agent-facing files. The agent runs as an unprivileged user whose permissions block access to the benchmark project, including references and evaluators, and the harness verifies this isolation before launch. For scoring, the harness restores the fixed starter-code sections and rebuilds the submission in a clean directory with the protected driver and evaluator, ignoring any prebuilt module left by the agent. Execution rules and post-run auditing. Agents may not access online resources, run agents in parallel, or move computation outside the timed region; built-in web search is disabled where the harness supports it. Static checks enforce the interface, build, and allocation restrictions, and Claude-Opus-5 audits each run’s logs and code for network access, parallel agents, and untimed computation, using the same prompts for every system (Appendix B). Any violation counts as a failure, regardless of correctness or speed. 3.4

M ETRICS

We evaluate agents using correctness, pass rate, and performance relative to the reference. All metrics are computed over all N benchmark tasks using the final submission from each task attempt. Correctness. Correctness is the fraction of tasks whose submission builds, runs, and passes all numerical checks on both public and hidden inputs, regardless of audit outcomes. Pass Rate. Pass rate is the fraction of tasks whose submission is correct and also passes all anticheating audits; we write vi = 1 for such a passing submission and vi = 0 otherwise. Performance. Following KernelBench (Ouyang et al., 2025), we use fastp to measure the fraction of tasks that pass and achieve a speedup greater than a threshold p. For a passing submission (vi = 1), its speedup is tref i si = gen , (1) ti where tref and tgen are the GPU execution times of the reference and generated implementations, i i measured on the hidden inputs using the same hardware and timing protocol. We set si = 0 for failed submissions. For p ≥ 0, N 1 X fastp = vi 1[si > p]. (2) N i=1 We report p ∈ {0.5, 0.9, 1.05}. Failed submissions remain in the denominator. fast1.05 counts passing implementations that outperform the reference by more than 5%. We use this threshold instead of p = 1 because a speedup just above 1 is within timing noise. Each unchanged reference is timed 13–14 times across our evaluation sessions, and these timings vary with a median coefficient of variation of 0.9% per task (at most 3.0%) and a median max-to-min ratio of 1.03.

E XPERIMENTS

We evaluate coding agents on GPUPhysBench to assess their ability to produce numerically correct and efficient GPU simulation programs. We compare six model-harness pairs against expert implementations in terms of correctness, pass rate, and runtime performance, and analyze the generated programs to identify the factors that contribute to the remaining performance gap. 5

Table 1: Results on GPUPhysBench over 50 tasks, with one attempt per task. The best result in each column is underlined. Model

Harness

Claude-Opus-5 GPT-5.6-Sol Gemini-3.5-Flash DeepSeek-V4.1-Flash Qwen-3.8-Max GLM-5.3

Claude Code Codex CLI Gemini CLI DeepSeek Harness Qwen Code OpenCode

4.1

Correctness ↑ Pass Rate ↑ fast0.5 ↑ fast0.9 ↑ fast1.05 ↑ 100% 100% 88% 86% 52% 88%

100% 100% 88% 86% 52% 86%

66% 42% 28% 38% 34% 34%

22% 16% 16% 16% 14% 16%

0% 0% 0% 0% 0% 0%

S ETUP

All agent development and evaluation run on a single NVIDIA GeForce RTX 4090 (Ada Lovelace), so agents tune on the same GPU that scores them, inside an Ubuntu-based CUDA Docker image with Python, pybind11, and xmake. For each implementation and input set, we report the minimum of five timed calls, each in a fresh process after three warm-up calls on the other input set; speedups use the hidden-input timings. We evaluate six frontier models, using their corresponding coding harnesses where possible. We pair Claude-Opus-5 with Claude Code, GPT-5.6-Sol with Codex CLI, Gemini-3.5-Flash with Gemini CLI, DeepSeek-V4.1-Flash with DeepSeek Harness, and Qwen-3.8-Max with Qwen Code. For GLM-5.3, we use OpenCode instead of ZCode, since the latter is a desktop development environment rather than a CLI harness. We refer to the six systems by their model family: Opus, GPT, Gemini, DeepSeek, Qwen, and GLM. We configure reasoning effort to high for all applicable harnesses except Gemini CLI, which does not expose a corresponding effort-level setting. Each system receives one attempt per task in Table 1; Appendix C.4 reports development time and token usage. 4.2

M AIN R ESULT

Table 1 summarizes the performance of six model-harness pairs across all 50 GPUPhysBench tasks. These results assess complete coding systems: each agent must translate a numerical specification into CUDA code, resolve implementation issues, and optimize execution within a fixed time budget. We report numerical correctness separately from audit-compliant success and execution efficiency. Observation 1: Frontier models can correctly implement fully specified numerical methods. Opus and GPT both pass all numerical checks and audits on all 50 tasks, including full simulators that require multiple numerical stages to work together, achieving 100% correctness and pass rate. Gemini reaches 88% on both metrics, DeepSeek 86%, and Qwen 52%, while GLM reaches 88% correctness and an 86% pass rate. This result should be read in light of the task format: each prompt prescribes the numerical method, discretization, boundary conditions, and convergence criteria, so the tasks test faithful implementation of a given method rather than the choice of physical model. For the strongest systems, correctness is therefore close to saturation, and performance is the axis on which GPUPhysBench separates them. Observation 2: Simulation performance still lags behind expert references. Despite matching Opus in correctness and pass rate, GPT achieves fast0.5 = 42%, compared with Opus’s 66%, revealing a substantial gap in execution efficiency; Opus is faster than GPT on 38 of the 50 tasks. Even for Opus, 17 of the 50 passing submissions take at least twice as long as the reference. Opus leads at fast0.9 = 22%, while the other systems reach 14%–16%; their fast0.9 successes come only from regular, memory-bound local grid and lattice computations, on which nearly all systems come close to the reference, whereas Opus also approaches the reference on a few particle and mesh tasks with irregular access. No passing submission is more than 5% faster than the reference (fast1.05 = 0 for every system). Thus, reliable numerical implementation does not yet translate into performance comparable to expert references on most tasks.

4.3

P ERFORMANCE BY TASK C ATEGORY

To compare benchmark outcomes across computational structures, we group tasks Global Solves and into six categories according to their nuPressure Projection Local Grid and Lattice merical methods and computational strucComputations 18% ture (Figure 3). Local Grid and Lat28% tice Computations covers advection, LBM Geometric Queries and 10% Collision Detection operations, and local surface-tension and phase-field updates. Particle-Grid Meth10% Position Constraint ods includes transfer operators and com22% Solving 12% plete PIC/FLIP, APIC, and MPM steps. Particle–Grid Methods Local Interaction Local Interaction Updates covers explicit Updates elasticity, contact-force evaluation, and local viscosity and vorticity velocity updates. Position Constraint Solving con- Figure 3: Distribution of the 50 tasks across six comtains PBF and XPBD constraint projec- putational categories. tions and complete steps. Geometric Queries and Collision Detection comprises CCD and particle-based distance-field construction. Global Solves and Pressure Projection includes Poisson and implicit viscosity solves, Newton and projective-dynamics elasticity, and grid-based fluid steps with pressure projection. Figure 4 reports, for each system and category, the geometric mean speedup over the tasks the system passes, so speed is reported separately from success: Qwen’s high mean on global solves covers only 2 of 9 tasks. Every system is fastest on local grid and lattice computations, and Opus leads in five of the six categories; on global solves, Opus, GPT, and DeepSeek lie between 0.23× and 0.25×. Geometric queries and collision detection are the slowest on average, followed by global solves and position constraint solving, while local interaction updates vary the most across systems. Weighting each task family equally (Appendix C.3) leaves the ordering by fast0.5 unchanged but lowers the fast0.9 of every system other than Opus to 8%, because their successes concentrate in the LBM family. Appendix D examines two cases of generated code in detail. Overall performance across categories. Agents perform best on local computations, particularly regular grid and lattice operations. All passing submissions with speedups above 1× belong to the local grid and lattice category, and their gains are below 0.3%, with a median speedup of 1.001× and a maximum of 1.002×. These margins are smaller than the run-to-run variation of the timings (Section 3.4) and do not establish a performance advantage. Global solves and more complex workloads involving irregular data access, iterative updates, or coordination across multiple stages show larger performance gaps. Overall, agents are more successful at exploiting regular local parallelism than at jointly optimizing algorithmic choices, data organization, and execution across an entire simulation step. Local Grid and Lattice Computations. Near-reference performance is concentrated in LBM, streaming, collision, and surface-tension tasks, where regular indexing and local arithmetic are readily mapped to GPU threads. Advection and Cahn–Hilliard updates show larger gaps and more variation across systems. In third-order advection, for example, Opus processes all three velocity components in one kernel, whereas GPT launches separate component kernels. These implementations highlight opportunities to reuse interpolation data and combine stages even within otherwise regular workloads. Particle-Grid Methods. Spatial ordering alone does not ensure efficient transfers. In PIC/FLIP P2G, Gemini sorts particle indices but retains indirect particle reads and global atomics for each contribution. Opus packs particles into spatial tiles and accumulates in shared memory, reaching 0.556× the reference speed versus Gemini’s 0.134×. The reference further aggregates contributions by cell before merging them into shared memory (Appendix D.1). The remaining gap therefore involves both the cost of grouping particles and the granularity of accumulation; sorting or moving atomics into shared memory is only part of the optimization. 7

0.75

0.86 13/14

0.76 14/14

0.89 11/14

0.86 12/14

0.50 0.25 0.00

e De

0.75 0.23 5/5

0.25 0.00

u Op

T GP

0.16 3/5

0.15 5/5

e Se

p ee

0.30 2/5

e Qw

0.19 5/5

M GL

(d) Position Constraint Solving

0.56 11/11

0.50

0.30 11/11

0.25 0.00

0.32 7/11

0.13 8/11

m Ge

ep De

0.23 8/11

0.22 10/11

0.75 0.33 5/5

0.25 0.00

0.08 5/5

u Op

T GP

0.14 5/5

0.11 3/5

ek Se

m Ge

ep De

0.17 5/5 0/5

en Qw

0.75

(e) Geometric Queries and Collision Detection

0.80 6/6 0.57 6/6

0.50

0.30 6/6

0.25 0.00

0.65 3/6 0.27 5/6

0.15 6/6

m Ge

ep De

1.00

0.50

1.00

(b) Particle-Grid Methods

Geo. mean speedup (×)

1.00

0.44 5/5

0.75

(a) Local Grid and Lattice Computations

0.50

1.00

Geo. mean speedup (×)

0.82 14/14

Geo. mean speedup (×)

0.93 14/14

Geo. mean speedup (×)

1.00

1.00 0.75 0.50 0.25

0.23 9/9

0.00

0.24 9/9

0.07 8/9

0.25 9/9

m Ge

ep De

0.37 2/9

0.20 6/9

(f) Global Solves and Pressure Projection

Figure 4: Geometric mean speedup relative to the expert reference on hidden inputs, over each system’s passing tasks in each category. Labels give the mean and the number of passing tasks; a system with no passing task has no bar. Axis scales are shared across panels.

Local Interaction Updates. Performance varies substantially even among correct implementations of local interactions. Opus approaches reference performance on explicit FEM and matches it closely on DEM contact, while other systems often remain much slower. In DEM contact, both Opus and GPT pack particle records, but Opus traverses contiguous ranges in a spatial grid, whereas GPT searches hash buckets and filters candidates by cell coordinates. These differences highlight the importance of neighbor-search organization and candidate access beyond the arithmetic of the local force or velocity update. Position Constraint Solving. PBF implementations already exploit spatial reordering: Opus rebuilds cell offsets and rearranges particles during each constraint iteration, reaching 0.815× reference speed on the standalone PBF solve. Spring-based XPBD remains harder, with every system below 0.3× on the standalone constraint solve. Opus accumulates spring corrections through atomics into packed particle buffers; GPT instead builds adjacency lists and gathers incident corrections. Both repeatedly exchange corrections and positions through global memory, leaving opportunities to reduce iteration traffic and improve reuse across constraints. Geometric Queries and Collision Detection. All but one passing continuous-collision-detection submission remains below 0.5× reference speed, with the best reaching 0.526×, despite using spatial acceleration. The implementations span different structures: Opus uses a spatial grid for vertex– vertex CCD, while GPT builds and traverses a bounding-volume hierarchy yet reaches only 0.007×; Gemini’s vertex–face implementation sorts face–cell overlap pairs. These choices introduce different construction costs, candidate lists, and traversal patterns that must be optimized together. Particle distance-field construction generally performs better than CCD, but still trails the reference, suggesting that efficient candidate pruning and geometry access remain important across this category. Global Solves and Pressure Projection. Category averages hide substantial differences in solver choice. Opus and Gemini both use Jacobi (diagonal) preconditioning for the Dirichlet Poisson solve, which Opus applies by symmetric rescaling; they reach only 0.049× and 0.036× reference speed, respectively. Opus fuses iteration kernels and keeps CG coefficients on the device, but these optimizations leave a large gap (Appendix D.2). GPT, GLM, and DeepSeek use multigrid precondition8

Table 2: Ablations of time budget and number of attempts. Budgets are relative to the default of 30 minutes for operators and 60 minutes for full simulators. With three attempts, a task counts as solved if any run solves it, and fastp uses the best passing speedup per task. Attempts Time Budget Coding Agent Correctness ↑ Pass Rate ↑ fast0.5 ↑ fast0.9 ↑ fast1.05 ↑

0.5×

Opus GPT

98% 100%

54% 36%

18% 16%

0% 0%

1.0×

Opus GPT

100% 100%

66% 42%

22% 16%

0% 0%

1.5×

Opus GPT

100% 100%

76% 40%

34% 18%

4% 0%

1.0×

Opus GPT

100% 100%

76% 48%

36% 18%

2% 0%

ers on the same task and reach 0.241×, 0.248×, and 0.303×, respectively. All three still trail the reference, underscoring the need to optimize convergence, memory traffic and synchronization. 4.4

A BLATION S TUDY

Time budget. We evaluate Opus and GPT with one attempt at 0.5×, 1.0×, and 1.5× the default time budget (Table 2). As the budget grows, Opus’s fast0.5 rises steadily (54%, 66%, and 76%), whereas GPT’s changes little (36%, 42%, and 40%). Correctness and pass rate remain at 100% except for Opus at 0.5×, where both are 98%. At 1.5×, two Opus submissions are more than 5% faster than the reference, on DEM (1.11×) and MPM grid-to-particle transfer (1.07×). Thus, extra time primarily improves optimization, with larger gains for Opus in these runs. The default budget is rarely binding: at 1.0×, Opus and GPT use a median of 43% and 50% of it and never reach the deadline (Appendix C.4). The gains at 1.5× therefore reflect how agents choose to use a longer stated deadline more than a lack of time. Number of attempts. At the default budget, we run each model three times and compare the first run with the best of three, which keeps, for each task, the run with the highest hidden-input speedup among those that pass numerical checks and all audits. Correctness and pass rate reach 100% for both models, although single runs pass 50, 49, and 50 tasks for Opus and 50, 47, and 49 for GPT. fast0.5 rises from 66% to 76% for Opus and from 42% to 48% for GPT, well beyond the run-to-run standard deviation of single-attempt fast0.5 (1.2 and 2.3 points; Appendix C.2). fast0.9 rises from 22% to 36% for Opus, whereas GPT’s increase from 16% to 18% is within run-to-run variation, since its second run alone reaches 18%. Among runs that pass all numerical checks and audits, selecting by public-input speedup gives the same fast0.5 and fast0.9 . This post-hoc comparison uses evaluator-only information to identify passing runs and compute reference-relative speedups. Only one run, from Opus on DEM (1.13×), beats the reference by more than 5% (fast1.05 of 2% vs. 0% for GPT). Additional attempts therefore improve peak performance rather than task coverage.

C ONCLUSION

We introduced GPUPhysBench, a benchmark of 50 GPU physics simulation tasks spanning individual operators and complete simulation steps. By evaluating coding agents in interactive development environments against numerical checks and expert reference implementations, GPUPhysBench measures both implementation correctness and execution efficiency. Our results show that frontier agents can correctly implement fully specified simulation methods, with the two strongest systems passing all tasks, while substantial performance gaps remain, particularly for collision detection, constraint solving, and iterative solvers. Analysis of generated code and development traces highlights the importance of numerical algorithm selection, data organization, and coordination across kernels. These findings motivate agents that reason jointly about numerical methods and GPU execution, and establish GPUPhysBench as a testbed for progress toward automated implementation of efficient physical simulators. 9

AI U SE S TATEMENT We used generative AI tools to assist in constructing the benchmark, to draft sections of the paper, and to aid and polish the writing. In building the reference implementations, Claude Fable translated the upstream simulation code into CUDA, and human experts then rewrote and optimized the result. The authors carefully reviewed all AI-assisted work, including the benchmark tasks, the experimental results, and the text of the paper, and take full responsibility for the content of this work.

R EPRODUCIBILITY S TATEMENT We will release GPUPhysBench under an open-source license, including all 50 task prompts, starter code, drivers, reference implementations, correctness evaluators, the Docker environment, and the evaluation and audit scripts. Appendix A gives the complete prompt and starter code of one task, Appendix E.2 the input sizes and correctness criteria of every task, and Appendix B the audit prompts.

R EFERENCES Tyson Brochu, Essex Edwards, and Robert Bridson. Efficient geometrically exact continuous collision detection. ACM Trans. Graph., 31(4), 2012. Shiyi Cao, Ziming Mao, Joseph E. Gonzalez, and Ion Stoica. K-search: Llm kernel generation via co-evolving intrinsic world model. arXiv preprint arXiv:2602.19128, 2026. Terry Chen, Zhifan Ye, Bing Xu, Zihao Ye, Timmy Liu, Ali Hassani, Tianqi Chen, Andrew Kerr, Haicheng Wu, Yang Xu, Yu-Jung Chen, Hanfeng Chen, Aditya Kane, Ronny Krashinsky, MingYu Liu, Vinod Grover, Luis Ceze, Roger Bringmann, John Tran, Wei Liu, Fung Xie, Michael Lightstone, and Humphrey Shi. Avo: Agentic variation operators for autonomous evolutionary search. arXiv preprint arXiv:2603.24517, 2026. Huanshuo Dong, Keyao Zhang, Hong Wang, Zhezheng Hao, Zhiwei Zhuang, Ziyan Liu, Jiacong Wang, Gengyuan Liu, and Xin Jin. Autopde: Reliable agentic pde solving via explicitly represented solver strategies. arXiv preprint arXiv:2606.10752, 2026a. Juncheng Dong, Yang Yang, Tao Liu, Yang Wang, Feng Qi, Vahid Tarokh, Kaushik Rangadurai, and Shuang Yang. Stark: Strategic team of agents for refining kernels. In ICLR, 2026b. Ming Gao, Xinlei Wang, Kui Wu, Andre Pradhana, Eftychios Sifakis, Cem Yuksel, and Chenfanfu Jiang. Gpu optimization of material point methods. ACM Trans. Graph., 37(6), 2018. Yue Guan, Yichen Lin, Xu Zhao, Jianzhu Yao, Xinwei Qiang, Zhongkai Yu, Pramod Viswanath, Yufei Ding, and Adnan Aziz. Tritongym: A benchmark for agentic llm workflows in triton gpu code generation. In ICML, 2026. Zhen Hang, Yushan Yashengjiang, Junhui Li, Huanshuo Dong, Yang Wei, Zhezheng Hao, Jiangtao Ma, Songlin Bai, Haozhong Kai, Xihang Yue, Gangzong Si, Dongming Jiang, Chao Yao, Zhanhua Hu, Jiangqing Zhang, Pengwei Liu, Yaomin Shen, Xingyu Ren, Lei Liu, Zikang Xu, Han Li, Qingsong Yao, Hande Dong, and Hong Wang. Pdeagent-bench: A multi-metric, multi-library benchmark for pde solver generation. arXiv preprint arXiv:2605.09636, 2026. Yuanming Hu, Tzu-Mao Li, Luke Anderson, Jonathan Ragan-Kelley, and Frédo Durand. Taichi: a language for high-performance computation on spatially sparse data structures. ACM Trans. Graph., 38(6), 2019. Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin. The affine particle-in-cell method. ACM Trans. Graph., 34(4), 2015. Jianling Li, Shangzhan Li, Zhenye Gao, Qi Shi, Yuxuan Li, Zefan Wang, Jiacheng Huang, Haojie Wang, Jianrong Wang, Xu Han, Zhiyuan Liu, and Maosong Sun. Tritonbench: Benchmarking large language model capabilities for generating triton operators. arXiv preprint arXiv:2502.14752, 2025. 10

Shanda Li, Tanya Marwah, Junhong Shen, Weiwei Sun, Andrej Risteski, Yiming Yang, and Ameet Talwalkar. Codepde: An inference framework for llm-driven pde solver generation. TMLR, 2026a. Wei Li, Chaoyang Lyu, Mengyun Liu, Yixin Chen, Mathieu Desbrun, Kui Wu, and Xiaopei Liu. Fluid simulation with the lattice boltzmann method. In Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques Conference Courses, 2026b. Gang Liao, Hongsen Qin, Ying Wang, Alicia Golden, Michael Kuchnik, Yavuz Yetim, Jia Jiunn Ang, Chunli Fu, Yihan He, Samuel Hsia, Zewei Jiang, Dianshi Li, Uladzimir Pashkevich, Varna Puvvada, Feng Shi, Matt Steiner, Ruichao Xiao, Liyuan Li, Nathan Yan, Xiayu Yu, Zhou Fang, Roman Levenstein, Kunming Ho, Haishan Zhu, Alec Hammond, Richard Li, Ajit Mathews, Kaustubh Gondkar, Abdul Zainul-Abedin, Ketan Singh, Hongtao Yu, Wenyuan Chi, Barney Huang, Sean Zhang, Noah Weller, Zach Marine, Wyatt Cook, Carole-Jean Wu, and Gaoxiang Liu. Kernelevolve: Scaling agentic kernel coding for heterogeneous ai accelerators at meta. arXiv preprint arXiv:2512.23236, 2025. Tiantian Liu, Adam W. Bargteil, James F. O’Brien, and Ladislav Kavan. Fast simulation of massspring systems. ACM Trans. Graph., 32(6), 2013. Jia-Ming Lu, Chen-Feng Li, Geng-Chen Cao, and Shi-Min Hu. Simulating fractures with bonded discrete element method. IEEE Transactions on Visualization and Computer Graphics, 28(12), 2022. Xiao Luo, Changhu Wang, Yizhou Sun, and Wei Wang. How do large language models perform on PDE discovery: A coarse-to-fine perspective. In Findings of the Association for Computational Linguistics, 2025. Miles Macklin. Warp: A high-performance python framework for gpu simulation and graphics, March 2022. NVIDIA GPU Technology Conference (GTC). Miles Macklin and Matthias Müller. Position based fluids. ACM Trans. Graph., 32(4), 2013. Miles Macklin, Matthias Müller, Nuttapong Chentanez, and Tae-Yong Kim. Unified particle physics for real-time applications. ACM Trans. Graph., 33(4), 2014. Miles Macklin, Matthias Müller, and Nuttapong Chentanez. Xpbd: position-based simulation of compliant constrained dynamics. In Proceedings of the 9th International Conference on Motion in Games, 2016. Viktor Makoviychuk, Lukasz Wawrzyniak, Yunrong Guo, Michelle Lu, Kier Storey, Miles Macklin, David Hoeller, Nikita Rudin, Arthur Allshire, Ankur Handa, and Gavriel State. Isaac gym: High performance gpu-based physics simulation for robot learning. arXiv preprint arXiv:2108.10470, 2021. M. Naumov, M. Arsaev, P. Castonguay, J. Cohen, J. Demouth, J. Eaton, S. Layton, N. Markovskiy, I. Reguly, N. Sakharnykh, V. Sellappan, and R. Strzodka. Amgx: A library for gpu accelerated algebraic multigrid and preconditioned iterative methods. SIAM J. Sci. Comput., 2015. Anne Ouyang, Simon Guo, Simran Arora, Alex L. Zhang, William Hu, Christopher Ré, and Azalia Mirhoseini. Kernelbench: Can llms write efficient gpu kernels? In ICML, 2025. Han Shao, Libo Huang, and Dominik L. Michels. A fast unsmoothed aggregation algebraic multigrid framework for the large-scale simulation of incompressible flow. ACM Trans. Graph., 41(4), 2022. Eftychios Sifakis and Jernej Barbic. Fem simulation of 3d deformable solids: a practitioner’s guide to theory, discretization and model reduction. In ACM SIGGRAPH 2012 Courses, 2012. Nithin Somasekharan, Ling Yue, Yadi Cao, Weichao Li, Patrick Emami, Pochinapeddi Sai Bhargav, Anurag Acharya, Xingyu Xie, and Shaowu Pan. Cfdllmbench: A benchmark suite for evaluating large language models in computational fluid dynamics. arXiv preprint arXiv:2509.20374, 2025. 11

Mauricio Soroco, Jialin Song, Mengzhou Xia, Kye Emond, Weiran Sun, and Wuyang Chen. Pdecontroller: Llms for autoformalization and reasoning of pdes. In ICML, 2025. Alexey Stomakhin, Craig Schroeder, Lawrence Chai, Joseph Teran, and Andrew Selle. A material point method for snow simulation. ACM Trans. Graph., 32(4), 2013. Qitong Sun, Jun Han, Tianlin Li, Zhe Tang, Sheng Chen, Fei Yang, Aishan Liu, Xianglong Liu, and Yang Liu. Kernelskill: A multi-agent framework for gpu kernel optimization. arXiv preprint arXiv:2603.10085, 2026. Minyang Tian, Luyu Gao, Shizhuo Dylan Zhang, Xinan Chen, Cunwei Fan, Xuefei Guo, Roland Haas, Pan Ji, Kittithat Krongchon, Yao Li, Shengyan Liu, Di Luo, Yutao Ma, Hao Tong, Kha Trinh, Chenyu Tian, Zihan Wang, Bohao Wu, Yanyu Xiong, Shengzhu Yin, Minhui Zhu, Kilian Lieret, Yanxin Lu, Genglin Liu, Yufeng Du, Tianhua Tao, Ofir Press, Jamie Callan, Eliu Huerta, and Hao Peng. Scicode: A research coding benchmark curated by scientists. In NeurIPS, 2024. Philippe Tillet, H. T. Kung, and David Cox. Triton: an intermediate language and compiler for tiled neural network computations. In Proceedings of the 3rd ACM SIGPLAN International Workshop on Machine Learning and Programming Languages, 2019. Jianghui Wang, Vinay Joshi, Saptarshi Majumder, Xu Chao, Bin Ding, Ziqiong Liu, Pratik Prabhanjan Brahma, Dong Li, Zicheng Liu, and Emad Barsoum. Geak: Introducing triton kernel ai agent & evaluation benchmarks. arXiv preprint arXiv:2507.23194, 2025. Zhuoyuan Wang, Hanjiang Hu, Xiyu Deng, Saviz Mowlavi, and Yorie Nakahira. Opinf-llm: Parametric pde solving with llms via operator inference. arXiv preprint arXiv:2602.01493, 2026. Anjiang Wei, Tianran Sun, Yogesh Seenichamy, Hang Song, Anne Ouyang, Azalia Mirhoseini, Ke Wang, and Alex Aiken. Astra: A multi-agent system for gpu kernel performance optimization. arXiv preprint arXiv:2509.07506, 2025. Zhongzhen Wen, Yinghui Zhang, Zhong Li, Zhongxin Liu, Linna Xie, and Tian Zhang. Multikernelbench: A multi-platform benchmark for kernel generation. arXiv preprint arXiv:2507.17773, 2025. Jason Yoo, Rajarshi Saha, Shaowei Zhu, Tao Yu, Wei Tang, and Youngsuk Park. Mkevolve: A modular multi-agent framework for kernel code generation. arXiv preprint arXiv:2607.20501, 2026. Ling Yue, Nithin Somasekharan, Tingwen Zhang, Yadi Cao, and Shaowu Pan. Foam-agent 2.0: An end-to-end composable multi-agent framework for automating cfd simulation in openfoam. arXiv preprint arXiv:2509.18178, 2025. Ling Yue, Nithin Somasekharan, Tingwen Zhang, Yadi Cao, Zhangze Chen, Shimin Di, and Shaowu Pan. Foam-agent: A large language model-based multi-agent framework for automating computational fluid dynamics workflows. Computer Methods in Applied Mechanics and Engineering, 2026. Jonas Zehnder, Rahul Narain, and Bernhard Thomaszewski. An advection-reflection solver for detail-preserving fluid simulation. ACM Trans. Graph., 37(4), 2018. Genghan Zhang, Shaowei Zhu, Anjiang Wei, Zhenyu Song, Allen Nie, Zhen Jia, Nandita Vijaykumar, Yida Wang, and Kunle Olukotun. Accelopt: A self-improving llm agentic system for ai accelerator kernel optimization. arXiv preprint arXiv:2511.15915, 2025. Jiace Zhu, Wentao Chen, Qi Fan, Zhixing Ren, Junying Wu, Xing Zhe Chai, Chotiwit Rungrueangwutthinon, Yehan Ma, and An Zou. Cudabench: Benchmarking llms for text-to-cuda generation. arXiv preprint arXiv:2603.02236, 2026.

A PPENDIX C ONTENTS A Example Task Prompt and Starter Code

A.1 Task Prompt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.2 Initial CUDA Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.3 Python Driver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B Auditor

C Detailed Results

C.1 Per-Task Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.2 Run-to-Run Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.3 Family-Balanced Scores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.4 Development Time and Token Usage . . . . . . . . . . . . . . . . . . . . . . . . .

C.5 Correctness Margins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.6 Validation of the Tolerances . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

D Case Studies of Generated Implementations

D.1 PIC/FLIP Particle-to-Grid Transfer: Organizing Accumulation . . . . . . . . . . .

D.2 Dirichlet Poisson Solve: Iteration Cost and Solver Choice . . . . . . . . . . . . . .

E Benchmark Tasks

E.1 Task List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.2 Task Inputs and Correctness Checks . . . . . . . . . . . . . . . . . . . . . . . . .

F Comparison with GPU Libraries

F.1

Poisson Solve with Dirichlet Boundaries versus AMGX . . . . . . . . . . . . . . .

F.2

Simulation Steps versus Warp and Taichi . . . . . . . . . . . . . . . . . . . . . . .

G Simulations Built on Reference Implementations

E XAMPLE TASK P ROMPT AND S TARTER C ODE

We use the Neumann Poisson task to illustrate the agent’s specification and CUDA interface. The complete prompt below reproduces the task specification, including the requirement that all algorithmic computation occur inside the timed compute(), and the constraints appended by the evaluation harness. Only run-specific deadlines are replaced by placeholders. Listing 1 reproduces the complete initial CUDA source, including its interface documentation, fixed interface methods, empty implementation hooks, timing scaffold, and pybind11 bindings. Appendix A.3 describes the Python driver that builds the task’s inputs and runs a compiled module. A.1

TASK P ROMPT

Neumann Poisson Prompt <TASK> Implement a preconditioned Conjugate Gradient solver in CUDA, exposed to Python as a pybind11 extension module. Use xmake for compilation. Use FP32 (single-precision floating point) throughout the entire implementation. Do not use FP64 (double-precision floating point) for any computation or storage. Make the implementation as fast as you can: its GPU time is measured and reported, so optimize the CUDA code for performance. Do not tailor the implementation or its optimizations to the example input that run poisson neumann.py builds: it must be correct and efficient for any valid input as specified below, and must not rely on properties that happen to hold for that particular example. 1. Linear System Solve the Poisson-type system A x = b on a uniform cell-centered grid of shape (res x, res y, res z); every input array (a diag, a x, a y, a z, b, is dof) and the output x have this shape. Only cells marked by the boolean array is dof are active degrees of freedom. Support arbitrary coefficients and active-cell configurations within the guarantees below, start from a zero initial guess, and iterate until the solution x you return – after the zero-mean shift described below – satisfies ||b - A x|| 2 / ||b|| 2 < tol, with A x taken over the active cells. This is the true residual of the returned x, not the residual carried along by the iteration’s recurrences, which can drift from it in FP32. The system has pure Neumann boundary conditions: the normal derivative is zero at both the grid boundary and interfaces with inactive cells. These conditions are already encoded in the input coefficients. Each active row sums to zero, up to FP32 rounding: its diagonal equals the sum of the magnitudes of its retained couplings. Removing a coupling also reduces the diagonal, so a missing neighbor means zero flux rather than a prescribed zero value. A valid input also satisfies the following. The off-diagonal entries are zero or negative. A coupling to an inactive cell or to a cell beyond the grid is stored as zero, and every entry of an inactive cell – a diag, a x, a y, a z and b – is zero. Every active cell’s diagonal is positive, though it may be very small. The tolerance tol is positive. No boundary fixes the solution value, so A is singular. The active cells form one component, connected through non-zero couplings, giving a one-dimensional null space spanned by the vector that is one on active cells and zero elsewhere. The right-hand side sums to zero over active cells up to floating-point rounding; prevent the resulting null-space component from growing during iteration. Since the solution is determined only up to an additive constant, return the representative with zero mean over active cells, and set inactive cells to zero. A right-hand side of zero has the solution zero. 2. Implementation and interface Fill in the implementation in the single file poisson neumann gen.cu and build it into the Python extension poisson neumann gen.so with the provided xmake.lua. The file already defines the PoissonNeumann class, its pybind11 bindings, and the module init; the comments above the class describe the input and output format – what each array holds and the coefficient conventions – and the comments inside it name each device buffer and member. The class already implements input(), exec() and output(): input() validates the arrays, stores the scalar parameters in members, copies each input array into a device buffer with the same layout, allocates a device buffer for the output, and then calls allocate(); exec() records a CUDA start event, calls compute(), synchronizes the device, records the stop event, and returns the elapsed GPU time in milliseconds; output() copies the output device buffer into the array it is given. Do not modify input(), exec(), output(), the destructor, the error helper, the bindings, or the members input() sets, and do not change what they do indirectly –

Neumann Poisson Prompt (continued) through macros, overloads or wrappers that redefine the CUDA calls or names they use, whether in the source file or through xmake.lua. The working directory also contains run poisson neumann.py, whose run(folder, name) builds the benchmark’s input system and drives one compiled module, so you can call run(".", "poisson neumann gen") to exercise your own build and read back its timing and output. You may change the compilation flags in xmake.lua (optimization, architecture or register options, for example), but not its defines, include paths, forced includes, source files, targets or build scripts. For scoring, the extension is rebuilt from your source file and xmake.lua in a clean directory, with the fixed code above restored from the original skeleton – so keep it, and the ”===== Fixed” and ”===== Implement below” marker comments, in place; a prebuilt .so you leave behind is not used. Each timed call runs in a fresh process, after warm-up calls on other inputs of the same kind, and the output of the timed call is the one checked. You may choose your own device-memory layout and internal data structures. You implement the three private methods marked in the file, and may add members, device functions and kernels: • allocate() allocates whatever extra memory compute() needs – device memory with cudaMalloc or its variants, or pinned host memory with cudaMallocHost – sized from the grid resolution and the scalar parameters. It may compute sizes and constants from those, and must do nothing else: no kernel launches, memory copies, memsets or other GPU work, no device, stream or kernel configuration (such as cudaFuncSetAttribute, cudaFuncSetCacheConfig, cudaDeviceSetLimit or an L2 access policy) and no creation of streams, events, CUDA graphs or library handles (do both in compute()), and nothing that reads the input data. Size queries that launch nothing (such as a CUB call with a null temporary buffer) and read-only queries of device or kernel properties are allowed. Storage whose size depends on the values in the input data, rather than only on the array shapes and scalar parameters, cannot be sized here: either reserve a bound computed from the shapes and scalar parameters, or allocate it inside compute(), where the allocation is timed. • compute() performs the entire solve. It reads the input device buffers, which it may overwrite, and writes the solution x into the output device buffer. • release() frees what allocate() allocated. Do not read or write any files: nothing is loaded from disk, and the solution must be written into the output device buffer, from which output() copies it, rather than saved to x.npy or anywhere else. Report CUDA errors as Python exceptions (e.g. throw std::runtime error, which pybind11 maps to a Python exception); the file’s ck() helper does this for a CUDA call. 3. Computation and timing Only the time exec() returns is scored, and it covers compute() alone, so everything you implement that belongs to the solve must run inside compute(), on every call – including building the preconditioner, the zero-mean shift of the solution, any conversion of the input buffers into another layout, and any zero-initialization of buffers. Do not perform any part of it anywhere else: not in allocate() or release(), not in constructors, static or global initializers or the module init, not on a host thread that outlives compute(), and not by reusing results from an earlier call. </TASK> <CONSTRAINT> 1. Do not perform any web searches or access any online resources. 2. Do not run agents in parallel. At most one agent may be working at any moment, counting yourself and any subagent you start: if you delegate, the subagent has to finish and hand back before you do anything else. Do not start two subagents in one step, do not start one while another is running, and do not put one in the background. 3. You have 30 minutes of wall-clock time, ending at <deadline iso> (Unix time <deadline epoch>). Run date +%s and compare it against that number whenever you want to know how much is left. At the deadline the run is killed mid-action and whatever is on disk is what gets evaluated – so get a working build in place early, and treat anything after that as optional improvement you can afford to lose. </CONSTRAINT>

A.2

I NITIAL CUDA C ODE

The agent implements allocate(), compute(), and release() while preserving the fixed interface. The fixed exec() method times compute() and synchronizes the device before recording the stop event. 1 2 3

#include <cuda runtime.h> #include <pybind11/numpy.h> #include <pybind11/pybind11.h>

4 5 6

#include <stdexcept> #include <string>

7 8

namespace py = pybind11;

9 10 11

namespace gen impl {

12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

// Solves a Poisson−type system A x = b on a uniform cell−centered grid of // shape (res x, res y, res z), with a preconditioned Conjugate Gradient // iteration starting from a zero initial guess. Only a subset of the cells // are active degrees of freedom; the rest carry no equation. // // −−− input format −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− // Every array is 3D of shape (res x, res y, res z), where [i, j, k] is the // cell at spatial index (x, y, z). The resolution is inferred from the // shapes; all six arrays must agree. a diag, a x, a y, a z and b are // float32; is dof is bool. // // a diag[i,j,k] the diagonal entry of A at cell (i, j, k) // a x[i,j,k] the matrix entry coupling (i, j, k) to (i+1, j, k) // a y[i,j,k] the matrix entry coupling (i, j, k) to (i, j+1, k) // a z[i,j,k] the matrix entry coupling (i, j, k) to (i, j, k+1) // b[i,j,k] the right−hand side // is dof[i,j,k] true when the cell is an active degree of freedom // // A is symmetric, so the entry coupling (i, j, k) to (i−1, j, k) is // a x[i−1, j, k], and likewise along y and z. A coupling whose neighbour // would fall outside the grid is absent from the system: its entry, at the // last index along that axis, is present in the array but zero. // // Sign convention: a x, a y and a z hold the actual signed matrix entries, // not positive coupling magnitudes. For the standard discrete Poisson // operator the off−diagonals are negative (e.g. −1 on a unit grid for an // interior face) and the diagonal is positive, so // // (A x)[I] = a diag[I] ∗ x[I] + sum over neighbours of a off ∗ x[nbr] // // with a plus sign in front of the off−diagonal sum. // // Cells where is dof is false hold no unknown: their row is not part of the // system, and they contribute nothing to an active cell’s equation. // // −−− pure Neumann boundary conditions −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− // The normal derivative is zero on every boundary, both at the edge of the // grid and at the interface with inactive cells. That makes every active // row sum to zero, up to FP32 rounding: a cell’s diagonal is the sum of the // magnitudes of the couplings it keeps, so on a unit grid a cell with six // active neighbours has a diag = 6 while a cell with three has a diag = 3. // Where a coupling is dropped the diagonal drops with it −− a dropped // coupling means no flux across that face, not a known value beyond it. // // Since no boundary prescribes a value, A is singular: it annihilates any // function that is constant on the active set. The active cells form a // single connected component, so that null space is exactly // one−dimensional, spanned by the vector that is 1 on active cells and 0 // elsewhere. Two things follow: // // − b is compatible: it sums to zero over the active cells, so a solution

64 65 66 67 68 69 70 71 72 73 74 75

// exists. It is only compatible to within floating−point rounding, and // the resulting component along the null space must not be allowed to // grow as the iteration proceeds. // − x is determined only up to an additive constant, which is why the // output below is pinned to a particular representative. // // −−− output format −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− // One float32 array of the same shape (res x, res y, res z) holding the // solution x, normalized to have zero mean over the active cells. The value // at every cell where is dof is false must be set to zero. Nothing is // written to disk: the caller supplies the destination array and output() // fills it.

76 77 78 79 80 81 82

// Throws a Python exception for a failed CUDA call. inline void ck(cudaError t e, const char ∗what) { if (e != cudaSuccess) throw std::runtime error(std::string(what) + ": " + cudaGetErrorString(e)); }

83 84 85 86 87 88

class PoissonNeumann { public: // ===== Fixed: do not modify input(), exec(), output(), the destructor, // ck(), the bindings or the members input() sets. =====

89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123

// (1) input: validate the arrays, keep the scalars, copy each input // array into a device buffer of the same layout, allocate a device // buffer for each output, then call allocate(). void input(py::array t<float, py::array::c style> a diag, py::array t<float, py::array::c style> a x, py::array t<float, py::array::c style> a y, py::array t<float, py::array::c style> a z, py::array t<float, py::array::c style> b, py::array t<bool, py::array::c style> is dof, float tol) { auto bd = a diag.request(), bx = a x.request(), by = a y.request(), bz = a z.request(), bb = b.request(), bm = is dof.request(); if (bd.ndim != 3) throw std::runtime error("a diag must have shape (res x, res y, res z)"); auto same = [&](const py::buffer info &o) { return o.ndim == 3 && o.shape[0] == bd.shape[0] && o.shape[1] == bd.shape[1] && o.shape[2] == bd.shape[2]; }; if (!same(bx)) throw std::runtime error("a x must have shape (res x, res y, res z)"); if (!same(by)) throw std::runtime error("a y must have shape (res x, res y, res z)"); if (!same(bz)) throw std::runtime error("a z must have shape (res x, res y, res z)"); if (!same(bb)) throw std::runtime error("b must have shape (res x, res y, res z)"); if (!same(bm)) throw std::runtime error("is dof must have shape (res x, res y, res z)"); res x = static cast<int>(bd.shape[0]); res y = static cast<int>(bd.shape[1]); res z = static cast<int>(bd.shape[2]); num cells = static cast<size t>(res x ) ∗ res y ∗ res z ; tol = tol;

124 125 126 127 128 129 130 131 132

const size t fbytes = sizeof(float) ∗ num cells ; const size t mbytes = sizeof(bool) ∗ num cells ; ck(cudaMalloc(&d a diag , fbytes), "cudaMalloc a diag"); ck(cudaMalloc(&d a x , fbytes), "cudaMalloc a x"); ck(cudaMalloc(&d a y , fbytes), "cudaMalloc a y"); ck(cudaMalloc(&d a z , fbytes), "cudaMalloc a z"); ck(cudaMalloc(&d b , fbytes), "cudaMalloc b"); ck(cudaMalloc(&d is dof , mbytes), "cudaMalloc is dof");

133 134 135 136 137 138 139

ck(cudaMalloc(&d x , fbytes), "cudaMalloc x"); ck(cudaMemcpy(d a diag , bd.ptr, fbytes, cudaMemcpyHostToDevice), "H2D a diag"); ck(cudaMemcpy(d a x , bx.ptr, fbytes, cudaMemcpyHostToDevice), "H2D a x"); ck(cudaMemcpy(d a y , by.ptr, fbytes, cudaMemcpyHostToDevice), "H2D a y"); ck(cudaMemcpy(d a z , bz.ptr, fbytes, cudaMemcpyHostToDevice), "H2D a z"); ck(cudaMemcpy(d b , bb.ptr, fbytes, cudaMemcpyHostToDevice), "H2D b"); ck(cudaMemcpy(d is dof , bm.ptr, mbytes, cudaMemcpyHostToDevice), "H2D is dof");

140

allocate(); ck(cudaDeviceSynchronize(), "input");

141 142 143

}

144 145 146 147 148 149 150 151 152 153

// (2) exec: run compute() between two CUDA events and return the GPU // time in milliseconds. The device is synchronized before the stop // event, so every piece of GPU work compute() issues is timed. float exec() { cudaEvent t start, stop; ck(cudaEventCreate(&start), "cudaEventCreate"); ck(cudaEventCreate(&stop), "cudaEventCreate"); ck(cudaEventRecord(start), "cudaEventRecord");

154

compute();

155 156

ck(cudaDeviceSynchronize(), "compute"); ck(cudaEventRecord(stop), "cudaEventRecord"); ck(cudaEventSynchronize(stop), "cudaEventSynchronize"); ck(cudaGetLastError(), "compute");

157 158 159 160 161

float elapsed ms = 0.0f; ck(cudaEventElapsedTime(&elapsed ms, start, stop), "cudaEventElapsedTime"); cudaEventDestroy(start); cudaEventDestroy(stop); return elapsed ms;

162 163 164 165 166 167

}

168 169 170 171 172 173 174 175 176 177 178 179

// (3) output: copy the output device buffer into the provided numpy // array, of shape (res x, res y, res z). void output(py::array t<float, py::array::c style> x) { auto bx = x.request(); if (bx.ndim != 3 || bx.shape[0] != res x || bx.shape[1] != res y || bx.shape[2] != res z ) throw std::runtime error("x must have shape (res x, res y, res z)"); ck(cudaMemcpy(bx.ptr, d x , sizeof(float) ∗ num cells , cudaMemcpyDeviceToHost), "D2H x"); }

180 181 182 183 184 185 186 187 188 189 190 191

~PoissonNeumann() { release(); cudaFree(d a diag ); cudaFree(d a x ); cudaFree(d a y ); cudaFree(d a z ); cudaFree(d b ); cudaFree(d is dof ); cudaFree(d x ); }

192 193 194 195 196 197

private: // −−−−− Set by input(): read them, do not reassign them. −−−−− int res x = 0, res y = 0, res z = 0; size t num cells = 0; // res x ∗ res y ∗ res z float tol = 0.0f;

198

// Input device buffers, laid out exactly as a diag, a x, a y, a z and // b: float32, (res x, res y, res z), row−major. compute() may // overwrite them. float ∗d a diag = nullptr; float ∗d a x = nullptr; float ∗d a y = nullptr; float ∗d a z = nullptr; float ∗d b = nullptr;

199 200 201 202 203 204 205 206 207

// Input device buffer, laid out exactly as is dof: bool (one byte per // cell), (res x, res y, res z), row−major. compute() may overwrite it. bool ∗d is dof = nullptr;

208 209 210 211

// Output device buffer, laid out exactly as x: float32, // (res x, res y, res z), row−major. compute() writes the result here, // zero at every non−DoF cell. float ∗d x = nullptr;

212 213 214 215 216

// ===== Implement below. You may add members, device functions and // kernels. =====

217 218 219

// allocate(): called once, at the end of input(). Allocate the extra // memory compute() needs (device memory with cudaMalloc or its // variants, pinned host memory with cudaMallocHost), sized from // res x , res y , res z and the scalar parameters. Allocation only: // no kernel launches, copies, memsets or other GPU work, no device, // stream or kernel configuration, no creation of streams, events, CUDA // graphs or library handles, and nothing that reads the input data. // Size queries that launch nothing and read−only property queries are // fine. Storage sized by the input data belongs in compute() (timed), // or reserve a bound here. void allocate() {}

220 221 222 223 224 225 226 227 228 229 230 231

// compute(): the whole solve, run between exec()’s CUDA events. Read // the input device buffers, write the output device buffer. All of // the algorithm’s work happens here, on every call −− the zero−mean // shift included. void compute() {}

232 233 234 235 236 237

// release(): free what allocate() allocated. Called by the destructor. void release() {}

238 239

};

240 241 242

} // namespace gen impl

243 244 245 246

PYBIND11 MODULE(poisson neumann gen, m) { m.doc() = "pybind11 + CUDA: preconditioned CG solver for a Poisson−type system";

247

py::class <gen impl::PoissonNeumann>(m, "PoissonNeumann") .def(py::init<>()) .def("input", &gen impl::PoissonNeumann::input, py::arg("a diag"), py::arg("a x"), py::arg("a y"), py::arg("a z"), py::arg("b"), py::arg("is dof"), py::arg("tol"), "Input the system arrays a diag, a x, a y, a z, the right−hand " "side b, the DoF mask is dof and the relative−residual tolerance," " copy the arrays to the device, keep the scalars, and allocate.") .def("exec", &gen impl::PoissonNeumann::exec, "Run the preconditioned CG solve on the GPU, leaving the solution " "normalized to zero mean over the active cells, and return its " "GPU time in milliseconds (measured with CUDA events around " "compute()).") .def("output", &gen impl::PoissonNeumann::output, py::arg("x"), "Output the solution back into a numpy array, zero at non−DoF " "cells.");

248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265

}

Listing 1: Unmodified CUDA starter code for the Neumann Poisson task.

A.3

P YTHON D RIVER

The agent’s working directory also contains the public driver run poisson neumann.py. It builds the benchmark system with Numba on a 2563 cell-centered grid: a solid cylinder along the z-axis, of radius 0.18 of the grid width, is embedded in the box, the active cells are those containing any fluid, each off-diagonal coefficient is the negated fluid fraction of the face between two active cells, and each diagonal is the sum of these fractions, so every active row sums to zero. The righthand side is b = Ax⋆ for a manufactured field x⋆ (a ramp plus seeded uniform noise, shifted to zero mean over the active cells), and the solver tolerance is 10−6 . Its run(folder, name) loads the compiled module and, on a fresh solver object each time, calls input() and exec(); it returns the fastest of five timed calls after three discarded warm-up calls, together with the last call’s solution, read back by output(). Agents use it to test and time their builds. For scoring, the harness uses the same driver code together with a hidden driver, which the agent cannot read. The hidden driver re-executes the public one with a different seed and a cylinder radius of 0.21, so the grid size and tolerance are unchanged but the active set, the cut-cell coefficients, and b differ. Each timed call runs in a fresh process after three warm-up calls on the other input set (Section 4.1). A submission is correct on an input set if its true relative residual ∥b − Ax∥2 /∥b∥2 is below the tolerance, the magnitude of its mean over the active cells is below 10−4 of its RMS, and it is exactly zero on inactive cells.

AUDITOR

Deterministic checks and controlled evaluation. The audit combines static checks with two LLM review sessions. Before compilation, static checks compare the protected CUDA interface against the starter code, ignoring comments and whitespace, and inspect the submitted build configuration. Changes to protected methods, members, bindings, or the error helper, attempts to redefine names used by the fixed code, and build changes beyond permitted compilation flags receive HARNESS MODIFIED. A separate check inspects allocate() and identified helper calls for prohibited operations, including kernel launches, copies, initialization, device configuration, stream or library-handle creation, and host-thread creation. Detected violations receive WORK OUTSIDE COMPUTE; allocation and recognized non-executing size queries are allowed. For evaluation, the harness restores the fixed code and rebuilds the submission in a clean directory, so the scored module uses the prescribed input, timing, and output interface. LLM audit configuration. We use Claude-Opus-5 through Claude Code for every evaluated system, with shared prompts and two separate read-only sessions. The transcript auditor receives only agent.log and checks online access and parallel agent execution. A deterministic scan of tool names omits the parallel-agent question when it detects no calls that could start an agent. The timer auditor receives only the submitted CUDA source and traces whether any algorithmic work occurs outside compute(), including work in constructors, initializers, cleanup, asynchronous host threads, or state reused across calls. The fixed exec() times compute(), so this check concerns the location and dependencies of the computation rather than rechecking the event scaffold. Initialization, input-dependent preparation, and layout conversion must all occur inside compute() on every call. Host synchronization, allocation, and device-buffer copies inside it are allowed. This division assigns explicit interface and API restrictions to deterministic checks and uses the LLM for contextual interpretation of execution traces and computation dependencies. A common auditor keeps the review procedure consistent across systems. Findings must cite log events or source lines; the auditor cannot access reference solutions, benchmark scores, or the other audit session. Audit explanations and transcripts are retained for inspection. Numerical correctness is evaluated separately. A run passes the audit only when all four checks—static, network, parallelagent, and timing—return passing verdicts; missing verdicts are not treated as passes. The following boxes reproduce the transcript and timer audit prompts from the evaluation harness. The transcript prompt shown includes both checks; when the tool-name scan detects no possible delegation, the harness omits Question 2 and requests only the network verdict. Transcript Auditor Prompt <TASK> The file agent.log in the current directory is a JSON event transcript (one JSON object per line) of an agent that solved a coding task under constraints. Report only what is present in the file: do not speculate about what the agent might have done, or about what its harness may do internally. QUESTION 1 – the network. The constraint was no web or network access of any kind. Determine whether the agent actually reached the network. Inspect every tool call and shell command it ran for anything that reaches out: curl, wget, pip/conda/apt install, git clone/fetch/pull, ssh, scp, nc, or python/node code using urllib, requests, httpx, socket or fetch. Also check whether any built-in web tool was used (names like WebSearch, WebFetch, or web search). Quote the exact command or event for anything you find, with its line number in agent.log. A file path that merely contains a URL-like string is not network access. QUESTION 2 – parallel agents. The constraint was at most one agent working at any moment, counting the main agent and any subagent it started. Delegating was allowed; overlapping was not. Find every tool call that starts an agent – the name varies by harness, and includes Task, agent, subagent, subagent fork, task, workflow and ralph, as well as anything whose arguments describe dispatching work to another agent – and judge whether any two of them overlapped. What overlap looks like in a transcript: - two or more agent-starting calls in the same assistant message or step, which is how a harness fans out in parallel - a second agent started before the first one’s result appeared in the transcript

- a call made in a background or asynchronous mode, whose result is collected later by a separate poll, wait or stop call - a single call that itself fans out, such as a workflow or a batch tool given a list of tasks to run at once One subagent at a time, each finishing before the next begins, satisfies the constraint. Parallel shell commands, background shell jobs and concurrent file reads are not agents and are out of scope. Quote the events you rely on, with their line numbers, and say which two agents you believe overlapped. End your reply with exactly two lines, in this order and nothing after them – one of each pair: VERDICT NETWORK: NO NETWORK ACCESS VERDICT NETWORK: NETWORK ACCESS FOUND VERDICT AGENTS: NO PARALLEL AGENTS VERDICT AGENTS: PARALLEL AGENTS FOUND </TASK> Timer Auditor Prompt <TASK> Read the CUDA source file in the current directory (* gen.cu). It defines one class whose input(), exec() and output() are fixed: input() uploads the arrays into device buffers and calls allocate(); exec() records a CUDA start event, calls compute(), synchronizes the device, records a stop event and returns the elapsed time; output() copies the results back. Only the number exec() returns is scored, so all of the algorithm must run inside compute(). The author wrote allocate(), compute() and release(), and anything else they added. Whether the fixed code is intact and whether allocate() only allocates have already been checked mechanically; do not repeat those checks. Answer one question: does any part of the algorithm run outside compute()? Look at everything the author added: constructors, static or global initializers, the module init, release(), helper functions and whoever calls them, host threads that keep running after compute() returns, and state kept from an earlier call – static or global variables, including device buffers, that let a later call skip work (for example a result cached against a hash or sample of the input). Trace where a suspicious result is computed and where it is consumed. Work inside compute() is allowed however it is organised, including host synchronisation, allocation, and copies between device buffers. Quote the relevant lines with line numbers for anything you find. Report only what the source shows. End your reply with exactly one of these two lines and nothing after it: VERDICT: TIMER COVERS ALL WORK VERDICT: WORK OUTSIDE COMPUTE </TASK>

Audit failure types. In the reported evaluation, all audit failures come from deterministic checks and receive HARNESS MODIFIED: submissions change or omit protected interface code or modify the build configuration beyond permitted compilation flags. The LLM audits report no timing, network-access, or parallel-agent violations. Example timing violation. In an earlier development run of GPT on poisson neumann, input() constructs a multigrid hierarchy on the CPU, and the timed solver reuses the uploaded coarse operators. This input-dependent computation is part of the solve and must be included in the measured time. In the current interface, input() is fixed, and the task prompt and starter code (Appendix A) explicitly require preconditioner construction and all other algorithmic work to run in compute(), within the timed exec() region.

D ETAILED R ESULTS

This section reports the per-task outcomes behind Table 1, their variation across repeated runs, scores that weight task families equally, the development cost of each system, and the margins by which the evaluated submissions pass or fail the numerical checks. C.1

P ER -TASK R ESULTS

gen Table 3 lists, for every task and system, the speedup si = tref of each passing submission on i /ti the hidden inputs, grouped by the categories of Figure 3.

All six systems pass 21 of the 50 tasks. Pass/fail differences across systems therefore arise mainly from the remaining tasks. The best speedup reaches at least 0.95 on nine tasks: the streaming, collision, and complete LBM steps, both surface-tension operators, and contact_dem. Most of these are single memory-bound passes that the reference already runs close to peak memory bandwidth, leaving little room to improve on it. The few submissions that exceed the reference all come from these tasks and do so by at most 0.2%. At the other end, no system reaches 0.5 on 16 tasks. These include both Poisson solves, the implicit viscosity solve, the Newton FEM step and the projective-dynamics steps, the semi-implicit MPM step, the XPBD tasks, and three of the four CCD queries. Most of these tasks involve an iterative solver, irregular contact or constraint processing, or a candidate search whose cost depends on the data. C.2

RUN - TO -RUN VARIATION

Table 4 reports the three independent default-budget runs of Opus and GPT used for the best-ofthree ablation (Section 4.4). Single-run results vary little. Over the three runs, Opus reaches fast0.5 of 64%–66% and GPT 38%–42%, and the geometric mean speedup over passing tasks is 0.52–0.53 for Opus and 0.29–0.33 for GPT. The difference between the two systems is therefore much larger than the variation within either. The pass rate varies more for GPT, whose runs pass 50, 47, and 49 tasks. Table 2 forms the best of three by hidden-input speedup, which an agent cannot observe. Among runs that pass all numerical checks and audits, selecting each task’s run by its public-input speedup instead gives the same fast0.5 , fast0.9 , and fast1.05 for both systems, and a geometric mean that differs by less than 0.002. Both selection rules are post-hoc comparisons: the agent has neither the full evaluation verdicts used to filter runs nor the reference timings used to compute speedups. Only one run passes the checks on one input set but not the other: the second GPT run on self_ collision_xpbd passes on the public inputs (0.88 of the tolerance) but fails on the hidden ones (1.38). Counting it as passing would not change any fastp , since its speedup is 0.11. The only task on which any run is more than 5% faster than the reference is dem, where the second Opus run reaches 1.13×. C.3

FAMILY-BALANCED S CORES

Some tasks are near-variants of one computation: the five advection schemes, the six LBM streaming, collision, and full-step tasks on two lattices, the four CCD primitive pairs, and the two Poisson solves. Table 5 compares the task-weighted metrics of Table 1 with scores in which each such family counts once. The family-weighted fast0.5 preserves the order of the systems, with Opus at 65% and GPT at 35%. The family-weighted fast0.9 is 8% for every system other than Opus, because their fast0.9 successes lie mostly in the LBM family. C.4

D EVELOPMENT T IME AND T OKEN U SAGE

Table 6 summarizes how each system uses its budget in the runs of Table 1. Opus, GPT, Gemini, and DeepSeek finish every task before the deadline, using a median of 25%–54% of the budget. GLM reaches the deadline on 4 tasks. Qwen reaches it on 19 tasks and uses a median of 92% of the budget. Token usage differs widely across systems: DeepSeek consumes the most input and output tokens, while Gemini produces the fewest output tokens and finishes fastest. 23

Table 3: Per-task results. Each entry is the hidden-input speedup si of a passing submission relative to the reference (higher is better; 1.00 matches the reference). “–” marks a submission that fails the numerical checks or produces no valid output, and “A” marks a numerically correct submission that fails an audit. Tasks marked ∗ are checked by the residual of the returned solution rather than by comparison with the reference. Task

Opus

GPT

Gemini

DeepSeek

Qwen

GLM

Local Grid and Lattice Computations advect_rk1 advect_rk2 advect_rk3 advect_rk4 advect_mc stream_d3q19 stream_d3q27 trt_collision_d3q19 trt_collision_d3q27 lbm_d3q19 lbm_d3q27 cahn_hilliard surface_tension_sdf surface_tension_phase_field

0.87 0.83 0.85 0.89 0.81 1.00 0.96 1.00 0.99 0.99 1.00 0.82 0.99 1.00

0.80 0.74 0.44 0.76 0.51 0.96 1.00 0.99 0.99 0.99 1.00 0.68 0.99 0.98

0.55 0.70 0.42 0.75 0.66 0.97 0.95 0.99 0.98 0.98 0.99 0.34 0.99 1.00

0.88 0.83 0.73 0.44 0.64 0.99 1.00 0.99 0.99 0.99 1.00 – 1.00 0.99

– 0.67 0.62 0.86 – – 1.00 1.00 0.99 0.99 1.00 0.84 0.98 1.00

– 0.82 – 0.63 0.54 1.00 1.00 0.99 0.99 0.99 1.00 0.64 0.95 1.00

Particle–Grid Methods p2g_pic_flip g2p_pic_flip p2g_apic g2p_apic p2g_mpm g2p_mpm velocity_gradient_mpm apic pic_flip mpm_explicit mpm_semi_implicit∗

0.56 0.76 0.77 0.66 0.92 0.37 0.46 0.51 0.44 0.68 0.37

0.56 0.76 0.10 0.76 0.69 0.37 0.40 0.21 0.16 0.36 0.05

0.13 0.11 0.12 0.67 0.12 0.16 0.16 – – – 0.03

0.13 0.37 – 0.62 0.54 – 0.46 0.17 – 0.26 –

0.15 0.62 0.08 0.65 0.10 0.14 0.24 – – 0.44 –

0.37 0.69 0.30 0.69 0.10 – 0.27 0.09 0.13 0.14 0.10

Local Interaction Updates mass_spring_explicit fem_explicit contact_dem dem viscosity_pbf vorticity_confinement_pbf

0.62 0.93 1.00 0.86 0.68 0.78

0.35 0.57 0.36 0.14 0.16 0.40

0.36 0.34 0.13 0.10 0.10 0.09

0.62 0.53 0.82 0.45 0.54 0.51

0.62 0.55 0.78 – – –

0.62 0.55 0.15 0.25 0.11 A

Position Constraint Solving solve_constraint_pbf pbf solve_constraint_xpbd self_collision_xpbd xpbd

0.82 0.57 0.27 0.41 0.32

0.56 0.62 0.11 0.08 0.22

0.10 – 0.20 – 0.20

0.23 0.02 0.28 0.18 0.26

– – 0.27 – 0.32

0.06 0.51 0.28 0.12 0.27

Geometric Queries and Collision Detection ccd_vv ccd_ve ccd_vf ccd_ee particle_sdf

0.21 0.40 0.38 0.22 0.54

0.01 0.10 0.16 0.04 0.57

0.03 0.11 0.53 0.12 0.26

0.06 0.05 – – 0.47

– – – – –

0.06 0.41 0.22 0.06 0.47

Global Solves and Pressure Projection poisson_dirichlet∗ poisson_neumann∗ viscosity_implicit mass_spring_newtonian_implicit∗ fem_newtonian_implicit∗ mass_spring_pd fem_pd stable_fluids mc_r

0.05 0.04 0.10 0.58 0.30 0.43 0.32 0.54 0.66

0.24 0.25 0.11 0.29 0.21 0.40 0.29 0.19 0.34

0.04 0.04 – 0.25 0.12 0.18 0.20 0.03 0.03

0.30 0.14 0.29 0.39 0.25 0.34 0.31 0.32 0.10

– – – – – 0.27 – – 0.52

0.25 – – 0.26 0.17 0.31 0.30 – 0.06

Table 4: Run-to-run variation at the default time limit. Runs 1–3 are independent single attempts; run 1 is the run in Table 1. Metrics are percentages over all 50 tasks; the geometric mean is the speedup over passing tasks. Best of 3 keeps, per task, one run that passes all numerical checks and audits: public selection picks the highest public-input speedup, and hidden selection picks the highest hidden-input speedup, as in Table 2. Both report the hidden-input speedup of the chosen run. These post-hoc selections require evaluator verdicts and reference timings unavailable to the agent during development. Model

Run

Pass Rate

fast0.5

fast0.9

fast1.05

Geo. mean

Claude-Opus-5 Run 1 100% 66% 22% 0% 0.53 Run 2 98% 64% 24% 2% 0.53 Run 3 100% 64% 28% 0% 0.52 Mean ± s.d. 99.3 ± 1.2 64.7 ± 1.2 24.7 ± 3.1 0.7 ± 1.2 0.53 ± 0.01 Best of 3, public selection 100% 76% 36% 2% 0.63 Best of 3, hidden selection 100% 76% 36% 2% 0.63 GPT-5.6-Sol

Run 1 100% 42% 16% 0% 0.33 Run 2 94% 38% 18% 0% 0.30 Run 3 98% 38% 16% 0% 0.29 Mean ± s.d. 97.3 ± 3.1 39.3 ± 2.3 16.7 ± 1.2 0.0 ± 0.0 0.31 ± 0.02 Best of 3, public selection 100% 48% 18% 0% 0.42 Best of 3, hidden selection 100% 48% 18% 0% 0.42

Table 5: Task-weighted / family-weighted metrics. The family-weighted score groups near-variant tasks into one family each: the five advection schemes, the six LBM streaming, collision and fullstep tasks, the four CCD primitive pairs, and the two Poisson solves. Every other task is its own family, giving 37 families of equal weight; within a family, tasks are weighted equally. Model

Pass Rate

Claude-Opus-5 100% / 100% GPT-5.6-Sol 100% / 100% Gemini-3.5-Flash 88% / 84% DeepSeek-V4.1-Flash 86% / 85% Qwen-3.8-Max 52% / 53% GLM-5.3 86% / 87%

C.5

fast0.5

fast0.9

66% / 65% 22% / 16% 42% / 35% 16% / 8% 28% / 14% 16% / 8% 38% / 29% 16% / 8% 34% / 28% 14% / 8% 34% / 26% 16% / 8%

C ORRECTNESS M ARGINS

Every numerical check records its error and its tolerance. For each submission, we take the worst error-to-tolerance ratio over the task’s checks on both the public and the hidden inputs; a ratio below 1 passes. Figure 5 shows the distribution of this ratio over the submissions that produce an output. It separates comparison-based checks from residual-based checks, whose ratio has a different meaning. Comparison-based checks. The 234 numerically correct submissions all lie well below the tolerance. Their median ratio is 0.027, 76% are below 0.1, and the largest is 0.60. The 27 failing submissions all lie above it. Twenty-five exceed their tolerance by more than 100×. The two closest failures still exceed it by 2.2× and 7.6×, and both are genuine deviations from the specification (Appendix C.6). No submission falls between 0.6 and 2.2 times its tolerance. The observed outcomes would therefore be unchanged by any threshold in that range. Residual-based checks. For the Poisson solves, the Newton steps and the semi-implicit MPM step, the ratio records where the submission’s own iteration stopped relative to the required residual. For the semi-implicit MPM step, the worst ratio of every passing submission (0.60) instead comes from a consistency check between the returned particle and grid states. Iterative solvers terminate once they meet the criterion, so passing submissions cluster just below 1, with a median of 0.73. This proximity reflects the stopping rule rather than a narrow margin. Any convergent solver can 25

Table 6: Agent development cost per task in the runs of Table 1, as medians over the 50 tasks. Time is the wall-clock development time, also given as a fraction of the task budget; deadline hits count runs stopped at the budget. Input tokens include cached prompt tokens, whose share over all tasks is given separately. Token counts are reported by each harness; for runs stopped at the deadline, some harnesses report only the usage recorded before termination. Model

Harness

Claude-Opus-5 GPT-5.6-Sol Gemini-3.5-Flash DeepSeek-V4.1-Flash Qwen-3.8-Max GLM-5.3

Claude Code Codex CLI Gemini CLI DeepSeek Harness Qwen Code OpenCode

Time (min) Budget used Deadline hits Input tok. Cached Output tok. 15.7 18.0 8.9 18.5 30.0 23.8

43% 50% 25% 54% 92% 66%

0 0 0 0 19 4

1.65M 2.72M 2.88M 7.34M 1.65M 2.96M

96% 97% 90% 99% 49% 97%

45k 33k 16k 116k 41k 66k

Comparison-based checks Passing (n=234) Failing (n=27)

no submission: 0.6–2.2 50 48 38 28 19

15 1

Residual-based checks 21

≤10

−5

−3

Passing (n=23) Failing (n=5)

no submission: 0.99–5 × 104

1 −1

1 5

10 10 10 10 10 10 Worst error / tolerance over the checks on both input sets (vertical line: tolerance)

1 7

Figure 5: Distribution of the worst error-to-tolerance ratio over the checks on both input sets, for every evaluated submission that produces an output, in half-decade bins. The vertical line marks the tolerance; the shaded band is the range between the largest passing and the smallest failing ratio, which contains no submission. Comparison checks measure the relative ℓ2 difference from the reference output; residual checks measure the true residual of a returned linear or Newton solution against the solver tolerance, and passing solvers stop just below it by design. Values are clipped to [10−5 , 108.5 ]; the rightmost residual failure is non-finite.

move further from the threshold by iterating longer, at a cost in speed. The five failing submissions exceed the required residual by at least 5 × 104 . Public and hidden inputs. No submission passes the checks on one input set and fails them on the other. The hidden inputs therefore rejected no submission that had passed on the public inputs in these runs. They remain a safeguard against implementations tailored to the public example. Failures without output. Eleven failing submissions produce no output to check, and they are not shown in Figure 5. Eight cannot be rebuilt from the submitted source: six do not compile, and two have removed the skeleton’s fixed sections, which the evaluator needs to restore. One fails at run time. The remaining two produce non-finite values, which fail the checks before any tolerance is applied. 26

C.6

VALIDATION OF THE T OLERANCES

The margins above describe the submissions that happen to have been evaluated. Three further experiments test the tolerances directly. The first measures how far a known-correct variant of each reference moves. The second measures how far plausible bugs move. The third inspects the submissions nearest the threshold. All three use the evaluator’s own checks on the public and hidden inputs, and report the worst error-to-tolerance ratio as in Appendix C.5. Known-correct variants. We rebuilt every reference with fused multiply-add contraction disabled (-fmad=false), changing nothing else. This changes the rounding of nearly every arithmetic expression, a choice a correct submission may equally make. All 45 comparison-based tasks pass, with a median ratio of 0.012 and a maximum of 0.18. The five residual-based tasks score 0.30–0.62, which again records where the solver stopped (or, for the semi-implicit MPM step, a consistency check between the returned particle and grid states) rather than any difference from the reference. Plausible implementation errors. We selected nine tasks spanning the categories of Figure 3. For each, we wrote four mutants of the reference, each a single small edit modeling a mistake an implementer could plausibly make from the specification. The mutants fall into five classes: omitting a term, flipping a sign, mishandling a boundary, reading a stage out of order or updating in place, and using a wrong coefficient or stopping criterion. Each task receives mutants from four of the five classes; the fifth is omitted where it has no meaningful single-edit form, such as a boundary rule in a purely per-cell collision step. • Detected. All 36 mutants fail. Their median ratio is about 1.6 × 103 , and 26 exceed the tolerance by more than 100×. • Closest to the threshold. Four mutants come within 10× of it: – In mc_r, each pressure solve stops at a relative residual of 10−2 instead of 10−4 . This mutant scores 1.07 on the public inputs and 3.4 on the hidden ones. – In pic_flip, gravity is applied with the wrong sign. It scores 1.7, because a single step’s gravity increment is small next to the pressure projection of the initial velocity field. – In ccd_ee, the bisection that locates each crossing stops after 15 steps instead of the specified 30. It scores 6.0. – In xpbd, positions are updated in place instead of double-buffered. It scores 8.7. Submissions nearest the threshold. We inspected the passing submissions with the largest ratios and the failing submissions with the smallest. • The two mc_r submissions at 0.56 and 0.60 are correct. Every mc_r submission scores at least 0.53 on the same hidden check. This floor is the reference’s own error: its pressure solves stop at the relative residual of 10−4 that the specification prescribes. Against a reference solved to 10−6 , the two submissions score 0.18 and 0.29. Re-solving one of them to 10−6 brings it to 0.001 of the tolerance, so FP32 rounding contributes almost nothing. On this task, the margin left to a correct solver is therefore about 1.7×, set by the solver tolerance rather than by floating-point variation. It could be widened by producing the reference output with a tighter solve. • The ccd_ee submission at 0.57 is also correct. All but one of its time-of-impact values agree with the reference to within one unit in the last place. The remaining edge crosses its partner about 10−6 from an endpoint. At that edge, the submission forms the difference between the two edges from absolute coordinates, and the reference forms it from a precomputed relative offset. The resulting cancellation rejects this crossing and returns the edge’s next impact instead. • The pic_flip submission that fails at 7.6× has a real bug. It updates the level set only in the cells within one cell of each particle, whereas the specified level set can be lowered by a particle at any cell centre within r + ∆x ≈ 1.87 cells, which spans two cells on each side. The missed cells are the air cells adjacent to the free surface, which set the free-surface fractions of the pressure system. Widening the three loops to two cells, with no other change, brings the submission to 0.16, in line with the passing submissions. 27

• The self_collision_xpbd submission that fails at 2.2× also has a real bug. Its hashed cell lookup can visit the same bucket twice, which records a few contact pairs twice, contrary to the specification that each pair appears once. Taken together, correct variation in these experiments stays below about 0.6 of the tolerance. Every planted bug exceeds it, most by orders of magnitude. The narrowest gap, on mc_r, comes from a solver tolerance inherited by the reference, not from floating-point variation.

C ASE S TUDIES OF G ENERATED I MPLEMENTATIONS

We examine two tasks that illustrate complementary optimization challenges: organizing particleto-grid accumulation, and reducing iteration costs in a global solve. The cases come from the same evaluation as Table 1 and were evaluated on the NVIDIA GeForce RTX 4090. Table 7 reports speedup as reference execution time divided by generated execution time on hidden inputs, consistent with the definition used for fastp ; all four submissions pass numerical checks and audits. The discussion combines evaluation results with inspection of final CUDA code and development traces. Table 7: Final performance of the selected implementations on hidden inputs. For passing submissions, higher tref /tgen is better, and 1× matches the reference. tref /tgen ↑

Task

Model

PIC/FLIP particle-to-grid transfer

Claude-Opus-5 Gemini-3.5-Flash

0.556× 0.134×

Dirichlet Poisson solve

Claude-Opus-5 Gemini-3.5-Flash

0.049× 0.036×

D.1

PIC/FLIP PARTICLE - TO -G RID T RANSFER : O RGANIZING ACCUMULATION

The p2g pic flip task transfers particle velocities onto a staggered MAC grid using trilinear weights. Each particle contributes mass and momentum to eight neighboring nodes for each velocity component, followed by normalization of momentum by mass. Nearby particles update overlapping grid nodes, making the organization of accumulation central to performance. Gemini sorts particle indices by grid-cell ID using CUB radix sort, then processes particles in that order. However, each thread still gathers particle data indirectly from the original arrays and performs global atomic additions for every particle-node contribution. Sorting improves spatial ordering without aggregating contributions before they reach global memory. Opus groups particles into 8 × 8 × 8 cell tiles through counting and prefix sums, then packs their positions and velocities into contiguous tile buffers. One CUDA block handles each tile, accumulating mass and momentum in shared memory over a 10 × 10 × 10 region that includes the halo. The block then merges these values into global accumulators. This reduces global atomic traffic, but each particle still performs atomic updates within the shared-memory tile. All grouping, packing, accumulation, and normalization occur inside the timed compute(). Opus and Gemini achieve speedups of 0.556× and 0.134× relative to the reference, respectively. The reference further sorts particles by cell within each tile and sums each cell’s particle contributions in registers before merging them into shared memory. Its fixed-capacity tile buckets also avoid a global counting-sort pass when capacity is sufficient, with a fallback for overflowing tiles. These differences illustrate why spatial grouping alone is insufficient: the granularity of aggregation and the cost of organizing particles both matter, even after most global atomics have been replaced by shared-memory accumulation. D.2

D IRICHLET P OISSON S OLVE : I TERATION C OST AND S OLVER C HOICE

This task solves a variable-coefficient Poisson system with homogeneous Dirichlet boundary conditions on a 2563 grid. The operator is symmetric positive definite on the active cells, and the solution must satisfy a relative residual tolerance of 10−6 . Unlike fixed-iteration updates, the total work depends on both the cost of each iteration and the convergence behavior of the chosen solver. Gemini implements conjugate gradients with diagonal Jacobi preconditioning. Each iteration uses separate kernels for the matrix-vector product, dot products, solution and residual updates, and search-direction update. Two scalar reductions are copied back to the host in a typical iteration to compute the CG coefficients, introducing repeated synchronization. Periodic convergence checks recompute the true residual before accepting the solution and restart CG if necessary. 29

Opus also uses diagonal preconditioning, applying it through symmetric scaling of the system inside compute(). This folds the preconditioner into an operator with unit diagonal. One main kernel fuses the deferred solution update, search-direction update, matrix-vector product, and dot-product accumulation; another updates the residual and accumulates its norm. Additional small reduction kernels keep the CG coefficients on the device. Host transfers are used for initialization, convergence checks every ten iterations, and true-residual verification, rather than for the coefficients of every iteration. The final implementations achieve speedups of 0.049× and 0.036× relative to the reference for Opus and Gemini, respectively. Opus reduces memory passes and host synchronization, but both implementations retain a diagonal preconditioner, whereas the expert reference uses multigridpreconditioned CG. Their large remaining gaps illustrate the limits of optimizing individual iterations without also improving the solver’s convergence behavior. Efficient global solves require numerical algorithm selection and GPU execution to be optimized together.

B ENCHMARK TASKS

E.1

TASK L IST

The benchmark covers fluid dynamics, deformable solids, and granular materials. The fluid tasks span Eulerian grid-based solvers (Zehnder et al., 2018), hybrid particle-grid methods such as particle-in-cell/fluid-implicit-particle (PIC/FLIP) and affine particle-in-cell (APIC) (Jiang et al., 2015), particle-based methods such as position-based fluids (PBF) (Macklin & Müller, 2013), and the lattice Boltzmann method (LBM) (Li et al., 2026b). They include both local operations, such as advection, particle-grid transfer, lattice collision and streaming, and phase-field evolution, and global or iterative computations, such as pressure projection and implicit viscosity. The solid and granular tasks cover mass-spring systems (Liu et al., 2013), the finite element method (FEM) (Sifakis & Barbic, 2012), the material point method (MPM) (Stomakhin et al., 2013), extended positionbased dynamics (XPBD) (Macklin et al., 2016), continuous collision detection (CCD) (Brochu et al., 2012), and the discrete element method (DEM) (Lu et al., 2022). Together, these tasks exercise explicit and implicit integration, projective and constraint-based solvers, contact handling, and coupled multi-stage simulation pipelines. Table 8 lists all 50 benchmark tasks under the six categories used in Section 4.3. Test names match the task directories in the benchmark codebase. Table 8: Complete GPUPhysBench task list. Descriptions summarize the computation specified in each task prompt. No. 1

Category Global Solves and Pressure Projection

Test name poisson_dirichlet

Global Solves and Pressure Projection

poisson_neumann

Global Solves and Pressure Projection Global Solves and Pressure Projection Global Solves and Pressure Projection Global Solves and Pressure Projection

viscosity_implicit

Global Solves and Pressure Projection

fem_pd

Global Solves and Pressure Projection

stable_fluids

Global Solves and Pressure Projection

mc_r

Position Constraint Solving Position Constraint Solving

solve_constraint_pbf

Position Constraint Solving Position Constraint Solving

solve_constraint_xpbd

4 5 6

12 13

mass_spring_newtonian_ implicit fem_newtonian_implicit mass_spring_pd

pbf

self_collision_xpbd

Description Solve a Poisson system with homogeneous Dirichlet boundaries using preconditioned CG. Solve a pure-Neumann Poisson system using preconditioned CG and return a zero-mean solution. Solve implicit viscous diffusion for velocity components on a MAC grid. Advance an implicit mass-spring step using Newton’s method. Advance an implicit tetrahedral FEM step using Newton’s method. Advance a mass-spring step using projective dynamics local-global iterations. Advance a tetrahedral FEM step using projective dynamics local-global iterations. Advance a MAC-grid fluid step with advection, implicit viscosity, and pressure projection. Advance a MAC-grid fluid step using MacCormack advection, reflection, and two pressure projections. Iteratively correct particle positions to satisfy PBF density constraints. Advance a full PBF step with prediction, density-constraint solving, velocity update, vorticity confinement, and XSPH viscosity. Solve XPBD spring constraints by iteratively correcting particle positions. Resolve cloth particle self-collisions with XPBD contact constraints.

Table 8: Complete GPUPhysBench task list (continued). No. 14

Category Position Constraint Solving

Test name xpbd