GAV: A scientific experimentation environment for PDE agents

Muheng Li; Wenlong Mou

cs.LGOct 21, 2026

GAV: A scientific experimentation environment for PDE agents

Authors: Muheng Li, Wenlong Mou

Organizations: Department of Statistical Sciences, University of Toronto · Vector Institute

Abstract

Reliable feedback is essential for coding agents to develop programs for partial differential equations (PDEs). Yet evaluating these programs on new PDE problems often requires expensive reference computations, while fixed datasets restrict the range of experiments. We introduce Generator-Agent-Verifier (GAV), a scientific experimentation environment that makes new PDE problems verifiable by construction. The Generator learns a distribution of continuous reference solutions using data from numerical simulations or experimental measurements to represent a target physical regime. Using the method of manufactured solutions, it works backward from the generated reference solutions to construct matching PDE problems without a separate numerical solve. The Agent designs experiments, investigates failures, and iteratively improves its program. The Verifier evaluates candidate programs against the reference solutions and provides numerical feedback on accuracy and computational cost. This feedback can also guide the generation of new problems that probe the current program's weaknesses. We provide a theoretical analysis of evaluation reliability and guided development. Experiments on Navier-Stokes and Vlasov-Poisson-type equations across different data sources and regimes demonstrate the effectiveness of GAV. In particular, on the numerical Navier-Stokes task, GAV reduces the mean test error of three existing PDE agents by more than an order of magnitude relative to their original setups. Code is available at https://github.com/PotatoJnny/PDE-GAV.

Figures & tables

Figure 1: GAV: a scientific experimentation environment for PDE agents. The Generator constructs PDE problems and reference solutions from scientific data. The Agent requests problems and input information, tests programs, and uses the Verifier’s feedback on error and runtime to guide further experiments and revisions. The diagram illustrates a numerical solution task.

Figure 2: Multi-turn setting. The selected program and updated memory log pass to the next round; research notes from both agents are retained regardless of the selection.

Method	Mean error (%)
COMSOL ( COMSOL AB, 2025 )	1.6352
FEniCSx ( Baratta et al., 2023 )	0.3510
Nektar++ ( Moxey et al., 2020 )	1.9042
GAV-Full (ours)	0.1046

Table 1: Numerical NS solver comparison. Mean relative velocity error (%) on the same test set; lower is better. Appendix C.2 describes the baselines.

Figure 3: NS program development with GAV-Full. Mean relative velocity error (%) on the numerical and PIV-derived NS test sets over five development rounds. Lower is better; the PIV panel uses a log scale.

Figure 4: Kinetic program development with GAV-Full. Mean transport score from Equation ( 10 ) for programs retained after rounds 1–6. Each task uses the same test set and evaluation queries across both models and all rounds. Scores are dimensionless; all three panels use logarithmic axes and lower is better. R6 is the final program. Strong-collision Astra retains its round-3 program in round 4.

Task	Model	R1 error	R6 error	Reduction (%)
ARTEMIS VP	Sol	0.03961	0.02956	25.4
ARTEMIS VP	Astra	0.01074	0.002854	73.4
Weak-collision VPFP	Sol	0.189	0.1352	28.5
Weak-collision VPFP	Astra	0.2668	0.1274	52.2
Strong-collision VPFP	Sol	$9.831\times 10^{-4}$	$9.781\times 10^{-5}$	90.1
Strong-collision VPFP	Astra	$4.999\times 10^{-6}$	$6.719\times 10^{-7}$	86.6

Table 2: First-round and final kinetic programs. Mean dimensionless transport score on each task’s test set, shared by Sol and Astra. Reduction is

100(1-\mathrm{error}_{\mathrm{R6}}/\mathrm{error}_{\mathrm{R1}})

, relative to the first retained program R1. The final program is R6 in all six runs.

Figure 5: Velocity error on a PIV-derived manufactured problem. Rounds 1 (top) and 5 (bottom) use the same generated reference at five times. Color shows $\|\mathbf{u}^{\mathrm{pred}}-\widetilde{\mathbf{u}}\|_{2}$ (m/s) on one full-window linear scale: white at zero, darker blue for larger errors.

Figure 6: Mean relative velocity error (%) on the numerical NS test set; lower is better.

Environment	Mean error (%)	Error / GAV
GAV	0.441	1.0 $\times$
w/o learning from data	1.468	3.3 $\times$
w/o MMS construction	0.988	2.2 $\times$
w/o continuous input queries	5.211	11.8 $\times$

Table 3: Mean relative velocity error (%) after three rounds on the same NS test set. Error / GAV is the ratio to GAV; lower is better.

Appendix figures & tables24 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 7: Navier–Stokes data sources. Examples of discrete solution data used to construct the environment. Color shows fluid speed at a single time. (a) Numerical simulation in a square domain with two solid obstacles (gray). (b) PIV measurements of flow behind a cylinder. The white rectangle marks the measurement region used to construct problems; gray areas indicate missing measurements. Each panel uses its own color scale; PIV speeds above $0.06$ m/s use the top color.

Setting	Numerical data	PIV data
Cosine basis functions per spatial axis	$17$
Chebyshev polynomial degree in time	$12$
Retained principal components	$256$	$32$
Batch size	$512$	$128$
Training steps	$30{,}000$	$5{,}000$
Checkpoint used (training step)	$15{,}000$	$5{,}000$

Appendix

Table 4: NS Generator settings.

Quantity	Accepted range or bound
Joint relative macroscopic change	$[0.002,0.20]$
Hermite-coefficient RMS change	$[0.005,0.35]$
Relative source change	$[10^{-4},0.80]$
Mean squared norm of added Hermite coefficients	$\geqslant 10^{-5}$
RMS change in coefficient time derivative	$\leqslant 4$
Initial relative macroscopic change	$\leqslant 2\times 10^{-6}$

Appendix

Table 5: Strong-collision VPFP perturbation acceptance settings. All listed checks must pass for an accepted child problem. Initial-state checks use the same spatial probes at

t=0

Error predictor $h$
Parameter	Value
Residual blocks	$6$
Hidden width	$512$
Embedding dimension for flow time $s$	$128$
Training steps	$5{,}000$
Batch size	$256$

Appendix

Table 6: NS guidance and perturbation settings. The numerical and PIV studies share these settings. Here

[t_{0},t_{1}]

is the physical time interval.

Study	Coding model	Rounds	Time budget (h)
Numerical and PIV-derived NS	Sol	5	4
Kinetic tasks	Sol	6	3
Kinetic tasks	Astra	6	3
NS ablations	Sol	3	2

Appendix

Table 7: Configurations of the shared multi-turn development setting. Time budgets are per agent per round and include tool calls and waiting. Each kinetic task has separate Sol and Astra runs.

Study	Use	Size	Composition
Numerical NS (Section 4.2 )	Selection, R2–R5	32	11 unguided, 11 guided, 10 perturbed
	Test	32	11 unguided, 11 guided, 10 perturbed
PIV-derived NS (Section 4.2 )	Selection per round	192	64 unguided, 64 guided, 64 perturbed
	Test	768	256 unguided, 256 guided, 256 perturbed
Existing PDE agents (Section 4.3 )	Initial development	64	Generated NS problems shared across environments
NS comparisons (Sections 4.3 and 4.4 )	Selection	24	8 unguided, 8 guided, 8 perturbed

Appendix

Table 8: Problem sets for initial development, program selection, and evaluation. Counts are per dataset. Unguided problems are sampled directly from the learned Generator. NS perturbations use unguided parents; kinetic test sets include both unguided and guided parents. Selection sets in Section 4.2 are refreshed each round; test sets are fixed. Numerical NS uses

192

selection problems in round one (

64

per source) and

32

in subsequent rounds. Selection sets for Sections 4.3 and 4.4 are fixed. The numerical NS test set is shared with COMSOL, FEniCSx, and Nektar++.

Method	Original	GAV	GAV-Full
CodePDE	11.77	5.59	8.52
AutoNumerics	13.22	9.71	8.61
PDE-SHARP	20.80	13.44	12.36

Appendix

Table 9: Total agent active time (hours).

Agent setting	Environment	Mean error (%)
Multi-turn setting	GAV-Full	0.1256
CodePDE	Original	3.1231
	GAV	0.1890
	GAV-Full	1.0555
AutoNumerics	Original	3.4970
	GAV	0.3086

Appendix

Table 10: Mean relative velocity error (%) of the final numerical NS programs from the main experiments on 96 interpolated test problems. Lower is better.

Environment	Mean error (%)
GAV	0.440
w/o learning from data	1.561
w/o MMS construction	0.864
w/o continuous input queries	5.915

Appendix

Table 11: Mean relative velocity error (%) for the environment comparison on 96 interpolated NS test problems. Lower is better.

Family	Model	S0	R1	R6	Change (%)
ARTEMIS VP	Sol	0.1044	0.007703	0.005259	-31.73
ARTEMIS VP	Astra	0.1044	$1.522\times 10^{-6}$	$8.662\times 10^{-6}$	+469.09
Weak VPFP	Sol	0.5652	0.08331	0.03437	-58.74
Weak VPFP	Astra	0.5652	0.3857	0.2211	-42.68
Strong VPFP	Sol	0.2054	0.001883	$4.233\times 10^{-4}$	-77.52
Strong VPFP	Astra	0.2054	$1.472\times 10^{-8}$	$2.176\times 10^{-9}$	-85.22

Appendix

Table 12: Mean dimensionless transport error on the independent fixed banks (six ARTEMIS cases and eight cases for each VPFP regime). The signed change is

100(\mathrm{error}_{R6}/\mathrm{error}_{R1}-1)

; positive values indicate higher error. R6 is the prespecified final round.

Round	Sol	Astra	Sol	Astra	Sol	Astra
	ARTEMIS VP		Weak VPFP		Strong VPFP
S0	0.1044	0.1044	0.5652	0.5652	0.2054	0.2054
R1	0.007703	$1.522\times 10^{-6}$	0.08331	0.3857	0.001883	$1.472\times 10^{-8}$
R2	0.007637	$1.725\times 10^{-6}$	0.03311	0.2143	0.001846	$2.346\times 10^{-9}$
R3	0.005257	$8.662\times 10^{-6}$	0.03297	0.2882	0.001845	$2.165\times 10^{-9}$
R4	0.005259	$8.662\times 10^{-6}$	0.03458	0.2882	0.001845	$2.165\times 10^{-9}$

Appendix

Table 13: Mean dimensionless transport error on the independent MMS banks: six ARTEMIS cases and eight cases for each VPFP regime. The same cases and evaluation queries are used for both models in every round. S0 is the initial program; R1–R6 are the retained development programs. All scheduled cases are included, and R6 is the prespecified final round.

Task	Model	$e_{\rho}$	$e_{\mathbf{E}}$	$e_{\mathbf{j}}$	$e_{\mathbf{P}}$	$e_{\mathbf{q}}$
ARTEMIS VP	Sol	0.002574	0.01625	0.008603	0.002552	$4.708\times 10^{-4}$
ARTEMIS VP	Astra	$3.062\times 10^{-6}$	$2.765\times 10^{-5}$	$9.302\times 10^{-6}$	$1.937\times 10^{-6}$	$8.055\times 10^{-6}$
Weak VPFP	Sol	0.001524	0.008831	0.009259	0.01128	0.1118
Weak VPFP	Astra	0.041	0.4014	0.2016	0.09325	0.3682
Strong VPFP	Sol	$2.574\times 10^{-5}$	$1.435\times 10^{-4}$	$1.896\times 10^{-4}$	$2.160\times 10^{-4}$	0.001203
Strong VPFP	Astra	$1.220\times 10^{-10}$	$1.092\times 10^{-9}$	$4.464\times 10^{-10}$	$5.534\times 10^{-10}$	$7.099\times 10^{-9}$

Appendix

Table 14: Final macroscopic errors on independent kinetic problems. Mean per-problem normalized RMS error for each field group in Equation ( 10 ), evaluated for R6 on six ARTEMIS-informed VP problems or eight problems per VPFP regime. These are dimensionless component errors, before weighting: density and electric field each have weight

1/8

, and particle flux, pressure, and heat flux each have weight

1/4

. Vector errors use the Euclidean norm and VPFP pressure uses the full tensor Frobenius norm. Lower is better.

Figure 8: Kinetic programs on independently constructed problems. Mean dimensionless transport score from Equation ( 10 ) for the initial program S0 and the programs retained after rounds 1–6. The VP set contains six MMS problems constructed from twelve unused ARTEMIS profiles; each VPFP regime has eight interpolated-reference MMS problems. Programs within each task share the same problems and hidden queries, using the frozen scoring calibration. VP is evaluated at

t=0.875

; VPFP uses space–time queries. All vertical axes are logarithmic, and lower is better. The independent sets are used for evaluation only; R6 remains the final program.

Figure 9: Percentage increase in mean error after guidance updates. Each comparison fixes the program and measures the joint relative

\ell_{2}

velocity error on the old and new problem sets.

Figure 10: Numerical NS program evolution over five rounds. Each round shows the Explorer’s and Refiner’s final submissions. Blue marks selected submissions; arrows connect each submission to its starting program. Round 5 retains the incumbent from round 4, shown at the bottom. Grid sizes are given per spatial axis.

Figure 11: ARTEMIS-informed VP program evolution with Astra. Each row shows the Explorer’s and Refiner’s final nominations; blue identifies the program selected on that round’s 48 problems, and arrows identify recorded starting programs. DG denotes discontinuous Galerkin, Gauss

k

denotes

k

quadrature nodes per velocity cell, and ETDRK4 is fourth-order exponential Runge–Kutta. Source-time counts are refinement targets; an earlier complete panel may be returned. R6 is the final program; common-test errors are shown in Figure 4 .

Figure 12: Weak-collision VPFP program evolution with Astra. The two columns contain the final nominations from each round, with blue marking the program retained on the common 48-problem selection bank. Arrows show recorded starting programs. OU denotes Ornstein–Uhlenbeck and RK4 denotes fourth-order Runge–Kutta. Grid dimensions count sampling nodes. The programs predict the required macroscopic quantities using different internal approximations. Figure 4 reports their errors on the separate test bank.

Figure 13: Strong-collision VPFP program evolution with Astra. Final Explorer and Refiner nominations over six rounds. Blue marks the program retained on each round’s 48 selection problems; arrows indicate recorded starting programs. Round 4 retains the round-3 Refiner program, shown by the bypass arrow and intervening box. Both round-5 candidates therefore start from R3. Grid sizes count nodes per spatial axis, Hermite degree specifies the velocity expansion, and RK4 is fourth-order Runge–Kutta. Test errors, including the R5 increase, are reported separately in Figure 4 .

Figure 14: PIV errors over the complete time window. Spatial relative

\ell_{2}

errors (%) for each velocity component, the joint velocity vector, and vorticity at all 193 query times. R2 coincides with R1. These are single-problem grid diagnostics.

Figure 15: Final-program fields on the PIV-derived example. Reference, R5 prediction, and signed error (prediction minus reference) for speed and vorticity at $t=0.307$ , $0.928$ , and $1.55$ s. Each field uses one scale for reference and prediction and a separate symmetric error scale, fixed over the full time window. Gray vorticity pixels lie outside the common interior-stencil support.

Figure 16: Velocity error across PIV development rounds. Pointwise vector error in m/s for the retained programs R1–R5 on the same manufactured problem at three times. All panels share one linear scale determined over all five programs and all 193 times. R2 repeats R1 because its program is unchanged.

Figure 17: Reference macroscopic fields for the strong-collision VPFP example. The first unguided test problem, sampled on a $64\times 64$ grid over $[0,2\pi)^{2}$ at the five displayed times. Horizontal and vertical directions are $x$ and $y$ . Rows show density, electric-field magnitude, particle-flux magnitude, central-pressure Frobenius norm, and central-heat-flux magnitude, in nondimensional model units. Each row uses a common linear scale across times. These are manufactured reference fields.

Panel of Figure 18 — Figure 18: R1 and R6 errors on the strong-collision VPFP example. Top: Sol density error $|\widehat{\rho}-\rho|$ . Bottom: Astra central-pressure error $\|\widehat{\mathbf{P}}-\mathbf{P}\|_{\mathrm{F}}$ , counting both off-diagonal tensor entries. Within each block, rows are R1 and R6 and columns are $t=0,0.25,0.5,0.75,1$ . Each image covers $[0,2\pi)^{2}$ on a $64\times 64$ grid, with $x$ horizontal and $y$ vertical. Each block has its own linear scale, fixed over both programs and all nine saved times. The errors are in nondimensional model units; the two blocks have different quantities and scales.

Figure 19: Weak-collision VPFP heat-flux errors across development rounds. Sol’s R1–R6 programs on the first unguided test problem. Panels show central-heat-flux error $\|\widehat{\mathbf{q}}-\mathbf{q}\|_{2}$ at $t=0,0.25,0.5,0.75,1$ on a $64\times 64$ grid over $[0,2\pi)^{2}$ ( $x$ horizontal, $y$ vertical). One linear color scale is shared by all rounds and nine saved times. Errors use nondimensional model units.

PDE Solving

PDE solving develops numerical and learned methods to approximate solutions of partial differential equations across scientific applications.

See all 250 papers

AutoPDE: Reliable Agentic PDE Solving via Explicitly Represented Solver Strategies
Jun 9, 2026Huanshuo Dong, Keyao Zhang, Hong Wang +6AI Coding Agents PDE Solving
University of Science and Technology of China · 2Zhejiang University · University of the Chinese Academy of Sciences +2
PDEAgent-Bench: A Multi-Metric, Multi-Library Benchmark for PDE Solver Generation
May 10, 2026Zhen Hang, Yushan Yashengjiang, Junhui Li +21PDE Solving Benchmark Design
University of Science and Technology of China · Beijing University of Posts and Telecommunications · Tencent +10
Reinforcement Learning with Verifiable Physics: Post-training LLMs with Continuous Rewards
Jul 11, 2026Pengfei Cai, Utkarsh Utkarsh, Alan Edelman +2RL for Code Generation Scientific Code Generation
Massachusetts Institute of Technology

Jun 9, 2026J/K move · Enter open · S save

AutoPDE: Reliable Agentic PDE Solving via Explicitly Represented Solver Strategies

Huanshuo Dong, Keyao Zhang, Hong Wang, Zhezheng Hao+5

University of Science and Technology of China · 2Zhejiang University · University of the Chinese Academy of Sciences+2

AI Coding Agents PDE Solving

Open paper AskPDF arXiv

Numerical solvers for partial differential equations (PDEs) are core computational tools in science and engineering. Building reliable PDE solvers requires not only executable code, but a numerical solver strategy, a set of decisions about discretization, stabilization, solver configuration, and resolution control, that matches the PDE structure. Recent LLM-based coding agents have begun to reduce the programming burden by generating and debugging solver implementations. However, they typically move directly from a PDE problem to solver code, leaving the solver strategy implicit in implementation details. Feedback from a failed solve is therefore routed back to code edits rather than to the underlying strategy, so numerical decisions remain hard to check before code is generated and hard to revise using numerical evidence when it fails. To address this limitation, we propose AutoPDE, a code agent that maintains the solver strategy as an explicitly represented object throughout the solving process: an independent, inspectable object that is built before any code is written and can be revised, using numerical evidence, whenever a solve fails. AutoPDE builds and maintains this object in three stages, all drawing from a library of reusable PDE-solving skills: PDE analysis identifies the equation type and algebraic structure; numerical method selection chooses a numerical method that matches the analysis result and commits to a discretization, stabilization, and linear solver accordingly; and adaptive tuning runs low-cost pilot solves to calibrate resolution and tolerances under the prescribed accuracy and runtime budget. We evaluate AutoPDE on the PDE Agent Bench, where experimental results show that AutoPDE achieves a pass rate of

54.5%

, improving over the strongest baseline by

14.2

percentage points.

Abstract

PDE Solving

Explore similar work

AutoPDE: Reliable Agentic PDE Solving via Explicitly Represented Solver Strategies

PDEAgent-Bench: A Multi-Metric, Multi-Library Benchmark for PDE Solver Generation

Reinforcement Learning with Verifiable Physics: Post-training LLMs with Continuous Rewards