Paper deep dive

LLM-Assisted Generation of Transparent, Open-Source Multiphysics Models of Electrochemical Devices

Sebastian Castro, Maya F. Schuchert, Spencer A. McCluskey, Eric W. Lees, Justin C. Bui

Year: 2026Venue: arXiv preprintArea: physics.chem-phType: PreprintEmbeddings: 77

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%

Last extracted: 10/11/2026, 4:24:20 AM

Summary

This paper presents a framework for generating transparent, open-source multiphysics models of electrochemical devices using Large Language Model (LLM) agents. By utilizing a human-specified, machine-readable modeling harness containing governing equations and numerical methods, the authors demonstrate that an AI agent can reproducibly construct a 1D continuum model for CO2 reduction in Julia. The agent-built models show high agreement (within 0.7%) with commercial COMSOL simulations and high reproducibility across independent builds, establishing a more accessible and transparent alternative to proprietary software.

Entities (8)

COMSOL · softwaretool · 99%Julia · programming-language · 99%LLM Agent · softwaretool · 98%CO2 Reduction · processreaction · 95%Gas Diffusion Electrode · device-component · 95%Justin C. Bui · person · 95%Sebastian Castro · person · 95%Modeling Harness · methodology · 90%

Relation Signals (5)

Sebastian Castro → authored → LLM-Assisted Generation of Transparent, Open-Source Multiphysics Models of Electrochemical Devices

confidence 95% · Sebastian Castro 1 , Maya F. Schuchert 1 , Spencer A. McCluskey 2 , Eric W. Lees 2,* , Justin C. Bui 1,*

LLM Agent → generates → Julia

confidence 95% · From this specification, the agent reproducibly constructs complete multiphysics models in open-source Julia.

LLM Agent → validatesagainst → COMSOL

confidence 95% · Independently built models, including fully autonomous agent-built models, agree with an equivalent COMSOL implementation to within 0.7% of the peak CO partial current density

Modeling Harness → guides → LLM Agent

confidence 90% · This implementation guide served as a harness that directed the AI agent through construction of the multiphysics model in open-source Julia

Julia → usedfor → CO2 Reduction

confidence 90% · we develop a machine-readable, human-specified modeling harness... From this specification, the agent reproducibly constructs complete multiphysics models in open-source Julia.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Multiphysics continuum models are powerful tools for studying electrochemical devices, enabling in silico reactor design and resolution of local pH, potential, and concentration fields that govern device performance but are difficult to measure experimentally. However, constructing such models requires substantial numerical expertise or reliance on proprietary software. Here, we show that frontier large language model agents can remove this implementation burden while keeping the underlying physics under researcher control. Using one-dimensional electrochemical CO2 reduction to CO in a porous gas diffusion electrode as a test case, we develop a machine-readable, human-specified modeling harness containing governing equations, parameters, numerical methods, logical build stages, and human-verifiable checkpoints. From this specification, the agent reproducibly constructs complete multiphysics models in open-source Julia. Independently built models, including fully autonomous agent-built models, agree with an equivalent COMSOL implementation to within 0.7% of the peak CO partial current density, and with one another to within 0.04%. Systematically planted errors demonstrate the importance of explicit specifications for reproducibility and reveal the agent's capabilities and limitations in debugging model physics. This framework establishes a more transparent approach to multiphysics modeling in which physical descriptions and governing equations, rather than specialized code, become the primary inputs for computational model development.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

73,650 characters extracted from source content.

Expand or collapse full text

arXiv–ResearchArticle LLM-AssistedGenerationofTransparent,Open-Source MultiphysicsModelsofElectrochemicalDevices Sebastian Castro 1 , Maya F. Schuchert 1 , Spencer A. McCluskey 2 , Eric W. Lees 2,* , Justin C. Bui 1,* 1 Department of Chemical and Biomolecular Engineering, New York University, Brooklyn, NY, USA 2 Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, BC, Canada * Corresponding authors. e-mail: justinbui@nyu.edu; eric.lees@ubc.ca 1 arXiv:2610.10320v1 [physics.chem-ph] 7 Oct 2026 arXiv–ResearchArticle Abstract Multiphysics continuum models are powerful tools for studying electrochemical devices, enabling in silico reactor design and resolution of local pH, potential, and concentration fields that govern device performance but are difficult to measure experimentally. However, constructing such models requires substantial numerical expertise or reliance on proprietary software. Here, we show that frontier large language model agents can remove this implementation burden while keeping the underlying physics under researcher control. Using one-dimensional electrochemical CO 2 reduction to CO in a porous gas diffusion electrode as a test case, we develop a machine-readable, human-specified modeling harness containing governing equations, parameters, numerical methods, logical build stages, and human-verifiable checkpoints. From this specification, the agent reproducibly constructs complete multiphysics models in open-source Julia. Independently built models, including fully autonomous agent-built models, agree with an equivalent COMSOL implementation to within 0.7% of the peak CO partial current density, and with one another to within 0.04%. Systematically planted errors demonstrate the importance of explicit specifications for reproducibility and reveal the agent’s capabilities and limitations in debugging model physics. This framework establishes a more transparent approach to multiphysics modeling in which physical descriptions and governing equations, rather than specialized code, become the primary inputs for computational model development. 2 arXiv–ResearchArticle Electrochemical reactors sit at the center of the energy transition: electrolyzers use renewable electricity to drive the synthesis of fuels and chemicals, while fuel cells and batteries store that renewable electricity in chemical form 1,2 . Each of these devices is a coupled, multi-component system for which performance is an emergent property of the interactions between material components (electrodes, membranes, and electrolytes) and operating conditions (temperature, pressure, flow rate). The spatially heterogeneous pH, potential, and reactant concentration profiles within these reactors directly govern selectivity, throughput, and durability, but these local chemical environments are largely inaccessible to direct measurement under realistic operating conditions 3,4 . Consequently, rational reactor design and optimization are difficult to achieve through experimentation alone. Continuum modeling can accelerate the optimization of electrochemical devices by translating their coupled physics (i.e., multiphysics) into a system of governing equations that can be solved numerically 5,6 . Conservation equations for mass, charge, and momentum, expressed with models of ion transport, electrochemical kinetics, and fluid flow in porous media, capture the relationship between local chemical environments and experimental observables. For instance, multiphysics models can resolve how competing reactions distribute current across an electrode 7 , how membrane chemistry impacts voltage losses 8 , and how spectator ions tune reaction selectivity 9 . Once parameterized and validated against experimental data, they also enable large design spaces to be explored at a fraction of the cost of experiments. The history of electrochemical continuum modeling is a history of shifting barriers. In the 1960s, building a model required simultaneous mastery of the physics of electrochemistry and the numerical methods needed to solve the resulting stiff, coupled partial differential equations (PDEs), often using Fortran 10,11 , which confined serious modeling to a small community. Commercial simulation platforms, principally COMSOL Multiphysics, placed those algorithms behind a graphical interface and a proprietary solver 12 , allowing practitioners to focus more effort on defining the physical equations and less on implementation. Physical intuition remained essential to choosing and evaluating the governing equations, constitutive relationships, and boundary conditions; the numerical machinery, however, receded into the background, and with it direct engagement with how models are solved. When a solver fails, diagnosing whether the failure originates from the governing equations, discretization, or numerical implementation can therefore be difficult. Generic messages such as “maximum iterations reached” rarely indicate which equation, term, or mesh element is responsible, leaving debugging to heuristic experience with mesh refinement, relaxation factors, and initial guesses. Commercial solvers are also less accessible and more expensive than open-source alternatives 13,14 , and proprietary formats complicate reproducibility, as model files are often tied to specific software versions and inaccessible to collaborators without licenses 15 . We contend that the emergence of large language model (LLM) agents is the next chapter in electrochemical continuum modeling. A researcher with deep domain knowledge but limited numerical training can instruct an artificial intelligence (AI) agent to implement a fully open-source continuum model in which every equation, solver step, and line of code is visible, modifiable, and explainable. The agent does not replace the physics expertise of the researcher but rather replaces the implementation burden. The researcher must still specify the 3 arXiv–ResearchArticle governing equations precisely, anticipate which physical mechanisms are relevant, design staged verification checks that catch errors before they propagate, and recognize when a converged solution is unphysical. Domain knowledge remains critical; what changes is that the auxiliary technical burden shifts from Fortran subroutines and COMSOL debugging toward specification and verification. Because the resulting model is open-source and every step of its numerics is inspectable, the practical barrier to continuum multiphysics modeling becomes the ability to state the physics precisely rather than access to a solver license or extensive training in numerical methods. Agentic AI has already shown potential for integration into the scientific process, including autonomy in designing experiments, executing code, and discovering new pathways for chemical synthesis 16,17 . LLM agents have been shown to autonomously write, execute, and debug finite-element code for continuum mechanics problems 18 , collaborate as multi-agent teams across a range of physics domains 19,20 , and generate PDE solvers that match expert-written code on benchmarks including compressible Navier–Stokes 21 . Nonetheless, this work has largely targeted weakly coupled or single-physics problems, leaving strongly coupled multiphysics systems comparatively unexplored. A recent multi-agent framework, for instance, runs finite-element simulations from natural language by retrieving domain knowledge from a curated knowledge graph and verifying the solver against closed-form analytical solutions 22 . Neither approach transfers directly here: a nonlinear, strongly coupled system of this kind generally possesses no analytical reference solution, while its governing physics must be specified through researcher-authored modeling choices specific to the simulated device rather than retrieved as a fixed set of facts. Electrochemical reactors are a particularly demanding case, with acid–base equilibria, interfacial kinetics, migration, and multicomponent gas transport entering one another's source terms and transport coefficients. The equations must be solved simultaneously as a single nonlinear system, with buffer equilibria that relax orders of magnitude faster than species diffuse, making that system stiff and difficult to solve robustly. To demonstrate the validity of our approach, we developed a detailed, machine-readable implementation guide for a one-dimensional (1D) continuum model (Fig. 1) of silver-catalyzed CO 2 reduction (CO 2 R) in a gas diffusion electrode (GDE). This implementation guide served as a harness that directed the AI agent through construction of the multiphysics model in open-source Julia 23 . Researchers of varying experience levels used the guide to instruct AI agents to build the model, without writing a line of code themselves. Because the complex numerical system has no analytical solution, validation of the approach was demonstrated against an equivalent model developed in COMSOL. Reproducibility was assessed by comparing independently generated model implementations across users. Beyond accuracy and reproducibility, the agentic workflow provides a transparent interface to the numerical implementation: the agent that built the code can inspect the Jacobian, trace convergence failures, and explain the numerical behavior in the language of physics directly to a human operator. Collectively, this framework establishes an approach to multiphysics modeling in which human-specified equation sets and natural-language instructions, rather than specialized code and numerical implementations, become the primary inputs for model development. In this paradigm, the specification harness—rather than any individual implementation—becomes the central scientific artifact that other practitioners can audit, reproduce, modify, and build upon. 4 arXiv–ResearchArticle Implementation Guide User creates and refines guide Defines physics parameters, methods Agent reads guide Collaborative Model Building Agent implements code and runs model stages Researcher verifies output Equations Parameters Numerical Methods Build Stages Researcher Provides physical insight, context, clarification Provides numerical insight and backend Continuum Model LLM Agent CO 2 (g) Carbon GDL Ionomer CLGDLElyte CO 2 Catalyst (b) 10 μm325 μm COH 2 H 2 O H 2 O =∇∙ CO H 2 H + OH – HCO – 3 CO 2– 3 CO 2 Cs + (aq) Stage 1 Liquid DOFs, diffusion only Dirichlet BCs Stage 2 Mass transfer at electrolyte BC Gas DOFs + phase transfer Stefan Maxwell diffusion Stage 3 Darcy pressure Finalize gas momentum balance Stage 4 CO 2 R: i 0 + λ HER: A HER + α c Kinetic Parameter Fit Potential DOFs Full Nernst–Planck migration Stage 5 (a) Gas Ag Catalyst Increasing Complexity Stop gates for human review Fig. 1. (a) Schematic of the agentic multiphysics modeling workflow employed here with a human-in-the-loop. (b) Schematic of the CO 2 reduction porous electrode model case study, along with the 5-stage implementation approach in Julia, with the primary physics additions/changes included in each stage. Elyte, electrolyte; CL, catalyst layer; GDL, gas diffusion layer; DOFs, degrees of freedom; BC, boundary condition; HER, hydrogen evolution reaction. Results Aspecification-driven,stagedbuild The modeling of electrochemical reactors is a prime example of a complex, multiphysics problem, and thus presents a suitable platform for establishing agentic multiphysics approaches. Every one of the nonlinear couplings introduced above, and the corresponding numerical strategy that makes the resulting stiff system tractable, must be stated explicitly before an agent can implement it. Hence, an implementation guide (i.e., machine-readable specifications, here markdown files) can serve as a harness for an agent building multiphysics models, fully describing the problem to be modeled: all relevant parameters, equations, and boundary conditions, together with a scaffolded build plan that the coding agent must follow. These guides link the scientific knowledge and intuition from the human to the numerical implementation of the agent. Precise and comprehensive specifications ensure consistency and accuracy between agents with probabilistic behavior, and the responsibility remains with the researcher to accurately define the problem. In other words, the agent's primary role in this workflow is the implementation of a specified guide, rather than the inference of a complete model from an incomplete specification. To provide the agent (and its human operator) with logical checkpoints at which to assess and correct mistakes in implementation, the guide decomposes model development into five stages of increasing physical complexity (Fig. 1b and Supplementary Fig. 1), from liquid Nernst–Planck transport, buffer reactions, and electrochemical kinetics at 400 degrees of freedom (DOFs) to the full electrostatically coupled system at 1,480 DOFs on the production mesh. Between stages, stop gates allow the human-in-the-loop to review quantitative physicality checks and diagnostic 5 arXiv–ResearchArticle plots prescribed by the guide itself, rather than leaving them to the agent's discretion, before authorizing the agent to advance to the next stage. These gates are solution-verification checks in the established sense 24 , applied by the human rather than the agent. Full stage definitions, pass/fail criteria, and details of the fitting step used to parameterize the model against the experimental data of Romiluyi et al. 25 are given in the Methods section. Our assessment of the AI-assisted build begins with the polarization and faradaic efficiency (FE) results, and how these outputs evolve as each stage of the implementation guide adds a layer of physics (Fig. 2a–c). Notably, the plots reproduce the canonical behavior of a CO 2 electrolyzer: FE CO peaks at intermediate potentials and declines rapidly at high potentials where the hydrogen evolution reaction (HER) outcompetes. Physically, this occurs because of the interplay between 1.01.52.02.53.0 −E [V vs SHE] cath 0 100 200 300 400 500 600 700 i [ m A / c m ] C O 2 a) Model Experiment 1.01.52.02.53.0 −E [V vs SHE] cath 10 −1 10 0 10 1 10 2 10 3 i [ m A / c m ] C O 2 b) 1.01.52.02.53.0 −E [V vs SHE] cath 0 20 40 60 80 100 F E [ % ] C O c) 0246810 x [μm] 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 pH d) 0246810 x [μm] 0.002 0.004 0.006 0.008 0.010 0.012 C [ m o l / L ] C O ( a q ) 2 e) 0246810 x [μm] 0.0 0.2 0.4 0.6 0.8 1.0 C [ m o l / L ] H C O 3 − f) 0246810 x [μm] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 C [ m o l / L ] C O 23 − g) −300−250−200−150−100−500 x [μm] 0.8600 0.8625 0.8650 0.8675 0.8700 0.8725 0.8750 0.8775 0.8800 p [ a t m ] C O 2 h) −300−250−200−150−100−500 x [μm] 0.000 0.005 0.010 0.015 0.020 p [ a t m ] C O i) Fig. 2. Stage-by-stage build progression and internal fields of final model. CO (a) linear-scale polarization, (b) log-scale polarization and (c) faradaic efficiency after each build stage and the fit. Converged Stage-5 profiles at five applied potentials: across the catalyst layer, (d) pH, (e) CO 2 , (f) HCO 3 – and (g) CO 3 2– concentrations; across the gas diffusion layer and catalyst layer, (h) CO 2 and (i) CO partial pressures. Markers represent experimental data. Voltages vs the standard hydrogen electrode (SHE). 6 arXiv–ResearchArticle CO 2 reduction kinetics, CO 2 mass transport, and the consumption of CO 2 by acid–base equilibria within the porous electrode 26 . The evolution of the curves coincides with the physics added at each stage: Stage 1 includes buffer equilibria, electrochemical kinetics, and liquid-phase species transport by diffusion alone, with fixed-concentration (Dirichlet) boundary conditions; the Stage 2 mass-transfer boundary condition allows CO 2 depletion and pH buildup and correspondingly reduces CO partial current, i CO ; Stage 3 phase transfer from the gas adds a liquid-phase CO 2 source and drives i CO back up; Stage 4 momentum balance in the gas phase adds a small but physical correction; the fit, performed on the Stage 4 model, brings i CO down to the experimental data; and Stage 5 migration raises i CO to its final value by increasing the flux of OH – and other anions away from the cathode, which mitigates the salting-out of CO 2 . Beyond examining the polarization data, the correctness of the implemented physics was complementarily verified at each stage by analyzing the species profiles to confirm the simulated device concentration fields are consistent with physical intuition. The most relevant field profiles are shown in Fig. 2d–i, and the full output of all the plotted profiles from every stage is shown in Supplementary Figs. 3–13. As expected, the simulated pH (Fig. 2d) is highest at the gas diffusion layer (GDL) side (x = 0) where the CO 2 consumption rate is highest, generating OH – anions via CO 2 R. Additionally, the dissolved CO 2 (Fig. 2e) is lowest near the GDL where CO 2 R rates are highest and increases toward the electrolyte, where CO 2 is slightly replenished by incoming (bi)carbonates and their equilibria. The (bi)carbonate field profiles (Fig. 2f–g) behave in accordance with their acid–base equilibria and CO 2 concentrations. In the gas phase (Fig. 2h–i), gaseous CO 2 is highest where it comes in from the flow channel, has a linear concentration gradient through the GDL where it undergoes diffusive transport only, and is consumed in the catalyst layer (CL). CO is generated in the CL and diffuses out into the gas channel. These profiles, which would be difficult to measure experimentally in the microscale porous electrode, show that the model is physically reasonable. AgreementwithCOMSOLandacrossbuilds To assess the efficacy of the specification guide in facilitating agentic multiphysics modeling, we leveraged the build to develop an experimentally calibrated continuum multiphysics model of a gas diffusion electrode used for CO 2 R. We consider an AI-assisted build successful if it yields consistent and physically reasonable results across independent builds, reproduces the simulation results of an established commercial tool reference (in this case, COMSOL) to within a few percent, and remains fully transparent end-to-end (i.e., all code is generated and available for examination by human evaluation or the audits of other agents). To confirm the accuracy of the agent’s multiphysics implementation, we benchmarked the three independent human-in-the-loop builds, made with Claude Opus 5, against a physically equivalent COMSOL model with the same parameters (Fig. 3). The Julia build outputs agree with the COMSOL model quantitatively, with the polarization and FE curves overlapping entirely. Evaluated at the 44 COMSOL sweep potentials that fall within the Julia sweep window, the average of the root mean square deviations (RMSDs) of each build against COMSOL is 2.27 mA cm −2 for i CO , 1.07 mA cm −2 for the H 2 partial current density, i H2 , and 0.14 percentage points for FE CO . The i CO RMSD represents 0.7% of the peak CO current density of about 330 mA cm −2 , and the Julia i CO 7 arXiv–ResearchArticle 1.01.52.02.53.0 −E [V vs SHE] cath 0 200 400 600 800 1000 1200 i [ m A / c m ] 2 a) i (Rep 1) CO i (Rep 1) H 2 i (Rep 2) CO i (Rep 2) H 2 i (Rep 3) CO i (Rep 3) H 2 i (COMSOL) CO i (COMSOL) H 2 i (exp) CO i (exp) H 2 1.01.52.02.53.0 −E [V vs SHE] cath 0 20 40 60 80 100 FE [%] b) 0246810 x [μm] 0.002 0.004 0.006 0.008 0.010 0.012 C [ m o l / L ] C O ( a q ) 2 c) Julia (Reps 1–3) COMSOL 1.01.52.02.53.0 −E [V vs SHE] cath −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Deviation from COMSOL [%] d) i CO i H 2 COMSOL 1.01.52.02.53.0 −E [V vs SHE] cath −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 F E d e v i a t i o n f r o m C O M S O L [ p p ] C O e) 0246810 x [μm] −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 C d e v i a t i o n f r o m C O M S O L [ % ] C O ( a q ) 2 f) Fig. 3. Agreement of three independently built Julia models with the COMSOL model. (a) CO and H 2 partial current densities, (b) faradaic efficiencies and (c) CO 2 (aq) concentration across the catalyst layer at five applied potentials, with the deviation of each from COMSOL below it: (d) currents, (e) FE CO , in percentage points, and (f) CO 2 (aq). Current and FE deviations are evaluated at the COMSOL sweep potentials. Markers represent experimental data. values sit systematically above the COMSOL values, by up to 2.2% at currents above 5 mA cm −2 and about 0.9% at the peak, a discrepancy likely due to slight differences in mesh generation or to numerical implementation differences between the solvers. Across the three builds evaluated on a common 89-point voltage sweep, the partial current densities also agreed with one another to within a pooled RMSD of 0.12 mA cm −2 for CO and 0.11 mA cm −2 for H 2 , and the CO faradaic efficiency agreed within 0.01 percentage points. The fitting approach was also remarkably consistent across the three builds. The CO 2 R Marcus–Hush–Chidsey 27 kinetic parameters were fit to literature data using a Nelder–Mead algorithm 28,29 , landing on an exchange current density (i 0 ) of 106.47 ± 0.11 A m −2 and a reorganization energy (λ) of 1.6111 ± 0.0004 eV. The small residual spread reflects the flatness of the fit objective; all three replicates used identical termination criteria and reached the same sum of squared errors to four significant figures (0.05608). These results confirm that independent implementations of the specification harness produce nearly identical multiphysics models both to each other and to an established commercial software package, providing confidence in both the methodological approach and the accuracy of the physical implementation. Solverimplementationandagentcost–benefitdeviation One benefit of commercial packages is that their plug-and-play numerical solvers have been long-optimized to ensure that solvers are properly matched to the systems that they solve. To 8 arXiv–ResearchArticle benchmark the quality and speed of the agent’s numerical implementation of CO 2 R against that of a commercial software package, we performed speed benchmarking against an equivalent COMSOL model with matched degrees of freedom. Notably, the agent’s Julia model demonstrates competitive solve speed to COMSOL with the added benefit of customization. Timed on the same machine at a matched compute budget and problem size (~20,000 DOFs), the full 89-potential sweep takes about 13 s with our production solver stack compared to 21 s for COMSOL’s 61-potential sweep (Fig. 4a). It is important to emphasize, however, that a direct comparison is not possible because of inherent differences in discretization (finite element method versus finite volume method), linear-solve library (PARDISO 30 vs UMFPACK 31 ), and solve environment (native Windows vs the Windows Subsystem for Linux). The solver ablation (Fig. 4b) reveals that speed gains in numerical solvers often come from simply tailoring the solver (e.g., its sparsity, the discretization method, etc.) to the problem solved. Agentic multiphysics implementations with frontier agents therefore benefit significantly from the cheapness of customizing solver parameters to the specific problem being solved. The gains come almost entirely from exploiting model sparsity, first in the linear solve and again in Jacobian assembly, each worth about two orders of magnitude against a dense alternative. We note that the components of the developed solver do not represent novel mathematics or solution techniques. However, with an open-source, modifiable solver, we no longer need to sacrifice performance for the solver generality that is often required of one-size-fits-all software packages. Full details of the solver stack are given in Supplementary Note 4. Despite agents mostly conforming to the specification, we observed deliberate, documented deviations during intermediate solver development. Each solver component was explicitly mandated in a reproducibility checklist, yet none of the builds implemented the analytical block COMSOL 32 cores Julia 32 threads 0 5 10 15 20 25 Solve time, ~20,000 DOFs [s] a) 21 s 12.8 s 10 1 10 2 10 3 10 4 Solve time, 20,038 DOFs [s] b) linear solve dense LU Jacobian dense AD predictor none linear solve UMFPACK, re-analyzed linear solve persistent MUMPS Jacobian analytic block reference production stack 42 min (198×) 31 min (147×) 23.4 s (1.83×) 19.6 s (1.53×) 19.3 s (1.51×) 13.0 s (1.02×) 12.8 s 10 2 10 3 catalyst-layer cells N CL 10 1 10 0 RMSD of i CO vs N CL limit [mA cm 2 ] c) N 1 CL 1.01.52.02.53.0 E cath (vs SHE) [V] 0 50 100 150 200 250 300 350 i CO [mA cm 2 ] d) 1.92.02.1 328 330 332 334 336 20 40 80 160 320 640 1000 N CL (coarse fine) Fig. 4. Solver performance of the Julia implementation. (a) Solve time for each implementation's full sweep (89 potentials in Julia, 61 in COMSOL), both run on the same machine at matched core budget (32 cores/threads, ≈20,000 DOFs). (b) Measured cost of Jacobian and linear-solver choices at 20,038 DOFs. Every configuration returns the same polarization curve to the precision written, so bars differ only in cost. Ratios are relative to the configuration specified in the guide (UMFPACK with pattern reuse and forward-mode automatic differentiation). LU, lower–upper factorization; AD, automatic differentiation. 9 arXiv–ResearchArticle Jacobian or installed MUMPS 32 , the open-source sparse direct solver that an earlier version of the guide specified, unless explicitly told to do so for speed benchmarking. The guide itself labeled these components performance multipliers, and the agents took the label at face value, opting for simpler alternatives and recording the reasoning in their build dialogue notes. The agent would weigh, for example, hand-deriving the analytical block Jacobian, which can yield silent calculation errors, against a forward-differentiation call 33 that takes a few lines of code and is correct to machine precision. Automatic differentiation has been adopted for the same reason in Newman-tradition porous electrode codes 34 , where hand-derived Jacobians are a known source of silent error. Similarly, the agents opted for the UMFPACK solver available in the Julia standard library, justifying the substitution as having similar performance while avoiding the extra effort required for installing MUMPS. Regardless, the difference in solve time for the solver items substituted by the agent was non-detrimental (Fig. 4b), while the implementation burden for the agent would have been intensive. Because both substitutions were at least as fast as the specified components, they were adopted in the final version of the solver stack. This behavior is nonetheless a cautionary tale for other users of this methodology: unless explicitly told and checked by a human operator, frontier agents often route numerics through simpler implementations that they deem equivalent and more efficient. Diagnostictransparencyduringmodeldevelopment Building a correct model depends first on the front-end work of developing and specifying the relevant physics and numerical methods within the implementation guide. Once the guide is sufficiently specified and could generate models, the process consisted of building, troubleshooting, refining, and iterating on the model harness. Intermediate builds using earlier versions of the implementation guide surfaced several critical errors, required substantial human input, and/or produced significant variation in build trajectories even from the same guide (Supplementary Note 7), reflecting the importance of the user in adequately specifying the model physics as to reduce variance in the LLM's probabilistic output 35 . As the guide was refined, these issues were progressively designed out and constrained, so that the final replicates required only minimal human intervention and had no critical failures. The challenges encountered throughout the development process for the implementation guide, along with the agents’ capacity to detect and address these failures, are shown in Fig. 5. An example of a physics bug that arose in the early builds occurred when switching from the Dirichlet to the Neumann mass-transfer boundary condition for CO 2 in Stage 2, where the model suddenly could not converge. The agent reasoned that the mass-transfer boundary condition was unable to sustain high currents, and, upon investigating the physics, found a single sign error in the mass-transfer flux term as the root cause, since it incorrectly caused CO 2 to leave the system rather than being supplied from the electrolyte side. The same sign error recurred independently in five of six early builds, across two agent models and four users, because the guide gave the flux without properly specifying its sign convention. In all runs the solver failed to converge, prompting investigation. In some cases, the agent found the bug as discussed above, while in others it could not identify the root cause without human intervention (e.g., prompting the agent to check the sign on the boundary condition). 10 arXiv–ResearchArticle ErrorSymptomAgent Handling Natural: early development Planted error study: mature guide CLElyte J i c i Case 1: Sign error in Stage 2 mass transfer boundary condition Case 3: Porous media diffusivity correction missing Case 2: Cs + omitted from ionic current sum residual iteration c i x No convergence at any voltage 10 attempts to resolve error Caught independently Not caught Caught with help Caught independently 1: Backtracking line search 2: 훼 BC blending (stalled at 0.97) 3: k MT scale ramp ... 9: Pseudo-arclength continuation echem ramp 10: PAC V-sweep from V = +1 Mass balance not satisfied Liquid potential drop too small Smoother profiles, but plausible results Pressure drop becomes negligible Identified mass transfer limitation from BC sign convention Compared to Nernst-Planck sign convention defined elsewhere Performed OH – physicality check D i,w D i,eff D i 3.6 × too high Silent during build; found only after cross-comparing builds Compared Ф L profiles vs reference Low Biot number fell within expected magnitude, produced plausible results and no convergence errors Nonconvergence exposed problem Not applicable Investigated Jacobian conditioning first, then caught physics bug after user nudge: “physics bug in the liquid block” ε 1.5 ⇄ Cs + i flux cond(J) ~10 23 Ф L x ! Parameters identical; traced root cause by reading Ф L assembly line by line residual iteration Tolerance too high Partial convergence; Kinks in polarization curve i CO –V Investigated after user observation; found high residual norm ‖F‖ ∞ ~ 10 –3 Tightened tolerance to 10 –6 and increased maximum iterations; converged with ‖F‖ ∞ ~ 10 –8 M [kg/mol]M [g/mol] FSG Correlation D gas D gas 30 × too high p i x Performed order-of-magnitude sanity check against known gas diffusivities Case 4: Incorrect units passed into FSG correlation tol. Observed higher CO currents and low CO 2 gas depletion at high overpotential. Found unit error and verified against literature diffusivity Fig. 5. How does the multiphysics agent handle errors? Errors encountered naturally during model development and planted systematically into a mature guide, their respective symptoms that the human can see and report to the agent, and the agent’s response to resolve the issue. The planted error study was conducted with Claude Opus 4.8. Agents were instructed not to consult outside resources, including other model builds, unless stated otherwise (see Supplementary Notes 6 and 8). Elyte, electrolyte; CL, catalyst layer; BC, boundary condition; α BC , boundary-condition blending factor; k MT , mass-transfer coefficient; PAC, pseudo-arclength continuation; FSG, Fuller–Schettler–Giddings; Φ L , liquid-phase potential; cond(J), condition number of the Jacobian. 11 arXiv–ResearchArticle Even when the physics is implemented correctly, convergence errors are common in numerical models and are the most common culprits that stall model development. In two independent late-stage builds from different users, spurious points in the polarization curve appeared that partially converged but were not on the correct solution path. In both cases, the human-in-the-loop identified the issue through visual plot inspection and communicated it to their agent, which found that loose convergence criteria were admitting partially converged points. The agents autonomously converged on the same fix, tightening the absolute residual tolerance. Nonetheless, the agents struggled to identify hidden physics bugs that do not block convergence, requiring human input to eliminate these physical inaccuracies. These silent errors are the most dangerous, since they can yield results that appear reasonable while being inaccurate. For example, an early guide included two expressions for the effective diffusivity of liquid-phase species, one with the Bruggeman correction for porous media and one without, resulting in effective diffusivities that differed by a factor of 3.54. In this case, none of the agents mentioned the ambiguity. Fortunately, when the agent was then instructed to cross-compare two replicates that had significant deviations in CO partial current density, it identified the difference in the parameter definition line by line. Notably, the ease of comparison was made possible by the standardization of the modular codebase in the specification guide. The agent could explain the physical implications of the differences in CO 2 diffusion coefficients, tracing the user's observation of higher CO partial currents to the fitting optimizer compensating for the decreased diffusivity by increasing exchange current density. It also located the two conflicting definitions in the guide itself regarding the use of the Bruggeman correction, which were reconciled in the next iteration of the specification. Because the agent writes every line of the solver, it does not merely return a number, exit flag, or convergence message; it can narrate the numerics in the language of the physics. During our builds, the agent could inspect the 1,480×1,480 sparse Stage-5 Jacobian and report its condition number, identify that the cell at the gas-side boundary carried the highest CO 2 R reaction rate while the electrolyte-side cell was HER-dominated, and explain why in terms of local CO 2 availability and pH, and verify current continuity against the boundary flux. This diagnostic ability represents a key benefit of having a specialized agent for each multiphysics problem as opposed to relying on a generalized software package. When using commercial software, errors can be difficult to self-identify, leading to reliance on reporting to a dedicated support team. Conversely, agents can point to specific variables and their residuals in a given mesh element, trace that to the exact physics, and tune solvers for convergence barriers such as stiffness and near-zero concentration terms. Errorhandlingwithmultiphysicsagents To test the agent's ability to detect and resolve issues in a controlled manner and establish the robustness of the guide itself, we conducted a planted-error study. Unlike the natural errors above, which were a result of and could be compounded by guide incompleteness or ambiguity, these were systematically placed in separate copies of a mature pre-release version of the guide to assess the model’s ability to identify and handle these errors. For four planted errors (referred to as Cases 1–4), the agent was not told which relation had been modified and was prompted to build without referring to any files besides the specification, just as any other build. A fifth build used the 12 arXiv–ResearchArticle unmodified guide as a control. All builds used Claude Opus 4.8 at maximum effort and received the same prompt. The four defects included a sign reversal in a boundary condition (Case 1, Fig. 5 row 1), an omitted term in a conservation law (Case 2, Fig. 5 row 2), an omitted correction factor in a constitutive relation (Case 3, Fig. 5 row 3), and a unit inconsistency in an empirical correlation (Case 4, Fig. 5 row 4). Each reproduces an error encountered naturally during development rather than one invented for the study (Supplementary Note 7). However, sandboxing these errors intentionally provides the cleanest possible study of how the agent addresses them, as opposed to allowing them to occur naturally due to unintentional guide ambiguities. Three of the four defects were resolved before the final stage of model development. Of these three, two were resolved without human intervention (Cases 1 and 4), and one was addressed only after a hint from the human-in-the-loop (Case 2). The fourth was never detected (Case 3). Detection tracked whether the guide contradicted itself, not the severity or even the visibility of the symptom. Case 1's sign reversal made the boundary-value problem unsolvable, and the sweep failed outright. Case 2 made the Jacobian structurally singular. Case 4 produced transport coefficients thirty-fold higher than their correct values and the agent rejected them on magnitude at Stage 3. In principle, Case 2 should have been the most consequential (breaking a conservation law) and was also the most conspicuous, but a visible symptom alone was not sufficient to be detected by the agent. Such errors often required a second textual signal, for instance, the modified guide contradicting itself. Case 1's reversed flux had the opposite sign to the Nernst–Planck boundary condition stated elsewhere in the same document, and Case 4's added units statement was incompatible with the Fuller–Schettler–Giddings correlation's own prefactor. Both errors were caught unaided. Case 2 left no such contradiction accessible to the agent and thus needed a hint from the operator. Of the four planted errors, Case 3 (omission of porous media correction for diffusivities) is the most informative failure, and the one that establishes the limits of the agents’ diagnostic capability. It is a plain constitutive error, no subtler in isolation than the others, but it was the only defect propagated consistently through every statement, example, and derived number in the guide. It was also numerically silent: the residuals converged at every stage, global charge conservation held to 2.4 × 10 –11 , and the one diagnostic that did register it, a liquid-potential drop at the bottom edge of the expected band, was explained away using the model's own inflated conductivity. Because the error could be compensated with the choice of kinetic parameters, the fitting step absorbed most of it, and the delivered model still reproduced the experimental polarization while carrying an exchange current density 41% below every other build. This is the error class that staged construction does not close on its own and requires discerning human rigor. This class of error is also quite common and is not specific to this multiphysics problem. A wrong constitutive relation transcribed faithfully and absorbed by a fit is a hazard of any model-building. The remedy is external reference, such as an independent implementation, an analytical limit, or a measurement the fit was not trained on. In this case, having the COMSOL model as an independent, orthogonal reference allowed the team to resolve silent errors in the implementation guide that caused deviations between the guide and the COMSOL implementation, which was treated as a “ground truth” or golden model. The study also shows that the build process is resilient to minor errors that would otherwise have major impacts. Beyond syntax errors, the agents caught common mistakes that hinder model 13 arXiv–ResearchArticle development and are tedious to find manually, including sign, unit, and math errors. Human mistakes are always a possibility regardless of experience level, so agents that can dissect code and dozens of parameters and equations in a fraction of the time provide a significant safety net. All syntax errors and most convergence errors arising from code implementation were caught by autonomous agent review, without human input. Fullyautonomousbuilds With human review eventually finding nothing else to catch during builds, the next step was pushing the bounds of agent autonomy. We gave the finalized implementation guide and the same prompt to three separate instances each of Claude Opus 5 and Sonnet 5 and instructed them to bypass the stop gates completely for a fully autonomous build with no human input. The agents were instructed to not refer to anything other than the guide, isolating the autonomous agents from previous builds and COMSOL data. The agent’s project memory, which carried a few results of earlier builds, was active in these builds as in all others (Methods); a tenth autonomous build made without it fitted i 0 = 106.59 A m −2 andλ= 1.6109 eV, inside the range of all nine builds. They checked their own plots and validated each stage's output against the physical and numerical checks provided in the guide. All six builds produced virtually the same output (Fig. 6 and Supplementary Note 8), replicating the canonical model with CO 2 R fitted parameters of i 0 = 106.56 ± 0.07 A m −2 andλ= 1.6113 ± 0.0002 eV for Opus builds and i 0 = 106.51 ± 0.01 A m −2 and λ= 1.6108 ± 0.0003 eV for Sonnet builds (Supplementary Table 9). These results show that LLM-assisted model building can culminate in fully autonomous builds once the implementation guide is fully developed. We emphasize that this was achievable only after substantial guide iteration. Nonetheless, the success of the Sonnet builds demonstrates a major payoff of the refinement, which is that less powerful models can do the same task at a cheaper rate. Usage costs can limit the ability of users to perform complex or large tasks within a certain time frame, but with well-defined tooling and harnesses (in our case, the specification), the scope of ambition can be made larger with cheaper agents. For instance, autonomy allows agents to run 1.01.52.02.53.0 −E [V vs SHE] cath 0 200 400 600 800 1000 i [mA/cm²] a) Opus 5, autonomous (n=3) Sonnet 5, autonomous (n=3) COMSOL i CO i H 2 i (exp) CO i (exp) H 2 1.01.52.02.53.0 −E [V vs SHE] cath 10 −1 10 0 10 1 10 2 10 3 i [mA/cm²] b) 1.01.52.02.53.0 −E [V vs SHE] cath −8 −6 −4 −2 0 2 4 6 8 i d e v i a t i o n v s C O M S O L [ % ] C O c) Opus 5 (n=3) Sonnet 5 (n=3) Fig. 6. Autonomous agent builds without a human-in-the-loop. Comparison of Claude Opus 5 and Sonnet 5 build replicate (a) linear-scale and (b) log-scale polarization curves to COMSOL and experimental data. (c) Deviation of each build's i CO from COMSOL at the COMSOL sweep potentials. 14 arXiv–ResearchArticle replicates or comparative studies in parallel without constant supervision, yet it can lead to fast usage depletion or high costs if used with stronger models like Opus and Fable. Cheaper agents, like Sonnet, require a fraction of the cost and can yield similar results depending on the task. While achieving autonomy for model building is a great milestone, it does not eliminate the risk that an autonomous build can silently inherit whatever is wrong with its specification guide, once again emphasizing the importance of the human in establishing a stable implementation guide to harness the multiphysics agent. Discussion This work extends LLM-assisted scientific computing beyond benchmark PDEs and weakly coupled or single-physics problems to a true multiphysics model of an experimentally relevant electrochemical device. Although the model is one-dimensional, it is not numerically simple: Nernst–Planck transport, carbonate buffer chemistry, electrochemical kinetics, Stefan–Maxwell gas transport, Darcy flow, and electrostatics are coupled through nonlinear source terms and transport coefficients in a stiff system containing 1,480 degrees of freedom. The resulting model reproduces expected physics and was calibrated to real experimental polarization and selectivity data and verified against an independently implemented COMSOL model. That provides a substantially stronger test of whether an agent can construct a model whose coupled physics actually reproduce the behavior of a real electrochemical system. More broadly, the results show that LLM-assisted multiphysics modeling can shift much of the technical burden from numerical implementation to physical specification and verification. The researcher still has to define the governing physics, choose the constitutive relationships and boundary conditions, recognize unphysical behavior, and validate the model. These decisions are where domain knowledge enters the model. For instance, we describe the porous electrode as a macrohomogeneous medium rather than a pore-resolved domain because that description resolves the polarization, selectivity, and concentration profiles studied here. The agent can implement either description, but the assumptions and simplifications appropriate to the question being asked remain the choice of the researcher. What the agent removes is the requirement that the researcher also translate that understanding into specialized numerical code. The implementation guide therefore becomes the central scientific artifact in this workflow. Supplementary Table 12 summarizes the components that most strongly determined whether independent agents built the intended model rather than a plausible approximation, including complete governing equations and boundary conditions, explicit unit, sign, and reference-frame conventions, numerical methods, staged verification criteria, and expected physical magnitudes. These requirements are not simply instructions for an AI agent. Together, they provide an explicit record of the assumptions, numerics, and physical reasoning underlying the model that can be inspected, revised, and shared independently of any particular implementation. The limitations are equally important: staged verification catches many implementation and convergence failures, but a physically incorrect assumption stated consistently throughout the specification can still be reproduced faithfully and even obscured by parameter fitting. Independent builds, external references, and human physical intuition therefore remain important safeguards. 15 arXiv–ResearchArticle The opportunity is therefore not to remove the researcher from continuum modeling, but to make the researcher’s understanding of the physics the primary interface to it. The models produced here are open source and fully inspectable, and the agent that builds them can trace numerical behavior back through residuals, local concentrations, governing equations, and physical mechanisms in the same language the researcher uses to reason about the system. Researchers who understand an electrochemical device but are not specialists in numerical methods can therefore build, troubleshoot, extend, and share rigorous multiphysics models without surrendering control of either the physics or the implementation. By replacing a barrier of specialized code and proprietary software with one of stating and evaluating the physics precisely, agentic AI can make transparent multiphysics modeling accessible to a much broader scientific community, expanding the role of multiphysics models in the future of science and engineering research. Methods Modelsystemandgoverningequations The model developed herein is based upon prior work on modeling electrochemical CO 2 reduction in 1D and, importantly, presents a stiff, highly nonlinear and coupled testbed with which to evaluate the efficacy and capabilities of our agentic multiphysics approach. Evaluating the effectiveness of building highly coupled and stiff numerical models with LLM agents requires an existing solution against which to verify the implementation. We therefore rebuilt a CO 2 R porous electrode model assembled from three previous COMSOL-based studies (Fig. 1b). The domain structure and GDE formulation follow Weng et al. 7 ; the Marcus–Hush–Chidsey (MHC) kinetics 27 follow Lees et al. 36 ; and the Ag-catalyzed CO 2 R model of Muhieddine Orfali et al. 37 serves as the COMSOL reference that we rebuilt for benchmarking. We resolve the cathode catalyst layer (CL) and gas diffusion layer (GDL) only and do not model the anode or membrane. The cathode alone is a sufficiently demanding multiphysics problem for a first demonstration, and full-cell models and models beyond 1D will be subject of future work. This system possesses each of the couplings described above. Two electrochemical reactions proceed at the catalyst surface, CO 2 + H 2 O + 2e – →CO + 2OH – (1) 2H 2 O + 2e – →H 2 + 2OH – (2) and both generate hydroxide, which is consumed by the carbonate buffer equilibria in the electrolyte, CO 2 + OH – ⇌HCO 3 – (3) HCO 3 – + OH – ⇌CO 3 2– + H 2 O(4) closing a loop in which current raises the local pH, the elevated pH consumes dissolved CO 2 , and the resulting CO 2 depletion suppresses the kinetic CO 2 R rate that produced the current as per mass transport limitations. Liquid-phase transport is described by Nernst–Planck 16 arXiv–ResearchArticle dilute-solution theory with acid–base and electrochemical source terms, gas transport by Stefan–Maxwell concentrated-solution theory with Henry’s law phase transfer, and gas-phase momentum by Darcy’s law; CO 2 R is modeled with MHC kinetics, evaluated with the closed-form approximation 38 , and the competing hydrogen evolution reaction (HER) with Butler–Volmer kinetics. The model is one-dimensional, steady-state, and isothermal at 50 °C, with fixed porosity and liquid saturation. The CL is treated as three phases, solid catalyst, ionomer/liquid, and gas-filled pore, each with fixed volume fractions defined by experimentally measured porous electrode structure. The GDL is treated as dry, and crossover to the anolyte is not modeled. The complete governing equations and boundary conditions are given in Supplementary Notes 1 and 2, the parameter set in Supplementary Note 9 and Supplementary Table 11, and the COMSOL build used for benchmarking in Supplementary Note 3. Theimplementationguide Critical to building well-developed, rigorous models are equally well-developed implementation guides for the agent. Implementation guides are machine-readable specifications (often in the form of markdown files) that fully describe and specify the problem to be modeled, including all relevant parameters, equations, boundary conditions, and include a scaffolded build plan that the coding agent must follow. Essentially, the implementation guide is a human-specified harness to constrain the multiphysics agent. These guides link the science from the human to the numerical implementation of the agent. Ensuring these guides are thorough and precise is necessary to prevent ambiguity for the agent that might lead to errors and inconsistencies in building the model. Hence, the human operator’s understanding of the physics and chemistry of a given multiphysics problem remains crucial to developing the model. In this workflow, the agent’s primary role is not inference but implementation, removing the numerical implementation barrier that otherwise hinders scientific modeling and has been largely inaccessible to chemists in the past. The guide used in this case study includes 18 sections covering everything needed for the full implementation of the multiphysics model: all equations, parameters, and solver information. Key structural elements include the complete physics of the system, full parameter tables, the numerical methods, and anticipated failure modes with solutions (Supplementary Note 10). The complete guide is available in the repository given under Code availability and is also provided as a supplementary file. In all, the guide's components are made to be as explicit and precise as possible to mitigate any inference or hallucination on the agent's part. This mirrors the prompting literature, where explicit, structured, example-bearing instructions outperform terse ones 39 . Importantly, a proper specification can also be read and audited by a human such that the human’s physical intuition properly guides the development of the guide and constrains the multiphysics problem. The physics lay out the governing equations explicitly, not citing by name such as “standard Nernst–Planck”. The parameter tables list every constant, including its name, symbol, value, units, and source. The numerical methods in this study include a finite-volume discretization, using a Newton solver with damping (Supplementary Table 5). Bootstrap protocols for preconditioning and methods for initializing solvers are included to improve convergence for stiff problems such as the CO 2 R problem employed herein. 17 arXiv–ResearchArticle The methods were specified from physical insight and user experience building models in COMSOL, described in sufficient detail to prevent the agent from making any critical algorithmic decisions. We note that the significant degree of scaffolding used in this workflow is required for the current state of frontier models due to their tendency to invoke unphysical shortcuts or to hallucinate parameters 40 . However, purpose-built agents 41,42 or future frontier models may be able to handle greater flexibility in specification. Notably, this methodology is also transferable to other systems: a fuel cell or a Li-ion pseudo-two-dimensional model 43,44 guide would follow the same structure, with only the equations and parameters needing changes. Regardless, no specification is immune to ambiguity, which is why staged verification and guide refinement are essential. Having a human-in-the-loop for the development of the guide, at this stage, appears essential. Because the guide is a plain-text markdown file, it can be version-controlled and developed collaboratively. In this work, the implementation guide was maintained in a shared GitHub repository, where co-authors iteratively refined its physics, parameters, and numerical methods and shared the next refinement targets. The same versioned guide could be given to multiple researchers who independently instructed agents to build the model, where each independent model build was treated as a replicate. Divergent outputs exposed ambiguities or errors in the guide (Supplementary Tables 8 and 10), while agreement across independent builds gave confidence in both the model and the specification. This ability to collaborate on a single human- and machine-readable artifact, rather than on an opaque, version-locked solver project, is itself a central advantage of the workflow. Stagedconstructionandstopgates Building the model in stages of sequentially increasing complexity is essential to enable human review of output, identify and resolve issues before they accumulate and compound, and subsequently iterate upon the guide (Fig. 1). The guide decomposes the model into five stages of increasing physical complexity. In Stage 1, liquid Nernst–Planck, buffer reactions, and electrochemical kinetics are included with Dirichlet boundary conditions (BCs), initializing the problem with 400 degrees of freedom (DOFs) at the production mesh of 80 CL cells and 120 GDL cells. In Stage 2, Sherwood–Reynolds mass-transfer (Robin) BCs replace the Dirichlet to capture finite mass transfer at the electrolyte boundary. In Stage 3, Stefan–Maxwell gas transport and Henry's law phase transfer are added, expanding the system to 1,200 DOFs. In Stage 4, the Darcy pressure equation is finalized, completing the gas-phase momentum balance. A fitting step then estimates the kinetic parameters for CO 2 R and HER against the experimental data of Romiluyi et al. 25 , exploiting the well-conditioned Stage 4 residual before the coupling of the charged species to electrostatic potential causes stiffness to grow further. Stage 5 then turns on the full electrochemistry for the migration term in the Nernst–Planck flux, solving for liquid potential φ L monolithically from the per-cell ionic charge balance in the CL and solid potentialφ s in the GDL and CL via Ohm's law. This brings the full system to 1,480 DOFs at the 80-CL/120-GDL cell production mesh (Supplementary Table 2). The kinetic parameters fitted at Stage 4 are held fixed in Stage 5, where migration raises the simulated CO partial current (Fig. 2a). Crucially, the stages are designed to include stop gates where the human-in-the-loop can review quantitative physicality checks (e.g., checking electroneutrality to machine precision, mass 18 arXiv–ResearchArticle conservation, current continuity, and Newton residual convergence) before authorizing the next stage. These checks are prescribed by the guide itself rather than left to the discretion of the agent. For each stage, the guide specifies the quantitative pass/fail criteria as well as diagnostic plots the agent must produce (Supplementary Note 5), to enable criteria-driven and reproducible human review across builders. These output plots force data visualization, allowing the user to directly evaluate the model and accelerate diagnostics. In our case, a physical understanding of the expected general trends, including those in polarization behavior, faradaic efficiencies (FEs), and concentration profiles, makes these plots more interpretable than a list of numbers. A visual inspection to verify that trends are accurate, such as pH increasing with cathodic current and CO 2 concentration dropping with increasing overpotential, provides evidence that the model is building properly. This framework is essential for revealing failure modes in isolation and facilitating independent human verification of all model outputs, rather than relying solely on the agent’s self-assessment. For instance, an error in the buffer equilibria implementation affects local pH, and thus CO 2 concentration, reaction rates, liquid phase potential, etc., contaminating practically every part of the highly coupled problem. Locating errors in a liquid-only 400-DOF sub-problem is much simpler than allowing those same errors to get buried inside the full-physics 1,480-DOF system. Establishing stop gates addresses the issue of compounding complexity and allows the user to more effectively troubleshoot. Errors during model building could be due to mistakes made by the agent, but these stop gates also reveal inconsistencies or ambiguities in the guide, allowing the user to refine and iterate to make it more comprehensive and reproducible across instances of the agent. We note that frontier AI models are improving rapidly, even within the timespan of this study, and could eventually check their own work or work in teams of agents 18,19 to develop these rigorous models. At this stage, however, we find that building a model specification of this scale from scratch in one prompt is not practical, and scaffolding makes the workflow reliable and generalizable. Agentsandmodels Throughout guide development, we used the most capable Claude model available at the time because of the task's complexity. In earlier builds from less mature guides, weaker Sonnet-tier models introduced more errors than Opus, and the additional time spent correcting them, at comparable token costs, outweighed any benefit. Opus 4.6 was used to initialize guide development and perform much of the front-end refinement. Opus 4.7 was then primarily used for finalizing the model, and Opus 5 was used to make the final builds for comparison to COMSOL. Opus 4.8 was used for the planted-error study, and Opus 5 and Sonnet 5 for the fully autonomous builds. The choice to use a frontier model (Claude) rather than training our own purpose-built agent was made to first demonstrate the feasibility of the methodology with a well-established and powerful agent that, in theory, possesses intrinsic knowledge of solving similar nonlinear PDE problems. These agents were used by multiple users across laboratories, with varying experience levels, to develop the same implementation guide. When building models, researchers tracked their build processes in build dialogue files (Supplementary Note 7 and Supplementary Table 8) and 19 arXiv–ResearchArticle instructed their agents to refer only to the implementation guide (apart from the agent’s project memory, described below), so that its quality could be evaluated on its own and without bias from other builds or documents in the repository. In practice, more resources and context are beneficial for agents, and comparing replicates built from the same guide proved helpful in compensating for the intrinsic probabilistic nature of LLMs. The coding agent’s persistent project memory was active in all nine builds reported here. Its index, loaded at the start of each session, quoted results of builds made from an earlier version of the guide, including i 0 ≈106.5 A m −2 , a Stage-5 peak i CO ≈334 mA cm −2 and a note that a COMSOL run at that i 0 agreed to about 1%. It held no value ofλor the HER parameters, no polarization curve and no COMSOL output, and no build opened another build’s files or COMSOL output before it was complete. The guide itself states expected values at its stop gates, about twenty of them measured on the first of the nine builds. The agreement between builds does not rest on these values: λ, the HER parameters and the full 89-point polarization curves, which had no stated target, agree across builds as closely as i 0 does, and each fit starts from the guide’s prescribed initial point and is deterministic. A tenth build made in a directory containing only the guide, with no project memory, fitted i 0 = 106.59 A m −2 andλ= 1.6109 eV (Supplementary Note 8). Planted-errorstudy To perform the planted error study, specific defects were planted in the specification guide, not in the code, reproducing the failure mode of practical interest. In this failure mode, the statement of the physics handed to the implementer is itself wrong and is transcribed faithfully, propagating the error from the implementation guide to the codebase. Detection therefore requires the agent to recognize a conflict between what it has been told to implement and either the physics or the remainder of the specification and to raise that conflict rather than deviate silently. Each build received the same instruction: implement the model exactly as the guide specifies, consult no other file, halt at each stage gate for human review, and maintain a running build log recording every error encountered, every error corrected without human input, and every discrepancy noticed in the guide. The agent was not told that a defect existed or where. At each stage gate, the human reviewer inspected the plotted stage output and approved or rejected progression without auditing the source code. A defect was scored as detected only if the build log identified it as a defect, named the affected relation, and proposed a correction. Where a modified guide contained explanatory or diagnostic text naming the affected quantity, that text was edited along with the defect, so detection could not follow from an isolated self-contradiction. Per-case results are given in Supplementary Note 6 and Supplementary Table 7. Solve-timebenchmarkingandmeshconvergence Both implementations were timed on the same workstation (Intel Xeon w7-3565X, 32 cores/64 threads, 191 GB RAM) at matched core budget, with COMSOL on 32 cores (20,033 DOFs; Supplementary Table 1) and Julia on 32 threads (20,038 DOFs; Supplementary Note 4 and Supplementary Table 6). COMSOL swept 61 potentials at 50 mV spacing and Julia 89 at 25 mV, and both converged at every potential (Fig. 4a–b). Reported Julia times are the fastest of three 20 arXiv–ResearchArticle repetitions, with run-to-run variability reported in Supplementary Table 4, and the COMSOL time was the same in both of its recorded runs; process start-up and just-in-time compilation, a further 15–20 s per invocation, are excluded from both. The Julia data were obtained under the Windows Subsystem for Linux rather than natively, because the MUMPS binary distribution, used for the solver ablation, has no working Windows-native build. Mesh convergence of the Stage 5 solution is documented in Supplementary Note 4 and Supplementary Fig. 2. Reportingsummary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Dataavailability The experimental partial current densities used for kinetic parameter estimation are taken from Romiluyi et al. and are reproduced in Supplementary Table 3. The simulated data underlying Figs. 2, 3 and 6, and the COMSOL results they are compared with, are included in the code repository (see Code availability); the timings in Fig. 4 are listed in Supplementary Tables 4 and 6. Codeavailability The implementation guide, the Julia model source for all five build stages, and the build dialogue logs are available at https://github.com/Sebobbit/CO2R-Bulk-Scale-Model, with the Julia source released under the MIT license. The four modified guides used in the planted-error study, together with the unmodified control guide, are provided as separate files in the same repository. The COMSOL reference model is included; opening and running it requires a COMSOL Multiphysics license (v6.4) with the Electrochemistry and CFD modules, and its exported results are provided as plain-text files so that the comparison in Fig. 3 can be reproduced without a license. Acknowledgements The authors acknowledge professors Miguel A. Modestino and Rohini B. Chandran for helpful discussions. Funding S.C. discloses support for the research of this work from the National Science Foundation Graduate Research Fellowship Program [DGE-2234660]. J.C.B. discloses support for the research of this work from Schmidt Science Fellows, in partnership with the Rhodes Trust. E.W.L. acknowledges funding support from the Peter Wall Legacy Award and the Natural Sciences and Engineering Research Council of Canada Discovery Grants program (RGPIN-2025-04804). Any opinions, findings, and 21 arXiv–ResearchArticle conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. Authorcontributions S.C.: Investigation, Visualization, Writing – original draft, Writing – review & editing. M.F.S.: Investigation, Software. S.M.: Investigation. E.W.L.: Investigation, Supervision, Writing – review & editing. J.C.B.: Investigation, Supervision, Writing – review & editing. Competinginterests J.C.B., S.C., M.F.S. and E.W.L. are inventors on a provisional patent application filed by New York University and the University of British Columbia covering the specification-driven modeling workflow described here (US Provisional Application no. 64/168,700). The authors declare they have no other conflicts of interest. 22 arXiv–ResearchArticle References 1. Chu, S. & Majumdar, A. Opportunities and challenges for a sustainable energy future. Nature 488, 294–303 (2012). 2. De Luna, P. et al. What would it take for renewably powered electrosynthesis to displace petrochemical processes? Science 364, eaav3506 (2019). 3. Burdyny, T. & Smith, W. A. CO2 reduction on gas-diffusion electrodes and why catalytic performance must be assessed at commercially-relevant conditions. Energy Environ. Sci. 12, 1442–1453 (2019). 4. Welch, A. J. et al. Operando Local pH Measurement within Gas Diffusion Electrodes Performing Electrochemical Carbon Dioxide Reduction. J. Phys. Chem. C 125, 20896–20904 (2021). 5. Bui, J. C. et al. Continuum Modeling of Porous Electrodes for Electrochemical Synthesis. Chem. Rev. 122, 11022–11084 (2022). 6. Newman, J. & Balsara, N. P. Electrochemical Systems. (John Wiley & Sons, Incorporated, Newark, 2021). 7. Weng, L.-C., Bell, A. T. & Weber, A. Z. Modeling gas-diffusion electrodes for CO2 reduction. Phys. Chem. Chem. Phys. 20, 16973–16984 (2018). 8. Weng, L.-C., Bell, A. T. & Weber, A. Z. Towards membrane-electrode assembly systems for CO2 reduction: a modeling study. Energy Environ. Sci. 12, 1950–1968 (2019). 9. Ringe, S. et al. Understanding cation effects in electrochemical CO2 reduction. Energy Environ. Sci. 12, 3001–3014 (2019). 10. Newman, J. Numerical Solution of Coupled, Ordinary Differential Equations. Ind. Eng. Chem. Fundam. 7, 514–517 (1968). 11. Newman, J. S. & Tobias, C. W. Theoretical Analysis of Current Distribution in Porous Electrodes. J. Electrochem. Soc. 109, 1183 (1962). 12. Dickinson, E. J. F., Ekström, H. & Fontes, E. COMSOL Multiphysics®: Finite element software for electrochemical analysis. A mini-review. Electrochem. Commun. 40, 71–74 (2014). 13. Sulzer, V., Marquis, S. G., Timms, R., Robinson, M. & Chapman, S. J. Python Battery Mathematical Modelling (PyBaMM). J. Open Res. Softw. 9, (2021). 14. Permann, C. J. et al. MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX 11, 100430 (2020). 15. Peng, R. D. Reproducible Research in Computational Science. Science 334, 1226–1227 (2011). 16. M. Bran, A. et al. Augmenting large language models with chemistry tools. Nat. Mach. Intell. 6, 525–535 (2024). 17. Boiko, D. A., MacKnight, R., Kline, B. & Gomes, G. Autonomous chemical research with large language models. Nature 624, 570–578 (2023). 18. Ni, B. & Buehler, M. J. MechAgents: Large language model multi-agent collaborations can solve mechanics problems, generate new data, and integrate knowledge. Extreme Mech. Lett. 67, 102131 (2024). 19. Park, D., Moon, H. & Ryu, S. A self-correcting multi-agent LLM framework for language-based physics simulation and explanation. Npj Artif. Intell. 2, 10 (2026). 20. Lupoiu, R. et al. A multi-agentic framework for real-time, autonomous freeform metasurface 23 arXiv–ResearchArticle design. Sci. Adv. 11, eadx8006 (2025). 21. Li, S. et al. CodePDE: An Inference Framework for LLM-driven PDE Solver Generation. Preprint at https://doi.org/10.48550/arXiv.2505.08783 (2026). 22. Adhikari, S., Noorsumar, G. & Jensen, Ø. PDE-Agents: An LLM-Orchestrated Multi-Agent Framework for Automated Finite Element Simulations with Knowledge Graph-Augmented Reasoning. Preprint at https://doi.org/10.48550/arXiv.2606.07850 (2026). 23. Bezanson, J., Edelman, A., Karpinski, S. & Shah, V. B. Julia: A Fresh Approach to Numerical Computing. SIAM Rev. 59, 65–98 (2017). 24. Oberkampf, W. L. & Trucano, T. G. Verification and validation in computational fluid dynamics. Prog. Aerosp. Sci. 38, 209–272 (2002). 25. Romiluyi, O., Danilovic, N., Bell, A. T. & Weber, A. Z. Membrane-electrode assembly design parameters for optimal CO2 reduction. Electrochem. Sci. Adv. 3, e2100186 (2023). 26. Corpus, K. R. M. et al. Coupling covariance matrix adaptation with continuum modeling for determination of kinetic parameters associated with electrochemical CO2 reduction. Joule 7, 1289–1307 (2023). 27. Chidsey, C. E. D. Free Energy and Temperature Dependence of Electron Transfer at the Metal-Electrolyte Interface. Science 251, 919–922 (1991). 28. Nelder, J. A. & Mead, R. A Simplex Method for Function Minimization. Comput. J. 7, 308–313 (1965). 29. Mogensen, P. K. & Riseth, A. N. Optim: A mathematical optimization package for Julia. J. Open Source Softw. 3, 615 (2018). 30. Schenk, O. & Gärtner, K. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Gener. Comput. Syst. 20, 475–487 (2004). 31. Davis, T. A. Algorithm 832: UMFPACK V4.3---an unsymmetric-pattern multifrontal method. ACM Trans. Math. Softw. TOMS 30, 196–199 (2004). 32. Amestoy, P. R., Duff, I. S., L’Excellent, J.-Y. & Koster, J. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. SIAM J. Matrix Anal. Appl. 23, 15–41 (2001). 33. Revels, J., Lubin, M. & Papamarkou, T. Forward-Mode Automatic Differentiation in Julia. Preprint at https://doi.org/10.48550/arXiv.1607.07892 (2016). 34. Brady, N. W., Mees, M., Vereecken, P. M. & Safari, M. Implementation of Dual Number Automatic Differentiation with John Newman’s BAND Algorithm. J. Electrochem. Soc. 168, 113501 (2021). 35. Ouyang, S., Zhang, J. M., Harman, M. & Wang, M. An Empirical Study of the Non-Determinism of ChatGPT in Code Generation. ACM Trans. Softw. Eng. Methodol. 34, 42:1-42:28 (2025). 36. Lees, E. W., Bui, J. C., Romiluyi, O., Bell, A. T. & Weber, A. Z. Exploring CO2 reduction and crossover in membrane electrode assemblies. Nat. Chem. Eng. 1, 340–353 (2024). 37. Muhieddine Orfali, D. et al. Pulsed Electrolysis Promotes Catalyst Activity in Dilute CO2 Streams. ACS Energy Lett. 11, 4988–4997 (2026). 38. Zeng, Y., Smith, R. B., Bai, P. & Bazant, M. Z. Simple formula for Marcus–Hush–Chidsey kinetics. J. Electroanal. Chem. 735, 77–83 (2014). 39. Schulhoff, S. et al. The Prompt Report: A Systematic Survey of Prompt Engineering Techniques. Preprint at https://doi.org/10.48550/arXiv.2406.06608 (2025). 24 arXiv–ResearchArticle 40. Huang, L. et al. A Survey on Hallucination in Large Language Models: Principles, Taxonomy, Challenges, and Open Questions. ACM Trans. Inf. Syst. 43, 42:1-42:55 (2025). 41. Zou, Y. et al. El Agente: An autonomous agent for quantum chemistry. Matter 8, 102263 (2025). 42. Zhang, Z. et al. El Agente Forjador: Task-Driven Agent Generation for Quantum Simulation. Preprint at https://doi.org/10.48550/arXiv.2604.14609 (2026). 43. Weber, A. Z. & Newman, J. Modeling Transport in Polymer-Electrolyte Fuel Cells. Chem. Rev. 104, 4679–4726 (2004). 44. Doyle, M., Fuller, T. F. & Newman, J. Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell. J. Electrochem. Soc. 140, 1526 (1993). 25

Instant research discovery